Merge pull request #2 from cgre-aachen/development

Development

gempy/DataManagement.py

@@ -233,7 +233,7 @@ def get_data(self, itype='all', numeric=False, verbosity=0):
             show_par_i = self.interfaces.columns
 
         if numeric:
-            show_par_f = ['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z']
+            show_par_f = ['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z', 'dip', 'azimuth', 'polarity']
             show_par_i = ['X', 'Y', 'Z']
             dtype = 'float'
         if itype == 'foliations':

gempy/Topology.py

@@ -20,9 +20,7 @@
 from skimage.future import graph
 from skimage.measure import label
 from skimage.measure import regionprops
-import networkx as nx
 import numpy as np
-import matplotlib.pyplot as plt
 
 # TODO: Across-fault edge identification
 # TODO: Across-unconformity edge identification

gempy/UncertaintyAnalysisPYMC2.py

@@ -5,6 +5,9 @@
 import theano
 import numpy as np
 import networkx as nx
+import gempy as gp
+import matplotlib.pyplot as plt
+import seaborn as sns
 
 
 class Posterior:
@@ -18,48 +21,224 @@ def __init__(self, dbname, topology=False, verbose=False):
         self.trace_names = self.db.trace_names[0]
         # get gempy block models
         try:
-            self.sols = self.db.gempy_model.gettrace()
+            self.lb, self.fb = self.db.gempy_model.gettrace()
         except AttributeError:
-            print("No GemPy block models tallied.")
-            self.sols = None
+            print("No GemPy model trace tallied.")
+            self.lb = None
+            self.fb = None
 
         if topology:
             # load graphs
             topo_trace = self.db.gempy_topo.gettrace()
-            self.graphs = topo_trace[:, 0]
+            self.topo_graphs = topo_trace[:, 0]
             # load centroids
-            self.centroids = topo_trace[:, 1]
+            self.topo_centroids = topo_trace[:, 1]
+            self.topo_labels_unique = topo_trace[:, 2]
+            self.topo_lith_to_labels_lot = topo_trace[:, 3]
+            self.topo_labels_to_lith_lot = topo_trace[:, 4]
             del topo_trace
 
         # load input data
         self.input_data = self.db.input_data.gettrace()
 
+        self.lith_prob = None
+        self.ie = None
+        self.ie_total = None
+
     def change_input_data(self, interp_data, i):
         """Changes input data in interp_data to posterior input data at iteration i."""
 
         # replace interface data
         interp_data.geo_data_res.interfaces[["X", "Y", "Z"]] = self.input_data[i][0]
         # replace foliation data
-        interp_data.geo_data_res.foliations[["G_x", "G_y", "G_z", "X", "Y", "Z"]] = self.input_data[i][1]
+        interp_data.geo_data_res.foliations[["G_x", "G_y", "G_z", "X", "Y", "Z", "dip", "azimuth", "polarity"]] = self.input_data[i][1]
         # do all the ugly updating stuff
-        # interp_data.interpolator.tg.final_potential_field_at_formations = theano.shared(np.zeros(
-        #     interp_data.interpolator.tg.n_formations_per_serie.get_value().sum(), dtype='float32'))
-        # interp_data.interpolator.tg.final_potential_field_at_faults = theano.shared(np.zeros(
-        #     interp_data.interpolator.tg.n_formations_per_serie.get_value().sum(), dtype='float32'))
-        # interp_data.update_interpolator()
+        interp_data.interpolator.tg.final_potential_field_at_formations = theano.shared(np.zeros(
+            interp_data.interpolator.tg.n_formations_per_serie.get_value().sum(), dtype='float32'))
+        interp_data.interpolator.tg.final_potential_field_at_faults = theano.shared(np.zeros(
+            interp_data.interpolator.tg.n_formations_per_serie.get_value().sum(), dtype='float32'))
+        interp_data.update_interpolator()
         if self.verbose:
             print("interp_data parameters changed.")
 
     def plot_topology_graph(self, i):
         # get centroid values into list
-        centroid_values = [triplet for triplet in self.centroids[i].values()]
+        centroid_values = [triplet for triplet in self.topo_centroids[i].values()]
         # unzip them into seperate lists of x,y,z coordinates
         centroids_x, centroids_y, centroids_z = list(zip(*centroid_values))
         # create new 2d pos dict for plot
         pos_dict = {}
         for j in range(len(centroids_x)):  # TODO: Change this directly to use zip?
             pos_dict[j + 1] = [centroids_x[j], centroids_y[j]]
         # draw
