Example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoir import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load sample data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('example_data/MackeyGlass_t17.txt')\n",
    "\n",
    "train = data[:4000].reshape(-1, 1)\n",
    "test = data[4000:6000].reshape(-1, 1)\n",
    "plt.title(\"Train data\")\n",
    "plt.plot(train)\n",
    "plt.show()\n",
    "plt.title(\"Zoom in on first 400 train points\")\n",
    "plt.plot(train[:400])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define noisy data for Bayesian Optimization\n",
    "snr = 4\n",
    "noise_std = np.sqrt(train.std(ddof=1)**2 / snr)\n",
    "\n",
    "print('Standard deviation of noise is:', noise_std)\n",
    "noisy_train = train + np.random.normal(0, noise_std, size=train.shape)\n",
    "noisy_test = test + np.random.normal(0, noise_std, size=test.shape)\n",
    "plt.plot(noisy_train[:400])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combine multiple series for evaluation if desired"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiple_train = np.hstack((train, noisy_train))\n",
    "multiple_test = np.hstack((test, noisy_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = {\n",
    "    'input_scaling': (0, 1),\n",
    "    'feedback_scaling': (0, 1),\n",
    "    'leaking_rate' : (0, 1), \n",
    "    'spectral_radius': (0, 1.25),\n",
    "    'regularization': (-12, 1),\n",
    "    'connectivity': (-3, 0),\n",
    "    'n_nodes':  (100, 1000)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set optimization parameters\n",
    "esn_cv = EchoStateNetworkCV(bounds=bounds,\n",
    "                            initial_samples=100,\n",
    "                            subsequence_length=250,\n",
    "                            eps=1e-3,\n",
    "                            cv_samples=3, \n",
    "                            max_iterations=1000, \n",
    "                            scoring_method='tanh',\n",
    "                            verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Optimize (this may take a while!)\n",
    "best_arguments = esn_cv.optimize(y=multiple_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train best model(s)\n",
    "esn = EchoStateNetwork(**best_arguments)\n",
    "esn.train(y=train)\n",
    "\n",
    "esn2 = EchoStateNetwork(**best_arguments)\n",
    "esn2.train(y=noisy_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test\n",
    "score = esn.test(y=test, scoring_method='rmse')\n",
    "score2 = esn.test(y=noisy_test, scoring_method='rmse')\n",
    "print(score)\n",
    "print(score2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diagnostic plots\n",
    "plt.plot(esn.predict(100), label='Predicted')\n",
    "plt.plot(test[:100], label='Ground truth')\n",
    "plt.title('Prediction on next 100 steps')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.plot(esn2.predict(100), label='Predicted')\n",
    "plt.plot(noisy_test[:100], label='Ground truth')\n",
    "plt.title('Prediction on next 100 steps')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manually tuned example, for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regular approach\n",
    "esn = EchoStateNetwork(n_nodes=1000, connectivity=0.01, input_scaling=0.5, feedback_scaling=0.5, leaking_rate=0.3, \n",
    "                       spectral_radius=1.25, regularization=1e-8, feedback=True)\n",
    "complete_data, y_used, burn_in = esn.train(y=train, burn_in=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prediction sample for ground truth\n",
    "plt.plot(test[:500], label='Test')\n",
    "plt.plot(esn.predict(500), label='Predicted')\n",
    "plt.legend()\n",
    "esn.test(y=test[:500], scoring_method='nrmse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stepwise prediction\n",
    "five_step_prediction = esn.predict_stepwise(test[:20], steps_ahead=5)\n",
    "plt.plot(test[:20], label='Test')\n",
    "plt.plot(five_step_prediction, label='Predicted')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "five_step_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}