Merge branch 'develop' into richard2

docs/Makefile

@@ -31,6 +31,6 @@ test : clean
 
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
-%: clean Makefile 
+%: Makefile 
 	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
 

docs/source/examples/basic_runs.rst

@@ -0,0 +1,29 @@
+Basic script to run the model
+=============================
+
+.. plot::
+    :context: reset
+
+    with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+        delta = pyDeltaRCM.DeltaModel(
+            out_dir=output_dir,
+            resume_checkpoint='../_resources/checkpoint')
+        delta.finalize()
+
+
+.. code::
+
+    delta = pyDeltaRCM.DeltaModel()
+
+    for _t in range(0, 1000):
+        delta.update()
+
+    delta.finalize()
+
+.. plot::
+    :context:
+    :include-source:
+
+    fig, ax = plt.subplots()
+    ax.imshow(delta.eta)
+    plt.show()

docs/source/examples/custom_class_preprocessor.rst

@@ -0,0 +1,104 @@
+Using custom classes with the Preprocessor
+==========================================
+
+Here, we set up three jobs to run as an ensemble of a single custom class.
+
+
+.. plot::
+    :context: reset
+    :include-source:
+
+    >>> import numpy as np
+    >>> import matplotlib.pyplot as plt
+
+    >>> import pyDeltaRCM
+
+    >>> class CustomRandomModel(pyDeltaRCM.DeltaModel):
+    ...     """
+    ...     A subclass of DeltaModel that runs fast for this example.
+    ... 
+    ...     Just for the sake of this example, we are implementing a 
+    ...     custom class that runs very quickly. We override the 
+    ...     `solve_water_and_sediment_timestep` method of the model to 
+    ...     simply add random gaussian blobs to the surface on each step.
+    ...     """
+    ...     def __init__(self, input_file=None, **kwargs):
+    ...     
+    ...         # inherit base DeltaModel methods
+    ...         super().__init__(input_file, **kwargs)
+    ...     
+    ...         self.noise_patch = int(25)
+    ...         self.noise_size = 5     # x y scale 
+    ...         self.noise_scale = 200  # z scale
+    ...         self._half_ns = self.noise_patch // 2
+    ...         self.noise_x, self.noise_y = np.meshgrid(
+    ...             np.linspace(-self._half_ns, self._half_ns, num=self.noise_patch),
+    ...             np.linspace(-self._half_ns, self._half_ns, num=self.noise_patch))
+    ...     
+    ...     def solve_water_and_sediment_timestep(self):
+    ...         """Overwrite method for documentation demonstration.
+    ...     
+    ...         This method now simply adds random gaussian noise on each step.
+    ...         """         
+    ...         # get a random x and y value
+    ...         #   important: use get_random_uniform for reproducibility!
+    ...         x, y, z = [pyDeltaRCM.shared_tools.get_random_uniform(1) for _ in range(3)]
+    ...         
+    ...         # rescale to fit inside domain
+    ...         x = int(x * (self.L - self.noise_patch))
+    ...         y = int(y * (self.W - self.noise_patch))
+    ...     
+    ...         # generate the blob
+    ...         sx = sy = self.noise_size
+    ...         exp = np.exp(-((self.noise_x)**2. / (2. * sx**2.) + (self.noise_y)**2. / (2. * sy**2.)))
+    ...         blob = (1. / (2. * np.pi * sx * sy) * exp * self.noise_scale)
+    ...         
+    ...         # place into domain
+    ...         self.eta[x:x+self.noise_patch, y:y+self.noise_patch] += blob
+
+
+Then, we pass this custom subclass to the :obj:`~pyDeltaRCM.preprocessor.Preprocessor.run_jobs` method.
+
+.. plot::
+    :context:
+
+    with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+        param_dict = dict(
+            out_dir=output_dir,
+            ensemble=3,
+            timesteps=50
+            )
+
+        # let the preprocessor set up the jobs for you from checkpoint
+        pp = pyDeltaRCM.Preprocessor(
+            param_dict)
+
+        # run the jobs
+        pp.run_jobs(DeltaModel=CustomRandomModel)
+
+.. code::
+
+    >>> # set up dictionary for parameters and create a `Preprocessor`
+    >>> param_dict = dict(
+    ...     ensemble=3,
+    ...     timesteps=50
+    ...     )
+
+    >>> # preprocessor set up the jobs
+    >>> pp = pyDeltaRCM.Preprocessor(
+    ...     param_dict)
+
+    >>> # run the jobs with custom class!
+    >>> pp.run_jobs(DeltaModel=CustomRandomModel)
+
+.. plot::
+    :context:
+    :include-source:
+
+    >>> fig, ax = plt.subplots(
+    ...     1, len(pp.job_list),
+    ...     figsize=(12, 4.8))
+    >>> for i in range(len(pp.job_list)):
+    ...     ax[i].imshow(pp.job_list[i].deltamodel.eta)
+    >>> plt.tight_layout()
+    >>> plt.show()