-        nx.draw_networkx(self.graphs[i], pos=pos_dict)
+        nx.draw_networkx(self.topo_graphs[i], pos=pos_dict)
+
+    def compute_posterior_model(self, interp_data, i):
+        self.change_input_data(interp_data, i)
+        return gp.compute_model(interp_data)
+
+    def plot_section(self, interp_data, i, dim, plot_data=False, plot_topo=False):
+        """Deprecated."""
+        self.change_input_data(interp_data, i)
+        lith_block, fault_block = gp.compute_model(interp_data)
+        #plt.imshow(lith_block[-1, 0,:].reshape(dim[0], dim[1], dim[2])[:, 0, :].T, origin="lower",
+        #           cmap=gp.colors.cmap, norm=gp.colors.norm)
+        gp.plot_section(interp_data.geo_data_res, lith_block[0], 0, plot_data=plot_data)
+
+        rs = interp_data.rescaling_factor
+        #if plot_data:
+        #    plt.scatter(interp_data.geo_data_res.interfaces["X"].values,
+        #                interp_data.geo_data_res.interfaces["Z"].values)
+
+        if plot_topo:
+            self.topo_plot_graph(i)
+
+    def topo_plot_graph(self, i):
+        pos_2d = {}
+        for key in self.topo_centroids[i].keys():
+            pos_2d[key] = [self.topo_centroids[i][key][0], self.topo_centroids[i][key][2]]
+        nx.draw_networkx(self.topo_graphs[i], pos=pos_2d)
+
+    def compute_posterior_models_all(self, interp_data, n=None, calc_fb=True):
+        """Computes block models from stored input parameters for all iterations."""
+        if self.lb is None:
+            # create the storage array
+            lb, fb = self.compute_posterior_model(interp_data, 1)
+            lb = lb[0]
+            fb = fb[0]
+            self.lb = np.empty_like(lb)
+            if calc_fb:
+                self.fb = np.empty_like(lb)
+
+            # compute model for every iteration
+            if n is None:
+                n = self.db.getstate()["sampler"]["_iter"]
+            for i in range(n):
+                if i == 0:
+                    lb, fb = self.compute_posterior_model(interp_data, i)
+                    self.lb = lb[0]
+                    self.fb = fb[0]
+                else:
+                    lb, fb = self.compute_posterior_model(interp_data, i)
+                    self.lb = np.vstack((self.lb, lb[0]))
+                    if calc_fb:
+                        self.fb = np.vstack((self.fb, fb[0]))
+        else:
+            print("self.lb already filled with something.")
+
+    def compute_entropy(self, interp_data):
+        """Computes the voxel information entropy of stored block models."""
+        if self.lb is None:
+            return "No models stored in self.lb, please run 'self.compute_posterior_models_all' to generate block" \
+                   " models for all iterations."
+
+        self.lith_prob = compute_prob_lith(self.lb)
+        self.ie = calcualte_ie_masked(self.lith_prob)
+        self.ie_total = calculate_ie_total(self.ie)
+        print("Information Entropy successfully calculated. Stored in self.ie and self.ie_total")
+
+
+def compute_prob_lith(lith_blocks):
+    """Blocks must be just the lith blocks!"""
+    lith_id = np.unique(lith_blocks)
+    lith_count = np.zeros_like(lith_blocks[0:len(lith_id)])
+    for i, l_id in enumerate(lith_id):
+        lith_count[i] = np.sum(lith_blocks == l_id, axis=0)
+    lith_prob = lith_count / len(lith_blocks)
+    return lith_prob
+
+
+def calcualte_ie_masked(lith_prob):
+    ie = np.zeros_like(lith_prob[0])
+    for l in lith_prob:
+        pm = np.ma.masked_equal(l, 0)  # mask where layer prob is 0
+        ie -= (pm * np.ma.log2(pm)).filled(0)
+    return ie
+
+
+def calculate_ie_total(ie, absolute=False):
+    if absolute:
+        return np.sum(ie)
+    else:
+        return np.sum(ie) / np.size(ie)
+
+
+def compare_graphs(G1, G2):
+    intersection = 0
+    union = G1.number_of_edges()
+
+    for edge in G1.edges_iter():
+        if G2.has_edge(edge[0], edge[1]):
+            intersection += 1
+        else:
+            union += 1
+
+    return intersection / union
+
+
+class Plane:
+    def __init__(self, group_id, data_obj):
+        self.group_id = group_id
+        self.data_obj = data_obj
+
+        # create dataframe bool filters for convenience
+        self.fol_f = self.data_obj.foliations["group_id"] == self.group_id
+        self.interf_f = self.data_obj.interfaces["group_id"] == self.group_id
+
+        # get indices for both foliations and interfaces
+        self.interf_i = self.data_obj.interfaces[self.interf_f].index
+        self.fol_i = self.data_obj.foliations[self.fol_f].index[0]
+
+        # normal
+        self.normal = None
+        # centroid
+        self.centroid = None
+        self.refresh()
+
+    # method: give dip, change interfaces accordingly
+    def interf_recalc(self, dip):
+        """Changes the dip of plane and recalculates Z coordinates for the points belonging to it."""
+        # modify the foliation
+        self.data_obj.foliations.set_value(self.fol_i, "dip", dip)
+        # get azimuth
+        az = float(self.data_obj.foliations.iloc[self.fol_i]["azimuth"])
+        # set polarity according to dip
+        if -90 < dip < 90:
+            polarity = 1
+        else:
+            polarity = -1
+        self.data_obj.foliations.set_value(self.fol_i, "polarity", polarity)
+        # modify gradient
+        self.data_obj.foliations.set_value(self.fol_i, "G_x",
+                                           np.sin(np.deg2rad(dip)) * np.sin(np.deg2rad(az)) * polarity)
+        self.data_obj.foliations.set_value(self.fol_i, "G_y",
+                                           np.sin(np.deg2rad(dip)) * np.cos(np.deg2rad(az)) * polarity)
+        self.data_obj.foliations.set_value(self.fol_i, "G_z", np.cos(np.deg2rad(dip)) * polarity)
+
+        # update normal
+        self.normal = self.get_normal()
+        # modify points (Z only so far)
+        a, b, c = self.normal
+        d = -a * self.centroid[0] - b * self.centroid[1] - c * self.centroid[2]
+        for i, row in self.data_obj.interfaces[self.interf_f].iterrows():
+            # iterate over each point and recalculate Z, set Z
+            # x, y, z = row["X"], row["Y"], row["Z"]
+            Z = (a * row["X"] + b * row["Y"] + d) / -c
+            self.data_obj.interfaces.set_value(i, "Z", Z)
+
+    def refresh(self):
+        # normal
+        self.normal = self.get_normal()
+        # centroid
+        self.centroid = [float(self.data_obj.foliations[self.fol_f]["X"]),
+                         float(self.data_obj.foliations[self.fol_f]["Y"]),
+                         float(self.data_obj.foliations[self.fol_f]["Z"])]
+
+    def get_normal(self):
+        """Just returns updated normal vector (values from dataframe)."""
+        normal = [float(self.data_obj.foliations.iloc[self.fol_i]["G_x"]),
+                  float(self.data_obj.foliations.iloc[self.fol_i]["G_y"]),
+                  float(self.data_obj.foliations.iloc[self.fol_i]["G_z"])]
+        return normal
 
 

