Merge pull request #54 from amoodie/longer_tests

Longer tests, logging, warnings

pyDeltaRCM/default.yml

@@ -1,9 +1,3 @@
-site_prefix:
-  type: 'str'
-  default: ''
-case_prefix:
-  type: 'str'
-  default: ''
 out_dir:
   type: 'str'
   default: 'deltaRCM_Output'

pyDeltaRCM/deltaRCM_tools.py

@@ -6,6 +6,7 @@
 import string
 import logging
 import time
+import warnings
 
 from math import floor, sqrt, pi
 import numpy as np
@@ -48,13 +49,16 @@ def run_one_timestep(self):
 
         timestep = self._time
 
+        self.logger.info('-' * 4 + ' Timestep ' +
+                         str(self._time) + ' ' + '-' * 4)
         if self.verbose > 0:
             print('-' * 20)
             print('Timestep: ' + str(self._time))
 
         if self._is_finalized:
             raise RuntimeError('Cannot update model, model already finalized!')
 
+        # model operations
         for iteration in range(self.itermax):
 
             self.init_water_iteration()
@@ -101,8 +105,10 @@ def expand_stratigraphy(self):
 
         """
 
-        if self.verbose:
-            self.logger.info('Expanding stratigraphy arrays')
+        _msg = 'Expanding stratigraphy arrays'
+        self.logger.info(_msg)
+        if self.verbose >= 2:
+            print(_msg)
 
         lil_blank = lil_matrix((self.L * self.W, self.n_steps),
                                dtype=np.float32)
@@ -143,8 +149,10 @@ def record_stratigraphy(self):
             if self.strata_eta.shape[1] <= timestep:
                 self.expand_stratigraphy()
 
+            _msg = 'Storing stratigraphy data'
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                self.logger.info('Storing stratigraphy data')
+                print(_msg)
 
             # ------------------ sand frac ------------------
             # -1 for cells with deposition volumes < vol_limit
@@ -211,8 +219,12 @@ def apply_subsidence(self):
             timestep = self._time
 
             if self.start_subsidence <= timestep:
+
+                _msg = 'Applying subsidence'
+                self.logger.info(_msg)
                 if self.verbose >= 2:
-                    self.logger.info('Applying subsidence')
+                    print(_msg)
+
                 self.eta[:] = self.eta - self.sigma
 
     def output_data(self):
@@ -252,34 +264,34 @@ def output_data(self):
                 plt.clim(self.clim_eta[0], self.clim_eta[1])
                 plt.colorbar()
                 plt.axis('equal')
-                self.save_figure(self.prefix + "eta_" + str(timestep))
+                self.save_figure(os.path.join(self.prefix, 'eta_' + str(timestep)))
 
             if self.save_stage_figs:
 
                 plt.pcolor(self.stage)
                 plt.colorbar()
                 plt.axis('equal')
-                self.save_figure(self.prefix + "stage_" + str(timestep))
+                self.save_figure(os.path.join(self.prefix, 'stage_' + str(timestep)))
 
             if self.save_depth_figs:
 
                 plt.pcolor(self.depth)
                 plt.colorbar()
                 plt.axis('equal')
-                self.save_figure(self.prefix + "depth_" + str(timestep))
+                self.save_figure(os.path.join(self.prefix, 'depth_' + str(timestep)))
 
             if self.save_discharge_figs:
 
                 plt.pcolor(self.qw)
                 plt.colorbar()
                 plt.axis('equal')
-                self.save_figure(self.prefix + "discharge_" + str(timestep))
+                self.save_figure(os.path.join(self.prefix, 'discharge_' + str(timestep)))
 
             if self.save_velocity_figs:
                 plt.pcolor(self.uw)
                 plt.colorbar()
                 plt.axis('equal')
-                self.save_figure(self.prefix + "velocity_" + str(timestep))
+                self.save_figure(os.path.join(self.prefix, 'velocity_' + str(timestep)))
 
             # ------------------ grids ------------------
             if self.save_eta_grids:
@@ -322,15 +334,18 @@ def output_strata(self):
 
         if self.save_strata:
 
+            _msg = 'Saving final stratigraphy to netCDF file'
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                self.logger.info('\nSaving final stratigraphy to netCDF file')
+                print(_msg)
 
             self.strata_eta = self.strata_eta[:, :self.strata_counter]
 
             shape = self.strata_eta.shape
             if shape[0] < 1:
                 raise RuntimeError('Stratigraphy are empty! '
-                                   'Are you sure you ran the model with `update()`?')
+                                   'Are you sure you ran the model at least '
+                                   'one timestep with `update()`?')
 
             total_strata_age = self.output_netcdf.createDimension(
                 'total_strata_age',

pyDeltaRCM/init_tools.py

@@ -4,7 +4,8 @@
 import re
 import string
 import logging
-import time
+# import time
+import warnings
 
 from math import floor, sqrt, pi
 import numpy as np
@@ -20,7 +21,7 @@
 import time as time_lib
 from scipy.sparse import lil_matrix, csc_matrix, hstack
 import logging
-import time
+# import time
 import yaml
 
 from . import shared_tools
@@ -30,22 +31,39 @@
 
 class init_tools(object):
 
+    def init_output_infrastructure(self):
+
+        # output directory config
+        self.prefix = self.out_dir
+        self.prefix_abspath = os.path.abspath(self.prefix)
+
+        # create directory if it does not exist
+        if not os.path.exists(self.prefix_abspath):
+            os.makedirs(self.prefix_abspath)
+            assert os.path.isdir(self.prefix_abspath)  # validate dir created
+
     def init_logger(self):
+        """Initialize a logger.
 
