clean up for pub

formatting in pipeline.py, grow_sim.py
cleaned up plotting notebook

pyutils/__init__.py

@@ -1,5 +1,5 @@
 # ====================================================================================== #
-# Module for forest data and simulations.
+# Module for sessile organism simulations.
 # Author : Eddie Lee, [email protected]
 # 
 #

pyutils/data.py

@@ -1,5 +1,5 @@
 # ====================================================================================== #
-# Useful routines for analysis of forest data.
+# Useful routines for analysis of observational data and simulation results.
 # Author : Eddie Lee, [email protected]
 # 
 #
@@ -25,13 +25,13 @@
 #     OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 #     SOFTWARE.
 # ====================================================================================== #
-from .utils import *
-from . import nearest_neighbor as nn
-import pandas as pd
 import matplotlib.pyplot as plt
 from matplotlib.patches import Circle
 from matplotlib.collections import PatchCollection
 
+from .utils import *
+from . import nearest_neighbor as nn
+
 
 
 def plot(xy, rbh, L,
@@ -92,12 +92,16 @@ def plot(xy, rbh, L,
         return fig
 
 def namibia_corr_fcn():
+    """Calculate correlation function for Namibian termate data set from Tarnita et al.
+    """
+
+    # load data
     df = pd.read_csv('../data/Tarnita/termite_mound_location_field_data/Namib_G1.txt',
                      sep='\t',
                      header=None)
     xy = df.values
 
-    # exclude out 100 meters from boundary of plot
+    # exclude area 100 meters in from boundaries of rectangular plot
     p, r = nn.pair_correlation(xy, np.linspace(0, 10, 50), (100, 100, 400, 400))
 
     return p, r

pyutils/grow_sim.py

@@ -1,5 +1,5 @@
 # ====================================================================================== #
-# Automata compartment model for forest growth.
+# Automata compartment model for sessile organism growth based on forests.
 # Author : Eddie Lee, [email protected]
 # 
 #
@@ -32,6 +32,7 @@
 from warnings import warn
 from misc.stats import PowerLaw
 from types import LambdaType
+
 from .utils import *
 
 
@@ -332,6 +333,7 @@ def sample(self, n_sample,
         return_trees : bool, False
         n_cpus : int, None
         **kwargs
+            These go into self.grow().
         
         Returns
         -------
@@ -341,6 +343,8 @@ def sample(self, n_sample,
             Time.
         ndarray
             Compartments r_k.
+        list of list of Tree
+            For each forest, i.e. outermost list length is given by n_forests.
         """
 
         if n_forests==1:

pyutils/nearest_neighbor.py

@@ -25,9 +25,10 @@
 #     OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 #     SOFTWARE.
 # ====================================================================================== #
-from .utils import *
 from warnings import warn
 
+from .utils import *
+
 
 
 def dist(xy):

pyutils/pipeline.py

@@ -1,5 +1,5 @@
 # ====================================================================================== #
-# Pipeline tools for analysis.
+# Pipeline tools for generating publication results.
 # Author : Eddie Lee, [email protected]
 # 
 #
@@ -25,11 +25,12 @@
 #     OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 #     SOFTWARE.
 # ====================================================================================== #
-from .utils import *
-from .grow_sim import Forest2D
 import pandas as pd
 from workspace.utils import save_pickle
 
+from .utils import *
+from .grow_sim import Forest2D
+
 
 
 def WEB_transience():
@@ -80,27 +81,29 @@ def mft_cutoff(nu=2.5, nForests=30):
     Parameters
     ----------
     nu : float, 2.5
+        Fatal fluctuation exponent.
     nForests : int, 30
         Number of indpt forests to run.
     """
 
-    basalRange = np.array([.8, .2, .05, .0125, .003125])
+    basalRange = np.array([.8, .2, .05, .0125, .003125])  # basal metabolism coeff
 
     # set up
-    r0 = 1  # smallest stem radius
+    r0 = 1     # smallest stem radius
     Abar = .5  # death rate coeff
-    cg = .3  # growth rate coeff
+    cg = .3    # growth rate coeff
 
     rRange = np.linspace(r0, 400, 2000)
     g0 = 500  # sapling introduction rate
-    L = 200  # length of boundary
+    L = 200   # length of boundary
     sampleSize = 5000  # number of samples to take
-    dt = .1  # sim time step
-
-    nk = {}
-    t = {}
-    rk = {}
+    dt = .1   # sim time step
 
+    nk = {}  # no. of trees per size class k (pop. no.)
+    t = {}   # time
+    rk = {}  # radius of size class k
+
+    # run forest sim over basal metabolism coeff
     for basal in basalRange:
         forest = Forest2D(L, g0, rRange, 
                           {'root':1,
@@ -124,8 +127,10 @@ def mft_cutoff(nu=2.5, nForests=30):
 
 def mft_cutoff_finite_size_checks(nu=2.5, run_smaller=True, run_larger=True):
     """Simulation of symmetric competition with changing cutoff modulated by adjusting
-    basal metabolic rate coefficient. This decreases and increases area generated by
-    .mft_cutoff() by a factor of 4 to test for finite size effects.
+    basal metabolic rate coefficient.
+
+    This allows for a factor of 16 difference in area generated by .mft_cutoff() to test
+    for finite size effects.
 
     Parameters
     ----------
@@ -134,25 +139,25 @@ def mft_cutoff_finite_size_checks(nu=2.5, run_smaller=True, run_larger=True):
     run_larger : bool, True
     """
 
-    basalRange = np.array([.8, .2, .05, .0125, .003125])
+    basalRange = np.array([.8, .2, .05, .0125, .003125])  # basal met coeff
 
     # set up
-    r0 = 1
-    Abar = .5
-    cg = .3
+    r0 = 1     # sapling radius
+    Abar = .5  # natural mortality coeff
+    cg = .3    # growth coeff
 
-    rRange = np.linspace(r0, 400, 2000)
-    sampleSize = 5000
-    dt = .1
-    nForests = 30
+    rRange = np.linspace(r0, 400, 2000)  # radii of size classes
+    sampleSize = 5000  # no. of samples to take
+    dt = .1            # time step size
+    nForests = 30      # no. of random forests to sim.
 
     # smaller system
     if run_smaller:
-        g0 = 500 / 4
-        L = 200 / 2
-        nk = {}
-        t = {}
-        rk = {}
+        g0 = 500 / 4  # incoming sapling rate
+        L = 200 / 2   # system length
+        nk = {}       # pop no. by size class k
+        t = {}        # time
+        rk = {}       # radiuso f size class k
 
         for basal in basalRange:
             forest = Forest2D(L, g0, rRange, 
@@ -213,20 +218,21 @@ def mft_cutoff_plot(nu=2.5):
     Parameters
     ----------
     nu : float, 2.5
+        Fatal fluctuation exponent.
     """
 
     basalRange = np.array([.8, .2, .05, .0125, .003125])
 
     # set up
-    r0 = 1  # smallest stem radius
+    r0 = 1     # smallest stem radius
     Abar = .5  # death rate coeff
-    cg = .3  # growth rate coeff
+    cg = .3    # growth rate coeff
 
     rRange = np.linspace(r0, 400, 2000)
     g0 = 500  # sapling introduction rate
-    L = 200  # length of boundary
+    L = 200   # length of boundary
     sampleSize = 4000  # number of samples to take
-    dt = .1  # sim time step
+    dt = .1   # sim time step
 
     def loop_wrapper(basal):
         forest = Forest2D(L, g0, rRange, 
@@ -253,17 +259,17 @@ def loop_wrapper(basal):
                 f'plotting/biomass_scaling_w_compet_{nu=}.p', True)
 
 def phase_space_scan_Abar():
-    """Scanning across natural mortality rate Abar as in Figure 3.
+    """Scanning across natural mortality rate Abar as in Figure 4.
     """
 
-    AbarRange = np.linspace(.75, 0, 5)  # keys to xy dict
-    areaDeathRateRange = np.logspace(-1, 2, 10)  # keys to dicts in xy
+    AbarRange = np.linspace(.75, 0, 5)           # death rate coeff
+    areaDeathRateRange = np.logspace(-1, 2, 10)  # comp attrition rate coeff
 
     # set up
-    r0 = 1
-    cg = .3
-    nu = 2.
-    basal = .05
+    r0 = 1       # sapling basal stem radius
+    cg = .3      # growth coeff
+    nu = 2.      # fatal fluc exponent
+    basal = .05  # basal met coeff
 
     rRange = np.linspace(r0, 800, 1600)  # growth saturates 
     g0 = 100
@@ -278,10 +284,12 @@ def phase_space_scan_Abar():
               'basal':basal,
               'sharing fraction':1,
               'resource efficiency':2}
-
+
+    # loop over natural mortality rate
     def loop_Abar(Abar):
         coeffs['death'] = Abar
-
+
+        # loop over competitive death rate
         def loop_wrapper(deathRate):
             coeffs['dep death rate'] = deathRate
             forest = Forest2D(L, g0, rRange, coeffs,
@@ -312,7 +320,7 @@ def loop_wrapper(deathRate):
 
         return xy, nk
 
-    xy = {}  # loop over mortality rates
+    xy = {}  # tree xy coord, indexed by Abar and areadeathrate
     nk = {}  # pop. number (can be used for equilibrium check)
     for Abar in AbarRange:
         xy[Abar], nk[Abar] = loop_Abar(Abar)
@@ -324,19 +332,19 @@ def loop_wrapper(deathRate):
 
 def phase_space_scan_abar():
     """Scanning across varying growth rate fixing natural mortality rate to 0 as in Figure
-    3.
+    4.
     """
 
     # for showing the spatial distributions
-    cgRange = np.logspace(np.log10(.5), -4, 4)
-    areaDeathRateRange = np.logspace(-1, 2, 10)  # keys to dicts in xy
+    cgRange = np.logspace(np.log10(.5), -4, 4)   # growth rate coeff
+    areaDeathRateRange = np.logspace(-1, 2, 10)  # comp attrition rate coeff
 
     # set up
-    r0 = 1
-    Abar = 0.
-    basal = .05
+    r0 = 1       # sapling radius
+    Abar = 0.    # death rate coeff
+    basal = .05  # basal met rate coeff
 
-    rRange = np.linspace(r0, 400, 800)  # growth saturates 
+    rRange = np.linspace(r0, 400, 800)  # growth saturates at radius=400
     g0 = 100
     L = 200
     burnIn = 400
@@ -348,7 +356,8 @@ def phase_space_scan_abar():
               'basal':basal,
               'sharing fraction':1,
               'resource efficiency':2}
-
+
+    # loop over sim parameters
     def loop_cg(cg):
         coeffs['grow'] = cg
 
@@ -395,17 +404,17 @@ def hex_packing():
 
     from .nearest_neighbor import pair_correlation
 
-    # for showing the spatial distributions
+    # this section for showing the spatial distributions
     areaDeathRateRange = np.logspace(-1, 3, 10)  # keys to dicts in xy
 
     # set up
     r0 = 1
     basal = 0
 
-    rRange = np.linspace(r0, 5, 5)  # growth saturates b/c max radius is 5
-    g0 = 100
-    L = 200
-    burnIn = 1_000  # in time steps
+    rRange = np.linspace(r0, 5, 5)  # growth saturates at max radius of 5
+    g0 = 100  # incoming sapling rate
+    L = 200   # system length
+    burnIn = 1_000  # time steps to ignore
     sampleSize = 1_000
     dt = .2
     coeffs = {'root':10,
@@ -441,7 +450,7 @@ def loop_wrapper(deathRate):
                  'g0','L','burnIn','sampleSize','dt','coeffs','xy','nk'],
                 'cache/packing_example.p', True)
 
-    # for plotting the correlation fcn
+    # this section for plotting the correlation fcn
     allxy = xy
     p = {}
     bins = np.linspace(0, 5, 40)  # this should be roughly aligned with the stats of the system

pyutils/utils.py

@@ -1,5 +1,5 @@
 # ====================================================================================== #
-# Forest analysis.
+# Useful functions for sessile package.
 # Author : Eddie Lee, [email protected]
 # 
 #
@@ -26,6 +26,8 @@
 #     SOFTWARE.
 # ====================================================================================== #
 import numpy as np
+import pandas as pd
+import os
 from netCDF4 import Dataset
 from scipy.spatial.distance import pdist, squareform
 from numpy import ma