Diagonal conceptor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.esn import ESN\n",
    "from src.utils import config, helper\n",
    "from src.conceptors import compute_conceptor, loading_ridge_report, compute_conceptor_diag, ridge_regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diagonal conceptor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prng = np.random.default_rng(1001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "esnConfig = config.ESNConfig(\n",
    "    input_size=1,\n",
    "    reservoir_size= N,\n",
    "    output_size=1,\n",
    "    feedback=False,\n",
    "    spectral_radius=1,\n",
    "    init_weights='normal',\n",
    "    init_weights_b='normal',\n",
    "    init_weights_in='normal',\n",
    "    init_weights_density=10/N,\n",
    "    init_weights_in_density=1.0,\n",
    "    init_weights__args={'loc': 0., 'scale': 1.},\n",
    "    init_weights_b__args={'loc': 0., 'scale': 0.2},\n",
    "    init_weights_in__args={'loc': 0., 'scale': 1.},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "esn = ESN(esnConfig, prng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_pattern = 5000\n",
    "ut = [\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(27)),\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(42)),\n",
    "    helper.n_periodic_jong(0, 5, T_pattern),\n",
    "    helper.n_periodic_jong(1, 5, T_pattern),\n",
    "    helper.n_sine(8.83, T_pattern, phase=0),\n",
    "    helper.n_sine(9.83, T_pattern, phase=0.)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_WASHOUT = 200\n",
    "T_ADAPT = 500\n",
    "T_LEARN = 5000 - 200 - 500\n",
    "LOADING_REGULARIZER = 1e-3\n",
    "APERTURE = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diagonal \"trick\": initial biais of the input driven reservoir with random conceptor\n",
    "C_diag = []\n",
    "for i in range (4):\n",
    "    diag = prng.uniform(size=(N))\n",
    "    C_diag.append(np.diag(diag))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptation (change C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_init = prng.uniform(low=-1., high=1., size=(esn.reservoir_size, 1))\n",
    "xt_conceptor, yt_conceptor = list(zip(*[\n",
    "    esn.harvest_states(ut[i][:T_WASHOUT+T_ADAPT].copy(), x_init=x_init.copy(), C=C_diag[i])\n",
    "    for i in range(len(ut))\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recompute the conceptor, but restrict it to a diagonale matrix (cf. eq 38 Jong (2021))\n",
    "Ci = [\n",
    "    \n",
    "    compute_conceptor_diag(xt_i[T_WASHOUT:, :], aperture=APERTURE)\n",
    "    for xt_i in xt_conceptor\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monitor the scaling values of C\n",
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(np.sort(Ci[i][:,0]))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning (Wout and W with C_diag on the loop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt_conceptor, yt_conceptor = list(zip(*[\n",
    "    esn.harvest_states(ut[i][T_WASHOUT+T_ADAPT:].copy(), x_init=xt_conceptor[i][[-1],:].T.copy(), C=Ci[i])\n",
    "    for i in range(len(ut))\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ut_learn = [ut[i][T_WASHOUT+T_ADAPT:] for i in range (len(ut))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = helper.concatenate_patterns(xt_conceptor, 1)\n",
    "U = helper.concatenate_patterns(ut_learn, 1)\n",
    "Y_T = U.copy() \n",
    "\n",
    "X_ = helper.concatenate_patterns(xt_conceptor, 1, shift=-1)\n",
    "B = np.repeat(esn.b, X_.shape[0], axis=1).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#W_loaded = loading_ridge_report(X, X_, B, regularizer=LOADING_REGULARIZER), good exercise to do it like that also! but more complex\n",
    "\n",
    "W_loaded = ridge_regression(X_, (np.dot(esn.w, X_.T) + np.dot(esn.w_in, U.T)).T, regularizer=LOADING_REGULARIZER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_before_loading = esn.w.copy()\n",
    "esn.w = W_loaded.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# W_out is not yet learned, the output is not reflecting the input (cf. scale)\n",
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(ut[i][:50], ls='--')\n",
    "    ax[i].plot(yt_conceptor[i][:50])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(ut[i][:30], ls='--')\n",
    "    ax[i].plot(xt_conceptor[i][:30][:,1:10])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "esn.update_weights(X, U, Y_T, alpha=1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ut_zero = np.zeros_like(ut[0])\n",
    "x_init = prng.uniform(low=-1., high=1., size=(esn.reservoir_size, 1))\n",
    "xt_trained, yt_trained = list(zip(*[\n",
    "    esn.harvest_states(ut_zero.copy(), x_init=x_init.copy(), C=Ci[i])\n",
    "    for i in range(len(ut))\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(ut[i][50:90], ls='--')\n",
    "    ax[i].plot(yt_trained[i][50:90])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(ut[i][:50], ls='--')\n",
    "    ax[i].plot(xt_trained[i][:50][:,10:20])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# full_procedure (in one step using trainingConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prng = np.random.default_rng(1001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 400\n",
    "T_WASHOUT = 200\n",
    "T_ADAPT = 500\n",
    "T_LEARN = 5000 - 200 - 500\n",
    "LOADING_REGULARIZER = 1e-3\n",
    "APERTURE = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trainingConfig = config.TrainingConfig(\n",
    "    washout=T_WASHOUT,\n",
    "    aperture=APERTURE,\n",
    "    adapt = T_ADAPT,\n",
    "    w_regularizer=LOADING_REGULARIZER,\n",
    "    wout_regularizer=1e-3\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "esnConfig = config.ESNConfig(\n",
    "    input_size=1,\n",
    "    reservoir_size= N,\n",
    "    output_size=1,\n",
    "    feedback=False,\n",
    "    spectral_radius=1,\n",
    "    init_weights='normal',\n",
    "    init_weights_b='normal',\n",
    "    init_weights_in='normal',\n",
    "    init_weights_density=10/N,\n",
    "    init_weights_in_density=1.0,\n",
    "    init_weights__args={'loc': 0., 'scale': 1.},\n",
    "    init_weights_b__args={'loc': 0., 'scale': 0.2},\n",
    "    init_weights_in__args={'loc': 0., 'scale': 1.},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_pattern = 5000\n",
    "ut = [\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(27)),\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(42)),\n",
    "    helper.n_periodic_jong(0, 5, T_pattern),\n",
    "    helper.n_periodic_jong(1, 5, T_pattern),\n",
    "    helper.n_sine(8.83, T_pattern, phase=0),\n",
    "    helper.n_sine(9.83, T_pattern, phase=0.)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "esn = ESN(esnConfig, prng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ci = esn.full_procedure_diag(ut, trainingConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ut_zero = np.zeros_like(ut[0])\n",
    "x_init = prng.uniform(low=-1., high=1., size=(esn.reservoir_size, 1))\n",
    "xt_trained, yt_trained = list(zip(*[\n",
    "    esn.harvest_states(ut_zero.copy(), x_init=x_init.copy(), C=Ci[i])\n",
    "    for i in range(len(ut))\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "for i in range(len(ut)):\n",
    "    ax[i].plot(ut[i][50:90], ls='--')\n",
    "    ax[i].plot(yt_trained[i][50:90])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hyperparameters search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_pattern = 5000\n",
    "ut = [\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(27)),\n",
    "    #helper.n_periodic(5, T_pattern, np.random.default_rng(42)),\n",
    "    helper.n_periodic_jong(0, 5, T_pattern),\n",
    "    helper.n_periodic_jong(1, 5, T_pattern),\n",
    "    helper.n_sine(8.83, T_pattern, phase=0),\n",
    "    helper.n_sine(9.83, T_pattern, phase=0.)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "Winscal = [1.2,1.4,1]\n",
    "SpecRadius = [0.6,1,1.4,2]\n",
    "Apertj = [1,6,8,10,14]\n",
    "\n",
    "b_scal = [0.1,0.2,0.4]\n",
    "Wout_reg = [0,1e-1,1e-2,1e-3]\n",
    "W_reg = [0,1e-1,1e-2,1e-3]\n",
    "\n",
    "\n",
    "T_WASHOUT = 200\n",
    "T_ADAPT = 500\n",
    "T_LEARN = 5000 - 200 - 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_RMSE = []\n",
    "def hyperparasearch (Winscal, b_scal, Wout_reg, SpecRadius, Apertj, N, ut, W_reg, T_WASHOUT, T_ADAPT, T_LEARN):\n",
    "    for b in b_scal:\n",
    "        for win in Winscal:\n",
    "            for rad in SpecRadius:\n",
    "                for wout in Wout_reg:\n",
    "                    for w in W_reg:\n",
    "                        for aj in Apertj: \n",
    "                            prng = np.random.default_rng(1001)\n",
    "                            print(\"aj: \"+str(aj)+\", b: \"+str(b)+\",wout:  \"+str(wout)+\", rad: \"+str(rad), \", w:\" + str(w) +\n",
    "                                   \",win : \"+str(win))\n",
    "                            # Initialization\n",
    "                            esnConfig = config.ESNConfig(\n",
    "                                input_size=1,\n",
    "                                reservoir_size= N,\n",
    "                                output_size=1,\n",
    "                                feedback=False,\n",
    "                                spectral_radius=rad,\n",
    "                                init_weights='normal',\n",
    "                                init_weights_b='normal',\n",
    "                                init_weights_in='normal',\n",
    "                                init_weights_density=10/N,\n",
    "                                init_weights_in_density=1.0,\n",
    "                                init_weights__args={'loc': 0., 'scale': 1.},\n",
    "                                init_weights_b__args={'loc': 0., 'scale': b},\n",
    "                                init_weights_in__args={'loc': 0., 'scale': win},\n",
    "                            )\n",
    "\n",
    "                            trainingConfig = config.TrainingConfig(\n",
    "                                washout=T_WASHOUT,\n",
    "                                aperture=aj,\n",
    "                                adapt = T_ADAPT,\n",
    "                                w_regularizer=w,\n",
    "                                wout_regularizer=wout\n",
    "                            )\n",
    "\n",
    "                            esn = ESN(esnConfig, prng)\n",
    "                            Ci = esn.full_procedure(ut, trainingConfig)\n",
    "\n",
    "                            ut_zero = np.zeros_like(ut[0])\n",
    "                            x_init = prng.uniform(low=-1., high=1., size=(esn.reservoir_size, 1))\n",
    "                            xt_trained, yt_trained = list(zip(*[\n",
    "                                esn.harvest_states(ut_zero.copy(), x_init=x_init.copy(), C=Ci[i])\n",
    "                                for i in range(len(ut))\n",
    "                            ]))\n",
    "\n",
    "\n",
    "                            fig, ax = plt.subplots(len(ut), 1, figsize=(12, len(ut)), sharex=True, sharey=True)\n",
    "                            for i in range(len(ut)):\n",
    "                                ax[i].plot(ut[i][50:90], ls='--')\n",
    "                                ax[i].plot(yt_trained[i][50:90])\n",
    "                            plt.tight_layout()\n",
    "                            plt.show()\n",
    "                            RMSE, _ = helper.testLRMSE(ut, yt_trained, 10)\n",
    "                            L_RMSE.append([min(RMSE),\"aj: \"+str(aj)+\", b: \"+str(b)+\",wout:  \"+str(wout)+\", rad: \"+str(rad), \", w:\" + str(w) +\n",
    "                                   \",win : \"+str(win)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "hyperparasearch (Winscal, b_scal, Wout_reg, SpecRadius, Apertj, N, ut,W_reg, T_WASHOUT, T_ADAPT, T_LEARN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "fac611bf71fb0da9434ed307100dbfd751275a58e894d436e9c2da9c708666f8"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.8.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "245.508px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}