-        if self.verbose >= 1:
+        The logger is initialized regardless of the value of ``self.verbose``.
+        The level of information printed to the log depends on the verbosity
+        setting. 
 
-            self.logger = logging.getLogger("driver")
-            self.logger.setLevel(logging.INFO)
+        """
 
-            # create the logging file handler
-            st = timestr = time.strftime("%Y%m%d-%H%M%S")
-            fh = logging.FileHandler("pyDeltaRCM_" + st + ".log")
-            formatter = logging.Formatter(
-                '%(asctime)s - %(name)s - %(levelname)s - %(message)s')
-            fh.setFormatter(formatter)
+        self.logger = logging.getLogger('driver')
+        self.logger.setLevel(logging.INFO)
 
-            # add handler to logger object
-            self.logger.addHandler(fh)
+        # create the logging file handler
+        st = timestr = time_lib.strftime('%Y%m%d-%H%M%S')
+        fh = logging.FileHandler(
+            self.prefix_abspath + '/pyDeltaRCM_' + st + '.log')
+        formatter = logging.Formatter(
+            '%(asctime)s - %(levelname)s - %(message)s')
+        fh.setFormatter(formatter)
+
+        # add handler to logger object
+        self.logger.addHandler(fh)
 
     def import_files(self):
 
@@ -97,13 +115,29 @@ def import_files(self):
         for k, v in list(input_file_vars.items()):
             setattr(self, k, v)
 
+    def determine_random_seed(self):
+        """Set the random seed if given.
+
+        If a random seed is specified, set the seed to this value.
+
+        Writes the seed to the log for record.
+        """
         if self.seed is not None:
+
+            _msg = 'Setting random seed to: %s ' % str(self.seed)
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                print("setting random seed to %s " % str(self.seed))
+                print(_msg)
+
             shared_tools.set_random_seed(self.seed)
 
+        # always write the seed to file for record and reproducability
+        self.logger.info('Random seed is: %s ' % str(self.seed))
+
     def set_constants(self):
 
+        self.logger.info('Setting model constants')
+
         self.g = 9.81   # (gravitation const.)
 
         sqrt2 = np.sqrt(2)
@@ -217,16 +251,6 @@ def create_other_variables(self):
 
         self.diffusion_multiplier = (self.dt / self.N_crossdiff * self.alpha
                                      * 0.5 / self.dx**2)
-
-        # output directory config
-        self.prefix = self.out_dir
-
-        if self.out_dir[-1] != '/':
-            self.prefix = self.out_dir + '/'
-        if self.site_prefix:
-            self.prefix += self.site_prefix + '_'
-        if self.case_prefix:
-            self.prefix += self.case_prefix + '_'
 
         self._is_finalized = False
 
@@ -235,8 +259,9 @@ def create_domain(self):
         Creates the model domain
         """
 
