# These two lines are necessary only if gempy is not installed
import sys, os
sys.path.append("../")
# Importing gempy
import gempy as gp
# Embedding matplotlib figures into the notebooks
%matplotlib inline
# Aux imports
import numpy as np
import matplotlib.pyplot as plt
geo_data = gp.create_data([696000-10000,747000 + 20600,6863000 - 20600,6950000 + 20600,-20000, 600],[16, 16, 16],
path_f = os.pardir+"/input_data/a_Foliations.csv",
path_i = os.pardir+"/input_data/a_Points.csv")
# Assigning series to formations as well as their order (timewise)
gp.set_data_series(geo_data, {"EarlyGranite_Series": 'EarlyGranite',
"BIF_Series":('SimpleMafic2', 'SimpleBIF'),
"SimpleMafic_Series":'SimpleMafic1'},
order_series = ["EarlyGranite_Series",
"BIF_Series",
"SimpleMafic_Series"], verbose=1)
geo_data.geo_data_type = 'none'
geo_data.interfaces['X_std'] = None
geo_data.interfaces['Y_std'] = 0
geo_data.interfaces['Z_std'] = 50
geo_data.foliations['X_std'] = None
geo_data.foliations['Y_std'] = 0
geo_data.foliations['Z_std'] = 0
geo_data.foliations['dip_std'] = 10
geo_data.foliations['azimuth_std'] = 10
geo_data.foliations.head()
X | Y | Z | azimuth | dip | polarity | formation | series | order_series | G_x | G_y | G_z | X_std | Y_std | Z_std | dip_std | azimuth_std | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 735082.0630 | 6879102.25 | 480.551436 | 276.153239 | 80.0 | 1 | EarlyGranite | EarlyGranite_Series | 1 | -0.979134 | 0.105560 | 0.173648 | None | 0 | 0 | 10 | 10 |
1 | 735326.7815 | 6875283.75 | 452.966456 | 114.812185 | 80.0 | 1 | EarlyGranite | EarlyGranite_Series | 1 | 0.893898 | -0.413270 | 0.173648 | None | 0 | 0 | 10 | 10 |
2 | 738704.4065 | 6877046.00 | 458.148740 | 170.056246 | 80.0 | 1 | EarlyGranite | EarlyGranite_Series | 1 | 0.170058 | -0.970014 | 0.173648 | None | 0 | 0 | 10 | 10 |
3 | 743403.6875 | 6877878.25 | 454.749952 | 169.802194 | 80.0 | 1 | EarlyGranite | EarlyGranite_Series | 1 | 0.174357 | -0.969250 | 0.173648 | None | 0 | 0 | 10 | 10 |
4 | 734151.9690 | 6874304.75 | 448.179887 | 206.518042 | 80.0 | 1 | EarlyGranite | EarlyGranite_Series | 1 | -0.439697 | -0.881200 | 0.173648 | None | 0 | 0 | 10 | 10 |
interp_data = gp.InterpolatorInput(geo_data, compile_theano=True, u_grade=[3,3,3])
I am here I am in the setting float32 I am here [2, 2]
interp_data.data.interfaces.head()
X | Y | Z | formation | series | order_series | X_std | Y_std | Z_std | formation number | |
---|---|---|---|---|---|---|---|---|---|---|
0 | 0.251624 | 0.2501 | 0.545482 | EarlyGranite | EarlyGranite_Series | 1 | NaN | 0.0 | 0.000515 | 1 |
1 | 0.280857 | 0.288407 | 0.545625 | EarlyGranite | EarlyGranite_Series | 1 | NaN | 0.0 | 0.000515 | 1 |
2 | 0.49254 | 0.399287 | 0.545915 | EarlyGranite | EarlyGranite_Series | 1 | NaN | 0.0 | 0.000515 | 1 |
3 | 0.584269 | 0.375096 | 0.545618 | EarlyGranite | EarlyGranite_Series | 1 | NaN | 0.0 | 0.000515 | 1 |
4 | 0.582253 | 0.403318 | 0.545733 | EarlyGranite | EarlyGranite_Series | 1 | NaN | 0.0 | 0.000515 | 1 |
sol = gp.compute_model(interp_data)
[3, 3, 3]
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 3.1069411932094226e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
gp.plot_section(geo_data, sol[-1,0,:], 7)
<gempy.Visualization.PlotData at 0x7fb40c0e18d0>
def precompute_grav(interp_data):
"""
model_res: number of voxels
interp_data: interp_data object with the specific resolution
"""
import sys
import gempy as gp
sys.path.append("../gempy")
import GeoPhysics
import theano.tensor as T
import theano
from importlib import reload
reload(GeoPhysics)
#inter_data = gp.InterpolatorInput(geo_data, compile_theano=False)
gpp = GeoPhysics.GeoPhysicsPreprocessing(interp_data,580, [7.050000e+05,747000,6863000,6925000,-20000, 200],
res_grav = [5, 5],
n_cells = 10, mode='range')
#res_grav = [125, 85],
#n_cells =1000)
print(gpp)
# Compile Theano function
x_1 = T.matrix()
x_2 = T.matrix()
sqd = T.sqrt(T.maximum(
(x_1 ** 2).sum(1).reshape((x_1.shape[0], 1)) +
(x_2 ** 2).sum(1).reshape((1, x_2.shape[0])) -
2 * x_1.dot(x_2.T), 0
))
eu = theano.function([x_1, x_2], sqd, allow_input_downcast=True)
# Init
i_0 = 0
model_res = interp_data.resolution[0] * interp_data.resolution[1] * interp_data.resolution[2]
b_all = np.zeros((0,model_res), dtype=bool)
# 25 is the size of the chunk
for i_1 in np.arange(25, gpp.airborne_plane.shape[0]+1, 25, dtype=int):#np.linspace(gpp.airborne_plane.shape[0]/25, gpp.airborne_plane.shape[0], 26, dtype=int):
airborne_plane_s = gpp.airborne_plane[i_0:i_1]
d = eu(interp_data.data.grid.grid.astype('float'), airborne_plane_s)
ab_g = interp_data.data.grid.grid[np.argmin(d, axis=0)]
ab_g[:,2] = ab_g[:,2] + 0.01
d2 = eu(ab_g, interp_data.data.grid.grid.astype('float'))
# Max range to select voxels
range_ = (interp_data.data.grid.grid[:, 2].max() - interp_data.data.grid.grid[:, 2].min()) *0.9
# Boolean of the selection
b = d2<range_
# Beoolean for all measurements
b_all = np.vstack((b_all, b))
# Selection of the x y z coordenates for each measurement and computed distance in those coordinates
s_gr_x = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, :25].T[b].reshape(25,-1) - ab_g[:,0].reshape(25,-1)).astype('float')
s_gr_y = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, 25:50].T[b].reshape(25,-1) - ab_g[:,1].reshape(25,-1)).astype('float')
s_gr_z = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, -25:].T[b].reshape(25,-1) - ab_g[:,2].reshape(25,-1)).astype('float')
# getting the coordinates of the corners of the voxel
x_cor = np.stack((s_gr_x - gpp.vox_size[0], s_gr_x + gpp.vox_size[0]), axis=2)
y_cor = np.stack((s_gr_y - gpp.vox_size[1], s_gr_y + gpp.vox_size[1]), axis=2)
z_cor = np.stack((s_gr_z - gpp.vox_size[2], s_gr_z + gpp.vox_size[2]), axis=2)
# ...and prepare them for a vectorial op
x_matrix = np.repeat(x_cor, 4, axis=2)
y_matrix = np.tile(np.repeat(y_cor, 2, axis=2), (1, 1, 2))
z_matrix = np.tile(z_cor, (1, 1, 4))
# Distances to each corner of the voxel
s_r = np.sqrt(x_matrix**2 + y_matrix**2 + z_matrix**2)
# Computing z component
from scipy.constants import G
mu = np.array([1, -1, -1, 1, -1, 1, 1, -1])
tz = np.sum(- G * mu * (
x_matrix * np.log(y_matrix + s_r) +
y_matrix * np.log(x_matrix + s_r) -
z_matrix * np.arctan(x_matrix * y_matrix /
(z_matrix * s_r))), axis=2)
# Stacking the precomputation
if i_0 == 0:
tz_all = tz
else:
tz_all = np.vstack((tz_all, tz))
i_0 = i_1
return tz_all, b_all
tz_, b_ = precompute_grav(interp_data)
0.1908919617392378 <GeoPhysics.GeoPhysicsPreprocessing object at 0x7fb474bef710>
import pymc
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'))
dip_u = pymc.Normal('dip_unc', interp_data.data.foliations['dip'].as_matrix().astype('float32'),
1./interp_data.data.foliations['dip_std'].as_matrix().astype('float32'))
azimuth_u = pymc.Normal('azimuth_unc', interp_data.data.foliations['azimuth'].as_matrix().astype('float32'),
1./interp_data.data.foliations['azimuth_std'].as_matrix().astype('float32'))
i = interp_data.get_input_data()
@pymc.deterministic
def grav(value=0, input_ = i, Z_rest_m = Z_rest, Z_ref_m = Z_ref, dip_m = dip_u, azimuth_m = azimuth_u):
# print(Z_rest_m)
input_[1] = dip_m.astype('float32')
input_[2] = azimuth_m.astype('float32')
input_[4][:, 2] = Z_ref_m
input_[5][:, 2] = Z_rest_m
sol = interp_data.th_fn(*input_)
#print("he", sol[0,0,:].sum())
# return solution
lith = sol[-1, 0, :]
lith_s = np.tile(lith, (25,1))[b_].reshape(25,-1)
grav = lith_s* tz_
return grav.sum(axis=1)
# return sol
# @pymc.deterministic
# def grav(value=0, model=gempy_model):
# lith = model[-1, 0, :]
# lith_s = np.tile(lith, (25,1))[b_].reshape(25,-1)
# grav = lith_s* tz_
# return grav.sum(axis=1)
[3, 3, 3]
np.random.seed(123456)
pymc_model = pymc.Model(input = [Z_ref, Z_rest, dip_u, azimuth_u, grav])
MCFS_RUN_16 = pymc.MCMC(pymc_model)#, db='pickle', dbname='MCMC16.pickle')
MCFS_RUN_16.sample(iter=100)
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.064899966100711e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
[------------- 35% ] 35 of 100 complete in 0.5 sec
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.571600425469114e-08 ' condition number: {}'.format(rcond), RuntimeWarning) /home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 4.8928015417004644e-08 ' condition number: {}'.format(rcond), RuntimeWarning) /home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.2049905718831724e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
[-----------------65%---- ] 65 of 100 complete in 1.0 sec
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.1274941625933934e-08 ' condition number: {}'.format(rcond), RuntimeWarning) /home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 4.182751922598982e-08 ' condition number: {}'.format(rcond), RuntimeWarning) /home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.1175184978546895e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
[-----------------93%--------------- ] 93 of 100 complete in 1.5 sec
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 5.718667139831268e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
[-----------------100%-----------------] 100 of 100 complete in 1.7 sec
/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 4.607261772093807e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
MCFS_RUN_16.trace('grav')[:][:, 8] - MCFS_RUN_16.trace('grav')[:][:, 8]
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., 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., 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.])
pymc.Matplot.plot(MCFS_RUN_8.trace('grav'))
Plotting grav_0 Plotting grav_1 Plotting grav_2 Plotting grav_3 Plotting grav_4 Plotting grav_5 Plotting grav_6 Plotting grav_7 Plotting grav_8 Plotting grav_9 Plotting grav_10 Plotting grav_11 Plotting grav_12 Plotting grav_13 Plotting grav_14 Plotting grav_15 Plotting grav_16 Plotting grav_17 Plotting grav_18 Plotting grav_19 Plotting grav_20 Plotting grav_21 Plotting grav_22 Plotting grav_23 Plotting grav_24
MCFS_RUN_8.trace('grav').stats()
/home/miguel/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py:224: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future return reshape(newshape, order=order)
{'95% HPD interval': array([[ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.12171418e-10, 3.31095837e-11, 6.26419005e-11, 6.26419005e-11, 6.48159254e-11, 5.46451502e-11, 6.13873423e-12, 6.91542367e-11, 6.91542367e-11, 7.70816107e-11, 2.02048424e-11, 0.00000000e+00, 3.20366493e-12, 3.20366493e-12, 3.20366493e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.15375083e-10, 1.66260573e-10, 1.91858801e-10, 1.91858801e-10, 2.12891696e-10, 2.11729590e-10, 6.02640924e-11, 2.12312001e-10, 2.12312001e-10, 2.39842258e-10, 1.24317727e-10, 0.00000000e+00, 6.40732986e-12, 6.40732986e-12, 6.40732986e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]), 'mc error': array([ 2.26429295e-15, 1.48340411e-15, 1.48340411e-15, 0.00000000e+00, 8.78814985e-15, 3.62284242e-13, 3.03915536e-13, 3.03915536e-13, 3.74529054e-13, 4.03935452e-13, 1.61350311e-13, 3.48668557e-13, 3.48668557e-13, 3.67599796e-13, 2.84254021e-13, 5.28052466e-14, 7.77665681e-14, 7.77665681e-14, 5.96847860e-14, 4.00230254e-14, 3.37692973e-14, 2.86210585e-14, 2.86210585e-14, 2.50164153e-14, 2.94905301e-14]), 'mean': array([ 1.10567730e-10, 1.10574970e-10, 1.10574970e-10, 1.10569845e-10, 1.12299885e-10, 9.45442974e-11, 1.22151815e-10, 1.22151815e-10, 1.36611500e-10, 1.43757082e-10, 3.52129071e-11, 1.32923797e-10, 1.32923797e-10, 1.52166459e-10, 6.20001626e-11, 4.30657457e-13, 4.90013288e-12, 4.90013288e-12, 5.78073620e-12, 4.99473076e-13, 9.97240459e-14, 1.04510031e-13, 1.04510031e-13, 7.24064340e-14, 8.20552659e-14]), 'n': 10000, 'quantiles': {2.5: array([ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.12171418e-10, 3.92743595e-11, 6.33794528e-11, 6.33794528e-11, 6.49037161e-11, 5.06680837e-11, 8.40731876e-12, 6.46170677e-11, 6.46170677e-11, 8.16187797e-11, 2.64371139e-11, 0.00000000e+00, 3.20366493e-12, 3.20366493e-12, 3.20366493e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 25: array([ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.12171418e-10, 7.57967951e-11, 1.09135500e-10, 1.09135500e-10, 1.14507490e-10, 1.19573808e-10, 2.48355487e-11, 1.08342865e-10, 1.08342865e-10, 1.19131714e-10, 3.72065544e-11, 0.00000000e+00, 3.20366493e-12, 3.20366493e-12, 3.20366493e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 50: array([ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.12171418e-10, 8.42406581e-11, 1.13772975e-10, 1.13772975e-10, 1.33863121e-10, 1.50652459e-10, 3.10678202e-11, 1.29284778e-10, 1.29284778e-10, 1.50023089e-10, 5.95970731e-11, 0.00000000e+00, 3.20366493e-12, 3.20366493e-12, 4.80549739e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 75: array([ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.12171418e-10, 1.10569519e-10, 1.36117739e-10, 1.36117739e-10, 1.58158371e-10, 1.72303704e-10, 4.32623804e-11, 1.57191250e-10, 1.57191250e-10, 1.77236998e-10, 7.80385840e-11, 0.00000000e+00, 4.80549739e-12, 4.80549739e-12, 6.40732986e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 97.5: array([ 1.10569585e-10, 1.10569845e-10, 1.10569845e-10, 1.10569845e-10, 1.15375083e-10, 1.83150105e-10, 1.92791660e-10, 1.92791660e-10, 2.13956295e-10, 2.09687746e-10, 7.11469536e-11, 2.10195184e-10, 2.10195184e-10, 2.48927323e-10, 1.34566627e-10, 0.00000000e+00, 6.40732986e-12, 6.40732986e-12, 6.40732986e-12, 6.40732986e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00])}, 'standard deviation': array([ 2.07384891e-13, 1.51876696e-13, 1.51876696e-13, 2.45569244e-25, 7.96366858e-13, 3.56143768e-11, 3.00626744e-11, 3.00626744e-11, 3.64739469e-11, 4.12920171e-11, 1.62213413e-11, 3.74268200e-11, 3.74268200e-11, 4.17495485e-11, 2.92291615e-11, 5.38072722e-12, 8.41123010e-12, 8.41123010e-12, 6.91166358e-12, 3.97433128e-12, 3.51588763e-12, 2.98574958e-12, 2.98574958e-12, 2.59181864e-12, 3.05998157e-12])}
pymc.Matplot.plot(MCFS_RUN_16.trace('grav'))
Plotting grav_0 Plotting grav_1 Plotting grav_2 Plotting grav_3 Plotting grav_4 Plotting grav_5 Plotting grav_6 Plotting grav_7 Plotting grav_8 Plotting grav_9 Plotting grav_10 Plotting grav_11 Plotting grav_12 Plotting grav_13 Plotting grav_14 Plotting grav_15 Plotting grav_16 Plotting grav_17 Plotting grav_18 Plotting grav_19 Plotting grav_20 Plotting grav_21 Plotting grav_22 Plotting grav_23 Plotting grav_24
MCFS_RUN_16.trace('grav').stats()
/home/miguel/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py:224: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future return reshape(newshape, order=order)
{'95% HPD interval': array([[ 2.08261741e-10, 2.29180491e-10, 2.28424031e-10, 2.28195964e-10, 2.28942868e-10, 6.55591626e-11, 1.75499614e-10, 2.34649796e-10, 2.38569940e-10, 1.91397284e-10, 1.22213502e-12, 8.45563834e-11, 3.16909883e-10, 2.94420915e-10, 6.77742005e-11, 0.00000000e+00, 3.62761143e-12, 5.93633010e-11, 2.75171395e-11, 1.51299609e-13, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 2.53748700e-10, 2.29180491e-10, 2.30386053e-10, 2.33132523e-10, 2.29180491e-10, 3.66883307e-10, 3.08708455e-10, 3.68880699e-10, 4.31454552e-10, 3.23452628e-10, 6.61590358e-11, 2.82424678e-10, 5.09034643e-10, 4.98972272e-10, 3.15359910e-10, 7.91654714e-13, 4.49444076e-11, 2.32832428e-10, 1.33458537e-10, 1.52027352e-11, 0.00000000e+00, 0.00000000e+00, 5.26937619e-12, 1.86193874e-12, 0.00000000e+00]]), 'mc error': array([ 1.32280300e-13, 2.52302460e-15, 6.10991769e-15, 1.62839910e-14, 1.40150209e-15, 1.03143467e-12, 3.30354254e-13, 2.98787924e-13, 4.89603552e-13, 4.05722508e-13, 2.65705530e-13, 4.84876510e-13, 4.33855357e-13, 6.00491229e-13, 7.77127127e-13, 1.19390193e-13, 1.88612003e-13, 5.09868760e-13, 3.07054711e-13, 9.48222087e-14, 1.09001836e-13, 9.43027635e-14, 8.62456843e-14, 6.22716685e-14, 4.60043759e-14]), 'mean': array([ 2.28390274e-10, 2.29183310e-10, 2.29296982e-10, 2.30746867e-10, 2.29199949e-10, 1.63487519e-10, 2.38909647e-10, 3.01904604e-10, 3.40856266e-10, 2.61871352e-10, 2.14441796e-11, 1.81357271e-10, 4.08300912e-10, 3.97357843e-10, 1.78844928e-10, 1.04108959e-12, 2.36025342e-11, 1.38476403e-10, 7.43944307e-11, 5.86322279e-12, 4.68140835e-13, 4.98094698e-13, 2.09434766e-12, 7.73394747e-13, 1.64335849e-13]), 'n': 10000, 'quantiles': {2.5: array([ 2.09347354e-10, 2.29180491e-10, 2.28778637e-10, 2.28195964e-10, 2.29180491e-10, 7.61758390e-11, 1.81044208e-10, 2.34649796e-10, 2.30869822e-10, 1.90003326e-10, 2.92852200e-12, 9.05001686e-11, 3.18375518e-10, 2.93809480e-10, 7.15716061e-11, 0.00000000e+00, 6.56863683e-12, 6.60488841e-11, 3.17572943e-11, 1.06535760e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 25: array([ 2.21242129e-10, 2.29180491e-10, 2.29180491e-10, 2.30008572e-10, 2.29180491e-10, 9.13216070e-11, 2.14663822e-10, 2.84211486e-10, 3.16057242e-10, 2.31635312e-10, 7.71491052e-12, 1.40643815e-10, 3.73541737e-10, 3.61330879e-10, 1.23887725e-10, 2.01732812e-13, 1.49314267e-11, 1.00644523e-10, 5.20595635e-11, 2.65341953e-12, 0.00000000e+00, 0.00000000e+00, 2.49613038e-13, 0.00000000e+00, 0.00000000e+00]), 50: array([ 2.22753379e-10, 2.29180491e-10, 2.29180491e-10, 2.30779939e-10, 2.29180491e-10, 1.13562012e-10, 2.35033889e-10, 3.02546816e-10, 3.45599104e-10, 2.69790447e-10, 1.32204846e-11, 1.78366538e-10, 4.06485266e-10, 3.97536387e-10, 1.67345765e-10, 2.01732812e-13, 2.11077823e-11, 1.32108841e-10, 6.93348134e-11, 3.98811194e-12, 0.00000000e+00, 0.00000000e+00, 1.10410889e-12, 0.00000000e+00, 0.00000000e+00]), 75: array([ 2.37118853e-10, 2.29180491e-10, 2.29180491e-10, 2.31619603e-10, 2.29180491e-10, 2.17325284e-10, 2.59093564e-10, 3.20178990e-10, 3.69559696e-10, 2.90139884e-10, 2.21570167e-11, 2.17055923e-10, 4.40232288e-10, 4.33518563e-10, 2.30306826e-10, 2.01732812e-13, 2.75975810e-11, 1.69402408e-10, 8.99523995e-11, 6.02515929e-12, 0.00000000e+00, 0.00000000e+00, 2.24353670e-12, 9.11446646e-13, 0.00000000e+00]), 97.5: array([ 2.55193226e-10, 2.29180491e-10, 2.31449870e-10, 2.33132523e-10, 2.29655737e-10, 4.17206162e-10, 3.18597294e-10, 3.68880699e-10, 4.25801335e-10, 3.22332753e-10, 1.09290003e-10, 2.93458994e-10, 5.10756084e-10, 4.98575415e-10, 3.23415917e-10, 2.03579624e-12, 5.61506315e-11, 2.46746166e-10, 1.45165187e-10, 2.27618075e-11, 0.00000000e+00, 0.00000000e+00, 6.52476776e-12, 1.86193874e-12, 0.00000000e+00])}, 'standard deviation': array([ 1.26323533e-11, 2.23414083e-13, 5.89942501e-13, 1.42999982e-12, 1.28823514e-13, 9.85194029e-11, 3.47252621e-11, 3.20307953e-11, 4.68242490e-11, 3.70783276e-11, 3.05773952e-11, 5.41270237e-11, 4.93982053e-11, 5.27084015e-11, 6.95314332e-11, 1.25009183e-11, 1.82905050e-11, 4.86057405e-11, 3.01496015e-11, 8.66728684e-12, 1.11647815e-11, 9.77116194e-12, 8.66561346e-12, 5.53211745e-12, 4.24558355e-12])}
pymc.Matplot.plot(MCFS_RUN_32.trace('grav'))
Plotting grav_0 Plotting grav_1 Plotting grav_2 Plotting grav_3 Plotting grav_4 Plotting grav_5 Plotting grav_6 Plotting grav_7 Plotting grav_8 Plotting grav_9 Plotting grav_10 Plotting grav_11 Plotting grav_12 Plotting grav_13 Plotting grav_14 Plotting grav_15 Plotting grav_16 Plotting grav_17 Plotting grav_18 Plotting grav_19 Plotting grav_20 Plotting grav_21 Plotting grav_22 Plotting grav_23 Plotting grav_24
MCFS_RUN_32.trace('grav').stats()
/home/miguel/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py:224: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future return reshape(newshape, order=order)
{'95% HPD interval': array([[ 2.41128359e-10, 2.66219621e-10, 2.64210043e-10, 2.64395469e-10, 2.66746266e-10, 6.54675421e-11, 2.15004651e-10, 2.76844585e-10, 3.07761580e-10, 1.96797564e-10, 6.86011132e-12, 1.63930456e-10, 3.77308708e-10, 3.55409560e-10, 5.78247115e-11, 0.00000000e+00, 1.31700445e-11, 1.97195186e-10, 3.25566082e-11, 1.24780951e-12, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 3.01463944e-10, 2.71813274e-10, 2.75877562e-10, 2.76909724e-10, 2.71260870e-10, 4.30915633e-10, 3.88248849e-10, 4.88037015e-10, 5.45232253e-10, 4.52212958e-10, 8.59399066e-11, 4.54601065e-10, 6.07463510e-10, 5.86945475e-10, 4.19539252e-10, 4.54843413e-12, 7.49238572e-11, 3.98917744e-10, 1.48596854e-10, 3.06827849e-11, 0.00000000e+00, 2.02189115e-12, 1.04845672e-11, 3.63835219e-12, 0.00000000e+00]]), 'mc error': array([ 1.68864305e-13, 1.48781968e-14, 2.73376353e-14, 2.42126182e-14, 1.03226712e-14, 1.12500468e-12, 4.56135708e-13, 5.15659773e-13, 5.46993964e-13, 7.72126255e-13, 3.26777836e-13, 7.93126525e-13, 5.54228377e-13, 5.71838787e-13, 1.14439546e-12, 1.40781739e-13, 2.51306017e-13, 5.05087697e-13, 3.09724394e-13, 1.42092613e-13, 9.98007139e-14, 9.88464287e-14, 1.06505638e-13, 7.33780521e-14, 8.05127312e-14]), 'mean': array([ 2.65290257e-10, 2.68458830e-10, 2.69402380e-10, 2.70428628e-10, 2.69468676e-10, 1.83717871e-10, 3.01388809e-10, 3.83407728e-10, 4.35400730e-10, 3.53446414e-10, 3.48815948e-11, 3.11026892e-10, 4.87911758e-10, 4.71685263e-10, 2.29442350e-10, 2.70644584e-12, 3.98064950e-11, 2.90790679e-10, 8.56491473e-11, 1.11611485e-11, 4.69282128e-13, 1.20853964e-12, 4.53213922e-12, 1.58876183e-12, 3.23097513e-13]), 'n': 10000, 'quantiles': {2.5: array([ 2.44282726e-10, 2.66291951e-10, 2.64480766e-10, 2.64296811e-10, 2.66934583e-10, 7.92477613e-11, 2.20336206e-10, 2.77381745e-10, 2.93930154e-10, 1.82783027e-10, 1.04262091e-11, 1.65490076e-10, 3.80905788e-10, 3.53391530e-10, 6.34883004e-11, 2.44462318e-13, 1.73212878e-11, 1.99216502e-10, 3.88293120e-11, 2.46048383e-12, 0.00000000e+00, 0.00000000e+00, 6.93671192e-13, 1.21696215e-13, 0.00000000e+00]), 25: array([ 2.53055681e-10, 2.68126585e-10, 2.68370705e-10, 2.69235210e-10, 2.68952656e-10, 1.03766210e-10, 2.72169374e-10, 3.51198388e-10, 4.03629665e-10, 2.96325632e-10, 1.91570027e-11, 2.58590975e-10, 4.45013910e-10, 4.29490839e-10, 1.33679957e-10, 7.56120119e-13, 2.75696129e-11, 2.58274737e-10, 6.23540637e-11, 5.22691954e-12, 0.00000000e+00, 1.17614215e-13, 1.76862130e-12, 3.59653395e-13, 0.00000000e+00]), 50: array([ 2.57854161e-10, 2.68268611e-10, 2.68918142e-10, 2.70242019e-10, 2.69740861e-10, 1.32305413e-10, 2.98686876e-10, 3.82221052e-10, 4.41844709e-10, 3.76102874e-10, 2.64736092e-11, 3.12461041e-10, 4.84339742e-10, 4.72760344e-10, 2.14873626e-10, 1.15688694e-12, 3.45363046e-11, 2.87944100e-10, 8.10139983e-11, 7.65127350e-12, 0.00000000e+00, 2.96794795e-13, 2.66698672e-12, 6.74765118e-13, 0.00000000e+00]), 75: array([ 2.77744806e-10, 2.68517110e-10, 2.70071373e-10, 2.71581641e-10, 2.70075476e-10, 2.34282592e-10, 3.28152644e-10, 4.16013085e-10, 4.75166259e-10, 4.20879469e-10, 3.77505924e-11, 3.62095207e-10, 5.27285097e-10, 5.14800896e-10, 3.26102230e-10, 1.91400131e-12, 4.40480537e-11, 3.19600861e-10, 1.02994479e-10, 1.14496794e-11, 0.00000000e+00, 6.54219625e-13, 4.05062466e-12, 1.21052298e-12, 0.00000000e+00]), 97.5: array([ 3.08347455e-10, 2.71947252e-10, 2.76211609e-10, 2.76864047e-10, 2.71500502e-10, 4.92177219e-10, 3.95685718e-10, 4.88872688e-10, 5.37161858e-10, 4.46605754e-10, 1.20213950e-10, 4.57109099e-10, 6.12110202e-10, 5.85978652e-10, 4.27691430e-10, 7.46535237e-12, 9.44205898e-11, 4.02102393e-10, 1.61886092e-10, 4.47680368e-11, 0.00000000e+00, 4.42617916e-12, 2.02681917e-11, 7.16953804e-12, 0.00000000e+00])}, 'standard deviation': array([ 1.70055923e-11, 1.50912046e-12, 2.67834944e-12, 2.80979729e-12, 1.13948541e-12, 1.13730483e-10, 4.34608467e-11, 5.22170804e-11, 5.99900892e-11, 7.87589601e-11, 3.29120365e-11, 7.55620053e-11, 5.99674547e-11, 6.02256128e-11, 1.09908629e-10, 1.43744719e-11, 2.55873463e-11, 5.06251314e-11, 3.25387740e-11, 1.42533264e-11, 9.79005165e-12, 1.02579171e-11, 1.09610389e-11, 7.69586455e-12, 8.46175593e-12])}
MCFS_RUN_8
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-12-597cccdee7a4> in <module>() ----> 1 MCFS_RUN_8 NameError: name 'MCFS_RUN_8' is not defined
plt.imshow(MCFS_RUN.trace("grav", chain=-1)[8].reshape(25,25),
cmap='viridis', origin='lower', extent=[7.050000e+05,747000,6863000,6950000])
<matplotlib.image.AxesImage at 0x7f15da0a3908>
MCFS_RUN.trace("grav", chain=-1)[0].sum(axis=1).reshape(25,25)
array([[ 3.19552082e-10, 3.07991459e-10, 2.92444738e-10, 2.87579093e-10, 2.84165329e-10, 2.81994144e-10, 2.80917395e-10, 2.80449998e-10, 2.80204608e-10, 2.80006593e-10, 2.79823936e-10, 2.79578611e-10, 2.79219379e-10, 2.78892797e-10, 2.78618848e-10, 2.78345920e-10, 2.78146705e-10, 2.78007418e-10, 2.77847951e-10, 2.77680208e-10, 2.77615464e-10, 2.77641539e-10, 2.77677996e-10, 2.77618973e-10, 2.77462932e-10], [ 3.48618204e-10, 3.32695833e-10, 3.10157260e-10, 3.02647577e-10, 2.97253676e-10, 2.93417443e-10, 2.90859408e-10, 2.89365574e-10, 2.88511098e-10, 2.87870467e-10, 2.87239086e-10, 2.86682567e-10, 2.86177215e-10, 2.85687192e-10, 2.85174714e-10, 2.84714551e-10, 2.84216459e-10, 2.83825571e-10, 2.83508315e-10, 2.83146926e-10, 2.82759727e-10, 2.82465553e-10, 2.82365292e-10, 2.82276505e-10, 2.82104501e-10], [ 3.85787409e-10, 3.65462854e-10, 3.36665405e-10, 3.26597117e-10, 3.18726584e-10, 3.12470831e-10, 3.08061197e-10, 3.04849587e-10, 3.02650740e-10, 3.01275243e-10, 3.00499858e-10, 2.99784066e-10, 2.99111748e-10, 2.98460501e-10, 2.97875129e-10, 2.97430461e-10, 2.97158031e-10, 2.96683721e-10, 2.95708563e-10, 2.94636045e-10, 2.93634050e-10, 2.92832879e-10, 2.91767281e-10, 2.91574354e-10, 2.91219237e-10], [ 4.28681874e-10, 4.02674847e-10, 3.66736331e-10, 3.54513581e-10, 3.44642802e-10, 3.36952182e-10, 3.31286726e-10, 3.26929139e-10, 3.23991673e-10, 3.22042009e-10, 3.20851593e-10, 3.19988679e-10, 3.19162360e-10, 3.18491879e-10, 3.18152251e-10, 3.18890826e-10, 3.22982908e-10, 3.26241383e-10, 3.21983651e-10, 3.16115503e-10, 3.12913533e-10, 3.10783333e-10, 3.08320464e-10, 3.07572252e-10, 3.07003403e-10], [ 4.75324034e-10, 4.40336428e-10, 3.95355683e-10, 3.81103432e-10, 3.70295932e-10, 3.62152218e-10, 3.56005607e-10, 3.51585034e-10, 3.48726651e-10, 3.47414295e-10, 3.47126206e-10, 3.47497260e-10, 3.47775677e-10, 3.47731230e-10, 3.47754264e-10, 3.50474970e-10, 3.63788520e-10, 3.77860800e-10, 3.70926827e-10, 3.55418988e-10, 3.45786179e-10, 3.40378850e-10, 3.35262077e-10, 3.33813683e-10, 3.32659130e-10], [ 5.29230838e-10, 4.81788186e-10, 4.23973968e-10, 4.08070915e-10, 3.96487306e-10, 3.87584186e-10, 3.81131387e-10, 3.77093112e-10, 3.75368560e-10, 3.75667710e-10, 3.77628013e-10, 3.80628592e-10, 3.84052651e-10, 3.86878322e-10, 3.88438173e-10, 3.91267133e-10, 4.05117305e-10, 4.24014716e-10, 4.23594771e-10, 4.08020394e-10, 3.94331377e-10, 3.85966880e-10, 3.77985797e-10, 3.75458956e-10, 3.73234676e-10], [ 5.75285298e-10, 5.23866926e-10, 4.55787172e-10, 4.44245087e-10, 4.32120682e-10, 4.18005049e-10, 4.07904966e-10, 4.03480101e-10, 4.03802238e-10, 4.09524600e-10, 4.17360272e-10, 4.22824821e-10, 4.29268721e-10, 4.36204023e-10, 4.39016875e-10, 4.37973157e-10, 4.40801821e-10, 4.50581085e-10, 4.57350783e-10, 4.53759812e-10, 4.45398227e-10, 4.38260870e-10, 4.28404587e-10, 4.24472386e-10, 4.21072906e-10], [ 5.59431086e-10, 5.38044442e-10, 4.91311179e-10, 4.90998717e-10, 4.80587507e-10, 4.59115543e-10, 4.41571692e-10, 4.34479784e-10, 4.36576489e-10, 4.48764297e-10, 4.65968301e-10, 4.79851258e-10, 4.92975390e-10, 5.06241415e-10, 5.10905951e-10, 5.05119537e-10, 4.94097045e-10, 4.87287955e-10, 4.86716363e-10, 4.83791179e-10, 4.76980243e-10, 4.69694147e-10, 4.56397478e-10, 4.50960924e-10, 4.46347808e-10], [ 4.41056578e-10, 4.60431278e-10, 4.83943349e-10, 5.02806735e-10, 5.06189615e-10, 4.92489857e-10, 4.73969286e-10, 4.63409220e-10, 4.64594134e-10, 4.77280984e-10, 4.98803401e-10, 5.24244176e-10, 5.51907013e-10, 5.74916095e-10, 5.83748080e-10, 5.74231853e-10, 5.51775620e-10, 5.28149729e-10, 5.11265773e-10, 4.97996575e-10, 4.84783543e-10, 4.72826227e-10, 4.52960123e-10, 4.45021239e-10, 4.38548257e-10], [ 2.75390695e-10, 3.01815945e-10, 3.78401706e-10, 4.24281839e-10, 4.60790493e-10, 4.78379227e-10, 4.75053141e-10, 4.64907555e-10, 4.63846245e-10, 4.74449867e-10, 4.93987619e-10, 5.21614052e-10, 5.53595538e-10, 5.82313258e-10, 5.98101446e-10, 5.91935015e-10, 5.66951146e-10, 5.37930884e-10, 5.13577189e-10, 4.93567702e-10, 4.73983770e-10, 4.54775263e-10, 4.22111645e-10, 4.09815369e-10, 4.00435072e-10], [ 1.54231127e-10, 1.70403548e-10, 2.33003763e-10, 2.88449758e-10, 3.56917081e-10, 4.10412115e-10, 4.31678646e-10, 4.31396403e-10, 4.31810519e-10, 4.43077008e-10, 4.62068776e-10, 4.85653907e-10, 5.11817457e-10, 5.37501046e-10, 5.56168129e-10, 5.56887111e-10, 5.40117598e-10, 5.17384968e-10, 4.95260372e-10, 4.74336583e-10, 4.52421063e-10, 4.25743348e-10, 3.60494601e-10, 3.34710290e-10, 3.19098681e-10], [ 8.65464689e-11, 9.67193013e-11, 1.32050966e-10, 1.75053618e-10, 2.44713602e-10, 3.15236713e-10, 3.58640854e-10, 3.73432207e-10, 3.79486970e-10, 3.93878940e-10, 4.15915013e-10, 4.39432641e-10, 4.61747514e-10, 4.81432998e-10, 4.95462425e-10, 4.99334496e-10, 4.94101456e-10, 4.84413523e-10, 4.69758478e-10, 4.48221461e-10, 4.19583876e-10, 3.77525520e-10, 2.66983353e-10, 2.27965559e-10, 2.07053805e-10], [ 4.60682531e-11, 5.25459083e-11, 7.50826975e-11, 1.00477587e-10, 1.47754875e-10, 2.10723487e-10, 2.64737545e-10, 2.96831600e-10, 3.15277443e-10, 3.35662081e-10, 3.61120205e-10, 3.87623366e-10, 4.11149707e-10, 4.28375530e-10, 4.38830214e-10, 4.43016468e-10, 4.44714592e-10, 4.45091742e-10, 4.33275228e-10, 3.99872960e-10, 3.48575343e-10, 2.84543259e-10, 1.66880701e-10, 1.35560437e-10, 1.18053706e-10], [ 2.24555480e-11, 2.69394251e-11, 4.09170505e-11, 5.45118370e-11, 7.80944223e-11, 1.15583200e-10, 1.59714044e-10, 1.98877367e-10, 2.30249960e-10, 2.59416665e-10, 2.93641823e-10, 3.33629833e-10, 3.66730986e-10, 3.86535608e-10, 3.96977878e-10, 4.01671223e-10, 4.02761226e-10, 3.96405842e-10, 3.66677001e-10, 3.08246427e-10, 2.38709297e-10, 1.75056984e-10, 9.60296614e-11, 7.77407745e-11, 6.61473105e-11], [ 8.18459705e-12, 1.11277425e-11, 2.07247895e-11, 2.90347005e-11, 4.05123204e-11, 5.73906518e-11, 8.11800075e-11, 1.08778280e-10, 1.38187532e-10, 1.71494641e-10, 2.18712087e-10, 2.82791644e-10, 3.40467049e-10, 3.70953688e-10, 3.82264962e-10, 3.84132414e-10, 3.74084064e-10, 3.41150351e-10, 2.77622926e-10, 2.02059807e-10, 1.41160949e-10, 1.00134492e-10, 5.57378342e-11, 4.34887270e-11, 3.51663132e-11], [ 2.23454009e-12, 3.75911328e-12, 1.03948435e-11, 1.59308372e-11, 2.33901368e-11, 3.30086615e-11, 4.56960566e-11, 6.26359490e-11, 8.58481706e-11, 1.19560171e-10, 1.75349643e-10, 2.57176738e-10, 3.34772882e-10, 3.76256440e-10, 3.85708340e-10, 3.74727173e-10, 3.41543337e-10, 2.77528201e-10, 1.93328100e-10, 1.24711385e-10, 8.46277618e-11, 6.06581112e-11, 3.20941941e-11, 2.37088840e-11, 1.76765991e-11], [ 7.01680276e-13, 1.37245619e-12, 6.25819131e-12, 1.01606726e-11, 1.53472914e-11, 2.24329731e-11, 3.20797228e-11, 4.59267403e-11, 6.76262935e-11, 1.05732963e-10, 1.73612145e-10, 2.63477665e-10, 3.35166238e-10, 3.66808217e-10, 3.65265196e-10, 3.37970417e-10, 2.83728111e-10, 2.05084349e-10, 1.28019325e-10, 7.93227499e-11, 5.39016553e-11, 3.81613244e-11, 1.87408935e-11, 1.25677595e-11, 8.54033736e-12], [ 1.82629224e-13, 5.25119306e-13, 3.95882622e-12, 7.00180923e-12, 1.09931608e-11, 1.65192700e-11, 2.44251257e-11, 3.65695605e-11, 5.76175717e-11, 1.00877612e-10, 1.77438742e-10, 2.61979228e-10, 3.13843328e-10, 3.26157987e-10, 3.07686447e-10, 2.63748792e-10, 1.99583995e-10, 1.30050497e-10, 7.84248633e-11, 5.00871419e-11, 3.43442977e-11, 2.41976378e-11, 1.07016392e-11, 6.45119190e-12, 3.66627817e-12], [ 3.84130435e-14, 1.62605648e-13, 2.63388801e-12, 4.84516684e-12, 7.90451228e-12, 1.21857590e-11, 1.84695814e-11, 2.83382994e-11, 4.62884274e-11, 8.60919812e-11, 1.55753656e-10, 2.25873052e-10, 2.61586069e-10, 2.56119358e-10, 2.18371760e-10, 1.64369500e-10, 1.11121499e-10, 7.00168082e-11, 4.50636322e-11, 3.08430864e-11, 2.19244441e-11, 1.57199273e-11, 6.36282017e-12, 3.53802083e-12, 1.64423592e-12], [ 5.82001170e-15, 1.00722161e-13, 1.84860124e-12, 3.34948735e-12, 5.50699588e-12, 8.70290435e-12, 1.31702235e-11, 2.00640391e-11, 3.18036802e-11, 5.57258557e-11, 1.00364316e-10, 1.52137337e-10, 1.77393359e-10, 1.61052465e-10, 1.19490338e-10, 8.01591577e-11, 5.29691861e-11, 3.62181865e-11, 2.59681359e-11, 1.89983408e-11, 1.38979447e-11, 9.82135454e-12, 3.60031991e-12, 1.69908050e-12, 6.05863353e-13], [ 0.00000000e+00, 4.14349520e-14, 1.14205836e-12, 2.27870164e-12, 3.73209312e-12, 5.83073330e-12, 8.75944603e-12, 1.29006264e-11, 1.88924842e-11, 2.85728228e-11, 4.59171206e-11, 7.01470926e-11, 8.33607595e-11, 7.20515258e-11, 5.05445690e-11, 3.52616853e-11, 2.63093524e-11, 2.01023039e-11, 1.52680936e-11, 1.16161596e-11, 8.52330002e-12, 5.73823653e-12, 1.66477153e-12, 6.26970629e-13, 1.58201204e-13], [ 0.00000000e+00, 2.53722083e-14, 5.40548627e-13, 1.26659072e-12, 2.36638825e-12, 3.69288270e-12, 5.36081629e-12, 7.67218304e-12, 1.05027481e-11, 1.41728678e-11, 1.86568926e-11, 2.43187895e-11, 2.78205610e-11, 2.51241274e-11, 2.04164343e-11, 1.68129133e-11, 1.37694869e-11, 1.12405977e-11, 8.83476650e-12, 6.60714486e-12, 4.59530264e-12, 2.85196555e-12, 5.76909462e-13, 1.57261522e-13, 3.19680440e-14], [ 0.00000000e+00, 0.00000000e+00, 2.09882898e-13, 5.56406182e-13, 1.22327564e-12, 1.95692515e-12, 2.95195755e-12, 4.12519317e-12, 5.53776097e-12, 7.03020904e-12, 8.50667705e-12, 9.67871532e-12, 1.02370087e-11, 1.00495316e-11, 9.26888447e-12, 8.22520479e-12, 7.08565683e-12, 5.69716558e-12, 4.42979040e-12, 3.13522197e-12, 2.02656462e-12, 1.06757921e-12, 1.18802241e-13, 6.09294653e-15, 0.00000000e+00], [ 0.00000000e+00, 0.00000000e+00, 6.19050598e-14, 1.82349262e-13, 4.63012005e-13, 7.73693259e-13, 1.30321267e-12, 1.92680377e-12, 2.57047086e-12, 3.17774365e-12, 3.71655111e-12, 4.08428321e-12, 4.24797988e-12, 4.18328121e-12, 3.95137488e-12, 3.60605077e-12, 3.10622861e-12, 2.39043406e-12, 1.75809184e-12, 1.11711362e-12, 6.08813920e-13, 2.33430758e-13, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 6.09289458e-15, 7.16611589e-14, 2.01023814e-13, 4.74061750e-13, 7.40613298e-13, 9.77019298e-13, 1.16198553e-12, 1.31455504e-12, 1.44331308e-12, 1.50199211e-12, 1.47523655e-12, 1.40197410e-12, 1.28113452e-12, 1.07424857e-12, 7.43729060e-13, 4.54497188e-13, 2.19371606e-13, 5.99110129e-14, 6.09294653e-15, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
(interp_data.data.grid.grid[:, 2].max() - interp_data.data.grid.grid[:, 2].min()) *0.9
0.18718532170546623
interp_data.resolution
array([50, 50, 50])
tz_, b_ = precompute_grav(interp_data)
0.1908919617392378 <GeoPhysics.GeoPhysicsPreprocessing object at 0x7f15ee7836a0>
%debug
> <ipython-input-52-0e9e43bb832d>(62)precompute_grav() 60 61 # Selection of the x y z coordenates for each measurement and computed distance in those coordinates ---> 62 s_gr_x = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, :25].T[b].reshape(25,-1) - ab_g[:,0].reshape(25,-1)).astype('float') 63 s_gr_y = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, 25:50].T[b].reshape(25,-1) - ab_g[:,1].reshape(25,-1)).astype('float') 64 s_gr_z = (np.repeat(interp_data.data.grid.grid, 25, axis=1)[:, -25:].T[b].reshape(25,-1) - ab_g[:,2].reshape(25,-1)).astype('float') ipdb> range <class 'range'> ipdb> range_ 0.18718532170546623 ipdb> i_0 0 ipdb> i_1 25 ipdb> ab_g array([[0.3372890540860558, 0.18611927452817562, 0.5528433277493487], [0.34800549097536115, 0.18611927452817562, 0.5528433277493487], [0.36943836475397174, 0.18611927452817562, 0.5528433277493487], [0.3908712385325824, 0.18611927452817562, 0.5528433277493487], [0.40158767542188767, 0.18611927452817562, 0.5528433277493487], [0.4230205492004983, 0.18611927452817562, 0.5528433277493487], [0.44445342297910895, 0.18611927452817562, 0.5528433277493487], [0.4551698598684143, 0.18611927452817562, 0.5528433277493487], [0.4766027336470249, 0.18611927452817562, 0.5528433277493487], [0.4980356074256355, 0.18611927452817562, 0.5528433277493487], [0.5087526878288355, 0.18611927452817562, 0.5528433277493487], [0.5301855616074462, 0.18611927452817562, 0.5528433277493487], [0.5516184353860568, 0.18611927452817562, 0.5528433277493487], [0.5730513091646675, 0.18611927452817562, 0.5528433277493487], [0.5837677460539727, 0.18611927452817562, 0.5528433277493487], [0.6052006198325834, 0.18611927452817562, 0.5528433277493487], [0.626633493611194, 0.18611927452817562, 0.5528433277493487], [0.6373499305004994, 0.18611927452817562, 0.5528433277493487], [0.65878280427911, 0.18611927452817562, 0.5528433277493487], [0.6802163215716154, 0.18611927452817562, 0.5528433277493487], [0.6909327584609206, 0.18611927452817562, 0.5528433277493487], [0.7123656322395313, 0.18611927452817562, 0.5528433277493487], [0.7337985060181419, 0.18611927452817562, 0.5528433277493487], [0.7445149429074472, 0.18611927452817562, 0.5528433277493487], [0.7659478166860578, 0.18611927452817562, 0.5528433277493487]], dtype=object) ipdb> exit
7523.4
import pymc3 as pm
# with pm.Model():
# Z = pm.Normal('Z_unc', interp_data.data.interfaces['Z'].as_matrix().astype('float'),
# interp_data.data.interfaces['Z_std'].as_matrix().astype('float'))
import theano
import theano.tensor as T
geomodel = theano.OpFromGraph( interp_data.interpolator.tg.input_parameters_list(),
[interp_data.interpolator.tg.whole_block_model(0)], on_unused_input='ignore',
)
input_data_P = interp_data.get_input_data()
[3, 3]
# This is the creation of the model
import pymc3 as pm
theano.config.compute_test_value = 'off'
theano.config.warn_float64 = 'warn'
model = pm.Model()
with model:
# Stochastic value
Z_rest = pm.Normal('Z_unc_rest',
interp_data.interpolator.pandas_rest_layer_points['Z'].as_matrix().astype('float32'),
interp_data.interpolator.pandas_rest_layer_points['Z_std'].as_matrix().astype('float32'),
dtype='float32', shape = (66))
Z_ref = pm.Normal('Z_unc_ref', interp_data.interpolator.ref_layer_points[:, 2].astype('float32'),
interp_data.interpolator.ref_layer_points[:, 2].astype('float32'),
dtype='float32', shape = (66))
# We convert a python variable to theano.shared
input_sh = []
for i in input_data_P:
input_sh.append(theano.shared(i))
# We add the stochastic value to the correspondant array
input_sh[4] = T.set_subtensor(
input_sh[4][:, 2], Z_ref)
input_sh[5] = T.set_subtensor(
input_sh[5][:, 2], Z_rest)
geo_model = pm.Deterministic('GeMpy', geomodel(input_sh[0], input_sh[1], input_sh[2],
input_sh[3], input_sh[4], input_sh[5]))
/home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:177: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. update_mapping=update_mapping) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. nw_x = x.type() /home/miguel/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py:242: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. tmp_replace = [(x, x.type()) for x, y in items] /home/miguel/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_op.py:660: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [t() for t in self.output_types]) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_opt.py:555: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. new_outer = outside_ins.dimshuffle(new_ord) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:1474: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. optimizer_profile = optimizer(fgraph) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_op.py:660: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [t() for t in self.output_types]) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py:96: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. [clone_d[i] for i in owner.inputs], strict=rebuild_strict) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:177: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. update_mapping=update_mapping) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/ops.py:55: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. return gof.Apply(self, [x], [x.type()]) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:1661: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. input_storage=input_storage_lists, storage_map=storage_map) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:1661: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. input_storage=input_storage_lists, storage_map=storage_map) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:1661: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. input_storage=input_storage_lists, storage_map=storage_map) /home/miguel/anaconda3/lib/python3.6/site-packages/theano/compile/function_module.py:1661: UserWarning: You are creating a TensorVariable with float64 dtype. You requested an action via the Theano flag warn_float64={ignore,warn,raise,pdb}. input_storage=input_storage_lists, storage_map=storage_map) /home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 4.2776672870559196e-08 ' condition number: {}'.format(rcond), RuntimeWarning)
theano.config.compute_test_value = 'ignore'
# This is the sampling
# BEFORE RUN THIS FOR LONG CHECK IN THE MODULE THEANOGRAF THAT THE FLAG THEANO OPTIMIZER IS IN 'fast_run'!!
with model:
# backend = pm.backends.ndarray.NDArray('geomodels')
step = pm.NUTS()
trace = pm.sample(30, tune=10, init=None, step=step, )
0%| | 0/40 [00:00<?, ?it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve Ill-conditioned matrix detected. Result is not guaranteed to be accurate. Reciprocal condition number: 4.2776672870559196e-08 ' condition number: {}'.format(rcond), RuntimeWarning) 100%|██████████| 40/40 [00:25<00:00, 2.03it/s]/home/miguel/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:418: UserWarning: Chain 0 contains only 30 samples. % (self._chain_id, n)) /home/miguel/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:448: UserWarning: Chain 0 reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize. 'reparameterize.' % self._chain_id)
input_data_T = interp_data.interpolator.tg.input_parameters_list()
select = interp_data.interpolator.pandas_rest_layer_points['formation'] == 'Reservoir'
import theano
import theano.tensor as T
geomodel = theano.OpFromGraph(input_data_T, [interp_data.interpolator.tg.whole_block_model(0)], on_unused_input='ignore')
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) ~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs) 624 try: --> 625 storage_map[ins] = [self._get_test_value(ins)] 626 compute_map[ins] = [True] ~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in _get_test_value(cls, v) 580 detailed_err_msg = utils.get_variable_trace_string(v) --> 581 raise AttributeError('%s has no test value %s' % (v, detailed_err_msg)) 582 AttributeError: Position of the dips has no test value Backtrace when that variable is created: File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2698, in run_cell interactivity=interactivity, compiler=compiler, result=result) File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2802, in run_ast_nodes if self.run_code(code, result): File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2862, in run_code exec(code_obj, self.user_global_ns, self.user_ns) File "<ipython-input-6-5483455a72af>", line 1, in <module> interp_data = gp.InterpolatorInput(geo_data, compile_theano=False, u_grade=[3,3]) File "../gempy/DataManagement.py", line 941, in __init__ self.interpolator = self.set_interpolator(**kwargs) File "../gempy/DataManagement.py", line 1094, in set_interpolator interpolator = self.InterpolatorClass(geo_data_in, geo_data_in.grid, *args, **kwargs) File "../gempy/DataManagement.py", line 1185, in __init__ self.tg = theanograf.TheanoGraph_pro(dtype=dtype, verbose=verbose,) File "../gempy/theanograf.py", line 94, in __init__ self.dips_position_all = T.matrix("Position of the dips") During handling of the above exception, another exception occurred: ValueError Traceback (most recent call last) <ipython-input-16-486a4325a379> in <module>() 1 import theano 2 import theano.tensor as T ----> 3 geomodel = theano.OpFromGraph(input_data_T, [interp_data.interpolator.tg.whole_block_model(0)], on_unused_input='ignore') ~/PycharmProjects/GeMpy/gempy/theanograf.py in whole_block_model(self, n_faults, compute_all) 1309 dict(input=self.len_series_f[n_faults:], taps=[0, 1]), 1310 dict(input=self.n_formations_per_serie[n_faults:], taps=[0, 1]), -> 1311 dict(input=self.u_grade_T[n_faults:], taps=[0])] 1312 # all_series_pf, updates3 = theano.scan( 1313 # fn=self.compute_a_series, ~/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan.py in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict, return_list) 771 # and outputs that needs to be separated 772 --> 773 condition, outputs, updates = scan_utils.get_updates_and_outputs(fn(*args)) 774 if condition is not None: 775 as_while = True ~/PycharmProjects/GeMpy/gempy/theanograf.py in compute_a_series(self, len_i_0, len_i_1, len_f_0, len_f_1, n_form_per_serie_0, n_form_per_serie_1, u_grade_iter, final_block) 1138 self.u_grade_T_op = u_grade_iter 1139 -> 1140 self.dips_position = self.dips_position_all[len_f_0: len_f_1, :] 1141 self.dips_position_tiled = T.tile(self.dips_position, (self.n_dimensions, 1)) 1142 ~/anaconda3/lib/python3.6/site-packages/theano/tensor/var.py in __getitem__(self, args) 577 self, *theano.tensor.subtensor.Subtensor.collapse( 578 args, --> 579 lambda entry: isinstance(entry, Variable))) 580 581 def take(self, indices, axis=None, mode='raise'): ~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs) 637 raise ValueError( 638 'Cannot compute test value: input %i (%s) of Op %s missing default value. %s' % --> 639 (i, ins, node, detailed_err_msg)) 640 elif config.compute_test_value == 'ignore': 641 # silently skip test ValueError: Cannot compute test value: input 0 (Position of the dips) of Op Subtensor{int64:int64:, ::}(Position of the dips, ScalarFromTensor.0, ScalarFromTensor.0) missing default value. Backtrace when that variable is created: File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2698, in run_cell interactivity=interactivity, compiler=compiler, result=result) File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2802, in run_ast_nodes if self.run_code(code, result): File "/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2862, in run_code exec(code_obj, self.user_global_ns, self.user_ns) File "<ipython-input-6-5483455a72af>", line 1, in <module> interp_data = gp.InterpolatorInput(geo_data, compile_theano=False, u_grade=[3,3]) File "../gempy/DataManagement.py", line 941, in __init__ self.interpolator = self.set_interpolator(**kwargs) File "../gempy/DataManagement.py", line 1094, in set_interpolator interpolator = self.InterpolatorClass(geo_data_in, geo_data_in.grid, *args, **kwargs) File "../gempy/DataManagement.py", line 1185, in __init__ self.tg = theanograf.TheanoGraph_pro(dtype=dtype, verbose=verbose,) File "../gempy/theanograf.py", line 94, in __init__ self.dips_position_all = T.matrix("Position of the dips")
Because now the GeMpy model is a theano operation and not a theano function, to call it we need to use theano variables (with theano functions we call them with python variables). This is very easy to modify, we just need to use theano shared to convert our python input data into theano variables.
The pymc3 objects are already theano variables (pm.Normal and so on). Now the trick is that using the theano function T.set_subtensor, we can change one deterministic value of the input arrays(the ones printed in the cell above) by a stochastic pymc3 object. Then with the new arrays we just have to call the theano operation and pymc will do the rest
# This is the creation of the model
import pymc3 as pm
theano.config.compute_test_value = 'off'
model = pm.Model()
with model:
# Stochastic value
reservoir = pm.Normal('reservoir', np.array([0], dtype='float64')
, sd=np.array([0.09], dtype='float64'), dtype='float64', shape=(1))
# We convert a python variable to theano.shared
ref = theano.shared(input_data_P[4])
rest = theano.shared(input_data_P[5])
# We add the stochastic value to the correspondant array
ref = pm.Deterministic('reference', T.set_subtensor(
ref[T.nonzero(T.cast(select.as_matrix(), "int8"))[0], 2],
ref[T.nonzero(T.cast(select.as_matrix(), "int8"))[0], 2]+reservoir))
rest = pm.Deterministic('rest', T.set_subtensor(
rest[T.nonzero(T.cast(select.as_matrix(), "int8"))[0], 2],
rest[T.nonzero(T.cast(select.as_matrix(), "int8"))[0], 2]+reservoir))#
geo_model = pm.Deterministic('GeMpy', geomodel(theano.shared(input_data_P[0]),
theano.shared(input_data_P[1]),
theano.shared(input_data_P[2]),
theano.shared(input_data_P[3]),
ref, rest))
theano.config.compute_test_value = 'ignore'
# This is the sampling
# BEFORE RUN THIS FOR LONG CHECK IN THE MODULE THEANOGRAF THAT THE FLAG THEANO OPTIMIZER IS IN 'fast_run'!!
with model:
# backend = pm.backends.ndarray.NDArray('geomodels')
step = pm.NUTS()
trace = pm.sample(30, init=None, step=step, )
gp.trace.get_values('GeMpy')[0][-1,0,:])
gp.plot_section(geo_data, trace.get_values('GeMpy')[0][-1, 0, :], 13,
direction='y', plot_data=False)
import matplotlib.pyplot as plt
for i in range(100):
gp.plot_section(geo_data, trace.get_values('GeMpy')[i][-1, 0, :], 13,
direction='y', plot_data=False)
plt.show()
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-17-6b10bc712b5f> in <module>() 1 import matplotlib.pyplot as plt 2 for i in range(100): ----> 3 gp.plot_section(geo_data, trace.get_values('GeMpy')[i][-1, 0, :], 13, 4 direction='y', plot_data=False) 5 plt.show() IndexError: index 30 is out of bounds for axis 0 with size 30
interp_data.interpolator.tg.u_grade_T.get_value()
array([3, 3])