forward_epicenter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "from numpy.linalg import inv\n",
    "import numpy.matlib\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')                  # do not show warnings\n",
    "\n",
    "# Python adaptation of forward_epicenter.m written by Carl Tape\n",
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "# Python coding by Nealey Sims\n",
    "#\n",
    "# This forward problem describes the straight-line ray path travel time\n",
    "# from a source with an unknown origin time to a receiver in a\n",
    "# homogeneous medium with unknown velocity V.\n",
    "#\n",
    "# NOTE: variables such as minitial, mtarget, and eobs may be over-written\n",
    "# within optimization.m when the user specifies multiple runs.\n",
    "#\n",
    "# NOTE: If there are constants in the Matlab in-line functions, then these\n",
    "# constants are assigned at the time the function is INITIALIZED, not at\n",
    "# the time the function is called.\n",
    "#\n",
    "# CALLED BY : hw_genlsq.ipynb\n",
    "#\n",
    "# Carl Tape, 2010-03-18\n",
    "#\n",
    "\n",
    "### Must be run within hw_genlsq.ipynb in order to obtain key variables ###\n",
    "\n",
    "run_forward_epicenter=False      #default=False; run without genlsq = True\n",
    "\n",
    "### be sure to change again to avoid changing variables in multiple locations\n",
    "\n",
    "if run_forward_epicenter==True:\n",
    "    #=========================================\n",
    "    # USER INPUT\n",
    "    \n",
    "    nsamples = 1000;\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",
    "    #=========================================\n",
    "    \n",
    "    inormalization = 1;\n",
    "    stnsamples = str(nsamples) + ' samples'\n",
    "    stlabS = ('Sd(m^n)','Sm(m^n)','S(m^n) = Sd + Sm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#==========================================================================\n",
    "# FORWARD PROBLEM \n",
    "\n",
    "ndata = 12;\n",
    "nparm = 4;\n",
    "\n",
    "# labels for parameters\n",
    "#mlabs = {'xs','ys','ts','v = ln(V/V0)'};\n",
    "mlabs = ['xs','ys','ts','v']\n",
    "ulabs = ['km','km','s','none']\n",
    "\n",
    "# velocity for prior model\n",
    "# v is the unitless logarithmic velocity, v = ln(V/V0)\n",
    "# V = V0*exp(v) has units km/s\n",
    "V = 5;      # km/s\n",
    "V0 = 1;     # arbitrary scale factor (V is defined below)\n",
    "\n",
    "# travel time computation (homogeneous velocity; straight ray paths)\n",
    "\n",
    "def d2(x,y,xr,yr):\n",
    "    d2x = (xr-x)**2 + (yr-y)**2\n",
    "    return d2x\n",
    "def d1(x,y,xr,yr):\n",
    "    d1x = np.sqrt(d2(x,y,xr,yr))\n",
    "    return d1x\n",
    "def tt(x,y,ts,v,xr,yr):\n",
    "    ttx = ts + (d1(x,y,xr,yr)/(V0*math.exp(v)))\n",
    "    return ttx\n",
    "\n",
    "# receiver locations for uniform grid\n",
    "xrecmin = 10; xrecmax = 80;\n",
    "yrecmin = 20; yrecmax = 90;\n",
    "xvec = np.linspace(xrecmin,xrecmax,4);\n",
    "yvec = np.linspace(yrecmin,yrecmax,3);\n",
    "[X,Y] = np.meshgrid(xvec,yvec);\n",
    "xrec = X.flatten(order='F')\n",
    "yrec = Y.flatten(order='F')\n",
    "\n",
    "# computation of predictions\n",
    "def d(m):\n",
    "    dm=np.array([tt(m[0],m[1],m[2],m[3],xrec[0],yrec[0]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[1],yrec[1]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[2],yrec[2]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[3],yrec[3]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[4],yrec[4]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[5],yrec[5]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[6],yrec[6]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[7],yrec[7]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[8],yrec[8]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[9],yrec[9]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[10],yrec[10]),\n",
    "              tt(m[0],m[1],m[2],m[3],xrec[11],yrec[11])],dtype='float')\n",
    "    return dm\n",
    "\n",
    "# ndata x nparm matrix of partial derivatives (differentiate tt with respect to each parameter)\n",
    "# evaluated at model m = (m(1),m(2),m(3),m(4))\n",
    "\n",
    "def G(m):\n",
    "    Gm=np.array([[-(d2(m[0],m[1],xrec[0],yrec[0]))**(-1/2)@(xrec[0]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[0],yrec[0]))**(-1/2)@(yrec[0]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[0],yrec[0])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[1],yrec[1]))**(-1/2)@(xrec[1]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[1],yrec[1]))**(-1/2)@(yrec[1]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[1],yrec[1])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[2],yrec[2]))**(-1/2)@(xrec[2]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[2],yrec[2]))**(-1/2)@(yrec[2]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[2],yrec[2])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[3],yrec[3]))**(-1/2)@(xrec[3]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[3],yrec[3]))**(-1/2)@(yrec[3]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[3],yrec[3])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[4],yrec[4]))**(-1/2)@(xrec[4]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[4],yrec[4]))**(-1/2)@(yrec[4]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[4],yrec[4])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[5],yrec[5]))**(-1/2)@(xrec[5]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[5],yrec[5]))**(-1/2)@(yrec[5]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[5],yrec[5])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[6],yrec[6]))**(-1/2)@(xrec[6]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[6],yrec[6]))**(-1/2)@(yrec[6]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[6],yrec[6])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[7],yrec[7]))**(-1/2)@(xrec[7]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[7],yrec[7]))**(-1/2)@(yrec[7]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[7],yrec[7])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[8],yrec[8]))**(-1/2)@(xrec[8]-m[0])/(V0*math.exp(m[3]))  ,   -(d2(m[0],m[1],xrec[8],yrec[8]))**(-1/2)@(yrec[8]-m[1])/(V0*math.exp(m[3]))  ,   1 , -d1(m[0],m[1],xrec[8],yrec[8])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[9],yrec[9]))**(-1/2)@(xrec[9]-m[0])/(V0*math.exp(m[3]))  ,  -(d2(m[0],m[1],xrec[9],yrec[9]))**(-1/2)@(yrec[9]-m[1])/(V0*math.exp(m[3]))  ,  1 , -d1(m[0],m[1],xrec[9],yrec[9])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[10],yrec[10]))**(-1/2)@(xrec[10]-m[0])/(V0*math.exp(m[3]))  ,  -(d2(m[0],m[1],xrec[10],yrec[10]))**(-1/2)@(yrec[10]-m[1])/(V0*math.exp(m[3]))  ,  1 , -d1(m[0],m[1],xrec[10],yrec[10])/(V0*math.exp(m[3]))],\n",
    "                 [-(d2(m[0],m[1],xrec[11],yrec[11]))**(-1/2)@(xrec[11]-m[0])/(V0*math.exp(m[3]))  ,  -(d2(m[0],m[1],xrec[11],yrec[11]))**(-1/2)@(yrec[11]-m[1])/(V0*math.exp(m[3]))  ,  1 , -d1(m[0],m[1],xrec[11],yrec[11])/(V0*math.exp(m[3]))]\n",
    "    ],dtype='float')\n",
    "    return Gm\n",
    "\n",
    "# nparm x nparm matrix of second partial derivatives (only used in full Newton method)\n",
    "# note: this contains the measurement index ii\n",
    "\n",
    "def G2(m,ii):\n",
    "    G2m=np.array([\n",
    "    [(d1(m[0],m[1],xrec[ii],yrec[ii]))**-3@(yrec[ii]-m[1])**2/(V0*math.exp(m[3])),               -(d1(m[0],m[1],xrec[ii],yrec[ii]))**-3@(xrec[ii]-m[0])*(yrec[ii]-m[1])/(V0*math.exp(m[3])),    0, (d1(m[0],m[1],xrec[ii],yrec[ii]))**-1@(xrec[ii]-m[0])/(V0*math.exp(m[3]))],\n",
    "    [-(d1(m[0],m[1],xrec[ii],yrec[ii]))**-3@(xrec[ii]-m[0])*(yrec[ii]-m[1])/(V0*math.exp(m[3])),  (d1(m[0],m[1],xrec[ii],yrec[ii]))**-3@(xrec[ii]-m[0])**2/(V0*math.exp(m[3])),                  0, (d1(m[0],m[1],xrec[ii],yrec[ii]))**-1@(yrec[ii]-m[1])/(V0*math.exp(m[3]))],\n",
    "    [0,                                                                                            0,                                                                                            0, 0],\n",
    "    [(d1(m[0],m[1],xrec[ii],yrec[ii]))**-1@(xrec[ii]-m[0])/(V0*math.exp(m[3])),                  (d1(m[0],m[1],xrec[ii],yrec[ii]))**-1@(yrec[ii]-m[1])/(V0*math.exp(m[3])),                      0, d1(m[0],m[1],xrec[ii],yrec[ii])/(V0*math.exp(m[3]))]])\n",
    "    return G2m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RANDOM VECTORS (for sampling covariance matrices)\n",
    "\n",
    "# Gaussian random vectors, each with mean = 0 and standard deviation = 1\n",
    "randn_vecs_m = np.random.randn(nparm,nsamples);   # model\n",
    "randn_vecs_d = np.random.randn(ndata,nsamples);   # data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PRIOR MODEL (MEAN MODEL) : xs, ys, ts, v\n",
    "\n",
    "# prior model\n",
    "# note: V and V0 are defined above\n",
    "mprior = np.transpose(np.array([[35, 45, 16, np.log(V/V0)]]))\n",
    "\n",
    "# prior model covariance matrix (assumed to be diagonal)\n",
    "sigma_prior = np.array([10, 10, 0.5, 0.2]).T       # standard deviations\n",
    "cprior0     = np.diag( np.square(sigma_prior) );   # diagonal covariance matrix\n",
    "if inormalization==1:\n",
    "    Cmfac = nparm;\n",
    "else:\n",
    "    Cmfac = 1;\n",
    "\n",
    "cprior   = np.dot(Cmfac, cprior0)             # WITH NORMALIZATION FACTOR\n",
    "icprior  = inv(cprior);                 # WITH NORMALIZATION FACTOR\n",
    "icprior0 = inv(cprior0);\n",
    "Lprior   = np.linalg.cholesky(cprior0).transpose()      # square-root (lower triangular)\n",
    "# sample the prior model distribution using the square-root UNNORMALIZED covariance matrix\n",
    "cov_samples_m  = np.dot(Lprior, randn_vecs_m)\n",
    "mprior_samples = np.matlib.repmat(mprior,1,nsamples) + cov_samples_m\n",
    "\n",
    "# compute the norm of each model sample using the inverse NORMALIZED covariance matrix\n",
    "norm2_mprior = np.zeros((nsamples,1))\n",
    "for xx in range(nsamples):\n",
    "    dm = np.array([mprior_samples[:,xx]]).T - mprior\n",
    "    norm2_mprior[xx] = np.dot(dm.T, icprior)@ dm\n",
    "#figure; plot(norm2_mprior,'.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INITIAL MODEL\n",
    "\n",
    "# minitial is DIFFERENT FOR EACH RUN, or you can fix it for testing purposes\n",
    "if irandom_initial_model == 1:\n",
    "    minitial = np.array([mprior_samples[:,0]]).T     # first sample (random)\n",
    "else:\n",
    "    minitial = np.array([                        # fixed\n",
    "        [46.5236],\n",
    "        [40.1182],\n",
    "        [15.3890],\n",
    "        [1.7748]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TARGET MODEL\n",
    "\n",
    "# mtarget is DIFFERENT FOR EACH RUN, or you can fix it for testing purposes\n",
    "if irandom_target_model == 1 : \n",
    "    mtarget = np.array([mprior_samples[:,-1]]).T    # last sample (random)\n",
    "else:\n",
    "    mtarget = np.array([                         # fixed  \n",
    "           [21.2922],\n",
    "           [46.2974],\n",
    "           [16.1314],\n",
    "            [2.0903]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TARGET DATA\n",
    "\n",
    "#dtarget = np.transpose(np.array(d(mtarget)))\n",
    "dtarget = np.array(d(mtarget))\n",
    "# data covariance matrix (assumed to be diagonal)\n",
    "tsigma = 0.5;                       # uncertainty in arrival time measurement, seconds\n",
    "sigma_obs = tsigma * np.ones((ndata)) # standard deviations\n",
    "cobs0     = np.diag( np.square(sigma_obs) )   # diagonal covariance matrix\n",
    "\n",
    "if inormalization==1:\n",
    "    Cdfac = ndata;\n",
    "else:\n",
    "    Cdfac = 1;\n",
    "\n",
    "cobs   = Cdfac* cobs0             # WITH NORMALIZATION FACTOR\n",
    "icobs  = inv(cobs);                 # WITH NORMALIZATION FACTOR\n",
    "icobs0 = inv(cobs0);\n",
    "Lcobs  = np.linalg.cholesky(cobs0).transpose()      # square-root (lower triangular)\n",
    "\n",
    "# sample the data distribution using the square-root UNNORMALIZED covariance matrix\n",
    "cov_samples_d = np.dot(Lcobs, randn_vecs_d)\n",
    "dobs_samples  = np.matlib.repmat(dtarget,1,nsamples) + cov_samples_d;\n",
    "\n",
    "# compute the norm of each data sample using the inverse NORMALIZED covariance matrix\n",
    "norm2_dobs = np.zeros((nsamples,1))\n",
    "for xx in range(nsamples):\n",
    "    dd = np.array([dobs_samples[:,xx]]).T - dtarget\n",
    "    norm2_dobs[xx] = np.dot(dd.T, icobs)@ dd;\n",
    "\n",
    "# Pick the uncertainties for the target data by simply choosing a sample\n",
    "#   from the realizations of the data covariance.\n",
    "# NOTE: eobs is DIFFERENT FOR EACH RUN, or you can fix it for testing purposes\n",
    "\n",
    "if idata_errors == 0:\n",
    "    eobs = np.zeros((ndata,1)).T          # no errors\n",
    "elif idata_errors == 1:\n",
    "    eobs = np.array([cov_samples_d[:,1]]);      # first sample (random)\n",
    "elif idata_errors == 2:\n",
    "    eobs = np.array([[                        # fixed\n",
    "            -0.8689,\n",
    "            -0.4666,\n",
    "            -0.0516,\n",
    "            0.2411,\n",
    "            0.2825,\n",
    "            -0.2301,\n",
    "            0.1977,\n",
    "            0.3291,\n",
    "            1.0063,\n",
    "            0.5674,\n",
    "            0.1348,\n",
    "            0.4603]])\n",
    "\n",
    "# \"true\" observations (includes added errors)\n",
    "\n",
    "dobs = np.add(dtarget, np.transpose(eobs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_epicenters(mprior_samples,mprior,minitial,mtarget,opts,mpost=[]):\n",
    "    #PLOT_EPICENTERS plots epicenters within the source-station geometry\n",
    "    #\n",
    "    # Carl Tape, 2010-02-24\n",
    "    \n",
    "    # input options\n",
    "    xrec = opts[0];\n",
    "    yrec = opts[1];\n",
    "    iray = opts[2];\n",
    "    ax0  = opts[3];\n",
    "    \n",
    "    nparm = np.shape(mprior_samples)[0]; nsample=np.shape(mprior_samples)[1]\n",
    "    ndata = len(xrec);\n",
    "    \n",
    "    # indices of xs and ys within model vector\n",
    "    ix = 0;\n",
    "    iy = 1;\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,8))\n",
    "    msizer = 20;    # receiver size\n",
    "    msizes = 15;    # source size\n",
    "    rthick = 1;     # receiver edge thickness\n",
    "    sthick = 2;     # source edge thickness\n",
    "    rfsize = 10;\n",
    "    \n",
    "    #xmin = ax0(1);\n",
    "    #xmax = ax0(2);\n",
    "    #ymin = ax0(3);\n",
    "    #ymax = ax0(4);\n",
    "    \n",
    "    # plot ray paths\n",
    "    if iray==1:\n",
    "        for ii in range(ndata):\n",
    "            plt.plot([minitial[ix], xrec[ii]],[minitial[iy], yrec[ii]],'k');\n",
    "            \n",
    "    # plot option depends on if a posterior model is passed\n",
    "    if len(mpost)>0:\n",
    "        p0 = plt.plot(mprior_samples[ix,:],mprior_samples[iy,:],'b.');\n",
    "        p1 = plt.plot(minitial[ix],minitial[iy],'o');\n",
    "        p2 = plt.plot(mpost[ix],mpost[iy],'co',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        pP = plt.plot(mprior[ix],mprior[iy],'bo',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        pT = plt.plot(mtarget[ix],mtarget[iy],'ro',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        plt.plot(xrec,yrec,'wv',markersize=msizer,markeredgecolor='k',linewidth=rthick)\n",
    "        for ii in range(ndata):\n",
    "            plt.text(xrec[ii],yrec[ii], str(ii+1),fontsize=rfsize,color='r',\n",
    "                     horizontalalignment='center',verticalalignment='center')\n",
    "\n",
    "        plt.legend([p0[0], p1[0], p2[0], pP[0], pT[0]],['Cprior sample','minitial','mpost','mprior','mtarget'],loc='upper right');\n",
    "        \n",
    "    else:\n",
    "        if len(mprior_samples)>1:\n",
    "            p0 = plt.plot(mprior_samples[ix,:],mprior_samples[iy,:],'b.');\n",
    "        \n",
    "        p1 = plt.plot(minitial[ix],minitial[iy],'ko',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        pP = plt.plot(mprior[ix],mprior[iy],'bo',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        pT = plt.plot(mtarget[ix],mtarget[iy],'ro',markersize=msizes,markeredgecolor='w',markeredgewidth='3',linewidth=sthick)\n",
    "        plt.plot(xrec,yrec,'wv',markersize=msizer,markeredgecolor='k',linewidth=rthick)\n",
    "        for ii in range(ndata):\n",
    "            plt.text(xrec[ii],yrec[ii],str(ii+1),fontsize=rfsize,color='r',\n",
    "                     horizontalalignment='center',verticalalignment='center')\n",
    "        if len(mprior_samples)>1:\n",
    "            plt.legend([p0[0], p1[0], pP[0], pT[0]],['Cprior sample','minitial','mprior','mtarget'],loc='upper right');\n",
    "            x = mprior_samples[ix,:]\n",
    "            y = mprior_samples[iy,:]\n",
    "            #iin = find( and( and( x >= ax0(1), x <= ax0(2)), and(y >= ax0(3), y <= ax0(4))) );\n",
    "            iin=0\n",
    "            for jj in range(len(x)):\n",
    "                if x[jj] >= ax0[0] and x[jj] <=ax0[1] and y[jj] >= ax0[2] and y[jj]<=ax0[3]:\n",
    "                    iin+=1\n",
    "            plt.title('%i/%i (%.2f) prior samples inside the region shown' %\n",
    "                (iin,nsample,iin/nsample))\n",
    "        else:\n",
    "            plt.legend([p1[0], pP[0], pT[0]],['minitial','mprior','mtarget'],loc='upper right');\n",
    "    \n",
    "    #set(gca,'xtick',[0:20:100],'ytick',[0:20:100]);\n",
    "    plt.gca().set_aspect('equal', adjustable='box')\n",
    "    plt.xlabel('X distance (km)'); plt.ylabel('Y distance (km)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# PLOTS\n",
    "\n",
    "axepi = [0, 100, 0, 100];\n",
    "\n",
    "# source-receiver geometry with ray paths\n",
    "plot_epicenters([[]],mprior,minitial,mtarget,[xrec,yrec,1,axepi]);\n",
    "\n",
    "# source-receiver geometry with ray paths\n",
    "plot_epicenters(mprior_samples,mprior,minitial,mtarget,[xrec,yrec,1,axepi]);\n",
    "\n",
    "# with prior samples (no ray paths)\n",
    "opts = [xrec,yrec,0,axepi];\n",
    "plot_epicenters(mprior_samples,mprior,minitial,mtarget,opts);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "seismo",
   "language": "python",
   "name": "seismo"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}