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('/kuafu/users/elibird/BishopAmbientNoise/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='/kuafu/users/elibird/BishopAmbientNoise/DispCurves/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('/kuafu/scratch/elibird/MammothSouth/CCOwensRiver/DispCurvesAlt/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='/kuafu/scratch/elibird/MammothSouth/CCOwensRiver/DispCurvesAlt/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('/home/elibird/200km/notebooks/xcor/BishopFinal/data/curve_fillmean.npy',curve_fillmean)
np.save('/home/elibird/200km/notebooks/xcor/BishopFinal/data/curve_fillstd.npy',curve_fillstd)
curve_fillmean=np.load('/home/elibird/200km/notebooks/xcor/BishopFinal/data/curve_fillmean.npy')
curve_fillstd=np.load('/home/elibird/200km/notebooks/xcor/BishopFinal/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')
def vs2vp(x):
return 0.9409+2.0947*x-0.8206*x**2+0.2683*x**3-0.0251*x**4
vs2vp(1.2)
2.68445104
# 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)
return lp + ldc + lTT
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]
if lcv==lcv:
return lp + ldc + lTT + lcv
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 [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_noamp.npy'%int(zb*1e3),thetapeaks)
np.save('data/MCMC/thetasall%s_r1_noamp.npy'%int(zb*1e3),thetasall)
np.save('data/MCMC/thetabests%s_r1_noamp.npy'%int(zb*1e3),thetabests)
np.save('data/MCMC/probsall%s_r1_noamp.npy'%int(zb*1e3),probsall)
np.save('data/MCMC/sbbools%s_r1_noamp.npy'%int(zb*1e3),sbbools)
/tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:14: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/943124656.py:12: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /tmp/ipykernel_2177305/943124656.py:13: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/943124656.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/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /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_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.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 [0.2]:
thetasall=np.load('data/MCMC/thetasall%s_r1_noamp.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r1_noamp.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_noamp.pck'%int(zb*1e3), 'wb') as file_handle:
pickle.dump(constvfun, file_handle)
for zb in [0.2]:
with open('data/MCMC/constvfun%s_noamp.pck'%int(zb*1e3), 'rb') as file_handle:
constvfun = pickle.load(file_handle)
thetasall=np.load('data/MCMC/thetasall%s_r1_noamp.npy'%int(zb*1e3))
probsall=np.load('data/MCMC/probsall%s_r1_noamp.npy'%int(zb*1e3))
thetabests=np.load('data/MCMC/thetabests%s_r1_noamp.npy'%int(zb*1e3))
sbbools=np.load('data/MCMC/sbbools%s_r1_noamp.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
if True:
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(runinv2_1step, np.arange(387))
thetasall=np.zeros((387,520,9))
probsall=np.zeros((387,520))
thetabests=np.zeros((387,9))
thetapeaks=np.empty((387,9))
thetapeaks[:]=np.nan
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']
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_r2_noamp.npy'%int(zb*1e3),thetapeaks)
np.save('data/MCMC/thetasall%s_r2_noamp.npy'%int(zb*1e3),thetasall)
np.save('data/MCMC/thetabests%s_r2_noamp.npy'%int(zb*1e3),thetabests)
np.save('data/MCMC/probsall%s_r2_noamp.npy'%int(zb*1e3),probsall)
np.save('data/MCMC/sbbools%s_r2_noamp.npy'%int(zb*1e3),sbbools)
/tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:18: 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/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:18: 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/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /home/elibird/miniconda3/envs/bird/lib/python3.9/site-packages/pysurf96/wrapper.py:99: RuntimeWarning: overflow encountered in cast error = surfdisp96( /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_2177305/526877977.py:16: RuntimeWarning: Mean of empty slice delaydata[0]=np.nanmean(meanpick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/526877977.py:17: RuntimeWarning: Mean of empty slice delaydata[1]=np.sqrt(np.nanmean(varpick[((startind*2)-25):((endind*2)-25)])) /tmp/ipykernel_2177305/526877977.py:18: RuntimeWarning: Mean of empty slice delaydata[2]=np.nanmean(ppick[((startind*2)-25):((endind*2)-25)]) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps) /tmp/ipykernel_2177305/3014339761.py:25: RuntimeWarning: invalid value encountered in arcsin i_s=np.arcsin(alphas*ps)
blockv_noamp=np.load('data/MCMC/thetasall%s_r2_noamp.npy'%int(zb*1e3))[:,:,7][sbbools.astype('bool')].flatten()
blockv_amp=np.load('data/MCMC/thetasall%s_r2.npy'%int(zb*1e3))[:,:,7][sbbools.astype('bool')].flatten()
plt.rcParams.update({'font.size': 30})
plt.hist(blockv_noamp,bins=100,alpha=0.5,density=True,label='Without Amplitude Information')
plt.hist(blockv_amp,bins=100,alpha=0.5,density=True,label='With Amplitude Information')
plt.xlim(1.2,3.0)
plt.ylabel('Density')
plt.xlabel('NBB Velocity (km/s)')
plt.legend(loc='upper left')
plt.tight_layout()
plt.savefig('figs/ampdist.png')