{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_11_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log\n",
    "import numpy as np\n",
    "from numpy import mod\n",
    "import scipy.sparse.linalg as las\n",
    "import scipy.interpolate as ip\n",
    "import scipy.spatial as sp\n",
    "from scipy import sparse\n",
    "import scipy.linalg as la\n",
    "import scipy.signal as sg\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "from matplotlib import patches as pch\n",
    "import timeit\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if( N==0 ):\n",
    "            continue; # skip vector of zero length\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "\n",
    "# sigma function, works with floats and Nx1 arrays\n",
    "def eda_sigma(w,x):\n",
    "    if( isinstance(x,np.ndarray) ):\n",
    "        N,M = np.shape(x);\n",
    "        y = np.zeros((N,1));\n",
    "        for i in range(N):\n",
    "            if( x[i,0] >= 0.0):\n",
    "                y[i,0] = 1/( 1 + exp(-w*x[i,0]) );\n",
    "            else:\n",
    "                t = exp(w*x[i,0]);\n",
    "                y[i,0] = t / (t + 1);\n",
    "    elif( isinstance(x,float) ):\n",
    "        y = 0.0;\n",
    "        if( x >= 0.0):\n",
    "            y = 1/( 1 + exp(-w*x) );\n",
    "        else:\n",
    "            t = exp(w*x);\n",
    "            y = t / (t + 1);\n",
    "    return(y);\n",
    "\n",
    "def eda_net_init(Lmax,Nmax):\n",
    "# initialize a network; produces only arrays of zeros\n",
    "# input:\n",
    "#        Lmax: number of layers\n",
    "#        Nmax: maximum number of neurons on any layer\n",
    "# output:\n",
    "#        N: neurons on a layer array\n",
    "#        w: weights\n",
    "#        b: biases\n",
    "#        a: activities\n",
    "\n",
    "    # N[0:Lmax,1]: column-vector with number of neurons in each layer;\n",
    "    N = np.zeros((Lmax,1),dtype=np.int8);\n",
    "\n",
    "    # biases: b[0:N[i,0],i) is a column-vector that gives the\n",
    "    # biases of the N(i) neurons on the i-th layer\n",
    "    b = np.zeros( (Nmax, Lmax) );\n",
    "\n",
    "    # weights: w[0:N[i,0],0:N[i-1,0],i) is a matrix that gives the\n",
    "    # weights between the N[i] neurons in the i-th layer and the\n",
    "    # N[i-1] neurons on the (i-1)-th layer. The weights assigned\n",
    "    # to the first layer are never used and should be set to zero\n",
    "    w = np.zeros( (Nmax, Nmax, Lmax) );\n",
    "\n",
    "    # activity: a[0:N[i,0],i) is a column-vector that gives the\n",
    "    # activitites of the N[i] neurons on the i-th layer\n",
    "    a = np.zeros( (Nmax, Lmax) );\n",
    "    \n",
    "    return N, w, b, a;\n",
    "\n",
    "def eda_net_view(N,w,b,fig,ax):\n",
    "# makes a picture of a neural net\n",
    "# you must set the figure before the call, e.g.\n",
    "# fig = plt.figure(1);\n",
    "# ax = plt.subplot(1,1,1);\n",
    "    Nmax = np.max(N);\n",
    "    Lmax, tmp = np.shape(N);\n",
    "\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig.set_size_inches(W,H);\n",
    "    ax.axis('off')\n",
    "    ax.axis([-0.5, Lmax-0.5, -0.5, Nmax-0.5]);\n",
    "\n",
    "    scale = np.zeros((Lmax,1));\n",
    "    for x in range(Lmax):\n",
    "        scale[x,0] =  Nmax/N[x,0];\n",
    "\n",
    "    for x in range(Lmax):\n",
    "        for y in range(N[x,0]):\n",
    "            lowerleft = ( x-0.2, y*scale[x,0]-0.4 );\n",
    "            width=0.4;\n",
    "            height = 0.8;\n",
    "            c=0.9;\n",
    "            ax.add_patch( pch.Rectangle( lowerleft, width, height, linewidth=2, fc=(c,c,c), ec=(0.0,0.0,0.0) ) );\n",
    "            strx = x;\n",
    "            stry = y*scale[x,0];\n",
    "            mystr = 'b = %.2f' % (b[y,x]);\n",
    "            ax.text(strx, stry, mystr, ha='center' );\n",
    "            if ( x != (Lmax-1) ):\n",
    "                for y2 in range(N[x+1,0]):\n",
    "                    if ( w[ y2 , y , x+1] != 0 ):\n",
    "                        ax.plot([x+0.2, x+0.8],[y*scale[x,0], y2*scale[x+1,0]],color=(0.7,0.7,0.7) );\n",
    "                        strx = x+0.5;\n",
    "                        stry = 0.07+(y*scale[x,0]+y2*scale[x+1,0])/2;\n",
    "                        mystr = 'w = %.2f' % (w[ y2 , y , x+1]);\n",
    "                        ax.text(strx, stry, mystr, ha='center' );\n",
    "                                               \n",
    "def eda_net_eval(N,w,b,a):\n",
    "# eda_eval_net evaluates a network with properites N, w and b \n",
    "# the \"input\" activities a[0:N[0,1],0] must be set\n",
    "# the output are the final activities and the final z's\n",
    "    Nmax = np.max(N);\n",
    "    Lmax, tmp = np.shape(N);\n",
    "    z = np.zeros( (Nmax, Lmax) );\n",
    "    for i in range(1,Lmax):\n",
    "        M = w[0:N[i,0],0:N[i-1,0],i];\n",
    "        v = a[0:N[i-1,0],i-1];\n",
    "        u = b[0:N[i,0],i];\n",
    "        zz = np.zeros((N[i,0],1));\n",
    "        zz[:,0] =  np.matmul( M, v ) + u;\n",
    "        z[0:N[i,0],i] = zz[:,0];\n",
    "        if i < (Lmax-1):\n",
    "            a[0:N[i,0],i] = eda_sigma(1.0,zz).ravel();\n",
    "        else:\n",
    "            a[0:N[i,0],i] = z[0:N[i,0],i];\n",
    "    return a, z;\n",
    "\n",
    "def eda_net_1dtower( Nc, slope, Xc, Dx, h, N, w, b ):\n",
    "# eda_net_linear\n",
    "# defines a network that computes a set of 1D towers\n",
    "# input Nc: number of towers\n",
    "#       slope: slope of all the towers\n",
    "#       Xc: centers of the towers\n",
    "#       Dx: width of the towers\n",
    "#       h: heights of the towaers\n",
    "#       N: layer array from eda_net_init()\n",
    "#       w: weights array from eda_net_init()\n",
    "#       b: bias array from from eda_net_init()\n",
    "# output: updated N, w, b arrays\n",
    "\n",
    "    N[:,0] = [1, Nc*2, 1];\n",
    "\n",
    "    for ncvalue in range(Nc):\n",
    "        w12 = 4*slope;\n",
    "        w22 = 4*slope;\n",
    "        w[ ((ncvalue)*2):((ncvalue)*2)+2, 0:N[0,0],  1] = np.reshape( [w12,w22], (2,1) );\n",
    "\n",
    "        w13 =  h[ncvalue,0];\n",
    "        w23 = -h[ncvalue,0];\n",
    "        w[ 0:N[2,0], ((ncvalue)*2):((ncvalue)*2)+2, 2] = np.reshape( [w13,w23], (1,2) );\n",
    "\n",
    "        b12 = -w12 * ( Xc[ncvalue,0] - Dx[ncvalue,0]/2 );\n",
    "        b22 = -w22 * ( Xc[ncvalue,0] + Dx[ncvalue,0]/2 );\n",
    "        b[ ((ncvalue)*2):((ncvalue)*2)+2, 1 ] = [b12, b22];\n",
    "\n",
    "    return N, w, b;\n",
    "\n",
    "\n",
    "def eda_net_2dtower( Nc, slope, Xc, Yc, Dx, Dy, h, N, w, b ):\n",
    "# defines a network that superimposes several towers\n",
    "# input:\n",
    "#    Nc: number of towers\n",
    "#    slope: slope of all the towers\n",
    "#    Xc, Yc: vector of the centers of the towers in x and y\n",
    "#    Dx, Dy: vector of the widths of the towers in x and y\n",
    "#    h: vector of heights of the towers \n",
    "#    N: layer array from eda_net_init()\n",
    "#    w: weights array from eda_net_init()\n",
    "#    b: bias array from from eda_net_init()\n",
    "# output: updated N, w, b arrays\n",
    "\n",
    "    N[:,0] = [2, Nc*4, Nc, 1];\n",
    "\n",
    "    for ncvalue in range(Nc):\n",
    "        w12 = 4*slope;\n",
    "        w22 = 4*slope;\n",
    "        w32 = 4*slope;\n",
    "        w42 = 4*slope;\n",
    "        w[ ((ncvalue)*4):((ncvalue)*4)+2,   0 ,1] = [w12, w22];\n",
    "        w[ ((ncvalue)*4)+2:((ncvalue)*4)+4, 1 ,1] = [w32, w42];\n",
    "    \n",
    "        w13 = 4*slope;\n",
    "        w23 = -4*slope;\n",
    "        w33 = 4*slope;\n",
    "        w43 = -4*slope;\n",
    "        w[ncvalue, ((ncvalue)*4):((ncvalue)*4)+4,2] = [w13, w23, w33, w43];\n",
    "    \n",
    "        w14 = h[ncvalue,0];\n",
    "        w[0, ncvalue, 3] = w14;\n",
    "\n",
    "        b12 = -w12 * ( Xc[ncvalue,0] - Dx[ncvalue,0]/2);\n",
    "        b22 = -w22 * ( Xc[ncvalue,0] + Dx[ncvalue,0]/2);\n",
    "        b32 = -w32 * ( Yc[ncvalue,0] - Dy[ncvalue,0]/2);\n",
    "        b42 = -w42 * ( Yc[ncvalue,0] + Dy[ncvalue,0]/2);\n",
    "        b[ ((ncvalue)*4):((ncvalue)*4)+4 ,1] = [b12, b22, b32, b42];\n",
    "    \n",
    "        b13 = -(3*(4*slope))/2;\n",
    "        b[ncvalue,2] = b13;\n",
    "\n",
    "    return N, w, b;\n",
    "\n",
    "def eda_net_map(N,w,b):\n",
    "# builds look-up tables to allow sequential ordering\n",
    "# of non-trivual biases and weights\n",
    "# used to build list of model parameters used in training\n",
    "# input:\n",
    "#        N: neurons on layer array\n",
    "#        w: weights\n",
    "#        b: biases\n",
    "# output\n",
    "#        Nb: number of biases\n",
    "#        Nw: Number of non-zero weights\n",
    "#            (use near-zero weights to fool it)\n",
    "#        layer_of_bias, vector(i) specifying layer of bias i\n",
    "#        neuron_of_bias, vector(i) specifying neuron of bias i\n",
    "#        index_of_bias, 2d array(i,j) specifying number of bias\n",
    "#            of layer j, neuron i\n",
    "#         layer_of_weight, vector(i) specifying layer of weight i\n",
    "#         r_neuron_of_weight, vector(i) specifying right neuron of weight i\n",
    "#         l_neuron_of_weight, vector(i) specifying left neuron of weight i\n",
    "#         index_of_weight, array(i,j,k) specifying number of weight(i,j,k)\n",
    "#         note: all the arrays are set to their maximum size, with unused indices set to (-1)\n",
    "\n",
    "    Nmax = np.max(N);\n",
    "    Lmax, tmp = np.shape(N);\n",
    "\n",
    "    layer_of_bias = -np.ones((Lmax*Nmax,1),dtype=np.int);\n",
    "    neuron_of_bias = -np.ones((Lmax*Nmax,1),dtype=np.int);\n",
    "    index_of_bias = -np.ones((Nmax,Lmax),dtype=np.int);\n",
    "    \n",
    "    Nb=0;  # number of nontrivial biases\n",
    "    for k in range(1,Lmax):\n",
    "        for j in range(N[k,0]):\n",
    "            layer_of_bias[Nb,0]=k;\n",
    "            neuron_of_bias[Nb,0]=j;\n",
    "            index_of_bias[j,k]=Nb;\n",
    "            Nb = Nb+1;\n",
    "\n",
    "    Nw=0;  # number of nontrivial weights\n",
    "    layer_of_weight = -np.ones((Lmax*Nmax*Nmax,1),dtype=np.int);\n",
    "    r_neuron_of_weight = -np.ones((Lmax*Nmax*Nmax,1),dtype=np.int);\n",
    "    l_neuron_of_weight = -np.ones((Lmax*Nmax*Nmax,1),dtype=np.int);\n",
    "    index_of_weight = -np.ones((Nmax,Nmax,Lmax),dtype=np.int);\n",
    "\n",
    "    for k in range(1,Lmax):\n",
    "        for i in range(N[k,0]):\n",
    "            for j in range(N[k-1,0]):\n",
    "                if( w[i,j,k] != 0.0 ):\n",
    "                    layer_of_weight[Nw,0]=k;\n",
    "                    r_neuron_of_weight[Nw,0]=i;\n",
    "                    l_neuron_of_weight[Nw,0]=j;\n",
    "                    index_of_weight[j,i,k]=Nw;\n",
    "                    Nw = Nw+1;\n",
    "\n",
    "    return Nb, Nw, layer_of_bias, neuron_of_bias, index_of_bias, layer_of_weight, r_neuron_of_weight, l_neuron_of_weight, index_of_weight; \n",
    "\n",
    "def eda_net_deriv( N,w,b,a ):\n",
    "# calculates the derivatives daLmax/dw and daLmax/dw\n",
    "# for a network with parameters N,w,b,a\n",
    "# and returns daLmaxdw, daLmaxdb\n",
    "\n",
    "    # N[i]:  number of neurons in layer i\n",
    "    \n",
    "    # Nmax: maximum number of neurons in any layer\n",
    "    Nmax = np.max(N);\n",
    "    \n",
    "    # Lmax number of layers\n",
    "    Lmax, tmp = np.shape(N);\n",
    "    \n",
    "    # b[i,j]: bias of neuron i in layer j\n",
    "    # b = np.zeros( (Nmax, Lmax) );\n",
    "\n",
    "    # w[i,j,k]: weights between neuron i in the k-th layer and\n",
    "    # neuron j in the (k-1)-th layer. The weights assigned\n",
    "    # to the first layer are never used and should be set to zero\n",
    "    # w = np.zeros( (Nmax, Nmax, Lmax) );\n",
    "\n",
    "    # activity: a[0:N[i,0],i) is a column-vector that gives the\n",
    "    # activitites of the N[i] neurons on the i-th layer\n",
    "    # a = np.zeros( (Nmax, Lmax) );\n",
    "    \n",
    "    # dada[i,j,k]: the change in the activity of the i-th neuron of the layer k\n",
    "    # due to a change in the activity of the j-th neuron in the layer k-1\n",
    "    dada = np.zeros((Nmax,Nmax,Lmax));\n",
    "\n",
    "    # daLmaxda[i,j,k]: the change in the activity of the i-th neuron in layer Lmax-1\n",
    "    # due to a change in the activity of the j-th neuron in the layer k-1\n",
    "    daLmaxda = np.zeros((Nmax,Nmax,Lmax));\n",
    "\n",
    "    # dadz[i,k]: the change in the activity of the i-th neuron of layer k\n",
    "    # due to a change in the input z of the same neuron\n",
    "    dadz = np.zeros((Nmax,Lmax));\n",
    "    \n",
    "    # daLmaxdb[i,j,k]: the change in the activity of the i-th neuron in layer Lmax-1\n",
    "    # due to a change in the bias of the of the j-th neuron in the layer k\n",
    "    daLmaxdb = np.zeros((Nmax,Nmax,Lmax));\n",
    "\n",
    "    # daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax-q\n",
    "    # due to a change in the weight connecting the j-th neuron in the layer l\n",
    "    # with the k-th neuron of layer l-1\n",
    "    daLmaxdw = np.zeros((Nmax,Nmax,Nmax,Lmax));\n",
    "\n",
    "    # da/dz\n",
    "    for i in range(Lmax):\n",
    "        if (i==(Lmax-1)):\n",
    "            # sigma function is sigma(z)=z for top layer\n",
    "            # d sigma / dz = 1\n",
    "            # MATLAB: dadz(1:N(i),i) = ones(N(i),1);\n",
    "            dadz[0:N[i,0],i] = np.ones((N[i,0],1)).ravel();\n",
    "        else:\n",
    "            # sigma function is sigma(z)=(1+exp(-z))**(-1) for other layers\n",
    "            # d sigma / dz = (-1)(-1) exp(-z)(1+exp(-z))^(-1)\n",
    "            #              = exp(-z) (a**2)\n",
    "            # and note exp(-z) = (1/a)-1\n",
    "            # so d sigma / dz = a-a**2\n",
    "            # MATLAB: dadz(1:N(i),i) = a(1:N(i),i) - (a(1:N(i),i).^2);\n",
    "            dadz[0:N[i,0],i] = a[0:N[i,0],i] - np.power( a[0:N[i,0],i], 2);  \n",
    "    \n",
    "    # da/da\n",
    "    for i in range(Lmax-1,0,-1): # MATLAB: [Lmax:-1:2]\n",
    "        if ( i == (Lmax-1) ):\n",
    "            # a=z in the top of the layer, and\n",
    "            # z = (sum) w a + b\n",
    "            # so da(layer i)/da(layer i-1) = dz/da(layer i-1) = w\n",
    "            # MATLAB: dada(1:N(i),1:N(i-1),i) = w(1:N(i),1:N(i-1),i);\n",
    "            dada[0:N[i,0],0:N[i-1,0],i] = w[0:N[i,0],0:N[i-1,0],i];\n",
    "            # initialize daLmaxda to prepare for backpropagation\n",
    "            # MATLAB: daLmaxda(1:N(i),1:N(i-1),i) = dada(1:N(i),1:N(i-1),i);\n",
    "            daLmaxda[0:N[i,0],0:N[i-1,0],i] = dada[0:N[i,0],0:N[i-1,0],i];\n",
    "        else:\n",
    "            # a=sigma(z) in the other layers, so use the chain rule\n",
    "            # da/da = da(layer)/dz * dz/da(layer i-1)\n",
    "            #    where as before dz/da(layer i-1) = w\n",
    "            # MATLAB: dada(1:N(i),1:N(i-1),i) = diag(dadz(1:N(i),i))*w(1:N(i),1:N(i-1),i);\n",
    "            V1 = np.zeros((N[i,0],1));\n",
    "            V1[:,0] = dadz[0:N[i,0],i];\n",
    "            M1 = np.diag(V1.ravel());\n",
    "            M2 = np.zeros((N[i,0],N[i-1,0]));\n",
    "            M2[:,:] = w[0:N[i,0],0:N[i-1,0],i];\n",
    "            dada[0:N[i,0],0:N[i-1,0],i] = np.matmul( M1, M2 );\n",
    "            # backpropagate\n",
    "            # MATLAB: daLmaxda(1:N(Lmax),1:N(i-1),i) =  ...\n",
    "            #                        daLmaxda(1:N(Lmax),1:N(i),i+1)*dada(1:N(i),1:N(i-1),i);\n",
    "            M1 = np.zeros((N[Lmax-1,0],N[i,0]));\n",
    "            M1[:,:] = daLmaxda[0:N[Lmax-1,0],0:N[i,0],i+1];\n",
    "            M2 = np.zeros((N[i,0],N[i-1,0]));\n",
    "            M2[:,:] = dada[0:N[i,0],0:N[i-1,0],i];\n",
    "            daLmaxda[0:N[Lmax-1,0],0:N[i-1,0],i] = np.matmul(M1, M2);\n",
    "        \n",
    "    # da/db\n",
    "    for i in range(Lmax):\n",
    "        if( i==(Lmax-1) ):\n",
    "            # a=z in the top layer, and z(layer i)= w *a(layer i-1) + b\n",
    "            # so da(layer i)/db = dz(layer i)/db = 1\n",
    "            daLmaxdb[0:N[Lmax-1,0],0:N[Lmax-1,0],Lmax-1] = np.identity(N[Lmax-1,0]);\n",
    "        else:\n",
    "            # a = sigma(z) in the other layers, so use the chain rule\n",
    "            # da(layer i)/db = da(layer i)/dz * dz/db\n",
    "            # where as before dzdb=1\n",
    "            dzdb = 1;\n",
    "            dadb = np.zeros((N[i,0],1));\n",
    "            dadb[:,0] = dadz[0:N[i,0],i] * dzdb;\n",
    "            # then apply the chain rule for activities\n",
    "            # da(layer Lmax)/db = da(layer Lmax)/da(layer i) * da(layer i)/db\n",
    "            # MATLAB: daLmaxdb(1:N(Lmax),1:N(i),i) = daLmaxda(1:N(Lmax),1:N(i),i+1) * diag(dadb);\n",
    "            M1 = np.zeros((N[Lmax-1,0],N[i,0]));\n",
    "            M1[:,:]=daLmaxda[0:N[Lmax-1,0],0:N[i,0],i+1];\n",
    "            M2 = np.diag(dadb.ravel());\n",
    "            daLmaxdb[0:N[Lmax-1,0],0:N[i,0],i] = np.matmul( M1, M2 );\n",
    "\n",
    "    # da/dw\n",
    "    # calculation of weight derivatives\n",
    "    # since z(layer i) = w*a(layer i-1) + b\n",
    "    # dz(layer i)/dw = a(layer i-1)\n",
    "    for l in range(1,Lmax): # MATLAB: [2:Lmax]\n",
    "        if( l==(Lmax-1) ):\n",
    "            # in the top layer, a=z, so da/dw = dz/dw = a\n",
    "            for j in range(N[l,0]):\n",
    "                daLmaxdw[j,j,0:N[l-1,0],l] = a[0:N[l-1,0],l-1];\n",
    "        else:\n",
    "            # in the other layers, use chain rule\n",
    "            # da/dw = da/dz * dz/dw\n",
    "            # and then the chain ruke again\n",
    "            # da(layer Lmax)/dw = da(layer Lmax)/da(layer i) * da(layer i)/dw\n",
    "            for j in range(N[l,0]):\n",
    "                dzdw = np.zeros((N[l-1,0],1));\n",
    "                dzdw[:,0] = a[0:N[l-1,0],l-1];\n",
    "                dadw = dadz[j,l] * dzdw;\n",
    "                # MATLAB: daLmaxdw(1:N(Lmax),j,1:N(l-1),l) = daLmaxda(1:N(Lmax),j,l+1) * dadw';\n",
    "                V1 = np.zeros((N[Lmax-1,0],1));\n",
    "                V1[:,0] = daLmaxda[0:N[Lmax-1,0],j,l+1];\n",
    "                M1 = np.matmul( V1, dadw.T );\n",
    "                daLmaxdw[0:N[Lmax-1,0],j,0:N[l-1,0],l] = M1;\n",
    "                \n",
    "    return daLmaxdw, daLmaxdb;\n",
    "\n",
    "def eda_net_linear( Nc, slope1, slope2, c, h, N, w, b ):\n",
    "# defines a network that cpmputes the linear function\n",
    "# d(x1, x2, ...) = c + h1*x1 + h2*x2 + h2*h3 ...\n",
    "# (good for implementing a filter)\n",
    "# input Nc: number of linear components\n",
    "#       slope1: slope of the first stage (say 0.01)\n",
    "#       slope2: slope of the second stage (say 0.02)\n",
    "#              should be different from slope1\n",
    "#       c: intercept\n",
    "#       h: slopes\n",
    "#       N: layer array from eda_net_init()\n",
    "#       w: weights array from eda_net_init()\n",
    "#       b: bias array from from eda_net_init()\n",
    "# output: updated N, w, b arrays\n",
    "\n",
    "    Nmax = np.max(N);\n",
    "    Lmax, tmp = np.shape(N);\n",
    "    N[:,0] = [Nc, Nc, Nc, 1];\n",
    "\n",
    "    b[0,3] = c;\n",
    "\n",
    "    for ncvalue in range(Nc):\n",
    "\n",
    "        w12 = 4*slope1;\n",
    "        w[ncvalue, ncvalue ,1] = w12;\n",
    "  \n",
    "        w13 = 4*slope2;\n",
    "        w[ncvalue, ncvalue ,2] = w13;\n",
    "    \n",
    "        w2w3 = w12*w13;\n",
    "    \n",
    "        w14 = 16*h[ncvalue,0]/w2w3;\n",
    "        w[0,ncvalue, 3] = w14;\n",
    "\n",
    "        b12 = 0;\n",
    "        b[ncvalue,1] = b12;\n",
    "    \n",
    "        b13 = -w13/2;\n",
    "        b[ncvalue,2] = b13;\n",
    "    \n",
    "        b14 = -8*h[ncvalue,0]/w2w3;\n",
    "        b[0,3] = b[0,3]+ b14;\n",
    "\n",
    "    return N, w, b;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28164 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_01, example of a Taylor Series\n",
    "# using the function\n",
    "# 1/(1-x) = 1 + x + x^2 + x^3 + ...\n",
    "\n",
    "# standard setup of time axis, t\n",
    "N = 101;       # number of points\n",
    "tmin = -1.2;   # start time\n",
    "tmax = 0.8;    # end time\n",
    "Dt = (tmax-tmin)/(N-1);  # sampling interval\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector\n",
    "\n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([tmin, tmax, 0, 2]);\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('y(t)');\n",
    "\n",
    "# true function\n",
    "y = 1.0/(1.0-t);\n",
    "plt.plot(t,y,'-',linewidth=5,color=(0.8,0.8,0.8));\n",
    "t0=0.0;\n",
    "y0=1./(1.0-t0);\n",
    "plt.plot(t0,y0,'ko',linewidth=2);\n",
    "\n",
    "# linear appriximation\n",
    "y1 = 1.0+t;\n",
    "plt.plot(t,y1,'k:',linewidth=2);\n",
    "\n",
    "# quadratic approximation\n",
    "y2 = 1.0+t+np.power(t,2);\n",
    "plt.plot(t,y2,'k--',linewidth=1.5);\n",
    "\n",
    "# cubic approximation\n",
    "y3 = 1+t+np.power(t,2)+np.power(t,3);\n",
    "plt. plot(t,y3,'k-',linewidth=3);\n",
    "           \n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21160 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_02, approximation for small distances on a sphere\n",
    "\n",
    "# degree--radian conversion factors\n",
    "RTOD = 180.0/pi;\n",
    "DTOR = pi/180.0;\n",
    "\n",
    "# the object is to compute distances between (lat1, lon1) and (lat2, lon2)\n",
    "# (lat1, lon1) is a single point\n",
    "# (lat2, lon2) are vectors with N points\n",
    "N = 101;  # number of points\n",
    "lat1 = 30.0*DTOR; # lat of single point\n",
    "lon1 = 30.0*DTOR; # lon of single point\n",
    "lat2min = 30.0*DTOR; # start lat of vector\n",
    "lat2max = 40.0*DTOR; # end lat of vector\n",
    "lat2 = eda_cvec( np.linspace(lat2min,lat2max,N) ); # vector of lats\n",
    "lon2 = lat2; # vector of lons\n",
    "latbar = (lat1 + lat2)/2; # vector of mean lat\n",
    "Dlon = lon2 - lon1;       # vector of difference in lat with respect to first endpoint\n",
    "Dlat = lat2 - lat1;       # vector of difference in lon with respect to first endpoint\n",
    "\n",
    "# law of cosines on a sphere\n",
    "cr = sin(lat1)*np.sin(lat2) + np.multiply( cos(lat1)*np.cos(lat2), np.cos(Dlon));\n",
    "r = np.arccos(cr);\n",
    "\n",
    "# flat earth approximation\n",
    "coslatbar2 = np.power( np.cos(latbar), 2);\n",
    "Dlat2 = np.power(Dlat,2);\n",
    "Dlon2 = np.power(Dlon,2);\n",
    "t1 = Dlat2 + np.multiply( Dlon2, coslatbar2 );\n",
    "r1 = np.sqrt( t1 );\n",
    "\n",
    "# higher order correction\n",
    "t2 = np.multiply( Dlat2, np.multiply( Dlon2, (0.25-(coslatbar2/6)) ) );\n",
    "r2 = np.sqrt( t1-t2 );\n",
    "\n",
    "# plot distances\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([30, 40, 0, 15]);\n",
    "plt.xlabel('latitude (deg)');\n",
    "plt.ylabel('r (deg)');\n",
    "plt.plot( RTOD*lat2, RTOD*r, '-', linewidth=7, color=(0.85,0.85,0.85) );\n",
    "plt.plot( RTOD*lat2, RTOD*r1, '-', linewidth=3, color=(0.7,0.7,0.7));\n",
    "plt.plot( RTOD*lat2, RTOD*r2, 'k-', linewidth=1 );\n",
    "\n",
    "# plot error in distances\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([30, 40, 0, 0.02]);\n",
    "plt.xlabel('latitude (deg)');\n",
    "plt.ylabel('error (deg)');\n",
    "plt.plot( RTOD*lat2, RTOD*(r1-r), 'k:', 'LineWidth', 2 );\n",
    "plt.plot( RTOD*lat2, RTOD*(r2-r), 'k-', 'LineWidth', 2 );\n",
    "\n",
    "plt.show();                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 16780 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_03, variance of nonlinear function\n",
    "# variance of T=1/f computed by:\n",
    "#        small number approximation\n",
    "#        numerical simulation\n",
    "\n",
    "# two time pi\n",
    "tpi = 2*pi;\n",
    "\n",
    "# w is angular frequency \"omega\"\n",
    "wbar = 1.3*tpi;       # mean\n",
    "sigmaw = 0.2;         # standard deviation\n",
    "sigmaw2 = sigmaw**2;  # variance\n",
    "\n",
    "# realizations of a Normally-distributed random process\n",
    "N = 10000;\n",
    "wr = np.random.normal(wbar,sigmaw,(N,1));\n",
    "\n",
    "# histogram of wr\n",
    "Nbins = 201;                # number of bins in histogram\n",
    "wbmin = 0;                  # start bin\n",
    "wbmax = 20;                 # end bin\n",
    "c, e = np.histogram(wr,Nbins,(wbmin,wbmax)); # make histogram\n",
    "Nc = len(c);                # length of count c\n",
    "Ne = len(e);                # length of edges e\n",
    "count = eda_cvec(c);        # count in each bin\n",
    "edges = eda_cvec(e[0:Ne-1]); # bin edge locations\n",
    "wb =  eda_cvec( 0.5*(e[0:Nbins]+e[1:Nbins+1]) ); # bin centers\n",
    "Dwb = wb[1,0]-wb[0,0]; # bin spacing\n",
    "\n",
    "# true Normal p.d.f.\n",
    "pw1 = (1/(sqrt(tpi)*sigmaw))*np.exp(-np.power(wb-wbar,2)/(2*sigmaw2));\n",
    "\n",
    "# approximate p.d.f from histogram\n",
    "pw2 = count / (Dwb*np.sum(count));\n",
    "\n",
    "# plot omega p.d.f.'s\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([0, 10, 0, 3]);\n",
    "plt.xlabel('w (radians/s)');\n",
    "plt.ylabel('p(w)');\n",
    "plt.plot( wb, pw1, '-', linewidth=4, color=(0.8,0.8,0.8) );\n",
    "plt.plot( wb, pw2, 'k:', linewidth=2 );\n",
    "\n",
    "# now consider period T = 2 pi / omega\n",
    "\n",
    "# approxmate mean and variance from linear transformation\n",
    "Tbar = tpi/wbar;\n",
    "sigmaT = sigmaw*tpi/(wbar**2);\n",
    "sigmaT2 = sigmaT**2;\n",
    "\n",
    "# convert realization to period\n",
    "Tr = np.divide( tpi, wr );\n",
    "\n",
    "# histogram or Tr\n",
    "Nbins = 201;                 # number of bins in histogram\n",
    "Tbmin = 0;                   # start bin\n",
    "Tbmax = 2;                   # end bin\n",
    "c, e = np.histogram(Tr,Nbins,(Tbmin,Tbmax)); # make histogram\n",
    "Nc = len(c);                 # length of counts c\n",
    "Ne = len(e);                 # length of edges e\n",
    "count = eda_cvec(c);         # count vector\n",
    "edges = eda_cvec(e[0:Ne-1]); # edges vector\n",
    "Tb = eda_cvec(0.5*(e[0:Nbins]+e[1:Nbins+1])); # bin centers\n",
    "DTb = Tb[1,0]-Tb[0,0];       # bin spacing\n",
    "\n",
    "# approximate p.d.f from histogram\n",
    "pT2 = count / (DTb*np.sum(count));\n",
    "\n",
    "# Normal p.d.f. based on approximate mean and variance\n",
    "pT1 = (1/(sqrt(tpi)*sigmaT))*np.exp(-np.power(Tb-Tbar,2)/(2*sigmaT2));\n",
    "\n",
    "# plot omega p.d.f.'s\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([0, 1, 0, 25]);\n",
    "plt.xlabel('T (s))');\n",
    "plt.ylabel('p(T)');\n",
    "plt.plot( Tb, pT1, '-', linewidth=4, color=(0.8,0.8,0.8) );\n",
    "plt.plot( Tb, pT2, 'k:', linewidth=2 );\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 63712 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18676 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations 4\n",
      "posterior error 5.619 (deg C)\n",
      "Estimated period 365.613 +/- 0.060 (95%) days\n"
     ]
    }
   ],
   "source": [
    "# edapy_11_04, fit of sinusoid of unknown frequency\n",
    "# to Black Rock Forest temperature data, using\n",
    "# fit of a sinusoid of unknown (but approximately annual)\n",
    "# iterative least squares method (Newton's method)\n",
    "\n",
    "# load Black Rock Forest temperature data\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "\n",
    "# clean the data by deleting outliers and zero (= no data) values\n",
    "test1 = (D[:,1]>-35);\n",
    "test2 = (D[:,1]<40);\n",
    "test3 = (D[:,1]!=0.0);\n",
    "r = np.where( test1 & test2 & test3);\n",
    "k = r[0];\n",
    "tmp = np.shape(k);\n",
    "N = tmp[0];\n",
    "t = D[k,0:1];    # cleaned time vector t of length N\n",
    "dobs = D[k,1:2]; # cleaned observation vector dobs of length N\n",
    "\n",
    "# some statistics\n",
    "dmean = np.mean(dobs);                  # mean\n",
    "amp = 0.5*(np.max(dobs)-np.min(dobs));  # amplitude\n",
    "\n",
    "# Annual period, frequency, angular frequency\n",
    "Ta = 365.0;\n",
    "fa = 1.0/Ta;\n",
    "wa = 2.0*pi*fa;\n",
    "\n",
    "# nonlinear theory and its derivative\n",
    "# note model parameters are scaled to be order unity\n",
    "# d = m1*dmean + amp*m2*sin(m4 wa t) + amp*m3*cos(m4 wa t)\n",
    "# dd/dm1 =  dmean\n",
    "# dd/dm2 =  amp*sin(m4 wa t)\n",
    "# dd/dm3 =  amp*cos(m4 wa t)\n",
    "# dd/dm4 =  amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)\n",
    "\n",
    "# inital guesses\n",
    "M=4;\n",
    "m = np.zeros((M,1));\n",
    "m[0,0] =  1.0;        # average level is mean\n",
    "m[1,0] =  0.05;       # phase with max in summer\n",
    "m[2,0] = -0.95;       # so only -cos() term\n",
    "m[3,0] =  1.0;        # annual frequency \n",
    "\n",
    "s = amp*np.sin(m[3,0]*wa*t);\n",
    "c = amp*np.cos(m[3,0]*wa*t);\n",
    "dpre = dmean*m[0,0] + m[1,0]*s + m[2,0]*c;\n",
    "\n",
    "# plot data and initial prediction\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([0,5000,-30,50]);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (deg C)');\n",
    "plt.plot( t, dobs, 'k-', linewidth=2 );\n",
    "plt.plot( t, dpre, '-', linewidth=2, color=(0.7,0.7,0.7) );\n",
    "\n",
    "# Newton's method\n",
    "s = amp*np.sin(m[3,0]*wa*t);\n",
    "c = amp*np.cos(m[3,0]*wa*t);\n",
    "dpre = dmean*m[0,0] + m[1,0]*s + m[2,0]*c;\n",
    "dd = dobs - dpre;\n",
    "E = np.matmul(dd.T,dd);\n",
    "Niter = 10;\n",
    "ic=0;\n",
    "for i in range(Niter):\n",
    "    ic=ic+1;\n",
    "    Eprev = E;\n",
    "    G = np.zeros((N,4));\n",
    "    G[:,0] = (dmean*np.ones((N,1))).ravel();\n",
    "    G[:,1] = s.ravel();\n",
    "    G[:,2] = c.ravel();\n",
    "    G[:,3] = (m[1,0]*wa*np.multiply(t,c)-m[2,0]*wa*np.multiply(t,s)).ravel();\n",
    "    GTG=np.matmul(G.T, G);\n",
    "    GTdd=np.matmul(G.T, dd);\n",
    "    dm = la.solve( GTG, GTdd );\n",
    "    m = m + dm;\n",
    "    s = amp*np.sin(m[3,0]*wa*t);\n",
    "    c = amp*np.cos(m[3,0]*wa*t);\n",
    "    dpre = m[0,0]*dmean + m[1,0]*s + m[2,0]*c;\n",
    "    dd = dobs - dpre;\n",
    "    E = np.matmul(dd.T,dd);\n",
    "    if( (abs(E-Eprev)/E) < 1.0e-6 ):\n",
    "        break\n",
    "\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([0,5000,-30,50]);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (deg C)');\n",
    "plt.plot( t, dobs, 'k-', linewidth=2 );\n",
    "plt.plot( t, dpre, '-', linewidth=2, color=(0.7,0.7,0.7) );\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# account for scalings and compute frequency and period\n",
    "ffinal = (wa*m[3,0])/(2*pi);\n",
    "Tfinal = (2*pi)/(wa*m[3,0]);\n",
    "                 \n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([320,420,0,1]);\n",
    "plt.xlabel('period (days)');\n",
    "plt.ylabel('normalized amplitude spectral density');\n",
    "plt.plot( [Tfinal, Tfinal], [0, 0.75], 'k-', linewidth=3 );\n",
    "\n",
    "# covariance\n",
    "# f = f0 + df   and   T = T0 + dT    and   T0 = 1 / f0\n",
    "# T = 1 / (f0 + df ) = (1/f0) (1+df/f0)^-1\n",
    "#   = (1/f0) (1-df/f0) = (1/f0) - df/(f0^2)\n",
    "# dT = (-1/(f0^2)) df\n",
    "# var(dT) = (-1/(f0^2))^2 var(df)\n",
    "\n",
    "# confidence intervals\n",
    "sigma2d = E/(N-M); # posterior variance\n",
    "GTGI = la.inv(GTG); # inverse of G_transpose times G\n",
    "sigma2m4 = sigma2d*GTGI[3,3]; # variance of mest[3]\n",
    "sigma2f = (wa**2)*sigma2m4/((2*pi)**2);  # variance of frequency f\n",
    "sigma2T = sigma2f / (ffinal**4); # variance of period T\n",
    "sigmaT = sqrt(sigma2T); # standard deviation of period T\n",
    "DT95 = 2*sigmaT; # 95% confidence\n",
    "plt.plot( [Tfinal-DT95, Tfinal+DT95], [0.7, 0.7], 'k-', linewidth=3 );\n",
    "\n",
    "# make evenly sampled time vectr\n",
    "Dt = t[1,0]-t[0,0];                 # sampling interval\n",
    "tstart = t[0,0];                    # start time\n",
    "tend = t[N-1,0];                    # end time\n",
    "N2 = floor((tend-tstart)/Dt)+1;     # N2 is lenth of evenly sampled time series\n",
    "N2 = 2*floor(N2/2); # insure even   # ensure N2 is even\n",
    "t2 = eda_cvec(np.linspace(tstart,tend,N2)); # evenly sampled time vector\n",
    "\n",
    "# standard frequency definitions\n",
    "fny = 1/(2*Dt);                      # Nyquist frequency\n",
    "Nf = floor(N2/2)+1;                  # number of non-negative frequencies\n",
    "Df = fny / floor(N2/2);              # frequency sampling\n",
    "f = eda_cvec(np.linspace(0,fny,Nf)); # vector of non-negative frequencies\n",
    "\n",
    "# interpolate data to even sampling\n",
    "# but note that data gaps filled in with stright lines\n",
    "d2f = ip.interp1d(t.ravel(),dobs.ravel(),kind='linear',fill_value='extrapolate');\n",
    "d2 = eda_cvec(d2f(t2));\n",
    "\n",
    "# demean, take FFT and compute amplitude spectral density\n",
    "d2 = d2-np.mean(d2);                # demean timeseries\n",
    "tmp = Dt*np.fft.fft(d2,axis=0);     # fourier trasnform\n",
    "d2tilda = eda_cvec(tmp[0:Nf,0]);    # keep only non-negative frequencies\n",
    "asd = np.abs(d2tilda);              # amplitude spectral density\n",
    "maxasd = np.max(asd);               # maximum value of a.s.d.\n",
    " \n",
    "plt.plot( np.divide(1,f[1:Nf]), 0.9*asd[1:Nf]/maxasd, 'k-' );\n",
    "plt.plot( np.divide(1,f[1:Nf]), 0.9*asd[1:Nf]/maxasd, 'ko' );\n",
    "\n",
    "plt.show();\n",
    "\n",
    "print('iterations %d' % (ic) );\n",
    "print('posterior error %.3f (deg C)' % (sqrt(sigma2d)) );\n",
    "print('Estimated period %.3f +/- %.3f (95%%) days' % (Tfinal, DT95) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 62784 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18676 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations 4347\n",
      "posterior error 5.619 (deg C)\n",
      "Estimated period 365.630 days\n"
     ]
    }
   ],
   "source": [
    "# edapy_11_05, sinusoidal fit using gradient descent method\n",
    "# fit of a sinusoid of unknown (but approximately annual)\n",
    "# frequency to Black Rock Forest temperature data using the\n",
    "# gradient descent method\n",
    "\n",
    "# load Black Rock Forest temperature data\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t');\n",
    "\n",
    "# clean the data by deleting outliers and zero (= no data) values\n",
    "test1 = (D[:,1]>-35);\n",
    "test2 = (D[:,1]<40);\n",
    "test3 = (D[:,1]!=0.0);\n",
    "r = np.where( test1 & test2 & test3 );\n",
    "k = r[0];\n",
    "tmp = np.shape(k);\n",
    "N = tmp[0];\n",
    "t = D[k,0:1];    # vector of cleaned times, length N\n",
    "dobs = D[k,1:2]; # vector of cleaned observations, length N\n",
    "\n",
    "# some statistics\n",
    "dmean = np.mean(dobs);                   # mean\n",
    "amp = 0.5*(np.max(dobs)-np.min(dobs));   # amplitude\n",
    "\n",
    "# Annual period, frequency, angular frequency\n",
    "Ta = 365.0;\n",
    "fa = 1.0/Ta;\n",
    "wa = 2.0*pi*fa;\n",
    "\n",
    "# nonlinear theory and its derivative\n",
    "# note model parameters are scaled to be order unity\n",
    "# d = m1*dmean + amp*m2*sin(m4 wa t) + amp*m3*cos(m4 wa t)\n",
    "# dd/dm1 =  dmean\n",
    "# dd/dm2 =  amp*sin(m4 wa t)\n",
    "# dd/dm3 =  amp*cos(m4 wa t)\n",
    "# dd/dm4 =  amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)\n",
    "\n",
    "# inital guesses\n",
    "M=4;\n",
    "m = np.zeros((M,1));\n",
    "m[0,0] =  1.0;        # average level is mean\n",
    "m[1,0] =  0.0;        # phase with max in summer\n",
    "m[2,0] = -0.75;       # so only -cos() term\n",
    "m[3,0] =  1.0;        # annual frequency \n",
    "\n",
    "Lm = sqrt(np.matmul(m.T,m)); # solution size\n",
    "\n",
    "s = amp*np.sin(m[3,0]*wa*t);\n",
    "c = amp*np.cos(m[3,0]*wa*t);\n",
    "dpre = dmean*m[0,0] + m[1,0]*s + m[2,0]*c;\n",
    "\n",
    "# plot data and initial prediction\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([0,5000,-30,50]);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (deg C)');\n",
    "plt.plot( t, dobs, 'k-', linewidth=2 );\n",
    "plt.plot( t, dpre, '-', linewidth=2, color=(0.7,0.7,0.7) );\n",
    "\n",
    "# error and its gradient at the trial solution\n",
    "mgo = m;\n",
    "s = amp*np.sin(mgo[3,0]*wa*t);\n",
    "c = amp*np.cos(mgo[3,0]*wa*t);\n",
    "dprego = dmean*mgo[0,0] + mgo[1,0]*s + mgo[2,0]*c;\n",
    "ego = dobs - dprego;\n",
    "Ego = np.matmul(ego.T,ego);\n",
    "Estart = Ego;\n",
    "Go = np.zeros((N,4));\n",
    "Go[:,0] = (dmean*np.ones((N,1))).ravel();\n",
    "Go[:,1] = s.ravel();\n",
    "Go[:,2] = c.ravel();\n",
    "Go[:,3] = (mgo[1,0]*wa*np.multiply(t,c)-mgo[2,0]*wa*np.multiply(t,s)).ravel();\n",
    "dEdmo = -2.0*np.matmul( Go.T, ego );\n",
    "\n",
    "# gradient methid parametres\n",
    "alpha = 0.1;\n",
    "c1 = 0.0001;\n",
    "c2 = 0.9;\n",
    "tau = 0.5;\n",
    "\n",
    "# gradient method\n",
    "Niter = 10000; # maximum number of iterations\n",
    "ic=0;          # iteration counter\n",
    "TOL = 1.0e-5;  # stop when fractional change in model parameters is less than\n",
    "for k in range(Niter):\n",
    "    \n",
    "    # downslope unit vector\n",
    "    v = -dEdmo / sqrt( np.matmul(dEdmo.T,dEdmo) );\n",
    "    \n",
    "    # backstep\n",
    "    NBACK = 100;  # maximum number of back-steps\n",
    "    for kk in range(NBACK):\n",
    "        ic=ic+1;                                        # increment iteration counter\n",
    "        mg = mgo+alpha*v;                               # update model\n",
    "        s = amp*np.sin(mg[3,0]*wa*t);                   # sines\n",
    "        c = amp*np.cos(mg[3,0]*wa*t);                   # cosines\n",
    "        dpreg = dmean*mg[0,0] + mg[1,0]*s + mg[2,0]*c;  # predicted data\n",
    "        eg = dobs - dpreg;                              # individual errors\n",
    "        Eg = np.matmul(eg.T,eg);                        # total error\n",
    "        G = np.zeros((N,4));                            # set up data kernel\n",
    "        G[:,0] = (dmean*np.ones((N,1))).ravel();\n",
    "        G[:,1] = s.ravel(); \n",
    "        G[:,2] = c.ravel();\n",
    "        G[:,3] = (mg[1,0]*wa*np.multiply(t,c)-mg[2,0]*wa*np.multiply(t,s)).ravel();\n",
    "        dEdm = -2.0*np.matmul(G.T,eg);                  # grdient of error\n",
    "        \n",
    "        RULE = 1;  # Rule=1, the simplest rule: error went down\n",
    "                   # RULE=2 is Armijo's rule\n",
    "                   # RULE=3 is RULE2 + Wolfe condition\n",
    "        if( RULE==1 ):\n",
    "            if( Eg<=Ego ):\n",
    "                break;\n",
    "        elif( RULE==2 ):\n",
    "            if( (Eg<=(Ego + c1*alpha*np.matmul(v.T,dEdmo))) ):\n",
    "                break;\n",
    "        elif( RULE==3 ):\n",
    "            test1 = (Eg<=(Ego + c1*alpha*np.matmul(v.T,dEdmo)));\n",
    "            test2 = ((-np.matmul(v.T,dEdm)) <= (-c2*np.matmul(v.T,dEdmo)));\n",
    "            if( test1 and test2 ):\n",
    "                break;\n",
    "        \n",
    "        # if the test failed, decrease the step size\n",
    "        alpha = tau*alpha;\n",
    "\n",
    "    # change in solution\n",
    "    Dmg = sqrt( np.matmul((mg-mgo).T,(mg-mgo)) );\n",
    "    \n",
    "    # update\n",
    "    mgo=mg;\n",
    "    dprego = dpreg;\n",
    "    ego = eg;\n",
    "    Ego = Eg;\n",
    "    Go = G;\n",
    "    dEdmo = dEdm;\n",
    "\n",
    "    if( (Dmg/Lm) < TOL ):\n",
    "        break;\n",
    "\n",
    "m = mg;         # estimated model is from last iteration\n",
    "dpre = dpreg;   # predicted datais from last iteration\n",
    "E = Ego;        # error from last iteration\n",
    "sigma2d = E/(N-M); # posterior variance\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([0,5000,-30,50]);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (deg C)');\n",
    "plt.plot( t, dobs, 'k-', linewidth=2 );\n",
    "plt.plot( t, dpre, '-', linewidth=2, color=(0.7,0.7,0.7) );\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# account for scalings and compute frequency and period\n",
    "ffinal = (wa*m[3,0])/(2*pi);\n",
    "Tfinal = (2*pi)/(wa*m[3,0]);\n",
    "                 \n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([320,420,0,1]);\n",
    "plt.xlabel('period (days)');\n",
    "plt.ylabel('normalized amplitude spectral density');\n",
    "plt.plot( [Tfinal, Tfinal], [0, 0.75], 'k-', linewidth=3 );\n",
    "\n",
    "# standard time & frequency definitions\n",
    "Dt = t[1,0]-t[0,0];\n",
    "tstart = t[0,0];\n",
    "tend = t[N-1,0];\n",
    "N2 = floor((tend-tstart)/Dt)+1;\n",
    "N2 = 2*floor(N2/2); # insure even\n",
    "t2 = np.zeros((N2,1));\n",
    "t2[:,0] = np.linspace(tstart,tend,N2);\n",
    "fny = 1/(2*Dt);\n",
    "Nf = floor(N2/2)+1;\n",
    "Df = fny / floor(N2/2);\n",
    "f = np.zeros((Nf,1));\n",
    "f[:,0] = np.linspace(0,fny,Nf);\n",
    "\n",
    "# interpolate data to even sampling\n",
    "# but note that data gaps filled in with stright lines\n",
    "d2f = ip.interp1d(t.ravel(),dobs.ravel(),kind='linear',fill_value='extrapolate');\n",
    "d2 = eda_cvec( d2f(t2) );\n",
    "\n",
    "# demean, take FFT and compute amplitude spectral density\n",
    "d2 = d2-np.mean(d2);               # demean\n",
    "tmp = Dt*np.fft.fft(d2,axis=0);    # fourier transform\n",
    "d2tilda = eda_cvec(tmp[0:Nf,0]);   # keep only non-negative frequencies\n",
    "asd = np.abs(d2tilda);             # amplitude spectral density\n",
    "maxasd = np.max(asd);              # maximum a.s.d.\n",
    "\n",
    "plt.plot( np.divide(1,f[1:Nf]), 0.9*asd[1:Nf]/maxasd, 'k-' );\n",
    "plt.plot( np.divide(1,f[1:Nf]), 0.9*asd[1:Nf]/maxasd, 'ko' );\n",
    "\n",
    "plt.show();\n",
    "\n",
    "print('iterations %d' % (ic) );\n",
    "print('posterior error %.3f (deg C)' % (sqrt(sigma2d)) );\n",
    "print('Estimated period %.3f days' % (Tfinal) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iterations: 9999\n",
      "mean: 9.961 amp: -1.498 -13.879\n",
      "Estimated period 367.022 error ratio 0.292\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 42160 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_06\n",
    "# stochastic gradient method, full error every\n",
    "# 100 iterations\n",
    "# non-linear least squares fit of a sinusoid\n",
    "# of unknown frequency, using the stochastic\n",
    "# gradient method.  The data are subsampled on\n",
    "# each iteration.  A fixed step size alpha is\n",
    "# used to adcance the solution anti-parallel to the\n",
    "# gradient. The solution is considered good enough\n",
    "# when its posterior variance is reduced to a\n",
    "# specififed value.\n",
    "\n",
    "# load Black Rock Forest temperature data\n",
    "D = np.genfromtxt('../Data/brf_temp.txt', delimiter='\\t')\n",
    "\n",
    "# clean the data by deleting outliers and zero (= no data) values\n",
    "test1 = (D[:,1]>-35);\n",
    "test2 = (D[:,1]<40);\n",
    "test3 = (D[:,1]!=0.0);\n",
    "r = np.where( test1 & test2 & test3);\n",
    "k = r[0];\n",
    "tmp = np.shape(k);\n",
    "N = tmp[0];\n",
    "t = D[k,0:1];    # cleaned time vector, length N\n",
    "dobs = D[k,1:2]; # cleaned observation vector, length N\n",
    "\n",
    "# subsampling information\n",
    "decimation_factor = 1000;\n",
    "Ns = floor(N/decimation_factor);  # number of points in subsampled data\n",
    "\n",
    "# some statistics\n",
    "dmean = np.mean(dobs);                   # mean\n",
    "amp = 0.5*(np.max(dobs)-np.min(dobs));   # amplitude\n",
    "\n",
    "# the solution considered good enough when\n",
    "# its posterior variance is smaller than\n",
    "low_enough_variance = 10.0;\n",
    "\n",
    "# Annual period, frequency, angular frequency\n",
    "Ta = 365.0;\n",
    "fa = 1.0/Ta;\n",
    "wa = 2.0*pi*fa;\n",
    "\n",
    "# the model parameters m in this problem have been\n",
    "# normalized so that we expect them to be of order\n",
    "# 1.  The DC offset is parameterized as a fraction\n",
    "# of the mean of the data, dmean*m0. The amplitudes\n",
    "# of the mean of the data, dmean*m0. The amplitudes\n",
    "# of the overall amplitude of the data, m1*amp and\n",
    "#  m2*amp. The frequency is parameterized as a\n",
    "# fraction of the annual frequecy, m3*wa.\n",
    "\n",
    "# The resulting nonlinear theory and its derivative\n",
    "# d = m0*dmean + m1*amp*sin(m3 wa t) + m2*amp*cos(m3 wa t)\n",
    "# dd/dm0 =  dmean\n",
    "# dd/dm1 =  amp*sin(m4 wa t)\n",
    "# dd/dm2 =  amp*cos(m4 wa t)\n",
    "# dd/dm3 =  amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)\n",
    "\n",
    "# inital guesses\n",
    "M=4;\n",
    "m = np.zeros((M,1));\n",
    "m[0,0] =  1.0;        # average level is mean\n",
    "m[1,0] =  0.0;        # phase with max in summer\n",
    "m[2,0] = -0.75;       # so only -cos() term\n",
    "m[3,0] =  1.0;        # annual frequency \n",
    "Lm = sqrt(np.matmul(m.T,m));\n",
    "\n",
    "# evaluate initial guess\n",
    "s = amp*np.sin(m[3,0]*wa*t);\n",
    "c = amp*np.cos(m[3,0]*wa*t);\n",
    "dpre = dmean*m[0,0] + m[1,0]*s + m[2,0]*c;\n",
    "\n",
    "# plot data and initial prediction\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([0,5000,-30,50]);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (deg C)');\n",
    "plt.plot( t, dobs, 'k-', linewidth=2 );\n",
    "plt.plot( t, dpre, '-', linewidth=2, color=(0.7,0.7,0.7) );\n",
    "\n",
    "# trial solution\n",
    "mg = m;\n",
    "mg0 = m;\n",
    "\n",
    "# the step size alpha is fixed\n",
    "alpha = 1e-4;\n",
    "\n",
    "myrange = 1;\n",
    "Niter = 10000; # number of iterations\n",
    "TOL = 1e-5;  # terminate iteration when change in model parameters drops below\n",
    "\n",
    "for k in range(Niter):\n",
    "    # subsample data\n",
    "    ns = np.random.randint(0,N,(Ns,1));  # Ns random integers in range 0 to N-1\n",
    "    ts = eda_cvec( t[ns,0] );            # subsampled time vector\n",
    "    dobss = eda_cvec( dobs[ns,0] );      # subsampled observations\n",
    "    \n",
    "    # compute predicted data, the error and the\n",
    "    # gradient of the error of the subsampled data\n",
    "    s = amp*np.sin(mg[3,0]*wa*ts);\n",
    "    c = amp*np.cos(mg[3,0]*wa*ts);\n",
    "    dpreg = dmean*mg[0,0] + mg[1,0]*s + mg[2,0]*c; # predicted daya\n",
    "    eg = dobss - dpreg;                            # individual errors\n",
    "    Eg = np.matmul(eg.T,eg);                       # total error\n",
    "    if( k==1 ):                                    # save starting error\n",
    "        Estart = Eg;\n",
    "    G = np.zeros((Ns,4));                          # set up data kernel\n",
    "    G[:,0] = (dmean*np.ones((Ns,1))).ravel();\n",
    "    G[:,1] = s.ravel();\n",
    "    G[:,2] = c.ravel();\n",
    "    G[:,3] = (m[1,0]*wa*np.multiply(ts,c)-m[2,0]*wa*np.multiply(ts,s)).ravel();                \n",
    "    dEdm = -2*np.matmul(G.T,eg);                   # gradient of error\n",
    "                           \n",
    "    # posterior variance of data\n",
    "    sigma2d =  Eg/(Ns-M);\n",
    "    \n",
    "    # compute just the error for the full dataset\n",
    "    # but only every 100 iterations\n",
    "    if( mod(k-1, 100) == 0 ):\n",
    "        sfull = amp*np.sin(mg[3,0]*wa*t); # sines\n",
    "        cfull = amp*np.cos(mg[3,0]*wa*t); # cosines\n",
    "        dpregfull = dmean*mg[0,0] + mg[1,0]*sfull + mg[2,0]*cfull; # predicte data\n",
    "        egfull = dobs - dpregfull; # individual errors\n",
    "        Egfull = np.matmul(egfull.T,egfull); # total error\n",
    "        # posterior variance of data\n",
    "        sigma2dfull =  Egfull/(N-M);\n",
    "        if( k==1 ):\n",
    "            Estartfull = Egfull;\n",
    "        # stop when the veriance is low enough\n",
    "        if( sigma2dfull <= low_enough_variance ):\n",
    "            break;\n",
    "    \n",
    "    # downslope unit vector\n",
    "    v = -dEdm / sqrt( np.matmul(dEdm.T,dEdm));\n",
    "    # update solution\n",
    "    mg = mg+alpha*v;\n",
    "    \n",
    "    # change in solution\n",
    "    Dmg = sqrt( np.matmul((mg-mgo).T,(mg-mgo)) );\n",
    "    if( (Dmg/Lm) < TOL ):\n",
    "        break;\n",
    "\n",
    "# estimated solution\n",
    "m = mg;                       # estimated model parametres from the last iteration\n",
    "s = amp*np.sin(m[3,0]*wa*t);  # sines\n",
    "c = amp*np.cos(m[3,0]*wa*t);  # cosines\n",
    "dpre = dmean*m[0,0] + m[1,0]*s + m[2,0]*c; # predicted data\n",
    "\n",
    "plt.plot(t,dpre,'-',linewidth=2,color=(0.9,0.9,0.9));\n",
    "\n",
    "ffinal = wa*m[3,0]/(2*pi);    # estimated freqiency\n",
    "Tfinal = (2*pi)/(wa*m[3,0]);  # estimated perior\n",
    "\n",
    "print('Iterations: %d' % (k) );\n",
    "print('mean: %.3f amp: %.3f %.3f' % (m[0,0]*dmean, m[1,0]*amp, m[2,0]*amp) );\n",
    "print('Estimated period %.3f error ratio %.3f' % (Tfinal, Egfull/Estartfull) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true:     m 0.5142 0.7524    E 100.490\n",
      "exact:    m 0.5100 0.7600    E 108.853 duration 7276 ms\n",
      "lookup:   m 0.5000 0.7700    E 120.737 duration 1204 ms\n",
      "speed improvement 83.45 percent\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 44864 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_07\n",
    "# example of a grid search that uses a table\n",
    "# lookup of the function\n",
    "# f(z) = sin(pi*sin(pi*z*sin(pi*(1-z^2))))\n",
    "\n",
    "# true values of the model parameters\n",
    "M=2;\n",
    "mtrue = np.zeros((M,1));\n",
    "mtrue[0,0] = 0.5142;\n",
    "mtrue[1,0] = 0.7524\n",
    "\n",
    "# standard setup of time axis, t\n",
    "N = 10001;   # length of time series\n",
    "tmin = 0.0;  # start time\n",
    "tmax = 1.0;  # end time\n",
    "Dt = (tmax-tmin)/(N-1); # sampling interval\n",
    "t = eda_cvec( np.linspace(tmin, tmax, N) ); # timve vector\n",
    "\n",
    "# subsamples indices for plotting purposes\n",
    "Nsub = 20;\n",
    "isub = range(0,N,Nsub);\n",
    "\n",
    "# the true data obey the formula\n",
    "# dtrue = f(z) with z = m1 + m2*t\n",
    "z = mtrue[0,0]+mtrue[1,0]*t;\n",
    "dtrue = np.sin( pi * np.sin( pi*np.multiply(z,np.sin(pi*(1-np.power(z,2)))))); # true data\n",
    "# synthetic estimated data is true data plus random noise\n",
    "sigmad = 0.1;\n",
    "dobs = dtrue + np.random.normal(0.0,sigmad,(N,1)); # observed = signal + noise\n",
    "etrue = dobs - dtrue; # individual errors\n",
    "Etrue = np.matmul( etrue.T, etrue ); # total error\n",
    "print('true:     m %.4f %.4f    E %.3f' % (mtrue[0,0], mtrue[1,0], Etrue) );\n",
    "\n",
    "# plot the true and observed data\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([tmin,tmax,-1.2,1.2]);\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.plot(t[isub,0],dobs[isub,0],'k.',linewidth=2);\n",
    "plt.plot(t[isub,0],dtrue[isub,0],'-',linewidth=7,color=(0.9,0.9,0.9));\n",
    "\n",
    "# range of m1 for grid search\n",
    "M1=101;       # m1-dimension of grid\n",
    "m1min = 0.0;  # minimum m1\n",
    "m1max = 1.0;  # maximum m1\n",
    "Dm1 = (m1max-m1min)/(M1-1); # sampling interval\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,M1) ); # vector of m1's\n",
    "\n",
    "# range of m2 for grid search\n",
    "M2=101;       # m2-dimension of grid\n",
    "m2min = 0.0;  # minimum m2\n",
    "m2max = 1.0;  # maximum m2 \n",
    "Dm2 = (m2max-m2min)/(M1-1); # sampling interval\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,M2) ); # vector of m2's\n",
    "\n",
    "tic = timeit.default_timer();\n",
    "\n",
    "# standard grid search that computes f(z) exactly\n",
    "# note use of tic/toc pair to time exectiion\n",
    "tic;\n",
    "first = 1;\n",
    "for m1trial in m1:\n",
    "    for m2trial in m2:\n",
    "        z = m1trial+m2trial*t; # evaluate z\n",
    "        dtrial = np.sin( pi * np.sin( pi*np.multiply(z,np.sin(pi*(1-np.power(z,2)))))); # evaluate d(z)\n",
    "        etrial = dobs - dtrial; # individual errors\n",
    "        Etrial = np.matmul(etrial.T,etrial); # total error\n",
    "        # keep solun with smallest error\n",
    "        if( first==1 ):\n",
    "            Ebest1 = Etrial;\n",
    "            m1best1 = m1trial;\n",
    "            m2best1 = m2trial;\n",
    "            first=0;\n",
    "        elif (Etrial < Ebest1 ):\n",
    "            Ebest1 = Etrial;\n",
    "            m1best1 = m1trial;\n",
    "            m2best1 = m2trial;\n",
    "toc = timeit.default_timer();\n",
    "duration1=1000*(toc-tic);\n",
    "\n",
    "# print and plot results\n",
    "print('exact:    m %.4f %.4f    E %.3f duration %.0f ms' % (m1best1, m2best1, Ebest1, duration1) );\n",
    "z = m1best1+m2best1*t;\n",
    "dest = np.sin( pi * np.sin( pi*np.multiply(z,np.sin(pi*(1-np.power(z,2))))));\n",
    "plt.plot(t[isub,0],dest[isub,0],':',linewidth=4,color=(0.6,0.6,0.6));\n",
    "\n",
    "# grid search with table lookup of the funtion f(z)\n",
    "tic=timeit.default_timer();\n",
    "\n",
    "# compute the table\n",
    "K = 1001;            # size of loolup table\n",
    "zmintable = -0.1;    # minimum z value\n",
    "zmaxtable =  2.1;    # maximum z value\n",
    "Dztable = (zmaxtable-zmintable)/(K-1); # sampling interval in z\n",
    "ztable = eda_cvec(np.linspace(zmintable,zmaxtable,K)); # vector of z's\n",
    "# corresponding vector of d's\n",
    "dtable = ((np.sin( pi * np.sin( pi*np.multiply(ztable,np.sin(pi*(1-np.power(ztable,2))))))));\n",
    "\n",
    "# grid search with table lookup\n",
    "first = 1;\n",
    "for m1trial in m1:\n",
    "    for m2trial in m2:\n",
    "        z = m1trial+m2trial*t; # evaluate z\n",
    "        nn = np.floor((z-zmintable)/Dztable)+1; # lookup-table index of that z\n",
    "        n = nn.astype(int);\n",
    "        dtrial = dtable[n,0]; # look up d(z) in table\n",
    "        etrial = dobs - dtrial; # indivisual errors\n",
    "        Etrial = np.matmul(etrial.T,etrial); # total error\n",
    "        # keep solun with smallest error\n",
    "        if( first==1 ):\n",
    "            Ebest2 = Etrial;\n",
    "            m1best2 = m1trial;\n",
    "            m2best2 = m2trial;\n",
    "            first=0;\n",
    "        elif (Etrial < Ebest2 ):\n",
    "            Ebest2 = Etrial;\n",
    "            m1best2 = m1trial;\n",
    "            m2best2 = m2trial;\n",
    "\n",
    "toc = timeit.default_timer();\n",
    "duration2=1000*(toc-tic);\n",
    "\n",
    "# print and plot estimated solution\n",
    "print('lookup:   m %.4f %.4f    E %.3f duration %.0f ms' % (m1best2, m2best2, Ebest2, duration2) );\n",
    "print('speed improvement %.2f percent' % (100*(duration1-duration2)/duration1) );\n",
    "z = m1best2+m2best2*t; # estimated z's\n",
    "dest = np.sin( pi * np.sin( pi*np.multiply(z,np.sin(pi*(1-np.power(z,2)))))); # estimated data\n",
    "plt.plot(t[isub,0],dest[isub,0],'k-',linewidth=2);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20492 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_08\n",
    "# plot sigma(x) = (1+exp(-wx))^-1\n",
    "# for a variety of w's.  Note that\n",
    "# dsigma/dx = w exp(-w) (1+exp(-wx))**(-2)\n",
    "# so the slope at x=0 is s=w/4\n",
    "\n",
    "Ns = 4;\n",
    "s = eda_cvec( [100, 3, 1, 0.1] );  # slopes\n",
    "\n",
    "# standard x axis setup\n",
    "Nx = 501;  # number of x's\n",
    "xmin = -2; # starting x\n",
    "xmax = 2;  # ending x\n",
    "x = eda_cvec( np.linspace(xmin,xmax,Nx) ); # vector of x's\n",
    "\n",
    "# plot the true and observed data\n",
    "fig1 = plt.figure(1,figsize=(8,8));\n",
    "\n",
    "for i in range(Ns):\n",
    "    w=4*s[i,0];                    # width\n",
    "    sigma = eda_sigma(w,x);        # evaluate sigma function\n",
    "    # make subplot\n",
    "    ax1 = plt.subplot(Ns, 1, i+1);\n",
    "    plt.axis([xmin,xmax,-1.5,1.5]);\n",
    "    plt.xlabel('x');\n",
    "    plt.ylabel('sigma(x)');\n",
    "    plt.plot(x,sigma,'k-',LineWidth=3);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4832 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 7780 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 4404 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 12844 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_09\n",
    "#\n",
    "# Use a neural net to approximate a step function\n",
    "\n",
    "# Step function definitions \n",
    "Xc = 2.0; # center of the step function\n",
    "h = 1.0; # height of the step function\n",
    "slope = 250.0; # slope of the step function\n",
    "\n",
    "# Neural Net definitions\n",
    "Lmax = 3; # number of layers\n",
    "Nmax = 2; # maximum number of neurons in any layer\n",
    "\n",
    "# initialize net to zeros\n",
    "N, w, b, a = eda_net_init(Lmax,Nmax);\n",
    "\n",
    "# set net values\n",
    "N[:,0] = [1, 1, 1]; # one neuron in each layer\n",
    "# weights\n",
    "w[0,0,1] = 4*slope; # (neuron 0 in layer 1) fr (neuron 0 in previous layer)\n",
    "w[0,0,2] = h;       # (neuron 0 in layer 2) fr (neuron 0 in previous layer)\n",
    "b[0,1] = -4*slope*Xc; # bias of (neuron 0 in layer 1)\n",
    "\n",
    "# plot network\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig1,ax1);\n",
    "plt.show();\n",
    "\n",
    "# standard setup of x-axis\n",
    "Nx= 100;\n",
    "xmin = 0.0;\n",
    "xmax = 4.0;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "\n",
    "# evaluate network for each x to get m(x)\n",
    "m = np.zeros((Nx,1));\n",
    "for kx in range(Nx):  # loop over x's\n",
    "    a[0:N[0,0],0] = x[kx,0]; # activity of (neuron 0 in layer 0) = x\n",
    "    a, z = eda_net_eval(N,w,b,a); # evaluate network\n",
    "    m[kx,0] = a[0:N[2,0],2]; # m(x) = activity of (neuron 0 in last layer)\n",
    "    \n",
    "# plot m(x)\n",
    "fig2 = plt.figure(2,figsize=(12,5));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,-1.5,1.5]);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('m(x)');\n",
    "plt.plot(x,m,'k-',LineWidth=2);\n",
    "\n",
    "# now do a different slope slope\n",
    "slope = 1.0;\n",
    "\n",
    "Lmax = 3;\n",
    "Nmax = 2;\n",
    "\n",
    "N, w, b, a = eda_net_init(Lmax,Nmax);\n",
    "N[:,0] = [1, 1, 1];\n",
    "w[0,0,1] = 4*slope;\n",
    "w[0,0,2] = h;\n",
    "b[0,1] = -4*slope*Xc;\n",
    "\n",
    "fig3 = plt.figure(1,figsize=(16,6));\n",
    "ax3 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig3,ax3);\n",
    "plt.show();\n",
    "\n",
    "Nx= 100;\n",
    "xmin = 0;\n",
    "xmax = 4;\n",
    "x = eda_cvec( np.linspace(xmin,xmax,Nx) );\n",
    "\n",
    "m = np.zeros((Nx,1));\n",
    "for kx in range(Nx):\n",
    "    a[0:N[0,0],0] = x[kx,0];\n",
    "    a, z = eda_net_eval(N,w,b,a);\n",
    "    m[kx,0] = a[0:N[2,0],2];\n",
    "    \n",
    "fig4 = plt.figure(2,figsize=(12,5));\n",
    "ax4 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,-1.5,1.5]);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('m(x)');\n",
    "plt.plot(x,m,'k-',LineWidth=2);\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14276 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 12156 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13308 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20584 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_10, neural net for 1D tower function\n",
    "\n",
    "# Tower function definitions\n",
    "Nc=1;                      # one tower\n",
    "Xc = 2.0*np.ones((Nc,1));  # center of the tower\n",
    "Dx = 1.0*np.ones((Nc,1));  # width of the tower\n",
    "h = 1.5*np.ones((Nc,1));   # the height of the tower function\n",
    "slope = 250.0; # slope of all the towers\n",
    "\n",
    "# neural net setup\n",
    "Lmax = 3; # number of layers\n",
    "Nmax = 2; # maximum number of neurons in any layer\n",
    "# initialize new to zeros\n",
    "N,w,b,a = eda_net_init(Lmax,Nmax);\n",
    "# initialize new to the tower\n",
    "N,w,b=eda_net_1dtower(Nc,slope,Xc,Dx,h,N,w,b);\n",
    "\n",
    "# plot the net\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig1,ax1);\n",
    "plt.show();\n",
    "\n",
    "# standars setup of an x-axis\n",
    "Nx= 100;\n",
    "xmin = 0.0;\n",
    "xmax = 4.0;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "\n",
    "# evaluate the network for each x to get m(x)\n",
    "m = np.zeros((Nx,1));\n",
    "for kx in range(Nx):  # loop over x\n",
    "    a[0:N[0,0],0] = x[kx,0]; # activity of (neoron 0 in layer 0) = x\n",
    "    a, z = eda_net_eval(N,w,b,a); # evaluate network\n",
    "    m[kx,0] = a[0:N[2,0],2]; # m(x) is activity of (neoron 0 in last layer)\n",
    "    \n",
    "# plot m(x)\n",
    "fig2 = plt.figure(2,figsize=(12,5));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,-1.0,2.0]);\n",
    "plt.xlabel('a[0,0]');\n",
    "plt.ylabel('a[2,0]');\n",
    "plt.title('tower function; slope = 999');\n",
    "plt.plot(x,m,'k-',LineWidth=2);\n",
    "plt.show();\n",
    "\n",
    "# repeat for different slope\n",
    "slope = 1.0;\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "[N,w,b]=eda_net_1dtower(Nc,slope,Xc,Dx,h,N,w,b);\n",
    "\n",
    "# plot the net\n",
    "fig3 = plt.figure(3,figsize=(16,6));\n",
    "ax3 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig3,ax3);\n",
    "plt.show();\n",
    "\n",
    "m = np.zeros((Nx,1));\n",
    "for kx in range(Nx):\n",
    "    a[0:N[0,0],0] = x[kx,0];\n",
    "    a, z = eda_net_eval(N,w,b,a);\n",
    "    m[kx,0] = a[0:N[2,0],2];\n",
    "    \n",
    "fig4 = plt.figure(4,figsize=(12,5));\n",
    "ax4 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,-1.0,2.0]);\n",
    "plt.xlabel('a[0,0]');\n",
    "plt.ylabel('a[2,0]');\n",
    "plt.title('tower function; slope = 999');\n",
    "plt.plot(x,m,'k-',LineWidth=2);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 71484 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20332 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edaps_11_11, approximate function with several towers\n",
    "\n",
    "# standard setup of an axis, xc, of tower cenetrs\n",
    "Nc= 6; # the number of centers in the approximation\n",
    "xcmin = 0.0; # maximum value of the centers\n",
    "xcmax = 1.0; # minimum value of the centers\n",
    "xc = eda_cvec(np.linspace(xcmax,xcmin,Nc)); # the centers\n",
    "Dx = np.ones((Nc,1))*(xcmax-xcmin)/(Nc-1); # widths\n",
    "# yc: the true value of the polynomial function\n",
    "y = 0.2 + 1.75*np.power(xc,2); # heights\n",
    "h = y; # height of the step functions\n",
    "slope = 250.0; # slopes all the same\n",
    "\n",
    "# ---Neural Net definitions---\n",
    "Lmax = 3; # number of layers;\n",
    "Nmax = Nc * 2; # maximum number of neurons in any layer\n",
    "# initialize net to zeros\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "# initialize net to towers\n",
    "[N,w,b] = eda_net_1dtower( Nc, slope, xc, Dx, h, N, w, b );\n",
    "\n",
    "# plot network\n",
    "fig1 = plt.figure(1,figsize=(14,10));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig1,ax1);\n",
    "plt.show();\n",
    "\n",
    "# standard setuo of x-axis\n",
    "Nx= 100;\n",
    "xmin = xcmin - 0.1;\n",
    "xmax = xcmax + 0.1;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "\n",
    "# evaluate the network for each x to get m(x)\n",
    "m = np.zeros((Nx,1));\n",
    "for kx in range(Nx):\n",
    "    a[0:N[0,0],0] = x[kx,0]; # activity of (neuron 0 in layer 0) = x\n",
    "    a, z = eda_net_eval(N,w,b,a); # evaluate network\n",
    "    m[kx,0] = a[0:N[2,0],2]; # m(x) = activity of (neuron 0 in last layer)\n",
    "    \n",
    "# true solution\n",
    "mtrue = 0.2 + 1.75*np.power(x,2);\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(12,5));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,-1.0,2.0]);\n",
    "plt.xlabel('a[0,0]');\n",
    "plt.ylabel('a[2,0]');\n",
    "plt.title('tower function; slope = 999');\n",
    "plt.plot(x,m,'k-',LineWidth=4);\n",
    "plt.plot(x,mtrue,'-',LineWidth=2,color=(0.8,0.8,0.8));\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24448 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAKq0lEQVR4nO3c32vld53H8dd7JhiZFCWZelG1iJZmxFZxOyo6gzgFBbVFwYtliwVFQQSHXsjCLuz+BSsIXoh644VQEfSi/sJ6oRWpIDMVJ+MPkPUHaLGgF22VOcmYpJ+9yLGkxU7T02NO3t3HAwL5fs/3fPJO+FycZ77JqTFGAAAA6OfYogcAAABgNoIOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoKmleS30yCOPjHmtdVgmk0kuXLiQjY2N7O7uPnl+jKP9rVRV1tfXc/bs2ayuri56nLm54YYbatEzPJuO+5yjpcM+BwD6cIcOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANLU0r4Umk8m8ljo0k8kkk8kkV65cye7u7qLHObCqenLu5eXlRY/DEffEE09kc3Mzm5ubGWMsepznZHl5OSsrKzl+/PiiRwEAOJLmFnQXLlyY11KH5sqVK7l48WIuXbrULugee+yxLC0tZW1tbdHjzM1NN9206BFekDY3N/Pggw/m4sWL2d7eXvQ4B1ZVueWWW3L77bfn+uuvX/Q4AABH0tyCbmNjY15LHZrJZJKNjY1cvnw5Ozs7ix7nwKoqW1tbWVlZeUEF3V133bXoEV6Qtra28tBDD+Xee+/N1atXFz3OgVVV7rjjjpw+fVrQAQA8g7kFXac7XH+3s7PzlI9Otre3s7u72/LnzuEaY2R7eztXr17N1tbWosc5sKrK9vZ2uz8TBQA4TN4UBQAAoClBBwAA0NTc/uTSn0Uthp87AAD8/+UOHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNLc1roaqa11KHpqqe/Ohk/9zdZmcxOu6XY8eO5dgxv3MCALiWuQXdqVOn5rXUoZlMJnn88ceztbWVnZ2dRY9zYFWVU6dOZX19PWtra4sehyNueXk5t956a+68885W+zxJTp8+neuuu27RYwAAHFlzC7ozZ87Ma6lDM5lMsrS0lBMnTmR3d3fR4xxYVWV9fT1nzpwRdDyrEydO5Ny5c7ntttsyxlj0OM/JysqKPQ4AcA1zC7rV1dV5LXVolpeXs7q6mrW1tacE3VF/0VtVT87d8efO4Tp+/HhOnjyZkydPLnoUAADmzD+oAAAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKApQQcAANCUoAMAAGhK0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmqoxxqJnAAAAYAbu0AEAADQl6AAAAJoSdAAAAE0JOgAAgKYEHQAAQFOCDgAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUE3J7Xn+1X1kunxu6vqV1X166r6z2d4zier6pdVdbmqvldVr5qef1lV3X+Y88NBzbjXP15VP6uqS1X1YFW9bnreXgcAeB4E3fy8N8nGGOMvVXU8yWeTvCfJ65Lc9fcXsE/z0yRvGmO8IcnXkvxPkowx/pzkkao6ezijw3Myy17/8hjj9WOMN2Zvn386sdcBAJ4vQTeD6d2GS9OP31XVA0k+mOTr00vekuTXY4zfjjH+luQrSd7/9HXGGA+MMSbTwx8neeW+h++brgkLU1X3VdVPquoXVfWx6elZ9vpf9h2uJBn7ju11AIAZCboZjDE+P73T8OYkD2fvbsPZJD+ZXvKKJH/Y95SHp+eu5aNJvrPv+KEkb5/LwDC7j4wxTid5U5J7qupkZtzrVfWJqvpN9u7Q3bPvIXsdAGBGgu75+UyS748xvplkbYzx1+n5+gfXjn9wbu/iqruz94L5U/tO/ynJy+c1KMzonqrayN4d5BuT3JwZ9/oY47NjjJuS/EeS/973kL0OADCjpUUP0FVVfTjJq5Kcn57aqapjY4wnsneX4sZ9l78yyR+fYZ13JvmvJO8YY1zd99CLk2zOe244qKo6l+SdSd42xphU1Q+yty9n2uv7fCXJ5/Yd2+sAADNyh24GVXU6yb8nuXv6ojZJfpXkNdPPLya5uapeXVUvSvJvSb4xfe75qjo//fxfknwhyfvGGH962pdZT/Lzf+53Atf00iSPTmPutUneOj0/y16/ed+6dyT5333H9joAwIzcoZvN+SRrSR6oqmTvf4C+neRc9t4gYmf6Qva7SY4n+eIY4xfT5742yY+mn38qyXVJvjpd5/djjPdNH7t9uiYsyv1JPl5Vl7MXcT+enp9lr5+f3o3eTvJokg/t+zr2OgDAjGqMZ/zXLp6DqrohyZfGGO96luu+leQD03cEvNZ1P0zy/jHGo3McE543ex0A4OgQdHNUVf+a5P6nvUX7LOu8LMnZMcZ985kM5steBwA4GgQdAABAU94UBQAAoClBBwAA0JSgAwAAaErQAQAANCXoAAAAmhJ0AAAATQk6AACApgQdAABAU4IOAACgKUEHAADQlKADAABoStABAAA0JegAAACaEnQAAABNCToAAICmBB0AAEBTgg4AAKCp/wNykvm7FHZbggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_12, newral net for 2D tower\n",
    "\n",
    "\n",
    "# set up towers (of which there are actually only 1)\n",
    "Nc = 1; # number of centers\n",
    "# center of the towers\n",
    "Xc = np.zeros((Nc,1));  # x position of center\n",
    "Xc[:,0] = [1.5];\n",
    "Yc = np.zeros((Nc,1));  # y position of center\n",
    "Yc[:,0] = [1.5];\n",
    "# widths of the towers\n",
    "Dx = np.zeros((Nc,1));\n",
    "Dx[:,0] = [1];\n",
    "Dy = np.zeros((Nc,1));\n",
    "Dy[:,0] = [1];\n",
    "# h: the height of the towers\n",
    "h = np.zeros((Nc,1));\n",
    "h[:,0] = [1];\n",
    "# slope: slope of all the towers \n",
    "slope = 250;\n",
    "\n",
    "# ---Neural Net definitions---\n",
    "Lmax = 4; # number of layers;\n",
    "Nmax = Nc * 4; # maximum number of neurons in any layer\n",
    "# initialize net to zeros\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "# initialize to towers\n",
    "N, w, b = eda_net_2dtower(Nc,slope,Xc,Yc,Dx,Dy,h,N,w,b);\n",
    "\n",
    "# plot network\n",
    "fig1 = plt.figure(1,figsize=(14,10));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig1,ax1);\n",
    "plt.show();\n",
    "\n",
    "# evaluate net on this (x,y) grid\n",
    "Nx = 50;\n",
    "Ny = 50;\n",
    "xmin = 0.0;\n",
    "ymin = 0.0;\n",
    "xmax = 3.1;\n",
    "ymax = 3.1;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "y = eda_cvec(np.linspace(ymin,ymax,Ny));\n",
    "m = np.zeros((Nx,Ny));\n",
    "zprev = np.zeros((Nx,Ny));\n",
    "\n",
    "# evauate the network for each x to get m(x)\n",
    "# and also save the intermediate output z02\n",
    "for kx in range(Nx):\n",
    "    for ky in range(Ny):\n",
    "        a[0:N[0,0],0] = [x[kx,0], x[ky,0] ]; \n",
    "        a, z = eda_net_eval(N,w,b,a);\n",
    "        zprev[kx,ky] = z[0,2];\n",
    "        m[kx,ky] = a[0:N[2,0],3];\n",
    "\n",
    "eda_draw('title z(0,2)',zprev,'   ','title a(0,3)', m);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained network\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22832 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------------\n",
      "Trained network\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23728 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 6860 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting error 847.81\n",
      "ending error   0.05\n"
     ]
    }
   ],
   "source": [
    "# eda11_13\n",
    "# simple least squares fit of gaussian with a single tower\n",
    "\n",
    "# grid of points in (x,y)\n",
    "Nx = 21;\n",
    "Ny = 21;\n",
    "xmax = 5;\n",
    "ymax = 5;\n",
    "xmin = 0;\n",
    "ymin = 0;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "y = eda_cvec(np.linspace(ymin,ymax,Ny));\n",
    "\n",
    "# (x,y) index arrays for easy access\n",
    "kofij=np.zeros((Nx,Ny),dtype=np.int);  # linear index k of a particular (i,j)\n",
    "iofk=np.zeros((Nx*Ny,1),dtype=np.int); # row index of a particular k\n",
    "jofk=np.zeros((Nx*Ny,1),dtype=np.int); # col index of a particular k\n",
    "Nxy=0;\n",
    "for i in range(Nx):\n",
    "    for j in range(Ny):\n",
    "        kofij[i,j]=Nxy;\n",
    "        iofk[Nxy,0]=i;\n",
    "        jofk[Nxy,0]=j;\n",
    "        Nxy=Nxy+1;\n",
    "\n",
    "# Gaussian\n",
    "dobs = np.zeros((Nx,Ny));\n",
    "dobsvec = np.zeros((Nxy,1));\n",
    "xbar = 2.4;\n",
    "ybar = 2.6;\n",
    "amp = 5.0;\n",
    "sigmax = 1.2;\n",
    "sigmay = 0.8;\n",
    "cx = 1/(2*sigmax*sigmax);\n",
    "cy = 1/(2*sigmay*sigmay);\n",
    "for i in range(Nx):\n",
    "    for j in range(Ny):\n",
    "        myd = amp*exp(-cx*((x[i,0]-xbar)**2) - cy*((y[j,0]-ybar)**2) );\n",
    "        dobs[i,j] = myd;\n",
    "        dobsvec[kofij[i,j],0] = myd;\n",
    "\n",
    "# build the test network for a single 2D tower\n",
    "Nc = 1; # One tower\n",
    "Xc = np.zeros((Nc,1));\n",
    "Xc[:,0] = [2.5]; # Xc: the center of the tower\n",
    "Yc = np.zeros((Nc,1));\n",
    "Yc[:,0] = [2.5];\n",
    "Dx  = np.zeros((Nc,1));\n",
    "Dx[:,0] = [2];   # Dx: the width of the tower\n",
    "Dy = np.zeros((Nc,1));\n",
    "Dy[:,0] = [2];\n",
    "h = np.zeros((Nc,1));\n",
    "h[:,0] = [1];    # h: the height of the tower\n",
    "slope = 5/4;     # slope of the component step function\n",
    "\n",
    "# ---Neural Net definitions---\n",
    "Lmax = 4;      # number of layers\n",
    "Nmax = Nc * 4; # maximum neurons in any layer\n",
    "\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "\n",
    "# daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax\n",
    "# due to a change in the bias of the of the j-th neuron in the layer k\n",
    "\n",
    "# daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax\n",
    "# due to a change in the weight connecting the j-th neuron in the layer l\n",
    "# with the k-th neuron of layer l-1\n",
    "\n",
    "N,w,b = eda_net_2dtower( Nc, slope, Xc, Yc, Dx, Dy, h, N, w, b );\n",
    "\n",
    "print('Untrained network');\n",
    "fig2 = plt.figure(2,figsize=(14,10));\n",
    "ax2 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig2,ax2);\n",
    "plt.show();\n",
    "\n",
    "Nb, Nw, layer_of_bias, neuron_of_bias, index_of_bias, \\\n",
    "layer_of_weight, r_neuron_of_weight, l_neuron_of_weight, index_of_weight \\\n",
    "= eda_net_map(N,w,b);\n",
    "\n",
    "# linearized data kernel\n",
    "M = Nb+Nw;\n",
    "G = np.zeros((Nxy,M));\n",
    "\n",
    "# predicted data\n",
    "dpre = np.zeros((Nx,Ny));\n",
    "dprevec = np.zeros((Nxy,1));\n",
    "\n",
    "# initial predicted data\n",
    "d0 = np.zeros((Nx,Ny));\n",
    "d0vec = np.zeros((Nxy,1));\n",
    "\n",
    "Niter=20;\n",
    "for iter in range(Niter):\n",
    "    # loop over (x,y)\n",
    "    for kx in range(Nx):\n",
    "        for ky in range(Ny):\n",
    "            # activities\n",
    "            a[0:N[0,0],0] = [x[kx,0],x[ky,0]];\n",
    "            a, z = eda_net_eval( N,w,b,a);\n",
    "            dpre[kx,ky] = a[0,3];\n",
    "            dprevec[kofij[kx,ky],0]=dpre[kx,ky];\n",
    "    \n",
    "            # calculate the derivatives\n",
    "            [daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );\n",
    "            # daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax\n",
    "            # due to a change in the bias of the of the j-th neuron in the layer k\n",
    "            # daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax\n",
    "            # due to a change in the weight connecting the j-th neuron in the layer l\n",
    "            # with the k-th neuron of layer l-1\n",
    "    \n",
    "            # biases first\n",
    "            for ib in range(Nb):\n",
    "                nob=neuron_of_bias[ib,0];\n",
    "                lob=layer_of_bias[ib,0];\n",
    "                G[kofij[kx,ky], ib ] = daLmaxdb[0, nob, lob ];\n",
    "\n",
    "            # weights second\n",
    "            for iw in range(Nw):\n",
    "                rnow = r_neuron_of_weight[iw,0];\n",
    "                lnow = l_neuron_of_weight[iw,0];\n",
    "                low = layer_of_weight[iw,0];\n",
    "                G[kofij[kx,ky], Nb+iw ] = daLmaxdw[0, rnow, lnow, low ];\n",
    "\n",
    "    Dd = dobsvec - dprevec;\n",
    "    if( iter==0 ):\n",
    "        d0 = np.copy(dpre);\n",
    "        d0vec = np.copy(dprevec);\n",
    "        error = np.matmul( Dd.T, Dd );\n",
    "        error0 = error;\n",
    "\n",
    "    # implement some damping to improve stability\n",
    "    # but relax it after a few iterations\n",
    "    if( iter<=5 ):\n",
    "        epsi=1.0;\n",
    "    elif (iter<=10):\n",
    "        epsi=0.1;\n",
    "    else:\n",
    "        epsi=0.001;\n",
    "\n",
    "    GTG = np.matmul(G.T,G) + epsi*np.identity(M);\n",
    "    GTd = np.matmul(G.T,Dd);\n",
    "    Dm = la.solve(GTG,GTd);\n",
    "\n",
    "    # add in bias perturbations\n",
    "    for ib in range(Nb):\n",
    "        nob=neuron_of_bias[ib,0];\n",
    "        lob=layer_of_bias[ib,0];\n",
    "        b[nob, lob] = b[nob, lob] + Dm[ib,0];\n",
    "\n",
    "    # add in weights perturbations\n",
    "    for iw in range(Nw):\n",
    "        rnow = r_neuron_of_weight[iw,0];\n",
    "        lnow = l_neuron_of_weight[iw,0];\n",
    "        low = layer_of_weight[iw,0];\n",
    "        w[rnow, lnow, low] = w[rnow, lnow, low] + Dm[Nb+iw, 0];\n",
    "        \n",
    "    # loop over (x,y)\n",
    "    for kx in range(Nx):\n",
    "        for ky in range(Ny):\n",
    "            # activities\n",
    "            a[0:N[0,0],0] = [x[kx,0], x[ky,0]];\n",
    "            a, z = eda_net_eval( N,w,b,a);\n",
    "            dpre[kx,ky] = a[0,3];\n",
    "            dprevec[kofij[kx,ky],0]=dpre[kx,ky];\n",
    "            \n",
    "\n",
    "print('------------------------------------------------------------------' );\n",
    "print('Trained network');\n",
    "fig3 = plt.figure(3,figsize=(14,10));\n",
    "ax3 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig3,ax3);\n",
    "plt.show();\n",
    "\n",
    "eda_draw( 'title d-obs', dobs, 'title d-start', d0, 'title d-pre',dpre);\n",
    "\n",
    "Dd = dobsvec - dprevec;\n",
    "error = np.matmul(Dd.T,Dd);\n",
    "print('starting error %.2f' %(error0) );\n",
    "print('ending error   %.2f' %(error) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6216 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24124 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda11_14\n",
    "# a network that implements the linear function y=cx+d\n",
    "# for small w112 and w113 so that the approximation\n",
    "# is reasonably accurate in the range -1<x<+1\n",
    "# (when c, d, are both or order 1)\n",
    "\n",
    "c = 0.5; # true slope\n",
    "d = 2.0; # true intercept\n",
    "\n",
    "# standard x axis setup\n",
    "Nx = 101;\n",
    "xmin = -1;\n",
    "xmax = 1;\n",
    "x = eda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "\n",
    "# the function we are trying o approximate\n",
    "ytrue = c*x+d;\n",
    "\n",
    "# ---Neural Net definitions---\n",
    "\n",
    "# Nmax: number of layers;\n",
    "Lmax = 4;\n",
    "\n",
    "# Lmax: maximum number of neurons in any layer\n",
    "Nmax = 1;\n",
    "\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "\n",
    "N[:,0] = [1, 1, 1, 1];\n",
    "\n",
    "# these two slopes are arbitrary, as long as they are small\n",
    "slope1 = 0.01;\n",
    "slope2 = 0.02;\n",
    "w[0,0,1] = 4*slope1;\n",
    "w[0,0,2] = 4*slope2;\n",
    "w2w3 = w[0,0,1]*w[0,0,2];\n",
    "b[0,1] = 0;\n",
    "\n",
    "# the details:\n",
    "# near x=0, sigma(wx) = (w/4)*x + (1/2)\n",
    "# z2 = w2*x + b2 = w2*x + 0\n",
    "# a2 = sigma(z2) = w2/4+(1/2)\n",
    "# z3 = w3*a2 + b3 = w3*(w2/4+(1/2)+b3 = (w2*w3/4)*x + w3/2 + b3\n",
    "# so with choice b3 = -w3/2\n",
    "# z3 = (w2*w3/4)*x\n",
    "# a3 = sigma(z3) = (w2*w3/16)*x + (1/2)\n",
    "# z4 = w4*a3+b4\n",
    "# a4 = c*x + d\n",
    "\n",
    "w[0,0,3] = 16*c/w2w3;\n",
    "\n",
    "b[0,1] = 0;\n",
    "b[0,2] = -w[0,0,2]/2;\n",
    "b[0,3] = d-8*c/w2w3;\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(14,10));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig1,ax1);\n",
    "plt.show();\n",
    "\n",
    "# forward propagate the network for each value of x\n",
    "yest = np.zeros((Nx,1));\n",
    "for kx in range(Nx):\n",
    "    a[0,0] = x[kx,0]; \n",
    "    a, z = eda_net_eval( N,w,b,a);\n",
    "    yest[kx,0]=a[0,3];\n",
    "\n",
    "# plot the results\n",
    "fig2 = plt.figure(2,figsize=(6,6));\n",
    "ax2 = plt.subplot(1,1,1);\n",
    "plt.axis( [xmin, xmax, d+c*xmin, d+c*xmax] );\n",
    "plt.plot( x, ytrue, '-', linewidth=7, color=(0.8,0.8,0.8) );\n",
    "plt.plot( x, yest, 'k-', linewidth=2);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y(x)');\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stage 1: 100 1.000000e-01; Stage 1: 100 5.000000e-03\n",
      "---------------------- Initial guess ------------------------\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 51208 characters]\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Untrained network\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 49544 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------- Stage 1 ------------------------\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 51160 characters]\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stage 1 Trained network\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 52704 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting     error 0.75\n",
      "intermediate error 0.21\n",
      "---------------------- Second Pass ------------------------\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 50960 characters]\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stage 2 Trained network\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 89692 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting     error 0.75\n",
      "intermediate error 0.21\n",
      "final        error 0.04\n",
      "---------------------- Application to New Data ------------------------\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 50588 characters]\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "untrained error 0.83\n",
      "trained   error 0.25\n"
     ]
    }
   ],
   "source": [
    "# eda11_15\n",
    "# network to match non-linear response funtion\n",
    "\n",
    "# least squares settngs\n",
    "Niter1=100; # Stage 1: just tune existing network\n",
    "epsi1 = 0.1;\n",
    "Niter2=100; # Stage 2: add links and retune \n",
    "epsi2 = 0.005;\n",
    "print('Stage 1: %d %.6e; Stage 1: %d %.6e' % (Niter1, epsi1, Niter2, epsi2) );\n",
    "\n",
    "# load NonLinear Response data\n",
    "D = np.genfromtxt('../Data/nonlinearresponse.txt', delimiter='\\t');\n",
    "N, tmp = np.shape(D);\n",
    "\n",
    "#  break out into separate variables, metric conversion\n",
    "traw = eda_cvec(D[:,0]); # time in days\n",
    "dobsraw = eda_cvec(D[:,1])/35.3146;  # discharge in cubic meters per second\n",
    "xraw = eda_cvec(D[:,2])*25.4; # precipitation in millimeters per day\n",
    "\n",
    "# in this example, we use only a Nx day portion of the\n",
    "# dataset, mainly so that when we plot it, storm events\n",
    "# are clearly visible\n",
    "\n",
    "Nx=501;\n",
    "istart = 0;\n",
    "t = eda_cvec(traw[istart:istart+Nx,0]);\n",
    "x = eda_cvec(xraw[istart:istart+Nx,0]);\n",
    "dobs = eda_cvec(dobsraw[istart:istart+Nx,0]);\n",
    "tmin = np.min(t);\n",
    "tmax = np.max(t);\n",
    "xnorm = 1/np.max(x);        # x normalization constant\n",
    "dobsnorm = 1/np.max(dobs);  # data normalization constant\n",
    "x = x * xnorm;              # normalied x\n",
    "dobs = dobs * dobsnorm;     # normalized data\n",
    "\n",
    "print('---------------------- Initial guess ------------------------' );\n",
    "\n",
    "# plot the discharge and precipitation data\n",
    "fig1 = plt.figure(1,figsize=(14,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, x, 'k-', linewidth=3);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('precip x (mm/day)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs (m3/s)');\n",
    "\n",
    "# build a guess for the filter\n",
    "Nc = 10;\n",
    "c0 = 2;\n",
    "To = 3; \n",
    "f0 = eda_cvec(np.exp(-np.linspace(0,Nc-1,Nc)/To)); # decaayng exponetial function\n",
    "f = np.zeros((Nc,1));\n",
    "f[0,0] = f0[0,0]/1000;    # manually reduce size of first element of filter\n",
    "f[1:Nc,0] = f0[0:Nc-1,0]; # rest of filter is decaying exponetial\n",
    "h = c0*f/np.sum(f);       # normalized fiter\n",
    "\n",
    "# ---Neural Net definitions---\n",
    "\n",
    "Lmax = 4; # number of layers\n",
    "Nmax = Nc; # maximum number of neurons in any layer\n",
    "[N,w,b,a] = eda_net_init(Lmax,Nmax);\n",
    "\n",
    "slope1 = 0.01;\n",
    "slope2 = 0.02;\n",
    "[N,w,b] = eda_net_linear( Nc, slope1, slope2, 0.0, h, N, w, b );\n",
    "\n",
    "# saved untrained network\n",
    "N_save = np.copy(N);\n",
    "w_save = np.copy(w);\n",
    "b_save = np.copy(b);\n",
    "a_save = np.copy(a);\n",
    "\n",
    "# predict data of untrained network\n",
    "d0 = np.zeros((Nx,1));\n",
    "for kx in range(Nx): # loop over starting point of filter\n",
    "    # set bottom activities for this kx\n",
    "    for j in range(N[0,0]): # loop over points in filter\n",
    "        k = kx-j;\n",
    "        if( k>=0 ):\n",
    "            a[j,0] = x[k,0];\n",
    "        else:\n",
    "            a[j,0] = 0;  \n",
    "        # evaluate network\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        d0[kx,0] = a[0,3];\n",
    "\n",
    "# plot prediction of untrained network\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, '-', linewidth=7, color=(0.8,0.8,0.8));\n",
    "plt.plot( t, d0, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs and d0 (m3/s)');\n",
    "plt.show();\n",
    "\n",
    "print('Untrained network');\n",
    "fig2 = plt.figure(2,figsize=(14,10));\n",
    "ax2 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig2,ax2);\n",
    "plt.show();\n",
    "\n",
    "# ---------------------   First Pass, no links added ----------------------------------\n",
    "\n",
    "Nb, Nw, layer_of_bias, neuron_of_bias, index_of_bias, \\\n",
    "layer_of_weight, r_neuron_of_weight, l_neuron_of_weight, index_of_weight \\\n",
    "= eda_net_map(N,w,b);\n",
    "\n",
    "# linearized data kernel\n",
    "M = Nb+Nw;\n",
    "G = np.zeros((Nx,M));\n",
    "\n",
    "# predicted data\n",
    "dpre = np.zeros((Nx,1));\n",
    "\n",
    "for iter in range(Niter1):\n",
    "    for kx in range(Nx): # loop over starting point of filter; also row of G\n",
    "        \n",
    "        # set bottom activities for this kx\n",
    "        for j in range(N[0,0]): # loop over points in filter\n",
    "            k = kx-j;\n",
    "            if( k>=0 ):\n",
    "                a[j,0] = x[k,0];\n",
    "            else:\n",
    "                a[j,0] = 0;\n",
    "                \n",
    "        # evaluate network\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        dpre[kx,0] = a[0,3];\n",
    "    \n",
    "        # calculate the derivatives\n",
    "        [daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );\n",
    "        # daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax\n",
    "        # due to a change in the bias of the of the j-th neuron in the layer k\n",
    "        # daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax\n",
    "        # due to a change in the weight connecting the j-th neuron in the layer l\n",
    "        # with the k-th neuron of layer l-1\n",
    "    \n",
    "        # biases first\n",
    "        for ib in range(Nb):\n",
    "            nob=neuron_of_bias[ib,0];\n",
    "            lob=layer_of_bias[ib,0];\n",
    "            G[kx, ib ] = daLmaxdb[0, nob, lob ];\n",
    "\n",
    "        # weights second\n",
    "        for iw in range(Nw):\n",
    "            rnow = r_neuron_of_weight[iw,0];\n",
    "            lnow = l_neuron_of_weight[iw,0];\n",
    "            low = layer_of_weight[iw,0];\n",
    "            G[kx, Nb+iw ] = daLmaxdw[0, rnow, lnow, low ];\n",
    "            \n",
    "        # end loop over rows of G\n",
    "\n",
    "    Dd = dobs - dpre;\n",
    "    if( iter==0 ):\n",
    "        d0 = np.copy(dpre);\n",
    "        error0 = np.matmul( Dd.T, Dd );\n",
    "\n",
    "    # implement some damping to improve stability\n",
    "    # but relax it after a few iterations\n",
    "    if( iter==50 ):\n",
    "        epsi1=epsi1/10;\n",
    "\n",
    "    GTG = np.matmul(G.T,G) + epsi1*np.identity(M);\n",
    "    GTd = np.matmul(G.T,Dd);\n",
    "    Dm = la.solve(GTG,GTd);\n",
    "\n",
    "    # add in bias perturbations\n",
    "    for ib in range(Nb):\n",
    "        nob=neuron_of_bias[ib,0];\n",
    "        lob=layer_of_bias[ib,0];\n",
    "        b[nob, lob] = b[nob, lob] + Dm[ib,0];\n",
    "\n",
    "    # add in weights perturbations\n",
    "    for iw in range(Nw):\n",
    "        rnow = r_neuron_of_weight[iw,0];\n",
    "        lnow = l_neuron_of_weight[iw,0];\n",
    "        low = layer_of_weight[iw,0];\n",
    "        w[rnow, lnow, low] = w[rnow, lnow, low] + Dm[Nb+iw, 0];\n",
    "        \n",
    "for kx in range(Nx): # loop over starting point of filter\n",
    "    # set bottom activities for this kx\n",
    "    for j in range(N[0,0]): # loop over points in filter\n",
    "        k = kx-j;\n",
    "        if( k>=0 ):\n",
    "            a[j,0] = x[k,0];\n",
    "        else:\n",
    "            a[j,0] = 0;  \n",
    "        # evaluate network\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        dpre[kx,0] = a[0,3];\n",
    "\n",
    "print('---------------------- Stage 1 ------------------------' );\n",
    "\n",
    "# plot the discharge and precipitation data\n",
    "fig3 = plt.figure(3,figsize=(14,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, x, 'k-', linewidth=3);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('precip x (mm/day)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs (m3/s)');\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, '-', linewidth=7, color=(0.8,0.8,0.8));\n",
    "plt.plot( t, dpre, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs and dpre (m3/s)');\n",
    "plt.show();\n",
    "\n",
    "print('Stage 1 Trained network');\n",
    "fig4 = plt.figure(4,figsize=(14,10));\n",
    "ax4 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig4,ax4);\n",
    "plt.show();\n",
    "\n",
    "Dd = dobs - dpre;\n",
    "error = np.matmul(Dd.T,Dd);\n",
    "error1 = error;\n",
    "print('starting     error %.2f' %(error0) );\n",
    "print('intermediate error %.2f' %(error1) );\n",
    "\n",
    "\n",
    "# ---------------------   Second Pass, links added  ----------------------------------\n",
    "\n",
    "i=0;\n",
    "for j in range(N[0,0]):\n",
    "    if( w[i,j,1] == 0 ):\n",
    "        w[i,j,1] = np.random.normal(0,0.0001,1);\n",
    "i=0\n",
    "for j in range(N[0,0]):\n",
    "    if( w[i,j,2] == 0 ):\n",
    "        w[i,j,2] = np.random.normal(0,0.0001,1);\n",
    "\n",
    "Nb, Nw, layer_of_bias, neuron_of_bias, index_of_bias, \\\n",
    "layer_of_weight, r_neuron_of_weight, l_neuron_of_weight, index_of_weight \\\n",
    "= eda_net_map(N,w,b);\n",
    "\n",
    "# linearized data kernel\n",
    "M = Nb+Nw;\n",
    "G = np.zeros((Nx,M));\n",
    "\n",
    "# predicted data\n",
    "dpre = np.zeros((Nx,1));\n",
    "\n",
    "for iter in range(Niter2):\n",
    "    for kx in range(Nx): # loop over starting point of filter; also row of G\n",
    "        \n",
    "        # set bottom activities for this kx\n",
    "        for j in range(N[0,0]): # loop over points in filter\n",
    "            k = kx-j;\n",
    "            if( k>=0 ):\n",
    "                a[j,0] = x[k,0];\n",
    "            else:\n",
    "                a[j,0] = 0;\n",
    "                \n",
    "        # evaluate network\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        dpre[kx,0] = a[0,3];\n",
    "    \n",
    "        # calculate the derivatives\n",
    "        [daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );\n",
    "        # daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax\n",
    "        # due to a change in the bias of the of the j-th neuron in the layer k\n",
    "        # daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax\n",
    "        # due to a change in the weight connecting the j-th neuron in the layer l\n",
    "        # with the k-th neuron of layer l-1\n",
    "    \n",
    "        # biases first\n",
    "        for ib in range(Nb):\n",
    "            nob=neuron_of_bias[ib,0];\n",
    "            lob=layer_of_bias[ib,0];\n",
    "            G[kx, ib ] = daLmaxdb[0, nob, lob ];\n",
    "\n",
    "        # weights second\n",
    "        for iw in range(Nw):\n",
    "            rnow = r_neuron_of_weight[iw,0];\n",
    "            lnow = l_neuron_of_weight[iw,0];\n",
    "            low = layer_of_weight[iw,0];\n",
    "            G[kx, Nb+iw ] = daLmaxdw[0, rnow, lnow, low ];\n",
    "            \n",
    "        # end loop over rows of G\n",
    "\n",
    "    Dd = dobs - dpre;\n",
    "\n",
    "    # implement some damping to improve stability\n",
    "    # but relax it after a few iterations\n",
    "    if( iter==50 ):\n",
    "        epsi1=epsi2/10;\n",
    "\n",
    "    GTG = np.matmul(G.T,G) + epsi2*np.identity(M);\n",
    "    GTd = np.matmul(G.T,Dd);\n",
    "    Dm = la.solve(GTG,GTd);\n",
    "\n",
    "    # add in bias perturbations\n",
    "    for ib in range(Nb):\n",
    "        nob=neuron_of_bias[ib,0];\n",
    "        lob=layer_of_bias[ib,0];\n",
    "        b[nob, lob] = b[nob, lob] + Dm[ib,0];\n",
    "\n",
    "    # add in weights perturbations\n",
    "    for iw in range(Nw):\n",
    "        rnow = r_neuron_of_weight[iw,0];\n",
    "        lnow = l_neuron_of_weight[iw,0];\n",
    "        low = layer_of_weight[iw,0];\n",
    "        w[rnow, lnow, low] = w[rnow, lnow, low] + Dm[Nb+iw, 0];\n",
    "        \n",
    "for kx in range(Nx): # loop over starting point of filter\n",
    "    # set bottom activities for this kx\n",
    "    for j in range(N[0,0]): # loop over points in filter\n",
    "        k = kx-j;\n",
    "        if( k>=0 ):\n",
    "            a[j,0] = x[k,0];\n",
    "        else:\n",
    "            a[j,0] = 0;  \n",
    "        # evaluate network\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        dpre[kx,0] = a[0,3];\n",
    "\n",
    "print('---------------------- Second Pass ------------------------' );\n",
    "\n",
    "# plot the discharge and precipitation data\n",
    "fig5 = plt.figure(3,figsize=(14,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, x, 'k-', linewidth=3);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('precip x (mm/day)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs (m3/s)');\n",
    "ax3 = plt.subplot(3,1,3)\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, '-', linewidth=7, color=(0.8,0.8,0.8));\n",
    "plt.plot( t, dpre, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs and dpre (m3/s)');\n",
    "plt.show();\n",
    "\n",
    "print('Stage 2 Trained network');\n",
    "fig6 = plt.figure(6,figsize=(14,10));\n",
    "ax6 = plt.subplot(1,1,1);\n",
    "eda_net_view(N,w,b,fig6,ax6);\n",
    "plt.show();\n",
    "\n",
    "Dd = dobs - dpre;\n",
    "error = np.matmul(Dd.T,Dd);\n",
    "print('starting     error %.2f' %(error0) );\n",
    "print('intermediate error %.2f' %(error1) );\n",
    "print('final        error %.2f' %(error) );\n",
    "\n",
    "# Now test against a new dataset\n",
    "Nx=501;\n",
    "istart = 500;\n",
    "t = np.zeros((Nx,1));\n",
    "t[:,0] = traw[istart:istart+Nx,0];\n",
    "x = np.zeros((Nx,1));\n",
    "x[:,0] = xraw[istart:istart+Nx,0];\n",
    "dobs = np.zeros((Nx,1))\n",
    "dobs[:,0] = dobsraw[istart:istart+Nx,0];\n",
    "tmin = np.min(t);\n",
    "tmax = np.max(t);\n",
    "x = x * xnorm;\n",
    "dobs = dobs * dobsnorm;\n",
    "\n",
    "for kx in range(Nx): # loop over starting point of filter\n",
    "    # set bottom activities for this kx\n",
    "    for j in range(N[0,0]): # loop over points in filter\n",
    "        k = kx-j;\n",
    "        if( k>=0 ):\n",
    "            a[j,0] = x[k,0];\n",
    "            a_save[j,0] = x[k,0];\n",
    "        else:\n",
    "            a[j,0] = 0;  \n",
    "            a_save[j,0] = 0;\n",
    "        # evaluate network\n",
    "        a_save, z = eda_net_eval( N_save,w_save,b_save,a_save);\n",
    "        a, z = eda_net_eval( N,w,b,a);\n",
    "        d0[kx,0] = a_save[0,3];\n",
    "        dpre[kx,0] = a[0,3];\n",
    "\n",
    "\n",
    "print('---------------------- Application to New Data ------------------------' );\n",
    "\n",
    "# plot the discharge and precipitation data\n",
    "fig10 = plt.figure(10,figsize=(14,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, x, 'k-', linewidth=3);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('precip x (mm/day)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, '-', linewidth=4, color=(0.8,0.8,0.8));\n",
    "plt.plot( t, d0, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs and d0 (m3/s)');\n",
    "ax3 = plt.subplot(3,1,3)\n",
    "plt.axis( [tmin, tmax, -0.5, 1.5 ] );\n",
    "plt.plot( t, dobs, '-', linewidth=4, color=(0.8,0.8,0.8));\n",
    "plt.plot( t, dpre, 'k-', linewidth=1);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', linewidth=1);\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('dobs and dpre (m3/s)');\n",
    "plt.show();\n",
    "\n",
    "Dd = dobs - d0;\n",
    "error = np.matmul(Dd.T,Dd);\n",
    "print('untrained error %.2f' %(error) );\n",
    "Dd = dobs - dpre;\n",
    "error = np.matmul(Dd.T,Dd);\n",
    "print('trained   error %.2f' %(error) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 64520 characters]\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_11_16, make data file of t, d(t), x(t) that\n",
    "# is ued in edapy_11_15 by using Eulers's method to\n",
    "# numerically integerate the nonlinear differential\n",
    "# equation developed in the text\n",
    "\n",
    "# Note: in order to prevent the overwriting of the\n",
    "# '../Data/nonlinearresponse.txt file of the standard\n",
    "# distribution, this script writes the file\n",
    "# '../Data/mynonlinearresponse.txt, instead.\n",
    "\n",
    "# stardard setup of a time axis, t\n",
    "N=1001;\n",
    "Dt = 1.0;\n",
    "t = eda_cvec(Dt*np.linspace(0,N-1,N));\n",
    "\n",
    "# main variables\n",
    "x = np.zeros((N,1));  # precipitation\n",
    "y = np.zeros((N,1));  # watershed volume\n",
    "d = np.zeros((N,1));  # discharge\n",
    "\n",
    "# precipitation events are randomly placed spikes\n",
    "Np = floor(N/10);\n",
    "# uniformly distributed times\n",
    "p = np.random.randint(low=0,high=N,size=(Np,1));\n",
    "cexp = 0.5;\n",
    "# exponentially-distributed amplitudes\n",
    "x[p,0] = np.random.exponential(scale=cexp,size=(Np,1));\n",
    "# plot precipitation, x(t)\n",
    "fig1 = plt.figure(1,figsize=(14,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [0, N, 0, 2 ] );\n",
    "plt.plot( t, x, 'k-', linewidth=2);\n",
    "plt.plot( [tmin,tmax], [0,0], 'k-', lw=1);\n",
    "plt.xlabel('time, t');\n",
    "plt.ylabel('precipitation, x');\n",
    "\n",
    "# quick and dirty integration of the differential\n",
    "# equation by Euler's method.  Since y(t) is only\n",
    "# exemplary, there is no special reason for the\n",
    "# solution to be especially accurate\n",
    "cf = 0.2;\n",
    "for i in range(1,N):  # loop over times\n",
    "    d[i,0] = cf*(y[i-1,0]**2);\n",
    "    dydt = x[i,0] - d[i,0];\n",
    "    y[i,0] = y[i-1,0] + dydt * Dt;\n",
    "    if( y[i,0] < 0.0 ): # negative volumes not allowed\n",
    "        y[i,0]=0.0;\n",
    "\n",
    "# plot discharge\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis( [0, N, 0, 1.1*max(d) ] );\n",
    "plt.plot(t,d,'k-',lw=2);\n",
    "plt.xlabel('time, t');\n",
    "plt.ylabel('discharge, d(t)');\n",
    "\n",
    "# plot volume\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis( [0, N, 0, 1.1*max(y) ] );\n",
    "plt.plot(t,y,'k-','LineWidth',2);\n",
    "plt.xlabel('time,t');\n",
    "plt.ylabel('volume, y(t)');\n",
    "plt.show();\n",
    "\n",
    "# concatenate time, discharge, precip into a matrix C\n",
    "C = np.concatenate( (t,d,x),axis=1); \n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(\"../Data/mynonlinearresponse.txt\",C, delimiter=\"\\t\");\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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
