ex3p2.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 ex3p2.ipynb\n",
    "- Aster Example 3.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run lib_header.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lib_inverse import *  # svdmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matlab code:\n",
    "# https://github.com/brianborchers/PEIP/blob/master/Examples/chap3/ex_3_2/ex_3_2.m\n",
    "#\n",
    "# CARL NOTE: The use of ndata/nparm and ndatad/nparmd is unfortunate.\n",
    "#            Maybe this can be fixed for future versions of the PEIP book code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set constants\n",
    "# The noise variance\n",
    "noise = 0.05\n",
    "\n",
    "# number of data and model parameters\n",
    "ndata = 210\n",
    "nparm = ndata\n",
    "\n",
    "# Discretizing values for nparm and ndata (210 data points)\n",
    "ndatad = ndata+1\n",
    "nparmd = nparm+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate time vector\n",
    "t = np.linspace(-5,100,nparmd)\n",
    "\n",
    "# Generate instrument impulse response as a critically-damped pulse\n",
    "# note: this is for plotting purposes only (only sigi and gmax are used)\n",
    "sigi = 10\n",
    "\n",
    "g = np.zeros(nparm)\n",
    "for i in range(nparm):\n",
    "    if (t[i] >= 0):\n",
    "        g[i] = t[i]*np.exp(-t[i]/sigi)\n",
    "\n",
    "# Normalize instrument response (i.e. max(g) = 1)\n",
    "g = np.array(g)\n",
    "gmax = np.max(g)    # 3.6788\n",
    "g = g/gmax\n",
    "\n",
    "# Plot of instrument response to unit area ground acceleration impulse.\n",
    "plt.figure(2)\n",
    "plt.plot(t[:nparm],g,'k')\n",
    "plt.axis('tight')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('V')\n",
    "print('Instrument response to unit ground acceleration impulse (fig. 2)')\n",
    "#print -deps2 c3fimp_resp.eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Populate G matrix \n",
    "# First the numerator which varies\n",
    "G = np.zeros((ndata,nparm))\n",
    "for i in range(1,ndatad):\n",
    "    for j in range(nparm):\n",
    "        tp = t[j]-t[i]\n",
    "        if (tp > 0):\n",
    "            G[i-1,j] = 0\n",
    "        else:\n",
    "            G[i-1,j] = -tp*np.exp(tp/sigi)\n",
    "\n",
    "# now divide everything by the denominator\n",
    "deltat =  t[1]-t[0]\n",
    "G = G/gmax * deltat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get SVD of G matrix\n",
    "[U,S,V] = svdmat(G)\n",
    "#[U, s, VH] = np.linalg.svd(G) \n",
    "#S=scipy.linalg.diagsvd(s,*G.shape)                  # vector of singular values\n",
    "#V=VH.T\n",
    "\n",
    "# Display image of G matrix with colorbar\n",
    "plt.figure(1)\n",
    "plt.imshow(G,cmap='gray')\n",
    "plt.colorbar\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.title('G matrix for convolution problem')\n",
    "print('Displaying image of G matrix (fig. 1)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True signal is two pulses of sig deviation\n",
    "sig = 2\n",
    "# Find unscaled true signal model  \n",
    "mtrue = np.exp(-(t[:nparm]-8)**2/(sig**2*2)).T + 0.5*np.exp(-(t[:nparm]-25)**2/(sig**2*2)).T\n",
    "# Rescale true signal model to have a max of 1\n",
    "mtrue = mtrue/np.max(mtrue)\n",
    "# Get true data without noise\n",
    "d = G@mtrue.reshape((len(mtrue),1))\n",
    "# Add random normal noise to the data\n",
    "dn = G@mtrue.reshape((len(mtrue),1)) + noise*np.random.randn(ndata,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using SVD with all 210 singular values\n",
    "nkeep = nparm\n",
    "# Find Up, Vp, Sp\n",
    "Up = U[:,:nkeep]\n",
    "Vp = V[:,:nkeep]\n",
    "Sp = S[:nkeep,:nkeep]\n",
    "\n",
    "# Generalized inverse solutions for noisy data (mn) \n",
    "# and noise-free data (mperf)\n",
    "Sinv = la.inv(Sp)\n",
    "mn    = Vp@Sinv@Up.T @ dn\n",
    "mperf = Vp@Sinv@Up.T @ d\n",
    "\n",
    "# Display semilog plot of singular values\n",
    "plt.figure(3)\n",
    "plt.semilogy(np.diag(S),'ko')\n",
    "plt.axis('tight')\n",
    "plt.xlabel('k')\n",
    "plt.ylabel(r'$s_k$')\n",
    "print('Displaying semilog plot of singular values (fig. 3)')\n",
    "\n",
    "# Plot true model\n",
    "plt.figure(4)\n",
    "plt.plot(t[:nparm],mtrue,'k')\n",
    "plt.xlim([-5,100])\n",
    "plt.ylim([0,1])\n",
    "plt.axis('tight') \n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Acceleration (m/$s^2$)')\n",
    "print('Displaying true model (fig. 4)')\n",
    "#print -deps2 c3fm_true.eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display predicted data using noise free model\n",
    "plt.figure(5)\n",
    "plt.plot(t[:nparm],d,'k')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('V')\n",
    "plt.axis('tight')\n",
    "print('Displaying predicted data from true model (without noise) (fig. 5)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display predicted data plus random independent noise\n",
    "plt.figure(6)\n",
    "plt.plot(t[:nparm],dn,'k')\n",
    "plt.xlim([-5,100])\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('V')\n",
    "plt.axis('tight')\n",
    "print('Displaying predicted data from true model plus independent noise (fig. 6)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display generalized inverse solution for noise-free data\n",
    "plt.figure(7)\n",
    "plt.plot(t[:nparm],mperf,'k')\n",
    "plt.xlim([-5,100])\n",
    "plt.ylim([0,1])\n",
    "plt.xlabel('Time (s)')\n",
    "plt.axis('tight')\n",
    "plt.ylabel('Acceleration (m/$s^2$)')\n",
    "print('Displaying generalized inverse solution for noise-free data (210 singular values) (fig. 7)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(8)\n",
    "plt.plot(t[:nparm],mn,'k')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.axis('tight')\n",
    "plt.ylabel('Acceleration (m/$s^2$)')\n",
    "print('Displaying generalized inverse solution for noisy data (210 singular values) (fig. 8)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Truncate SVD to 26 singular values\n",
    "nkeep = 26\n",
    "Up = U[:,:nkeep]\n",
    "Vp = V[:,:nkeep]\n",
    "Sp = S[:nkeep,:nkeep]\n",
    "\n",
    "# Get model for truncated SVD (m2) with noisy data\n",
    "Sinv = la.inv(Sp)\n",
    "m2 = Vp@Sinv@Up.T @ dn\n",
    "\n",
    "# Display generalized inverse solution for noisy data\n",
    "plt.figure(9)\n",
    "plt.plot(t[:nparm],m2,'k')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.axis('tight')\n",
    "plt.ylabel('Acceleration (m/$s^2$)')\n",
    "print('Displaying generalized inverse solution for noisy data (26 singular values) (fig. 9)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get resolution matrix\n",
    "Rm = Vp@Vp.T\n",
    "# Display image of resolution matrix for truncated solution with colorbar\n",
    "plt.figure(10)\n",
    "\n",
    "# use tight bounds on the colorbar instead of automatically expanded bounds\n",
    "plt.imshow(Rm,cmap='gray')\n",
    "plt.xlabel('j')\n",
    "plt.ylabel('i')\n",
    "plt.colorbar()\n",
    "print('Displaying image of resolution matrix for truncated SVD solution (26 singular values) (fig. 10)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(11)\n",
    "plt.plot(t[:nparm], Rm[79,:],'k')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.axis('tight')\n",
    "plt.ylabel('Element Value')\n",
    "print(['Displaying column of model resolution matrix for truncated SVD solution (fig. 11)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#show successive TSVD solutions\n",
    "print('Animating TSVD solutions as p increases (fig. 12) [THIS MAY TAKE AWHILE]')\n",
    "\n",
    "# the maximum p to use, the fit model and the diagonal elements\n",
    "p = np.linalg.matrix_rank(G)\n",
    "m = np.zeros((nparm,1))\n",
    "ss = np.diag(S)\n",
    "r = []\n",
    "mnorm = []\n",
    "\n",
    "prun = p\n",
    "#prun = 30\n",
    "\n",
    "plt.figure(12)\n",
    "for i in range(prun):\n",
    "    # adjust the predicted model to have p singular values\n",
    "    m = m + (U[:,i].T@dn/ss[i])*V[:,i].reshape((nparm,1))\n",
    "    # keep track of the residuals for each p\n",
    "    r.append(np.linalg.norm(G@m-dn))\n",
    "    # keep track of the model norm for each p\n",
    "    mnorm.append(np.linalg.norm(m))\n",
    "    \n",
    "    # plot the newly fit model\n",
    "    plt.cla()\n",
    "    plt.plot(m)\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Acceleration (m/$s^2$)')\n",
    "    plt.axis([0,p+1,-1,1.5])\n",
    "    # plot the true model\n",
    "    plt.plot(mtrue,'r-.',lw=2)\n",
    "    plt.title(i+1)\n",
    "    #plt.draw()\n",
    "    plt.pause(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPTION A: see code above\n",
    "# OPTION B: use tsvd(), which returns norm-squared calculations\n",
    "f_r, rss, f_r_ss = tsvd(dn,G,np.arange(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Examine the trade-off curve (collected in the loop above)\n",
    "# (The L curve is introduced in Chapter 4.)\n",
    "plt.figure(13)\n",
    "plt.plot(r,mnorm,'b.-')\n",
    "plt.plot(np.sqrt(rss),np.sqrt(f_r_ss),'r--')  # note square root\n",
    "plt.plot(r[nkeep-1],mnorm[nkeep-1],'ro')\n",
    "plt.ylabel(r'$||m||_2$')\n",
    "plt.xlabel(r'$||Gm-d||_2$')\n",
    "plt.text(r[nkeep-1], mnorm[nkeep-1], f'{nkeep} singular values in TSVD')\n",
    "print('Displaying TSVD trade-off curve (fig. 13)')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "7f341b321178e62fbbfacfb91167ed733f9907f6f0692efdceb8394e8303ba73"
  },
  "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
}