lab_iter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "# Applied Seismology, GEOS 626\n",
    "# Inverse Problems and Parameter Estimation, GEOS 627\n",
    "\n",
    "- genlsq = Generalized Least Squares\n",
    "- template script for the iterative quasi-Newton method for a 4-parameter inversion for epicenter (xs, ys),\n",
    "origin time (ts), and velocity (V).\n",
    "- the algorithm employs generalized least squares, whereby both data covariances and model covariances are\n",
    "used.\n",
    "- calls forward_epicenter.ipynb\n",
    "- background reading: Tarantola book (2005), Ch. 3 and Appendix 6.22"
   ]
  },
  {
   "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": [
    "import matplotlib\n",
    "import math\n",
    "import scipy.linalg as la\n",
    "from scipy.linalg import inv\n",
    "#from numpy.linalg import inv\n",
    "import numpy.matlib\n",
    "#from scipy.stats.distributions import chi2\n",
    "\n",
    "from lib_inverse import *  # plot_ellipse\n",
    "from lib_peip import *     # chi2inv\n",
    "from lib_geos import *     # plot_histo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# USER INPUT\n",
    "nsamples = 1000                # for visualizing covariance matrices (Cm and Cd)\n",
    "irandom_initial_model = 0      # 0(fixed), 1(random)\n",
    "irandom_target_model = 0       # 0(fixed), 1(random)\n",
    "idata_errors = 2               # 0(none),  1(random), 2(fixed)\n",
    "ifig = 1                       # 0,1\n",
    "\n",
    "inormalization = 1\n",
    "stnsamples = str(nsamples) + ' samples'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FORWARD PROBLEM: EXPLORE THIS SEPARATE NOTEBOOK   \n",
    "%run ./forward_epicenter.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ifig==1:\n",
    "    # plot different histograms of properties of the prior model covariance samples\n",
    "    fig = plt.figure(figsize=(9,9))\n",
    "    for kk in range(nparm):\n",
    "        sigma = sigma_prior[kk]\n",
    "        edges = np.arange(-4*sigma,4*sigma,sigma/2)\n",
    "        etemp = cov_samples_m[kk,:]\n",
    "        plt.subplot(2,2,kk+1)\n",
    "        plt.grid()\n",
    "        plot_histo(etemp,edges)\n",
    "        #plt.bar(etemp)\n",
    "        plt.ylim([0, 0.4])\n",
    "        plt.title('mprior samples: Model parameter '+str(kk+1)+' ('+str(mlabs[kk])+')\\n'+'mean = %.5f; std = %.5f' % (np.mean(etemp),np.std(etemp)))\n",
    "        plt.tight_layout()\n",
    "    # plot different histograms of properties of the data covariance samples\n",
    "    fig2 = plt.figure(figsize=(9,11))\n",
    "    for ii in range(ndata):\n",
    "        sigma = sigma_obs[ii]\n",
    "        edges = np.arange(-4*sigma,4*sigma,sigma/2)\n",
    "        etemp = cov_samples_d[ii,:]\n",
    "        plt.subplot(4,3,ii+1)\n",
    "        plot_histo(etemp,edges)\n",
    "        plt.ylim([0, 0.4])\n",
    "        plt.title('Data index ' + str(ii+1) +'\\n'+'mean = %.5f; std = %.5f' % (np.mean(etemp),np.std(etemp))) \n",
    "        plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# predictions for prior and initial models (not necessary)\n",
    "dprior   = d(mprior)\n",
    "dinitial = d(minitial)\n",
    "\n",
    "if ifig==1:\n",
    "    fig3 = plt.figure(figsize=(10,10))\n",
    "    ivecd = np.arange(1,ndata+1)\n",
    "    plt.plot(ivecd,dobs_samples,'.-')\n",
    "    #p1 = plt.plot(ivecd,dprior,'bo-',linewidth=2,markersize=10,markerfacecolor='b',markeredgecolor='w')\n",
    "    p2 = plt.plot(ivecd,dinitial,'ko-',linewidth=2,markersize=10,markerfacecolor='k',markeredgecolor='w')\n",
    "    p3 = plt.plot(ivecd,dtarget,'ro--',linewidth=2,markersize=10,markerfacecolor='r',markeredgecolor='w')\n",
    "    p4 = plt.plot(ivecd,dobs,'ro-',linewidth=2,markersize=10,markerfacecolor='r',markeredgecolor='w')\n",
    "    #plt.legend([p1[0],p2[0],p3[0],p4[0]],['g(mprior)','g(minitial)','g(mtarget)','g(mtarget) + errors','location','northwest'])\n",
    "    plt.legend([p2[0],p3[0],p4[0]],['d(minitial)','d(mtarget)','d(mtarget) + errors','location','northwest'])\n",
    "    plt.title('BLACK = d(minitial); RED DASHED = d(mtarget); RED = d(mtarget) + errors');\n",
    "    plt.xlim([0.5, ndata+0.5])  #set(gca,'xtick',[1:ndata]);\n",
    "    plt.xlabel('Data index')\n",
    "    plt.ylabel('Prediction value, d = g(m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# function to display model parameters\n",
    "def printm(nn,niter,mprior,mnew,mtarget):\n",
    "    print('%i/%i : %10s%16s%16s' % (nn,niter,'prior', 'current', 'target'))\n",
    "    for ii in range(len(mprior)):\n",
    "        print('%16.6f%16.6f%16.6f' % (mprior[ii],mnew[ii],mtarget[ii]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MISFIT FUNCTION: least squares, Tarantola (2005), Eq. 6.251\n",
    "# This calls the function d(m) to compute the predictions.\n",
    "\n",
    "# data misfit\n",
    "def Sd(m,dobs,icobs):\n",
    "    sd = 0.5 * (d(m)-dobs).T @ icobs @ (d(m)-dobs)\n",
    "    return sd\n",
    "# model misfit (related to regularization)\n",
    "def Sm(m,mprior,icprior):\n",
    "    sm = 0.5 * (m-mprior).T @ icprior @ (m-mprior)\n",
    "    return sm\n",
    "# total misfit\n",
    "def S(m,dobs,mprior,icobs,icprior):\n",
    "    s = Sd(m,dobs,icobs) + Sm(m,mprior,icprior)\n",
    "    return s\n",
    "\n",
    "# initial model\n",
    "#mnew = mprior     # prior model\n",
    "mnew = minitial\n",
    "\n",
    "Sd_0 = Sd(mnew,dobs,icobs)\n",
    "Sm_0 = Sm(mnew,mprior,icprior)\n",
    "S_0  = S(mnew,dobs,mprior,icobs,icprior)\n",
    "stS0 = 'S(m0) = %.3f = %.3f(D) + %.3f(M)' % (S_0,Sd_0,Sm_0)\n",
    "print(stS0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "niter = 0\n",
    "#niter = input(' Select the number of iterations (< 10) or 0 to exit: ')\n",
    "#niter = int(niter)\n",
    "if int(niter) == 0:\n",
    "    raise Exception('no iterations requested: exiting')\n",
    "    #sys.exit('no iterations requested: exiting')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# display information for initial model\n",
    "nn = 0\n",
    "printm(nn,niter,mprior,mnew,mtarget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize arrays\n",
    "iter_vec = np.transpose(np.arange(0,niter+1))\n",
    "Sd_vec = np.zeros((niter+1,1))\n",
    "Sm_vec = np.zeros((niter+1,1))\n",
    "S_vec = np.zeros((niter+1,1))\n",
    "\n",
    "# misfit for initial model\n",
    "Sd_vec[0] = Sd_0\n",
    "Sm_vec[0] = Sm_0\n",
    "S_vec[0]  = S_0\n",
    "\n",
    "for nn in range(1,niter+1):\n",
    "    #///////////////////////////////\n",
    "    # CODE HERE for quasi-Newton algorithm\n",
    "    print('code here')\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    # misfit function for new model\n",
    "    # note: book-keeping only -- not used within the algorithm above\n",
    "    #Sd_vec[nn] = \n",
    "    #Sm_vec[nn] = \n",
    "    #S_vec[nn] = \n",
    "\n",
    "    printm(nn,niter,mprior,mnew,mtarget)\n",
    "    \n",
    "    #///////////////////////////////\n",
    "\n",
    "# misfit function values\n",
    "print('summary of misfit function:');\n",
    "print('%8s%16s%16s%16s' % ('iter','Sd','Sm','S = Sm + Sd'))\n",
    "for nn in range(niter):\n",
    "    print('%8i%16.10f%16.10f%16.10f' % (iter_vec[nn],Sd_vec[nn],Sm_vec[nn],S_vec[nn]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot convergence curve\n",
    "ya = math.floor(np.min(np.log10([Sm_vec, Sd_vec, S_vec])))\n",
    "yb = math.ceil(np.max(np.log10([Sm_vec, Sd_vec, S_vec])))\n",
    "ylims = [10**ya, 10**yb]\n",
    "\n",
    "plt.plot(iter_vec,np.log10(S_vec),'k.-',iter_vec,np.log10(Sm_vec),'b.-',iter_vec,np.log10(Sd_vec),'r.-',\n",
    "    linewidth=2,markersize=20)\n",
    "plt.legend(('S(mⁿ) = Sd + Sm','Sm(mⁿ)','Sd(mⁿ)'),fontsize=14)\n",
    "plt.xlim([-0.5, niter+0.5])\n",
    "#plt.ylim([np.log10(ylims[0]),np.log10(ylims[0])])\n",
    "\n",
    "plt.locator_params(axis=\"x\", integer=True, tight=True)\n",
    "plt.xlabel('n, iteration')\n",
    "plt.ylabel(' log10[ S(mⁿ) ], misfit function')\n",
    "plt.title(str(niter) + ' iterations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exit early\n",
    "#raise SystemExit(\"stop here\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the correlation matrix from the covariance matrix\n",
    "def corrcov(C):\n",
    "    nx,ny = C.shape\n",
    "    if nx != ny:\n",
    "        return\n",
    "    sigma = np.sqrt(np.diag(C))\n",
    "    outer_v = np.outer(sigma,sigma)\n",
    "    Crho = C / outer_v\n",
    "    Crho[C == 0] = 0\n",
    "    return Crho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#///////////////////////////////\n",
    "# COMPUTE THE FOLLOWING\n",
    "# mpost       posterior model (\"final\" model)\n",
    "# dpost       predictions for mpost\n",
    "# Gpost       partial derivatives matrix at mpost\n",
    "# cpost0      posterior covariance matrix (use icobs0 and icprior0)\n",
    "# sigma_post  variances of the posterior covariance matrix\n",
    "# rho_post    posterior correlation matrix (hint: see Tarantola, Section 3.3)\n",
    "# CODE HERE\n",
    "\n",
    "\n",
    "\n",
    "#///////////////////////////////"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a priori model correlations (for comparison)\n",
    "rho_prior = corrcov(cprior0)\n",
    "\n",
    "print('mpost \\n',mpost, '\\n', 'cpost0 \\n', cpost0, '\\n','rho_post \\n', rho_post)\n",
    "print('variances:')\n",
    "print(np.diag(cpost0))\n",
    "print('sigma_post \\n',sigma_post)\n",
    "\n",
    "# posterior data covariance matrix (e.g., Tarantola Eq. 3.44)\n",
    "\n",
    "cpost0_d     = Gpost @ cpost0 @ Gpost.T\n",
    "cpost0_d     = (cpost0_d + cpost0_d.T)/2   # force to be symmetric\n",
    "sigma_post_d = np.sqrt(np.diag(cpost0_d))  # IGNORING OFF-DIAGONAL ELEMENTS\n",
    "rho_post_d   = corrcov(cpost0_d)           # posterior correlation matrix\n",
    "rho_prior_d  = corrcov(cobs0)              # prior, for comparison\n",
    "\n",
    "# probably cannot get a cholesky decomposition for Cpost_d\n",
    "# Lpost_d = chol(cpost0_d,'lower');\n",
    "# dpost_samples = zeros(ndata,nsamples);\n",
    "# dcov_samples = zeros(ndata,nsamples);\n",
    "# for xx=1:nsamples, randn_vecs_d(:,xx) = randn(ndata,1); end\n",
    "# dcov_samples  = Lpost_d * randn_vecs_d;\n",
    "# dpost_samples = repmat(dpost,1,nsamples) + dcov_samples;\n",
    "# std_samples_d = std(dpost_samples');\n",
    "\n",
    "#format long\n",
    "print('model summary (%i iterations):' % (niter))\n",
    "print('%16s%16s%16s%16s' % ('prior', 'initial','posterior', 'target'))\n",
    "for ii in range(nparm):\n",
    "    print('%16s%16s%16s%16s' % (str(mprior[ii]),str(minitial[ii]),str(mpost[ii]),str(mtarget[ii])))\n",
    "print('data summary (%i observations):' % (ndata))\n",
    "print('%16s%16s%16s%16s%16s' % ('prior', 'initial','posterior', 'target', 'actual'))\n",
    "for ii in range(ndata):\n",
    "    print('%16s%16s%16s%16s%16s' % (str(dprior[ii]),str(dinitial[ii]),str(dpost[ii]), \n",
    "                                    str(dtarget[ii]), str(dobs[ii])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cholesky decomposition to obtain the square-root of cpost0\n",
    "# NOTE: for large problems, this is not possible due to poor\n",
    "#       conditioning of cpost0 or the inability to compute cpost0\n",
    "Lpost = np.linalg.cholesky(cpost0)\n",
    "\n",
    "# samples of the posterior distribution\n",
    "mpost_samples = np.zeros((nparm,nsamples))\n",
    "mcov_samples  = np.zeros((nparm,nsamples))\n",
    "for xx in range(nsamples): \n",
    "    randn_vecs_m[:,xx] = np.random.randn(nparm)\n",
    "mcov_samples  = Lpost @ randn_vecs_m\n",
    "mpost_samples = np.matlib.repmat(mpost,1,nsamples) + mcov_samples\n",
    "\n",
    "# compare the standard deviation with sigma_post\n",
    "std_samples = np.std(mpost_samples.T, axis=0)\n",
    "\n",
    "# compare posterior model distribution with prior\n",
    "print('  ')\n",
    "print(' Compare model uncertainties : ')\n",
    "#print('              model parameter: ' + ''.join([f\"{_:>13}\" for _ in mlabs]))  \n",
    "print('             model parameter :',''.join(['%13s' % m for m in mlabs]))\n",
    "print('                       units :',''.join(['%13s' % u for u in ulabs]))\n",
    "print('                 sigma_prior =',''.join(['%13.5s' % u for u in sigma_prior]))\n",
    "print('                  sigma_post =',''.join(['%13.5s' % u for u in sigma_post]))\n",
    "print('   std(%6.0f mpost_samples) =' % nsamples,''.join(['%13.5s' % u for u in std_samples]))\n",
    "print('      sigma_post/sigma_prior =',''.join(['%13.5s' % u for u in np.divide(sigma_post,sigma_prior)]))\n",
    "print('1 - (sigma_post/sigma_prior) =',''.join(['%13.5s' % u for u in 1 - np.divide(sigma_post,sigma_prior)]))           \n",
    "\n",
    "# compute the predictions associated with the posterior samples,\n",
    "# then compare std_d_samples with sigma_post_d\n",
    "d_samples = np.zeros((ndata,nsamples))\n",
    "for xx in range(nsamples):\n",
    "    ms = mpost_samples[:,xx]\n",
    "    d_samples[:,xx] = d(ms).flatten()\n",
    "\n",
    "covd_samples  = np.cov(d_samples)\n",
    "rhod_samples  = corrcov(covd_samples)\n",
    "std_d_samples = np.sqrt(np.diag(covd_samples))\n",
    "#std_d_samples = std(d_samples')\n",
    "\n",
    "print('  ')\n",
    "print('Compare data uncertainties : ')\n",
    "print('%16s %10s %10s %10s' % ('prior','post','samples','post/prior'))\n",
    "call = [sigma_obs, sigma_post_d, std_d_samples, np.divide(sigma_post_d,sigma_obs)]\n",
    "for ii in range(ndata):\n",
    "    print('%6i%10.4f %10.4f %10.4f %10.4f' % (ii+1,call[0][ii],call[1][ii],call[2][ii],call[3][ii]))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_covsamples(msamples,rho,tlab,sym1,msamples2,rho2,tlab2,sym2,mlabs):\n",
    "    # Python adaptation of plot_covsamples\n",
    "    #PLOT_COVMSAMPLES generates plots for samples of covariance matrix\n",
    "    #\n",
    "    # INPUT\n",
    "    #    msamples   n x p matrix of vector samples\n",
    "    #    rho        n x n 'analytical' correlation matrix\n",
    "    #    tlab       label for plot\n",
    "    #    sym1       plotting symbol\n",
    "    #    msamples2  optional: 2nd set of samples ([] for none)\n",
    "    #    rho2       optional: 2nd 'analytical' correlation matrix ([] for none)\n",
    "    #    tlab2      optional: label for plot\n",
    "    #    sym2       optional: plotting symbol\n",
    "    #    mlabs      optional: labels for each variable ([] for default)\n",
    "    #\n",
    "    # EXAMPLE: \n",
    "    #    plot_covsamples(mpost_samples,rho_post,'mpost','b.',[],[],[],[],mlabs);\n",
    "    #\n",
    "    # NOTE: We could alternatively estimate the covariance matrix\n",
    "    # (and correlation matrix) directly from the input samples.\n",
    "    #\n",
    "    # Carl Tape, 2012-01-01\n",
    "    #\n",
    "    NMAX = 6   # max number to make into scatterplot\n",
    "    \n",
    "    n = np.shape(msamples)[0]\n",
    "    p = np.shape(msamples)[1]\n",
    "    print('plot_covsamples: n = %i, p = %i' % (n,p))\n",
    "    \n",
    "    if len(mlabs)==0:\n",
    "        #mlabs = strtrim(cellstr(num2str([1:n]')))\n",
    "        mlabs = []    #np.matlib.repmat(str(''),n,1)\n",
    "        for ii in range(n):\n",
    "            mlabs.append( 'i'+str(ii+1))        \n",
    "    # whether to plot a second set of samples\n",
    "    if len(msamples2)!=0 and len(rho2)!=0:\n",
    "        iplot_second_set = 1\n",
    "    else:\n",
    "        iplot_second_set = 0\n",
    "    \n",
    "    # kk=0: correlation matrices from Cpost\n",
    "    # kk=1: correlation matrices based on input SAMPLES\n",
    "    fig = plt.figure(figsize=(8,8))\n",
    "    nr = 1+iplot_second_set\n",
    "    nc = 2\n",
    "    for kk in [0,1]:\n",
    "        if kk==0:\n",
    "            F1 = rho\n",
    "            if iplot_second_set==1:\n",
    "                F2 = rho2\n",
    "            stag = ''\n",
    "        else:\n",
    "            F1 = np.corrcoef(msamples)\n",
    "            if iplot_second_set==2:\n",
    "                F2 = np.corrcoef(msamples2)\n",
    "            stag = 'sample'\n",
    "        \n",
    "        # first matrix plot\n",
    "        pind = kk+1+iplot_second_set*(kk)\n",
    "        ax = fig.add_subplot(nr,nc,pind)\n",
    "        p1=ax.imshow(F1,cmap='jet',vmin=-1,vmax=1)\n",
    "        ax.xaxis.tick_top()\n",
    "        if iplot_second_set==0:\n",
    "            fig.colorbar(p1,shrink=0.3,aspect=10,ticks=[-1,-0.5,0,0.5,1])\n",
    "        else:\n",
    "            fig.colorbar(p1,shrink=0.7,aspect=10,ticks=[-1,-0.5,0,0.5,1])\n",
    "        if len(mlabs)==4:\n",
    "            plt.xticks([0,1,2,3],mlabs)\n",
    "            plt.yticks([0,1,2,3],mlabs)\n",
    "        else:\n",
    "            plt.xticks(np.arange(0,n,1),mlabs)\n",
    "            plt.yticks(np.arange(0,n,1),mlabs)\n",
    "        plt.title(str(stag) + ' correlation matrix for ' + str(tlab));\n",
    "        \n",
    "        # second matrix plot (OPTIONAL)\n",
    "        if iplot_second_set==1:\n",
    "            ax2 = fig.add_subplot(nr,nc,pind+1)\n",
    "            p2 = ax2.imshow(F2,cmap='jet',vmin=-1,vmax=1); #caxis([-1 1]), colorbar\n",
    "            ax2.xaxis.tick_top()\n",
    "            if len(mlabs)==4:\n",
    "                plt.xticks([0,1,2,3],mlabs)\n",
    "                plt.yticks([0,1,2,3],mlabs)\n",
    "            else:\n",
    "                plt.xticks(np.arange(0,n,1),mlabs)\n",
    "                plt.yticks(np.arange(0,n,1),mlabs)\n",
    "            plt.title(str(stag) + ' correlation matrix for ' + str(tlab2))\n",
    "            fig.colorbar(p2,shrink=0.7,aspect=10,ticks=[-1,-0.5,0,0.5,1])\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    # scatterplots\n",
    "    if n > NMAX:\n",
    "        print('n = %i is > %i, so no scatterplots made' % (n,NMAX))\n",
    "    else:\n",
    "        fig = plt.figure(figsize=(8,8))\n",
    "        nr = n-1\n",
    "        nc = n-1\n",
    "        for ii in range(n-1):\n",
    "            jj=ii+1\n",
    "            while jj<n: \n",
    "                px = np.array([msamples[ii,:]])\n",
    "                py = np.array([msamples[jj,:]])\n",
    "                iplot = nc*(ii) + jj;\n",
    "                #disp([ii jj iplot]);\n",
    "                ax = fig.add_subplot(nr,nc,iplot)\n",
    "                ax.plot(px,py,sym1,markersize=2)\n",
    "                ax.xaxis.tick_bottom()\n",
    "                plt.xlabel(mlabs[ii])\n",
    "                plt.ylabel(mlabs[jj])\n",
    "                st1 = 'corr(%s) = %.3f (%.3f)' % (tlab,np.corrcoef(px,py,ddof=0)[0,1],rho[ii,jj])\n",
    "                if iplot_second_set==1:\n",
    "                    px = msamples2[ii,:]\n",
    "                    py = msamples2[jj,:]\n",
    "                    plt.plot(px,py,sym2,markersize=2);\n",
    "                    st2 = 'corr(%s) = %.3f (%.3f)' % (tlab2,np.corrcoef(px,py,ddof=0)[0,1],rho2[ii,jj])\n",
    "                    plt.title(str(st1)+'\\n'+ str(st2))\n",
    "                else:\n",
    "                    plt.title(str(st1))\n",
    "                jj+=1\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display distributions for each model parameter (nparm ROWS of cov_samples_m)\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "nr = 2\n",
    "nc = 2\n",
    "for kk in range(nparm):\n",
    "    sigma = sigma_post[kk]\n",
    "    edges = np.arange(-4*sigma,4*sigma,sigma/2)\n",
    "    etemp = mcov_samples[kk,:]\n",
    "    plt.subplot(nr,nc,kk+1)\n",
    "    plot_histo(etemp,edges)\n",
    "    plt.ylim([0,0.4])\n",
    "    plt.grid\n",
    "    stl1 = 'mpost samples'\n",
    "    stl2 = 'Model parameter ' + str(kk+1) + ' (' + str(mlabs[kk]) + ')'\n",
    "    stl3 = 'mean = %.5s; std = %.5s' % (np.mean(etemp),np.std(etemp))\n",
    "    if kk==0: \n",
    "        plt.title(str(stl1) + ': ' + str(stl2) + '\\n' + str(stl3))\n",
    "    else: \n",
    "        plt.title(str(stl2) + '\\n' + str(stl3))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# correlation matrices and scatterplots\n",
    "#plot_covsamples(mprior_samples,rho_prior,'mprior','b.',[],[],[],[],mlabs);\n",
    "plot_covsamples(mpost_samples,rho_post,'mpost','r.',[],[],[],[],mlabs)\n",
    "plot_covsamples(mprior_samples,rho_prior,'mprior','b.',mpost_samples,rho_post,'mpost','r.',mlabs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 'physical view' of estimated posterior data uncertainties\n",
    "# note: plot either sigma_post_d (from Cpost_d) or std_d_samples (from d(Cpost_samples))\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.subplot(aspect=1)\n",
    "plt.plot(mpost_samples[0,:],mpost_samples[1,:],'c.')\n",
    "plt.plot(mpost[0],mpost[1],'o',markersize=10,markerfacecolor='c',markeredgecolor='w')\n",
    "#plot(mprior(1),mprior(2),'o','markersize',10,'markerfacecolor','b','markeredgecolor','w');\n",
    "#scatter(xrec,yrec,16^2,sigma_post_d,'filled','V'); title('estimated uncertainties for posterior predictions'); \n",
    "plt.scatter(xrec,yrec,16**2,std_d_samples,marker='v',edgecolors='k',cmap='hot')\n",
    "plt.xlabel('X distance (km)')\n",
    "plt.ylabel('Y distance (km)')\n",
    "plt.colorbar(shrink=0.8)\n",
    "plt.title('uncertainties for posterior predictions, computed from samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot predictions for samples of the posterior\n",
    "plot_covsamples(d_samples,rho_post_d,'dpost','k.',[],[],[],[],[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plot solution\n",
    "plot_epicenters(mprior_samples,mprior,minitial,mtarget,[xrec,yrec,0,axepi],mpost)\n",
    "# plot the cpost0 samples and re-plot the two markers\n",
    "plt.plot(mpost_samples[0,:],mpost_samples[1,:],'c.')\n",
    "plt.plot(mpost[0],mpost[1],'o',markersize=10,markerfacecolor='c',markeredgecolor='w')\n",
    "plt.plot(mtarget[0],mtarget[1],'o',markersize=10,markerfacecolor='r',markeredgecolor='w')\n",
    "plt.plot(minitial[0],minitial[1],'o',markersize=10,markerfacecolor='k',markeredgecolor='w')\n",
    "\n",
    "# CODE HERE FOR PLOTTING THE ELLIPSE\n",
    "# chi2inv() in lib_peip.py\n",
    "# plot_ellipse() in lib_geos.py\n",
    "\n",
    "\n",
    "\n",
    "plt.title('samples of prior (blue) and posterior (cyan)')\n",
    "#plt.xlim([0,100])\n",
    "#plt.ylim([0,100])\n",
    "#plt.axis('equal')\n",
    "if xx==1:\n",
    "    plt.axis(axepi)\n",
    "#plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "7f341b321178e62fbbfacfb91167ed733f9907f6f0692efdceb8394e8303ba73"
  },
  "kernelspec": {
   "display_name": "base",
   "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
}