(FEATURE) added logging to discrete interpolator

for DOF and can now specify element_volume when adding a feature

Author: lachlangrose

LoopStructural/interpolators/discrete_interpolator.py

@@ -52,7 +52,8 @@ def __init__(self, support):
         self.eq_const_c_ = 0
         self.constraints = {}
         self.interpolation_weights= {}
-
+        logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.nx))
+        
     def set_property_name(self, propertyname):
         """
         Set the property name attribute, this is usually used to
@@ -89,7 +90,8 @@ def set_region(self, region=None):
         self.region_map[self.region] = np.array(
             range(0, len(self.region_map[self.region])))
         self.nx = len(self.support.nodes[self.region])
-        
+        logger.info("Interpolation now uses region and has {} degrees of freedom".format(self.nx))
+
     def set_interpolation_weights(self, weights):
         """
         Set the interpolation weights dictionary

LoopStructural/modelling/core/geological_model.py

@@ -36,7 +36,7 @@
 from LoopStructural.utils.helper import (all_heading, gradient_vec_names,
                                          strike_dip_vector)
 
-from LoopStructural.utils import getLogger
+from LoopStructural.utils import getLogger, log_to_file
 logger = getLogger(__name__)
 if not surfe:
     logger.warning("Cannot import Surfe")
@@ -86,7 +86,7 @@ class GeologicalModel:
 
     """
     def __init__(self, origin, maximum, rescale=True, nsteps=(40, 40, 40),
-                 reuse_supports=False):
+                 reuse_supports=False, logfile=None, loglevel='info'):
         """
         Parameters
         ----------