-        # ---- empty arrays ----
+        self.logger.info('Creating model domain')
 
+        # ---- empty arrays ----
         self.x, self.y = np.meshgrid(
             np.arange(0, self.W), np.arange(0, self.L))
 
@@ -326,11 +351,11 @@ def init_stratigraphy(self):
             self.n_steps = 5 * self.save_dt
 
             self.strata_sand_frac = lil_matrix((self.L * self.W, self.n_steps),
-                                               dtype = np.float32)
+                                               dtype=np.float32)
 
             self.init_eta = self.eta.copy()
             self.strata_eta = lil_matrix((self.L * self.W, self.n_steps),
-                                         dtype = np.float32)
+                                         dtype=np.float32)
 
     def init_output_grids(self):
         """Creates a netCDF file to store output grids.
@@ -348,22 +373,20 @@ def init_output_grids(self):
                 self.save_velocity_grids or
                 self.save_strata):
 
+            _msg = 'Generating netCDF file for output grids'
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                self.logger.info('Generating netCDF file for output grids...')
+                print(_msg)
 
             directory = self.prefix
             filename = 'pyDeltaRCM_output.nc'
 
-            if not os.path.exists(directory):
-                if self.verbose >= 2:
-                    self.logger.info('Creating output directory')
-                os.makedirs(directory)
-
             file_path = os.path.join(directory, filename)
-
             if os.path.exists(file_path):
+                _msg = 'Replacing existing netCDF file'
+                self.logger.warning(_msg)
                 if self.verbose >= 2:
-                    self.logger.info('*** Replaced existing netCDF file ***')
+                    warnings.warn(UserWarning(_msg))
                 os.remove(file_path)
 
             self.output_netcdf = Dataset(file_path, 'w',
@@ -421,8 +444,10 @@ def init_output_grids(self):
                                                              ('total_time', 'length', 'width'))
                 velocity.units = 'meters per second'
 
+            _msg = 'Output netCDF file created'
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                self.logger.info('Output netCDF file created.')
+                print(_msg)
 
     def init_subsidence(self):
         """

pyDeltaRCM/model.py

@@ -56,16 +56,20 @@ def __init__(self, input_file=None):
         self.default_file = os.path.join(_src_dir, 'default.yml')
         self.import_files()
 
-        self.create_other_variables()
-
+        self.init_output_infrastructure()
         self.init_logger()
 
+        self.create_other_variables()
+
+        self.determine_random_seed()
         self.create_domain()
 
         self.init_subsidence()
         self.init_stratigraphy()
         self.init_output_grids()
 
+        self.logger.info('Model initialization complete')
+
     def update(self):
         """Run the model for one full instance
 
@@ -113,12 +117,16 @@ def finalize(self):
 
         """
 
+        self.logger.info('Finalize model run')
+
         self.output_strata()
 
         try:
             self.output_netcdf.close()
-            if self.verbose >= 1:
-                print('Closed output netcdf file.')
+            _msg = 'Closed output netcdf file'
+            self.logger.info(_msg)
+            if self.verbose >= 2:
+                print(_msg)
         except Exception:
             pass
 

pyDeltaRCM/water_tools.py

@@ -194,9 +194,10 @@ def check_size_of_indices_matrix(self, it):
             Once it reaches it > self.itmax/2 once, make the size
             self.iter for all further timesteps
             """
-
+            _msg = 'Increasing size of self.indices'
+            self.logger.info(_msg)
             if self.verbose >= 2:
-                self.logger.info('Increasing size of self.indices')
+                print(_msg)
 
             indices_blank = np.zeros(
                 (np.int(self.Np_water), np.int(self.itmax / 4)), dtype=np.int)

run_pyDeltaRCM.py

@@ -3,7 +3,7 @@
 
 if __name__ == '__main__':
 
-    delta = DeltaModel(input_file = os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+    delta = DeltaModel()
 
     for time in range(0, 1):
         delta.update()