gempy/theanograf.py

@@ -28,7 +28,7 @@
 import numpy as np
 import sys
 
-theano.config.optimizer = 'fast_compile'
+theano.config.optimizer = 'fast_run'
 theano.config.exception_verbosity = 'high'
 theano.config.compute_test_value = 'off'
 theano.config.floatX = 'float32'

input_data/tutorial_ch3_foliations

@@ -0,0 +1,3 @@
+,G_x,G_y,G_z,X,X_std,Y,Y_std,Z,Z_std,annotations,azimuth,azimuth_std,dip,dip_std,formation,formation number,group_id,isFault,order_series,polarity,series
+0,-0.5169921826496945,-0.00855947322267713,0.8559473222677055,500.0,,100.0,,1148.0,,"${\bf{x}}_{\beta \,{\bf{1}},0}$",269.051481039,,31.1354510371,,Layer 2,1,l2_a,False,1,1,Default serie
+1,0.5161215666197044,-0.014273273413155542,0.8563964047893331,2500.0,,100.0,,1147.33333333,,"${\bf{x}}_{\beta \,{\bf{1}},1}$",91.5841034228,,31.0856523502,,Layer 2,1,l2_b,False,1,1,Default serie

input_data/tutorial_ch3_interfaces

@@ -0,0 +1,25 @@
+,X,X_std,Y,Y_std,Z,Z_std,annotations,formation,formation number,group_id,isFault,order_series,series
+0,250,0.0,0,0.0,996,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},0}$",Layer 2,1,l2_a,False,1,Default serie
+1,2500,0.0,200,0.0,1149,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},1}$",Layer 2,1,l2_b,False,1,Default serie
+2,2250,0.0,100,0.0,1298,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},2}$",Layer 2,1,l2_b,False,1,Default serie
+3,2750,0.0,0,0.0,995,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},3}$",Layer 2,1,l2_b,False,1,Default serie
+4,500,0.0,200,0.0,1149,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},4}$",Layer 2,1,l2_a,False,1,Default serie
+5,750,0.0,100,0.0,1299,0.0,"${\bf{x}}_{\alpha \,{\bf{1}},5}$",Layer 2,1,l2_a,False,1,Default serie
+6,750,0.0,100,0.0,699,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},0}$",Layer 5,2,l5_a,False,1,Default serie
+7,250,0.0,0,0.0,396,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},1}$",Layer 5,2,l5_a,False,1,Default serie
+8,2250,0.0,100,0.0,698,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},2}$",Layer 5,2,l5_b,False,1,Default serie
+9,2500,0.0,200,0.0,549,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},3}$",Layer 5,2,l5_b,False,1,Default serie
+10,500,0.0,200,0.0,549,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},4}$",Layer 5,2,l5_a,False,1,Default serie
+11,2750,0.0,0,0.0,395,0.0,"${\bf{x}}_{\alpha \,{\bf{2}},5}$",Layer 5,2,l5_b,False,1,Default serie
+12,250,0.0,0,0.0,796,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},0}$",Layer 3,3,l3_a,False,1,Default serie
+13,2750,0.0,0,0.0,795,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},1}$",Layer 3,3,l3_b,False,1,Default serie
+14,500,0.0,200,0.0,949,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},2}$",Layer 3,3,l3_a,False,1,Default serie
+15,750,0.0,100,0.0,1099,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},3}$",Layer 3,3,l3_a,False,1,Default serie
+16,2250,0.0,100,0.0,1098,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},4}$",Layer 3,3,l3_b,False,1,Default serie
+17,2500,0.0,200,0.0,949,0.0,"${\bf{x}}_{\alpha \,{\bf{3}},5}$",Layer 3,3,l3_b,False,1,Default serie
+18,2750,0.0,0,0.0,595,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},0}$",Layer 4,4,l4_b,False,1,Default serie
+19,250,0.0,0,0.0,596,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},1}$",Layer 4,4,l4_a,False,1,Default serie
+20,2500,0.0,200,0.0,749,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},2}$",Layer 4,4,l4_b,False,1,Default serie
+21,2250,0.0,100,0.0,898,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},3}$",Layer 4,4,l4_b,False,1,Default serie
+22,500,0.0,200,0.0,749,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},4}$",Layer 4,4,l4_a,False,1,Default serie
+23,750,0.0,100,0.0,899,0.0,"${\bf{x}}_{\alpha \,{\bf{4}},5}$",Layer 4,4,l4_a,False,1,Default serie

setup.py

@@ -2,7 +2,7 @@
 
 setup(
     name='gempy',
-    version='0.99',
+    version='0.991',
     packages=['gempy'],
     install_requires=[
         'numpy',
@@ -13,7 +13,7 @@
         'seaborn'
     ],
     url='https://github.com/cgre-aachen/gempy',
-    download_url='https://github.com/cgre-aachen/gempy/archive/0.99.tar.gz',
+    download_url='https://github.com/cgre-aachen/gempy/archive/0.991.tar.gz',
     license='MIT',
     author='Miguel de la Varga, Alexander Schaaf',
     author_email='[email protected]',