Fastslow initval

pyrealm/pmodel/scaler.py

@@ -407,6 +407,7 @@ def get_daily_means(
     def fill_daily_to_subdaily(
         self,
         values: NDArray,
+        previous_value: NDArray | None = None,
         update_point: str = "max",
         kind: str = "previous",
         fill_from: np.timedelta64 | None = None,
@@ -419,45 +420,53 @@ def fill_daily_to_subdaily(
         axis of the `values` must be the same length as the number of days used to
         create the instance.
 
+        The update point defaults to the maximum time of day during the acclimation
+        window. It can also be set to the mean time of day, but note that this implies
+        that the plant predicts the daily values between the mean and max observation
+        time. The ``fill_from`` argument can be used to set the update point to an
+        arbitrary time of day.
+
         Two interpolation kinds are currently implemented:
 
         * ``previous`` interpolates the daily value as a constant, until updating to the
           next daily value. This option will fill values until the end of the time
-          series.
+          series. The
         * ``linear`` interpolates linearly between the update points of the daily
           values. The interpolated values are held constant for the first day and then
           interpolated linearly: this is to avoid plants adapting optimally to future
           conditions.
 
-        The update point defaults to the maximum time of day during the acclimation
-        window. It can also be set to the mean time of day, but note that this implies
-        that the plant predicts the daily values between the mean and max observation
-        time. The ``fill_from`` argument can be used to set the update point to an
-        arbitrary time of day.
+        Subdaily observations before the update point on the first day of the time
+        series are filled with ``np.nan``. The ``previous_value`` argument can be used
+        to provide an alternative value, allowing time series to be processed in blocks,
+        but this option is only currently implemented for the ``previous``
+        interpolation method.
 
         Args:
             values: An array with the first dimension matching the number of days in the
-              instances :class:`~pyrealm.pmodel.scaler.SubdailyScaler` object.
+                instances :class:`~pyrealm.pmodel.scaler.SubdailyScaler` object.
+            previous_value: An array with dimensions equal to a slice across the first
+                axis of the values array.
             update_point: The point in the acclimation window at which the plant updates
-              to the new daily value: one of 'mean' or 'max'.
+                to the new daily value: one of 'mean' or 'max'.
             kind: The kind of interpolation to use to fill between daily values: one of
-              'previous' or 'linear',
+                'previous' or 'linear',
             fill_from: As an alternative to ``update_point``, an
-              :class:`numpy.timedelta64` value giving the time of day from which to fill
-              values forward.
+                :class:`numpy.timedelta64` value giving the time of day from which to
+                fill values forward.
         """
 
         if values.shape[0] != self.n_days:
-            raise ValueError("Values is not of length n_days on its first axis.")
+            raise ValueError("Values is not of length n_days on its first axis")
 
         if fill_from is not None:
             if not isinstance(fill_from, np.timedelta64):
-                raise ValueError("The fill_from argument must be a timedelta64 value.")
+                raise ValueError("The fill_from argument must be a timedelta64 value")
 
             # Convert to seconds and check it is in range
             _fill_from = fill_from.astype("timedelta64[s]")
             if not (_fill_from >= 0 and _fill_from < 24 * 60 * 60):
-                raise ValueError("The fill_from argument is not >= 0 and < 24 hours.")
+                raise ValueError("The fill_from argument is not >= 0 and < 24 hours")
 
             update_time = self.observation_dates + _fill_from
 
@@ -468,8 +477,23 @@ def fill_daily_to_subdaily(
         else:
             raise ValueError("Unknown update point")
 
-        # Note that interp1d cannot handle datetime64 inputs, so need to interpolate
-        # using datetimes cast to integer types
+        # Check previous value settings - only allow with previous interpolation and
+        # check the previous value shape matches
+        if previous_value is not None:
+            if kind == "linear":
+                raise NotImplementedError(
+                    "Using previous value with kind='linear' is not implemented"
+                )
+
+            # Use np.broadcast_shapes here to handle checking array shapes. This is
+            # mostly to catch the fact that () and (1,) are equivalent.
+            try:
+                np.broadcast_shapes(previous_value.shape, values.shape)
+            except ValueError:
+                raise ValueError(
+                    "The input to previous_value is not congruent with "
+                    "the shape of the observed data"
+                )
 
         # Use fill_value to handle extrapolation before or after update point:
         # - previous will fill the last value forward to the end of the time series,
@@ -479,8 +503,16 @@ def fill_daily_to_subdaily(
         #   value until _after_ the update point.
 
         if kind == "previous":
-            fill_value = (None, values[-1])
+            # The fill values here are used to extend the last daily value out to the
+            # end of the subdaily observations but also to fill any provided previous
+            # values for subdaily observations _before_ the first daily value. If
+            # the default previous value of None is supplied, this inserts np.nan as
+            # expected.
+            fill_value = (previous_value, values[-1])
         elif kind == "linear":
+            # Shift the values forward by a day, inserting a copy of the first day at
+            # the start. This then avoids plants seeing the future and provides values
+            # up until the last observation.
             values = np.insert(values, 0, values[0], axis=0)
             update_time = np.append(
                 update_time, update_time[-1] + np.timedelta64(1, "D")
@@ -489,13 +521,16 @@ def fill_daily_to_subdaily(
         else:
             raise ValueError("Unsupported interpolation option")
 
+        # Note that interp1d cannot handle datetime64 inputs, so need to interpolate
+        # using datetimes cast to integer types
         interp_fun = interp1d(
             update_time.astype("int"),
             values,
             axis=0,
             kind=kind,
             bounds_error=False,
             fill_value=fill_value,
+            assume_sorted=True,
         )
 
         # TODO - The kind "previous" might be replaceable with bottleneck.push

pyrealm/pmodel/subdaily.py

@@ -41,7 +41,10 @@
 
 
 def memory_effect(
-    values: NDArray, alpha: float = 0.067, allow_holdover: bool = False
+    values: NDArray,
+    previous_values: NDArray | None = None,
+    alpha: float = 0.067,
+    allow_holdover: bool = False,
 ) -> NDArray:
     r"""Apply a memory effect to a variable.
 
@@ -89,6 +92,8 @@ def memory_effect(
 
     Args:
         values: The values to apply the memory effect to.
+        previous_values: Last available realised value used if model is fitted in
+            chunks and value at t=0 is not optimal.
         alpha: The relative weight applied to the most recent observation.
         allow_holdover: Allow missing values to be filled by holding over earlier
             values.
@@ -105,7 +110,10 @@ def memory_effect(
     # Initialise the output storage and set the first values to be a slice along the
     # first axis of the input values
     memory_values = np.empty_like(values, dtype=np.float32)
-    memory_values[0] = values[0]
+    if previous_values is None:
+        memory_values[0] = values[0]
+    else:
+        memory_values[0] = previous_values * (1 - alpha) + values[0] * alpha
 
     # Handle the data if there are no missing data,
     if not nan_present:
@@ -177,6 +185,10 @@ class SubdailyPModel:
       more rapid acclimation: :math:`\alpha=1` results in immediate acclimation and
       :math:`\alpha=0` results in no acclimation at all, with values pinned to the
       initial estimates.
+    * By default, the initial realised value :math:`R_1` for each of the three slowly
+      acclimating variables is assumed to be the first optimal value :math:`O_1`, but
+      the `previous_realised` argument can be used to provide values of :math:`R_0` from
+      which to calculate :math:`R_{1} = R_{0}(1 - \alpha) + O_{1} \alpha`.
     * The realised values are then filled back onto the original subdaily timescale,
       with :math:`V_{cmax}` and :math:`J_{max}` then being calculated from the slowly
       responding :math:`V_{cmax25}` and :math:`J_{max25}` and the actual subdaily
@@ -221,10 +233,12 @@ class SubdailyPModel:
         allow_partial_data: Should estimates of daily optimal conditions be calculated
           with missing values in the acclimation window.
         reference_kphio: An optional alternative reference value for the quantum yield
-            efficiency of photosynthesis (:math:`\phi_0`, -) to be passed to the kphio
-            calculation method.
+          efficiency of photosynthesis (:math:`\phi_0`, -) to be passed to the kphio
+          calculation method.
         fill_kind: The approach used to fill daily realised values to the subdaily
           timescale, currently one of 'previous' or 'linear'.
+        previous_realised: A tuple of previous realised values of three NumPy arrays
+          (xi_real, vcmax25_real, jmax25_real).
     """
 
     def __init__(
@@ -241,6 +255,7 @@ def __init__(
         allow_holdover: bool = False,
         allow_partial_data: bool = False,
         fill_kind: str = "previous",
+        previous_realised: tuple[NDArray, NDArray, NDArray] | None = None,
     ) -> None:
         # Warn about the API
         warn(
@@ -379,25 +394,73 @@ def __init__(
             1 / calc_ftemp_arrh(tk_acclim, self.env.pmodel_const.subdaily_jmax25_ha)
         )
 
+        """Instantaneous optimal :math:`x_{i}`, :math:`V_{cmax}` and :math:`J_{max}`"""
+        # Check the shape of previous realised values are congruent with a slice across
+        # the time axis
+        if previous_realised is not None:
+            if fill_kind != "previous":
+                raise NotImplementedError(
+                    "Using previous_realised is only implemented for "
+                    "fill_kind = 'previous'"
+                )
+
+            # All variables should share the shape of a slice along the first axis of
+            # the environmental forcings
+            expected_shape = self.env.tc[0].shape
+            if not (
+                (previous_realised[0].shape == expected_shape)
+                and (previous_realised[1].shape == expected_shape)
+                and (previous_realised[2].shape == expected_shape)
+            ):
+                raise ValueError(
+                    "`previous_realised` entries have wrong shape in Subdaily PModel"
+                )
+            else:
+                previous_xi_real, previous_vcmax25_real, previous_jmax25_real = (
+                    previous_realised
+                )
+        else:
+            previous_xi_real, previous_vcmax25_real, previous_jmax25_real = [
+                None,
+                None,
+                None,
+            ]
+
         # 5) Calculate the realised daily values from the instantaneous optimal values
         self.xi_real: NDArray = memory_effect(
-            self.pmodel_acclim.optchi.xi, alpha=alpha, allow_holdover=allow_holdover
+            self.pmodel_acclim.optchi.xi,
+            previous_values=previous_xi_real,
+            alpha=alpha,
+            allow_holdover=allow_holdover,
         )
         r"""Realised daily slow responses in :math:`\xi`"""
         self.vcmax25_real: NDArray = memory_effect(
-            self.vcmax25_opt, alpha=alpha, allow_holdover=allow_holdover
+            self.vcmax25_opt,
+            previous_values=previous_vcmax25_real,
+            alpha=alpha,
+            allow_holdover=allow_holdover,
         )
         r"""Realised daily slow responses in :math:`V_{cmax25}`"""
         self.jmax25_real: NDArray = memory_effect(
-            self.jmax25_opt, alpha=alpha, allow_holdover=allow_holdover
+            self.jmax25_opt,
+            previous_values=previous_jmax25_real,
+            alpha=alpha,
+            allow_holdover=allow_holdover,
         )
+
         r"""Realised daily slow responses in :math:`J_{max25}`"""
 
         # 6) Fill the realised xi, jmax25 and vcmax25 from daily values back to the
         # subdaily timescale.
-        self.subdaily_vcmax25 = fs_scaler.fill_daily_to_subdaily(self.vcmax25_real)
-        self.subdaily_jmax25 = fs_scaler.fill_daily_to_subdaily(self.jmax25_real)
-        self.subdaily_xi = fs_scaler.fill_daily_to_subdaily(self.xi_real)
+        self.subdaily_xi = fs_scaler.fill_daily_to_subdaily(
+            self.xi_real, previous_value=previous_xi_real
+        )
+        self.subdaily_vcmax25 = fs_scaler.fill_daily_to_subdaily(
+            self.vcmax25_real, previous_value=previous_vcmax25_real
+        )
+        self.subdaily_jmax25 = fs_scaler.fill_daily_to_subdaily(
+            self.jmax25_real, previous_value=previous_jmax25_real
+        )
 
         # 7) Adjust subdaily jmax25 and vcmax25 back to jmax and vcmax given the
         #    actual subdaily temperatures.

tests/unit/pmodel/test_memory_effect.py

@@ -54,6 +54,43 @@ def test_memory_effect(inputs, alpha):
     assert np.allclose(result, expected)
 
 
+@pytest.mark.parametrize(
+    argnames="inputs_whole",
+    argvalues=[
+        pytest.param(np.arange(0, 10), id="1D"),
+        pytest.param(
+            np.column_stack([np.arange(0, 10)] * 4) + np.arange(4),
+            id="2D",
+        ),
+        pytest.param(
+            np.dstack([np.column_stack([np.arange(0, 10)] * 4)] * 4)
+            + np.arange(16).reshape(4, 4),
+            id="3D",
+        ),
+    ],
+)
+@pytest.mark.parametrize(argnames="alpha", argvalues=(0.0, 0.5, 1.0))
+def test_memory_effect_chunked(inputs_whole, alpha):
+    """Test that the memory effect works when chunking the time series up.
+
+    This compares the output of `test_memory_effect` with the output of
+    `memory_effect` which gets the two time chunks fed sequentially.,
+    """
+    from pyrealm.pmodel import memory_effect
+
+    result_whole = memory_effect(inputs_whole, alpha=alpha)
+
+    [inputs_chunk1, inputs_chunk2] = np.split(inputs_whole, [5], axis=0)
+
+    result_chunk1 = memory_effect(inputs_chunk1, alpha=alpha)
+
+    result_chunk2 = memory_effect(
+        inputs_chunk2, previous_values=result_chunk1[-1], alpha=alpha
+    )
+
+    assert np.allclose(result_whole[-1], result_chunk2[-1])
+
+
 @pytest.mark.parametrize(
     argnames="inputs,allow_holdover,context_manager,expected",
     argvalues=[