Note
Go to the end to download the full example code.
Dual data types (low and high moment seperately)¶
try:
# Check if the code is running in an IPython kernel (which includes Jupyter notebooks)
get_ipython()
# If the above line doesn't raise an error, it means we are in a Jupyter environment
# Execute the magic commands using IPython's run_line_magic function
get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')
except:
# If get_ipython() raises an error, we are not in a Jupyter environment
# # # #%load_ext autoreload
# # # #%autoreload 2
pass
import integrate as ig
import h5py
import numpy as np
import matplotlib.pyplot as plt
hardcopy=True
case = 'DAUGAARD'
files = ig.get_case_data(case=case)
f_data_h5 = files[0]
f_data_h5 = 'DAUGAARD_AVG.h5'
file_gex= ig.get_gex_file_from_data(f_data_h5)
print("Using data file: %s" % f_data_h5)
print("Using GEX file: %s" % file_gex)
1. Setup the prior model ($rho(mathbf{m},mathbf{d})$¶
In this example a simple layered prior model will be considered
1a. first, a sample of the prior model parameters, $rho(mathbf{m})$, will be generated¶
N=25000
# Layered model
f_prior_h5 = ig.prior_model_layered(N=N,lay_dist='chi2', NLAY_deg=3, RHO_min=1, RHO_max=3000)
f_prior_data_h5 = ig.prior_data_em(f_prior_h5, file_gex, showInfo=0)
ig.integrate_update_prior_attributes(f_prior_data_h5)
into TWO data sets, LOW and HIGH moment.
ig.plot_data_prior(f_prior_data_h5, f_data_h5)
Read D_obs from f_data_h5
with h5py.File(f_data_h5, 'r') as f:
D_obs = f['D1/d_obs'][:]
D_std = f['D1/d_std'][:]
# Alternatively, use the ig.load_data function
#D_obs = ig.load_data(f_data_h5, id=1, showInfo=1)['d_obs'][0]
#D_std = ig.load_data(f_data_h5, id=1, showInfo=1)['d_std'][0]
with h5py.File(f_prior_data_h5, 'r') as f:
D = f['/D1'][:]
# Alternatively, use the ig.load_prior_data function
#D = ig.load_prior_data(f_prior_data_h5)[0][0]
# Now splot the into low and high moment data sets
# The low moment data set will be the first 14 columns, and the high moment data set will be the last columns.
nd = D_obs.shape[1]
n_low = 14
n_high = nd - n_low
# set i low to 0:n_low-1
i_low = range(n_low)
i_high = range(n_low, nd)
# Split prior data
D_low = D[:,i_low]
D_high = D[:,i_high]
# Split observed data
D_obs_low = D_obs[:,i_low]
D_std_low = D_std[:,i_low]*2
D_obs_high = D_obs[:,i_high]
D_std_high = D_std[:,i_high]*2
plt.semilogy(D_obs[0],'k.',markersize=10, label='D_obs')
plt.semilogy(D_obs_low[0],'.', color='gray', markersize=6, label='D_obs_low')
plt.semilogy(D_obs_high[0],'.', color='gray',markersize=4, label='D_obs_high')
plt.semilogy(D[0],'b.',markersize=10, label='D_obs')
plt.semilogy(D_low[0],'.', color='yellow', markersize=6, label='D_obs_low')
plt.semilogy(D_high[0],'.', color='yellow', markersize=6, label='D_obs_low')
plt.legend()
f_data_dual_h5 = 'DAUGAARD_AVG_dual.h5'
useOldMethod = False
if useOldMethod:
ig.copy_hdf5_file(f_data_h5,f_data_dual_h5)
# Delete D1
with h5py.File(f_data_dual_h5, 'a') as f:
# show groups in f['']
if 'D1' in f.keys():
del(f['D1'])
# Update D1 and D2
with h5py.File(f_data_dual_h5, 'a') as f:
# remove 'D1'
#del f['D1']
f.create_dataset('D1/d_obs', data=D_obs_low)
f.create_dataset('D1/d_std', data=D_std_low)
f['D1'].attrs['noise_model'] = 'gaussian'
f.create_dataset('D2/d_obs', data=D_obs_high)
f.create_dataset('D2/d_std', data=D_std_high)
f['D2'].attrs['noise_model'] = 'gaussian'
else:
# Alternatively, use the ig.save_data_gaussian function
ig.copy_hdf5_file(f_data_h5,f_data_dual_h5)
ig.save_data_gaussian(D_obs_low, D_std = D_std_low, f_data_h5 = f_data_dual_h5, id=1, showInfo=0)
ig.save_data_gaussian(D_obs_high, D_std = D_std_high, f_data_h5 = f_data_dual_h5, id=2, showInfo=0)
f_prior_data_dual_h5 = 'PRIOR_dual.h5'
ig.copy_hdf5_file(f_prior_data_h5,f_prior_data_dual_h5)
if useOldMethod:
with h5py.File(f_prior_data_dual_h5, 'a') as f:
# show groups in f['']
if 'D1' in f.keys():
del(f['D1'])
f.create_dataset('D1', data=D_low)
f.create_dataset('D2', data=D_high)
else:
# Alternatively, use the ig.save_data_gaussian function
ig.save_prior_data(f_prior_data_dual_h5, D_low, id=1, force_delete=True)
ig.save_prior_data(f_prior_data_dual_h5, D_high, id=2, force_delete=False)
ig.plot_data_prior(f_prior_data_dual_h5, f_data_dual_h5, id=1)
ig.plot_data_prior(f_prior_data_dual_h5, f_data_dual_h5, id=2)
Sample the posterior $sigma(mathbf{m})$¶
The posterior distribution is sampling using the extended rejection sampler.
N_use = 100000 #%N
N_cpu = 8
f_post_arr = []
updatePostStat=False
showInfo = 1
import time
t_inversion = []
for itype in [0,1,2,3]:
t_start = time.time()
if itype == 0:
# LOW AND HIGH MOMENT AS ONE DATA SET - THE ORIGINAL METHOD
f_post_h5 = ig.integrate_rejection(f_prior_data_h5,
f_data_h5,
f_post_h5='POST_type%d.h5' % itype,
N_use = N_use, Ncpu=N_cpu,
showInfo=showInfo,
updatePostStat=updatePostStat)
elif itype == 1:
# LOW MOMENT ONLY
f_post_h5 = ig.integrate_rejection(f_prior_data_dual_h5,
f_data_dual_h5,
f_post_h5='POST_type%d.h5' % itype,
N_use = N_use, Ncpu=N_cpu,
showInfo=showInfo,
updatePostStat=updatePostStat,
id_use = [1])
elif itype == 2:
# HIGH MOMENT ONLY
f_post_h5 = ig.integrate_rejection(f_prior_data_dual_h5,
f_data_dual_h5,
# f_post_h5='POST_type%d.h5' % itype,
N_use = N_use, Ncpu=N_cpu,
showInfo=showInfo,
updatePostStat=updatePostStat,
id_use = [2])
elif itype == 3:
# JOINT INVERSION USING BOTH LOW AND HIGH MOMENT
f_post_h5 = ig.integrate_rejection(f_prior_data_dual_h5,
f_data_dual_h5,
f_post_h5='POST_type%d.h5' % itype,
N_use = N_use, Ncpu=N_cpu,
showInfo=showInfo,
updatePostStat=updatePostStat,
id_use = [1,2])
t_inversion.append(time.time() - t_start)
f_post_arr.append(f_post_h5)
for f_post_h5 in f_post_arr:
ig.integrate_posterior_stats(f_post_h5)
for f_post_h5 in f_post_arr:
# % Plot Profiles
ig.plot_profile(f_post_h5, i1=1000, i2=2000, im=1, hardcopy=hardcopy)
for f_post_h5 in f_post_arr:
ig.plot_T_EV(f_post_h5, pl='T', hardcopy=hardcopy)
plt.show()
#ig.plot_T_EV(f_post_h5, pl='EV', hardcopy=hardcopy)
for f_post_h5 in f_post_arr:
# % Plot prior, posterior, and observed data
ig.plot_data_prior_post(f_post_h5, i_plot=100, hardcopy=hardcopy)
#ig.plot_data_prior_post(f_post_h5, i_plot=0, hardcopy=hardcopy)
with h5py.File(f_data_h5, 'r') as f:
nobs = f['D1/d_obs'].shape[0]
ni=len(f_post_arr)
T = np.zeros((ni,nobs))
EV = np.zeros((ni,nobs))
EV_post = np.zeros((ni,nobs))
N_UNIQUE = np.zeros((ni,nobs))
i=-1
for f_post_h5 in f_post_arr:
with h5py.File(f_post_h5, 'r') as f:
i=i+1
T[i] = f['/T'][:]
EV[i] = f['/EV'][:]
N_UNIQUE[i] = f['/N_UNIQUE'][:]
EV_post[i] = f['/EV_post'][:]
fig, ax = plt.subplots(4,1,figsize=(10,4))
ax[0].semilogy(T.T,'-', linewidth=1)
ax[0].set_ylabel('T')
ax[1].plot(EV.T,'-', linewidth=1)
ax[1].set_ylabel('EV')
ax[2].plot(EV_post.T,'-', linewidth=1)
ax[2].set_ylabel('EV_post')
ax[3].plot(N_UNIQUE.T,'-', linewidth=1)
ax[3].set_ylabel('N_UNIQUE')
fig=plt.figure(figsize=(6,6))
plt.semilogx(T,EV,'.')
plt.xlabel('T')
plt.ylabel('EV')
fig=plt.figure(figsize=(6,6))
plt.semilogx(T,N_UNIQUE,'.')
plt.xlabel('T')
plt.ylabel('N_UNIQUE')
ig.post_to_csv(f_post_h5)