Source code for theanograf

"""
    This file is part of gempy.

    gempy is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    gempy is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with gempy.  If not, see <http://www.gnu.org/licenses/>.


DEP-- I need to update this string
Function that generates the symbolic code to perform the interpolation. Calling this function creates
 both the theano functions for the potential field and the block.

Returns:
    theano function for the potential field
    theano function for the block
"""
import theano
import theano.tensor as T
import numpy as np
import sys

theano.config.openmp_elemwise_minsize = 50000
theano.config.openmp = True

theano.config.optimizer = 'fast_run'
theano.config.floatX = 'float64'

theano.config.exception_verbosity = 'high'
theano.config.compute_test_value = 'off'
theano.config.profile_memory = False
theano.config.scan.debug = False
theano.config.profile = False

[docs]class TheanoGraph_pro(object): """ This class is used to help to divide the construction of the graph into sensical parts. All its methods build a part of the graph. Every method can be seen as a branch and collection of branches until the last method that will be the whole tree. Every part of the graph could be compiled separately but as we increase the complexity the input of each of these methods is more and more difficult to provide (if you are in a branch close to the trunk you need all the results of the branches above) """ def __init__(self, output='geology', verbose=[0], dtype='float32'): """ In the init we need to create all the symbolic parameters that are used in the process. Most of the variables are shared parameters initialized with random values. At this stage we only care about the type and shape of the parameters. After we have the graph built we can update the value of these shared parameters to our data (in the interpolatorClass). Args: u_grade: grade of the drift to compile the right graph. I found out that we can make a graph that takes this as variable so this argument will be deprecated soon verbose (list): name of the nodes you want to print dtype (str): type of float either 32 or 64 """ # Pass the verbose list as property self.verbose = verbose self.compute_all = False theano.config.floatX = dtype # Creation of symbolic parameters # ============= # Constants # ============= # Arbitrary values to get the same results that GeoModeller. These parameters are a mystery for me yet. I have # to ask. In my humble opinion they weight the contribution of the interfaces against the # foliations. self.i_reescale = theano.shared(np.cast[dtype](4.)) self.gi_reescale = theano.shared(np.cast[dtype](2.)) # Number of dimensions. Now it is not too variable anymore self.n_dimensions = 3 # ====================== # INITIALIZE SHARED # ==================== self.u_grade_T = theano.shared(np.arange(2, dtype='int64'), "Grade of the universal drift") self.a_T = theano.shared(np.cast[dtype](1.), "Range") self.c_o_T = theano.shared(np.cast[dtype](1.), 'Covariance at 0') self.nugget_effect_grad_T = theano.shared(np.cast[dtype](0.01)) # -DEP- # self.c_resc = theano.shared(np.cast[dtype](1), "Rescaling factor") self.grid_val_T = theano.shared(np.cast[dtype](np.zeros((2, 3))), 'Coordinates of the grid ' 'points to interpolate') # Shape is 9x2, 9 drift funcitons and 2 points self.universal_grid_matrix_T = theano.shared(np.cast[dtype](np.zeros((9, 2)))) # -DEP- Now I pass it as attribute when I create that part of the graph #self.n_faults = theano.shared(0, 'Number of faults to compute')# self.final_block = theano.shared(np.cast[dtype](np.zeros((1, 3))), "Final block computed") #self.yet_simulated = theano.shared(np.ones(3, dtype='int64'), "Points to be computed yet") # This parameters give me the shape of the different groups of data. I pass all data together and I threshold it # using these values to the different potential fields and formations self.len_series_i = theano.shared(np.arange(2, dtype='int64'), 'Length of interfaces in every series') self.len_series_f = theano.shared(np.arange(2, dtype='int64'), 'Length of foliations in every series') self.n_formations_per_serie = theano.shared(np.arange(3, dtype='int64'), 'List with the number of formations') self.n_formation = theano.shared(np.arange(2,5, dtype='int64'), "Value of the formation") self.number_of_points_per_formation_T = theano.shared(np.zeros(3, dtype='int64')) # ====================== # VAR #======================= self.dips_position_all = T.matrix("Position of the dips") self.dip_angles_all = T.vector("Angle of every dip") self.azimuth_all = T.vector("Azimuth") self.polarity_all = T.vector("Polarity") self.ref_layer_points_all = T.matrix("Reference points for every layer") self.rest_layer_points_all = T.matrix("Rest of the points of the layers") # Tiling dips to the 3 spatial coordinations self.dips_position = self.dips_position_all self.dips_position_tiled = T.tile(self.dips_position, (self.n_dimensions, 1)) # These are subsets of the data for each series. I initialized them as the whole arrays but then they will take # the data of every potential field self.dip_angles = self.dip_angles_all self.azimuth = self.azimuth_all self.polarity = self.polarity_all self.ref_layer_points = self.ref_layer_points_all self.rest_layer_points = self.rest_layer_points_all # Block model out of the faults. It is initialized with shape(0, grid+ ref+rest) so if I do not change it during # the computation it does not have any effect self.u_grade_T_op = theano.shared(0) self.len_points = self.rest_layer_points_all.shape[0] # self.fault_matrix_init = T.zeros((0, self.grid_val_T.shape[0] + 2*self.len_points)) self.len_i_0 = 0 self.len_i_1 = 1 self.potential_field_at_interfaces_values = T.vector('potential_field_at_interfaces_values') # TODO tidy up this initializations # self.final_potential_field_at_formations = T.zeros(self.n_formations_per_serie.get_value()[-1]) # self.pot_value = theano.shared(np.zeros(self.n_formations_per_serie.sum(), dtype='float64'), 'average potential field') self.final_potential_field_at_formations = theano.shared(np.zeros(self.n_formations_per_serie.get_value().sum(), dtype=dtype))#np.array([], ndmin=1, dtype=dtype)) #T.vector('Final value of potential fields at interfacse', dtype=dtype) self.final_potential_field_at_faults = theano.shared(np.zeros(self.n_formations_per_serie.get_value().sum(), dtype=dtype)) self.final_potential_field_at_formations_op = self.final_potential_field_at_formations self.final_potential_field_at_faults_op = self.final_potential_field_at_faults #self.potential_field_at_interfaces_value = theano.shared(np.cast[np.int64](np.zeros((2, 2))), 'Potential field value at each interface',) if output is 'gravity': self.densities = theano.shared(np.cast[dtype](np.zeros(3)), "List with the densities") self.tz = theano.shared(np.cast[dtype](np.zeros((1, 3))), "Component z") self.select = theano.shared(np.cast['int8'](np.zeros(3)), "Select nearby cells") # Init fault relation matrix self.fault_relation = theano.shared(np.array([[0, 1, 0, 1], [0, 0, 1, 1], [0, 0, 0, 1], [0, 0, 0, 0]]), 'relation matrix')
[docs] def input_parameters_list(self): """ Create a list with the symbolic variables to use when we compile the theano function Returns: list: [self.dips_position_all, self.dip_angles_all, self.azimuth_all, self.polarity_all, self.ref_layer_points_all, self.rest_layer_points_all] """ ipl = [self.dips_position_all, self.dip_angles_all, self.azimuth_all, self.polarity_all, self.ref_layer_points_all, self.rest_layer_points_all] return ipl
[docs] @staticmethod def squared_euclidean_distances(x_1, x_2): """ Compute the euclidian distances in 3D between all the points in x_1 and x_2 Args: x_1 (theano.tensor.matrix): shape n_points x number dimension x_2 (theano.tensor.matrix): shape n_points x number dimension Returns: theano.tensor.matrix: Distancse matrix. shape n_points x n_points """ # T.maximum avoid negative numbers increasing stability 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.000001 )) return sqd
[docs] def matrices_shapes(self): """ Get all the lengths of the matrices that form the covariance matrix Returns: length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C """ # Calculating the dimensions of the length_of_CG = self.dips_position_tiled.shape[0] length_of_CGI = self.rest_layer_points.shape[0] # if self.u_grade_T.get_value() == 0: # length_of_U_I = 0*self.u_grade_T # else: # length_of_U_I = 3**self.u_grade_T length_of_U_I = self.u_grade_T_op # Self fault matrix contains the block and the potential field (I am not able to split them). Therefore we need # to divide it by 2 length_of_faults = T.cast(self.fault_matrix.shape[0]/2, 'int64') length_of_C = length_of_CG + length_of_CGI + length_of_U_I + length_of_faults if 'matrices_shapes' in self.verbose: length_of_CG = theano.printing.Print("length_of_CG")(length_of_CG) length_of_CGI = theano.printing.Print("length_of_CGI")(length_of_CGI) length_of_U_I = theano.printing.Print("length_of_U_I")(length_of_U_I) length_of_faults = theano.printing.Print("length_of_faults")(length_of_faults) length_of_C = theano.printing.Print("length_of_C")(length_of_C) return length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C
[docs] def cov_interfaces(self): """ Create covariance function for the interfaces Returns: theano.tensor.matrix: covariance of the interfaces. Shape number of points in rest x number of points in rest """ # Compute euclidian distances sed_rest_rest = self.squared_euclidean_distances(self.rest_layer_points, self.rest_layer_points) sed_ref_rest = self.squared_euclidean_distances(self.ref_layer_points, self.rest_layer_points) sed_rest_ref = self.squared_euclidean_distances(self.rest_layer_points, self.ref_layer_points) sed_ref_ref = self.squared_euclidean_distances(self.ref_layer_points, self.ref_layer_points) # Covariance matrix for interfaces C_I = (self.c_o_T * self.i_reescale * ( (sed_rest_rest < self.a_T) * # Rest - Rest Covariances Matrix (1 - 7 * (sed_rest_rest / self.a_T) ** 2 + 35 / 4 * (sed_rest_rest / self.a_T) ** 3 - 7 / 2 * (sed_rest_rest / self.a_T) ** 5 + 3 / 4 * (sed_rest_rest / self.a_T) ** 7) - ((sed_ref_rest < self.a_T) * # Reference - Rest (1 - 7 * (sed_ref_rest / self.a_T) ** 2 + 35 / 4 * (sed_ref_rest / self.a_T) ** 3 - 7 / 2 * (sed_ref_rest / self.a_T) ** 5 + 3 / 4 * (sed_ref_rest / self.a_T) ** 7)) - ((sed_rest_ref < self.a_T) * # Rest - Reference (1 - 7 * (sed_rest_ref / self.a_T) ** 2 + 35 / 4 * (sed_rest_ref / self.a_T) ** 3 - 7 / 2 * (sed_rest_ref / self.a_T) ** 5 + 3 / 4 * (sed_rest_ref / self.a_T) ** 7)) + ((sed_ref_ref < self.a_T) * # Reference - References (1 - 7 * (sed_ref_ref / self.a_T) ** 2 + 35 / 4 * (sed_ref_ref / self.a_T) ** 3 - 7 / 2 * (sed_ref_ref / self.a_T) ** 5 + 3 / 4 * (sed_ref_ref / self.a_T) ** 7)))) + 1e-7 # Add name to the theano node C_I.name = 'Covariance Interfaces' if str(sys._getframe().f_code.co_name) in self.verbose: C_I = theano.printing.Print('Cov interfaces')(C_I) return C_I
[docs] def cov_gradients(self, verbose=0): """ Create covariance function for the gradiens Returns: theano.tensor.matrix: covariance of the gradients. Shape number of points in dip_pos x number of points in dip_pos """ # Euclidean distances sed_dips_dips = self.squared_euclidean_distances(self.dips_position_tiled, self.dips_position_tiled) # Cartesian distances between dips positions h_u = T.vertical_stack( T.tile(self.dips_position[:, 0] - self.dips_position[:, 0].reshape((self.dips_position[:, 0].shape[0], 1)), self.n_dimensions), T.tile(self.dips_position[:, 1] - self.dips_position[:, 1].reshape((self.dips_position[:, 1].shape[0], 1)), self.n_dimensions), T.tile(self.dips_position[:, 2] - self.dips_position[:, 2].reshape((self.dips_position[:, 2].shape[0], 1)), self.n_dimensions)) # Transpose h_v = h_u.T # Perpendicularity matrix. Boolean matrix to separate cross-covariance and # every gradient direction covariance (block diagonal) perpendicularity_matrix = T.zeros_like(sed_dips_dips) # Cross-covariances of x perpendicularity_matrix = T.set_subtensor( perpendicularity_matrix[0:self.dips_position.shape[0], 0:self.dips_position.shape[0]], 1) # Cross-covariances of y perpendicularity_matrix = T.set_subtensor( perpendicularity_matrix[self.dips_position.shape[0]:self.dips_position.shape[0] * 2, self.dips_position.shape[0]:self.dips_position.shape[0] * 2], 1) # Cross-covariances of z perpendicularity_matrix = T.set_subtensor( perpendicularity_matrix[self.dips_position.shape[0] * 2:self.dips_position.shape[0] * 3, self.dips_position.shape[0] * 2:self.dips_position.shape[0] * 3], 1) # Covariance matrix for gradients at every xyz direction and their cross-covariances C_G = T.switch( T.eq(sed_dips_dips, 0), # This is the condition 0, # If true it is equal to 0. This is how a direction affect another ( # else, following Chiles book (h_u * h_v / sed_dips_dips ** 2) * (( (sed_dips_dips < self.a_T) * # first derivative (-self.c_o_T * ((-14 / self.a_T ** 2) + 105 / 4 * sed_dips_dips / self.a_T ** 3 - 35 / 2 * sed_dips_dips ** 3 / self.a_T ** 5 + 21 / 4 * sed_dips_dips ** 5 / self.a_T ** 7))) + (sed_dips_dips < self.a_T) * # Second derivative self.c_o_T * 7 * (9 * sed_dips_dips ** 5 - 20 * self.a_T ** 2 * sed_dips_dips ** 3 + 15 * self.a_T ** 4 * sed_dips_dips - 4 * self.a_T ** 5) / (2 * self.a_T ** 7)) - (perpendicularity_matrix * (sed_dips_dips < self.a_T) * # first derivative self.c_o_T * ((-14 / self.a_T ** 2) + 105 / 4 * sed_dips_dips / self.a_T ** 3 - 35 / 2 * sed_dips_dips ** 3 / self.a_T ** 5 + 21 / 4 * sed_dips_dips ** 5 / self.a_T ** 7))) ) # Setting nugget effect of the gradients # TODO: This function can be substitued by simply adding the nugget effect to the diag if I remove the condition C_G = T.fill_diagonal(C_G, -self.c_o_T * (-14 / self.a_T ** 2) + self.nugget_effect_grad_T) # Add name to the theano node C_G.name = 'Covariance Gradient' if verbose > 1: theano.printing.pydotprint(C_G, outfile="graphs/" + sys._getframe().f_code.co_name + ".png", var_with_name_simple=True) if str(sys._getframe().f_code.co_name) in self.verbose: C_G = theano.printing.Print('Cov Gradients')(C_G) return C_G
[docs] def cov_interface_gradients(self): """ Create covariance function for the gradiens Returns: theano.tensor.matrix: covariance of the gradients. Shape number of points in rest x number of points in dip_pos """ # Euclidian distances sed_dips_rest = self.squared_euclidean_distances(self.dips_position_tiled, self.rest_layer_points) sed_dips_ref = self.squared_euclidean_distances(self.dips_position_tiled, self.ref_layer_points) # Cartesian distances between dips and interface points # Rest hu_rest = T.vertical_stack( (self.dips_position[:, 0] - self.rest_layer_points[:, 0].reshape( (self.rest_layer_points[:, 0].shape[0], 1))).T, (self.dips_position[:, 1] - self.rest_layer_points[:, 1].reshape( (self.rest_layer_points[:, 1].shape[0], 1))).T, (self.dips_position[:, 2] - self.rest_layer_points[:, 2].reshape( (self.rest_layer_points[:, 2].shape[0], 1))).T ) # Reference point hu_ref = T.vertical_stack( (self.dips_position[:, 0] - self.ref_layer_points[:, 0].reshape( (self.ref_layer_points[:, 0].shape[0], 1))).T, (self.dips_position[:, 1] - self.ref_layer_points[:, 1].reshape( (self.ref_layer_points[:, 1].shape[0], 1))).T, (self.dips_position[:, 2] - self.ref_layer_points[:, 2].reshape( (self.ref_layer_points[:, 2].shape[0], 1))).T ) # Cross-Covariance gradients-interfaces C_GI = self.gi_reescale * ( (hu_rest * (sed_dips_rest < self.a_T) * # first derivative (- self.c_o_T * ((-14 / self.a_T ** 2) + 105 / 4 * sed_dips_rest / self.a_T ** 3 - 35 / 2 * sed_dips_rest ** 3 / self.a_T ** 5 + 21 / 4 * sed_dips_rest ** 5 / self.a_T ** 7))) - (hu_ref * (sed_dips_ref < self.a_T) * # first derivative (- self.c_o_T * ((-14 / self.a_T ** 2) + 105 / 4 * sed_dips_ref / self.a_T ** 3 - 35 / 2 * sed_dips_ref ** 3 / self.a_T ** 5 + 21 / 4 * sed_dips_ref ** 5 / self.a_T ** 7))) ).T # Add name to the theano node C_GI.name = 'Covariance gradient interface' if str(sys._getframe().f_code.co_name)+'_g' in self.verbose: theano.printing.pydotprint(C_GI, outfile="graphs/" + sys._getframe().f_code.co_name + ".png", var_with_name_simple=True) return C_GI
[docs] def universal_matrix(self): """ Create the drift matrices for the potential field and its gradient Returns: theano.tensor.matrix: Drift matrix for the interfaces. Shape number of points in rest x 3**degree drift (except degree 0 that is 0) theano.tensor.matrix: Drift matrix for the gradients. Shape number of points in dips x 3**degree drift (except degree 0 that is 0) """ # # Init # U_I = None # U_G = None # # if self.u_grade_T.get_value() == 1: # # ========================== # # Condition of universality 1 degree # # # Gradients # n = self.dips_position.shape[0] # U_G = T.zeros((n * self.n_dimensions, self.n_dimensions)) # # x # U_G = T.set_subtensor( # U_G[:n, 0], 1) # # y # U_G = T.set_subtensor( # U_G[n:n * 2, 1], 1 # ) # # z # U_G = T.set_subtensor( # U_G[n * 2: n * 3, 2], 1 # ) # # # Interface # # Cartesian distances between reference points and rest # hx = T.stack( # (self.rest_layer_points[:, 0] - self.ref_layer_points[:, 0]), # (self.rest_layer_points[:, 1] - self.ref_layer_points[:, 1]), # (self.rest_layer_points[:, 2] - self.ref_layer_points[:, 2]) # ).T # # U_I = - hx * self.gi_reescale # elif self.u_grade_T.get_value() == 2: # ========================== # Condition of universality 2 degree # Gradients n = self.dips_position.shape[0] U_G = T.zeros((n * self.n_dimensions, 3 * self.n_dimensions)) # x U_G = T.set_subtensor(U_G[:n, 0], 1) # y U_G = T.set_subtensor(U_G[n * 1:n * 2, 1], 1) # z U_G = T.set_subtensor(U_G[n * 2: n * 3, 2], 1) # x**2 U_G = T.set_subtensor(U_G[:n, 3], 2 * self.gi_reescale * self.dips_position[:, 0]) # y**2 U_G = T.set_subtensor(U_G[n * 1:n * 2, 4], 2 * self.gi_reescale * self.dips_position[:, 1]) # z**2 U_G = T.set_subtensor(U_G[n * 2: n * 3, 5], 2 * self.gi_reescale * self.dips_position[:, 2]) # xy U_G = T.set_subtensor(U_G[:n, 6], self.gi_reescale * self.dips_position[:, 1]) # This is y U_G = T.set_subtensor(U_G[n * 1:n * 2, 6], self.gi_reescale * self.dips_position[:, 0]) # This is x # xz U_G = T.set_subtensor(U_G[:n, 7], self.gi_reescale * self.dips_position[:, 2]) # This is z U_G = T.set_subtensor(U_G[n * 2: n * 3, 7], self.gi_reescale * self.dips_position[:, 0]) # This is x # yz U_G = T.set_subtensor(U_G[n * 1:n * 2, 8], self.gi_reescale * self.dips_position[:, 2]) # This is z U_G = T.set_subtensor(U_G[n * 2:n * 3, 8], self.gi_reescale * self.dips_position[:, 1]) # This is y # Interface U_I = - T.stack( (self.gi_reescale * (self.rest_layer_points[:, 0] - self.ref_layer_points[:, 0]), self.gi_reescale * (self.rest_layer_points[:, 1] - self.ref_layer_points[:, 1]), self.gi_reescale * (self.rest_layer_points[:, 2] - self.ref_layer_points[:, 2]), self.gi_reescale ** 2 * (self.rest_layer_points[:, 0] ** 2 - self.ref_layer_points[:, 0] ** 2), self.gi_reescale ** 2 * (self.rest_layer_points[:, 1] ** 2 - self.ref_layer_points[:, 1] ** 2), self.gi_reescale ** 2 * (self.rest_layer_points[:, 2] ** 2 - self.ref_layer_points[:, 2] ** 2), self.gi_reescale ** 2 * ( self.rest_layer_points[:, 0] * self.rest_layer_points[:, 1] - self.ref_layer_points[:, 0] * self.ref_layer_points[:, 1]), self.gi_reescale ** 2 * ( self.rest_layer_points[:, 0] * self.rest_layer_points[:, 2] - self.ref_layer_points[:, 0] * self.ref_layer_points[:, 2]), self.gi_reescale ** 2 * ( self.rest_layer_points[:, 1] * self.rest_layer_points[:, 2] - self.ref_layer_points[:, 1] * self.ref_layer_points[:, 2]), )).T if 'U_I' in self.verbose: U_I = theano.printing.Print('U_I')(U_I) if 'U_G' in self.verbose: U_G = theano.printing.Print('U_G')(U_G) if str(sys._getframe().f_code.co_name)+'_g' in self.verbose: theano.printing.pydotprint(U_I, outfile="graphs/" + sys._getframe().f_code.co_name + "_i.png", var_with_name_simple=True) theano.printing.pydotprint(U_G, outfile="graphs/" + sys._getframe().f_code.co_name + "_g.png", var_with_name_simple=True) # Add name to the theano node if U_I: U_I.name = 'Drift interfaces' U_G.name = 'Drift foliations' return U_I[:, :self.u_grade_T_op], U_G[:, :self.u_grade_T_op]
[docs] def faults_matrix(self): """ This function creates the part of the graph that generates the faults function creating a "block model" at the references and the rest of the points. Then this vector has to be appended to the covariance function Returns: theano.tensor.matrix: Drift matrix for the interfaces. Shape number of points in rest x n faults. This drif is a simple addition of an arbitrary number theano.tensor.matrix: Drift matrix for the gradients. Shape number of points in dips x n faults. For discrete values this matrix will be null since the derivative of a constant is 0 """ length_of_CG, length_of_CGI, length_of_U_I, length_of_faults = self.matrices_shapes()[:4] # self.fault_matrix contains the faults volume of the grid and the rest and ref points. For the drift we need # to make it relative to the reference point interface_loc = self.fault_matrix.shape[1] - 2*self.len_points fault_matrix_at_interfaces_rest = self.fault_matrix[::2, interface_loc + self.len_i_0: interface_loc + self.len_i_1] fault_matrix_at_interfaces_ref = self.fault_matrix[::2, interface_loc+self.len_points+self.len_i_0: interface_loc+self.len_points+self.len_i_1] F_I = (fault_matrix_at_interfaces_ref - fault_matrix_at_interfaces_rest)+0.0001 # As long as the drift is a constant F_G is null F_G = T.zeros((length_of_faults, length_of_CG)) + 0.0001 if str(sys._getframe().f_code.co_name) in self.verbose: F_I = theano.printing.Print('Faults interfaces matrix')(F_I) F_G = theano.printing.Print('Faults gradients matrix')(F_G) return F_I, F_G
[docs] def covariance_matrix(self): """ Set all the previous covariances together in the universal cokriging matrix Returns: theano.tensor.matrix: Multivariate covariance """ # Lengths length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C = self.matrices_shapes() # Individual matrices C_G = self.cov_gradients() C_I = self.cov_interfaces() C_GI = self.cov_interface_gradients() U_I, U_G = self.universal_matrix() F_I, F_G = self.faults_matrix() # ================================= # Creation of the Covariance Matrix # ================================= C_matrix = T.zeros((length_of_C, length_of_C)) # First row of matrices # Set C_G C_matrix = T.set_subtensor(C_matrix[0:length_of_CG, 0:length_of_CG], C_G) # Set CGI C_matrix = T.set_subtensor(C_matrix[0:length_of_CG, length_of_CG:length_of_CG + length_of_CGI], C_GI.T) # Set UG #if not self.u_grade_T.get_value() == 0: C_matrix = T.set_subtensor(C_matrix[0:length_of_CG, length_of_CG+length_of_CGI:length_of_CG+length_of_CGI+length_of_U_I], U_G) # Set FG. I cannot use -index because when is -0 is equivalent to 0 C_matrix = T.set_subtensor(C_matrix[0:length_of_CG, length_of_CG+length_of_CGI+length_of_U_I:], F_G.T) # Second row of matrices # Set C_IG C_matrix = T.set_subtensor(C_matrix[length_of_CG:length_of_CG + length_of_CGI, 0:length_of_CG], C_GI) # Set C_I C_matrix = T.set_subtensor(C_matrix[length_of_CG:length_of_CG + length_of_CGI, length_of_CG:length_of_CG + length_of_CGI], C_I) # Set U_I #if not self.u_grade_T.get_value() == 0: C_matrix = T.set_subtensor(C_matrix[length_of_CG:length_of_CG + length_of_CGI, length_of_CG+length_of_CGI:length_of_CG+length_of_CGI+length_of_U_I], U_I) # Set F_I C_matrix = T.set_subtensor(C_matrix[length_of_CG:length_of_CG + length_of_CGI, length_of_CG+length_of_CGI+length_of_U_I:], F_I.T) # Third row of matrices # Set U_G # if not self.u_grade_T.get_value() == 0: C_matrix = T.set_subtensor(C_matrix[length_of_CG+length_of_CGI:length_of_CG+length_of_CGI+length_of_U_I, 0:length_of_CG], U_G.T) # Set U_I C_matrix = T.set_subtensor(C_matrix[length_of_CG+length_of_CGI:length_of_CG+length_of_CGI+length_of_U_I, length_of_CG:length_of_CG + length_of_CGI], U_I.T) # Fourth row of matrices # Set F_G C_matrix = T.set_subtensor(C_matrix[length_of_CG+length_of_CGI+length_of_U_I:, 0:length_of_CG], F_G) # Set F_I C_matrix = T.set_subtensor(C_matrix[length_of_CG+length_of_CGI+length_of_U_I:, length_of_CG:length_of_CG + length_of_CGI], F_I) # Add name to the theano node C_matrix.name = 'Block Covariance Matrix' if str(sys._getframe().f_code.co_name) in self.verbose: C_matrix = theano.printing.Print('cov_function')(C_matrix) return C_matrix
[docs] def b_vector(self, verbose=0): """ Creation of the independent vector b to solve the kriging system Args: verbose: -deprecated- Returns: theano.tensor.vector: independent vector """ length_of_C = self.matrices_shapes()[-1] # ===================== # Creation of the gradients G vector # Calculation of the cartesian components of the dips assuming the unit module G_x = T.sin(T.deg2rad(self.dip_angles)) * T.sin(T.deg2rad(self.azimuth)) * self.polarity G_y = T.sin(T.deg2rad(self.dip_angles)) * T.cos(T.deg2rad(self.azimuth)) * self.polarity G_z = T.cos(T.deg2rad(self.dip_angles)) * self.polarity G = T.concatenate((G_x, G_y, G_z)) # Creation of the Dual Kriging vector b = T.zeros((length_of_C,)) # G = T.tile(G, (1, 1) b = T.set_subtensor(b[0:G.shape[0]], G) if verbose > 1: theano.printing.pydotprint(b, outfile="graphs/" + sys._getframe().f_code.co_name + "_i.png", var_with_name_simple=True) if str(sys._getframe().f_code.co_name) in self.verbose: b = theano.printing.Print('b vector')(b) # Add name to the theano node b.name = 'b vector' return b
[docs] def solve_kriging(self): """ Solve the kriging system. This has to get substituted by a more efficient and stable method QR decomposition in all likelihood Returns: theano.tensor.vector: Dual kriging parameters """ C_matrix = self.covariance_matrix() b = self.b_vector() # Solving the kriging system import theano.tensor.slinalg b2 = T.tile(b, (1, 1)).T DK_parameters = theano.tensor.slinalg.solve(C_matrix, b2) DK_parameters = DK_parameters.reshape((DK_parameters.shape[0],)) # Add name to the theano node DK_parameters.name = 'Dual Kriging parameters' if str(sys._getframe().f_code.co_name) in self.verbose: DK_parameters = theano.printing.Print(DK_parameters.name)(DK_parameters) return DK_parameters
[docs] def x_to_interpolate(self, verbose=0): """ here I add to the grid points also the references points(to check the value of the potential field at the interfaces). Also here I will check what parts of the grid have been already computed in a previous series to avoid to recompute. Returns: theano.tensor.matrix: The 3D points of the given grid plus the reference and rest points """ #yet_simulated = self.yet_simulated_func() # Removing points no simulated pns = (self.grid_val_T * self.yet_simulated.reshape((self.yet_simulated.shape[0], 1))).nonzero_values() # Adding the rest interface points grid_val = T.vertical_stack(pns.reshape((-1, 3)), self.rest_layer_points_all) # Adding the ref interface points grid_val = T.vertical_stack(grid_val, self.ref_layer_points_all) # Removing points no simulated # grid_val = (grid_val* self.yet_simulated.reshape((self.yet_simulated.shape[0], 1))).nonzero_values() if verbose > 1: theano.printing.pydotprint(grid_val, outfile="graphs/" + sys._getframe().f_code.co_name + ".png", var_with_name_simple=True) if 'grid_val' in self.verbose: grid_val = theano.printing.Print('Points to interpolate')(grid_val) return grid_val
[docs] def extend_dual_kriging(self): """ Tile the dual kriging vector to cover all the points to interpolate.So far I just make a matrix with the dimensions len(DK)x(grid) but in the future maybe I have to try to loop all this part so consume less memory Returns: theano.tensor.matrix: Matrix with the Dk parameters repeated for all the points to interpolate """ grid_val = self.x_to_interpolate() DK_parameters = self.solve_kriging() # Creation of a matrix of dimensions equal to the grid with the weights for every point (big 4D matrix in # ravel form) # TODO IMP: Change the tile by a simple dot op DK_weights = T.tile(DK_parameters, (grid_val.shape[0], 1)).T return DK_weights
[docs] def gradient_contribution(self): """ Computation of the contribution of the foliations at every point to interpolate Returns: theano.tensor.vector: Contribution of all foliations (input) at every point to interpolate """ weights = self.extend_dual_kriging() length_of_CG = self.matrices_shapes()[0] grid_val = self.x_to_interpolate() # Cartesian distances between the point to simulate and the dips hu_SimPoint = T.vertical_stack( (self.dips_position[:, 0] - grid_val[:, 0].reshape((grid_val[:, 0].shape[0], 1))).T, (self.dips_position[:, 1] - grid_val[:, 1].reshape((grid_val[:, 1].shape[0], 1))).T, (self.dips_position[:, 2] - grid_val[:, 2].reshape((grid_val[:, 2].shape[0], 1))).T ) # Euclidian distances sed_dips_SimPoint = self.squared_euclidean_distances(self.dips_position_tiled, grid_val) # Gradient contribution sigma_0_grad = T.sum( (weights[:length_of_CG, :] * self.gi_reescale * (-hu_SimPoint * (sed_dips_SimPoint < self.a_T) * # first derivative (- self.c_o_T * ((-14 / self.a_T ** 2) + 105 / 4 * sed_dips_SimPoint / self.a_T ** 3 - 35 / 2 * sed_dips_SimPoint ** 3 / self.a_T ** 5 + 21 / 4 * sed_dips_SimPoint ** 5 / self.a_T ** 7)))), axis=0) # Add name to the theano node sigma_0_grad.name = 'Contribution of the foliations to the potential field at every point of the grid' return sigma_0_grad
[docs] def interface_contribution(self): """ Computation of the contribution of the interfaces at every point to interpolate Returns: theano.tensor.vector: Contribution of all interfaces (input) at every point to interpolate """ weights = self.extend_dual_kriging() length_of_CG, length_of_CGI = self.matrices_shapes()[:2] grid_val = self.x_to_interpolate() # Euclidian distances sed_rest_SimPoint = self.squared_euclidean_distances(self.rest_layer_points, grid_val) sed_ref_SimPoint = self.squared_euclidean_distances(self.ref_layer_points, grid_val) # Interface contribution sigma_0_interf = (T.sum( -weights[length_of_CG:length_of_CG + length_of_CGI, :] * (self.c_o_T * self.i_reescale * ( (sed_rest_SimPoint < self.a_T) * # SimPoint - Rest Covariances Matrix (1 - 7 * (sed_rest_SimPoint / self.a_T) ** 2 + 35 / 4 * (sed_rest_SimPoint / self.a_T) ** 3 - 7 / 2 * (sed_rest_SimPoint / self.a_T) ** 5 + 3 / 4 * (sed_rest_SimPoint / self.a_T) ** 7) - ((sed_ref_SimPoint < self.a_T) * # SimPoint- Ref (1 - 7 * (sed_ref_SimPoint / self.a_T) ** 2 + 35 / 4 * (sed_ref_SimPoint / self.a_T) ** 3 - 7 / 2 * (sed_ref_SimPoint / self.a_T) ** 5 + 3 / 4 * (sed_ref_SimPoint / self.a_T) ** 7)))), axis=0)) # Add name to the theano node sigma_0_interf.name = 'Contribution of the interfaces to the potential field at every point of the grid' return sigma_0_interf
[docs] def universal_drift_contribution(self): """ Computation of the contribution of the universal drift at every point to interpolate Returns: theano.tensor.vector: Contribution of the universal drift (input) at every point to interpolate """ weights = self.extend_dual_kriging() length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C = self.matrices_shapes() grid_val = self.x_to_interpolate() # Universal drift contribution # Universal terms used to calculate f0 # Here I create the universal terms for rest and ref. The universal terms for the grid are done in python # and append here. The idea is that the grid is kind of constant so I do not have to recompute it every # time _universal_terms_interfaces_rest = T.horizontal_stack( self.rest_layer_points_all, (self.rest_layer_points_all ** 2), T.stack((self.rest_layer_points_all[:, 0] * self.rest_layer_points_all[:, 1], self.rest_layer_points_all[:, 0] * self.rest_layer_points_all[:, 2], self.rest_layer_points_all[:, 1] * self.rest_layer_points_all[:, 2]), axis=1)) _universal_terms_interfaces_ref = T.horizontal_stack( self.ref_layer_points_all, (self.ref_layer_points_all ** 2), T.stack((self.ref_layer_points_all[:, 0] * self.ref_layer_points_all[:, 1], self.ref_layer_points_all[:, 0] * self.ref_layer_points_all[:, 2], self.ref_layer_points_all[:, 1] * self.ref_layer_points_all[:, 2]), axis=1), ) # I append rest and ref to grid universal_grid_interfaces_matrix = T.horizontal_stack( (self.universal_grid_matrix_T * self.yet_simulated).nonzero_values().reshape((9, -1)), T.vertical_stack(_universal_terms_interfaces_rest, _universal_terms_interfaces_ref).T) # T.tile(_universal_terms_interfaces, (1, 2))) # These are the magic terms to get the same as geomodeller gi_rescale_aux = T.repeat(self.gi_reescale, 9) gi_rescale_aux = T.set_subtensor(gi_rescale_aux[:3], 1) _aux_magic_term = T.tile(gi_rescale_aux[:self.u_grade_T_op], (grid_val.shape[0], 1)).T # Drif contribution f_0 = (T.sum( weights[length_of_CG + length_of_CGI:length_of_CG + length_of_CGI + length_of_U_I, :] * self.gi_reescale * _aux_magic_term * universal_grid_interfaces_matrix[:self.u_grade_T_op] , axis=0)) if not type(f_0) == int: f_0.name = 'Contribution of the universal drift to the potential field at every point of the grid' if str(sys._getframe().f_code.co_name) in self.verbose: f_0 = theano.printing.Print('Universal terms contribution')(f_0) return f_0
[docs] def faults_contribution(self): """ Computation of the contribution of the faults drift at every point to interpolate. To get these we need to compute a whole block model with the faults data Returns: theano.tensor.vector: Contribution of the faults drift (input) at every point to interpolate """ weights = self.extend_dual_kriging() length_of_CG, length_of_CGI, length_of_U_I, length_of_faults, length_of_C = self.matrices_shapes() grid_val = self.x_to_interpolate() # Contribution # selecting the voxels that are not computed in the previous pot. field # interface_loc = self.fault_matrix.shape[1] - 2*self.len_points # # # # fault_matrix_at_interfaces_rest = self.fault_matrix[::2, # interface_loc + self.len_i_0: interface_loc + self.len_i_1] # fault_matrix_at_interfaces_ref = self.fault_matrix[::2, # interface_loc+self.len_points+self.len_i_0: interface_loc+self.len_points+self.len_i_1] # # aux_ones = T.ones([2 * self.len_points]) faults_select = T.concatenate((self.yet_simulated, aux_ones)) nfc = faults_select.nonzero_values().shape[0] if 'nfc' in self.verbose: nfc = theano.printing.Print('Faults nfc')(nfc) fault_matrix_selection = (self.fault_matrix[::2, :]+1) * faults_select if 'fault matrix selection' in self.verbose: fault_matrix_selection = theano.printing.Print('f_sle')(fault_matrix_selection) fault_matrix_selection_non_zero0 = fault_matrix_selection.nonzero_values() fault_matrix_selection_non_zero = fault_matrix_selection_non_zero0.reshape((length_of_faults, nfc)) f_1 = T.sum( weights[length_of_CG + length_of_CGI + length_of_U_I:, :] * fault_matrix_selection_non_zero, axis=0) # Add name to the theano node f_1.name = 'Faults contribution' if str(sys._getframe().f_code.co_name) in self.verbose: f_1 = theano.printing.Print('Faults contribution')(f_1) return f_1
[docs] def potential_field_at_all(self): """ Compute the potential field at all the interpolation points, i.e. grid plus rest plus ref Returns: theano.tensor.vector: Potential fields at all points """ sigma_0_grad = self.gradient_contribution() sigma_0_interf = self.interface_contribution() f_0 = self.universal_drift_contribution() f_1 = self.faults_contribution() Z_x = (sigma_0_grad + sigma_0_interf + f_0 + f_1) # Add an arbitrary number at the potential field to get unique values for each of them Z_x += T.repeat(T.cast(1000 - 50*self.n_formation_op[0], "float32"), Z_x.shape[0]) Z_x.name = 'Value of the potential field at every point' if str(sys._getframe().f_code.co_name) in self.verbose: Z_x = theano.printing.Print('Potential field at all points')(Z_x) return Z_x
[docs] def potential_field_at_interfaces(self): """ Potential field at interfaces. To avoid errors I take all the points of rest that belong to one interface and make the average Returns: theano.tensor.vector: Potential field values at the interfaces of a given series """ sigma_0_grad = self.gradient_contribution() sigma_0_interf = self.interface_contribution() f_0 = self.universal_drift_contribution() f_1 = self.faults_contribution() potential_field_interfaces = (sigma_0_grad + sigma_0_interf + f_0 + f_1)[-2*self.len_points: -self.len_points] npf = T.cumsum(T.concatenate((T.stack(0), self.number_of_points_per_formation_T))) # Loop to obtain the average Zx for every intertace def average_potential(dim_a, dim_b, pfi): """ Function to make the average of the potential field at an interface Args: dim: size of the rest values vector per formation pfi: the values of all the rest values potentials Return: theano.tensor.vector: average of the potential per formation """ average = pfi[T.cast(dim_a, "int32"): T.cast(dim_b, "int32")].sum() / T.cast((dim_b - dim_a), 'float32') return average # , {self.pot_value: T.stack([average])} potential_field_interfaces_unique, updates1 = theano.scan( fn=average_potential, outputs_info=None, sequences=dict(input=npf, taps=[0, 1]), non_sequences=potential_field_interfaces) # Add name to the theano node potential_field_interfaces_unique.name = 'Value of the potential field at the interfaces' if str(sys._getframe().f_code.co_name) in self.verbose: potential_field_interfaces_unique = theano.printing.Print(potential_field_interfaces_unique.name)\ (potential_field_interfaces_unique) return potential_field_interfaces_unique
[docs] def block_series(self): """ Compute the part of the block model of a given series (dictated by the bool array yet to be computed) Returns: theano.tensor.vector: Value of lithology at every interpolated point """ # Graph to compute the potential field Z_x = self.potential_field_at_all() # Max and min values of the potential field. # TODO this may be expensive because I guess that is a sort algorithm. We just need a +inf and -inf... I guess max_pot = T.max(Z_x) #T.max(potential_field_unique) + 1 min_pot = T.min(Z_x) #T.min(potential_field_unique) - 1 # Value of the potential field at the interfaces of the computed series self.potential_field_at_interfaces_values = T.sort(self.potential_field_at_interfaces()[self.n_formation_op-1])[::-1] # 1000 and 50 are used here to give values to the scalar fields far enough to not interfere self.potential_field_at_interfaces_values += T.repeat(T.cast(1000 - 50*self.n_formation_op[0], "float32"), self.potential_field_at_interfaces_values.shape[0]) # A tensor with the values to segment potential_field_iter = T.concatenate((T.stack([max_pot]), self.potential_field_at_interfaces_values, T.stack([min_pot]))) if "potential_field_iter" in self.verbose: potential_field_iter = theano.printing.Print("potential_field_iter")(potential_field_iter) # Loop to segment the distinct lithologies def compare(a, b, n_formation, Zx): """ Treshold of the points to interpolate given 2 potential field values. TODO: This function is the one we need to change for a sigmoid function Args: a (scalar): Upper limit of the potential field b (scalar): Lower limit of the potential field n_formation (scalar): Value given to the segmentation, i.e. lithology number Zx (vector): Potential field values at all the interpolated points Returns: theano.tensor.vector: segmented values """ return T.le(Zx, a) * T.ge(Zx, b) * n_formation partial_block, updates2 = theano.scan( fn=compare, outputs_info=None, sequences=[dict(input=potential_field_iter, taps=[0, 1]), self.n_formation_op], non_sequences=Z_x) # For every formation we get a vector so we need to sum compress them to one dimension partial_block = partial_block.sum(axis=0) # Add name to the theano node partial_block.name = 'The chunk of block model of a specific series' if str(sys._getframe().f_code.co_name) in self.verbose: partial_block = theano.printing.Print(partial_block.name)(partial_block) return partial_block
[docs] def compute_a_fault(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, fault_matrix ): """ Function that loops each fault, generating a potential field for each on them with the respective block model Args: len_i_0: Lenght of rest of previous series len_i_1: Lenght of rest for the computed series len_f_0: Lenght of dips of previous series len_f_1: Length of dips of the computed series n_form_per_serie_0: Number of formations of previous series n_form_per_serie_1: Number of formations of the computed series Returns: theano.tensor.matrix: block model derived from the faults that afterwards is used as a drift for the "real" data """ # THIS IS THE FAULTS BLOCK. # ================== # Preparing the data # ================== # compute the youngest fault and consecutively the others # Slice the matrices for the corresponding series # Theano shared self.number_of_points_per_formation_T_op = self.number_of_points_per_formation_T[n_form_per_serie_0: n_form_per_serie_1] self.n_formation_op = self.n_formation[n_form_per_serie_0: n_form_per_serie_1] if 'n_formation' in self.verbose: self.n_formation_op = theano.printing.Print('n_formation_fault')(self.n_formation_op) self.u_grade_T_op = u_grade_iter self.dips_position = self.dips_position_all[len_f_0: len_f_1, :] self.dips_position_tiled = T.tile(self.dips_position, (self.n_dimensions, 1)) # Theano Var self.dip_angles = self.dip_angles_all[len_f_0: len_f_1] self.azimuth = self.azimuth_all[len_f_0: len_f_1] self.polarity = self.polarity_all[len_f_0: len_f_1] self.ref_layer_points = self.ref_layer_points_all[len_i_0: len_i_1, :] self.rest_layer_points = self.rest_layer_points_all[len_i_0: len_i_1, :] # Updating the interface points involved in the iteration. This is important for the fault drift self.len_i_0 = len_i_0 self.len_i_1 = len_i_1 # Extracting a the subset of the fault matrix to the scalar field of the current iterations # self.fault_relation = theano.shared(np.array([[0,1,0,1], # [0,0,1,1], # [0,0,0,1], # [0,0,0,0]]), 'relation matrix') faults_relation_op = self.fault_relation[:, T.cast(self.n_formation_op-1, 'int8')] faults_relation_rep = T.repeat(faults_relation_op, 2) if 'faults_relation' in self.verbose: faults_relation_rep = theano.printing.Print('SELECT')(faults_relation_rep) self.fault_matrix = fault_matrix[T.nonzero(T.cast(faults_relation_rep, "int8"))[0], :] if 'fault_matrix_loop' in self.verbose: self.fault_matrix = theano.printing.Print('self fault matrix')(self.fault_matrix) # ==================== # Computing the fault scalar field # ==================== faults_matrix = self.block_series() aux_ones = T.ones([2*self.len_points]) faults_select = T.concatenate((self.yet_simulated, aux_ones)) # Update the block matrix block_matrix = T.set_subtensor( final_block[0, :], T.cast(T.cast(faults_matrix, 'bool'), 'int8')) # Update the potential field matrix if self.compute_all: potential_field_values = self.potential_field_at_all() block_matrix = T.set_subtensor( block_matrix[1, T.nonzero(T.cast(faults_select, "int8"))[0]], potential_field_values) # Store the potential field at the interfaces self.final_potential_field_at_faults_op = T.set_subtensor(self.final_potential_field_at_faults_op[self.n_formation_op-1], self.potential_field_at_interfaces_values) aux_ind = T.max(self.n_formation_op, 0) # Setting the values of the fault matrix computed in the current iteration fault_matrix = T.set_subtensor(fault_matrix[(aux_ind-1)*2:aux_ind*2, :], block_matrix) return block_matrix, fault_matrix, self.final_potential_field_at_faults_op,
[docs] def 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, fault_block): """ Function that loops each series, generating a potential field for each on them with the respective block model Args: len_i_0: Lenght of rest of previous series len_i_1: Lenght of rest for the computed series len_f_0: Lenght of dips of previous series len_f_1: Length of dips of the computed series n_form_per_serie_0: Number of formations of previous series n_form_per_serie_1: Number of formations of the computed series Returns: theano.tensor.matrix: final block model """ # Setting the fault contribution to kriging from the previous loop self.fault_matrix = fault_block # THIS IS THE FINAL BLOCK. (DO I NEED TO LOOP THE FAULTS FIRST? Yes you do) # ================== # Preparing the data # ================== # Vector that controls the points that have been simulated in previous iterations self.yet_simulated = T.eq(final_block[0, :-2 * self.len_points], 0) self.yet_simulated.name = 'Yet simulated LITHOLOGY node' # Theano shared self.number_of_points_per_formation_T_op = self.number_of_points_per_formation_T[n_form_per_serie_0: n_form_per_serie_1] self.n_formation_op = self.n_formation[n_form_per_serie_0: n_form_per_serie_1] self.u_grade_T_op = u_grade_iter self.dips_position = self.dips_position_all[len_f_0: len_f_1, :] self.dips_position_tiled = T.tile(self.dips_position, (self.n_dimensions, 1)) # Theano Var self.dip_angles = self.dip_angles_all[len_f_0: len_f_1] self.azimuth = self.azimuth_all[len_f_0: len_f_1] self.polarity = self.polarity_all[len_f_0: len_f_1] self.ref_layer_points = self.ref_layer_points_all[len_i_0: len_i_1, :] self.rest_layer_points = self.rest_layer_points_all[len_i_0: len_i_1, :] # For the contribution of the faults I did not find a better way self.len_i_0 = len_i_0 self.len_i_1 = len_i_1 # Printing if 'yet_simulated' in self.verbose: self.yet_simulated = theano.printing.Print(self.yet_simulated.name)(self.yet_simulated) if 'n_formation' in self.verbose: self.n_formation_op = theano.printing.Print('n_formation_series')(self.n_formation_op) # ==================== # Computing the series # ==================== aux_ones = T.ones([2 * self.len_points]) lith_select = T.concatenate((self.yet_simulated, aux_ones)) potential_field_contribution = self.block_series() #[:-2*self.len_points] # Updating the block model with the lithology block final_block = T.set_subtensor( final_block[0, T.nonzero(T.cast(lith_select, "int8"))[0]], potential_field_contribution) # Store the potential field at the interfaces self.final_potential_field_at_formations_op = T.set_subtensor( self.final_potential_field_at_formations_op[self.n_formation_op - 1], self.potential_field_at_interfaces_values) # Update the potential field matrix if self.compute_all: potential_field_values = self.potential_field_at_all()#[:-2*self.len_points] final_block = T.set_subtensor( final_block[1, T.nonzero(T.cast(lith_select, "int8"))[0]], potential_field_values) return final_block, self.final_potential_field_at_formations_op
def compute_geological_model(self, n_faults=0, compute_all=True): # TODO This part should go to the init eventually # Init all if compute_all: # Change the flag to extend the graph in the compute fault and compute series function self.compute_all = True # Init faults block. Here we store the block and potential field results of one iteration fault_block_init = T.zeros((2, self.grid_val_T.shape[0] + 2 * self.len_points)) fault_block_init.name = 'final block of faults init' # Here we store the value of the potential field at interfaces pfai_fault = T.zeros((0, self.n_formations_per_serie[-1])) # Init lithology block. Here we store the block and potential field results lith_block_init = T.zeros((2, self.grid_val_T.shape[0] + 2 * self.len_points)) lith_block_init.name = 'final block of lithologies init' lith_matrix = T.zeros((0, 0, self.grid_val_T.shape[0] + 2 * self.len_points)) pfai_lith = T.zeros((0, self.n_formations_per_serie[-1])) # Init to matrix which contains the block and scalar field of every fault self.fault_matrix = T.zeros((n_faults*2, self.grid_val_T.shape[0] + 2 * self.len_points)) else: # TODO Check if still makes sense to have to condition and if yes debugging # Change the flag to extend the graph in the compute fault and compute series function self.compute_all = False # Init faults block. Here we store the block and potential field results fault_block_init = T.zeros((1, self.grid_val_T.shape[0] + 2 * self.len_points)) fault_block_init.name = 'final block of faults init' fault_matrix = T.zeros((0, self.grid_val_T.shape[0] + 2 * self.len_points)) # Here we store the value of the potential field at interfaces pfai_fault = T.zeros((0, self.n_formations_per_serie[-1])) # Init lithology block. Here we store the block and potential field results lith_block_init = T.zeros((1, self.grid_val_T.shape[0] + 2 * self.len_points)) lith_block_init.name = 'final block of lithologies init' lith_matrix = T.zeros((0, self.grid_val_T.shape[0] + 2 * self.len_points)) pfai_lith = T.zeros((0, self.n_formations_per_serie[-1])) self.fault_matrix = T.zeros((n_faults * 2, self.grid_val_T.shape[0] + 2 * self.len_points)) # Compute Faults if n_faults != 0: # --DEP--? Initialize yet simulated self.yet_simulated = T.eq(fault_block_init[0, :-2 * self.len_points], 0) # Looping fault_loop, updates3 = theano.scan( fn=self.compute_a_fault, outputs_info=[fault_block_init, dict(initial=self.fault_matrix, taps=[-1]), None], # This line may be used for the faults network sequences=[dict(input=self.len_series_i[:n_faults + 1], taps=[0, 1]), dict(input=self.len_series_f[:n_faults + 1], taps=[0, 1]), dict(input=self.n_formations_per_serie[:n_faults + 1], taps=[0, 1]), dict(input=self.u_grade_T[:n_faults + 1], taps=[0])], return_list=True, ) # We return the last iteration of the fault matrix self.fault_matrix = fault_loop[1][-1] # For this we return every iteration since is each potential field at interface pfai_fault = fault_loop[2] # Check if there are lithologies to compute if len(self.len_series_f.get_value()) - 1 > n_faults: # Compute Lithologies lith_loop, updates2 = theano.scan( fn=self.compute_a_series, outputs_info=[lith_block_init, None], sequences=[dict(input=self.len_series_i[n_faults:], taps=[0, 1]), dict(input=self.len_series_f[n_faults:], taps=[0, 1]), dict(input=self.n_formations_per_serie[n_faults:], taps=[0, 1]), dict(input=self.u_grade_T[n_faults:], taps=[0])], non_sequences=[self.fault_matrix], name='Looping interfaces', profile=False, return_list=True ) lith_matrix = lith_loop[0][-1] pfai_lith = lith_loop[1] pfai = T.vertical_stack(pfai_fault, pfai_lith) return [lith_matrix[:, :-2 * self.len_points], self.fault_matrix[:, :-2 * self.len_points], pfai] # ================================== # Geophysics # ================================== def switch_densities(self, n_formation, density, density_block): density_block = T.switch(T.eq(density_block, n_formation), density, density_block) return density_block def compute_forward_gravity(self, n_faults=0, compute_all=True): # densities, tz, select, # TODO: Assert outside that densities is the same size as formations (minus faults) # Compute the geological model lith_matrix, fault_matrix, pfai = self.compute_geological_model(n_faults=n_faults, compute_all=compute_all) if n_faults == 0: formations = T.concatenate([self.n_formation[::-1], T.stack([0])]) else: formations = T.concatenate([self.n_formation[:n_faults-1:-1], T.stack([0])]) formations = theano.printing.Print('formations')(formations) # Substitue lithologies by its density density_block_loop, updates4 = theano.scan(self.switch_densities, outputs_info=[lith_matrix[0]], sequences=[formations, self.densities], return_list = True ) if True: density_block_loop = theano.printing.Print('density block')(density_block_loop[-1]) n_measurements = self.tz.shape[0] # Tiling the density block for each measurent and picking just the closer to them. This has to be possible to # optimize densities_rep = T.tile(density_block_loop[-1], n_measurements) densities_selected = densities_rep[T.nonzero(T.cast(self.select, "int8"))[0]] densities_selected_reshaped = densities_selected.reshape((n_measurements, -1)) # # # density times the component z of gravity grav = densities_selected_reshaped * self.tz #return [lith_matrix, self.fault_matrix, pfai, grav.sum(axis=1)] return [lith_matrix, self.fault_matrix, pfai, grav.sum(axis=1)]