TS_nonlinear_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import matplotlib\n",
    "from sklearn import preprocessing\n",
    "import os\n",
    "import numpy as np\n",
    "import random as rd\n",
    "import tensorflow as tf \n",
    "\n",
    "from scipy.integrate import odeint\n",
    "import pickle\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from scipy.signal import savgol_filter\n",
    "import scipy\n",
    "import joblib\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.compat.v1.disable_eager_execution()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def toggle_switch(yz, t, beta_y, beta_z, gamma, n_x, n_y, x):\n",
    "    \"\"\"\n",
    "    Right hand side for cascade X -> Y -> Z.  Return dy/dt and dz/dt.\n",
    "    \"\"\"\n",
    "    # Unpack y and z\n",
    "    y, z = yz\n",
    "    \n",
    "    # Compute dy/dt\n",
    "    dy_dt = beta_y * x**n_x / (1 + x**n_x) - y\n",
    "    \n",
    "    # Compute dz/dt\n",
    "    dz_dt = gamma * (beta_z * y**n_y / (1 + y**n_y) - z)\n",
    "    \n",
    "    # Return the result as a NumPy array\n",
    "    return np.array([dy_dt, dz_dt])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def general_toggle_switch(yz, t, *args):\n",
    "    \"\"\"\n",
    "    Right hand side for cascade X -> Y -> Z.  Return dy/dt and dz/dt.\n",
    "    \"\"\"\n",
    "    args = (np.array(args).reshape(2, 6))\n",
    "    args1 = args[0]\n",
    "    args2 = args[1]\n",
    "    # Unpack y and z\n",
    "    y, z = yz\n",
    "    \n",
    "    # Compute dy/dt\n",
    "    dy_dt = generalized_hill_function(yz, args1) - y\n",
    "    \n",
    "    # Compute dz/dt\n",
    "    dz_dt = generalized_hill_function(yz, args2) - z\n",
    "    \n",
    "    # Return the result as a NumPy array\n",
    "    return np.array([dy_dt, dz_dt])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generalized_hill_function(y, args):\n",
    "    a1    = args[0]\n",
    "    a2    = args[1]\n",
    "    b1    = args[2]\n",
    "    b2    = args[3]\n",
    "    n1    = args[4]\n",
    "    n2    = args[5]\n",
    "    x1, x2= y\n",
    "    return (1 + a1*x1**n1 + a2*x2**n2)/(1 + b1*x1**n1 + b2*x2**n2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of time points we want for the solutions\n",
    "n = 400\n",
    "\n",
    "# Time points we want for the solution\n",
    "t = np.linspace(0, 10, n)\n",
    "\n",
    "# Initial condition\n",
    "yz_0 = np.array([0.0, 0.0])\n",
    "\n",
    "# Parameters\n",
    "beta_y = 1.0\n",
    "beta_z = 1.0\n",
    "gamma = 1.0\n",
    "n_x = 2\n",
    "n_y = 2\n",
    "x_0 = 2.0\n",
    "\n",
    "a1_node1_param = 1\n",
    "a2_node1_param = 1\n",
    "b1_node1_param = 1\n",
    "b2_node1_param = 1\n",
    "n1_node1_param = 1\n",
    "n2_node1_param = 1\n",
    "\n",
    "a1_node2_param = 1\n",
    "a2_node2_param = 1\n",
    "b1_node2_param = 1\n",
    "b2_node2_param = 1\n",
    "n1_node2_param = 1 \n",
    "n2_node2_param = 1\n",
    "\n",
    "# Package parameters into a tuple\n",
    "args = (beta_y, beta_z, gamma, n_x, n_y, x_0)\n",
    "\n",
    "# Integrate ODES\n",
    "yz_og = scipy.integrate.odeint(toggle_switch, yz_0, t, args=args).T\n",
    "\n",
    "#y_og, z_og = yz[0], yz[1]\n",
    "\n",
    "#plt.plot(t, y)\n",
    "#plt.plot(t, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /Users/dennisjoshy/anaconda3/lib/python3.7/site-packages/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "If using Keras pass *_constraint arguments to layers.\n"
     ]
    }
   ],
   "source": [
    "#x1=tf.compat.v1.placeholder(tf.compat.v1.double, shape = (len(y)))\n",
    "\n",
    "yz_var=tf.compat.v1.placeholder(tf.compat.v1.double, shape = (2, len(yz_og[0])))\n",
    "\n",
    "#n = tf.compat.v1.Variable(np.random.randn(), name = 'n', dtype = tf.double)\n",
    "#n1_node1 = tf.compat.v1.get_variable('n1',\n",
    "#                   dtype=tf.double,\n",
    "#                   shape=(),\n",
    "#                   initializer=tf.random_uniform_initializer(minval=1., maxval=10.),\n",
    "#                   constraint=lambda z: tf.clip_by_value(z, 1, 10))\n",
    "\n",
    "a1_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "a2_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "b1_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "b2_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "n1_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "n2_node1 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "\n",
    "a1_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "a2_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "b1_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "b2_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "n1_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "n2_node2 = tf.Variable(tf.compat.v1.truncated_normal((1, ), mean=3,stddev=0.3,dtype=tf.double))\n",
    "\n",
    "#n1_node1 = tf.Variable(1.0, dtype = tf.double)\n",
    "\n",
    "#Kx = tf.Variable(tf.compat.v1.truncated_normal((3, 3), mean=0.0,stddev=0.1,dtype=tf.double));\n",
    "#np.abs(Y - W*b)\n",
    "\n",
    "\n",
    "#last_col = tf.constant(np.zeros(shape=(3, 1)), dtype=tf.dtypes.double)\n",
    "\n",
    "#last_col = tf.concat([last_col, [[1.]]], axis=0)\n",
    "\n",
    "#Kx = tf.concat([Kx, last_col], axis=1)  \n",
    "#print(Kx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess = tf.compat.v1.InteractiveSession();\n",
    "sess.run(tf.compat.v1.global_variables_initializer());"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#args = tuple(sess.run([a1_node1, a2_node1, b1_node1, b2_node1, n1_node1, n2_node1, a1_node2, a2_node2, b1_node2, b2_node2, n1_node2, n2_node2]))\n",
    "cost = tf.reduce_sum(tf.pow(yz_og - scipy.integrate.odeint(general_toggle_switch, yz_0, t, args = tuple(sess.run([a1_node1, a2_node1, b1_node1, b2_node1, n1_node1, n2_node1, a1_node2, a2_node2, b1_node2, b2_node2, n1_node2, n2_node2]))).T, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "R2 = (1 - tf.divide(tf.math.reduce_sum(tf.math.square(tf.concat([Xf, [Xf[0]**n1]], axis = 0) - tf.math.reduce_mean(tf.concat([Xf, [Xf[0]**n1]], axis = 0), axis=0))), tf.math.reduce_sum(tf.math.square(tf.concat([Xf, [Xf[0]**n1]], axis = 0) - tf.matmul(Kx, tf.concat([Xf, [Xf[0]**n1]], axis = 0)))))) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "No gradients provided for any variable, check your graph for ops that do not support gradients, between variables [\"<tf.Variable 'Variable:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_1:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_2:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_3:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_4:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_5:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_6:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_7:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_8:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_9:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_10:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_11:0' shape=(1,) dtype=float64>\"] and loss Tensor(\"Sum_2:0\", shape=(), dtype=float64).",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-3aec1c618b61>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mlearning_rate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.00005\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0merror_threshold\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.000000001\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0moptimizer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mAdamOptimizer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlearning_rate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlearning_rate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.9\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta2\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.99\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-08\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0muse_locking\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Adam'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcost\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0minit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mglobal_variables_initializer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/tensorflow_core/python/training/optimizer.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(self, loss, global_step, var_list, gate_gradients, aggregation_method, colocate_gradients_with_ops, name, grad_loss)\u001b[0m\n\u001b[1;32m    408\u001b[0m           \u001b[0;34m\"No gradients provided for any variable, check your graph for ops\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    409\u001b[0m           \u001b[0;34m\" that do not support gradients, between variables %s and loss %s.\"\u001b[0m \u001b[0;34m%\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 410\u001b[0;31m           ([str(v) for _, v in grads_and_vars], loss))\n\u001b[0m\u001b[1;32m    411\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    412\u001b[0m     return self.apply_gradients(grads_and_vars, global_step=global_step,\n",
      "\u001b[0;31mValueError\u001b[0m: No gradients provided for any variable, check your graph for ops that do not support gradients, between variables [\"<tf.Variable 'Variable:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_1:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_2:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_3:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_4:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_5:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_6:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_7:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_8:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_9:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_10:0' shape=(1,) dtype=float64>\", \"<tf.Variable 'Variable_11:0' shape=(1,) dtype=float64>\"] and loss Tensor(\"Sum_2:0\", shape=(), dtype=float64)."
     ]
    }
   ],
   "source": [
    "training_epochs = 85000\n",
    "learning_rate = 0.00005\n",
    "error_threshold = 0.000000001\n",
    "optimizer = tf.compat.v1.train.AdamOptimizer(learning_rate=learning_rate, beta1=0.9, beta2=0.99, epsilon=1e-08, use_locking=False, name='Adam').minimize(cost)\n",
    "init = tf.compat.v1.global_variables_initializer()\n",
    "c = 100\n",
    "epoch = 0\n",
    "\n",
    "with tf.compat.v1.Session() as sesh:    \n",
    "    sesh.run(init)    \n",
    "    print(\"Initial n1\", sesh.run(n1_node1))\n",
    "    while epoch < training_epochs and c > error_threshold:\n",
    "        args = tuple(sess.run([a1_node1, a2_node1, b1_node1, b2_node1, n1_node1, n2_node1, a1_node2, a2_node2, b1_node2, b2_node2, n1_node2, n2_node2]))\n",
    "        yz = scipy.integrate.odeint(general_toggle_switch, yz_0, t, args = args).T\n",
    "        c = sesh.run(cost, feed_dict = {yz_var: yz})\n",
    "        print(c)\n",
    "        if epoch % 5000 == 0:\n",
    "            print(\"Epoch:\", epoch, \"{:.5f}\".format(c))\n",
    "            print(\"Exponent\", sesh.run(n1))\n",
    "            #print(\"R2\", sesh.run(R2, feed_dict = {Xp: np.array(Xp_data).T, Xf: np.array(Xf_data).T}))\n",
    "        sesh.run(optimizer, feed_dict = {yz_var: yz})\n",
    "        epoch+=1\n",
    "    #KxT_num = sesh.run(Kx)\n",
    "    #sesh.close()\n",
    "    #print(sesh.run(n1))\n",
    "    print(sesh.run(cost, feed_dict = {yz_var: yz}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_name = 'CFS_data'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler().fit(Xp)\n",
    "Xp_ref_scaled = np.array(Xp)#scaler.transform(Xp)\n",
    "Xf_ref_scaled = np.array(Xf)#scaler.transform(Xf)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,10))\n",
    "\n",
    "ax1.scatter(np.linspace(0, len(Xp_ref_scaled[:, 0]), len(Xp_ref_scaled[:, 0])), Xp_ref_scaled[:, 0], color = 'r')\n",
    "ax1.scatter(np.linspace(1, len(Xf_ref_scaled[:, 0])+1, len(Xf_ref_scaled[:, 0])), Xf_ref_scaled[:, 0], color = 'b')\n",
    "\n",
    "ax1.set_xlim([-2, 22])\n",
    "ax1.set_ylim([-5, 5])\n",
    "\n",
    "ax2.scatter(np.linspace(0, len(Xp_ref_scaled[:, 1]), len(Xp_ref_scaled[:, 1])), Xp_ref_scaled[:, 1], color = 'r')\n",
    "ax2.scatter(np.linspace(1, len(Xf_ref_scaled[:, 1])+1, len(Xf_ref_scaled[:, 1])), Xf_ref_scaled[:, 1], color = 'b')\n",
    "\n",
    "#ax2.set_xlim([-2, 22])\n",
    "#ax2.set_ylim([-5, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scaler_filename = \"CFS_scaler.save\"\n",
    "joblib.dump(scaler, scaler_filename) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PsiXp = np.vstack([Xp_ref_scaled.T, Xp_ref_scaled.T[0]**2.0])\n",
    "PsiXf = np.vstack([Xf_ref_scaled.T, Xf_ref_scaled.T[0]**2.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kx = np.matmul(PsiXf, np.matmul(PsiXp.T, np.linalg.inv(np.matmul(PsiXp, PsiXp.T))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum((PsiXf - np.matmul(Kx, PsiXp))**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import joblib\n",
    "\n",
    "scaler_filename = 'CFS_scaler.save'\n",
    "scaler = joblib.load(scaler_filename) \n",
    "\n",
    "P = np.diag(scaler.scale_)\n",
    "b = scaler.mean_\n",
    "\n",
    "K11 = np.array([[0.7, 0], [0.8, 0.9]])\n",
    "K12 = np.array([[0], [0.6]])\n",
    "K21 = np.array([[0, 0]])\n",
    "K22 = np.array([[0.49]])\n",
    "\n",
    "Ks11 = np.matmul(np.matmul(P, K11), np.linalg.inv(P))\n",
    "Ks12 = np.matmul(P, K12)\n",
    "Ks13 = np.matmul(np.eye(2) - Ks11, b)\n",
    "Ks21 = np.matmul(K21, np.linalg.inv(P))\n",
    "Ks22 = K22\n",
    "Ks23 = np.matmul(Ks21, b)\n",
    "\n",
    "Ks1 = np.concatenate([Ks11, Ks12, Ks13.reshape(-1, 1)], axis = 1)\n",
    "Ks2 = np.concatenate([Ks21, Ks22, Ks23.reshape(-1, 1)], axis = 1)\n",
    "Ks = np.concatenate([Ks1, Ks2, np.array([[0, 0, 0, 1]])], axis = 0)\n",
    "\n",
    "Ks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xp_final_learned = []\n",
    "Xf_final_learned = []\n",
    "Xp_final_actual = []\n",
    "Xf_final_actual = []\n",
    "#t = np.array([i for i in range(0, N+1)])\n",
    "x_learned = np.zeros(2)\n",
    "x_actual = np.zeros(2)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,10))\n",
    "#fig.suptitle('Horizontally stacked subplots')\n",
    "T = 100\n",
    "t = np.linspace(0, T, T+1)\n",
    "\n",
    "for ic in ICs:\n",
    "    #ic_s = scaler.transform(np.array([ic]))\n",
    "    x_learned = np.array(ic)\n",
    "    x_actual[0] = ic[0]\n",
    "    x_actual[1] = ic[1]\n",
    "    X_learned = []\n",
    "    #print(\"x_actual\", x_actual)\n",
    "    #print(\"x_learned\", x_learned)\n",
    "    X_actual = []\n",
    "    X_learned.append([x_learned[0], x_learned[1]]) ## ICs\n",
    "    X_actual.append(np.array([x_actual[0], x_actual[1]])) \n",
    "    for k in range(0, T):\n",
    "        y_learned = np.matmul(KxT_num, np.vstack([x_learned[0], x_learned[1], x_learned[0]**2]))[0:2].T[0]\n",
    "        x_learned = y_learned#np.array([[y_learned[0], y_learned[1]]])\n",
    "        X_learned.append([x_learned[0], x_learned[1]]) \n",
    "        y_actual = model(x_actual)\n",
    "        x_actual = y_actual\n",
    "        X_actual.append(x_actual)\n",
    "    ax1.scatter(t, np.array(X_actual).T[0])\n",
    "    ax1.plot(t, np.array(X_learned).T[0])\n",
    "    ax2.scatter(t, np.array(X_actual).T[1])\n",
    "    ax2.plot(t, np.array(X_learned).T[1])"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}