covC_plot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "33eb2e052a6246fb",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "# Inverse Problems and Parameter Estimation, GEOS 627/427, University of Alaska Fairbanks\n",
    "\n",
    "- script covC_plot.ipynb\n",
    "- plotting covariance functions used in GEOS 627"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f77600021dbe2974",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70e0083",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from lib_inverse import covC\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd51924a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample distances\n",
    "sigma  = 4\n",
    "Lprime = 22\n",
    "dmax   = 2.5*Lprime\n",
    "n      = 1e4\n",
    "dcont  = np.linspace(0,dmax,num=int(n))\n",
    "\n",
    "# some reference heights (horizontal lines) of e-folding depths\n",
    "print('reference lines at y = exp(-efold)/sigma^2):')\n",
    "#efold = 2**np.arange(start=-3,stop=1.5,step=0.5)   # LFACTOR = 1 (Tarantola)\n",
    "efold = 2**np.arange(start=-1,stop=1.5,step=0.5)   # LFACTOR = 2 (preferred)\n",
    "print(np.exp(-efold),'\\n')\n",
    "\n",
    "# some reference lengths (vertical lines)\n",
    "print('reference lines at x = dref:')\n",
    "yfold = sigma**2*np.exp(-efold)\n",
    "#[xmat,ymat] = horzlines(yfold,axtemp(1),axtemp(2));\n",
    "dref = Lprime*np.array([1, 1/np.sqrt(2), 1/2])\n",
    "print(dref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6401b00d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "# cyan line shows L as defined in Tarantola 2005\n",
    "# red line shows L' = 2L\n",
    "nr = 2\n",
    "nc = 2\n",
    "\n",
    "plt.subplots(nrows=nr,ncols=nc,figsize=(10,12))\n",
    "\n",
    "stcovs = ['Gaussian','Exponential','Circular']\n",
    "axtemp = [0, dmax, sigma**2*-0.1, sigma**2*1.1]\n",
    "\n",
    "for icov in range(3):\n",
    "    Cd = covC(dcont,[icov+1,Lprime,sigma])\n",
    "    \n",
    "    plt.subplot(nr,nc,icov+1)\n",
    "    plt.plot(axtemp[0:2],[0,0],'k')\n",
    "    plt.plot([dref[0],dref[0]], [axtemp[2],axtemp[3]],'r')\n",
    "    plt.plot([dref[1],dref[1]], [axtemp[2],axtemp[3]],'g')\n",
    "    plt.plot([dref[2],dref[2]], [axtemp[2],axtemp[3]],'c')\n",
    "    \n",
    "    Cplot = covC(dref[0],[icov+1,Lprime,sigma])\n",
    "    plt.plot(axtemp[1:3],[Cplot,Cplot],'r')\n",
    "    \n",
    "    Cplot = covC(dref[1],[icov+1,Lprime,sigma])\n",
    "    plt.plot(axtemp[1:3],[Cplot,Cplot],'g')\n",
    "    \n",
    "    Cplot = covC(dref[2],[icov+1,Lprime,sigma])\n",
    "    plt.plot(axtemp[1:3],[Cplot,Cplot],'c')\n",
    "    \n",
    "    plt.hlines(yfold,axtemp[0],axtemp[1],colors='k',linestyle='dashed')\n",
    "    \n",
    "    #plt.plot(axtemp[0:3],[sigma**2*.np.exp(-2),sigma**2*.np.exp(-2)],'k--')   # check\n",
    "    plt.plot(dcont,Cd,'b',lw=2)\n",
    "    \n",
    "    if icov == 0:\n",
    "        plt.plot(dcont,covC(dcont,[2,Lprime,sigma]),'b--',lw=1)\n",
    "        \n",
    "        for ii in range(len(yfold)):\n",
    "            txt = 'exp(-%.2f)' % efold[ii]\n",
    "            plt.text(0.75*dmax,yfold[ii]+0.5,txt)\n",
    "            \n",
    "    elif icov == 1:\n",
    "        plt.plot(dcont,covC(dcont,[1,Lprime,sigma]),'b--',lw=1)\n",
    "        \n",
    "    plt.xlabel('Distance')\n",
    "    plt.title('%s covariance' % stcovs[icov])\n",
    "    plt.axis(axtemp)\n",
    "    \n",
    "# Matern covariance   \n",
    "# note: nu=0.5 is exponential; nu->infty is Gaussian\n",
    "# note: there may be numerical artifacts associated with besselk and gamma\n",
    "\n",
    "icov = 4\n",
    "nuvec = [0.25, 0.5, 1.5, 100]\n",
    "stit = 'Matern covariance\\nnu = %.1f, %.1f, %.1f, %.1f' % (nuvec[0],nuvec[1],nuvec[2],nuvec[3])\n",
    "\n",
    "plt.subplot(nr,nc,4)\n",
    "for ii in range(len(nuvec)):\n",
    "    nu = nuvec[ii]\n",
    "    Cd = covC(dcont,[icov,Lprime,sigma,nu])\n",
    "    \n",
    "    plt.plot(axtemp[0:2],[0,0],'k')\n",
    "    plt.plot([Lprime, Lprime],[axtemp[2],axtemp[3]],'r')\n",
    "    ##plt.plot([Lprime/np.sqrt(2), Lprime/np.sqrt(2)],[axtemp[2],axtemp[3]],'g')\n",
    "    ##plt.plot([Lprime/2, Lprime/2],[axtemp[2],axtemp[3]],'c')\n",
    "    \n",
    "    Cplot = covC(Lprime,[icov,Lprime,sigma,nu])\n",
    "    plt.plot(axtemp[0:2],[Cplot, Cplot],'r')\n",
    "    \n",
    "    ##Cplot = covC(Lprime/np.sqrt(2),[icov,Lprime,sigma,nu])\n",
    "    ##plt.plot(axtemp[0:2],[Cplot, Cplot],'g')\n",
    "    \n",
    "    ##Cplot = covC(Lprime/2,[icov,Lprime,sigma,nu])\n",
    "    ##plt.plot(axtemp[0:2],[Cplot, Cplot],'c')\n",
    "    \n",
    "    #plt.plot(axtemp[0:2],[sigma**2*np.exp(-2),sigma**2*np.exp(-2),'k--')    # check\n",
    "    plt.hlines(yfold,axtemp[0],axtemp[1],colors='k',linestyle='dashed') #check\n",
    "    plt.plot(dcont,Cd,'b',lw=1)\n",
    "    plt.xlabel('Distance')\n",
    "    plt.axis(axtemp)\n",
    "    plt.title(stit)\n",
    "    \n",
    "plt.subplots_adjust(wspace=0.3,hspace=0.3)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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": 5
}