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/DispCurves/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 False:
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')
ccurve_fillstd=np.load('data/curve_fillstd.npy')
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 False:
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')
#curve_fillmean[:,:5][curve_fillmean[:,:5]<1]=np.nan
plt.pcolormesh(frqs,np.arange(1574),curve_fillmean)
plt.colorbar()
plt.xscale('log')
#Plotting Example Dispersion Curve
fig,ax=plt.subplots(1,figsize=(20,20))
ich=630
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.3,0.6,1,2,3,4))
ax.set_xlabel('Period (s)')
ax.set_ylabel('Rayleigh Phase Velocity (m/s)')
plt.savefig('figs/egdcurve.png')
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[6], line 4 2 fig,ax=plt.subplots(1,figsize=(20,20)) 3 ich=630 ----> 4 dstack=(data_stacks[ich]/np.max(data_stacks[ich],axis=0)) 5 ax.pcolormesh(1/frqs,vels,dstack,cmap='Spectral_r') 6 #ax.plot(1/frqs,curve_fillmean[ich]*1e3,'k',linewidth=8) NameError: name 'data_stacks' is not defined
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.nansum(((y[fs<2.5][::5]-ymod[fs<2.5][::5])/yerr[fs<2.5][::5])**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<3.0 and 1.5<vs4<4.0 and 1.2<vsb<3.0 and 0.1<=d1<1.0 and 0.1<=d2<2.0 and 0.1<=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):
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta
vp1,vp2,vp3,vp4,vpb=vs2vp(np.array((vs1,vs2,vs3,vs4,vsb)))
lp = lnphysvalsprior(theta,delaydata,zb)
if np.isfinite(lp):
ldc = lnlike(theta, dispcurve, dispcurve_err,delaydata,zb)
lTT = lnTTprior(theta,delaydata)
lnA=lnamp(theta,delaydata)
return lp + ldc + lTT + lnA
else:
return -np.inf
def lnprob_step2(theta, dispcurve, dispcurve_err,delaydata,zb):
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = theta
vp1,vp2,vp3,vp4,vpb=vs2vp(np.array((vs1,vs2,vs3,vs4,vsb)))
lp = lnphysvalsprior(theta,delaydata,zb)
if np.isfinite(lp):
ldc = lnlike(theta, dispcurve, dispcurve_err,delaydata,zb)
lTT = lnTTprior(theta,delaydata)
lcv=constvfun(vs1,vs2,vs3,vs4,vsb)[0]
lnA=lnamp(theta,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
from multiprocessing import Pool
for zb in (2**np.linspace(-3,2,6)*0.2):
def runinv1_1step(i):
thetapeak=np.empty(9)
thetapeak[:]=np.nan
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)
sbbool=False
if delaydata[0]==delaydata[0]:
sbbool=True
cfm_use=curve_fillmean[startind:endind]
cfs_use=curve_fillstd[startind:endind]
cfm_use_masked = np.ma.masked_array(cfm_use, np.isnan(cfm_use))
cfs_use_masked = np.ma.masked_array(cfs_use, np.isnan(cfs_use))
dispcurve=np.average(cfm_use_masked, axis=0, weights=1/cfs_use**2)
dispcurve_err=np.sqrt(np.average(cfs_use_masked**2,axis=0,weights=1/cfs_use**2))
data = (dispcurve, dispcurve_err, delaydata, zb)
nwalkers=40
niter=6500
p0 = np.array([np.random.uniform(0.4,0.8,nwalkers*10),
np.random.uniform(0.4,1.8,nwalkers*10),
np.random.uniform(0.8,3.0,nwalkers*10),
np.random.uniform(1.5,4.0,nwalkers*10),
np.random.uniform(0.02,1.0,nwalkers*10),
np.random.uniform(0.02,2.0,nwalkers*10),
np.random.uniform(0.02,3.0,nwalkers*10),
np.random.uniform(1.2,3.0,nwalkers*10),
np.random.uniform(0.0,2.0,nwalkers*10)]).T
p0=p0[(p0[:,3]>p0[:,2])&(p0[:,2]>p0[:,1])&(p0[:,1]>p0[:,0])][:nwalkers]
ndim = 9
def main(p0,nwalkers,niter,ndim,lnprob,data):
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=data)
p0, _, _ = sampler.run_mcmc(p0, 1500)
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
thetas=samples.reshape(-1,40,9)[::500,:,:].reshape(-1,9)
probs=np.zeros(thetas.shape[0])
for j in range(thetas.shape[0]):
if lnphysvalsprior(thetas[j], delaydata,zb)==0:
probs[j]=np.exp(lnprob_step1(thetas[j], dispcurve, dispcurve_err, delaydata,zb))
else:
probs[j]=0
thetabest=samples[np.argmax(sampler.flatlnprobability)]
if False:
if sbbool:
thetapeak=findpeak(thetas,sbbool)
else:
thetapeak[:7]=findpeak(thetas,sbbool)
return dict(thetas=thetas,probs=probs,thetabest=thetabest,thetapeak=thetapeak,
sbbool=sbbool,dispcurve=dispcurve,dispcurve_err=dispcurve_err,delaydata=delaydata)
if __name__ == '__main__':
with Pool(20) as p:
allres=p.map(runinv1_1step, np.arange(60))
thetasall=np.zeros((60,520,9))
probsall=np.zeros((60,520))
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 ires,res in enumerate(allres):
thetasall[ires]=res['thetas']
probsall[ires]=res['probs']
thetabests[ires]=res['thetabest']
thetapeaks[ires]=res['thetapeak']
sbbools[ires]=res['sbbool']
dispcurves[ires]=res['dispcurve']
dispcurve_errs[ires]=res['dispcurve_err']
delaydata_all[ires]=res['delaydata']
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)
/tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] /tmp/ipykernel_2979128/1133263168.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1133263168.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1133263168.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
import matplotlib.pyplot as plt
import tqdm as tqdm
import numpy as np
import scipy.stats
import pickle
for zb in tqdm.tqdm((2**np.linspace(-3,2,6)*0.2)):
thetasall=np.load('data/MCMC/thetasall%s_r1.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r1.npy'%int(zb*1e3))
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.4:1.8:nbins*1j,0.8:3.0:nbins*1j, 1.5:4.0:nbins*1j]
vg = np.linspace(0.4, 0.8, nbins)
wg = np.linspace(0.4, 1.8, nbins)
xg = np.linspace(0.8, 3.0, nbins)
yg = np.linspace(1.5, 4.0, 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])
with open('data/MCMC/constvfun%s.pck'%int(zb*1e3), 'wb') as file_handle:
pickle.dump(constvfun, file_handle)
for zb in (2**np.linspace(-3,2,6)*0.2):
with open('data/MCMC/constvfun%s.pck'%int(zb*1e3), 'rb') as file_handle:
constvfun = pickle.load(file_handle)
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))
def runinv2_1step(i):
thetapeak=np.empty(9)
thetapeak[:]=np.nan
startind=4*i
midind=4*i+12
endind=4*i+25
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
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=40
niter=6500
p0 = np.array([np.random.uniform(0.4,0.8,nwalkers*10),
np.random.uniform(0.4,1.8,nwalkers*10),
np.random.uniform(0.8,3.0,nwalkers*10),
np.random.uniform(1.5,4.0,nwalkers*10),
np.random.uniform(0.02,1.0,nwalkers*10),
np.random.uniform(0.02,2.0,nwalkers*10),
np.random.uniform(0.02,3.0,nwalkers*10),
np.random.uniform(1.2,3.0,nwalkers*10),
np.random.uniform(0.0,2.0,nwalkers*10)]).T
p0=p0[(p0[:,3]>p0[:,2])&(p0[:,2]>p0[:,1])&(p0[:,1]>p0[:,0])][:nwalkers]
ndim = 9
def main(p0,nwalkers,niter,ndim,lnprob,data):
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=data)
p0, _, _ = sampler.run_mcmc(p0, 1500)
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)]
thetas=samples.reshape(-1,40,9)[::500,:,:].reshape(-1,9)
probs=np.zeros(thetas.shape[0])
for j in range(thetas.shape[0]):
if lnphysvalsprior(thetas[j], delaydata,zb)==0:
probs[j]=np.exp(lnprob_step2(thetas[j], dispcurve, dispcurve_err, delaydata,zb))
else:
probs[j]=0
return dict(thetas=thetas,probs=probs,thetabest=thetabest,
sbbool=sbbool,dispcurve=dispcurve,dispcurve_err=dispcurve_err,delaydata=delaydata)
if __name__ == '__main__':
with Pool(20) as p:
allres=p.map(runinv2_1step, np.arange(387))
thetasall=np.zeros((387,520,9))
probsall=np.zeros((387,520))
thetabests=np.zeros((387,9))
sbbools=np.zeros(387)
dispcurves=np.zeros((387,50))
dispcurve_errs=np.zeros((387,50))
delaydata_all=np.zeros((387,3))
for ires,res in enumerate(allres):
thetasall[ires]=res['thetas']
probsall[ires]=res['probs']
thetabests[ires]=res['thetabest']
sbbools[ires]=res['sbbool']
dispcurves[ires]=res['dispcurve']
dispcurve_errs[ires]=res['dispcurve_err']
delaydata_all[ires]=res['delaydata']
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)
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((387,50))
dcurve_obs_err=np.zeros((387,50))
for chan in range(387):
thetabest=thetabests[chan].copy()
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
startind=4*chan
midind=4*chan+12
endind=4*chan+25
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)
np.nanmean(probs)
/tmp/ipykernel_2979128/1740600676.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1740600676.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1740600676.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
probsums=np.zeros(6)
for zbi,zb in enumerate((2**np.linspace(-3,2,6)*0.2)):
probs=np.zeros(387)
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))
probsums[zbi]=np.mean(np.log10(np.mean(probsall,axis=1)))
probsums=np.zeros(6)
for zbi,zb in enumerate((2**np.linspace(-3,2,6)*0.2)):
probs=np.zeros(387)
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
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=np.nanmean(thetasall,axis=1)
#thetabests=np.load('data/MCMC/thetabests%s_r2.npy'%int(zb*1e3))
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((387,50))
dcurve_obs_err=np.zeros((387,50))
for chan in range(387):
thetabest=thetabests[chan].copy()
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
startind=4*chan
midind=4*chan+12
endind=4*chan+25
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_step2(thetabest, dcurve_obs[chan], dcurve_obs_err[chan],delaydata,zb)
probsums[zbi]=np.nanmean(probs)
/tmp/ipykernel_2979128/2727449843.py:21: RuntimeWarning: Mean of empty slice thetabests=np.nanmean(thetasall,axis=1) /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /tmp/ipykernel_2979128/2727449843.py:40: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/2727449843.py:41: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/2727449843.py:42: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
plt.plot((2**np.linspace(-3,2,6)*0.2),np.exp(probsums)/np.max(np.exp(probsums)))
plt.xscale('log')
plt.xticks(ticks=np.array((0.05,0.1,0.2,0.4)),labels=np.array((0.05,0.1,0.2,0.4)),fontsize=40)
plt.yticks(fontsize=40)
plt.ylim(0,1.1)
plt.ylabel('average model\n probability',fontsize=40)
plt.xlabel('slumped block thickness',fontsize=40)
plt.tight_layout()
plt.savefig('figs/modprobs.png')
zb=(2**np.linspace(-3,2,6)*0.2)[3]
#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))
thetameds=np.nanmedian(thetasall,axis=1)
thetameds[:,4:7][thetameds[:,4:7]<0.1]=0.1
thetameds[:,4][thetameds[:,4]>1]=1
thetameds[:,5][thetameds[:,5]>2]=2
thetameds[:,6][thetameds[:,6]>3]=3
thetameds[(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
thetameds_smooth=np.apply_along_axis(moving_average,0,thetameds)
#find array channels associated with model columns and smooth result
midinds_fullarr=(4*(np.arange(0,387))+12)*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((373,250))
velmodmatr=np.zeros((387,250))
dep_ax=np.linspace(0,5,250)
#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((373,250))
modstdmatr=np.zeros((387,250))
for chan in range(387):
thetabest=thetameds[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(373):
thetabest=thetameds_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,5,250)
dcurve_mod=np.zeros((387,50))
t_residuals=np.repeat(np.nan,387)
t_obs=np.repeat(np.nan,387)
t_mod=np.repeat(np.nan,387)
amp_obs=np.repeat(np.nan,387)
amp_err=np.repeat(np.nan,387)
amp_mod=np.repeat(np.nan,387)
ps=np.repeat(np.nan,387)
dcurve_obs=np.zeros((387,50))
dcurve_obs_err=np.zeros((387,50))
t_err=np.repeat(np.nan,387)
probs=np.repeat(np.nan,387)
f1, f2 = 0.25, 4.0
periods=1/np.logspace(np.log10(f1),np.log10(f2),50)
for chan in range(387):
thetabest=thetameds[chan].copy()
vs1,vs2,vs3,vs4,d1,d2,d3,vsb,db = thetabest
ele=elevations[chan]
startind=4*chan
midind=4*chan+12
endind=4*chan+25
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_2979128/1793960487.py:19: RuntimeWarning: All-NaN slice encountered thetameds=np.nanmedian(thetasall,axis=1) /tmp/ipykernel_2979128/1793960487.py:32: RuntimeWarning: invalid value encountered in scalar divide return np.nansum(a*w)/((~np.isnan(a))*w).sum() /tmp/ipykernel_2979128/1793960487.py:32: RuntimeWarning: invalid value encountered in scalar divide return np.nansum(a*w)/((~np.isnan(a))*w).sum() /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /tmp/ipykernel_2979128/1793960487.py:157: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/1793960487.py:158: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2979128/1793960487.py:159: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2979128/3684714611.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,373)
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,373)
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((373,50))
ps_smooth=np.apply_along_axis(moving_average,0,ps)
probs_smooth=np.repeat(np.nan,373)
for chan in range(373):
sbt=sbts_smooth[chan]
sbbool=sbt>0
thetabest=thetameds_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_2979128/1793960487.py:32: RuntimeWarning: invalid value encountered in scalar divide return np.nansum(a*w)/((~np.isnan(a))*w).sum() /tmp/ipykernel_2979128/3684714611.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
-158.28853703703786
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'),loc='lower right')
for axi in range(4):
ax=list((ax1,ax2,ax3,ax4))[axi]
axlab=list(('c)','b)','d)','a)'))[axi]
ax.annotate(
axlab,
xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
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,label='Observed')
ax5.plot(invdist,amp_mod,color='k',linewidth=6,label='Predicted')
ax5.fill_between(invdist,amp_obs-amp_err,amp_obs+amp_err,alpha=0.25)
plt.legend()
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)')
ax4.legend(('Observed','Modeled'),loc='lower right')
plt.tight_layout()
for axi in range(4):
ax=list((ax1,ax2,ax3,ax4))[axi]
axlab=list(('c)','b)','d)','a)'))[axi]
ax.annotate(
axlab,
xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
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)
thetameds=np.nanmedian(thetasall,axis=1)
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((4,0))
ax1.set_aspect(2)
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((4,0))
ax2.set_aspect(2)
fig.colorbar(velstd,label='Velocity s.d./\n Velocity (-)')
ax2.annotate(
'c)', xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
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,520).reshape(-1,520,1),thetasall[:,:,4:7]),axis=2)
deps4plt_cum=np.cumsum(deps4plt,axis=2)
deps_plt=np.cumsum(np.concatenate(([elevations_smooth],thetameds_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=(thetameds_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((4,0))
ax3.set_ylabel('depth (km)')
ax3.set_xlabel('Along-fiber Distance (km)')
plt.legend()
ax3.set_aspect(2)
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((4,0))
fig.colorbar(velmodr,label='Velocity (km/s)')
ax4.annotate(
'a)', xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
ax4.set_aspect(2)
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((4,0))
fig.colorbar(velstd,label='Velocity s.d./\n Velocity (-)')
ax5.annotate(
'b)', xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
ax5.set_aspect(2)
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,520).reshape(-1,520,1),thetasall[:,:,4:7]),axis=2)
deps4plt_cum=np.cumsum(deps4plt,axis=2)
deps_plt=np.cumsum(np.concatenate(([elevations],thetameds[:,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=(thetameds[:,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((4,0))
ax6.set_ylabel('depth (km)')
ax6.set_xlabel('Along-fiber Distance (km)')
ax6.set_aspect(2)
plt.savefig('figs/dstdr.png')
/tmp/ipykernel_2979128/1250971831.py:1: RuntimeWarning: All-NaN slice encountered thetameds=np.nanmedian(thetasall,axis=1) /tmp/ipykernel_2979128/1793960487.py:32: RuntimeWarning: invalid value encountered in scalar divide 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=thetameds[:,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<3:
axes[i,j].set_ylabel('layer %s thickness (km)'%(plot_i+1))
else:
axes[i,j].legend(loc='lower left')
axes[i,j].set_ylabel('slumped block depth (km)')
axlab=list(('a)','b)','c)','d)'))[plot_i]
axes[i,j].annotate(
axlab,
xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
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<3:
axes[i,j].set_ylabel('layer %s thickness (km)'%(plot_i+1))
else:
axes[i,j].legend(loc='lower left')
axes[i,j].set_ylabel('slumped block depth (km)')
axlab=list(('a)','b)','c)','d)'))[plot_i]
axes[i,j].annotate(
axlab,
xy=(0, 1), xycoords='axes fraction',
xytext=(+0.5, -0.5), textcoords='offset fontsize',
fontsize='medium', verticalalignment='top', fontfamily='serif',
bbox=dict(facecolor='0.7', edgecolor='none', pad=3.0))
plt.savefig('figs/dlayerr.png')
/home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1563: RuntimeWarning: All-NaN slice encountered return function_base._ureduce(a, /tmp/ipykernel_2979128/1793960487.py:32: RuntimeWarning: invalid value encountered in scalar divide return np.nansum(a*w)/((~np.isnan(a))*w).sum()
len(sbbools)
387
onedmod=(np.nanmean(thetameds,axis=0)*1e3).round(-1).astype('int')
onedmoderr=(np.nanstd(thetameds,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(thetameds[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(thetasall[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)
thetabests
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])
str_res
'Model residuals show the greatest reduction for a slumped block of 200 meter thickness. The block takes on a velocity of 2370 m/s with a standard deviation of 40 m/s. In order of increasing depth, the remaining model layers take on mean velocities of 630, 1270, 2350, and 3350 m/s, with respective standard deviations of 40, 170, 120, and 160 m/s. As the array crosses onto the proposed location of the NBB, the inversion indicates the top of the slumped block is located $560\\pm110$ meters beneath the ground surface. 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. 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 $930\\pm220$. '
import corner
corner.corner(thetasall[:,:,:7].reshape(-1,7))