{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_05_00 clear all variables and import vatious modules\n",
    "# this cell must be exectuted first, before any of the others\n",
    "\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import scipy.sparse.linalg as las\n",
    "from scipy import sparse\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\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",
    "\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",
    "# function to perform the multiplication FT (F v)\n",
    "def FTFmul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the geneneralized least\n",
    "    # squares problem Fm=f.  Note that \"F\" must be called \"edaFsparse\".\n",
    "    global edaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    temp = edaFsparse*vv;\n",
    "    return edaFsparse.transpose()*temp;\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1: prior [10.] estimate [13.]\n",
      "m2: prior [20.] estimate [15.]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 6240 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_01: application of Bayes theorem to least squares\n",
    "\n",
    "# set up vectors m1 and m2\n",
    "L = 40;\n",
    "Dm = 1.0;\n",
    "m1 = eda_cvec( np.linspace(0,Dm*(L-1),L) );\n",
    "m2 = eda_cvec( np.linspace(0,Dm*(L-1),L) );\n",
    "\n",
    "# make normal distribution pp(m)\n",
    "m1bar=10;\n",
    "m2bar=20;\n",
    "Cpm=[[6,0],[0,10]];\n",
    "CpmI=la.inv(Cpm);\n",
    "norm=1/(2*pi*sqrt(la.det(Cpm)));\n",
    "Pp=np.zeros((L,L));\n",
    "dm = np.zeros((2,1));\n",
    "for i in range(L):\n",
    "    for j in range(L):\n",
    "        dm = eda_cvec( [m1[i,0]-m1bar, m2[j,0]-m2bar] );\n",
    "        Pp[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dm.T, CpmI ), dm ));\n",
    "\n",
    "# search grid for maximum \n",
    "k = np.argmax( Pp, axis=0 );\n",
    "Pptemp = np.max( Pp, axis=0 );\n",
    "j = np.argmax(Pptemp);\n",
    "i=k[j];\n",
    "m1bar=m1[i];\n",
    "m2bar=m2[j];\n",
    "\n",
    "# make normal distribution pdgm(m)\n",
    "# with a one-row data kernel that the\n",
    "# difference of two model parameters is d1\n",
    "\n",
    "G = np.zeros((1,2));\n",
    "G[0,:] = [1, -1];\n",
    "d1 = 0;\n",
    "Cd = 3;\n",
    "CdI = 1/Cd;\n",
    "norm = 1/sqrt(2*pi*Cd);\n",
    "Pdgm = np.zeros((L,L));\n",
    "m = np.zeros((2,1));\n",
    "for i in range(L):\n",
    "    for j in range(L):\n",
    "        m[:,0] = [m1[i,0], m2[j,0]];\n",
    "        dd = d1 - np.matmul( G, m );\n",
    "        Pdgm[i,j]=norm*np.exp( -0.5 * CdI * np.matmul( dd.T,dd ) );\n",
    "\n",
    "# multiply the two to get Pmgd\n",
    "Pmgd = np.multiply( Pp , Pdgm);\n",
    "                \n",
    "# search grid for maximum \n",
    "k = np.argmax( Pmgd, axis=0 );     # row indices of Pmgd with the largest value in each col\n",
    "Pmgdtemp = np.max( Pmgd, axis=0 ); # maximums in each column\n",
    "j = np.argmax(Pmgdtemp);           # indice to the column with the overall maximum\n",
    "i=k[j];                            # row of the overall maximum\n",
    "m1est=m1[i];                       # m1 position on the grid of the maximum\n",
    "m2est=m2[j];                       # m2 position on the grid of the maximum\n",
    "                \n",
    "print('m1: prior', m1bar, 'estimate', m1est);\n",
    "print('m2: prior', m2bar, 'estimate', m2est);\n",
    "\n",
    "# plot result\n",
    "eda_draw('title pp(d)',Pp,'title p(d|m)',Pdgm,'title p(m|d)',Pmgd);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5292 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_02: illustrate the product of twos Normal p.d.f.s\n",
    "\n",
    "# set up vectors m1 and m2\n",
    "L = 40;\n",
    "Dd = 1.0;\n",
    "d1 = eda_cvec( np.linspace(0,Dd*(L-1),L) );\n",
    "d2 = eda_cvec( np.linspace(0,Dd*(L-1),L) );\n",
    "\n",
    "# first Normal p.d.f.\n",
    "d1bar=15; # d1 mean\n",
    "d2bar=15; # d2 mean\n",
    "C=[[10,8],[8,10]]; # covariance\n",
    "CI=la.inv(C);      # its inverse\n",
    "norm=1/(2*pi*sqrt(la.det(C)));  # normalization factor for 2D Normal pdf\n",
    "pa=np.zeros((L,L)); # matrix to hold p.d.f.\n",
    "dd = np.zeros((2,1)); # vector to hold difference between value and mean\n",
    "for i in range(L):\n",
    "    for j in range(L):\n",
    "        dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];\n",
    "        pa[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd ));\n",
    "        \n",
    "# second Normal p.d.f.\n",
    "d1bar=15;\n",
    "d2bar=25;\n",
    "C=[[10,-8],[-8,10]];\n",
    "CI=la.inv(C);\n",
    "norm=1/(2*pi*sqrt(la.det(C)));\n",
    "pb=np.zeros((L,L));\n",
    "dd = np.zeros((2,1));\n",
    "for i in range(L):\n",
    "    for j in range(L):\n",
    "        dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];\n",
    "        pb[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd ));\n",
    "\n",
    "# product pc = pa * pb\n",
    "pc=np.zeros((L,L));\n",
    "pc = np.multiply(pa, pb);\n",
    "\n",
    "# normalize to unit area\n",
    "norm = (Dd**2)*np.sum(pc);  # sum is approximation for integral\n",
    "pc = pc/norm;\n",
    "\n",
    "# plot results\n",
    "eda_draw('title pA',pa,'title pB',pb,'title pC', pc);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28172 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_03: straight-line fit with unequal variance\n",
    "# illustrate least squares fitting of a line\n",
    "# when the data have unequal variance\n",
    "# this version uses G and d\n",
    "\n",
    "# set of vector of x's\n",
    "N = 51;\n",
    "N2 = floor((N-1)/2);\n",
    "Dx = 2;\n",
    "xmin = 0.0;\n",
    "xmax = xmin+Dx*(N-1);\n",
    "x = eda_cvec( np.linspace(xmin, xmax, N) );\n",
    "\n",
    "# prepare some test data hat consists of\n",
    "# two sections, that differ slightly in\n",
    "# interrcept and slope\n",
    "M=2;\n",
    "a = 1.0;\n",
    "b = 2.0;\n",
    "dobs = a + b*x;\n",
    "Da = 1;\n",
    "Db = 0.3;\n",
    "dobs[N2:N,0:1] = dobs[N2:N,0:1] + Da + Db*x[N2:N,0:1];\n",
    "\n",
    "# Part 1: variance, high for thr first section low for the second section\n",
    "sd = np.zeros((N,1));\n",
    "sd[0:N2,0:1] = 20.0;\n",
    "sd[N2:N,0:1] = 2.0;\n",
    "sd2i = np.power(sd,-2);\n",
    "Cdi = np.diag(sd2i.ravel());\n",
    "\n",
    "# set up the data kernel\n",
    "G=np.zeros((N,2));\n",
    "G[:,0]=np.ones((N,1)).ravel();\n",
    "G[:,1]=x.ravel();\n",
    "\n",
    "# solve for estimated model parameters\n",
    "FTF = np.matmul( G.T , np.matmul( Cdi, G ) );\n",
    "FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );\n",
    "mest = la.solve(FTF,FTf);\n",
    "# predict data\n",
    "dpre = mest[0] + mest[1] * x;\n",
    "             \n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([x[0,0], x[N-1,0], 0, 300])\n",
    "plt.plot(x,dobs,'ko');\n",
    "# 95% error bars\n",
    "for p in range(N):\n",
    "    plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "\n",
    "# Part 2: variance, high for thr first section low for the second section\n",
    "sd = np.zeros((N,1));\n",
    "sd[0:N2,0:1] = 2.0;\n",
    "sd[N2:N,0:1] = 20.0;\n",
    "sd2i = np.power(sd,-2);\n",
    "Cdi = np.diag(sd2i.ravel());\n",
    "\n",
    "# set up the data kernel\n",
    "G[0:N,0:1]=np.ones((N,1));\n",
    "G[0:N,1:2]=x;\n",
    "\n",
    "# solve for estimated model parameters\n",
    "FTF = np.matmul( G.T , np.matmul( Cdi, G ) );\n",
    "FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );\n",
    "mest = la.solve(FTF,FTf);\n",
    "# predict data\n",
    "dpre = mest[0,0] + mest[1,0] * x;\n",
    "             \n",
    "#  plot data\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([x[0,0], x[N-1,0], 0, 300])\n",
    "plt.plot(x,dobs,'ko');\n",
    "for p in range(N):\n",
    "    plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14524 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_04: Straight-line fit with diagonal variance\n",
    "# illustrate least squares fitting of a line\n",
    "# when the data have unequal variance\n",
    "# this version uses F and f\n",
    "\n",
    "# set of vector of x's\n",
    "N = 51;\n",
    "N2 = floor((N-1)/2);\n",
    "Dx = 2;\n",
    "xmin = 0.0;\n",
    "xmax = xmin+Dx*(N-1);\n",
    "x = eda_cvec( np.linspace(xmin, xmax, N) );\n",
    "\n",
    "M=2; # number of model parametres\n",
    "\n",
    "# prepare some test data hat consists of\n",
    "# two sections, that differ slightly in\n",
    "# interrcept and slope\n",
    "a = 1.0;\n",
    "b = 2.0;\n",
    "dobs = a + b*x;\n",
    "Da = 1;\n",
    "Db = 0.3;\n",
    "dobs[N2:N,0] = np.copy( dobs[N2:N,0] + Da + Db*x[N2:N,0] );\n",
    "\n",
    "# Part 1: variance, high for thr first section low for the second section\n",
    "sd = np.zeros((N,1));\n",
    "sd[0:N2,0] = 20;\n",
    "sd[N2:N,0] = 2;\n",
    "\n",
    "# set up the data kernel\n",
    "G=np.zeros((N,2));\n",
    "G[0:N,0:1]=np.ones((N,1));\n",
    "G[0:N,1:2]=x;\n",
    "\n",
    "# set up F and f\n",
    "F = np.zeros((N,M));\n",
    "for i in range(N):\n",
    "    F[i:i+1,0:M] = G[i:i+1,0:M] / sd[i,0];\n",
    "f = np.divide( dobs, sd );\n",
    "\n",
    "# solve for estimated model parameters\n",
    "FTF = np.matmul( F.T, F );\n",
    "FTf = np.matmul( F.T, f );\n",
    "mest = la.solve(FTF,FTf);\n",
    "\n",
    "# predict data\n",
    "dpre = mest[0] + mest[1] * x;\n",
    "             \n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([x[0,0], x[N-1,0], 0, 300])\n",
    "plt.plot(x,dobs,'ko');\n",
    "# 95% error barx\n",
    "for p in range(N):\n",
    "    plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', 'LineWidth', 2 );\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "built 141 rows out of 141\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15752 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_05: using smoothness to fill in data gaps\n",
    "# simple example of smoothness applied\n",
    "# to the problem of filling in data gaps\n",
    "# this version uses non-sparse arrays\n",
    "\n",
    "# the model parameters are the values of a sinusoidal\n",
    "# function evenly spaced along the x axis\n",
    "M=101;\n",
    "Dx=1.0;\n",
    "Dx2 = Dx**2;\n",
    "xmin=0.0;\n",
    "xmax = xmin+Dx*(M-1);\n",
    "x = eda_cvec( np.linspace(xmin, xmax, M) );\n",
    "\n",
    "# the data are the values the function measured on\n",
    "# only a subset of points along the x axis\n",
    "N=40;\n",
    "\n",
    "# randomly choose N of the M points to have data\n",
    "rowindex = np.sort( np.random.randint(M,size=N) );  # N random integers between 1 and M\n",
    "xd = eda_cvec( x[rowindex,:] );  # xd are the x's of the N data\n",
    "\n",
    "# data is sinusoidal\n",
    "freq=0.02;\n",
    "d = np.sin(2*pi*freq*xd);\n",
    "\n",
    "# prior variances\n",
    "# 1% error of signal level of 1 i 0.1\n",
    "sigmadi =  1.0/0.01;\n",
    "# maximum slope of about (1-0)/10 = 0.1\n",
    "sigmaD2i = 1.0/0.1;\n",
    "# maximum 2nd der about (0-2*1-0)/(20*20) = 0.005\n",
    "sigmaD1i = 1.0/0.005;\n",
    "\n",
    "# build F and f\n",
    "K=M;\n",
    "F=np.zeros((N+K,M));\n",
    "fobs=np.zeros((N+K,1));\n",
    "nextrow = 0;\n",
    "\n",
    "# the model parameter equals the observed data\n",
    "for p in range(N):\n",
    "    F[nextrow,rowindex[p]] = sigmadi;\n",
    "    fobs[nextrow,0]=sigmadi*d[p,0];\n",
    "    nextrow=nextrow+1;\n",
    "\n",
    "# smoothness (small second derivative) in the interior\n",
    "for p in range(M-2):\n",
    "    F[nextrow,p] = sigmaD2i/Dx2;\n",
    "    F[nextrow,p+1] = -2*sigmaD2i/Dx2;\n",
    "    F[nextrow,p+2] = sigmaD2i/Dx2;\n",
    "    fobs[nextrow,0]=0.0;\n",
    "    nextrow=nextrow+1;\n",
    "\n",
    "# flatness (small first derivative) at the left end\n",
    "F[nextrow,0]=-sigmaD1i*Dx;\n",
    "F[nextrow,1]=sigmaD1i*Dx;\n",
    "fobs[nextrow,0]=0;\n",
    "nextrow=nextrow+1;\n",
    "\n",
    "# flatness (small first derivative) at the right end\n",
    "F[nextrow,M-2]=-sigmaD1i*Dx;\n",
    "F[nextrow,M-1]=sigmaD1i*Dx;\n",
    "fobs[nextrow,0]=0;\n",
    "nextrow=nextrow+1;\n",
    "\n",
    "print(\"built %d rows out of %d\" % (nextrow,N+K));\n",
    "\n",
    "# solve generalized least squares equation\n",
    "FTF = np.matmul( F.T, F );\n",
    "FTf = np.matmul( F.T, fobs );\n",
    "mest = np.zeros((M,1));\n",
    "mest[0:M,0:1] = la.solve(FTF,FTf);\n",
    "\n",
    "# error\n",
    "e = d-mest[rowindex];      # individual error\n",
    "E = np.matmul( e.T, e );   # overall error\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])\n",
    "plt.plot(xd,d,'ko');\n",
    "plt.plot(x,mest,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 2. 0. 0. 0.]\n",
      " [3. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 4. 0.]\n",
      " [0. 0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "# eda05_06: sparse matrix from list of non-zero elements\n",
    "\n",
    "# build lists\n",
    "L=7;\n",
    "M=6;\n",
    "Frow = eda_cvec( 0, 2, 3, 5);\n",
    "Fcol = eda_cvec( 0, 2, 0, 4);\n",
    "Fval = eda_cvec( 1.0, 2.0, 3.0, 4.0);\n",
    "KK, tmp = np.shape(Fval);\n",
    "\n",
    "# build sparse matrix\n",
    "edaFsparse = sparse.coo_matrix((Fval[0:KK,0],\n",
    "                (Frow[0:KK,0],Fcol[0:KK,0])), shape=(L,M) );\n",
    "\n",
    "M=sparse.csr_matrix.todense(edaFsparse);\n",
    "print(M);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "341 non-zero-elements, expected 341\n",
      "141 rows, expected 141\n",
      "E: [[3.1975708e-07]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15676 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_07: sparse matrix implementation of Fm=f\n",
    "# simple example of smoothness applied\n",
    "# to the problem of filling in data gaps\n",
    "# this version uses sparse arrays\n",
    "\n",
    "# the model parameters are the values of a sinusoidal\n",
    "# function evenly spaced along the x axis\n",
    "M=101;\n",
    "Dx=1.0;\n",
    "Dx2 = Dx**2;\n",
    "xmin=0.0;\n",
    "xmax = xmin+Dx*(M-1);\n",
    "x = eda_cvec( np.linspace(xmin, xmax, M) );\n",
    "\n",
    "# the data are the values the function measured on\n",
    "# only a subset of points along the x axis\n",
    "N=40;\n",
    "\n",
    "# randomly choose N of the M points to have data\n",
    "rowindex = np.sort( np.random.randint(M,size=N) ); # N random integers between 1 and M\n",
    "xd = eda_cvec( x[rowindex,:] ); # x's of the data\n",
    "\n",
    "# data is sinusoidal\n",
    "freq=0.02;\n",
    "dobs = np.sin(2*pi*freq*xd);\n",
    "\n",
    "# prior variances\n",
    "# 1% error of signal level of 1 i 0.1\n",
    "sigmadi =  1.0/0.01;\n",
    "# maximum slope of about (1-0)/10 = 0.1\n",
    "sigmaD2i = 1.0/0.1;\n",
    "# maximum 2nd der about (0-2*1-0)/(20*20) = 0.005\n",
    "sigmaD1i = 1.0/0.005;\n",
    "\n",
    "# set up f\n",
    "L=N+M;  # rows in F, N data, and N priors\n",
    "fobs=np.zeros((L,1));\n",
    "\n",
    "# set up F.  In this impletation, the non-zero\n",
    "# elements of F are put first in a vector Fval, with\n",
    "# two auxilary arrays Frow and Fcol\n",
    "K = int;\n",
    "K = N + 3*(M-2) + 2*2 + 100;  # needs to be at least as large as the number of non-zero elments in F\n",
    "Frow = np.zeros((K,1),dtype=np.int32);  # the row index vector is integer\n",
    "Fcol = np.zeros((K,1),dtype=np.int32);  # the col index vector is integer\n",
    "Fval = np.zeros((K,1),dtype=np.double); # the value vector is double precision\n",
    "\n",
    "# a trick for keeping tracks of rows is to use a row counter\n",
    "# an increment it after each row is added.  The counters\n",
    "# are defined as integers.\n",
    "nextrow = int;  # next row in F; note is an integer\n",
    "nextrow = 0;    # starting value of zero\n",
    "KK = int;  # next element in Fval; note is an integer\n",
    "KK = 0;    # starting value of zero\n",
    "\n",
    "# the first N rows of F are that the model parameter\n",
    "# equals the observed data\n",
    "for p in range(N):\n",
    "    # F[nextrow,rowindex[p]] = 1;\n",
    "    Frow[KK,0]=nextrow;\n",
    "    Fcol[KK,0]=rowindex[p];\n",
    "    Fval[KK,0]=sigmadi;\n",
    "    KK=KK+1;\n",
    "    # f[nextrow]=d[p];\n",
    "    fobs[nextrow,0]=sigmadi*dobs[p,0];\n",
    "    nextrow=nextrow+1;\n",
    "\n",
    "# implement the smoothness constraints in the interior\n",
    "# of the x interval\n",
    "for p in range(M-2):\n",
    "    # F[nextrow,p] = shi/Dx2;\n",
    "    Frow[KK,0]=nextrow;\n",
    "    Fcol[KK,0]=p;\n",
    "    Fval[KK,0]=sigmaD2i/Dx2;\n",
    "    KK=KK+1;\n",
    "    # F[nextrow,p+1] = -2*shi/Dx2;\n",
    "    Frow[KK,0]=nextrow;\n",
    "    Fcol[KK,0]=p+1;\n",
    "    Fval[KK,0]=-2*sigmaD2i/Dx2;\n",
    "    KK=KK+1;\n",
    "    # F[nextrow,p+2] = shi/Dx2;\n",
    "    Frow[KK,0]=nextrow;\n",
    "    Fcol[KK,0]=p+2;\n",
    "    Fval[KK,0]=sigmaD2i/Dx2;\n",
    "    KK=KK+1;\n",
    "    # f[nextrow,0]=0.0;\n",
    "    fobs[nextrow,0]=0.0;\n",
    "    nextrow=nextrow+1;\n",
    "\n",
    "# implement flatness constraint at left end of the x interval\n",
    "# F[nextrow,0]=-shi*Dx;\n",
    "Frow[KK,0]=nextrow;\n",
    "Fcol[KK,0]=0;\n",
    "Fval[KK,0]=-sigmaD1i*Dx;\n",
    "KK=KK+1;\n",
    "# F[nextrow,1]=shi*Dx;\n",
    "Frow[KK,0]=nextrow;\n",
    "Fcol[KK,0]=1;\n",
    "Fval[KK,0]=sigmaD1i*Dx;\n",
    "KK=KK+1;\n",
    "# f[nextrow,0]=0;\n",
    "fobs[nextrow,0]=0;\n",
    "nextrow=nextrow+1;\n",
    "\n",
    "# implement flatness constraint at right end of the x interval\n",
    "# F[nextrow,M-2]=-shi*Dx;\n",
    "Frow[KK,0]=nextrow;\n",
    "Fcol[KK,0]=M-2;\n",
    "Fval[KK,0]=-sigmaD1i*Dx;\n",
    "KK=KK+1;\n",
    "# F[nextrow,M-1]=shi*Dx;\n",
    "Frow[KK,0]=nextrow;\n",
    "Fcol[KK,0]=M-1;\n",
    "Fval[KK,0]=sigmaD1i*Dx;\n",
    "KK=KK+1;\n",
    "# f[nextrow]=0;\n",
    "fobs[nextrow]=0;\n",
    "nextrow=nextrow+1;\n",
    "\n",
    "print(\"%d non-zero-elements, expected %d\" % (KK, N+3*(M-2)+2*2) );\n",
    "print(\"%d rows, expected %d\" % (nextrow, N+M) );\n",
    "# actual rows in F\n",
    "\n",
    "# the \"F\" matrix must be called \"edaFsparse\" in order to correctly interface\n",
    "# witk the FYFmul() function called by the biconjugate gradient solver\n",
    "# Note the reduced length of the arrays Fval, Frow, Fcol\n",
    "edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0],\n",
    "                            Fcol[0:KK,0])), shape=((N+M),M) );\n",
    "del Fval;  # delete no loger needed arrays\n",
    "del Frow;\n",
    "del Fcol; \n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),\n",
    "                      matvec=FTFmul,rmatvec=FTFmul);\n",
    "\n",
    "# solve for eusing congugate gradient algorithm\n",
    "mest = np.zeros((M,1));\n",
    "q=las.bicg(LO,edaFsparse.T*fobs,tol=1e-12,\n",
    "           maxiter=(3*(N+M)+100));\n",
    "mest[:,0] = q[0];\n",
    "\n",
    "# error\n",
    "e = dobs-mest[rowindex];  # individual errors\n",
    "E = np.matmul( e.T, e );    # overall error\n",
    "print('E:',E);\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])\n",
    "plt.plot(xd,dobs,'ko');\n",
    "plt.plot(x,mest,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\n",
      " [[ 1.  5.  9.]\n",
      " [ 2.  6. 10.]\n",
      " [ 3.  7. 11.]\n",
      " [ 4.  8. 12.]]\n",
      "\n",
      "Part 1: with a formula\n",
      "m\n",
      " [[ 1.]\n",
      " [ 5.]\n",
      " [ 9.]\n",
      " [ 2.]\n",
      " [ 6.]\n",
      " [10.]\n",
      " [ 3.]\n",
      " [ 7.]\n",
      " [11.]\n",
      " [ 4.]\n",
      " [ 8.]\n",
      " [12.]]\n",
      "B\n",
      " [[ 1.  5.  9.]\n",
      " [ 2.  6. 10.]\n",
      " [ 3.  7. 11.]\n",
      " [ 4.  8. 12.]]\n",
      "\n",
      "Part 2: with index arrays\n",
      "m\n",
      " [[ 1.]\n",
      " [ 5.]\n",
      " [ 9.]\n",
      " [ 2.]\n",
      " [ 6.]\n",
      " [10.]\n",
      " [ 3.]\n",
      " [ 7.]\n",
      " [11.]\n",
      " [ 4.]\n",
      " [ 8.]\n",
      " [12.]]\n",
      "B\n",
      " [[ 1.  5.  9.]\n",
      " [ 2.  6. 10.]\n",
      " [ 3.  7. 11.]\n",
      " [ 4.  8. 12.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_05_08: folding out a matrix into a vector\n",
    "\n",
    "I = 4;   # number of rows of a matrix A\n",
    "J = 3;   # number of coluns of a matrix A\n",
    "M = I*J; # total number of elements\n",
    "\n",
    "# make a matrix A\n",
    "A=np.zeros((I,J));\n",
    "A[:,:] = [ [1, 5, 9], [2, 6, 10], [3, 7, 11], [4, 8, 12]];\n",
    "print('A\\n', A);\n",
    "\n",
    "# make a Mx1 vector m\n",
    "m=np.zeros((M,1));\n",
    "\n",
    "# PART 1: using a formula\n",
    "print('\\nPart 1: with a formula');\n",
    "\n",
    "# copy elements of A[i,j] into m[k]\n",
    "for i in range(I):\n",
    "    for j in range(J):\n",
    "        k = i*J+j; # row-wise unfold\n",
    "        m[k,0] = A[i,j];\n",
    "print('m\\n', m);\n",
    "\n",
    "# copy the elements of m[k] into a matrix B[i,j]\n",
    "B=np.zeros((I,J));\n",
    "for k in range(M):\n",
    "    i = floor(k/J); # row-wise fold\n",
    "    j = k-i*J;\n",
    "    B[i,j]=m[k];\n",
    "print('B\\n',B);\n",
    "\n",
    "# PART 2: with index arrays\n",
    "print('\\nPart 2: with index arrays')\n",
    "\n",
    "# the idea is to define a set of \"pointer arrays\"\n",
    "# that translate between (i,j) and k.  The arrays\n",
    "# are built during one traverse of (i,j).  One\n",
    "# advanage is that no \"formulas\" are needed.\n",
    "i_of_k  = np.zeros((M,1), dtype=np.int32);  # i given k\n",
    "j_of_k =  np.zeros((M,1), dtype=np.int32);  # j given k\n",
    "k_of_ij = np.zeros((I,J), dtype=np.int32);  # k given (i,j)\n",
    "k=0; # initialize counter\n",
    "for i in range(I):\n",
    "    for j in range(J):\n",
    "        i_of_k[k,0]=i;\n",
    "        j_of_k[k,0]=j;\n",
    "        k_of_ij[i,j]=k;\n",
    "        k=k+1;\n",
    "M=k;\n",
    "\n",
    "# copy A[i,j] to m[k]\n",
    "m=np.zeros((M,1));\n",
    "for i in range(I):\n",
    "    for j in range(J):\n",
    "        k = k_of_ij[i,j];\n",
    "        m[k,0] = A[i,j];\n",
    "print('m\\n', m);\n",
    "  \n",
    "# copy m[k] to B[i,j]\n",
    "B=np.zeros((I,J));\n",
    "for k in range(M):\n",
    "    i = i_of_k[k,0];\n",
    "    j = j_of_k[k,0];\n",
    "    B[i,j]=m[k,0];\n",
    "print('B\\n',B);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8093 non-zero elements\n",
      "1849 rows, expected 1849\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 8264 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior errror in data 2.3415950502294525\n"
     ]
    }
   ],
   "source": [
    "# edapy_05_09: fill in data gaps in 2D using diff eqn\n",
    "# use generealized least square to fill in\n",
    "# data gaps for data that satisfy Laplace's\n",
    "# equation in the in the interior of a grid,\n",
    "# and the normal derivative is zero on the\n",
    "# edges of the grid\n",
    "\n",
    "# load the data.  These data are only 'simulated'\n",
    "# (also called \"synthetic\") and were created by the\n",
    "# script eda05_09\n",
    "\n",
    "# load (x, y, pressure) data from file\n",
    "D = np.genfromtxt('../data/pressure.txt', delimiter='\\t')\n",
    "[N, K]=D.shape;\n",
    "xobs = eda_cvec( D[:,0] );\n",
    "yobs = eda_cvec( D[:,1] );\n",
    "dobs = eda_cvec( D[:,2] );\n",
    "\n",
    "# assign variance of data\n",
    "sigmad = 0.1;\n",
    "rsigmad = 1/sigmad;\n",
    "rsigmad2 = rsigmad**2;\n",
    "\n",
    "# rows of grid are x\n",
    "I = 41;\n",
    "Dx = 1;\n",
    "xmin = 0;\n",
    "xmax = Dx*(I-1);\n",
    "x = eda_cvec( np.linspace(xmin, xmax, I) );\n",
    "rDx = 1/Dx;\n",
    "rDx2 = rDx**2;\n",
    "\n",
    "# columns of grid are y\n",
    "J = 41;\n",
    "Dy = 1;\n",
    "ymin = 0;\n",
    "ymax = Dy*(J-1);\n",
    "y = eda_cvec( np.linspace(ymin, ymax, J) );\n",
    "rDy = 1/Dy;\n",
    "rDy2 = rDy**2;\n",
    "\n",
    "# associate observations with points on the grid\n",
    "# by scaling the range (xmin, xmax) to (1, N)\n",
    "rowindex = np.zeros( (N,1), dtype=np.int32 );\n",
    "rowindex[:,0] = np.floor(I*(xobs-xmin)/(xmax-xmin)+0.5).ravel();\n",
    "colindex = np.zeros( (N,1), dtype=np.int32 );\n",
    "colindex[:,0] = np.floor(J*(yobs-ymin)/(ymax-ymin)+0.5).ravel();\n",
    "\n",
    "# don't let any data fall outside the grid.  In this\n",
    "# case, we just put it on the edge of the grid.  The\n",
    "# alternative is to throw it away.\n",
    "for i in range(N):\n",
    "    if( rowindex[i,0]<0 ):\n",
    "        rowindex[i,0]=1;\n",
    "    elif( rowindex[i,0]>=I ):\n",
    "        rowindex[i,0]=I-1;\n",
    "    if( colindex[i,0]<0 ):\n",
    "        colindex[i,0]=0;\n",
    "    elif( colindex[i,0]>=J ):\n",
    "        colindex[i,0]=J-1;\n",
    "\n",
    "# copy the data to the grid A for later use\n",
    "Aobs = np.zeros((I,J));\n",
    "for i in range(N):\n",
    "    Aobs[rowindex[i,0],colindex[i,0]] = dobs[i,0];\n",
    "\n",
    "# one model parameter per grid point\n",
    "M = I*J;\n",
    "\n",
    "# Fm=f is LxM with N rows of data plus M rows of prior information\n",
    "L = N+M;\n",
    "\n",
    "# F will have no more than K elements\n",
    "K = N + 5*M + 100;\n",
    "\n",
    "# r.h.s. of Fm=F\n",
    "fobs = np.zeros((L,1));\n",
    "\n",
    "# reminders on how to move from grid to vector & vice-versa\n",
    "# m(k) = A(i,j);\n",
    "# k = i*J+j;\n",
    "# i = floor(k/J);\n",
    "# j = k-i*J;\n",
    "# B(i,j)=m(k);\n",
    "\n",
    "# in this implementation, I build vectors on zon-zero elements and\n",
    "# then change them into a sparse matrix at the end\n",
    "Frow = np.zeros((K,1),dtype=np.int32);\n",
    "Fcol = np.zeros((K,1),dtype=np.int32);\n",
    "Fval = np.zeros((K,1),dtype=np.double);\n",
    "\n",
    "# build F and f\n",
    "\n",
    "# nextrow is used to keep track of the current row of F and f\n",
    "nextrow = int;\n",
    "nextrow = 0;\n",
    "\n",
    "# KK is used to keep track of the current row of Fvec, Fcol, Frow\n",
    "KK = int;\n",
    "KK = 0;\n",
    "\n",
    "# k=dindex(i): data i corresponds to row k of m\n",
    "dindex = np.zeros((N,1),dtype=np.int32);\n",
    "\n",
    "# first part of Fm=f: the data\n",
    "for i in range(N):\n",
    "    k = rowindex[i]*J+colindex[i];\n",
    "    # F[q,k] = rsigmad;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = rsigmad;\n",
    "    KK = KK+1;\n",
    "    fobs[nextrow,0] = rsigmad*dobs[i,0];\n",
    "    dindex[i:0] = k;\n",
    "    nextrow = nextrow+1;\n",
    "\n",
    "# error of prior information\n",
    "sigmah = 1; # was 1000\n",
    "rsigmah = 1/sigmah;\n",
    "\n",
    "# second part, laplaces equation in interior\n",
    "# second derivative is a (1,-2,1)/Dx^2 in x plus a (1,-2,1)/Dy^2 in y\n",
    "for i in range(1,I-1): # Matlab [2:I-1]\n",
    "    for j in range(1,J-1): # Matlab [2:J-1]\n",
    "        # k = (i-1)*J+j; and F(q,k) = rsigmah*(-2*rDx2-2*rDy2);\n",
    "        k = i*J+j;\n",
    "        Frow[KK,0] = nextrow;\n",
    "        Fcol[KK,0] = k;\n",
    "        Fval[KK,0] = rsigmah*(-2*rDx2-2*rDy2);\n",
    "        KK=KK+1;\n",
    "        # k = ((i-1)-1)*J+j; and F(q,k) = rsigmah*rDx2;\n",
    "        k = (i-1)*J+j;\n",
    "        Frow[KK,0] = nextrow;\n",
    "        Fcol[KK,0] = k;\n",
    "        Fval[KK,0] = rsigmah*rDx2\n",
    "        KK=KK+1;\n",
    "        # k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx2;\n",
    "        k = (i+1)*J+j;\n",
    "        Frow[KK,0] = nextrow;\n",
    "        Fcol[KK,0] = k;\n",
    "        Fval[KK,0] = rsigmah*rDx2;\n",
    "        KK=KK+1;\n",
    "        #  k = (i-1)*J+(j-1); and F(q,k) = rsigmah*rDy2;\n",
    "        k = i*J+(j-1);\n",
    "        Frow[KK,0] = nextrow;\n",
    "        Fcol[KK,0] = k;\n",
    "        Fval[KK,0] = rsigmah*rDy2;\n",
    "        KK=KK+1;\n",
    "        # k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy2;\n",
    "        k = i*J+(j+1);\n",
    "        Frow[KK,0] = nextrow;\n",
    "        Fcol[KK,0] = k;\n",
    "        Fval[KK,0] = rsigmah*rDy2;\n",
    "        KK=KK+1;\n",
    "        fobs[nextrow,0] = 0.0;\n",
    "        nextrow = nextrow+1;\n",
    "\n",
    "# third part: top edge\n",
    "for j in range(1,J-1): # MATLAB [2:J-1]\n",
    "    i = 0;\n",
    "    # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDx;\n",
    "    k = i*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = -rsigmah*rDx;\n",
    "    KK=KK+1;\n",
    "    # k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx;\n",
    "    k = (i+1)*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = rsigmah*rDx;\n",
    "    KK=KK+1;\n",
    "    fobs[nextrow,0] = 0.0;\n",
    "    nextrow = nextrow+1;\n",
    "    \n",
    "# fourth part: bottom edge\n",
    "for j in range(1,J-1): # MATLAB [2:J-1]\n",
    "    i = I-1;\n",
    "    # k = (i-1)*J+j; and F(q,k) = rsigmah*rDx;\n",
    "    k = i*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = rsigmah*rDx;\n",
    "    KK=KK+1;\n",
    "    # k = ((i-1)-1)*J+j; and F(q,k) = -rsigmah*rDx;\n",
    "    k = (i-1)*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = -rsigmah*rDx;\n",
    "    KK=KK+1;\n",
    "    fobs[nextrow,0] = 0.0;\n",
    "    nextrow = nextrow+1;\n",
    "\n",
    "# fifth part: left edge\n",
    "for i in range(1,I-1): # MATLAB [2:I-1]:\n",
    "    j = 0;\n",
    "    # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDy; \n",
    "    k = i*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = -rsigmah*rDy;\n",
    "    KK=KK+1;\n",
    "    # k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy;\n",
    "    k = i*J+(j+1);\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = rsigmah*rDy;\n",
    "    KK=KK+1;\n",
    "    fobs[nextrow,0] = 0.0;\n",
    "    nextrow = nextrow+1;\n",
    "    \n",
    "# sixth part: right edge\n",
    "for i in range(1,I-1): # MATLAB [2:I-1]\n",
    "    j = J-1;\n",
    "    # k = (i-1)*J+j; and F(q,k) = rsigmah*rDy;\n",
    "    k = i*J+j;\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = rsigmah*rDy;\n",
    "    KK=KK+1;\n",
    "    # k = (i-1)*J+(j-1); and F(q,k) = -rsigmah*rDy;\n",
    "    k = i*J+(j-1);\n",
    "    Frow[KK,0] = nextrow;\n",
    "    Fcol[KK,0] = k;\n",
    "    Fval[KK,0] = -rsigmah*rDy;\n",
    "    KK=KK+1;\n",
    "    fobs[nextrow,0] = 0.0;\n",
    "    nextrow = nextrow+1;\n",
    "\n",
    "# seventh part: four corners\n",
    "# top left\n",
    "rDz = 1/sqrt((Dx^2)+(Dy^2));\n",
    "i=0;\n",
    "j=0;\n",
    "# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;\n",
    "k = i*J+j;\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = -rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "# k = ((i+1)-1)*J+(j+1); F(q,k) = rsigmah*rDz;\n",
    "k=(i+1)*J+(j+1);\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "fobs[nextrow,0]=0.0;\n",
    "nextrow=nextrow+1;\n",
    "# top right\n",
    "i=0;\n",
    "j=J-1;\n",
    "# k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;\n",
    "k = i*J+j;\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "# k = ((i+1)-1)*J+(j-1); F(q,k) = -rsigmah*rDz;\n",
    "k=(i+1)*J+(j-1);\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = -rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "fobs[nextrow,0]=0.0;\n",
    "nextrow=nextrow+1;\n",
    "# bottom left\n",
    "i=I-1;\n",
    "j=0;\n",
    "# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;\n",
    "k = i*J+j;\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = -rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "# k = ((i-1)-1)*J+(j+1); and F(q,k) = rsigmah*rDz;\n",
    "k = (i-1)*J+(j+1);\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "fobs[nextrow,0]=0.0;\n",
    "nextrow=nextrow+1;\n",
    "# bottom right\n",
    "i=I-1;\n",
    "j=0;\n",
    "# k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;\n",
    "k = i*J+j;\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "# k = ((i-1)-1)*J+(j-1); and F(q,k) = -rsigmah*rDz;\n",
    "k = (i-1)*J+(j-1);\n",
    "Frow[KK,0] = nextrow;\n",
    "Fcol[KK,0] = k;\n",
    "Fval[KK,0] = -rsigmah*rDz;\n",
    "KK=KK+1;\n",
    "fobs[nextrow,0]=0.0;\n",
    "nextrow=nextrow+1;\n",
    "\n",
    "# actual rows in F\n",
    "LL = int;\n",
    "LL = nextrow;\n",
    "\n",
    "print(\"%d non-zero elements\" % (KK));\n",
    "print(\"%d rows, expected %d\" % (LL, (N+M)));\n",
    "\n",
    "# the \"F\" matrix must be called \"edaFsparse\" in order to correctly interface\n",
    "# witk the FYFmul() function called by the biconjugate gradient solver\n",
    "# Note the reduced length of the arrays Fval, Frow, Fcol\n",
    "edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0], Fcol[0:KK,0])), shape=(LL,M) );\n",
    "del Fval;  # delete the arrays used to build the sparse matrix\n",
    "del Frow;\n",
    "del Fcol; \n",
    "\n",
    "# define linear operator needed for conjugate gradient solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=FTFmul,rmatvec=FTFmul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "mest = np.zeros((M,1));\n",
    "q=las.bicg(LO,edaFsparse.transpose()*fobs,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "mest[:,0] = q[0];\n",
    "\n",
    "# reorganize estimated model parameters back into grid\n",
    "Apre = np.zeros((I,J));\n",
    "for k in range(M):\n",
    "    i = floor(k/J);\n",
    "    j = k-i*J;\n",
    "    Apre[i,j] = mest[k];\n",
    "\n",
    "# plot results\n",
    "eda_draw('title p-obs',Aobs,'title p-pre',Apre);\n",
    "\n",
    "# compute weighted prediction error Ed\n",
    "dpre=np.zeros((N,1));\n",
    "for i in range(N):\n",
    "    k=dindex[i];\n",
    "    dpre[i,0] = mest[k,0];\n",
    "e = dobs-dpre;\n",
    "Ed = np.matmul( e.T, e );\n",
    "sigmadpos = sqrt(Ed/N);\n",
    "print('posterior errror in data',sigmadpos);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5564 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_05_10: make simulated (or \"synthetic\") pressure data\n",
    "\n",
    "# note: in order not to accidentally overwrite the data file\n",
    "# in the distribution, this code only writes the data to\n",
    "# mypressure.txt not pressure.txt; you must rename it to\n",
    "# have is used by eda05_08\n",
    "\n",
    "# rows of grid are x\n",
    "I = 41;\n",
    "Dx = 1;\n",
    "xmin = 0;\n",
    "xmax = Dx*(I-1);\n",
    "xobs = eda_cvec( np.linspace(xmin, xmax, I) );\n",
    "\n",
    "# columns of grid are y\n",
    "J = 41;\n",
    "Dy = 1;\n",
    "ymin = 0;\n",
    "ymax = Dy*(J-1);\n",
    "yobs = eda_cvec( np.linspace(ymin, ymax, J) );\n",
    "\n",
    "# randomly choose N of the I*J points to have data\n",
    "N = floor(I*J/10);\n",
    "rowindex = np.random.randint(I,size=N);\n",
    "xd = eda_cvec( xobs[rowindex,:] ); # x coords of data\n",
    "colindex = np.random.randint(J,size=N);\n",
    "yd = eda_cvec( yobs[colindex,:] ); # y coords of data\n",
    "\n",
    "# synthetoc data is chosen to be a solution to Laplace's equation\n",
    "kappa = 0.1;\n",
    "dtrue = eda_cvec(\n",
    "    10*np.multiply(np.sin(kappa*xd),np.exp(-kappa*yd))  );\n",
    "\n",
    "# synthetic observed data are true data plus random noise\n",
    "sigmad = 0.1;\n",
    "noise = np.random.normal(loc=0.0, scale=sigmad, size=(N,1));\n",
    "dobs = dtrue + noise;\n",
    "\n",
    "# populate a matrix A with data and display it\n",
    "A = np.zeros((I,J));\n",
    "for i in range(N):\n",
    "    A[rowindex[i],colindex[i]] = dobs[i];\n",
    "eda_draw('title pressure',A);\n",
    "\n",
    "# save matrix to file of tab-separated values\n",
    "# note that the filename has been changed to \"pressure2.txt\"\n",
    "# so that the original file \"pressure.txr\" is not overwritten\n",
    "# If you want to investigate the performance of eda05_08 with\n",
    "# different choices of data, change the filename in eda05_08 to\n",
    "# \"pressure2.txt\"\n",
    "DD = np.concatenate( (xd,yd,dobs),axis=1);  # concatenate (x,y,data) vectors into matrix\n",
    "np.savetxt(\"../data/mypressure.txt\",DD, delimiter=\"\\t\"); # write the matrix to a file\n",
    "\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
}
