import numpy as np
import matplotlib.pyplot as plt
import emcee
import corner
import math
import utm
from pysurf96 import surf96
%matplotlib inline
plt.rcParams['figure.figsize'] = (20,10)
import os
import sys
import time
import multiprocessing as mp
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import obspy
from scipy.signal import iirpeak, tukey, butter, filtfilt, hilbert
from func_disp import *
import glob
f1, f2 = 0.25, 4.0
def extr_disp(input):
'''
# Extract dispersion curve from a frequency-velocity image
# Input:
# input - a list of the parameters below:
# f - frequency
# v - phase velocity
# disp - amplitude at each (f,v) point (nv*nf)
'''
f, v, disp = input
nonzero_inds=np.sum(disp,axis=0)!=0
if np.sum(nonzero_inds)==0:
return np.repeat(np.nan,len(f))
zero_inds=np.sum(disp,axis=0)==0
fmin,fmax=np.min(f[nonzero_inds]),np.max(f[nonzero_inds])
f_ref_all = [2.0,1.0,0.5,0.4,0.26] # give the starting point of dispersion picking, a list
vmax_all = [600,1000,2500,3400,3450]# corresponding to f_ref_set, max velocity at that frequency
f_ref_set = []
vmax_set = []
for i in range(len(f_ref_all)):
if (f_ref_all[i]<fmax)&(f_ref_all[i]>fmin):
f_ref_set.append(f_ref_all[i])
vmax_set.append(vmax_all[i])
if len(f_ref_set)>0:
xpt = []
ypt = []
for k in range(len(f_ref_set)):
f_ref = f_ref_set[k]
vmax = vmax_set[k]
index = np.where(f >= f_ref)[0][0]
f_ref = f[index]
disp_ref = disp[:, index]
v_ref = v[np.argmax(disp_ref[v<vmax])]
ypt.append(np.abs(v - v_ref).argmin())
xpt.append(np.abs(f - f_ref).argmin())
ArrPt = AutoSearch(ypt[0], xpt[0], disp, step=5)
#ArrPt = AutoSearchMultiplePoints(np.array(ypt), np.array(xpt), disp)
voutput = v[ArrPt]
voutput[zero_inds]=np.nan
voutput[voutput==v[-1]]=np.nan
voutput[voutput==v[0]]=np.nan
return voutput
else:
return np.repeat(np.nan,len(f))
from scipy import interpolate
filenames=glob.glob('data/dcurves/Shot*')
filenames.sort()
x=np.load(filenames[0])
vels=x['vels']
frqs=x['freq']
vel_std_lb=interpolate.interp1d(vels[:-1]/1000,2*np.diff(vels)/1000)
if True:
import glob
import tqdm
picks=list()
for i in range(1574):
picks.append(list())
data_stacks=list()
for i in range(1574):
data_stacks.append(np.zeros((len(vels),len(frqs))))
for i_file in tqdm.tqdm(range(len(filenames))):
fname='data/dcurves/Shot_%s.npz'%i_file
x=np.load(fname)
data=x['res']
if (data!=data).any():
print(fname)
recs=x['rec']
for i in range(len(data)):
picks[recs[i]].append(extr_disp(list((frqs,vels,data[i]))))
data_stacks[recs[i]]+=data[i]
curve_fillmed=np.zeros((len(picks),len(frqs)))
curve_fillmean=np.zeros((len(picks),len(frqs)))
curve_fillstd=np.zeros((len(picks),len(frqs)))
for i in tqdm.tqdm(range(len(picks))):
if len(picks[i])>0:
curve_fillmean[i]=np.nanmean(np.array(picks[i]),axis=0)/1000
curve_fillstd[i]=np.nanstd(np.array(picks[i]),axis=0)/1000
rep_inds=np.where(vel_std_lb(curve_fillmean)>curve_fillstd)
for repi in tqdm.tqdm(range(len(rep_inds[0]))):
curve_fillstd[rep_inds[0][repi],rep_inds[1][repi]]=vel_std_lb(curve_fillmean[rep_inds[0][repi],rep_inds[1][repi]])
# np.save('data/curve_fillmean.npy',curve_fillmean)
# np.save('data/curve_fillstd.npy',curve_fillstd)
#curve_fillmean=np.load('data/curve_fillmean.npy')
#curve_fillstd=np.load('data/curve_fillstd.npy')
100%|███████████████████████████████████████| 1574/1574 [12:25<00:00, 2.11it/s] 100%|██████████████████████████████████████| 1574/1574 [00:02<00:00, 672.19it/s] 100%|████████████████████████████████████| 2105/2105 [00:00<00:00, 68753.73it/s]
300*2+5999
6599
#Plotting Example Dispersion Curve
fig,ax=plt.subplots(1,figsize=(20,20))
ich=300
dstack=(data_stacks[ich]/np.max(data_stacks[ich],axis=0))
ax.pcolormesh(1/frqs,vels,dstack,cmap='Spectral_r')
ax.plot(1/frqs,curve_fillmean[ich]*1e3,'k',linewidth=8)
ax.fill_between(1/frqs,curve_fillmean[ich]*1e3-curve_fillstd[ich]*1e3,curve_fillmean[ich]*1e3+curve_fillstd[ich]*1e3,alpha=0.2,color='black')
ax.set_xlim((0.4,4))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_yticks((500,1000,2000))
ax.set_yticklabels((500,1000,2000))
ax.set_xticks((0.4,0.6,1.0,2.0,3.0,4.0))
ax.set_xticklabels((0.4,0.6,1,2,3,4))
ax.set_xlabel('Period (s)')
ax.set_ylabel('Rayleigh Phase Velocity (m/s)')
plt.savefig('figs/egdcurve.png')
nan_arr = np.empty((25,1000,3))
nan_vec=np.empty((1000))
nan_arr[:] = np.nan
nan_vec[:] = np.nan
meanpick=np.concatenate((nan_vec,np.load('data/meanpick.npy')))
stdpick=np.concatenate((nan_vec,np.load('data/stdpick.npy')))
ppick=np.concatenate((nan_vec,np.load('data/ppick.npy')))
ppickmin=np.nanpercentile(ppick,5)
ppickmax=np.nanpercentile(ppick,95)
ppick[ppick<ppickmin]=ppickmin
ppick[ppick>ppickmax]=ppickmax
meanratpick=np.concatenate((nan_vec,np.load('data/meanratpick.npy')))
stdratpick=np.concatenate((nan_vec,np.load('data/stdratpick.npy')))
stdratpick[stdratpick<0.05]=0.05
stdpick[stdpick<0.035]=0.035
varpick=stdpick**2
varratpick=stdratpick**2
rat,rat_err,p_overall=np.load('data/ratdat.npy')
# function to empirically convert vs to vp (see Brocher, 2005)
def vs2vp(x):
return 0.9409+2.0947*x-0.8206*x**2+0.2683*x**3-0.0251*x**4
pi=np.pi
#function to predict rayleigh dispersion curve for our specific model setup
def pred_dcurve(theta,sb=False,zb=0.4):
# inputs:
# theta - vs model parameterization in form (vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db)
# #vs1-vs4: s velocity in layers 1-4
# #d1-d3: thickness of layers 1-3
# #vsb: s velocity in slumped block
# #db: depth of slumped block beneath the surface
# sb - assume slumped block in forward model? (boolean)
# zb - block thickness (km)
f1, f2 = 0.25, 4.0
periods= 1/np.logspace(np.log10(f1),np.log10(f2),50)[::-1]
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db=theta
depth_cum=np.cumsum((d1,d2,d3,5))
if (sb==True):
#find layer block starts in
startblock=np.where((depth_cum>db))[0][0]
#find layer block ends in
endblock=np.where((depth_cum>(db+zb)))[0][0]
betas_t=np.array((vs1,vs2,vs3,vs4))
betas=np.concatenate((betas_t[:(startblock+1)],[vsb],betas_t[endblock:]))
#empirically go from vs to vp and rho
alphas=vs2vp(betas)
rhos = (alphas**0.25)*1.74
#get layer thicknesses with block imposed over a portion of the section
depths=np.diff(np.concatenate(([0],depth_cum[:startblock],[db,db+zb],depth_cum[endblock:])))
else:
betas=np.array((vs1,vs2,vs3,vs4))
#empirically go from vs to vp and rho
alphas=vs2vp(betas)
rhos = (alphas**0.25)*1.74
depths=np.array([d1,d2,d3,5.0])
#forward model of dispersion curve
try:
disp_curve = surf96(
depths,
alphas,
betas,
rhos,
periods,
wave="rayleigh",
mode=1,
velocity="phase",
flat_earth=True,
ddc0=0.005,
h0=0.005)
return disp_curve[::-1]
except:
return np.repeat(np.nan,50)
def calc_trans_coeff(ps,alpha1,alpha2,beta1,beta2,rho1,rho2,ends=(0,0)):
# calculate transmission coefficient across layer boundary
# formulas taken from Aki and Richards
# inputs:
## ps - horizontal slowness
## alpha1 - p velocity in starting medium
## alpha2 - p velocity in ending medium
## beta1 - s velocity in starting medium
## beta2 - s velocity in ending medium
## rho1 - density in starting medium
## rho2 - density in ending medium
## ends - (starting phase: 0 or 1, ending phase: 0 or 1) 0 is p, 1 is s
### e.g. if we are considering a p to s transmission coefficient, we would take ends = (0,1)
# outputs:
# tcoeff - transmission coefficient
start=ends[0]
end=ends[1]
cos=np.cos
i1=np.arcsin(alpha1*ps)
i2=np.arcsin(alpha2*ps)
j1=np.arcsin(beta1*ps)
j2=np.arcsin(beta2*ps)
a=rho2*(1-2*(beta2*ps)**2)-rho1*(1-2*(beta1*ps)**2)
b=rho2*(1-2*(beta2*ps)**2)+2*rho1*(beta1*ps)**2
c=rho1*(1-2*(beta1*ps)**2)+2*rho2*(beta2*ps)**2
d=2*(rho2*beta2**2-rho1*beta1**2)
E=b*cos(i1)/alpha1+c*cos(i2)/alpha2
F=b*cos(j1)/beta1+c*cos(j2)/beta2
G=a-d*cos(i1)*cos(j2)/(alpha1*beta2)
H=a-d*cos(i2)*cos(j1)/(alpha2*beta1)
D=E*F+G*H*ps**2
if start==0:
if end==0:
tcoeff=(2*rho1*cos(i1)/alpha1)*F*alpha1/(alpha2*D)
if end==1:
tcoeff=(2*rho1*cos(i1)/alpha1)*H*ps*alpha1/(beta2*D)
if start==1:
if end==0:
tcoeff=(-2*rho1*cos(j1)/beta1)*G*ps*beta1/(alpha2*D)
if end==1:
tcoeff=(2*rho1*cos(j1)/beta1)*F*beta1/(beta2*D)
return tcoeff
def dir_sens(p,alpha,beta,pors):
# calculate coefficient to calculate along-fiber motion
# inputs:
## p - horizontal slowness
## alpha - p velocity in top layer
## beta - s velocity in top layer
## pors - P or S arrival? P is 0, S is 1
# output: coefficient to calculate along-fiber motion
if pors==0:
return p*alpha
if pors==1:
eta=np.sqrt(1/beta-p)
return eta*beta
def predamprat(ps,alphas,betas,rhos,direct=(0,0,0),conv=(0,1,1)):
# predict amplitude ratio of direct and converted phases traveling through three layers
# inputs:
## ps - horizontal slowness
## alphas - vector of length n, p velocity in top n layers, top to bottom
## betas - vector of length n, s velocity in top n layers, top to bottom
## rhos - vector of length n, density in top n layers
## direct - phase (p or s) the direct arrival takes on in each of the n layers,
### bottom to top, p is 0 and s is 1, e.g. a direct p arrival is (0,0,0) for three layers
## converted - phase (p or s) the converted arrival takes on in each of the n layers,
### bottom to top, p is 0 and s is 1
### e.g. a Ps arrival converted below the second layer is (0,1,1) for three layers
# output: amplitude ratio between direct and converted phases
n=len(alphas)
d_amp=1
c_amp=1
for layer_i in range(n-1):
d_amp=d_amp*calc_trans_coeff(ps,alphas[layer_i+1],alphas[layer_i],betas[layer_i+1],
betas[layer_i],rhos[layer_i+1],rhos[layer_i],
ends=direct[(n-2-layer_i):(n-layer_i)])
c_amp=c_amp*calc_trans_coeff(ps,alphas[layer_i+1],alphas[layer_i],betas[layer_i+1],
betas[layer_i],rhos[layer_i+1],rhos[layer_i],
ends=conv[(n-2-layer_i):(n-layer_i)])
dir_sens_direct=dir_sens(ps,alphas[0],betas[0],direct[-1])
dir_sens_conv=dir_sens(ps,alphas[0],betas[0],conv[-1])
return (dir_sens_direct*d_amp)/(c_amp*dir_sens_conv)
#see if model fits data
f1, f2 = 0.25, 4.0
fs=np.logspace(np.log10(f1),np.log10(f2),50)
def lnlike(theta, y, yerr,delaydata,zb):
if delaydata[0]==delaydata[0]:
ymod=pred_dcurve(theta,sb=True,zb=zb)
else:
ymod=pred_dcurve(theta,sb=False,zb=zb)
if np.all(ymod==ymod):
LnLike = -np.sum(((y[fs<2.5][::4]-ymod[fs<2.5][::4])/yerr[fs<2.5][::4])**2)/2
return LnLike
else:
return -np.inf
def delaytime(theta,ps):
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta
depth_cum=np.cumsum((d1,d2,d3,5))
startblock=np.where((depth_cum>db))[0][0]
depths=np.diff(np.concatenate(([0],depth_cum[:startblock],[db])))
betas=np.array((vs1,vs2,vs3,vs4))
alphas=vs2vp(betas)
i_s=np.arcsin(alphas*ps)
j_s=np.arcsin(betas*ps)
TTs=0
for i_layer,depth in enumerate(depths):
TTs+=(depth*np.cos(j_s[i_layer]))/betas[i_layer]-(depth*np.cos(i_s[i_layer]))/alphas[i_layer]
return TTs
#check if TT difference is consistent with observation
def lnTTprior(theta,delaydata):
if delaydata[0]==delaydata[0]:
delays_obs=delaydata[0]
delays_obs_err=delaydata[1]
ps=delaydata[2]
TTs=delaytime(theta,ps)
lnTT=-(((TTs-delays_obs)/delays_obs_err)**2)/2
else:
lnTT=0
return lnTT
#check if values are reasonable (e.g., reasonable range of velocities, depth of block isn't
#outside the model), returns -np.inf if not
def lnphysvalsprior(theta,delaydata,zb):
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta
#get maximum p velocity in block and overlying layers, we require that this is less than 1/hoz. slowness
depth_cum=np.cumsum((d1,d2,d3,5))
#require that block isn't outside the model
if (delaydata[0]==delaydata[0]):
ps=delaydata[2]
if (db+zb)>depth_cum[-1]:
return -np.inf
else:
startblock=np.where((depth_cum>db))[0][0]
vpmax=vs2vp(np.max(np.concatenate((np.array((vs1,vs2,vs3,vs4))[:startblock+1],[vsb]))))
if (vpmax>1/ps)|(vpmax>1/p_overall):
return -np.inf
#require reasonable velocities
if 0.4<vs1<0.8 and 0.4<vs2<1.8 and 0.8<vs3<2.5 and 1.5<vs4<3.5 and 1.2<vsb<3.0 and 0.02<=d1<1.0 and 0.02<=d2<2.0 and 0.02<=d3<3.0 and 2>db>0 and vs1<vs2<vs3<vs4:
return 0.0
else:
return -np.inf
def lnamp(theta,delaydata):
if delaydata[0]==delaydata[0]:
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta
depth_cum=np.cumsum((d1,d2,d3,5))
#find layer block starts in
startblock=np.where((depth_cum>db))[0][0]
betas_t=np.array((vs1,vs2,vs3,vs4))
betas=np.concatenate((betas_t[:(startblock+1)],[vsb]))
n=len(betas)
alphas=vs2vp(betas)
rhos=(alphas**0.25)*1.74
direct=np.zeros(n)
conv=np.ones(n)
conv[0]=0
amprat_pred=predamprat(p_overall,alphas,betas,rhos,direct=direct,conv=conv)
return -(((amprat_pred-rat)/rat_err)**2)/2
else:
return 0
def lnprob_step1(theta, dispcurve, dispcurve_err,delaydata,zb):
theta_alt=theta.copy()
theta_alt[4:7][theta_alt[4:7]<0.1]=0.1
if theta_alt[4]>1:
theta_alt[4]=1
if theta_alt[5]>2:
theta_alt[5]=2
if theta_alt[6]>3:
theta_alt[6]=3
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta_alt
vp1,vp2,vp3,vp4,vpb=vs2vp(np.array((vs1,vs2,vs3,vs4,vsb)))
lp = lnphysvalsprior(theta_alt,delaydata,zb)
if np.isfinite(lp):
ldc = lnlike(theta_alt, dispcurve, dispcurve_err,delaydata,zb)
lTT = lnTTprior(theta_alt,delaydata)
lnA=lnamp(theta_alt,delaydata)
return lp + ldc + lTT + lnA
else:
return -np.inf
def lnprob_step2(theta, dispcurve, dispcurve_err,delaydata,zb):
theta_alt=theta.copy()
theta_alt[4:7][theta_alt[4:7]<0.1]=0.1
if theta_alt[4]>1:
theta_alt[4]=1
if theta_alt[5]>2:
theta_alt[5]=2
if theta_alt[6]>3:
theta_alt[6]=3
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta_alt
vp1,vp2,vp3,vp4,vpb=vs2vp(np.array((vs1,vs2,vs3,vs4,vsb)))
lp = lnphysvalsprior(theta_alt,delaydata,zb)
if np.isfinite(lp):
ldc = lnlike(theta_alt, dispcurve, dispcurve_err,delaydata,zb)
lTT = lnTTprior(theta_alt,delaydata)
lcv=constvfun(vs1,vs2,vs3,vs4,vsb)[0]
lnA=lnamp(theta_alt,delaydata)
if lcv==lcv:
return lp + ldc + lTT + lcv + lnA
else:
return -np.inf
else:
return -np.inf
import scipy.stats
from geneticalgorithm import geneticalgorithm as ga
import sys, os
def blockPrint():
sys.stdout = open(os.devnull, 'w')
def enablePrint():
sys.stdout = sys.__stdout__
def findpeak(thetasall_point,sbbool):
if sbbool:
k_t=scipy.stats.gaussian_kde(thetasall_point.T)
def k(x):
return -k_t(x)
lbs=np.array((0.4,0.4,0.8,1.5,0.1,0.1,0.1,1.2,0))
ubs=np.array((0.8,1.8,2.5,3.5,1,2,3,3,2))
varbound=np.concatenate(([lbs],[ubs])).T
model=ga(function=k,dimension=9,variable_type='real',variable_boundaries=varbound,convergence_curve =False, progress_bar= False)
blockPrint()
model.run()
enablePrint()
else:
k_t=scipy.stats.gaussian_kde(thetasall_point[:,:7].T)
def k(x):
return -k_t(x)
lbs=np.array((0.4,0.4,0.8,1.5,0.1,0.1,0.1))
ubs=np.array((0.8,1.8,2.5,3.5,1,2,3))
varbound=np.concatenate(([lbs],[ubs])).T
model=ga(function=k,dimension=7,variable_type='real',variable_boundaries=varbound,convergence_curve =False, progress_bar= False)
blockPrint()
model.run()
enablePrint()
return model.output_dict['variable']
import tqdm
for zb in np.linspace(0.1,0.6,6):
if True:
thetasall=np.zeros((60,500,9))
probsall=np.zeros((60,500))
thetabests=np.zeros((60,9))
thetapeaks=np.empty((60,9))
thetapeaks[:]=np.nan
sbbools=np.zeros(60)
dispcurves=np.zeros((60,50))
dispcurve_errs=np.zeros((60,50))
delaydata_all=np.zeros((60,3))
for i in tqdm.tqdm(np.arange(0,60)):
startind=25*i
midind=(25*i)+50
endind=(25*i)+100
if (0<=(startind*2-25))&(((endind*2)-25)<len(meanpick)):
delaydata=np.zeros(3)
delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)])
delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)]))
delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)])
else:
delaydata=np.repeat(np.nan,3)
delaydata_all[i]=delaydata
sbbool=False
if delaydata[0]==delaydata[0]:
sbbool=True
sbbools[i]=sbbool
cfm_use=curve_fillmean[startind:endind]
cfs_use=curve_fillstd[startind:endind]
dispcurve=np.average(cfm_use, axis=0, weights=1/cfs_use**2)
dispcurve_err=np.sqrt(np.average(cfs_use**2,axis=0,weights=1/cfs_use**2))
dispcurves[i]=dispcurve
dispcurve_errs[i]=dispcurve_err
data = (dispcurve, dispcurve_err, delaydata, zb)
nwalkers=128
niter=500
initial=np.array((0.6, 0.7, 1.4, 2.9, 0.2, 0.2,0.2,2.4,0.6))
ndim = len(initial)
p0 = [np.array(initial) + 1e-1 * np.random.randn(ndim) for k in range(nwalkers)]
def main(p0,nwalkers,niter,ndim,lnprob,data):
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=data)
p0, _, _ = sampler.run_mcmc(p0, 100)
sampler.reset()
pos, prob, state = sampler.run_mcmc(p0, niter)
return sampler, pos, prob, state
sampler, pos, prob, state = main(p0,nwalkers,niter,ndim,lnprob_step1,data)
samples = sampler.flatchain
size=500
thetas=samples[np.random.randint(len(samples), size=size)]
thetasall[i,:,:]=thetas
probs=np.zeros(size)
for j in range(size):
if lnphysvalsprior(thetas[j], delaydata,zb)==0:
probs[j]=np.exp(lnprob_step1(thetas[j], dispcurve, dispcurve_err, delaydata,zb))
else:
probs[j]=0
probsall[i,:]=probs
thetabest=samples[np.argmax(sampler.flatlnprobability)]
thetabests[i,:]=thetabest
if False:
if sbbool:
thetapeaks[i,:]=findpeak(thetas,sbbool)
else:
thetapeaks[i,:7]=findpeak(thetas,sbbool)
np.save('data/MCMC/thetapeaks%s_r1.npy'%int(zb*1e3),thetapeaks)
np.save('data/MCMC/thetasall%s_r1.npy'%int(zb*1e3),thetasall)
np.save('data/MCMC/thetabests%s_r1.npy'%int(zb*1e3),thetabests)
np.save('data/MCMC/probsall%s_r1.npy'%int(zb*1e3),probsall)
np.save('data/MCMC/sbbools%s_r1.npy'%int(zb*1e3),sbbools)
plt.close()
thetasall=np.load('data/MCMC/thetasall%s_r1.npy'%int(zb*1e3))
probsall=np.load('data/MCMC/probsall%s_r1.npy'%int(zb*1e3))
thetabests=np.load('data/MCMC/thetabests%s_r1.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r1.npy'%int(zb*1e3))
if True:
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
nbins=20
# create data
v = thetasall[:,:,0].flatten()
w = thetasall[:,:,1].flatten()
x = thetasall[:,:,2].flatten()
y = thetasall[:,:,3].flatten()
# Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
k = scipy.stats.gaussian_kde([v,w,x,y])
vi, wi, xi, yi = np.mgrid[0.4:0.8:nbins*1j,0.5:1.8:nbins*1j,0.8:2.5:nbins*1j, 1.5:3.5:nbins*1j]
vg = np.linspace(0.4, 0.8, nbins)
wg = np.linspace(0.4, 1.8, nbins)
xg = np.linspace(0.8, 2.5, nbins)
yg = np.linspace(1.5, 3.5, nbins)
zi = k(np.vstack([vi.flatten(), wi.flatten(),xi.flatten(), yi.flatten()]))
zi=np.log(zi)
normval=np.max(zi)
zi=zi-normval
constvfun_4layer=interpolate.RegularGridInterpolator((vg,wg,xg,yg),zi.reshape((nbins,nbins,nbins,nbins)),fill_value=np.nan)
block = thetasall[sbbools.astype('bool'),:,7].flatten()
kblock = scipy.stats.gaussian_kde([block])
gblock = [np.linspace(1.2, 3.0, nbins)]
ziblock=kblock(gblock)
ziblock=np.log(ziblock)
normvalblock=np.max(ziblock)
ziblock=ziblock-normvalblock
constvfun_block=interpolate.RegularGridInterpolator(gblock,ziblock,fill_value=np.nan)
def constvfun(vs1,vs2,vs3,vs4,vsb):
return constvfun_block([vsb])+constvfun_4layer([vs1,vs2,vs3,vs4])
if True:
thetasall=np.zeros((384,500,9))
probsall=np.zeros((384,500))
thetabests=np.zeros((384,9))
thetapeaks=np.empty((384,9))
thetapeaks[:]=np.nan
sbbools=np.zeros(384)
for i in tqdm.tqdm(np.arange(0,384)):
startind=4*i
midind=4*i+20
endind=4*i+40
if (0<=(startind*2-25))&(((endind*2)-25)<len(meanpick)):
delaydata=np.zeros(3)
delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)])
delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)]))
delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)])
else:
delaydata=np.repeat(np.nan,3)
sbbool=False
if delaydata[0]==delaydata[0]:
sbbool=True
sbbools[i]=sbbool
cfm_use=curve_fillmean[startind:endind]
cfs_use=curve_fillstd[startind:endind]
dispcurve=np.average(cfm_use, axis=0, weights=1/cfs_use**2)
dispcurve_err=np.sqrt(np.average(cfs_use**2,axis=0,weights=1/cfs_use**2))
data = (dispcurve, dispcurve_err,delaydata,zb)
nwalkers=128
niter=500
initial=np.array((0.6, 0.7, 1.4, 2.9, 0.2, 0.2,0.2,2.4,0.6))
ndim = len(initial)
p0 = [np.array(initial) + 1e-1 * np.random.randn(ndim) for k in range(nwalkers)]
def main(p0,nwalkers,niter,ndim,lnprob,data):
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=data)
p0, _, _ = sampler.run_mcmc(p0, 100)
sampler.reset()
pos, prob, state = sampler.run_mcmc(p0, niter)
return sampler, pos, prob, state
sampler, pos, prob, state = main(p0,nwalkers,niter,ndim,lnprob_step2,data)
samples = sampler.flatchain
thetabest=samples[np.argmax(sampler.flatlnprobability)]
size=500
thetas=samples[np.random.randint(len(samples), size=size)]
thetasall[i,:,:]=thetas
probs=np.zeros(size)
for j in range(size):
if lnphysvalsprior(thetas[j],delaydata,zb)==0:
probs[j]=np.exp(lnprob_step2(thetas[j], dispcurve, dispcurve_err,delaydata,zb))
else:
probs[j]=0
probsall[i,:]=probs
thetabests[i,:]=thetabest
if True:
if sbbool:
thetapeaks[i,:]=findpeak(thetas,sbbool)
else:
thetapeaks[i,:7]=findpeak(thetas,sbbool)
np.save('data/MCMC/thetapeaks%s_r2.npy'%int(zb*1e3),thetapeaks)
np.save('data/MCMC/thetasall%s_r2.npy'%int(zb*1e3),thetasall)
np.save('data/MCMC/thetabests%s_r2.npy'%int(zb*1e3),thetabests)
np.save('data/MCMC/probsall%s_r2.npy'%int(zb*1e3),probsall)
np.save('data/MCMC/sbbools%s_r2.npy'%int(zb*1e3),sbbools)
plt.close()
thetasall=np.load('data/MCMC/thetasall%s_r2.npy'%int(zb*1e3))
probsall=np.load('data/MCMC/probsall%s_r2.npy'%int(zb*1e3))
thetabests=np.load('data/MCMC/thetabests%s_r2.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r2.npy'%int(zb*1e3))
0%| | 0/60 [00:00<?, ?it/s]/home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in double_scalars lnpdiff = f + nlp - state.log_prob[j] 2%|▋ | 1/60 [00:22<22:24, 22.78s/it]/tmp/ipykernel_3248220/4131823394.py:19: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/4131823394.py:20: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_3248220/4131823394.py:21: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) 58%|█████████████████████████ | 35/60 [13:44<09:32, 22.90s/it]/tmp/ipykernel_3248220/3802658083.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) 100%|███████████████████████████████████████████| 60/60 [31:12<00:00, 31.21s/it] 1%|▍ | 4/384 [03:28<5:29:12, 51.98s/it]/tmp/ipykernel_3248220/4131823394.py:124: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/4131823394.py:125: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_3248220/4131823394.py:126: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) 6%|██▍ | 23/384 [19:45<5:09:12, 51.39s/it]
probsums=np.zeros(6)
for zbi,zb in enumerate(np.linspace(0.1,0.6,6)):
thetasall=np.load('data/MCMC/thetasall%s_r2.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r2.npy'%int(zb*1e3))
thetasall[(1-sbbools).astype('bool'),:,7:]=np.nan
probsall=np.load('data/MCMC/probsall%s_r2.npy'%int(zb*1e3))
thetabests=np.zeros((thetasall.shape[0],thetasall.shape[2]))
thetasall[:,:,4:7][thetasall[:,:,4:7]<0.1]=0.1
thetasall[:,:,4][thetasall[:,:,4]>1]=1
thetasall[:,:,5][thetasall[:,:,5]>2]=2
thetasall[:,:,6][thetasall[:,:,6]>3]=3
thetabests=np.load('data/MCMC/thetapeaks%s_r2.npy'%int(zb*1e3))
for i in range(len(thetabests)):
if thetabests[i][0]!=0:
close_i=np.argmin(np.nansum(abs(thetasall[i]-thetabests[i]),axis=1))
thetabests[i]=thetasall[i,close_i]
thetabests[:,4:7][thetabests[:,4:7]<0.1]=0.1
thetabests[:,4][thetabests[:,4]>1]=1
thetabests[:,5][thetabests[:,5]>2]=2
thetabests[:,6][thetabests[:,6]>3]=3
thetabests[(1-sbbools).astype('bool'),7:]=np.nan
dcurve_obs=np.zeros((384,50))
dcurve_obs_err=np.zeros((384,50))
for chan in range(384):
thetabest=thetabests[chan].copy()
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
startind=4*chan
midind=4*chan+20
endind=4*chan+40
if (0<=(startind*2-25))&(((endind*2)-25)<len(meanpick)):
delaydata=np.zeros(3)
delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)])
delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)]))
delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)])
else:
delaydata=np.repeat(np.nan,3)
sbbool=False
if delaydata[0]==delaydata[0]:
sbbool=True
#predict dispersion curves along model
cfm_use=curve_fillmean[startind:endind]
cfs_use=curve_fillstd[startind:endind]
dcurve_obs[chan,:]=np.average(cfm_use, axis=0,weights=1/cfs_use**2)
dcurve_obs_err[chan,:]=np.sqrt(np.average(cfs_use**2,axis=0,weights=1/cfs_use**2))
if sbbool==False:
thetabest[thetabest!=thetabest]=1.5
probs[chan]=lnprob_step1(thetabest, dcurve_obs[chan], dcurve_obs_err[chan],delaydata,zb)
probsums[zbi]=np.nansum(probs)
/tmp/ipykernel_3248220/2934882799.py:36: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/2934882799.py:37: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_3248220/2934882799.py:38: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/3802658083.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
zb=np.linspace(0.1,0.6,6)[np.argmax(probsums)]
zb
0.2
zb=np.linspace(0.1,0.6,6)[np.argmax(probsums)]
#pick highest likelihood model
np.sum(probs)
import tqdm
import scipy.stats
#read in results
thetasall=np.load('data/MCMC/thetasall%s_r2.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r2.npy'%int(zb*1e3))
thetasall[(1-sbbools).astype('bool'),:,7:]=np.nan
probsall=np.load('data/MCMC/probsall%s_r2.npy'%int(zb*1e3))
thetabests=np.zeros((thetasall.shape[0],thetasall.shape[2]))
thetasall[:,:,4:7][thetasall[:,:,4:7]<0.1]=0.1
thetasall[:,:,4][thetasall[:,:,4]>1]=1
thetasall[:,:,5][thetasall[:,:,5]>2]=2
thetasall[:,:,6][thetasall[:,:,6]>3]=3
thetabests=np.load('data/MCMC/thetapeaks%s_r2.npy'%int(zb*1e3))
for i in range(len(thetabests)):
if thetabests[i][0]!=0:
close_i=np.argmin(np.nansum(abs(thetasall[i]-thetabests[i]),axis=1))
thetabests[i]=thetasall[i,close_i]
thetabests[:,4:7][thetabests[:,4:7]<0.1]=0.1
thetabests[:,4][thetabests[:,4]>1]=1
thetabests[:,5][thetabests[:,5]>2]=2
thetabests[:,6][thetabests[:,6]>3]=3
thetabests[(1-sbbools).astype('bool'),7:]=np.nan
#where converted phase is not present, set model to np.nan
#define a moving average function which ignores nans and outliers
def nanaverage(a):
w=np.sin(np.linspace(0,np.pi,len(a)))**2
return np.nansum(a*w)/((~np.isnan(a))*w).sum()
'''
def my_ave(a):
sd=np.nanstd(a)
m=np.nanmean(a)
w=np.sin(np.linspace(0,np.pi,len(a)))**2
return np.nanmean(a[(a<(m+2*sd))&(a>(m-2*sd))],weights=w)
'''
def moving_average(a, n=15) :
return np.array([nanaverage(a[i:i+n])for i in range(len(a)-n+1)])
#make vector of slumped block thicknesses, which goes to zeros where the block is absent
sbts=np.zeros(len(sbbools))
sbts[(sbbools).astype('bool')]=zb
sbts_smooth=moving_average(sbts)
#smoothed model
thetabests_smooth=np.apply_along_axis(moving_average,0,thetabests)
#find array channels associated with model columns and smooth result
midinds_fullarr=(4*(np.arange(0,384))+20)*2+5999
midinds_fullarr_smooth=moving_average(midinds_fullarr)
#read array location to get elevations & smooth elevations
chan_csv=pd.read_csv('data/Mammoth_south_100km.csv')
#get along-fiber distance of channels from channel 5999
ch_lat=np.array(chan_csv[chan_csv['status']=='good']['latitude'])
ch_lon=np.array(chan_csv[chan_csv['status']=='good']['longitude'])
_,_,zno,zlt=utm.from_latlon(np.mean(ch_lat),np.mean(ch_lon))
ch_x,ch_y,_,_=utm.from_latlon(ch_lat,ch_lon,zno,zlt)
along_fiber_dist=np.cumsum(np.sqrt(np.concatenate(([0],np.diff(ch_x)))**2+np.concatenate(([0],np.diff(ch_y)))**2))/1e3
distfrom5999=(along_fiber_dist-along_fiber_dist[5999])
#along -fiber distance in plotted model
invdist=distfrom5999[midinds_fullarr]
invdist_smooth=moving_average(invdist)
elevations=np.array(chan_csv[chan_csv['status']=='good']['elevation'])[midinds_fullarr]/1000
elevations=-elevations+np.max(elevations)
elevations_smooth=moving_average(elevations)
#define blank matrices to plot models velmodmat is smooth and velmodmatr is rough
velmodmat=np.zeros((370,200))
velmodmatr=np.zeros((384,200))
dep_ax=np.linspace(0,3.5,200)
#standard deviation in model parameters
thetastds=np.std(thetasall,axis=1)
thetastds_smooth=np.sqrt(np.apply_along_axis(moving_average,0,thetastds**2))
#define blank matrices to plot velocity uncertainties, modstdmat is smooth, modstdmatr is rough
modstdmat=np.zeros((370,200))
modstdmatr=np.zeros((384,200))
for chan in range(384):
thetabest=thetabests[chan]
thetastd=thetastds[chan]
ele=elevations[chan]
sbt=sbts[chan]
sbbool=bool(sbbools[chan])
depsums=np.cumsum(np.concatenate(([ele],thetabest[4:7])))
velmodmatr[chan,dep_ax<=depsums[0]]=np.nan
velmodmatr[chan,(dep_ax>depsums[0])&(dep_ax<=depsums[1])]=thetabest[0]
velmodmatr[chan,(dep_ax>depsums[1])&(dep_ax<=depsums[2])]=thetabest[1]
velmodmatr[chan,(dep_ax>depsums[2])&(dep_ax<=depsums[3])]=thetabest[2]
velmodmatr[chan,dep_ax>depsums[3]]=thetabest[3]
modstdmatr[chan,dep_ax<=depsums[0]]=np.nan
modstdmatr[chan,(dep_ax>depsums[0])&(dep_ax<=depsums[1])]=thetastd[0]
modstdmatr[chan,(dep_ax>depsums[1])&(dep_ax<=depsums[2])]=thetastd[1]
modstdmatr[chan,(dep_ax>depsums[2])&(dep_ax<=depsums[3])]=thetastd[2]
modstdmatr[chan,dep_ax>depsums[3]]=thetastd[3]
if sbbool:
velmodmatr[chan,(dep_ax>thetabest[8]+ele)&(dep_ax<thetabest[8]+sbt+ele)]=thetabest[7]
modstdmatr[chan,(dep_ax>thetabest[8]+ele)&(dep_ax<thetabest[8]+sbt+ele)]=thetastd[7]
#make smooth model and uncertainty matrix for plotting
for chan in range(370):
thetabest=thetabests_smooth[chan]
thetastd=thetastds_smooth[chan]
ele=elevations_smooth[chan]
sbt=moving_average(sbts)[chan]
depsums=np.cumsum(np.concatenate(([ele],thetabest[4:7])))
velmodmat[chan,dep_ax<=depsums[0]]=np.nan
velmodmat[chan,(dep_ax>depsums[0])&(dep_ax<=depsums[1])]=thetabest[0]
velmodmat[chan,(dep_ax>depsums[1])&(dep_ax<=depsums[2])]=thetabest[1]
velmodmat[chan,(dep_ax>depsums[2])&(dep_ax<=depsums[3])]=thetabest[2]
velmodmat[chan,dep_ax>depsums[3]]=thetabest[3]
modstdmat[chan,dep_ax<=depsums[0]]=np.nan
modstdmat[chan,(dep_ax>depsums[0])&(dep_ax<=depsums[1])]=thetastd[0]
modstdmat[chan,(dep_ax>depsums[1])&(dep_ax<=depsums[2])]=thetastd[1]
modstdmat[chan,(dep_ax>depsums[2])&(dep_ax<=depsums[3])]=thetastd[2]
modstdmat[chan,dep_ax>depsums[3]]=thetastd[3]
if sbt>0:
velmodmat[chan,(dep_ax>thetabest[8]+ele)&(dep_ax<thetabest[8]+sbt+ele)]=thetabest[7]
modstdmat[chan,(dep_ax>thetabest[8]+ele)&(dep_ax<thetabest[8]+sbt+ele)]=thetastd[7]
#now determine the data fit of the rough model
dep_ax=np.linspace(0,3.5,200)
dcurve_mod=np.zeros((384,50))
t_residuals=np.repeat(np.nan,384)
t_obs=np.repeat(np.nan,384)
t_mod=np.repeat(np.nan,384)
amp_obs=np.repeat(np.nan,384)
amp_err=np.repeat(np.nan,384)
amp_mod=np.repeat(np.nan,384)
ps=np.repeat(np.nan,384)
dcurve_obs=np.zeros((384,50))
dcurve_obs_err=np.zeros((384,50))
t_err=np.repeat(np.nan,384)
probs=np.repeat(np.nan,384)
f1, f2 = 0.25, 4.0
periods=1/np.logspace(np.log10(f1),np.log10(f2),50)
for chan in range(384):
thetabest=thetabests[chan].copy()
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
ele=elevations[chan]
startind=4*chan
midind=4*chan+20
endind=4*chan+40
if (0<=(startind*2-25))&(((endind*2)-25)<len(meanpick)):
delaydata=np.zeros(3)
delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)])
delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)]))
delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)])
else:
delaydata=np.repeat(np.nan,3)
sbbool=False
if delaydata[0]==delaydata[0]:
sbbool=True
#predict dispersion curves along model
dcurve_mod[chan,:]=pred_dcurve(thetabest,sb=sbbool,zb=zb)
cfm_use=curve_fillmean[startind:endind]
cfs_use=curve_fillstd[startind:endind]
dcurve_obs[chan,:]=np.average(cfm_use, axis=0,weights=1/cfs_use**2)
dcurve_obs_err[chan,:]=np.sqrt(np.average(cfs_use**2,axis=0,weights=1/cfs_use**2))
if dcurve_mod[chan,0]!=dcurve_mod[chan,0]:
print(chan)
if sbbool:
t_obs[chan]=delaydata[0]
t_err[chan]=delaydata[1]
ps[chan]=delaydata[2]
#predict time delay along model
t_mod[chan]=delaytime(thetabest,ps[chan])
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
depth_cum=np.cumsum((d1,d2,d3,5))
#find layer block starts in
startblock=np.where((depth_cum>db))[0][0]
betas_t=np.array((vs1,vs2,vs3,vs4))
betas=np.concatenate((betas_t[:(startblock+1)],[vsb]))
n=len(betas)
alphas=vs2vp(betas)
rhos=(alphas**0.25)*1.74
direct=np.zeros(n)
conv=np.ones(n)
conv[0]=0
#predict amplitude ratio along model
amp_mod[chan]=predamprat(p_overall,alphas,betas,rhos,direct=direct,conv=conv)
if amp_mod[chan]!=amp_mod[chan]:
print(chan)
amp_obs[chan]=rat
amp_err[chan]=rat_err
else:
thetabest[thetabest!=thetabest]=1.5
probs[chan]=lnprob_step1(thetabest, dcurve_obs[chan], dcurve_obs_err[chan],delaydata,zb)
/tmp/ipykernel_3248220/2658180458.py:35: RuntimeWarning: invalid value encountered in double_scalars return np.nansum(a*w)/((~np.isnan(a))*w).sum() /tmp/ipykernel_3248220/2658180458.py:35: RuntimeWarning: invalid value encountered in double_scalars return np.nansum(a*w)/((~np.isnan(a))*w).sum() /tmp/ipykernel_3248220/2658180458.py:160: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/2658180458.py:161: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_3248220/2658180458.py:162: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_3248220/3802658083.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
#here, we are determining the data fit of the smoothed model
t_obs_smooth=np.apply_along_axis(moving_average,0,dcurve_obs)
t_mod_smooth=np.repeat(np.nan,370)
t_obs_smooth=np.apply_along_axis(moving_average,0,t_obs)
t_err_smooth=np.apply_along_axis(moving_average,0,t_err)
amp_obs_smooth=np.apply_along_axis(moving_average,0,amp_obs)
amp_mod_smooth=np.repeat(np.nan,370)
amp_err_smooth=np.apply_along_axis(moving_average,0,amp_err)
dcurve_obs_smooth=np.apply_along_axis(moving_average,0,dcurve_obs)
dcurve_obs_err_smooth=np.apply_along_axis(moving_average,0,dcurve_obs_err)
dcurve_mod_smooth=np.zeros((370,50))
ps_smooth=np.apply_along_axis(moving_average,0,ps)
probs_smooth=np.repeat(np.nan,370)
for chan in range(370):
sbt=sbts_smooth[chan]
sbbool=sbt>0
thetabest=thetabests_smooth[chan].copy()
#predict dispersion curve along model
dcurve_mod_smooth[chan,:]=pred_dcurve(thetabest,sb=sbbool,zb=sbt)
if sbbool:
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
depth_cum=np.cumsum((d1,d2,d3,5))
#find layer block starts in
startblock=np.where((depth_cum>db))[0][0]
betas_t=np.array((vs1,vs2,vs3,vs4))
betas=np.concatenate((betas_t[:(startblock+1)],[vsb]))
#predict delay along model
t_mod_smooth[chan]=delaytime(thetabest,ps_smooth[chan])
n=len(betas)
alphas=vs2vp(betas)
rhos=(alphas**0.25)*1.74
direct=np.zeros(n)
conv=np.ones(n)
conv[0]=0
#predict amplitude along model
amp_mod_smooth[chan]=predamprat(p_overall,alphas,betas,rhos,direct=direct,conv=conv)
delaydata=np.array((t_obs_smooth[chan],t_err_smooth[chan],ps_smooth[chan]))
thetabest[thetabest!=thetabest]=1.5
probs_smooth[chan]=lnprob_step1(thetabest, dcurve_obs_smooth[chan], dcurve_obs_err_smooth[chan],delaydata,zb)
#pick highest probability model:
#in picking the highest probability smoothed model,
#we ignore the transition between
#the areas of the model with and without the block
np.sum(probs_smooth[(sbts_smooth==0)|(sbts_smooth==zb)])
/tmp/ipykernel_3248220/2658180458.py:35: RuntimeWarning: invalid value encountered in double_scalars return np.nansum(a*w)/((~np.isnan(a))*w).sum() /tmp/ipykernel_3248220/3802658083.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
-219.87545387159744
plt.rcParams.update({'font.size': 44})
dcurve_res_smooth=abs(dcurve_obs_smooth-dcurve_mod_smooth)/dcurve_obs_err_smooth
fig, ax1=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2)
mod=ax1.pcolormesh(invdist_smooth,periods,dcurve_mod_smooth.T,vmax=2.5,cmap='Spectral')
fig.colorbar(mod,label='Modeled Velocity (km/s)')
ax1.set_yscale('log')
ax1.set_ylim((0.4,4))
ax1.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax1.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_mod.png')
fig, ax2=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2)
obs=ax2.pcolormesh(invdist_smooth,periods,dcurve_obs_smooth.T,vmax=2.5,cmap='Spectral')
ax2.set_ylim((0.4,4))
fig.colorbar(obs,label='Observed Velocity (km/s)')
ax2.set_yscale('log')
ax2.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax2.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_obs.png')
fig, ax3=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2)
res=ax3.pcolormesh(invdist_smooth,periods,dcurve_res_smooth.T,vmax=2,cmap='Spectral')
fig.colorbar(res,label='Normalized Residual (in s.d.)')
ax3.set_ylim((0.4,4))
ax3.set_yscale('log')
ax3.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax3.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_res.png')
fig, ax4=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax4.plot(invdist_smooth,t_obs_smooth,linewidth=6)
ax4.plot(invdist_smooth,t_mod_smooth,linewidth=6,color='k')
ax4.fill_between(invdist_smooth,t_obs_smooth-t_err_smooth,t_obs_smooth+t_err_smooth,alpha=0.25)
ax4.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax4.set_ylabel('Delay between P and converted phase (s)')
plt.legend(('Observed','Modeled'))
plt.savefig('figs/tfit.png')
fig, ax5=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax5.plot(invdist_smooth,amp_obs_smooth,linewidth=6)
ax5.plot(invdist_smooth,amp_mod_smooth,color='k',linewidth=6)
ax5.fill_between(invdist_smooth,amp_obs_smooth-amp_err_smooth,amp_obs_smooth+amp_err_smooth,alpha=0.25)
ax5.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax5.set_ylabel('Amplitude Ratio between P and Ps')
plt.savefig('figs/ampfit.png')
fig, ((ax4,ax2),(ax1,ax3))=plt.subplots(2,2,figsize=(40,32))
plt.tight_layout()
mod=ax1.pcolormesh(invdist_smooth,periods,dcurve_mod_smooth.T,vmax=2.5,cmap='Spectral')
fig.colorbar(mod,label='Modeled Velocity (km/s)')
ax1.set_yscale('log')
ax1.set_ylim((0.4,4))
ax1.set_xlabel('Along-fiber Distance (km)')
ax1.set_ylabel('Period (s)')
obs=ax2.pcolormesh(invdist_smooth,periods,dcurve_obs_smooth.T,vmax=2.5,cmap='Spectral')
ax2.set_ylim((0.4,4))
fig.colorbar(obs,label='Observed Velocity (km/s)')
ax2.set_yscale('log')
ax2.set_xlabel('Along-fiber Distance (km)')
ax2.set_ylabel('Period (s)')
res=ax3.pcolormesh(invdist_smooth,periods,dcurve_res_smooth.T,vmax=2,cmap='Spectral')
fig.colorbar(res,label='Normalized Residual (in s.d.)')
ax3.set_ylim((0.4,4))
ax3.set_yscale('log')
ax3.set_xlabel('Along-fiber Distance (km)')
ax3.set_ylabel('Period (s)')
ax4.plot(invdist_smooth,t_obs_smooth,linewidth=6)
ax4.plot(invdist_smooth,t_mod_smooth,linewidth=6,color='k')
ax4.fill_between(invdist_smooth,t_obs_smooth-t_err_smooth,t_obs_smooth+t_err_smooth,alpha=0.25)
ax4.set_xlim(np.min(invdist_smooth[t_mod_smooth==t_mod_smooth]),np.max(invdist_smooth[t_mod_smooth==t_mod_smooth]))
ax4.set_xlabel('Along-fiber Distance (km)')
ax4.set_ylabel('Delay between P and converted phase (s)')
ax4.legend(('Observed','Modeled'))
plt.tight_layout()
plt.savefig('figs/datafit.png')
plt.rcParams.update({'font.size': 42})
dcurve_res=abs(dcurve_obs-dcurve_mod)/dcurve_obs_err
fig, ax1=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2,bottom=0.2)
mod=ax1.pcolormesh(invdist,periods,dcurve_mod.T,vmax=2.5,cmap='Spectral')
fig.colorbar(mod,label='Modeled Velocity (km/s)')
ax1.set_yscale('log')
ax1.set_ylim((0.4,4))
ax1.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax1.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_modr.png')
fig, ax2=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2)
obs=ax2.pcolormesh(invdist,periods,dcurve_obs.T,vmax=2.5,cmap='Spectral')
ax2.set_ylim((0.4,4))
fig.colorbar(obs,label='Observed Velocity (km/s)')
ax2.set_yscale('log')
ax2.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax2.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_obsr.png')
fig, ax3=plt.subplots(1,1,figsize=(24,16))
plt.subplots_adjust(left=0.2)
res=ax3.pcolormesh(invdist,periods,dcurve_res.T,vmax=2,cmap='Spectral')
fig.colorbar(res,label='Normalized Residual (in s.d.)')
ax3.set_ylim((0.4,4))
ax3.set_yscale('log')
ax3.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax3.set_ylabel('Period (s)')
plt.savefig('figs/dcurve_resr.png')
fig, ax4=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax4.plot(invdist,t_obs,linewidth=6)
ax4.plot(invdist,t_mod,linewidth=6,color='k')
ax4.fill_between(invdist,t_obs-t_err,t_obs+t_err,alpha=0.25)
ax4.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax4.set_ylabel('Delay between P and converted phase (s)')
plt.legend(('Observed','Modeled'))
plt.savefig('figs/tfitr.png')
fig, ax5=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax5.plot(invdist,amp_obs,linewidth=6)
ax5.plot(invdist,amp_mod,color='k',linewidth=6)
ax5.fill_between(invdist,amp_obs-amp_err,amp_obs+amp_err,alpha=0.25)
ax5.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax5.set_ylabel('Amplitude Ratio between P and Ps')
plt.savefig('figs/ampfitr.png')
fig, ax5=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax5.plot(invdist_smooth,amp_obs_smooth,linewidth=6)
ax5.plot(invdist_smooth,amp_mod_smooth,color='k',linewidth=6)
ax5.fill_between(invdist_smooth,amp_obs_smooth-amp_err_smooth,amp_obs_smooth+amp_err_smooth,alpha=0.25)
ax5.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax5.set_ylabel('Amplitude Ratio between P and Ps')
plt.savefig('figs/ampfit.png')
fig, ax6=plt.subplots(1,1,figsize=(20,16))
plt.subplots_adjust(left=0.2)
ax6.plot(invdist,amp_obs,linewidth=6,label='observation')
ax6.plot(invdist,amp_mod,color='k',linewidth=6,label='prediction')
ax6.plot(invdist_smooth,amp_mod_smooth,color='r',linewidth=6,label='prediction (smoothed)')
ax6.fill_between(invdist,amp_obs-amp_err,amp_obs+amp_err,alpha=0.25)
ax6.set_xlabel('Along-fiber Distance from Channel 5999 (km)')
ax6.set_ylabel('Amplitude Ratio between P and Ps')
plt.legend()
plt.savefig('figs/ampfitfull.png')
fig, ((ax4,ax2),(ax1,ax3))=plt.subplots(2,2,figsize=(40,32))
mod=ax1.pcolormesh(invdist,periods,dcurve_mod.T,vmax=2.5,cmap='Spectral')
fig.colorbar(mod,label='Modeled Velocity (km/s)')
ax1.set_yscale('log')
ax1.set_ylim((0.4,4))
ax1.set_xlabel('Along-fiber Distance (km)')
ax1.set_ylabel('Period (s)')
obs=ax2.pcolormesh(invdist,periods,dcurve_obs.T,vmax=2.5,cmap='Spectral')
ax2.set_ylim((0.4,4))
fig.colorbar(obs,label='Observed Velocity (km/s)')
ax2.set_yscale('log')
ax2.set_xlabel('Along-fiber Distance (km)')
ax2.set_ylabel('Period (s)')
res=ax3.pcolormesh(invdist,periods,dcurve_res.T,vmax=2,cmap='Spectral')
fig.colorbar(res,label='Normalized Residual (in s.d.)')
ax3.set_ylim((0.4,4))
ax3.set_yscale('log')
ax3.set_xlabel('Along-fiber Distance (km)')
ax3.set_ylabel('Period (s)')
ax4.plot(invdist,t_obs,linewidth=6)
ax4.plot(invdist,t_mod,linewidth=6,color='k')
ax4.fill_between(invdist,t_obs-t_err,t_obs+t_err,alpha=0.25)
ax4.set_xlim(np.min(invdist[t_mod==t_mod]),np.max(invdist[t_mod==t_mod]))
ax4.set_xlabel('Along-fiber Distance (km)')
ax4.set_ylabel('Delay between P and converted phase (s)')
plt.legend(('Observed','Modeled'))
plt.tight_layout()
plt.savefig('figs/datafitr.png')
def forceAspect(ax,aspect=1):
im = ax.get_images()
extent = im[0].get_extent()
ax.set_aspect(abs((extent[1]-extent[0])/(extent[3]-extent[2]))/aspect)
plt.rcParams.update({'font.size': 60})
xl_smooth=invdist_smooth[-1]-invdist_smooth[0]
xl=invdist[-1]-invdist[0]
fig, ax1=plt.subplots(1,1,figsize=(50,16))
plt.tight_layout()
velmod=ax1.pcolormesh(invdist_smooth,dep_ax,velmodmat.T,cmap='Spectral')
ax1.set_ylabel('depth (km)')
ax1.set_xlabel('Along-fiber Distance (km)')
ax1.set_ylim((2.5,0))
ax1.set_aspect(4)
fig.colorbar(velmod,label='Velocity (km/s)')
plt.savefig('figs/vmod.png')
fig, ax2=plt.subplots(1,1,figsize=(50,16))
plt.tight_layout()
velstd=ax2.pcolormesh(invdist_smooth,dep_ax,(modstdmat/velmodmat).T,cmap='Spectral')
ax2.set_ylabel('depth (km)')
ax2.set_xlabel('Along-fiber Distance (km)')
ax2.set_ylim((2.5,0))
ax2.set_aspect(4)
fig.colorbar(velstd,label='Velocity s.d./\n Velocity (-)')
plt.savefig('figs/vstd.png')
fig, ax3=plt.subplots(1,1,figsize=(50,16))
plt.subplots_adjust(left=0.2,bottom=0.2)
deps4plt=np.concatenate((np.repeat(elevations,500).reshape(-1,500,1),thetasall[:,:,4:7]),axis=2)
deps4plt_cum=np.cumsum(deps4plt,axis=2)
deps_plt=np.cumsum(np.concatenate(([elevations_smooth],thetabests_smooth[:,4:7].T),axis=0),axis=0)
deps_pltm=np.apply_along_axis(moving_average,0,np.percentile(deps4plt_cum,20,axis=1)).T
deps_pltp=np.apply_along_axis(moving_average,0,np.percentile(deps4plt_cum,80,axis=1)).T
blockdeps_plt=(thetabests_smooth[:,8]+elevations_smooth).T
blockdeps_pltm=np.apply_along_axis(moving_average,0,np.percentile(thetasall[:,:,8],20,axis=1)+elevations).T
blockdeps_pltp=np.apply_along_axis(moving_average,0,np.percentile(thetasall[:,:,8],80,axis=1)+elevations).T
for i in range(4):
ax3.plot(invdist_smooth,deps_plt[i,:],color='blue',linewidth=6)
if i==0:
ax3.plot(invdist_smooth,deps_plt[i,:],color='blue',linewidth=6,label='layer boundary')
ax3.fill_between(invdist_smooth,deps_pltm[i,:],deps_pltp[i,:],color='blue',alpha=0.2)
ax3.plot(invdist_smooth,blockdeps_plt,color='red',linewidth=6,label='top of slumped block')
ax3.fill_between(invdist_smooth,blockdeps_pltm,blockdeps_pltp,color='red',alpha=0.2)
ax3.set_ylim((2.5,0))
ax3.set_ylabel('depth (km)')
ax3.set_xlabel('Along-fiber Distance (km)')
plt.legend()
ax3.set_aspect(4)
plt.savefig('figs/dstd.png')
fig, ax4=plt.subplots(1,1,figsize=(50,16))
plt.tight_layout()
velmodr=ax4.pcolormesh(invdist,dep_ax,velmodmatr.T,cmap='Spectral')
ax4.set_ylabel('depth (km)')
ax4.set_xlabel('Along-fiber Distance (km)')
ax4.set_ylim((2.5,0))
fig.colorbar(velmodr,label='Velocity (km/s)')
ax4.set_aspect(4)
plt.savefig('figs/vmodr.png')
fig, ax5=plt.subplots(1,1,figsize=(50,16))
plt.tight_layout()
velstd=ax5.pcolormesh(invdist,dep_ax,(modstdmatr/velmodmatr).T,cmap='Spectral')
ax5.set_ylabel('depth (km)')
ax5.set_xlabel('Along-fiber Distance (km)')
ax5.set_ylim((2.5,0))
fig.colorbar(velstd,label='Velocity s.d./\n Velocity (-)')
ax5.set_aspect(4)
plt.savefig('figs/vstdr.png')
fig, ax6=plt.subplots(1,1,figsize=(50,16))
plt.subplots_adjust(left=0.2,bottom=0.2)
deps4plt=np.concatenate((np.repeat(elevations,500).reshape(-1,500,1),thetasall[:,:,4:7]),axis=2)
deps4plt_cum=np.cumsum(deps4plt,axis=2)
deps_plt=np.cumsum(np.concatenate(([elevations],thetabests[:,4:7].T),axis=0),axis=0)
deps_pltm=np.percentile(deps4plt_cum,20,axis=1).T
deps_pltp=np.percentile(deps4plt_cum,80,axis=1).T
blockdeps_plt=(thetabests[:,8]+elevations).T
blockdeps_pltm=(np.percentile(thetasall[:,:,8],20,axis=1)+elevations).T
blockdeps_pltp=(np.percentile(thetasall[:,:,8],80,axis=1)+elevations).T
for i in range(4):
ax6.plot(invdist,deps_plt[i,:],color='blue',linewidth=6)
if i==0:
ax6.plot(invdist,deps_plt[i,:],color='blue',linewidth=6,label='layer boundary')
ax6.fill_between(invdist,deps_pltm[i,:],deps_pltp[i,:],color='blue',alpha=0.2)
ax6.plot(invdist,blockdeps_plt,color='red',linewidth=6,label='top of slumped block')
ax6.fill_between(invdist,blockdeps_pltm,blockdeps_pltp,color='red',alpha=0.2)
ax6.set_ylim((2.5,0))
ax6.set_ylabel('depth (km)')
ax6.set_xlabel('Along-fiber Distance (km)')
ax6.set_aspect(4)
plt.savefig('figs/dstdr.png')
/tmp/ipykernel_3248220/2658180458.py:35: RuntimeWarning: invalid value encountered in double_scalars return np.nansum(a*w)/((~np.isnan(a))*w).sum()
plt.rcParams.update({'font.size': 44})
depth_pars=thetasall[:,:,np.array((4,5,6,8))]
depth5=np.nanpercentile(depth_pars,5,axis=1)
depth25=np.nanpercentile(depth_pars,25,axis=1)
depth50=np.nanpercentile(depth_pars,50,axis=1)
depth75=np.nanpercentile(depth_pars,75,axis=1)
depth95=np.nanpercentile(depth_pars,95,axis=1)
depths=thetabests[:,np.array((4,5,6,8))]
depth5_smooth=np.apply_along_axis(moving_average,0,depth5)
depth25_smooth=np.apply_along_axis(moving_average,0,depth25)
depth50_smooth=np.apply_along_axis(moving_average,0,depth50)
depth75_smooth=np.apply_along_axis(moving_average,0,depth75)
depth95_smooth=np.apply_along_axis(moving_average,0,depth95)
depths_smooth=np.apply_along_axis(moving_average,0,depths)
axes=np.zeros((2,2))
fig,axes = plt.subplots(2,2,figsize=(40,30))
for plot_i in range(4):
i=np.array((0,0,1,1))[plot_i]
j=np.array((0,1,0,1))[plot_i]
axes[i,j].plot(invdist_smooth,depths_smooth[:,plot_i],'k',label='best fit',linewidth=6)
axes[i,j].plot(invdist_smooth,depth50_smooth[:,plot_i],'b',label='median',linewidth=6)
axes[i,j].fill_between(invdist_smooth,depth5_smooth[:,plot_i],depth25_smooth[:,plot_i],color='b',alpha=0.2,label='5-95 percentile')
axes[i,j].fill_between(invdist_smooth,depth25_smooth[:,plot_i],depth75_smooth[:,plot_i],color='b',alpha=0.4,label='25-75 percentile')
axes[i,j].fill_between(invdist_smooth,depth75_smooth[:,plot_i],depth95_smooth[:,plot_i],color='b',alpha=0.2)
axes[i,j].set_xlabel('Along-fiber Distance (km)')
plt.xlim(np.min(invdist)-0.5,np.max(invdist)+0.5)
if plot_i==0:
axes[i,j].legend(loc='upper left')
if plot_i<3:
axes[i,j].set_ylabel('layer %s thickness (km)'%(plot_i+1))
else:
axes[i,j].set_ylabel('slumped block depth (km)')
plt.savefig('figs/dlayer.png')
axes=np.zeros((2,2))
fig,axes = plt.subplots(2,2,figsize=(40,30))
for plot_i in range(4):
i=np.array((0,0,1,1))[plot_i]
j=np.array((0,1,0,1))[plot_i]
axes[i,j].plot(invdist,depths[:,plot_i],'k',label='best fit',linewidth=6)
axes[i,j].plot(invdist,depth50[:,plot_i],'b',label='median',linewidth=6)
axes[i,j].fill_between(invdist,depth5[:,plot_i],depth25[:,plot_i],color='b',alpha=0.2,label='5-95 percentile')
axes[i,j].fill_between(invdist,depth25[:,plot_i],depth75[:,plot_i],color='b',alpha=0.4,label='25-75 percentile')
axes[i,j].fill_between(invdist,depth75[:,plot_i],depth95[:,plot_i],color='b',alpha=0.2)
axes[i,j].set_xlabel('Along-fiber Distance (km)')
plt.xlim(np.min(invdist)-0.5,np.max(invdist)+0.5)
if plot_i==0:
axes[i,j].legend(loc='upper left')
if plot_i<3:
axes[i,j].set_ylabel('layer %s thickness (km)'%(plot_i+1))
else:
axes[i,j].set_ylabel('slumped block depth (km)')
plt.savefig('figs/dlayerr.png')
/home/elibird/.local/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1584: RuntimeWarning: All-NaN slice encountered result = np.apply_along_axis(_nanquantile_1d, axis, a, q, /tmp/ipykernel_3248220/2658180458.py:35: RuntimeWarning: invalid value encountered in double_scalars return np.nansum(a*w)/((~np.isnan(a))*w).sum()
onedmod=(np.nanmean(thetabests,axis=0)*1e3).round(-1).astype('int')
onedmoderr=(np.nanstd(thetabests,axis=0)*1e3).round(-1).astype('int')
str_res='Model residuals show the greatest reduction for a slumped block of %s meter thickness. '%int(zb*1e3)
str_res+='The block takes on a velocity of %s m/s with a standard deviation of %s m/s. '%(onedmod[7],onedmoderr[7])
str_res+='In order of increasing depth, the remaining model layers take on mean velocities of '
for i in range(3):
str_res+='%s, '%(onedmod[i])
str_res+='and %s m/s, '%(onedmod[3])
str_res+='with respective standard deviations of '
for i in range(3):
str_res+='%s, '%(onedmoderr[i])
str_res+='and %s m/s. '%(onedmoderr[3])
startblockinds=np.where(sbbools.astype('bool'))[0][:10]
startdep=int(round(np.mean(thetabests[startblockinds,8]*1e3),-1))
startdeperr=int(round(np.std(thetasall[startblockinds,:,8]*1e3),-1))
str_res+='As the array crosses onto the proposed location of the NBB, the inversion indicates the top of the slumped block is located $%s\pm%s$ meters beneath the ground surface. '%(startdep,startdeperr)
str_res+=' This depth slightly exceeds that proposed for the block on the basis of resistivity data. This disagreement could be explained by lateral heterogeneity not accounted for in the one-dimensional resistivity-based estimate of block depth. '
pastblockinds=(midinds_fullarr>8500)&sbbools.astype('bool')
pastdep=int(round(np.mean(thetabests[pastblockinds,8]*1e3),-1))
pastdeperr=int(round(np.std(thetasall[pastblockinds,:,8]*1e3),-1))
str_res+='As the array moves to the south, off of the proposed location, delays between the {\it P} and {\it Ps} arrivals increase, and in turn the modeled depth of the converting boundary increases to $%s\pm%s$. '%(pastdep,pastdeperr)