1.general_framework.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "754c1d64",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2023-11-10T10:40:49.539850Z"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "from joblib import Parallel, delayed\n",
    "from tqdm import tqdm\n",
    "from importlib import reload\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from seaborn import heatmap, color_palette"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb44384e-7cbc-4611-a015-e1c073f61aa9",
   "metadata": {},
   "source": [
    "# Datasets loading\n",
    "\n",
    "Lots of different on availabale : https://towardsdatascience.com/a-data-lakes-worth-of-audio-datasets-b45b88cd4ad\n",
    "\n",
    "Review: \n",
    "https://arxiv.org/abs/2012.02974\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00291787-4082-48ea-9002-54d777cac5b6",
   "metadata": {},
   "source": [
    "## Prediction ahead\n",
    "\n",
    "Datasets available :\n",
    "\n",
    "* ReservoirPy: MackeyGlass, Lorenz\n",
    "* Custom: Sunspot"
   ]
  },
  {
   "cell_type": "raw",
   "id": "099fb604-5450-495e-8e6b-b6ff9eb78ff4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-20T23:34:32.438023Z",
     "iopub.status.busy": "2024-11-20T23:34:32.437367Z",
     "iopub.status.idle": "2024-11-20T23:34:32.885775Z",
     "shell.execute_reply": "2024-11-20T23:34:32.885504Z",
     "shell.execute_reply.started": "2024-11-20T23:34:32.437978Z"
    }
   },
   "source": [
    "from datasets.load_forecasting import load_dataset_forecasting\n",
    "is_instances_classification = False\n",
    "dataset_name = \"Lorenz\"\n",
    "step_ahead=5\n",
    "\n",
    "is_multivariate, sampling_rate, X_train_raw, X_test_raw, Y_train_raw, Y_test = load_dataset_forecasting(dataset_name, step_ahead, visualize=True)\n",
    "use_spectral_representation = False\n",
    "spectral_representation = None # Can be None, \"stft\" or \"mfcc\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ae6a7ebe247f27",
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a741331a-93ba-4132-8449-0e9521c4ddcc",
   "metadata": {},
   "source": [
    "Datasets available :\n",
    "\n",
    "* Custom :  FSDD, HAART, JapaneseVowels\n",
    "* Aeon : SpokenArabicDigits, CatsDogs, LSST\n",
    "* Torchaudio: SPEECHCOMMANDS\n",
    "\n",
    "More on https://www.timeseriesclassification.com/dataset.php or https://pytorch.org/audio/stable/datasets.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1e39882-bed2-4c3a-97f9-962a9e7247f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets.load_classification import load_dataset_classification\n",
    "is_instances_classification = True\n",
    "spectral_representation = \"mfcc\"  # Can be None, \"stft\" or \"mfcc\"\n",
    "\n",
    "dataset_name = \"FSDD\"\n",
    "\n",
    "use_spectral_representation, is_multivariate, sampling_rate, X_train_raw, X_test_raw, Y_train_raw, Y_test, groups = load_dataset_classification(dataset_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d5ddde-f921-468c-abfc-cb56d1255be3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets.preprocessing import plot_classes_distribution\n",
    "\n",
    "# Plot data distribution\n",
    "if is_instances_classification:\n",
    "    plot_classes_distribution(Y_train_raw, Y_test, val = False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af55bde5-5bb3-41ea-ae0d-10d7bc7401f4",
   "metadata": {},
   "source": [
    "# Cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90237f4e-2418-48a9-a067-aece8c922912",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold, TimeSeriesSplit, StratifiedGroupKFold\n",
    "from datasets.preprocessing import flexible_indexing, plot_classes_distribution\n",
    "\n",
    "# CROSS-VALIDATION METHODS\n",
    "# SEED\n",
    "SEED = 49387\n",
    "\n",
    "use_cross_validation = False\n",
    "\n",
    "n_splits=2\n",
    "\n",
    "if use_cross_validation: # we split train between train and val\n",
    "    if is_instances_classification:\n",
    "        if groups is None:\n",
    "            splits = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=SEED).split(X_train_raw, np.argmax(Y_train_raw, axis=1))\n",
    "        else:\n",
    "            splits = StratifiedGroupKFold(n_splits=n_splits, shuffle=True, random_state=SEED).split(X_train_raw, np.argmax(Y_train_raw, axis=1), groups)\n",
    "    else: #prediction\n",
    "        splits = TimeSeriesSplit(n_splits=n_splits).split(X_train_raw)\n",
    "    \n",
    "    for i, (train_index, val_index) in enumerate(splits):\n",
    "        X_train = flexible_indexing(X_train_raw, train_index)\n",
    "        X_val = flexible_indexing(X_train_raw, val_index)\n",
    "        Y_train = flexible_indexing(Y_train_raw, train_index)\n",
    "        Y_val = flexible_indexing(Y_train_raw, val_index)\n",
    "        # SPLITS\n",
    "        if is_multivariate:\n",
    "            X_train_band, X_val_band, X_test_band = X_train, X_val, X_test_raw\n",
    "            del X_train, X_val\n",
    "        else:\n",
    "            X_test = X_test_raw\n",
    "\n",
    "    if is_instances_classification:\n",
    "        plot_classes_distribution(Y_train, Y_val, val = True)\n",
    "else: # then we use the test dataset\n",
    "    if is_multivariate:\n",
    "        X_train_band, X_test_band = X_train_raw, X_test_raw\n",
    "        X_val_band = None\n",
    "    else:\n",
    "        X_test, X_train = X_test_raw, X_train_raw\n",
    "        X_val, X_val_band = None, None\n",
    "    Y_train = Y_train_raw\n",
    "#del Y_train_raw, X_train_raw, X_test_raw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a468ca5c",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "# Preprocessing\n",
    "\n",
    "In this part we convert (when needed) the univariatedatasets into multivariate ones\n",
    "\n",
    "Spectrograms_vs_Cochleagrams : https://www.researchgate.net/publication/340510607_Speech_recognition_using_very_deep_neural_networks_Spectrograms_vs_Cochleagrams"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "766304bc-3c47-4813-9828-91d8fed5845b",
   "metadata": {},
   "source": [
    "## Multivariate generation (if not multivariate) and train_validation split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "380ad0f4-4dbc-4b97-8d28-8ffb357715d7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-20T11:11:18.377724Z",
     "start_time": "2023-10-20T11:11:18.143305Z"
    }
   },
   "outputs": [],
   "source": [
    "from datasets.multivariate_generation import generate_multivariate_dataset, extract_peak_frequencies\n",
    "\n",
    "freq_train_data = X_train_band if is_multivariate else X_train\n",
    "flat_train_data = np.concatenate(freq_train_data, axis=0) if is_instances_classification else freq_train_data\n",
    "\n",
    "WINDOW_LENGTH = 10\n",
    "\n",
    "smoothed_peaks_freqs = extract_peak_frequencies(flat_train_data, sampling_rate, smooth=True, window_length=WINDOW_LENGTH, threshold=1e-5, nperseg=1024, visualize=True)\n",
    "nb_smoothed = len(smoothed_peaks_freqs[0]) if is_multivariate else len(smoothed_peaks_freqs)\n",
    "print(\"Smoothed peaks: \", nb_smoothed)\n",
    "\n",
    "peak_freqs = extract_peak_frequencies(flat_train_data, sampling_rate, smooth=True, window_length=1, threshold=1e-5, nperseg=1024, visualize=True)\n",
    "nb_peaks = len(peak_freqs[0]) if is_multivariate else len(peak_freqs)\n",
    "print(\"Non-smoothed peaks: \", nb_peaks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc3488bf-2624-44aa-9ca4-7d96c975c392",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets.multivariate_generation import generate_multivariate_dataset, extract_peak_frequencies\n",
    "\n",
    "# if it has use_spectral_representation, then it is multivariate\n",
    "if use_spectral_representation == True: \n",
    "    if is_multivariate==False:\n",
    "        raise ValueError(\"Cannot use spectral representation if it's not multivariate !\")\n",
    "\n",
    "if not is_multivariate:\n",
    "    print(f\"Converting single variate to {spectral_representation}\")\n",
    "    X_train_band = generate_multivariate_dataset(\n",
    "        X_train, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "    )\n",
    "    if use_cross_validation:\n",
    "        X_val_band = generate_multivariate_dataset(\n",
    "            X_val, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "        )\n",
    "    X_test_band = generate_multivariate_dataset(\n",
    "            X_test, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "    )\n",
    "        \n",
    "elif is_multivariate and not use_spectral_representation:\n",
    "    print(f\"Converting multivariate to {spectral_representation}\") # if None we convert to temporal\n",
    "    X_train_band = generate_multivariate_dataset(\n",
    "        X_train_band, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "    )\n",
    "    if use_cross_validation:\n",
    "        X_val_band = generate_multivariate_dataset(\n",
    "            X_val_band, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "        )\n",
    "    X_test_band = generate_multivariate_dataset(\n",
    "        X_test_band, sampling_rate, is_instances_classification, peak_freqs, spectral_representation, hop=100\n",
    "    )\n",
    "else:\n",
    "    print(\"Data is already spectral, nothing to do\")\n",
    "\n",
    "if spectral_representation is not None:\n",
    "    use_spectral_representation = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6b859e5-c6dd-4086-ac56-42a7de3ed25e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Standardizing the amplitudes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb45e001-4920-47e1-8dbb-6681cf9b32e9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-20T11:12:00.856649Z",
     "start_time": "2023-10-20T11:11:39.230512Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from datasets.preprocessing import scale_data\n",
    "\n",
    "X_train_band_not_scale = X_train_band\n",
    "scaler_multi = MinMaxScaler(feature_range=(0, 1))\n",
    "X_train_band, X_val_band, X_test_band = scale_data(X_train_band, X_val_band, X_test_band, scaler_multi, is_instances_classification)\n",
    "            \n",
    "if not is_multivariate:\n",
    "    X_train_not_scale = X_train\n",
    "    scaler_x_uni = MinMaxScaler(feature_range=(0, 1))\n",
    "    X_train, X_val, X_test = scale_data(X_train, X_val, X_test, scaler_multi, is_instances_classification)       \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ddcde6-215a-464a-8634-624f11da4380",
   "metadata": {},
   "source": [
    "## Visualisation checks\n",
    "\n",
    "test if the filtering happened correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87795848-f154-419a-b7d4-35400767f51d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from seaborn import color_palette\n",
    "import librosa\n",
    "import librosa.display\n",
    "from matplotlib.patches import ConnectionPatch\n",
    "\n",
    "# Define different color palettes for each step\n",
    "colors_step1 = color_palette(\"tab20\")    # For the first plot\n",
    "colors_step2 = color_palette(\"Set2\")     # For the second plot\n",
    "colors_step3 = color_palette(\"Dark2\")     # For the third plot\n",
    "fontsize = 14\n",
    "\n",
    "# Create a figure with three subplots side by side\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "\n",
    "# First Plot: Original Time Series\n",
    "if is_instances_classification:\n",
    "    train_data = np.concatenate(X_train, axis=0).T\n",
    "else:\n",
    "    train_data = X_train_raw.T\n",
    "\n",
    "START = 0\n",
    "END = 500\n",
    "DIFF = END - START\n",
    "\n",
    "ax = axes[0]\n",
    "for idx, i in enumerate(range(train_data.shape[0])):\n",
    "    ax.plot(range(DIFF), train_data[i, START:END], color=colors_step1[idx % len(colors_step1)])\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "ax.set_xlabel('Time', fontsize=fontsize)\n",
    "#ax.set_ylabel('Value', fontsize=fontsize)\n",
    "ax.set_title('Original Time Series Data', fontsize=fontsize)\n",
    "\n",
    "# Second Plot: Transformed Data\n",
    "if is_instances_classification:\n",
    "    train_band_data = np.concatenate(X_train_band_not_scale, axis=0).T\n",
    "else:\n",
    "    train_band_data = X_train_band_not_scale.T\n",
    "    \n",
    "ax = axes[1]\n",
    "win_length = 50\n",
    "freqs = librosa.fft_frequencies(sr=sampling_rate, n_fft=win_length)\n",
    "\n",
    "if use_spectral_representation:\n",
    "    END = 200\n",
    "    DIFF = END - START\n",
    "    S_db = librosa.amplitude_to_db(train_band_data, ref=np.max)[:,START:END]\n",
    "\n",
    "    # Plot the spectrogram using imshow for manual control\n",
    "    img = ax.imshow(\n",
    "        S_db,\n",
    "        aspect='auto',\n",
    "        origin='lower',\n",
    "        interpolation='none',\n",
    "        cmap='plasma'\n",
    "    )\n",
    "    \n",
    "    # Adjust the yticks to show one out of every two labels\n",
    "    num_freq_bins = S_db.shape[0]\n",
    "    ytick_labels = [f\"{int(freqs[i])}\" if i % 2 == 0 else \"\" for i in range(num_freq_bins)]\n",
    "    ax.set_yticks(np.arange(num_freq_bins))  # Keep all the ticks\n",
    "    ax.set_yticklabels(ytick_labels, fontsize=fontsize)  # Set custom labels\n",
    "    ax.tick_params(axis='both', labelsize=fontsize)\n",
    "\n",
    "    ax.set_xlabel('Steps', fontsize=fontsize)\n",
    "    ax.set_title('Transformed spectral Data', fontsize=fontsize)\n",
    "    fig.colorbar(img, ax=ax, format=\"%+2.0f dB\")\n",
    "else:\n",
    "    random_indexes = np.random.randint(0, train_band_data.shape[0], size=3)\n",
    "    for idx, i in enumerate(random_indexes):\n",
    "        ax.plot(range(DIFF), train_band_data[i, START:END],\n",
    "                color=colors_step2[idx % len(colors_step2)], label=f'Feature {i}')\n",
    "\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.tick_params(axis='both', labelsize=fontsize)\n",
    "    ax.set_xlabel('Time', fontsize=fontsize)\n",
    "    ax.set_title('Transformed spectral Data', fontsize=fontsize)\n",
    "    ax.legend(fontsize=fontsize, loc='upper right')\n",
    "\n",
    "# Third Plot: Normalized Data\n",
    "if is_instances_classification:\n",
    "    train_band_data = np.concatenate(X_train_band, axis=0).T\n",
    "else:\n",
    "    train_band_data = X_train_band.T\n",
    "ax = axes[2]\n",
    "if 'random_indexes' not in locals():\n",
    "    random_indexes = np.random.randint(0, train_band_data.shape[0], size=3)\n",
    "    \n",
    "for idx, i in enumerate(random_indexes):\n",
    "    ax.plot(range(DIFF), train_band_data[i, START:END], color=colors_step3[idx % len(colors_step3)],\n",
    "            label=f'Feature {i}')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "ax.set_xlabel('Steps', fontsize=fontsize)\n",
    "#ax.set_ylabel('Normalized Value', fontsize=fontsize)\n",
    "ax.set_title('Normalized spectral Data', fontsize=fontsize)\n",
    "ax.legend(fontsize=fontsize, loc='upper right')\n",
    "\n",
    "# Add arrows between subplots to indicate transformation\n",
    "fig.tight_layout()\n",
    "\n",
    "# Draw arrows using ConnectionPatch for better control\n",
    "def add_arrow(axes_from, axes_to, x_from, x_to):\n",
    "    con = ConnectionPatch(\n",
    "        xyA=(x_from, 0.5), xyB=(x_to, 0.5),\n",
    "        coordsA='axes fraction', coordsB='axes fraction',\n",
    "        axesA=axes_from, axesB=axes_to,\n",
    "        arrowstyle=\"-|>\", color=\"k\", shrinkB=5\n",
    "    )\n",
    "    axes_from.add_artist(con)\n",
    "\n",
    "add_arrow(axes[0], axes[1], 0.94, -0.13)\n",
    "add_arrow(axes[1], axes[2], 1.25, -0.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bde977ab-e2ee-43a7-a916-3ae2f4a95780",
   "metadata": {},
   "source": [
    "## Noizing\n",
    "\n",
    "This is to optionnally add noize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b49a01ef-7b49-4493-aa96-3ea65b18b118",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define noise parameter\n",
    "noise_std = 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97747959-8de6-4fdc-85e8-f8e6425b3399",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets.preprocessing import add_noise, duplicate_data\n",
    "\n",
    "#Train/Val/Test\n",
    "if is_instances_classification:\n",
    "    # UNI\n",
    "    if not is_multivariate:\n",
    "        X_train_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_train, desc=\"TRAIN UNI\")]\n",
    "        if X_val is not None:\n",
    "            X_val_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_val, desc=\"VAL UNI\")]\n",
    "        X_test_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_test, desc=\"TEST UNI\")]\n",
    "        \n",
    "    # MULTIX\n",
    "    X_train_band_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_train_band, desc=\"TRAIN MULTI\")]\n",
    "    if X_val_band is not None:\n",
    "        X_val_band_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_val_band, desc=\"VAL MULTI\")]\n",
    "    X_test_band_noisy = [add_noise(instance, noise_std) for instance in tqdm(X_test_band, desc=\"TEST MULTI\")]\n",
    "\n",
    "else:  #if prediction\n",
    "    # UNI\n",
    "    if not is_multivariate:\n",
    "        X_train_noisy = add_noise(X_train, noise_std)\n",
    "        X_test_noisy = add_noise(X_test, noise_std)\n",
    "        if X_val is not None:\n",
    "            X_val_noisy = add_noise(X_val, noise_std)\n",
    "\n",
    "    # MULTI\n",
    "    X_train_band_noisy = add_noise(X_train_band, noise_std)\n",
    "    if X_val_band is not None:\n",
    "        X_val_band_noisy = add_noise(X_val_band, noise_std)\n",
    "    X_test_band_noisy = add_noise(X_test_band, noise_std)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25d6f51b-0c6c-4cec-bbd5-1fff0223efed",
   "metadata": {},
   "source": [
    "## Pretrain\n",
    "\n",
    "We use the train dataset as pretraining (pretraining is an unsupervised task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f5fc6e-c2cf-4441-817a-6ba38f14e2de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the number of instances you want to select\n",
    "x_size = len(X_train_band) if is_multivariate else len(X_train)\n",
    "num_samples_for_pretrain = 500 if x_size >= 500 else x_size\n",
    "if is_instances_classification:\n",
    "    indices = np.random.choice(x_size, num_samples_for_pretrain, replace=False)\n",
    "else:\n",
    "    indices = range(x_size)\n",
    "    \n",
    "if not is_multivariate:\n",
    "    X_pretrain_uni = np.array(X_train, dtype=object)[indices]\n",
    "    X_pretrain_noisy = np.array(X_train_noisy, dtype=object)[indices]\n",
    "X_pretrain_band = np.array(X_train_band, dtype=object)[indices]\n",
    "X_pretrain_band_noisy = np.array(X_train_band_noisy, dtype=object)[indices]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd93f662-f197-4b0a-a140-090f3c0909d2",
   "metadata": {},
   "source": [
    "# Generating reservoirs\n",
    "\n",
    "We are interrested in two technique to genereate reservoir. \n",
    "* One is called mean-HAG (hadsp), was studied in previous paper, and recombines inputs based on their activity to reach a given activation target.\n",
    "* The other called variance-HAG (desp) recombines inputs in order to reach a given standard deviation of activation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b33224-cc12-4edd-948c-235f415a4ab8",
   "metadata": {},
   "source": [
    "## Data format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aa63393-5c4e-4cce-9ed0-7caf47af8589",
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_instances_classification:\n",
    "    common_index = 1\n",
    "    print(\"Common index for multivariate classification should be 1\")\n",
    "    print(\"\\nCheck it ! \\nFirst array \", X_train_band[1].shape, \" and second array\", X_train_band[2].shape)\n",
    "    print(\"\\nCheck it ! \\nFirst array \", X_test_band[1].shape, \" and second array\", X_test_band[2].shape)\n",
    "    common_size = X_train_band[0].shape[common_index]\n",
    "else:\n",
    "    common_index = 1\n",
    "    print(\"Common index for multivariate prediction should be 1\")\n",
    "    print(\"\\nCheck it ! \\nFirst array \", X_train_band.shape, \" and second array\", X_train_band.shape)\n",
    "    common_size = X_train_band.shape[common_index]\n",
    "\n",
    "print(\"Common size:\", common_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24988d03-a16b-46f7-b3f1-372452272905",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Min window size to get all the dynamics ? \n",
    "min_window_size = sampling_rate/np.max(np.hstack(peak_freqs))\n",
    "max_window_size = sampling_rate/np.min(np.hstack(peak_freqs))\n",
    "\n",
    "print(min_window_size)\n",
    "print(max_window_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d6435e2-6f19-4987-9309-ab74694678fd",
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Shared parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd8001c-ba74-457f-9518-c0afa04f3606",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Reservoir parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "615a38d2df41a727",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from reservoir.activation_functions import tanh, heaviside, sigmoid, softplus, heaviside\n",
    "\n",
    "# the activation function choosen for the rest of the experiment\n",
    "# activation_function = lambda x : sigmoid(2*(x-0.5))tanh(x)\n",
    "activation_function = lambda x : tanh(x)\n",
    "\n",
    "x=np.linspace(0, 2, 100)\n",
    "plt.plot(x, activation_function(x))\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b4dc28-34d5-4a98-a047-7fa6e4e9da47",
   "metadata": {},
   "source": [
    "**common_size** : the number of different dimensions in the input data\n",
    "\n",
    "**K** : the number of neurons that will receive a particular time serie as input "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2433d6d3-acc2-4bbc-8373-12d692eda119",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math \n",
    "\n",
    "RESERVOIR_SIZE = 500\n",
    "\n",
    "# We want the size of the reservoir to be at least RESERVOIR_SIZE\n",
    "K = math.ceil(RESERVOIR_SIZE / common_size)\n",
    "n = common_size * K\n",
    "print(\"Dimension of our reservoir :\", n)\n",
    "print(\"Copy of each time serie :\", K)\n",
    "print(\"Number of multivariate inputs :\", common_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "970f3072-3e8d-49f1-a890-5c3a28df3d9e",
   "metadata": {},
   "source": [
    "### Hyperparameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bbc2c01-de54-4bf6-8fa0-1dbde5fe21c2",
   "metadata": {},
   "source": [
    "#### Customs"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a6c7bf77-e19f-4790-b401-a40d219d8e7c",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "time_increment = 82 # int(max_window_size) or int(min_window_size+1)\n",
    "time_increment_span = 14\n",
    "MAX_TIME_INCREMENT = TIME_INCREMENT + time_increment_span #int(max_window_size) or None  or TIME_INCREMENT\n",
    "weight_increment = 0.030000000000000002\n",
    "\n",
    "target_rate = 0.58\n",
    "RATE_SPREAD = 0.34\n",
    "\n",
    "min_variance = 0.014000000000000002\n",
    "variance_spread = 0.003\n",
    "\n",
    "bias_scaling = 0.0\n",
    "input_scaling = 0.060000000000000005\n",
    "\n",
    "intrinsic_saturation = 0.92\n",
    "intrinsic_coef = 0.8\n",
    "\n",
    "max_partners= 20\n",
    "use_full_instance = False\n",
    "\n",
    "# Doesn't move :\n",
    "input_connectivity = 1\n",
    "connectivity = 0\n",
    "leaky_rate = 1\n",
    "if int(max_window_size) < TIME_INCREMENT or TIME_INCREMENT < min_window_size:\n",
    "    raise ValueError(f\"INCREMENT must be greater than {min_window_size} and smaller than {max_window_size}. Current INCREMENT is {TIME_INCREMENT}.\")\n",
    "\n",
    "TIME_INCREMENT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e7ba684-539a-4f94-99f0-9bf242717c8b",
   "metadata": {},
   "source": [
    "#### Optuna's bests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfae9536-b10e-43b5-9c28-3e66dc8d10ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import optuna\n",
    "from performances.utility import retrieve_best_model\n",
    "\n",
    "function_name = \"hadsp\"  # \"desp\" ou \"hadsp\" or \"random\"\n",
    "study = retrieve_best_model(function_name, dataset_name, is_multivariate, variate_type = \"multi\", data_type = \"normal\")\n",
    "\n",
    "# Print the best hyperparameters\n",
    "print(\"Best hyperparameters: \")\n",
    "for param_name, param_value in study.best_trial.params.items():\n",
    "    globals()[param_name] = param_value\n",
    "    print(param_name, param_value)\n",
    "\n",
    "if not is_instances_classification:\n",
    "    use_full_instance = None\n",
    "\n",
    "if function_name in [\"hadsp\", \"desp\"]:\n",
    "    max_increment_span = int(max_window_size) if int(max_window_size) - 100 < 0 else int(max_window_size) - 100\n",
    "    MAX_TIME_INCREMENT = time_increment + time_increment_span #int(max_window_size) or None or TIME_INCREMENT\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e19fd1fc57ed9157",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Function to initialise and generate reservoir\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab39af15",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from reservoir.reservoir import init_matrices\n",
    "from connexion_generation.hag import run_algorithm\n",
    "from scipy import sparse\n",
    "target_rate = 0.7\n",
    "rate_spread = 0.1\n",
    "\n",
    "def initialise_and_run(input_scaling, n, input_connectivity, connectivity, K, bias_scaling, seed, training_set):\n",
    "    \n",
    "    Win, W, bias = init_matrices(n, input_connectivity, connectivity,  K, seed=seed)\n",
    "    bias *= bias_scaling\n",
    "    Win *= input_scaling\n",
    "\n",
    "    if function_name == \"hadsp\":\n",
    "        W, (state_history, delta_z_history, W_history) = run_algorithm(W, Win, bias, leaky_rate, activation_function, training_set, time_increment, weight_increment,\n",
    "                                target_rate, rate_spread, function_name, is_instance=is_instances_classification, use_full_instance = use_full_instance, \n",
    "                                max_increment=MAX_TIME_INCREMENT, max_partners=max_partners, method = \"pearson\", \n",
    "                                n_jobs = 12, visualize=False, record_history=True)\n",
    "\n",
    "    elif function_name == \"desp\":\n",
    "        W, (state_history, delta_z_history, W_history) = run_algorithm(W, Win, bias, leaky_rate, activation_function, training_set, time_increment, weight_increment,\n",
    "                            min_variance, variance_spread, function_name, is_instance=is_instances_classification, use_full_instance = use_full_instance, \n",
    "                            max_increment=MAX_TIME_INCREMENT, max_partners=max_partners, method = \"pearson\", \n",
    "                            intrinsic_saturation=intrinsic_saturation, intrinsic_coef=intrinsic_coef, n_jobs = 12, visualize=False, record_history=True)\n",
    "    else:\n",
    "        raise ValueError(f\"Invalid function_name: {function_name}, should be 'hadsp' or 'desp'\")\n",
    "    connectivity =  np.count_nonzero(W) / (W.shape[0] * W.shape[1])\n",
    "    eigen = sparse.linalg.eigs(W, k=1, which=\"LM\", maxiter=W.shape[0] * 20, tol=0.1, return_eigenvectors=False)\n",
    "    sr = np.max(np.abs(eigen))\n",
    "\n",
    "    \n",
    "    \n",
    "    return Win, W, bias, connectivity, sr, state_history, delta_z_history, W_history\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c4de89f833989d0",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## HAG algorithm generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54472b03-ca39-4f13-92cf-e71e6d63c135",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "(Win_hag_multi, \n",
    " W_hag_multi, \n",
    " bias_hag_multi, \n",
    " connectivity_band, \n",
    " sr_hag_multi, \n",
    " state_history_hag,\n",
    " delta_z_history,\n",
    " W_history) = initialise_and_run(input_scaling, n, input_connectivity, connectivity, K, bias_scaling, SEED, X_pretrain_band)\n",
    "\n",
    "\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "custom_colormap = ListedColormap(np.vstack((plt.cm.cividis(0.0), plt.cm.cividis(np.linspace(0.5, 1, 128)))))\n",
    "vmin = 0\n",
    "vmax = max(np.max(Win_hag_multi), np.max(W_hag_multi))\n",
    "fig, axs = plt.subplots(ncols=3, gridspec_kw=dict(width_ratios=[0.5,6,0.2]))\n",
    "heatmap(Win_hag_multi, cmap=custom_colormap, cbar=False, square=True, ax=axs[0], vmin=vmin, vmax=vmax)\n",
    "heatmap(W_hag_multi, cmap=custom_colormap, yticklabels=False, cbar=False, ax=axs[1], vmin=vmin, vmax=vmax)\n",
    "fig.colorbar(axs[1].collections[0], cax=axs[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c886405-da86-4b5e-858b-fabeff6f5473",
   "metadata": {},
   "source": [
    "## Random matrice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99d44149-8fd2-4528-9cdf-1e384eff50c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoir.reservoir import init_matrices\n",
    "from scipy import sparse\n",
    "\n",
    "# To set particular values\n",
    "function_name = \"random\"  # \"desp\" ou \"hadsp\" or \"random\"\n",
    "study = retrieve_best_model(function_name, dataset_name, is_multivariate, variate_type = \"multi\", data_type = \"normal\")\n",
    "\n",
    "# Print the best hyperparameters\n",
    "print(\"Best hyperparameters: \")\n",
    "random_esn = {}\n",
    "for param_name, param_value in study.best_trial.params.items():\n",
    "    random_esn[param_name] = param_value\n",
    "    print(param_name, param_value)\n",
    "\n",
    "if not is_instances_classification:\n",
    "    use_full_instance = None\n",
    "\n",
    "Win_random_multi, W_random_multi, bias_random_multi = init_matrices(n, 1, random_esn[\"connectivity\"], K, seed=19823)\n",
    "bias_random_multi= bias_random_multi*random_esn['bias_scaling']\n",
    "Win_random_multi= Win_random_multi*random_esn['input_scaling']\n",
    "\n",
    "\n",
    "eigen_random_multi = sparse.linalg.eigs(W_random_multi, k=1, which=\"LM\", maxiter=W_random_multi.shape[0] * 20, tol=0.1, return_eigenvectors=False)\n",
    "sr_random_multi = np.max(np.abs(eigen_random_multi))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21b9a1da-cc49-4b49-aa7c-935d8c12cb2c",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9206b894-d2b6-43a2-8f23-8b62678c44d7",
   "metadata": {},
   "source": [
    "## Algorithm Dynamics\n",
    "\n",
    "### ∆z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22be4078-9172-48f0-b39f-9b9c634c136e",
   "metadata": {},
   "outputs": [],
   "source": [
    "if delta_z_history is not None:\n",
    "    fig, ax = plt.subplots(figsize=(16, 5))\n",
    "    print(f'Dynamics')\n",
    "    \n",
    "    # NEURON ACTIVITY PLOT\n",
    "    random_neurons_indices =  np.sort(np.random.randint(RESERVOIR_SIZE, size=7)) #Size max is 19 because there is not enough colors\n",
    "    random_neurons_indices = np.append(random_neurons_indices, 51)\n",
    "    colors = color_palette(\"tab20\")\n",
    "    # NUMBER_OF_STEP_TO_WATCH\n",
    "    WATCH_FROM = 0\n",
    "    WATCH_TO = 3000\n",
    "    neurons_evolution = np.array(delta_z_history)[WATCH_FROM:WATCH_TO]\n",
    "    x =range(len(delta_z_history))[WATCH_FROM:WATCH_TO]\n",
    "    j = 0\n",
    "    for i in random_neurons_indices:\n",
    "        ax.plot(x, neurons_evolution[: ,i], label=str(i), color=colors[j])\n",
    "        j += 1\n",
    "    fontsize=18\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.tick_params(axis='both', labelsize=fontsize)\n",
    "    plt.xlabel('Time', size=fontsize)\n",
    "    plt.ylabel('Value', size=fontsize)\n",
    "    plt.legend(fontsize=fontsize)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3466f83e-ce7c-4625-93a2-a8630dcbd57c",
   "metadata": {},
   "source": [
    "### Rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eef5f04-6a58-46eb-a958-3937fa36947a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "function_name = \"hadsp\"\n",
    "# NEURON ACTIVITY PLOT \n",
    "random_neurons_indices =  np.sort(np.random.randint(RESERVOIR_SIZE, size=4)) #Size max is 19 because there is not enough colors\n",
    "colors = color_palette(\"tab20\")\n",
    "# NUMBER_OF_STEP_TO_WATCH \n",
    "WATCH_FROM = 0\n",
    "WATCH_TO = 10000\n",
    "neurons_evolution = np.array(state_history_hag)[WATCH_FROM:WATCH_TO]\n",
    "x =range(len(state_history_hag))[WATCH_FROM:WATCH_TO]\n",
    "j = 0\n",
    "for i in random_neurons_indices:\n",
    "    ax.plot(x, neurons_evolution[: ,i], label=str(i), color=colors[j])\n",
    "    j += 1\n",
    "if function_name == \"hadsp\":\n",
    "    ax.plot(x, np.ones(len(x))*(target_rate-rate_spread), color=\"red\", linestyle = \"--\")\n",
    "    ax.plot(x, np.ones(len(x))*(target_rate+rate_spread), color=\"red\", linestyle = \"--\")\n",
    "    ax.text(WATCH_TO, target_rate - rate_spread, r'$\\rho_r - \\beta_r$', color='red', fontsize=fontsize, va='center', ha='left')\n",
    "    ax.text(WATCH_TO, target_rate + rate_spread, r'$\\rho_r + \\beta_r$', color='red', fontsize=fontsize, va='center', ha='left')\n",
    "\n",
    "# target_rate, rate_spread\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Time', size=fontsize)\n",
    "plt.ylabel('Neuron Activity', size=fontsize)\n",
    "plt.legend(title=\"Neurons :\",fontsize=fontsize, title_fontsize=fontsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da0eae33-4a27-49ce-b7b5-1be0dadd9c7e",
   "metadata": {},
   "source": [
    "## Pearson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24170c46-7cc0-4b20-85bf-e29b22589da1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from analysis.richness import pearson\n",
    "num_windows = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e041190-1f47-4252-bfe4-6c55e894ba11",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_correlations_hag, std_correlations_hag = pearson(state_history_hag, num_windows=num_windows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f696a2a0-a60f-4704-b947-605e58b7879f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoir.reservoir import update_reservoir\n",
    "\n",
    "states_history_multi = []\n",
    "neurons_state = np.random.uniform(0, 1, bias_hag_multi.size)\n",
    "inputs = np.concatenate(X_pretrain_band, axis=0) if is_instances_classification else X_pretrain_band\n",
    "for input_value in inputs:\n",
    "    neurons_state = update_reservoir(W_hag_multi, Win_hag_multi, input_value, neurons_state, leaky_rate, bias_hag_multi, activation_function)\n",
    "    states_history_multi.append(neurons_state)\n",
    "\n",
    "mean_correlations_multi, std_correlations_multi = pearson(states_history_multi, num_windows=num_windows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "386591fb-1869-4f7a-8a0b-fd12acee7aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(num_windows)\n",
    "plt.plot(time_windows, mean_correlations_multi, marker='.', linestyle='-', color='b',  label=\"with same input after training\")\n",
    "plt.fill_between(time_windows, np.array(mean_correlations_multi) - np.array(std_correlations_multi),\n",
    "                 np.array(mean_correlations_multi) + np.array(std_correlations_multi), color='blue', alpha=0.2)\n",
    "\n",
    "plt.plot(time_windows, mean_correlations_hag, marker='.', linestyle='-', color='red', label=\"during training\")\n",
    "plt.fill_between(time_windows, np.array(mean_correlations_hag) - np.array(std_correlations_hag),\n",
    "                 np.array(mean_correlations_hag) + np.array(std_correlations_hag), color='red', alpha=0.2)\n",
    "plt.ylim(0, 1)\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Average Correlation', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d42acc5c-4ab9-42a2-bf63-194210f62198",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoir.reservoir import update_reservoir\n",
    "\n",
    "states_history_random = []\n",
    "neurons_state = np.random.uniform(0, 1, bias_random_multi.size)\n",
    "inputs = np.concatenate(X_pretrain_band, axis=0) if is_instances_classification else np.array(X_pretrain_band)\n",
    "for input_value in inputs:\n",
    "    neurons_state = update_reservoir(W_random_multi, Win_random_multi, input_value, neurons_state, leaky_rate, bias_random_multi, activation_function)\n",
    "    states_history_random.append(neurons_state)\n",
    "\n",
    "mean_correlations_random, std_correlations_random = pearson(np.array(states_history_random), num_windows=num_windows)\n",
    "\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(num_windows)\n",
    "plt.plot(time_windows, mean_correlations_random, marker='.', linestyle='-', color='g',  label=\"with same input for random matrice\")\n",
    "plt.fill_between(time_windows, np.array(mean_correlations_random) - np.array(std_correlations_random),\n",
    "                 np.array(mean_correlations_random) + np.array(std_correlations_random), color='g', alpha=0.2)\n",
    "plt.ylim(0, 1)\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Average Correlation', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03592584-13e0-4abb-b06e-b68cc6b3a33f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mean_correlations_random, std_correlations_random = pearson(inputs, num_windows=num_windows)\n",
    "\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(num_windows)\n",
    "plt.plot(time_windows, mean_correlations_random, marker='.', linestyle='-', color='y',  label=\"inputs correlation\")\n",
    "plt.fill_between(time_windows, np.array(mean_correlations_random) - np.array(std_correlations_random),\n",
    "                 np.array(mean_correlations_random) + np.array(std_correlations_random), color='y', alpha=0.2)\n",
    "plt.ylim(0, 1)\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Average Correlation', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42557f32-cb89-45ce-810c-15a35cb5d775",
   "metadata": {},
   "source": [
    "## IPC"
   ]
  },
  {
   "cell_type": "raw",
   "id": "358478ca-5834-43ea-affb-9e94d0a84e48",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-09-09T11:11:47.398677Z",
     "iopub.status.busy": "2024-09-09T11:11:47.398377Z",
     "iopub.status.idle": "2024-09-09T11:12:00.366036Z",
     "shell.execute_reply": "2024-09-09T11:12:00.362075Z",
     "shell.execute_reply.started": "2024-09-09T11:11:47.398655Z"
    },
    "scrolled": true
   },
   "source": [
    "from reservoir.reservoir import update_reservoir\n",
    "import analysis.capacities as CAP\n",
    "\n",
    "MAX_DEG = 8\n",
    "MAX_DEL = 100\n",
    "\n",
    "def cap2vec(capacities,maxdel=MAX_DEL,maxdeg=MAX_DEG):\n",
    "    vec = np.zeros((maxdel,maxdeg))\n",
    "    for idx in range(len(capacities)):\n",
    "        delay=capacities[idx]['delay']\n",
    "        degree=capacities[idx]['degree']\n",
    "        if (delay<=maxdel) and (degree<=maxdeg):\n",
    "            vec[delay-1,degree-1]+=capacities[idx]['score']\n",
    "    return vec\n",
    "\n",
    "Win_capa = Win_hag_multi[Win_hag_multi != 0]\n",
    "Win_capa = Win_capa.reshape(Win_capa.size, -1)\n",
    "\n",
    "\n",
    "V_adsp = []\n",
    "for j, W in enumerate(W_history[0:30:4]):    \n",
    "    warmup_steps = 500\n",
    "    warmups = 2.0*np.random.rand(warmup_steps,1)-1.0\n",
    "\n",
    "    # Important: capacity theory with Legendre polynomials assumes uniform random inputs in [-1,1]\n",
    "    steps = 1000000\n",
    "    inputs = 2.0*np.random.rand(steps, 1)-1.0\n",
    "\n",
    "    state = np.random.uniform(0, 1, n)\n",
    "    for i in range(warmup_steps):\n",
    "        state = update_reservoir(W , Win_capa, warmups[i], state, leaky_rate, bias_hag_multi, activation_function)\n",
    "\n",
    "    capa_history = []\n",
    "    for i in range(steps):\n",
    "        state = update_reservoir(W, Win_capa, inputs[i], state, leaky_rate, bias_hag_multi, activation_function)\n",
    "        capa_history.append(state)\n",
    "    capa_history = np.array(capa_history)\n",
    "\n",
    "    # Measallcapsure capacities on inputs and states after removing a \"warmup period\"\n",
    "    # If you require more extensive printed output of individual capacities:\n",
    "    # set verbose = 1\n",
    "    Citer=CAP.capacity_iterator(maxdeg=MAX_DEG, maxdel=MAX_DEL)#, verbose = 1)\n",
    "    totalcap_adsp,allcaps_adsp,numcaps_adsp,nodes = Citer.collect(inputs,capa_history)\n",
    "\n",
    "    print(\"\\nMeasured \",numcaps_adsp,\" capacities above threshold.\\nTotal capacity = \",totalcap_adsp)\n",
    "    V_adsp.append(cap2vec(allcaps_adsp,maxdel = MAX_DEL, maxdeg = MAX_DEG))\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c2244fa0-fb6b-4dc1-aa43-112e2ac98047",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-09-09T11:12:01.961425Z",
     "iopub.status.busy": "2024-09-09T11:12:01.960646Z",
     "iopub.status.idle": "2024-09-09T11:12:01.997951Z",
     "shell.execute_reply": "2024-09-09T11:12:01.997259Z",
     "shell.execute_reply.started": "2024-09-09T11:12:01.961368Z"
    }
   },
   "source": [
    "\n",
    "# Assuming V_adsp is defined somewhere and MAX_DEG is set\n",
    "V_adsp_array = np.array(V_adsp)  # This should be a 3D array where each slice [i,:,:] corresponds to a 'degree'\n",
    "delrange = np.arange(1, MAX_DEG + 1)\n",
    "\n",
    "# Create a figure for the plot\n",
    "plt.figure(figsize=(10, 6))\n",
    "\n",
    "# This will keep track of the cumulative height of the bars (starting at zero)\n",
    "cumulative_height = np.zeros(8)\n",
    "\n",
    "# Loop through each component and plot\n",
    "for component in range(V_adsp_array.shape[2]):\n",
    "    component_cap = np.sum(V_adsp_array,axis=1)[:,component]\n",
    "    plt.bar(range(0,30,4), component_cap, bottom=cumulative_height, label=f'{component + 1}')\n",
    "    cumulative_height += component_cap\n",
    "\n",
    "plt.xlabel('step')\n",
    "plt.ylabel('Information procssing capacity')\n",
    "plt.title('Information procssing capacity over time')\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "272473cc-1309-4a35-978e-d6197c3a8f31",
   "metadata": {},
   "source": [
    "## Cumulative explained variance\n",
    "Or squared_uncoupled_dynamics\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "013b6cd7-d060-4ba4-add7-13e0a492ff5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import squared_uncoupled_dynamics, squared_uncoupled_dynamics_alternative\n",
    "\n",
    "uds_hag = squared_uncoupled_dynamics_alternative(np.array(state_history_hag), size_window=500, num_windows=num_windows, A=0.99)\n",
    "uds_multi = squared_uncoupled_dynamics_alternative(np.array(states_history_multi), size_window=500, num_windows=num_windows, A=0.99)\n",
    "\n",
    "max_value = np.max([np.max(uds_hag), np.max(uds_hag)])\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(len(uds_hag))\n",
    "plt.plot(range(len(uds_hag)), uds_hag, marker='.', linestyle='-', color='r',  label=\"during training\")\n",
    "plt.plot(range(len(uds_multi)), uds_multi, marker='.', linestyle='-', color='b',  label=\"with same input after training\")\n",
    "plt.ylim(0, max_value)\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('CEV', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8580e57-d01b-4f38-b535-6b47e2850277",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import squared_uncoupled_dynamics\n",
    "\n",
    "uds_random = squared_uncoupled_dynamics_alternative(np.array(states_history_random), size_window=500, num_windows=num_windows, A=0.99)\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "plt.plot(range(len(uds_random)), uds_random, marker='.', linestyle='-', color='g',  label=\"with same input for random matrice\")\n",
    "\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('CEV', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fc60097-34c7-4c5e-9ff9-a4e50bd8cee0",
   "metadata": {},
   "source": [
    "## Linear uncoupled dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08a6a8db-77c0-461e-b65a-319a9fa362b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import linear_uncoupled_dynamics\n",
    "num_windows = 150\n",
    "\n",
    "LUDs_hag = linear_uncoupled_dynamics(np.array(state_history_hag), size_window=500, num_windows=num_windows, theta=0.9)\n",
    "LUDs_multi = linear_uncoupled_dynamics(np.array(states_history_multi), size_window=500, num_windows=num_windows, theta=0.9)\n",
    "\n",
    "max_value = np.max([np.max(LUDs_hag), np.max(LUDs_multi)])\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(len(LUDs_hag))\n",
    "plt.plot(range(len(uds_hag)), LUDs_hag, marker='.', linestyle='-', color='r',  label=\"during training\")\n",
    "plt.plot(range(len(LUDs_multi)), LUDs_multi, marker='.', linestyle='-', color='b',  label=\"with same input after training\")\n",
    "plt.ylim(0, max_value)\n",
    "fontsize=18\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Linear Uncoupled Dynamics', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4086ba68-ffa6-404b-81c4-3c553e92c135",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import linear_uncoupled_dynamics\n",
    "\n",
    "LUDs_random = linear_uncoupled_dynamics(np.array(states_history_random), size_window=500, num_windows=num_windows, theta=0.9)\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(len(LUDs_random))\n",
    "plt.plot(time_windows, LUDs_random, marker='.', linestyle='-', color='g',  label=\"with same input for random matrice\")\n",
    "fontsize=18\n",
    "max_value = np.max([np.max(LUDs_random)])\n",
    "plt.ylim(0, max_value)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Linear Uncoupled Dynamics', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea32b731-f24b-4da3-9273-aff069b3b9e5",
   "metadata": {},
   "source": [
    "## Condition number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74141fa4-b0cf-427e-b534-a0fed799ff27",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import condition_number\n",
    "num_windows = 150\n",
    "\n",
    "CNs_hag = condition_number(np.array(state_history_hag), size_window=500, num_windows=num_windows)\n",
    "CNs_multi = condition_number(np.array(states_history_multi), size_window=500, num_windows=num_windows)\n",
    "\n",
    "max_value = np.max([np.max(CNs_hag), np.max(CNs_hag)])\n",
    "min_value = np.min([np.min(CNs_hag), np.min(CNs_hag)])\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(len(CNs_hag))\n",
    "plt.plot(range(len(uds_hag)), CNs_hag, marker='.', linestyle='-', color='r',  label=\"during training\")\n",
    "plt.plot(range(len(uds_multi)), CNs_multi, marker='.', linestyle='-', color='b',  label=\"with same input after training\")\n",
    "plt.ylim(min_value, max_value)\n",
    "fontsize=18\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Condition Numbers', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e956ac90-411a-4b4e-8a8a-ef80ffb3f075",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.richness\n",
    "reload(analysis.richness)\n",
    "from analysis.richness import condition_number\n",
    "\n",
    "CNs_random = condition_number(np.array(states_history_random), size_window=500, num_windows=num_windows)\n",
    "\n",
    "# Plotting the mean correlations with the standard deviation area\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "\n",
    "time_windows = range(len(CNs_random))\n",
    "plt.plot(time_windows, CNs_random, marker='.', linestyle='-', color='g',  label=\"with same input for random matrice\")\n",
    "fontsize=18\n",
    "max_value = np.max([np.max(CNs_random)])\n",
    "min_value = np.min([np.min(CNs_random)])\n",
    "\n",
    "plt.ylim(min_value, max_value)\n",
    "ax.set_yscale('log')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.tick_params(axis='both', labelsize=fontsize)\n",
    "plt.xlabel('Steps', size=fontsize)\n",
    "plt.ylabel('Linear Uncoupled Dynamics', size=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.legend(fontsize=fontsize)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db9b3d83",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-09T13:36:57.885053Z",
     "start_time": "2023-10-09T13:36:57.882050Z"
    },
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Univariate case"
   ]
  },
  {
   "cell_type": "raw",
   "id": "70e02ca4-fea0-464f-b21b-29d141aa79b8",
   "metadata": {
    "editable": true,
    "execution": {
     "iopub.execute_input": "2024-07-22T12:57:10.271291Z",
     "iopub.status.busy": "2024-07-22T12:57:10.271110Z"
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "if not is_multivariate:\n",
    "    # hag + uni\n",
    "    (Win_hag_uni, \n",
    "     W_hag_uni, \n",
    "     bias_hag_uni, \n",
    "     connectivity_hag_uni, \n",
    "     sr_hag_uni,\n",
    "     state_history_hag_uni,\n",
    "     _,\n",
    "     _) = initialise_and_train(input_scaling, n, input_connectivity, connectivity, n, bias_scaling, SEED, X_pretrain_uni)\n",
    "\n",
    "    from matplotlib.colors import ListedColormap\n",
    "    \n",
    "    custom_colormap = ListedColormap(np.vstack((plt.cm.cividis(0.0), plt.cm.cividis(np.linspace(0.5, 1, 128)))))\n",
    "    vmin = 0\n",
    "    vmax = max(np.max(Win_hag_uni), np.max(W_hag_uni))\n",
    "    fig, axs = plt.subplots(ncols=3, gridspec_kw=dict(width_ratios=[0.5,6,0.2]))\n",
    "    heatmap(Win_hag_uni, cmap=custom_colormap, cbar=False, square=True, ax=axs[0], vmin=vmin, vmax=vmax)\n",
    "    heatmap(W_hag_uni, cmap=custom_colormap, yticklabels=False, cbar=False, ax=axs[1], vmin=vmin, vmax=vmax)\n",
    "    fig.colorbar(axs[1].collections[0], cax=axs[2])"
   ]
  },
  {
   "cell_type": "raw",
   "id": "73fe1a2c-c54e-4ae8-80fe-f57447c8dbe9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-07-22T12:41:26.032884Z",
     "iopub.status.busy": "2024-07-22T12:41:26.032101Z",
     "iopub.status.idle": "2024-07-22T12:43:14.065679Z",
     "shell.execute_reply": "2024-07-22T12:43:14.065383Z",
     "shell.execute_reply.started": "2024-07-22T12:41:26.032821Z"
    }
   },
   "source": [
    "from analysis.richness import pearson\n",
    "\n",
    "pearson(state_history_hag_uni)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "6d205ad0-0c4e-4e76-9a21-d69d1b963c54",
   "metadata": {},
   "source": [
    "if not is_multivariate:    \n",
    "    # random + uni\n",
    "    Win_normal, W_normal, bias_normal = init_matrices(n, 1, connectivity_hag_uni, n, sr_hag_uni)\n",
    "    bias_normal= bias_normal*bias_scaling\n",
    "    Win_normal= Win_normal*input_scaling\n",
    "    \n",
    "    eigen_normal = sparse.linalg.eigs(W_normal, k=1, which=\"LM\", maxiter=W_normal.shape[0] * 20, tol=0.1, return_eigenvectors=False)\n",
    "    sr_normal = np.max(np.abs(eigen_normal))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "084e5dfe-8d4d-4e7e-89af-01fbe9e436d1",
   "metadata": {},
   "source": [
    "print(connectivity_band)\n",
    "if not is_multivariate:\n",
    "    print(connectivity_hag_uni)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5b447f6-e964-4f6c-9116-597cd29c1755",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-09T13:13:26.845564Z",
     "start_time": "2023-10-09T13:13:26.821527Z"
    },
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f4ffd9c-ffc5-47ea-885d-f2f2a3e549fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "from sklearn.metrics import mutual_info_score\n",
    "\n",
    "N_JOBS = -1\n",
    "RIDGE_COEF= 1e-9 if ridge is None else 10**ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89a95dbc-2f36-485e-b1c4-85f6ac6f30ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "from performances.esn_model_evaluation import init_nvar_model, init_reservoir_model, init_ip_reservoir_model\n",
    "\n",
    "reservoir_hag_multi, readout_hag_multi = init_reservoir_model(W_hag_multi, Win_hag_multi, bias_hag_multi, leaky_rate, activation_function, ridge_coef=RIDGE_COEF)\n",
    "reservoir_random_multi, readout_random_multi = init_reservoir_model(W_random_multi, Win_random_multi, bias_random_multi, leaky_rate, activation_function, ridge_coef=RIDGE_COEF)\n",
    "\n",
    "if False: # TODO: change condition for when we want to plot univariate data\n",
    "    reservoir_hag_uni, readout_hag_uni  = init_reservoir_model(W_hag_uni, Win_hag_uni, bias_hag_uni, leaky_rate, activation_function, ridge_coef=RIDGE_COEF)\n",
    "    reservoir_random_uni, readout_random_uni = init_reservoir_model(W_random_uni, Win_random_uni, bias_random_uni, leaky_rate, activation_function, ridge_coef=RIDGE_COEF)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e27380a6-3976-442d-8bb0-34bfc142985a",
   "metadata": {},
   "source": [
    "## Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08537a67-9097-4930-b5af-29659c8cfdc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not is_instances_classification:\n",
    "    raise ValueError(\"This is not the right Classification section.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f3b3a4e-b06e-450a-b7c5-68802be8013d",
   "metadata": {},
   "source": [
    "### Classification for multivariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3534f99-25b5-4653-8dcf-7073183e0eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"X_train_band_duplicated example shape :\", X_train_band[1].shape)     \n",
    "print(\"We should have :\", X_train_band[0].shape[1], \"==\", common_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00b9130e-a834-4256-8f29-e42add2eee19",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from performances.esn_model_evaluation import train_model_for_classification, predict_model_for_classification, compute_score\n",
    "# To remember : \n",
    "#  For reservoirpy   pre_s = W @ r + Win @ (u + noise_gen(dist=dist, shape=u.shape, gain=g_in)) + bias\n",
    "\n",
    "train_data_multi = X_train_band # X_train_band_noisy or X_train_band\n",
    "Y_data = Y_val if use_cross_validation else Y_test\n",
    "\n",
    "mode=\"sequence-to-vector\"\n",
    "train_model_for_classification(reservoir_hag_multi, readout_hag_multi, train_data_multi, Y_train, n_jobs = N_JOBS, mode=mode)\n",
    "train_model_for_classification(reservoir_random_multi, readout_random_multi, train_data_multi, Y_train, N_JOBS, mode=mode)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d7e344e-de91-4709-bfce-467b1d049546",
   "metadata": {},
   "source": [
    "#### noisy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd3dc743-36c8-4821-948d-f14f43af8104",
   "metadata": {},
   "outputs": [],
   "source": [
    "from performances.esn_model_evaluation import predict_model_for_classification, compute_score\n",
    "\n",
    "test_data_multi_noisy = X_test_band_noisy # X_test_band_noisy or X_test_band\n",
    "\n",
    "Y_pred_hag_multi = predict_model_for_classification(reservoir_hag_multi, readout_hag_multi, test_data_multi_noisy, N_JOBS)\n",
    "score = compute_score(Y_pred_hag_multi, Y_data, is_instances_classification, function_name + \" multi\", verbosity=1)\n",
    "\n",
    "Y_pred_random_multi = predict_model_for_classification(reservoir_random_multi, readout_random_multi, test_data_multi_noisy, N_JOBS)\n",
    "score = compute_score(Y_pred_random_multi, Y_data, is_instances_classification, \"random multi\", verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6287d16-2f61-4d2a-97e3-0a0844918879",
   "metadata": {},
   "source": [
    "#### normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49735fb0-8701-4a4f-8d0c-a319c142f674",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "test_data_multi = X_test_band # X_test_band_noisy or X_test_band\n",
    "\n",
    "Y_pred_hag_multi = predict_model_for_classification(reservoir_hag_multi, readout_hag_multi, test_data_multi, -1)\n",
    "score = compute_score(Y_pred_hag_multi, Y_data, is_instances_classification, function_name + \" multi\", verbosity=1)\n",
    "\n",
    "Y_pred_random_multi = predict_model_for_classification(reservoir_random_multi, readout_random_multi, test_data_multi, N_JOBS)\n",
    "score = compute_score(Y_pred_random_multi, Y_data, is_instances_classification, \"random multi\", verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a5f3553-1d28-4340-b4b9-f3dc673bdcac",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-09T13:59:05.097521Z",
     "start_time": "2023-10-09T13:59:04.992489Z"
    }
   },
   "source": [
    "### Classification for univariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db50e2ca-93d7-48bb-bcff-4f885e8d5a6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:  # TODO: change condition for when we want to plot univariate data\n",
    "    # Create a list to store the arrays with the same shape as the expected input of the reservoir\n",
    "\n",
    "    train_data_uni = [ts.reshape(-1, 1) for ts in X_train]\n",
    "    test_data_uni = [ts.reshape(-1, 1) for ts in X_test]\n",
    "\n",
    "    print(\"number of instances in train_data_uni :\", len(train_data_uni), \"should be equal to\", len(X_train))     \n",
    "    print(\"example of train_data_uni train shape :\", train_data_uni[0].shape)     \n",
    "    print(\"We should have :\", train_data_uni[0].shape[1], \"==\", 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b91d01e8-8763-4bf3-ae0a-4c4569231dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "if False: # TODO: change condition for when we want to plot univariate data\n",
    "    reservoir_hag_uni, readout_hag_uni = train_model_for_classification(W_hag_uni, Win_hag_uni, bias_hag_uni, leaky_rate, activation_function, train_data_uni, Y_train, N_JOBS, RIDGE_COEF, mode=\"sequence-to-vector\")\n",
    "\n",
    "    reservoir_random_uni, readout_random_uni = train_model_for_classification(W_normal, Win_normal, bias_normal, leaky_rate, activation_function, train_data_uni, Y_train, N_JOBS, RIDGE_COEF, mode=\"sequence-to-vector\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c5f4028-f457-46b4-a219-7f5fed575af6",
   "metadata": {},
   "outputs": [],
   "source": [
    "if False: # TODO: change condition for when we want to plot univariate data\n",
    "    Y_pred_hag_uni = predict_model_for_classification(reservoir_hag_uni, readout_hag_uni, test_data_uni, N_JOBS)\n",
    "    score = compute_score(Y_pred_hag_uni, Y_test, is_instances_classification, function_name + \" uni\", verbosity=1)\n",
    "    \n",
    "    Y_pred_normal = predict_model_for_classification(reservoir_random_uni, readout_random_uni, test_data_uni, N_JOBS)\n",
    "    score = compute_score(Y_pred_normal, Y_test, is_instances_classification, \"random uni\", verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dc7ab1a-e020-489a-8ff9-9bbd689de15e",
   "metadata": {},
   "source": [
    "## Prediction ahead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f25c06e5-71c1-4f72-ac0f-7e4f1f15afd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "if is_instances_classification:\n",
    "    raise ValueError(\"This is not the right Classification section.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d11c5fc-45b0-4180-9953-bb13800da95b",
   "metadata": {},
   "source": [
    "### Plot datasets\n",
    "Noisy or normal dataset can be ploted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b6d45a7-49c9-4e1e-a473-652bc90d64eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concatenate train and test arrays for plotting\n",
    "combined_data = np.concatenate((X_train_band, X_val_band), axis=0)\n",
    "\n",
    "# noisy version\n",
    "combined_data_noisy = np.concatenate((X_train_band, X_val_band_noisy), axis=0)\n",
    "combined_Y =np.concatenate((Y_train, Y_val), axis=0)\n",
    "\n",
    "# Calculate the merge point index\n",
    "merge_point_index = X_train_band.shape[0]\n",
    "\n",
    "# Define the range around the merge point to plot\n",
    "start_index = merge_point_index - 100\n",
    "end_index = merge_point_index + 100\n",
    "\n",
    "# Plot for a subset N features within a range arround transition from train to test\n",
    "N = 3\n",
    "plt.figure(figsize=(16, 5))\n",
    "for i in [1, 13, 17]: \n",
    "    plt.plot(range(start_index, end_index), combined_data_noisy[start_index:end_index, i], label=f'Feature {i}')\n",
    "plt.plot(range(start_index, end_index), combined_Y[start_index:end_index], label=\"Prediction\")\n",
    "plt.title('Feature Values Around Merge Point')\n",
    "plt.xlabel('Sample Index')\n",
    "plt.ylabel('Feature Value')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f1bcdbf",
   "metadata": {},
   "source": [
    "### Training\n",
    "Noisy or normal dataset can be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d012581-4559-4924-8ba1-9d1b6dea2ad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from performances.esn_model_evaluation import train_model_for_prediction, compute_score\n",
    "from performances.esn_model_evaluation import init_nvar_model, init_reservoir_model, init_ip_reservoir_model\n",
    "\n",
    "if not is_multivariate:\n",
    "    train_data_uni = X_train # X_train_noisy or X_train\n",
    "    \n",
    "    # Training random + MG\n",
    "    esn_random_uni = train_model_for_prediction(reservoir, readout, train_data_uni, Y_train, RIDGE_COEF)\n",
    "    \n",
    "    # Training for hag + MG\n",
    "    esn_hag_uni = train_model_for_prediction(reservoir, readout, train_data_uni, Y_train, RIDGE_COEF)\n",
    "\n",
    "    \n",
    "train_data_multi = X_train_band # X_train_band_noisy or train_band_inputs\n",
    "\n",
    "# Training random + bandfilter\n",
    "esn_random_multi = train_model_for_prediction(reservoir, readout, bias_random_multi, leaky_rate, activation_function, train_data_multi, Y_train, RIDGE_COEF)\n",
    "\n",
    "# Training output HASDP + bandfilter\n",
    "esn_hag_multi = train_model_for_prediction(reservoir, readout, leaky_rate, activation_function, train_data_multi, Y_train, RIDGE_COEF)\n",
    "                                                       "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11fcfe07-1ad0-4d8f-827f-bb9776e8086e",
   "metadata": {},
   "source": [
    "### Prediction\n",
    "Noisy or normal dataset can be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a65dac84",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "if not is_multivariate:\n",
    "    test_data_uni = X_val # X_val_noisy or X_val\n",
    "\n",
    "    # Prediction for random + MG\n",
    "    y_pred_random_uni = esn_random_uni.run(test_data_uni, reset=False) \n",
    "\n",
    "    # Prediction for hag + MG\n",
    "    y_pred_hag_uni = esn_hag_uni.run(test_data_uni, reset=False) \n",
    "\n",
    "\n",
    "test_data_multi = X_val_band # X_test_band_noisy_duplicated or X_test_band_duplicated\n",
    "\n",
    "# Prediction for random + bandfilter\n",
    "y_pred_random_multi = esn_random_multi.run(test_data_multi, reset=False)\n",
    "\n",
    "# Prediction for hag + bandfilter\n",
    "y_pred_hag_multi = esn_hag_multi.run(test_data_multi, reset=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61165568",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from plots.plots import plot_results\n",
    "from performances.esn_model_evaluation import compute_score\n",
    "\n",
    "START_STEP = 30\n",
    "END_STEP = 500\n",
    "slice_range = slice(START_STEP, END_STEP)\n",
    "\n",
    "if not is_multivariate:\n",
    "    print(\"nrmse normal         :\", compute_score(Y_val[slice_range], y_pred_random_uni[slice_range], is_instances_classification))\n",
    "    print(\"nrmse hag          :\", compute_score(Y_val[slice_range], y_pred_hag_uni[slice_range], is_instances_classification))\n",
    "\n",
    "print(\"nrmse random + band  :\", compute_score(Y_val[slice_range], y_pred_random_multi[slice_range], is_instances_classification))\n",
    "print(\"nrmse hag + band   :\", compute_score(Y_val[slice_range], y_pred_hag_multi[slice_range], is_instances_classification))\n",
    "\n",
    "plot_results(y_pred_hag_multi, Y_test, 0, 300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2645bd84-e26a-497c-9053-6f6dd4e69845",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not is_multivariate:\n",
    "    test_data_uni = X_val_noisy # X_val_noisy or X_val\n",
    "\n",
    "    # Prediction for random + MG\n",
    "    y_pred_random_uni = esn_random_uni.run(test_data_uni, reset=False) \n",
    "\n",
    "    # Prediction for hag + MG\n",
    "    y_pred_hag_uni = esn_hag_uni.run(test_data_uni, reset=False) \n",
    "\n",
    "test_data_multi = X_val_band_noisy # X_test_band_noisy_duplicated or X_test_band_duplicated\n",
    "\n",
    "# Prediction for random + bandfilter\n",
    "y_pred_random_multi = esn_random_multi.run(test_data_multi, reset=False)\n",
    "\n",
    "# Prediction for hag + bandfilter\n",
    "y_pred_hag_multi = esn_hag_multi.run(test_data_multi, reset=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26d132a1-0766-497f-bd98-635048b14b44",
   "metadata": {},
   "outputs": [],
   "source": [
    "from plots.plots import plot_results\n",
    "\n",
    "START_STEP = 30\n",
    "END_STEP = 500\n",
    "slice_range = slice(START_STEP, END_STEP)\n",
    "\n",
    "if not is_multivariate:\n",
    "    print(\"nrmse normal         :\", compute_score(y_pred_random_uni[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "    print(\"nrmse hag          :\", compute_score(y_pred_hag_uni[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "print(\"nrmse random + band  :\", compute_score(y_pred_random_multi[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "print(\"nrmse hag + band   :\", compute_score(y_pred_hag_multi[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "\n",
    "plot_results(y_pred_hag_multi, Y_test, 0, 300)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89fb0591-4213-425e-92ff-a5da80e718dc",
   "metadata": {},
   "source": [
    "#### Moving average "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0b5bf2d-368f-4e55-b068-c42c36efe2cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from performances.esn_model_evaluation import compute_score\n",
    "\n",
    "# moving average of the y\n",
    "span=7\n",
    "pad_width = span // 2\n",
    "\n",
    "if not is_multivariate:\n",
    "    ave_y_random_uni = np.convolve(np.pad(y_pred_random_uni.flatten(), pad_width, mode='edge'), np.ones(span), 'valid') / span\n",
    "    ave_y_hag_uni = np.convolve(np.pad(y_pred_hag_uni.flatten(), pad_width, mode='edge') , np.ones(span), 'valid') / span\n",
    "ave_y_random_multi = np.convolve(np.pad(y_pred_random_multi.flatten(), pad_width, mode='edge'), np.ones(span), 'valid') / span\n",
    "ave_y_hag_multi = np.convolve(np.pad(y_pred_hag_multi.flatten(), pad_width, mode='edge'), np.ones(span), 'valid') / span\n",
    "\n",
    "if not is_multivariate:\n",
    "    print(\"nrmse normal         :\", compute_score(ave_y_random_uni[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "    print(\"nrmse hag          :\", compute_score(ave_y_hag_uni[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "print(\"nrmse random + band  :\", compute_score(ave_y_random_multi[slice_range], Y_val[slice_range], is_instances_classification))\n",
    "print(\"nrmse hag + band   :\", compute_score(ave_y_hag_multi[slice_range], Y_val[slice_range], is_instances_classification))\n",
    " \n",
    "plot_results(ave_y_hag_multi.reshape(-1,1), Y_test, 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bc19b97",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nrmse_array_random_uni = []\n",
    "nrmse_array_hag_uni = []\n",
    "nrmse_array_random_multi = []\n",
    "nrmse_array_hag_multi = []\n",
    "\n",
    "for i in range(len(Y_val)-100- step_ahead):\n",
    "    Y_val_i = Y_val[i:100+i]\n",
    "    nrmse_array_random_uni.append(compute_score(Y_val_i, y_pred_random_uni[i:100+i], is_instances_classification))\n",
    "    nrmse_array_hag_uni.append(compute_score(Y_val_i, y_pred_hag_uni[i:100+i], is_instances_classification))\n",
    "    nrmse_array_random_multi.append(compute_score(Y_val_i, y_pred_random_multi[i:100+i], is_instances_classification))\n",
    "    nrmse_array_hag_multi.append(compute_score(Y_val_i, y_pred_hag_multi[i:100+i], is_instances_classification))\n",
    "    \n",
    "log10_nrmse_random_uni= np.log10(nrmse_array_random_uni)\n",
    "log10_nrmse_hag_uni = np.log10(nrmse_array_hag_uni)\n",
    "log10_nrmse_random_multi = np.log10(nrmse_array_random_multi)\n",
    "log10_nrmse_hag_multi = np.log10(nrmse_array_hag_multi)\n",
    "plt.figure()\n",
    "plt.plot(log10_nrmse_random_uni[:1000])\n",
    "plt.plot(log10_nrmse_hag_uni[:1000])\n",
    "plt.plot(log10_nrmse_random_multi[:1000])\n",
    "plt.plot(log10_nrmse_hag_multi[:1000])\n",
    "\n",
    "plt.xlabel('Time steps')\n",
    "plt.ylabel('Log10 NRMSE')\n",
    "plt.legend([\"hag+band\", \"random\", \" random + bandfilter\", \"hag\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b27da37",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cb2478c-bb79-4bc0-855f-b028109c8022",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.cluster.hierarchy import linkage, leaves_list\n",
    "\n",
    "scaler = MinMaxScaler(feature_range=(-1, 1))\n",
    "\n",
    "# pretrain\n",
    "# Be really carefull of the column order here !\n",
    "df_data = scaler.fit_transform(X_pretrain_band.T)\n",
    "df_data = df_data.T\n",
    "df = pd.DataFrame(df_data.T)\n",
    "# Initialize a progress bar for total number of series\n",
    "progress_bar = tqdm(total=df.shape[1]**2, position=0, leave=True)\n",
    "\n",
    "# Initialize an empty correlation matrix\n",
    "correlation_matrix = pd.DataFrame(index=df.columns, columns=df.columns)\n",
    "\n",
    "# Calculate the correlation for each pair of series\n",
    "for col1 in df.columns:\n",
    "    progress_bar.set_description(f\"Processing {col1}\")\n",
    "    for col2 in df.columns:\n",
    "        correlation_matrix.loc[col1, col2] = df[col1].corr(df[col2], method='pearson', min_periods=5)\n",
    "\n",
    "        progress_bar.update(1)  # Update the progress bar after processing each series\n",
    "    \n",
    "progress_bar.close()\n",
    "\n",
    "# Convert correlation_matrix to numeric as it is stored as objects due to tqdm\n",
    "correlation_matrix = correlation_matrix.apply(pd.to_numeric)\n",
    "\n",
    "# Perform hierarchical clustering\n",
    "linked = linkage(correlation_matrix, 'single')\n",
    "\n",
    "# Get the order of rows/columns after hierarchical clustering\n",
    "row_order = leaves_list(linked)\n",
    "\n",
    "# Reorder the correlation matrix\n",
    "sorted_corr_matrix = correlation_matrix.iloc[row_order, :].iloc[:, row_order]\n",
    "\n",
    "# Visualize the sorted correlation matrix with a heatmap\n",
    "plt.figure(figsize=(10,7))\n",
    "sns.heatmap(sorted_corr_matrix, annot=False, cmap='vlag', vmin=-1, vmax=1)\n",
    "plt.title('Clustered Pairwise Correlation of Time Series')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10ec25bb-f3e8-4cdb-87ee-9772460a8cc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import coo_matrix\n",
    "\n",
    "row_order_r = np.array([i + k for i in row_order*K for k in range(K)])\n",
    "\n",
    "# Convert the sparse matrix to a dense format (if memory allows)\n",
    "dense_matrix = W_hag_multi.toarray()\n",
    "\n",
    "# Reorder the dense matrix using the repeated ordering\n",
    "reordered_matrix = dense_matrix[np.ix_(row_order_r, row_order_r)]\n",
    "\n",
    "# Convert the reordered dense matrix back to a sparse format if needed\n",
    "sparse_reordered_matrix = coo_matrix(reordered_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "901730be-d407-4d6d-973d-9fc24951bab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(sparse_reordered_matrix.todense(), cmap=color_palette(\"vlag\", as_cmap=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66e46145-b737-4e62-a8ed-9860be34eb2e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-10-10T12:52:00.989091Z",
     "start_time": "2023-10-10T12:52:00.984122Z"
    }
   },
   "source": [
    "## Motifs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e78de13-7b2d-4ea0-9854-c7ef05c9780e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import analysis.topology\n",
    "reload(analysis.topology)\n",
    "from analysis.topology import motif_distribution, draw_motifs_distribution\n",
    "\n",
    "motifs_count = motif_distribution(W_hag_multi.A)\n",
    "draw_motifs_distribution(motifs_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68f7d48e-892f-44a0-9f33-a0b561e3a1c2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68e4480d-5e4f-4603-90c0-e2f0576e0288",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import poisson, binom\n",
    "\n",
    "def analyze_connectivity_matrix(matrix):\n",
    "    # Extract weights from the matrix (ignoring the diagonal and zeros)\n",
    "    weights = matrix.flatten()\n",
    "    weights = weights[weights != 0]\n",
    "    bin_centers, counts = np.unique(weights, return_counts=True)\n",
    "    \n",
    "    # Calculate the difference for all centers\n",
    "    diffs = np.diff(bin_centers)\n",
    "    # Add the last difference for the last bin\n",
    "    diffs = np.append(diffs, diffs[-1])\n",
    "    \n",
    "    # Calculate the bin edges based on bin centers and differences\n",
    "    bin_edges = bin_centers - diffs/2\n",
    "    # Add the last bin edge\n",
    "    bin_edges = np.append(bin_edges, bin_centers[-1] + diffs[-1]/2)\n",
    "    \n",
    "    # Plot histogram\n",
    "    plt.bar(bin_centers, counts, align='center', alpha=0.6, width=np.diff(bin_centers).min())\n",
    "    \n",
    "    # Fit to Poisson distribution\n",
    "    lambda_est = np.mean(weights)\n",
    "    plt.plot(bin_centers, poisson.pmf(range(len(bin_centers)), lambda_est)*counts[0], 'r-', label='Poisson fit')\n",
    "    \n",
    "    # Fit to Binomial distribution using derived relations\n",
    "    mean = np.mean(weights)\n",
    "    variance = np.var(weights)\n",
    "    \n",
    "    # Calculate p and n estimates\n",
    "    p_est = mean ** 2 / (n * mean - variance) if (n * mean - variance) != 0 else 0\n",
    "    n_est = int(round(mean / p_est)) if p_est != 0 else 0  # n should be integer\n",
    "\n",
    "    # Check parameter validity\n",
    "    if not(0 < p_est < 1):\n",
    "        print(\"Estimated parameters are not valid for the Binomial distribution.\")\n",
    "    else:\n",
    "        x_vals = range(len(bin_centers))\n",
    "        plt.plot(bin_centers, binom.pmf(x_vals, n_est, p_est) * counts[0], 'g-', label='Binomial fit')\n",
    "\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "    return {\"Poisson\": lambda_est, \"Binomial\": (n_est, p_est)}\n",
    "\n",
    "\n",
    "# Assuming W_hag_multi.A is your connectivity matrix\n",
    "analyze_connectivity_matrix(W_hag_multi.A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d995e386-8847-4da0-8197-c73ff68179bb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5ca5733-cd7f-44df-b321-5146b6df7a39",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6298ca7-92fe-4f3b-a3fc-6a94d33ab0f5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "029ae5c6-b0d8-431b-aef0-a7d5606eb4bf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54442bc3-7602-4ab1-b008-deb139afde1f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccced642-5be0-459c-a8b9-d0c1f2de5e1f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d57ae5f-70d7-4250-9044-0e217865a4b5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hag_env",
   "language": "python",
   "name": "hag_env"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}