{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0",
   "metadata": {},
   "source": [
    "# Advanced Computation\n",
    "\n",
    "This notebook demonstrates advanced signal processing techniques including anomaly detection, real-time monitoring, harmonic-percussive separation, and Empirical Mode Decomposition (EMD) using vitalDSP."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-1",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = \"sphinx_gallery\"\n",
    "from plotly import graph_objects as go\n",
    "from vitalDSP.notebooks import load_sample_ecg_small, load_sample_ppg, plot_trace\n",
    "\n",
    "fs = 128\n",
    "signal_col, date_col = load_sample_ecg_small()\n",
    "signal_col = np.array(signal_col)\n",
    "\n",
    "print(f\"ECG signal: {len(signal_col)} samples at {fs} Hz\")\n",
    "plot_trace(signal_col)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-3",
   "metadata": {},
   "source": [
    "## Anomaly Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.anomaly_detection import AnomalyDetection\n",
    "\n",
    "ad = AnomalyDetection(signal_col)\n",
    "\n",
    "# Z-score based anomaly detection\n",
    "anomalies_z = ad.detect_anomalies(method='z_score', threshold=3.0)\n",
    "print(f\"Z-score anomalies detected: {len(anomalies_z)}\")\n",
    "\n",
    "# IQR based\n",
    "anomalies_iqr = ad.detect_anomalies(method='iqr')\n",
    "print(f\"IQR anomalies detected    : {len(anomalies_iqr)}\")\n",
    "\n",
    "# Visualize anomalies\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=signal_col, mode='lines', name='ECG Signal', opacity=0.7))\n",
    "if len(anomalies_z) > 0:\n",
    "    fig.add_trace(go.Scatter(\n",
    "        x=anomalies_z,\n",
    "        y=signal_col[anomalies_z],\n",
    "        mode='markers', marker=dict(color='red', size=8, symbol='circle-open'),\n",
    "        name='Anomalies (Z-score)'\n",
    "    ))\n",
    "fig.update_layout(title='ECG Anomaly Detection', xaxis_title='Sample', yaxis_title='Amplitude')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-5",
   "metadata": {},
   "source": [
    "## Real-Time Anomaly Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.real_time_anomaly_detection import RealTimeAnomalyDetection\n",
    "\n",
    "# Initialize with a sliding window size\n",
    "rtad = RealTimeAnomalyDetection(window_size=64)\n",
    "\n",
    "# Train autoencoder on the first 512 samples as normal baseline\n",
    "training_data = signal_col[:512]\n",
    "rtad.train_autoencoder(training_data, encoding_dim=8)\n",
    "\n",
    "# Stream test samples and detect anomalies\n",
    "test_segment = signal_col[512:640]\n",
    "anomaly_flags = [rtad.detect_autoencoder(np.array([s]), threshold=0.1) for s in test_segment]\n",
    "stat_flags    = [rtad.detect_statistical(np.array([s]), method='z_score', threshold=2.5) for s in test_segment]\n",
    "\n",
    "print(f\"Autoencoder anomalies in test segment: {sum(anomaly_flags)} / {len(test_segment)} samples\")\n",
    "print(f\"Statistical anomalies                : {sum(stat_flags)} / {len(test_segment)} samples\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-7",
   "metadata": {},
   "source": [
    "## Harmonic-Percussive Source Separation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.harmonic_percussive_separation import HarmonicPercussiveSeparation\n",
    "\n",
    "hps = HarmonicPercussiveSeparation(signal_col)\n",
    "harmonic, percussive = hps.separate()\n",
    "\n",
    "print(f\"Harmonic component  : shape={harmonic.shape}, power={np.mean(harmonic**2):.4f}\")\n",
    "print(f\"Percussive component: shape={percussive.shape}, power={np.mean(percussive**2):.4f}\")\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='Original', opacity=0.4))\n",
    "fig.add_trace(go.Scatter(y=harmonic[:512], mode='lines', name='Harmonic (smooth)'))\n",
    "fig.add_trace(go.Scatter(y=percussive[:512], mode='lines', name='Percussive (transients)'))\n",
    "fig.update_layout(\n",
    "    title='Harmonic-Percussive Separation',\n",
    "    xaxis_title='Sample', yaxis_title='Amplitude'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-9",
   "metadata": {},
   "source": [
    "## Empirical Mode Decomposition (EMD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.emd import EMD\n",
    "\n",
    "# Use a shorter segment for EMD (computationally intensive)\n",
    "short_signal = signal_col[:512]\n",
    "\n",
    "emd = EMD(short_signal)\n",
    "imfs = emd.emd()\n",
    "\n",
    "print(f\"Number of IMFs extracted: {len(imfs)}\")\n",
    "for i, imf in enumerate(imfs):\n",
    "    print(f\"  IMF {i+1}: mean={np.mean(imf):.4f}, std={np.std(imf):.4f}\")\n",
    "\n",
    "# Plot IMFs\n",
    "n_plot = min(4, len(imfs))\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=short_signal, mode='lines', name='Original', opacity=0.4))\n",
    "for i in range(n_plot):\n",
    "    fig.add_trace(go.Scatter(y=imfs[i], mode='lines', name=f'IMF {i+1}'))\n",
    "fig.update_layout(\n",
    "    title='EMD: Intrinsic Mode Functions',\n",
    "    xaxis_title='Sample', yaxis_title='Amplitude'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-11",
   "metadata": {},
   "source": [
    "## Advanced Kalman Filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.kalman_filter import KalmanFilter\n",
    "\n",
    "# 1-D Kalman filter with identity transition / measurement matrices\n",
    "kf = KalmanFilter(\n",
    "    initial_state=np.array([[signal_col[0]]]),\n",
    "    initial_covariance=np.array([[1.0]]),\n",
    "    process_covariance=np.array([[1e-3]]),\n",
    "    measurement_covariance=np.array([[0.1]])\n",
    ")\n",
    "\n",
    "H = np.array([[1.0]])  # measurement matrix\n",
    "F = np.array([[1.0]])  # state transition matrix\n",
    "\n",
    "kf_filtered = kf.filter(signal_col, measurement_matrix=H, transition_matrix=F)\n",
    "\n",
    "plot_trace(signal_col, kf_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-13",
   "metadata": {},
   "source": [
    "## Sparse Signal Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.sparse_signal_processing import SparseSignalProcessing\n",
    "\n",
    "ssp = SparseSignalProcessing(signal_col)\n",
    "\n",
    "# Build a DCT basis for sparse representation\n",
    "N = len(signal_col)\n",
    "k = np.arange(N)\n",
    "basis = np.cos(np.pi * k[:, None] * (2 * k[None, :] + 1) / (2 * N)) * np.sqrt(2 / N)\n",
    "basis[0] *= 1 / np.sqrt(2)\n",
    "\n",
    "sparse_repr = ssp.sparse_representation(basis)\n",
    "\n",
    "# Threshold: keep only 5% of coefficients\n",
    "threshold = np.percentile(np.abs(sparse_repr), 95)\n",
    "sparse_thresh = ssp.thresholding(sparse_repr, threshold=threshold)\n",
    "\n",
    "reconstructed = ssp.reconstruction(sparse_thresh, inverse_basis=basis.T)\n",
    "print(f\"Sparse reconstruction MSE (5% kept): {np.mean((signal_col - reconstructed)**2):.6f}\")\n",
    "\n",
    "plot_trace(signal_col, reconstructed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Process Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.bayesian_analysis import GaussianProcess\n\n# Fit a GP to a short ECG segment as a smooth interpolator\nn_train = 40\nX_train = np.arange(n_train).reshape(-1, 1).astype(float)\ny_train = signal_col[:n_train]\n\ngp = GaussianProcess(length_scale=5.0, noise=0.05)\ngp.update(X_train, y_train)\n\nX_test = np.linspace(0, n_train - 1, 200).reshape(-1, 1)\nmean_pred, var_pred = gp.predict(X_test)\nstd_pred = np.sqrt(var_pred)\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(x=X_train.ravel(), y=y_train,\n                         mode='markers', name='Training points'))\nfig.add_trace(go.Scatter(x=X_test.ravel(), y=mean_pred,\n                         mode='lines', name='GP mean prediction'))\nfig.add_trace(go.Scatter(\n    x=np.concatenate([X_test.ravel(), X_test.ravel()[::-1]]),\n    y=np.concatenate([mean_pred + 2*std_pred, (mean_pred - 2*std_pred)[::-1]]),\n    fill='toself', opacity=0.2, name='95% confidence'\n))\nfig.update_layout(title='Gaussian Process Regression on ECG Segment',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.bayesian_analysis import BayesianOptimization\n\n# Optimize a synthetic objective (e.g., finding best filter cutoff)\ndef signal_snr_objective(cutoff):\n    from vitalDSP.filtering.signal_filtering import SignalFiltering\n    filtered = SignalFiltering(signal_col).lowpass(cutoff=cutoff, fs=fs, order=4)\n    noise = signal_col - filtered\n    snr = 10 * np.log10(np.mean(filtered**2) / (np.mean(noise**2) + 1e-8))\n    return snr\n\nbo = BayesianOptimization(func=signal_snr_objective, bounds=(5.0, 45.0))\nbest_cutoff, best_snr = bo.optimize(n_iter=8)\nprint(f'Optimal lowpass cutoff: {float(best_cutoff):.2f} Hz (SNR={best_snr:.2f} dB)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generative Signal Synthesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.generative_signal_synthesis import GenerativeSignalSynthesis\n\ngss = GenerativeSignalSynthesis()\n\nsynth_ecg  = gss.generate(method='random_noise', length=512)\nsynth_ar   = gss.generate(method='autoregressive', length=512, coeffs=[0.9, -0.2])\n\nprint(f'Synthetic noise signal : {synth_ecg.shape}')\nprint(f'Synthetic sine signal  : {synth_ar.shape}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=synth_ecg, mode='lines', name='Random noise synthesis'))\nfig.add_trace(go.Scatter(y=synth_ar, mode='lines', name='Autoregressive synthesis'))\nfig.update_layout(title='Generative Signal Synthesis',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multimodal Signal Fusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.multimodal_fusion import MultimodalFusion\nfrom vitalDSP.notebooks import load_sample_ppg\n\nppg_col, _ = load_sample_ppg()\nppg_col = np.array(ppg_col)\nn = min(len(signal_col), len(ppg_col))\n\n# Normalize both signals before fusion\necg_norm = (signal_col[:n] - signal_col[:n].mean()) / (signal_col[:n].std() + 1e-8)\nppg_norm = (ppg_col[:n] - ppg_col[:n].mean()) / (ppg_col[:n].std() + 1e-8)\n\nmf = MultimodalFusion(signals=[ecg_norm, ppg_norm])\n\nfused_ws  = mf.fuse(strategy='weighted_sum', weights=[0.6, 0.4])\nfused_pca = mf.fuse(strategy='pca', n_components=1)\nfused_max = mf.fuse(strategy='maximum')\n\nprint(f'Weighted sum fusion : {fused_ws.shape}')\nprint(f'PCA fusion          : {fused_pca.shape}')\nprint(f'Maximum fusion      : {fused_max.shape}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=ecg_norm[:512], mode='lines', name='ECG (norm)', opacity=0.5))\nfig.add_trace(go.Scatter(y=ppg_norm[:512], mode='lines', name='PPG (norm)', opacity=0.5))\nfig.add_trace(go.Scatter(y=fused_ws[:512], mode='lines', name='Weighted Sum (0.6/0.4)'))\nfig.add_trace(go.Scatter(y=fused_pca[:512], mode='lines', name='PCA Fusion'))\nfig.update_layout(title='Multimodal ECG-PPG Fusion',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural Network Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.advanced_computation.neural_network_filtering import NeuralNetworkFiltering\n\n# Train a feedforward network to learn the signal structure\n# and apply it as a learned filter\nnnf = NeuralNetworkFiltering(\n    signal_col[:512],\n    network_type='feedforward',\n    hidden_layers=[32, 16],\n    epochs=5,\n    batch_size=32\n)\nnnf.train()\nfiltered_nn = nnf.apply_filter().ravel()\n\nprint(f'NN filtered signal shape: {filtered_nn.shape}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=signal_col[1:512], mode='lines', name='Original ECG', opacity=0.5))\nfig.add_trace(go.Scatter(y=filtered_nn, mode='lines', name='NN Filtered'))\nfig.update_layout(title='Neural Network Filtering',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}