import sys
sys.path.append("../")
import gempy as gp
%matplotlib inline
import pymc
import numpy as np
import math
geo_data = gp.read_pickle("../Tutorial/BasicFault.pickle")
geo_data.n_faults = 1
geo_data.interfaces.head()
X | Y | Z | formation | series | order_series | isFault | |
---|---|---|---|---|---|---|---|
0 | 800.0 | 1000.0 | -1600.0 | MainFault | fault | 1 | True |
1 | 1200.0 | 1000.0 | -400.0 | MainFault | fault | 1 | True |
2 | 1100.0 | 1000.0 | -700.0 | MainFault | fault | 1 | True |
3 | 900.0 | 1000.0 | -1300.0 | MainFault | fault | 1 | True |
4 | 1000.0 | 1000.0 | -1000.0 | MainFault | fault | 1 | True |
geo_data.geo_data_type = 'whatever'
gp.plot_data(geo_data, direction='y')
<gempy.Visualization.PlotData at 0x7f5ab768a6d8>
# stdev for x,y,z interface coordinates
geo_data.interfaces['X_std'] = None
geo_data.interfaces['Y_std'] = None
geo_data.interfaces['Z_std'] = 100
geo_data.interfaces['dist_type'] = "Normal" # should be initially None
geo_data.interfaces['X_dist'] = None
geo_data.interfaces['Y_dist'] = None
geo_data.interfaces['Z_dist'] = None
#geo_data.foliations['X_std'] = None
#geo_data.foliations['Y_std'] = 0
#geo_data.foliations['Z_std'] = 0
geo_data.interfaces["Z_std"] = 500
geo_data.interfaces.head()
X | Y | Z | formation | series | order_series | isFault | X_std | Y_std | Z_std | dist_type | X_dist | Y_dist | Z_dist | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 800.0 | 1000.0 | -1600.0 | MainFault | fault | 1 | True | None | None | 500 | Normal | None | None | None |
1 | 1200.0 | 1000.0 | -400.0 | MainFault | fault | 1 | True | None | None | 500 | Normal | None | None | None |
2 | 1100.0 | 1000.0 | -700.0 | MainFault | fault | 1 | True | None | None | 500 | Normal | None | None | None |
3 | 900.0 | 1000.0 | -1300.0 | MainFault | fault | 1 | True | None | None | 500 | Normal | None | None | None |
4 | 1000.0 | 1000.0 | -1000.0 | MainFault | fault | 1 | True | None | None | 500 | Normal | None | None | None |
interp_data = gp.InterpolatorInput(geo_data,
u_grade = [3, 3],
compile_theano=True)
float32 [2, 2]
sol = gp.compute_model(interp_data)
[3, 3]
adding dataframe columns for uncertainty analysis
Z_rest = pymc.Normal('Z_unc_rest', interp_data.interpolator.pandas_rest_layer_points['Z'].as_matrix().astype('float32'),
1./interp_data.interpolator.pandas_rest_layer_points['Z_std'].as_matrix().astype('float32'))
Z_ref = pymc.Normal('Z_unc_ref', interp_data.interpolator.pandas_ref_layer_points['Z'].as_matrix().astype('float32'),
1./interp_data.interpolator.pandas_ref_layer_points['Z_std'].as_matrix().astype('float32'))
Z_rest.get_stoch_value()
array([[ 1.46208002, -0.16228111, 0.15947797, 0.27908407, -0.06837028, 0.11290613, 1.56334182, 1.21185192, 0.67252936, 0.33573741, 0.24979466, 0.20912485, 0.46076347, 0.64697701, 0.7478793 , 0.25877077, 0.31715363, 0.83621002, 0.38905729, 0.76213303, 0.91305479, 0.58656948, 1.35064574, 0.62597685, 0.97317111, 0.2814209 , 0.58052269, 0.34527976, 0.44034878, 0.5655428 , -0.11403995, 0.4879321 , 0.38160715, 0.54497289]])
i = interp_data.get_input_data()
@pymc.deterministic
def gempy_model(value=0, input_ = i, Z_rest_m = Z_rest, Z_ref_m = Z_ref):
#xp = X_priors, yp = Y_priors, zp = Z_priors, ref = data_ref):
Z_ref_m_rep = np.apply_over_axes(lambda x,y:
np.repeat(x, interp_data.interpolator.len_interfaces - 1),
Z_ref_m, axes = [0])
input_[4][:, 2] = Z_ref_m_rep
input_[5][:, 2] = Z_rest_m
sol, pot = interp_data.th_fn(*input_)
print(sol[0,0,:])
# return solution
return sol
[3, 3] [ 0. 0. 0. ..., 0. 0. 0.]
# set number of iterations
iterations = 10
# set model
pymc_model = pymc.Model(Z_ref, Z_rest, gempy_model)
MCFS_RUN = pymc.MCMC(pymc_model)
MCFS_RUN.sample(iter=iterations)
[ 2. 2. 2. ..., 2. 2. 2.] [------- 20% ] 2 of 10 complete in 2.3 sec[ 0. 0. 0. ..., 2. 2. 2.] [----------- 30% ] 3 of 10 complete in 4.6 sec[ 2. 2. 2. ..., 4. 4. 4.] [--------------- 40% ] 4 of 10 complete in 6.8 sec[ 2. 2. 2. ..., 5. 5. 5.] [-----------------50% ] 5 of 10 complete in 9.0 sec
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
def p(i):
ax.imshow(MCFS_RUN.trace("gempy_model", chain=-1)[i][0, 0, :].reshape(50,50,50)[:,24,:].T, origin="lower")
fig, ax = plt.subplots()
interact(p, i=widgets.IntSlider(min=0, max=9, step=1, value=0));
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-29-373a0f7b5c39> in <module>() 1 fig, ax = plt.subplots() ----> 2 interact(p, i=widgets.IntSlider(min=0, max=9, step=1, value=0)); NameError: name 'interact' is not defined
plt.imshow(MCFS_RUN.trace("gempy_model", chain=-1)[3][0, 0, :].reshape(50,50,50)[:,24,:].T, origin="lower")
<matplotlib.image.AxesImage at 0x7fc37822bc18>