@@ -120,6 +120,11 @@ def __init__(self, origin, maximum, rescale=True, nsteps=(40, 40, 40),
 
 
         """
+        if logfile:
+            self.logfile = logfile
+            log_to_file(logfile,loglevel)
+        
+        logger.info('Initialising geological model')
         self.features = []
         self.feature_name_index = {}
         self.data = None
@@ -128,19 +133,23 @@ def __init__(self, origin, maximum, rescale=True, nsteps=(40, 40, 40),
         # we want to rescale the model area so that the maximum length is
         # 1
         self.origin = np.array(origin).astype(float)
+        logger.info('Model origin: {} {} {}'.format(self.origin[0],self.origin[1],self.origin[2]))
 
         self.maximum = np.array(maximum).astype(float)
+        logger.info('Model maximum: {} {} {}'.format(self.maximum[0],self.maximum[1],self.maximum[2]))
         lengths = self.maximum - self.origin
         self.scale_factor = 1.
         self.bounding_box = np.zeros((2, 3))
         self.bounding_box[1, :] = self.maximum - self.origin
         self.bounding_box[1, :] = self.maximum - self.origin
         if rescale:
             self.scale_factor = np.max(lengths)
+            logger.info('Rescaling model using scale factor {}'.format(self.scale_factor))
 
         self.bounding_box /= self.scale_factor
         self.support = {}
         self.reuse_supports = reuse_supports
+        logger.info('Reusing interpolation supports: {}'.format(self.reuse_supports))
         self.stratigraphic_column = None
         self.parameters = {'features': [], 'model': {'bounding_box': self.origin.tolist() + self.maximum.tolist(),
                                                      'rescale': rescale,
@@ -167,6 +176,7 @@ def from_map2loop_directory(cls, m2l_directory,**kwargs):
             the created geological model and a dictionary of the map2loop data
         """
         from LoopStructural.utils import build_model, process_map2loop
+        logger.info('LoopStructural model initialised from m2l directory: {}'.format(m2l_directory))
         m2lflags = kwargs.pop('m2lflags',{})
         m2l_data = process_map2loop(m2l_directory,m2lflags)
         return build_model(m2l_data,**kwargs), m2l_data
@@ -192,11 +202,12 @@ def from_file(cls, file):
             return None
         model = pickle.load(open(file,'rb'))
         if type(model) == GeologicalModel:
+            logger.info('GeologicalModel initialised from file')
             return model
         else:
             logger.error('{} does not contain a geological model'.format(file))
             return None
-
+        
     def __getitem__(self, feature_name):
         """Accessor for feature in features using feature_name_index
 
@@ -227,6 +238,7 @@ def to_file(self, file):
                         "pip install dill")
             return
         try:
+            logger.info('Writing GeologicalModel to: {}'.format(file))
             pickle.dump(self,open(file,'wb'))
         except pickle.PicklingError:
             logger.error('Error saving file')
@@ -285,7 +297,7 @@ def set_model_data(self, data):
                 data = pd.read_csv(data)
             except:
                 logger.error("Could not load pandas data frame from data")
-
+        logger.info('Adding data to GeologicalModel with {} data points'.format(len(data)))
         self.data = data.copy()
         self.data['X'] -= self.origin[0]
         self.data['Y'] -= self.origin[1]
@@ -294,7 +306,7 @@ def set_model_data(self, data):
         self.data['Y'] /= self.scale_factor
         self.data['Z'] /= self.scale_factor
         if 'type' in self.data:
-            logger.warning("'type' is being replaced with 'feature_name' \n")
+            logger.warning("'type' is depreciated replace with 'feature_name' \n")
             self.data.rename(columns={'type':'feature_name'},inplace=True)
         for h in all_heading():
             if h not in self.data:
@@ -305,6 +317,7 @@ def set_model_data(self, data):
                     self.data[h] = 0
 
         if 'strike' in self.data and 'dip' in self.data:
+            logger.info('Converting strike and dip to vectors')
             mask = np.all(~np.isnan(self.data.loc[:, ['strike', 'dip']]),
                           axis=1)
             self.data.loc[mask, gradient_vec_names()] = strike_dip_vector(
@@ -359,6 +372,7 @@ def set_stratigraphic_column(self, stratigraphic_column,cmap='tab20'):
         """
         # if the colour for a unit hasn't been specified we can just sample from 
         # a colour map e.g. tab20
+        logger.info('Adding stratigraphic column to model')
         random_colour = True
         n_units=0
         for g in stratigraphic_column.keys():
@@ -403,7 +417,7 @@ def create_from_feature_list(self, features):
                 self.create_and_add_folded_foliation(f)
 
     def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
-                         buffer=0.2, **kwargs):
+                         buffer=0.2, element_volume = None, **kwargs):
         """
         Returns an interpolator given the arguments, also constructs a
         support for a discrete interpolator
@@ -435,11 +449,12 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
         bb[1, :] += buffer  # *(bb[1,:]-bb[0,:])
         box_vol = (bb[1, 0]-bb[0, 0]) * (bb[1, 1]-bb[0, 1]) * (bb[1, 2]-bb[0, 2])
         if interpolatortype == "PLI":
-            nelements /= 5
-            ele_vol = box_vol / nelements
+            if element_volume is None:
+                nelements /= 5
+                element_volume = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
-            step_vector[:] = ele_vol ** (1. / 3.)
+            step_vector[:] = element_volume ** (1. / 3.)
             # step_vector /= np.array([1,1,2])
             # number of steps is the length of the box / step vector
             nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
@@ -457,10 +472,11 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
 
         if interpolatortype == 'FDI':
             # find the volume of one element
-            ele_vol = box_vol / nelements
+            if element_volume is None:
+                element_volume = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
-            step_vector[:] = ele_vol ** (1. / 3.)
+            step_vector[:] = element_volume ** (1. / 3.)
             # number of steps is the length of the box / step vector
             nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             if np.any(np.less(nsteps, 3)):
@@ -478,11 +494,12 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
             return FDI(grid)
 
         if interpolatortype == "DFI":  # "fold" in kwargs:
-            nelements /= 5
-            ele_vol = box_vol / nelements
+            if element_volume is None:
+                nelements /= 5
+                element_volume = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
-            step_vector[:] = ele_vol ** (1. / 3.)
+            step_vector[:] = element_volume ** (1. / 3.)
             # number of steps is the length of the box / step vector
             nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             # create a structured grid using the origin and number of steps