In [1]:
# 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
In [2]:
geo_data = gp.read_pickle('../input_data/sandstone.pickle')
In [3]:
# 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)
Out[3]:
In [4]:
geo_data.interfaces['X_std'] = None
geo_data.interfaces['Y_std'] = 0
geo_data.interfaces['Z_std'] = 100
geo_data.foliations['X_std'] = None
geo_data.foliations['Y_std'] = 0
geo_data.foliations['Z_std'] = 0
In [5]:
geo_data.foliations['dip_std'] = 10
geo_data.foliations['azimuth_std'] = 10
geo_data.foliations.head()
Out[5]:
In [6]:
interp_data = gp.InterpolatorInput(geo_data, compile_theano=False, u_grade=[3,3])
In [7]:
interp_data.data.interfaces.head()
Out[7]:
In [8]:
import pymc3 as pm
In [9]:
# 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'))
In [ ]:
In [10]:
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',
)
In [11]:
input_data_P = interp_data.get_input_data()
In [12]:
# 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]))
In [13]:
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, )
In [15]:
input_data_T = interp_data.interpolator.tg.input_parameters_list()
select = interp_data.interpolator.pandas_rest_layer_points['formation'] == 'Reservoir'
In [16]:
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')
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
In [ ]:
# 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))
In [ ]:
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, )
In [ ]:
gp.trace.get_values('GeMpy')[0][-1,0,:])
In [ ]:
gp.plot_section(geo_data, trace.get_values('GeMpy')[0][-1, 0, :], 13,
direction='y', plot_data=False)
In [17]:
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()
In [18]:
interp_data.interpolator.tg.u_grade_T.get_value()
Out[18]:
In [ ]: