IPython notebook code base

.gitignore

@@ -2,9 +2,46 @@
 *.zip
 *.rtf
 *.eps
-Python_Codebase_Notes.txt
-UTILS\*
+*.pyc
+*.mdl
+*.slx
 *.m
-*\.ipynb_checkpoints\*
-*\__pychache__
-*\.DS_Store
+*.doc#
+
+Python_Codebase_Notes.txt
+UTILS/*
+auxiliary book code files/*
+*/.ipynb_checkpoints/*
+*/*/.ipynb_checkpoints/*
+*/__pycache__/*
+*/*/__pycache__/*
+*/.DS_Store
+*/*/.DS_Store
+
+CH01/dog.jpg
+CH01/None0000000.png
+CH01/temp.mat
+CH02/CH02_SEC04_2_SpectrogramBeethoven/*
+CH02/None0000000.png
+CH03/CC2.mat
+CH03/jelly.jpg
+CH03/mustache.jpg
+CH05/catData.mat
+CH05/catData_w.mat
+CH05/dogData.mat
+CH05/dogData_w.mat
+CH05/census1994.csv
+CH06/catData_w.mat
+CH06/dogData_w.mat
+CH06/CH06_SEC01_1_NN/*
+CH06/CH06_SEC06_1_NNLorenz/*
+CH06/extra/*
+CH08/CH08_SEC07_3_LQG_Simulink/*
+CH09/*.mat
+CH09/CH09_SEC02_2_BalancedTruncation/*
+CH09/CH09_z_extra_BT_ERA_OKID/*
+CH10/CH10_SEC02_GA_PID/*
+CH10/CH10_SEC03_ESC_Krstic_ex1p3/*
+CH10/CH10_SEC03_ESC_Krstic_ex1p3/NEWkrstic_example_1p3_wnoise.slx
+CH12/old_extra/*
+

CH04/CH04_SEC07_2_RegressAIC_BIC.ipynb

@@ -0,0 +1,94 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import rcParams\n",
+    "from statsmodels.tsa import arima_process, arima_model\n",
+    "# Using the StatsModels module available at\n",
+    "# https://www.statsmodels.org/dev/install.html\n",
+    "\n",
+    "\n",
+    "rcParams.update({'font.size': 18})\n",
+    "plt.rcParams['figure.figsize'] = [12, 12]\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "AIC: [175.44759775 175.44759775   0.        ]\n",
+      "BIC: [180.65793813 180.65793813   0.        ]\n"
+     ]
+    }
+   ],
+   "source": [
+    "np.random.seed(123) # for random data reproducibility\n",
+    "T = 100 # Sample size\n",
+    "# DGP = sm.ARIMA(x,order=(0,0,0))\n",
+    "# arparams = np.array([.2, 0.5])\n",
+    "# maparams = np.array([-4])\n",
+    "# arparams = np.r_[1, arparams]\n",
+    "\n",
+    "arparams = np.array([-4, .2, 0.5])\n",
+    "maparams = np.array([1])\n",
+    "\n",
+    "\n",
+    "\n",
+    "arma_process = sm.tsa.arima_process.ArmaProcess(arparams, maparams)\n",
+    "y = arma_process.generate_sample(T,scale=2)\n",
+    "\n",
+    "logL = np.zeros(3) # log likelihood vector\n",
+    "aic = np.zeros(3)  # AIC vector\n",
+    "bic = np.zeros(3)  # BIC vector\n",
+    "\n",
+    "for j in range(2):\n",
+    "    model_res = sm.tsa.arima_model.ARMA(y, (0,0)).fit(trend='c', disp=0,start_ar_lags=j+1,method='mle')\n",
+    "    logL[j] = model_res.llf\n",
+    "    aic[j] = model_res.aic\n",
+    "    bic[j] = model_res.bic\n",
+    "\n",
+    "print('AIC: {:}'.format(aic))\n",
+    "print('BIC: {:}'.format(bic))"
+   ]
+  },
+  {
+   "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": 2
+}

CH07/CH07_SEC03_SINDY_Lorenz.ipynb

@@ -0,0 +1,175 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import rcParams\n",
+    "from mpl_toolkits.mplot3d import Axes3D\n",
+    "from scipy import integrate\n",
+    "\n",
+    "\n",
+    "rcParams.update({'font.size': 18})\n",
+    "plt.rcParams['figure.figsize'] = [12, 12]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Simulate the Lorenz System\n",
+    "\n",
+    "dt = 0.01\n",
+    "T = 50\n",
+    "t = np.arange(dt,T+dt,dt)\n",
+    "beta = 8/3\n",
+    "sigma = 10\n",
+    "rho = 28\n",
+    "n = 3\n",
+    "\n",
+    "def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):\n",
+    "    x, y, z = x_y_z\n",
+    "    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]\n",
+    "\n",
+    "np.random.seed(123)\n",
+    "x0 = (-8,8,27)\n",
+    "\n",
+    "x = integrate.odeint(lorenz_deriv, x0, t,rtol=10**(-12),atol=10**(-12)*np.ones_like(x0))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Compute Derivative\n",
+    "dx = np.zeros_like(x)\n",
+    "for j in range(len(t)):\n",
+    "    dx[j,:] = lorenz_deriv(x[j,:],0,sigma,beta,rho)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 68,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## SINDy Function Definitions\n",
+    "\n",
+    "def poolData(yin,nVars,polyorder):\n",
+    "    n = yin.shape[0]\n",
+    "    yout = np.zeros((n,1))\n",
+    "    \n",
+    "    # poly order 0\n",
+    "    yout[:,0] = np.ones(n)\n",
+    "    \n",
+    "    # poly order 1\n",
+    "    for i in range(nVars):\n",
+    "        yout = np.append(yout,yin[:,i].reshape((yin.shape[0],1)),axis=1)\n",
+    "    \n",
+    "    # poly order 2\n",
+    "    if polyorder >= 2:\n",
+    "        for i in range(nVars):\n",
+    "            for j in range(i,nVars):\n",
+    "                yout = np.append(yout,(yin[:,i]*yin[:,j]).reshape((yin.shape[0],1)),axis=1)\n",
+    "                \n",
+    "    # poly order 3\n",
+    "    if polyorder >= 3:\n",
+    "        for i in range(nVars):\n",
+    "            for j in range(i,nVars):\n",
+    "                for k in range(j,nVars):\n",
+    "                    yout = np.append(yout,(yin[:,i]*yin[:,j]*yin[:,k]).reshape((yin.shape[0],1)),axis=1)\n",
+    "    \n",
+    "    return yout\n",
+    "\n",
+    "def sparsifyDynamics(Theta,dXdt,lamb,n):\n",
+    "    Xi = np.linalg.lstsq(Theta,dXdt,rcond=None)[0] # Initial guess: Least-squares\n",
+    "    \n",
+    "    for k in range(10):\n",
+    "        smallinds = np.abs(Xi) < lamb # Find small coefficients\n",
+    "        Xi[smallinds] = 0                          # and threshold\n",
+    "        for ind in range(n):                       # n is state dimension\n",
+    "            biginds = smallinds[:,ind] == 0\n",
+    "            # Regress dynamics onto remaining terms to find sparse Xi\n",
+    "            Xi[biginds,ind] = np.linalg.lstsq(Theta[:,biginds],dXdt[:,ind],rcond=None)[0]\n",
+    "            \n",
+    "    return Xi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 71,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[  0.           0.           0.        ]\n",
+      " [-10.          28.           0.        ]\n",
+      " [ 10.          -1.           0.        ]\n",
+      " [  0.           0.          -2.66666667]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           1.        ]\n",
+      " [  0.          -1.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]\n",
+      " [  0.           0.           0.        ]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "Theta = poolData(x,n,3) # Up to third order polynomials\n",
+    "lamb = 0.025 # sparsification knob lambda\n",
+    "Xi = sparsifyDynamics(Theta,dx,lamb,n)\n",
+    "\n",
+    "print(Xi)"
+   ]
+  },
+  {
+   "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": 2
+}

CH07/CH07_SEC05_HAVOK_Lorenz.ipynb

@@ -0,0 +1,140 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import rcParams\n",
+    "from mpl_toolkits.mplot3d import Axes3D\n",
+    "from scipy import integrate\n",
+    "\n",
+    "\n",
+    "rcParams.update({'font.size': 18})\n",
+    "plt.rcParams['figure.figsize'] = [12, 12]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Simulate the Lorenz System\n",
+    "\n",
+    "dt = 0.01\n",
+    "T = 50\n",
+    "t = np.arange(0,T+dt,dt)\n",
+    "beta = 8/3\n",
+    "sigma = 10\n",
+    "rho = 28\n",
+    "\n",
+    "\n",
+    "def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):\n",
+    "    x, y, z = x_y_z\n",
+    "    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]\n",
+    "\n",
+    "np.random.seed(123)\n",
+    "x0 = (-8,8,27)\n",
+    "\n",
+    "x = integrate.odeint(lorenz_deriv, x0, t,rtol=10**(-12),atol=10**(-12)*np.ones_like(x0))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Eigen-time delay coordinates\n",
+    "stackmax = 10 # Number of shift-stacked rows\n",
+    "r = 10        # rank of HAVOK model\n",
+    "H = np.zeros((stackmax,x.shape[0]-stackmax))\n",
+    "\n",
+    "for k in range(stackmax):\n",
+    "    H[k,:] = x[k:-(stackmax-k),0]\n",
+    "    \n",
+    "U,S,VT = np.linalg.svd(H,full_matrices=0)\n",
+    "V = VT.T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Compute Derivatives (4th Order Central Difference)\n",
+    "# dV = np.zeros((V.shape[0]-5,r))\n",
+    "# for i in range(2,V.shape[0]-3):\n",
+    "#     for k in range(r):\n",
+    "#         dV[i-1,k] = (1/(12*dt))\n",
+    "\n",
+    "dV = (1/(12*dt)) * (-V[4:,:] + 8*V[3:-1,:] - 8*V[1:-3,:] + V[:-4,:])\n",
+    "\n",
+    "# trim first and last two that are lost in derivative\n",
+    "V = V[2:-2]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 46,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Build HAVOK Regression Model on Time Delay Coordinates\n",
+    "Xi = np.linalg.lstsq(V,dV,rcond=None)[0]\n",
+    "A = Xi[:(r-1),:(r-1)].T\n",
+    "B = Xi[-1,:(r-1)].T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0.16666666666666666\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(1/2/3)"
+   ]
+  },
+  {
+   "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.6.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

CH08/CH08_SEC07_2b_Obsv.ipynb

@@ -0,0 +1,124 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import rcParams\n",
+    "from control.matlab import *\n",
+    "import slycot\n",
+    "from scipy import integrate\n",
+    "from scipy.linalg import schur\n",
+    "# Python control toolbox available at https://python-control.readthedocs.io/\n",
+    "\n",
+    "plt.rcParams['figure.figsize'] = [8, 8]\n",
+    "plt.rcParams.update({'font.size': 18})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Observability matrix:\n",
+      "[[ 0.    0.    1.    0.  ]\n",
+      " [ 0.    0.    0.    1.  ]\n",
+      " [ 0.    0.1  -6.    0.  ]\n",
+      " [ 0.   -0.02  0.2  -6.  ]]\n",
+      "Observability matrix determinant: 0.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "m = 1\n",
+    "M = 5\n",
+    "L = 2\n",
+    "g = -10\n",
+    "d = 1\n",
+    "\n",
+    "b = -1 # pendulum down (b = -1)\n",
+    "\n",
+    "A = np.array([[0,1,0,0],\\\n",
+    "              [0,-d/M,b*m*g/M,0],\\\n",
+    "              [0,0,0,1],\\\n",
+    "              [0,-b*d/(M*L),-b*(m+M)*g/(M*L),0]])\n",
+    "\n",
+    "B = np.array([0,1/M,0,b/(M*L)]).reshape((4,1))\n",
+    "\n",
+    "C = np.array([0,0,1,0]) # only observable if x measured... because x can't be\n",
+    "\n",
+    "print('Observability matrix:\\n{}'.format(obsv(A,C)))\n",
+    "print('Observability matrix determinant: {}'.format(np.linalg.det(obsv(A,C))))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Observability matrix:\n",
+      "[[ 0.   1.   0. ]\n",
+      " [ 0.   0.   1. ]\n",
+      " [ 0.1 -6.   0. ]]\n",
+      "Gramian determinant: 0.031250000000001\n"
+     ]
+    }
+   ],
+   "source": [
+    "## Which measurements are best if we omit \"x\"\n",
+    "Ah = A[1:,1:]\n",
+    "Bh = B[1:]\n",
+    "# Ch = np.array([1,0,0])\n",
+    "Ch = np.array([0,1,0])\n",
+    "# Ch = np.array([0,0,1])\n",
+    "\n",
+    "print('Observability matrix:\\n{}'.format(obsv(Ah,Ch)))\n",
+    "\n",
+    "Ch = Ch.reshape((1,len(Ch)))\n",
+    "Dh = np.zeros((Ch.shape[0],Bh.shape[1]))\n",
+    "sys = ss(Ah,Bh,Ch,Dh)\n",
+    "print('Gramian determinant: {}'.format(np.linalg.det(gram(sys,'o'))))\n"
+   ]
+  },
+  {
+   "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": 2
+}