docs/source/examples/index.rst

@@ -2,6 +2,15 @@
 Examples of using and working with pyDeltaRCM
 =============================================
 
+Running the model
+-----------------
+
+.. toctree::
+   :maxdepth: 1
+
+   basic_runs
+   resume_from_checkpoint
+
 
 Modifying initial conditions
 ----------------------------
@@ -27,13 +36,11 @@ Modifying boundary conditions
 Modifying internal computations
 -------------------------------
 
-* none
-
+.. toctree::
+   :maxdepth: 1
 
-Modifying behavior of the Preprocessor
---------------------------------------
+   overwrite_topo_diffusion
 
-* none
 
 Modifying Model Input/Output
 ----------------------------
@@ -43,3 +50,12 @@ Modifying Model Input/Output
 
    custom_yaml
    custom_saving
+
+
+Working with the Preprocessor
+-----------------------------
+
+.. toctree::
+   :maxdepth: 1
+
+   custom_class_preprocessor

docs/source/examples/overwrite_topo_diffusion.rst

@@ -0,0 +1,154 @@
+Changing topographic diffusion to represent tree throw
+======================================================
+
+Here, we demonstrate how to overwrite an existing method of the `DeltaModel` to achieve custom behavior during model runtime.
+
+.. note::
+
+    This example demonstrates several best-practices, including using yaml parameters specifications, custom figure and grid saving setup, and using :obj:`~pyDeltaRCM.shared_tools.get_random_uniform` for random number generation.
+
+In implementing custom model subclasses, it is common to want to change runtime behavior of the model. 
+This can often be achieved by using hooks alone, but sometimes a combination of hooks and overwriting existing methods is necessary.
+
+In this example, we calculate a diffusion multiplier to represent tree throw. 
+In this super simple and **likely way over-exaggerated** representation of this process, we assume:
+
+* there are trees everywhere the elevation of a cell has been above `self.dry_depth` for 10 consecutive timesteps
+* there is a probability threshold of tree throw occurring anywhere trees exist
+* tree throw makes the diffusion multiplier for that cell on that timestep really big!
+
+.. important::
+
+    There are all kinds of problems with the assumptions in this example. Don't read into it too much. It's an example to show how to modify code, not how to represent tree throw.
+
+.. plot::
+    :context: reset
+    :include-source:
+
+    class TreeThrowModel(pyDeltaRCM.DeltaModel):
+        """Implementation of tree throw.
+        """
+
+        def __init__(self, input_file=None, **kwargs):
+
+            # inherit from base model
+            super().__init__(input_file, **kwargs)
+            self.hook_after_create_domain()
+
+        def hook_import_files(self):
+            """Define the custom YAML parameters."""
+            # whether to run vegetation
+            self.subclass_parameters['tree_throw'] = {
+                'type': 'bool', 'default': False
+            }
+
+            # tree throw multiplier
+            self.subclass_parameters['p_tt_mult'] = {
+                'type': ['int', 'float'], 'default': 100
+            }
+
+            # tree throw establish timesteps
+            self.subclass_parameters['p_tt_iter'] = {
+                'type': ['int', 'float'], 'default': 10
+            }
+
+            # tree throw prob threshold
+            self.subclass_parameters['p_tt_prob'] = {
+                'type': ['int', 'float'], 'default': 0.2
+            }
+
+        def hook_init_output_file(self):
+            """Add non-standard grids, figures and metadata to be saved."""
+            # save a figure of the active layer each save_dt
+            self._save_fig_list['tree'] = ['tree']
+
+            # save the active layer grid each save_dt w/ a short name
+            self._save_var_list['tree'] = ['tree', 'boolean',
+                                           'i4', ('time',
+                                                  'x', 'y')]
+
+        def hook_after_create_domain(self):
+            """Add fields to the model for all tree parameterizations.
+            """
+            self.tree = np.zeros(self.depth.shape, dtype=np.int64)
+            self.dry_count = np.zeros(self.depth.shape, dtype=np.int64)
+
+            self.tree_multiplier = np.ones_like(self.depth)
+
+        def hook_after_route_sediment(self):
+            """Apply vegetation growth/death rules.
+            """
+            # determine cells dry and increment counter
+            _where_dry = (self.depth < self.dry_depth)
+            self.dry_count[~_where_dry] = 0  # any wet gets reset
+            self.dry_count[_where_dry] += 1  # any dry gets incremented
+
+            # if tree_throw is on, run the tree placing routine
+            if self.tree_throw:
+
+                # trees die anywhere wet
+                self.tree[~_where_dry] = int(0)
+
+                # trees go anywhere dry for more than threshold
+                _where_above_thresh = (self.dry_count >= self.p_tt_iter)
+
+                self.tree[_where_above_thresh] = int(1)
+
+                # determine the multiplier field
+                _rand = np.array([get_random_uniform(1) for i in np.arange(self.depth.size)]).reshape(self.depth.shape)
+                _thrown = np.logical_and((_rand < self.p_tt_prob), self.tree)
+
+                # ignore the strip of land
+                _thrown[self.cell_type == -2] = 0
+
+                # set to ones everywhere, then overwrite with multiplier
+                self.tree_multiplier[:] = 1
+                self.tree_multiplier[_thrown] = self.p_tt_mult
+
+        def topo_diffusion(self):
+            """Overwrite with new behavior.
+
+            This method is very similar to the base DeltaModel code, but we add an
+            additional multiplier to represent tree throw.
+            """
+            for _ in range(self.N_crossdiff):
+
+                a = ndimage.convolve(self.eta, self.kernel1, mode='constant')
+                b = ndimage.convolve(self.qs, self.kernel2, mode='constant')
+                c = ndimage.convolve(self.qs * self.eta, self.kernel2,
+                                     mode='constant')
+
+                self.cf = (self.tree_multiplier * self.diffusion_multiplier *
+                           (self.qs * a - self.eta * b + c))
+
+                self.cf[self.cell_type == -2] = 0
+                self.cf[0, :] = 0
+
+                self.eta += self.cf
+
+And the model is then instantiated with:
+
+.. plot::
+    :context:
+
+    with pyDeltaRCM.shared_tools._docs_temp_directory() as output_dir:
+        mdl = TreeThrowModel(
+            out_dir=output_dir,
+            tree_throw=True)
+
+.. code:: python
+
+    mdl = TreeThrowModel(
+        tree_throw=True)
+
+We don't actually run this model at all in this example.
+Let's plot the ``.tree`` field just to see that the subclass was instantiated correctly.
+
+.. plot::
+    :context:
+    :include-source:
+
+    fig, ax = plt.subplots()
+    im = ax.imshow(mdl.tree)
+    plt.colorbar(im, ax=ax, shrink=0.5)
+    plt.show()