ex3p1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Inverse Problems and Parameter Estimation, GEOS 627/427, University of Alaska Fairbanks\n",
    "\n",
    "- script ex3p1.ipynb\n",
    "- Aster Example 3.1\n",
    "- see also lib_geos.py functions svdmat() and svdall()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run lib_header.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from numpy.linalg import inv\n",
    "import scipy.linalg as la\n",
    "\n",
    "from lib_inverse import *  # plot_matrix, svdmat, svdall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# design matrix for linear forward model\n",
    "s2 = np.sqrt(2)\n",
    "G  = np.array([[1,0,0,1,0,0,1,0,0],\n",
    "            [0,1,0,0,1,0,0,1,0],\n",
    "            [0,0,1,0,0,1,0,0,1],\n",
    "            [1,1,1,0,0,0,0,0,0],\n",
    "            [0,0,0,1,1,1,0,0,0],\n",
    "            [0,0,0,0,0,0,1,1,1],\n",
    "            [s2,0,0,0,s2,0,0,0,s2],\n",
    "            [0,0,0,0,0,0,0,0,s2]],dtype='float')\n",
    "showmat(G,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find dimensions of G\n",
    "[ndata,nparm] = G.shape\n",
    "print('m = ndata = %i' % ndata)\n",
    "print('n = nparm = %i' % nparm)\n",
    "print('G is %i x %i' % (ndata,nparm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the singular values for the system matrix\n",
    "[U,svec,VH] = la.svd(G)\n",
    "S = la.diagsvd(svec,*G.shape)\n",
    "V = VH.T\n",
    "\n",
    "# Display singular values\n",
    "print('Singular values:')\n",
    "showmat([svec],3)\n",
    "\n",
    "# Find and display system rank\n",
    "print('System rank:')\n",
    "p = np.linalg.matrix_rank(G)\n",
    "print(p)\n",
    "\n",
    "# Display null space vectors \n",
    "print('Model null space vectors')\n",
    "showmat(V[:,p:nparm],2)\n",
    "\n",
    "# Display null space vectors reshaped to match tomography example geometry\n",
    "print('Model null space vectors reshaped into matrices')\n",
    "m01 = V[:,p].reshape((3,3))\n",
    "m02 = V[:,p+1].reshape((3,3))\n",
    "showmat(m01,2)\n",
    "showmat(m02,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot_matrix(m01,True,'gray')\n",
    "plt.clim([-0.5,0.5])\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "plt.title('Reshaped Null Space Model v8')\n",
    "print('Displaying image of null space model v8 (fig. 1)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display image of null space model V.,9\n",
    "plt.figure(2)\n",
    "plot_matrix(m02,True,'gray')\n",
    "plt.clim([-0.5,0.5])\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "plt.title('Reshaped Null Space Model v9')\n",
    "print('Displaying image of null space model v9 (fig. 2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display data space null vector\n",
    "print('Data null space vector')\n",
    "\n",
    "# Find and display model resolution matrix\n",
    "Vp = V[:,:p]\n",
    "Rm = Vp @ Vp.T\n",
    "\n",
    "# Display full model resolution matrix and colorbar\n",
    "plt.figure(3)\n",
    "plot_matrix(Rm,True,'gray')\n",
    "plt.clim([-0.1,1])\n",
    "plt.xlabel('k')\n",
    "plt.ylabel(\"k'\")\n",
    "plt.colorbar()\n",
    "#plt.xticks(range(9))\n",
    "#plt.yticks(range(9))\n",
    "print('Displaying full model resolution matrix (fig. 3)')\n",
    "plt.title('Model Resolution Matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rmdiag = np.diag(Rm).reshape((3,3))\n",
    "showmat(Rmdiag,2)\n",
    "\n",
    "plt.figure(4)\n",
    "plot_matrix(Rmdiag,True,'gray')\n",
    "plt.clim([-0.1,1])\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "plt.title('Reshaped Model Resolution Matrix Diagonal Elements')\n",
    "print('Displaying diagonal elements of model resolution matrix (fig. 4)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spike resolution test\n",
    "# Construct spike model\n",
    "mtest = np.zeros((nparm,1))\n",
    "mtest[4] = 1\n",
    "\n",
    "# Get noise free data for the spike model (forward problem)\n",
    "dtest = G @ mtest\n",
    "\n",
    "# Let python compute the psuedoinverse for us with pinv\n",
    "mdagger = np.linalg.pinv(G) @ dtest\n",
    "\n",
    "# Display reshaped noise free spike model with color bar\n",
    "plt.figure(5)\n",
    "plot_matrix(mtest.reshape((3,3)),True,'gray')\n",
    "plt.clim([-0.1,1])\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "plt.title('Spike Test Model')\n",
    "print('Displaying spike test model (fig. 5)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(6)\n",
    "plot_matrix(mdagger.reshape((3,3)),True,'gray')\n",
    "plt.clim([-0.1,1])\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "plt.title('Recovered Model from Spike Test Model Data')\n",
    "print('Displaying recovered model from noise-free spike test data (fig. 6)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a4e78903152c5c8c0a1611fb4b22aa78baa308a7734308b75819bd49ff5f4ec4"
  },
  "kernelspec": {
   "display_name": "inverse [conda env:inverse]",
   "language": "python",
   "name": "conda-env-inverse-py"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}