1818)
1919
2020
21- def memory_effect (values : NDArray , alpha : float = 0.067 ) -> NDArray :
21+ def memory_effect (
22+ values : NDArray , alpha : float = 0.067 , handle_nan : bool = False
23+ ) -> NDArray :
2224 r"""Apply a memory effect to a variable.
2325
2426 Three key photosynthetic parameters (:math:`\xi`, :math:`V_{cmax25}` and
@@ -27,37 +29,91 @@ def memory_effect(values: NDArray, alpha: float = 0.067) -> NDArray:
2729 average to apply a lagged response to one of these parameters.
2830
2931 The estimation uses the paramater `alpha` (:math:`\alpha`) to control the speed of
30- convergence of the estimated values (:math:`E `) to the calculated optimal values
32+ convergence of the realised values (:math:`R `) to the calculated optimal values
3133 (:math:`O`):
3234
3335 .. math::
3436
35- E_ {t} = E_ {t-1}(1 - \alpha) + O_{t} \alpha
37+ R_ {t} = R_ {t-1}(1 - \alpha) + O_{t} \alpha
3638
37- For :math:`t_{0}`, the first value in the optimal values is used so :math:`E_ {0} =
39+ For :math:`t_{0}`, the first value in the optimal values is used so :math:`R_ {0} =
3840 O_{0}`.
3941
4042 The ``values`` array can have multiple dimensions but the first dimension is always
4143 assumed to represent time and the memory effect is calculated only along the first
4244 dimension.
4345
46+ By default, the ``values`` array must not contain missing values (`numpy.nan`).
47+ However, :math:`V_{cmax}` and :math:`J_{max}` are not estimable in some conditions
48+ (namely when :math:`m \le c^{\ast}`, see
49+ :class:`~pyrealm.pmodel.pmodel.CalcOptimalChi`) and so missing values in P Model
50+ predictions can arise even when the forcing data is complete, breaking the recursion
51+ shown above. When ``handle_nan=True``, this function fills missing data as follow:
52+
53+ +-------------------+--------+-------------------------------------------------+
54+ | | | Current optimal (:math:`O_{t}`) |
55+ +-------------------+--------+-----------------+-------------------------------+
56+ | | | NA | not NA |
57+ +-------------------+--------+-----------------+-------------------------------+
58+ | Previous | NA | NA | O_{t} |
59+ | realised +--------+-----------------+-------------------------------+
60+ | (:math:`R_{t-1}`) | not NA | :math:`R_{t-1}` | :math:`R_{t-1}(1-a) + O_{t}a` |
61+ +-------------------+--------+-----------------+-------------------------------+
62+
63+ Initial missing values are kept, and the first observed optimal value is accepted as
64+ the first realised value (as with the start of the recursion above). After this, if
65+ the current optimal value is missing, then the previous estimate of the realised
66+ value is held over until it can next be updated from observed data.
67+
4468 Args:
4569 values: The values to apply the memory effect to.
46- alpha: The relative weight applied to the most recent observation
70+ alpha: The relative weight applied to the most recent observation.
71+ handle_nan: Allow missing values to be handled.
4772
4873 Returns:
4974 An array of the same shape as ``values`` with the memory effect applied.
5075 """
5176
77+ # Check for nan and nan handling
78+ nan_present = np .any (np .isnan (values ))
79+ if nan_present and not handle_nan :
80+ raise ValueError ("Missing values in data passed to memory_effect" )
81+
5282 # Initialise the output storage and set the first values to be a slice along the
5383 # first axis of the input values
5484 memory_values = np .empty_like (values , dtype = np .float32 )
5585 memory_values [0 ] = values [0 ]
5686
57- # Loop over the first axis, in each case taking slices through the first axis of the
58- # inputs. This handles arrays of any dimension.
87+ # Handle the data if there are no missing data,
88+ if not nan_present :
89+ # Loop over the first axis, in each case taking slices through the first axis of
90+ # the inputs. This handles arrays of any dimension.
91+ for idx in range (1 , len (memory_values )):
92+ memory_values [idx ] = (
93+ memory_values [idx - 1 ] * (1 - alpha ) + values [idx ] * alpha
94+ )
95+
96+ return memory_values
97+
98+ # Otherwise, do the same thing but handling missing data at each step.
5999 for idx in range (1 , len (memory_values )):
60- memory_values [idx ] = memory_values [idx - 1 ] * (1 - alpha ) + values [idx ] * alpha
100+ # Need to check for nan conditions:
101+ # - the previous value might be nan from an initial nan or sequence of nans, in
102+ # which case the current value is accepted without weighting - it could be nan
103+ # itself to extend a chain of initial nan values.
104+ # - the current value might be nan, in which case the previous value gets
105+ # held over as the current value.
106+ prev_nan = np .isnan (memory_values [idx - 1 ])
107+ curr_nan = np .isnan (values [idx ])
108+ memory_values [idx ] = np .where (
109+ prev_nan ,
110+ values [idx ],
111+ np .where (
112+ curr_nan ,
113+ memory_values [idx - 1 ],
114+ memory_values [idx - 1 ] * (1 - alpha ) + values [idx ] * alpha ,
115+ ),
116+ )
61117
62118 return memory_values
63119
@@ -92,7 +148,8 @@ class FastSlowPModel:
92148 * The :meth:`~pyrealm.pmodel.subdaily.memory_effect` function is then used to
93149 calculate realised slowly responding values for :math:`\xi`, :math:`V_{cmax25}`
94150 and :math:`J_{max25}`, given a weight :math:`\alpha \in [0,1]` that sets the speed
95- of acclimation.
151+ of acclimation. The ``handle_nan`` argument is passed to this function to set
152+ whether missing values in the optimal predictions are permitted and handled.
96153 * The realised values are then filled back onto the original subdaily timescale,
97154 with :math:`V_{cmax}` and :math:`J_{max}` then being calculated from the slowly
98155 responding :math:`V_{cmax25}` and :math:`J_{max25}` and the actual subdaily
@@ -107,6 +164,8 @@ class FastSlowPModel:
107164 fapar: The :math:`f_{APAR}` for each observation.
108165 ppfd: The PPDF for each observation.
109166 alpha: The :math:`\alpha` weight.
167+ handle_nan: Should the :func:`~pyrealm.pmodel.subdaily.memory_effect` function
168+ be allowe to handle missing values.
110169 kphio: The quantum yield efficiency of photosynthesis (:math:`\phi_0`, -).
111170 fill_kind: The approach used to fill daily realised values to the subdaily
112171 timescale, currently one of 'previous' or 'linear'.
@@ -120,6 +179,7 @@ def __init__(
120179 fapar : NDArray ,
121180 alpha : float = 1 / 15 ,
122181 kphio : float = 1 / 8 ,
182+ handle_nan : bool = False ,
123183 fill_kind : str = "previous" ,
124184 ) -> None :
125185 # Warn about the API
@@ -192,11 +252,17 @@ def __init__(
192252 )
193253
194254 # Calculate the realised values from the instantaneous optimal values
195- self .xi_real : NDArray = memory_effect (self .pmodel_acclim .optchi .xi , alpha = alpha )
255+ self .xi_real : NDArray = memory_effect (
256+ self .pmodel_acclim .optchi .xi , alpha = alpha , handle_nan = handle_nan
257+ )
196258 r"""Realised daily slow responses in :math:`\xi`"""
197- self .vcmax25_real : NDArray = memory_effect (self .vcmax25_opt , alpha = alpha )
259+ self .vcmax25_real : NDArray = memory_effect (
260+ self .vcmax25_opt , alpha = alpha , handle_nan = handle_nan
261+ )
198262 r"""Realised daily slow responses in :math:`V_{cmax25}`"""
199- self .jmax25_real : NDArray = memory_effect (self .jmax25_opt , alpha = alpha )
263+ self .jmax25_real : NDArray = memory_effect (
264+ self .jmax25_opt , alpha = alpha , handle_nan = handle_nan
265+ )
200266 r"""Realised daily slow responses in :math:`J_{max25}`"""
201267
202268 # Fill the daily realised values onto the subdaily scale
@@ -292,6 +358,8 @@ class FastSlowPModel_JAMES:
292358 fapar: The :math:`f_{APAR}` for each observation.
293359 ppfd: The PPDF for each observation.
294360 alpha: The :math:`\alpha` weight.
361+ handle_nan: Should the :func:`~pyrealm.pmodel.subdaily.memory_effect` function
362+ be allowe to handle missing values.
295363 kphio: The quantum yield efficiency of photosynthesis (:math:`\phi_0`, -).
296364 vpd_scaler: An alternate
297365 :class:`~pyrealm.pmodel.fast_slow_scaler.FastSlowScaler` instance used to
@@ -310,6 +378,7 @@ def __init__(
310378 ppfd : NDArray ,
311379 fapar : NDArray ,
312380 alpha : float = 1 / 15 ,
381+ handle_nan : bool = False ,
313382 kphio : float = 1 / 8 ,
314383 vpd_scaler : Optional [FastSlowScaler ] = None ,
315384 fill_from : Optional [np .timedelta64 ] = None ,
@@ -388,9 +457,13 @@ def __init__(
388457 )
389458
390459 # Calculate the realised values from the instantaneous optimal values
391- self .vcmax25_real : NDArray = memory_effect (self .vcmax25_opt , alpha = alpha )
460+ self .vcmax25_real : NDArray = memory_effect (
461+ self .vcmax25_opt , alpha = alpha , handle_nan = handle_nan
462+ )
392463 r"""Realised daily slow responses in :math:`V_{cmax25}`"""
393- self .jmax25_real : NDArray = memory_effect (self .jmax25_opt , alpha = alpha )
464+ self .jmax25_real : NDArray = memory_effect (
465+ self .jmax25_opt , alpha = alpha , handle_nan = handle_nan
466+ )
394467 r"""Realised daily slow responses in :math:`J_{max25}`"""
395468
396469 # Calculate the daily xi value, which does not have a slow reponse in this
0 commit comments