{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from math import exp, pi, sin, cos, tan, atan, sqrt, floor, log10, nan, log   # math functions\n",
    "import numpy as np\n",
    "import scipy.linalg as la             # linear algebra functions\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt  # general plotting \n",
    "from matplotlib.colors import ListedColormap\n",
    "import time;\n",
    "from matplotlib import animation, rc\n",
    "from IPython.display import HTML\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 gda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = gda_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 = gda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0:1] = v;  # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=gda_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=gda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def gda_cvec1(v):\n",
    "    if isinstance(v, int) or isinstance(v, np.int32):\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(\"gda_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) or isinstance(vi, np.int32): #patch v->vi 20230418\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(\"gda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(v); # patch v -> np.array(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) or isinstance(vi, np.int32): #patch v->vi 20230418\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(\"gda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"gda_cvec: %s not supported\" % type(v));\n",
    "\n",
    "        \n",
    "# gda_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 gda_draw(*argv):\n",
    "    DOCOLOR=True;\n",
    "    if( DOCOLOR ):\n",
    "        bwcmap = matplotlib.colormaps['jet'];\n",
    "    else:\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# useful tidbits\n",
    "\n",
    "# don't run this!\n",
    "xxx\n",
    "\n",
    "# read file\n",
    "myfile = \"planedata.txt\";\n",
    "D = np.genfromtxt(myfile, delimiter=\"\\t\");\n",
    "Nrows, Ncols = np.shape(D);\n",
    "\n",
    "# write file\n",
    "myfile = \"planedata.txt\";\n",
    "np.savetxt(myfile,np.concatenate((xobs,yobs,dobs),axis=1),delimiter=\"\\t\");\n",
    "\n",
    "# plot image\n",
    "mycmap = matplotlib.colormaps['jet'];\n",
    "plt.imshow(Dobs.T,cmap=mycmap,vmin=0.0,vmax=50.0);\n",
    "\n",
    "# create random integers\n",
    "ixr = np.random.randint(low=0, high=Nx, size=(Nr,1) );\n",
    "iyr = np.random.randint(low=0, high=Ny, size=(Nr,1) );\n",
    "\n",
    "# plot colored dots\n",
    "plt.scatter(xobs,yobs,c=dobs,cmap=mycmap,vmin=0.0,vmax=50.0,edgecolors=[0,0,0]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mtrue: 1.000  2.000  3.000\n",
      "dogmad: 5.000\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14016 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wrote file planedata.txt\n"
     ]
    }
   ],
   "source": [
    "# make test data for planar fit problem\n",
    "\n",
    "# x-axis\n",
    "Nx = 11;\n",
    "xmin = 0.0;\n",
    "xmax = 10.0;\n",
    "Dx = (xmax-xmin)/(Nx-1);\n",
    "x = gda_cvec( np.linspace(xmin,xmax,Nx) );\n",
    "# print(x);\n",
    "\n",
    "# y-axis\n",
    "Ny = 11;\n",
    "ymin = 0.0;\n",
    "ymax = 10.0;\n",
    "Dy = (ymax-ymin)/(Ny-1);\n",
    "y = gda_cvec( np.linspace(ymin,ymax,Ny) );\n",
    "# print(y);\n",
    "\n",
    "# choose model paramaters for d = m0 + m1 x + m2 y and print them\n",
    "mtrue = gda_cvec( 1.0, 2.0, 3.0 );\n",
    "print(\"mtrue: %.3f  %.3f  %.3f\" % (mtrue[0,0], mtrue[1,0], mtrue[2,0] ) );\n",
    "\n",
    "# evaluate plane on an (x,y) grid\n",
    "X = np.zeros((Nx,Ny));\n",
    "Y = np.zeros((Nx,Ny));\n",
    "Dtrue = np.zeros((Nx,Ny));\n",
    "for ix in range(Nx):\n",
    "    for iy in range(Ny):\n",
    "        X[ix,iy] = x[ix,0];\n",
    "        Y[ix,iy] = y[iy,0];\n",
    "        Dtrue[ix,iy] = mtrue[0,0] + mtrue[1,0]*x[ix,0] + mtrue[2,0]*y[iy,0];\n",
    "        \n",
    "# add noise to data\n",
    "sigmad = 5.0;\n",
    "print(\"dogmad: %.3f\" % (sigmad) );\n",
    "Dobs = Dtrue + np.random.normal( loc=0.0, scale=sigmad, size=(Nx,Ny) );\n",
    "\n",
    "# grab Nr random data points as the \"observed data\"\n",
    "Nr = 20;\n",
    "ixr = np.random.randint(low=0, high=Nx, size=(Nr,1) );\n",
    "iyr = np.random.randint(low=0, high=Ny, size=(Nr,1) );\n",
    "xobs = np.zeros((Nr,1));\n",
    "yobs = np.zeros((Nr,1));\n",
    "dobs = np.zeros((Nr,1));\n",
    "for ir in range(Nr):\n",
    "    xobs[ir,0] = X[ ixr[ir], iyr[ir] ];\n",
    "    yobs[ir,0] = Y[ ixr[ir], iyr[ir] ];\n",
    "    dobs[ir,0] = Dobs[ ixr[ir], iyr[ir] ];\n",
    "\n",
    "# plot plane\n",
    "fig1 = plt.figure();\n",
    "ax = plt.subplot(1,1,1);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.axis( [xmin, xmax, ymin, ymax] );\n",
    "mycmap = matplotlib.colormaps['jet'];\n",
    "if( 0 ): # checks orientation on 1\n",
    "    Dobs = np.zeros((Nx,Ny));\n",
    "    Dobs[0,0] = 1.0;\n",
    "    Dobs[Nx-1,0] = 2.0;\n",
    "    Dobs[0,Ny-1] = 3.0;\n",
    "plt.imshow(Dobs.T,cmap=mycmap,vmin=0.0,vmax=50.0);\n",
    "\n",
    "# overlay observed data locations\n",
    "plt.plot( xobs, yobs, 'ko' );\n",
    "plt.show();  \n",
    "\n",
    "# write random data points to file\n",
    "myfile = \"planedata.txt\";\n",
    "np.savetxt(myfile,np.concatenate((xobs,yobs,dobs),axis=1),delimiter=\"\\t\");\n",
    "print(\"wrote file %s\" % (myfile) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 20 nt 3 table from file planedata.txt\n",
      "estimated sigmad 4.427\n",
      "estimated solution\n",
      "m0 = 0.835 +/- 4.410 (95 pct)\n",
      "m1 = 1.976 +/- 0.723 (95 pct)\n",
      "m2 = 2.777 +/- 0.766 (95 pct)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 16168 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# least squares fit of plane to data\n",
    "\n",
    "# read (x,y,d) data from file\n",
    "myfile = \"planedata.txt\";\n",
    "D = np.genfromtxt(myfile, delimiter=\"\\t\");\n",
    "Nrows, Ncols = np.shape(D);\n",
    "print(\"read %d nt %d table from file %s\" % (Nrows,Ncols,myfile) );\n",
    "\n",
    "# rename columns xobs, yobs, dobs\n",
    "N = Nrows;\n",
    "M = 3; # 3 unkownns\n",
    "xobs = D[0:N,0:1];\n",
    "yobs = D[0:N,1:2];\n",
    "dobs = D[0:N,2:3];\n",
    "# print( np.concatenate( (xobs,yobs,dobs), axis=1 ) );\n",
    "\n",
    "# create G for planar model\n",
    "# di = m0 + xi * m1 + yi * m2\n",
    "G = np.concatenate( (np.ones((N,1)), xobs, yobs), axis=1 );\n",
    "# print(G);\n",
    "\n",
    "# standard least squares, mest = inv(GT G) GT dobs\n",
    "GTG = np.matmul(G.T,G);\n",
    "GTd = np.matmul(G.T,dobs);\n",
    "mest = la.solve( GTG, GTd );\n",
    "# print(mest);\n",
    "\n",
    "# predicted data dpre = G m\n",
    "dpre = np.matmul( G, mest );\n",
    "\n",
    "# error\n",
    "e = dobs - dpre;\n",
    "E = np.matmul( e.T, e );\n",
    "\n",
    "# estimate of standard deviation of data\n",
    "sigmad2 = E/(N-M);\n",
    "sigmad = sqrt(sigmad2);\n",
    "print(\"estimated sigmad %.3f\" % (sigmad) );\n",
    "\n",
    "# covariance\n",
    "Cm = sigmad2 * la.inv(GTG);\n",
    "\n",
    "print('estimated solution');\n",
    "for i in range(M):\n",
    "    print(\"m%d = %.3f +/- %.3f (95 pct)\" % (i,mest[i,0],2*sqrt(Cm[i,i])) );            \n",
    "                \n",
    "# x-axis\n",
    "Nx = 11;\n",
    "xmin = 0.0;\n",
    "xmax = 10.0;\n",
    "Dx = (xmax-xmin)/(Nx-1);\n",
    "x = gda_cvec( np.linspace(xmin,xmax,Nx) );\n",
    "# print(x);\n",
    "\n",
    "# y-axis\n",
    "Ny = 11;\n",
    "ymin = 0.0;\n",
    "ymax = 10.0;\n",
    "Dy = (ymax-ymin)/(Ny-1);\n",
    "y = gda_cvec( np.linspace(ymin,ymax,Ny) );\n",
    "# print(y);\n",
    "\n",
    "# evaluate dpre on (x,y) grid\n",
    "Dpre = np.zeros((Nx,Ny));\n",
    "for ix in range(Nx):\n",
    "    for iy in range(Ny):\n",
    "        Dpre[ix,iy] = mest[0,0] + mest[1,0]*x[ix,0] + mest[2,0]*y[iy,0];\n",
    "\n",
    "# plot predicted data as image\n",
    "fig1 = plt.figure();\n",
    "ax = plt.subplot(1,1,1);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.axis( [xmin, xmax, ymin, ymax] );\n",
    "mycmap = matplotlib.colormaps['jet'];\n",
    "if( 0 ): # checks orientation on 1\n",
    "    Dobs = np.zeros((Nx,Ny));\n",
    "    Dobs[0,0] = 1.0;\n",
    "    Dobs[Nx-1,0] = 2.0;\n",
    "    Dobs[0,Ny-1] = 3.0;\n",
    "plt.imshow(Dpre.T,cmap=mycmap,vmin=0.0,vmax=50.0);\n",
    "\n",
    "# overlay observed data\n",
    "plt.scatter(xobs,yobs,c=dobs,cmap=mycmap,vmin=0.0,vmax=50.0,edgecolors=[0,0,0]);\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
