lactationcurve.characteristics.lactation_curve_characteristics

  1"""
  2Lactation Curve Characteristics (LCC) utilities.
  3
  4This module derives and evaluates **lactation curve characteristics** (LCCs) for
  5common lactation models using symbolic calculus (via SymPy) and fast numeric
  6evaluation (via NumPy). Characteristics include:
  7
  8- **time_to_peak**: day in milk (DIM) at which the model reaches its maximum
  9  (derived where the first derivative is zero).
 10- **peak_yield**: the yield at the time of peak.
 11- **cumulative_milk_yield**: total yield over a lactation horizon (default 305 days),
 12  computed as the definite integral of the model over time.
 13- **persistency**: the average slope after peak until the end of lactation (derived),
 14  or literature-based formulas (for Wood and MilkBot).
 15
 16Features
 17--------
 18- **Symbolic derivation** (derivative/solve and integration) **cached** per
 19  (model, characteristic) to avoid recomputation.
 20- **Fallback numeric methods** for robustness when the symbolic expression is
 21  not suitable for lambdification or yields invalid values.
 22- Works with frequentist- or Bayesian-fitted parameters via the
 23  `lactationcurve.fitting` API.
 24
 25Notes
 26-----
 27- Units: DIM in days, milk in kg or lbs.
 28- Symbolic expressions can be complex; a light safety check is applied before
 29  lambdification (`is_valid_sympy_expr`) and large/invalid expressions are rejected.
 30- The default laction curve model is the Wood model, but other models are supported as well.
 31
 32
 33Author: Meike van Leerdam
 34Last update: 21-May-2026
 35"""
 36
 37import numpy as np
 38from sympy import (
 39    Integral,
 40    diff,
 41    exp,
 42    integrate,
 43    lambdify,
 44    ln,
 45    log,
 46    nan,
 47    oo,
 48    simplify,
 49    solve,
 50    symbols,
 51    zoo,
 52)
 53
 54from lactationcurve.fitting import (
 55    bayesian_fit_milkbot_single_lactation,
 56    fit_lactation_curve,
 57    get_lc_parameters,
 58)
 59from lactationcurve.preprocessing import validate_and_prepare_inputs
 60
 61
 62# safe guard for extreme expressions or weird results
 63def is_valid_sympy_expr(expr) -> bool:
 64    """Check whether a SymPy expression looks safe to lambdify.
 65
 66    The validation rejects expressions that contain infinities, NaNs, or an
 67    excessive number of operations.
 68
 69    Args:
 70        expr: A SymPy expression to validate.
 71
 72    Returns:
 73        True if the expression appears safe to lambdify; False otherwise.
 74
 75    Notes:
 76        - Expressions containing `oo`, `-oo`, `zoo`, or `nan` are rejected.
 77        - Expressions with `count_ops() > 5000` are rejected as too large/complex.
 78    """
 79    try:
 80        if expr.has(oo, -oo, zoo, nan):
 81            return False
 82        # reject absurdly large expressions
 83        if expr.count_ops() > 5000:
 84            return False
 85        return True
 86    except Exception:
 87        return False
 88
 89
 90# Store derived function in cache so they do not need to be calculated all the time again.
 91_LCC_CACHE = {}
 92
 93
 94def lactation_curve_characteristic_function(
 95    model="wood", characteristic=None, lactation_length=305
 96) -> tuple:
 97    """Build (or fetch from cache) a symbolic expression and fast numeric function for an LCC.
 98
 99    This function derives the requested **lactation curve characteristic** for a given
100    model using SymPy (derivative / root finding / integration). It returns the symbolic
101    expression, the tuple of parameter symbols (argument order), and a lambdified
102    numeric function that can be evaluated with numerical parameters.
103
104    The symbolic derivation and integration are done only once per (model, characteristic)
105    and then **cached** for reuse.
106
107    Args:
108        model (str): Model name. Options:
109            'milkbot', 'wood', 'wilmink', 'ali_schaeffer', 'fischer',
110            'brody', 'sikka', 'nelder', 'dhanoa', 'emmans', 'hayashi',
111            'rook', 'dijkstra', 'prasad'.
112        characteristic (str | None): Desired characteristic. Options:
113            'time_to_peak', 'peak_yield', 'cumulative_milk_yield', 'persistency'.
114            If `None` or unrecognized, a dict of all available characteristics is returned
115            (with `persistency` possibly `None` if derivation is not feasible).
116        lactation_length (int): Length of lactation in days used in persistency
117            and cumulative milk yield computation (default 305).
118
119    Returns:
120        tuple:
121            expr: SymPy expression (or dict of expressions if `characteristic` is None).
122            params: Tuple of SymPy symbols for model parameters (argument order).
123            func: Lambdified numeric function `f(*params)` (or dict of functions).
124
125    Raises:
126        Exception: If the model is unknown, or if no positive real solution for
127            peak timing/yield exists where required. In this case None is returned
128             for the characteristic,
129             and the function can be evaluated with numeric methods as a fallback.
130    """
131    # check the cache
132    storage = (model, characteristic)
133    if storage in _LCC_CACHE:
134        return (
135            _LCC_CACHE[storage]["expr"],
136            _LCC_CACHE[storage]["params"],
137            _LCC_CACHE[storage]["func"],
138        )
139
140    # make sure model is all lowercase
141    model = model.lower()
142
143    # define functions
144    if model == "brody":
145        # === BRODY 1 ===
146        a, b, k1, k2, t = symbols("a b k1 k2 t", real=True, positive=True)
147        function = a * exp(-k1 * t) - b * exp(-k2 * t)
148
149    elif model == "sikka":
150        # === SIKKA ===
151        a, b, c, t = symbols("a b c t", real=True, positive=True)
152        function = a * exp(b * t - c * t**2)
153
154    elif model == "fischer":
155        # === FISCHER ===
156        a, b, c, t = symbols("a b c t", real=True, positive=True)
157        function = a - b * t - a * exp(-c * t)
158
159    elif model == "nelder":
160        # === NELDER ===
161        a, b, c, t = symbols("a b c t", real=True, positive=True)
162        function = t / (a + b * t + c * t**2)
163
164    elif model == "wood":
165        # === WOOD ===
166        a, b, c, t = symbols("a b c t", real=True, positive=True)
167        function = a * t**b * exp(-c * t)
168
169    elif model == "dhanoa":
170        # === DHANOA ===
171        a, b, c, t = symbols("a b c t", real=True, positive=True)
172        function = a * t ** (b * c) * exp(-c * t)
173
174    elif model == "emmans":
175        # === EMMANS ===
176        a, b, c, d, t = symbols("a b c d t", real=True, positive=True)
177        function = a * exp(-exp(d - b * t)) * exp(-c * t)
178
179    elif model == "ali_schaeffer":
180        # ====ALI=====
181        a, b, c, d, k, t = symbols("a b c d k t", real=True)
182
183        t1 = t / 305
184        w1 = log(305 / t)
185
186        function = a + b * t1 + c * t1**2 + d * w1 + k * w1**2
187
188    elif model == "wilmink":
189        # === WILMINK ===
190        a, c, k, t = symbols("a c k t", real=True, positive=True)
191        b = symbols("b", real=True)  # b can be negative in Wilmink
192        function = a + b * exp(-k * t) - c * t
193
194    elif model == "hayashi":
195        # === HAYASHI ===
196        a, b, c, t = symbols("a b c t", real=True, positive=True)
197        function = b * (exp(-t / c) - exp(-t / (a * c)))
198
199    elif model == "rook":
200        # === ROOK ===
201        a, b, c, d, t = symbols("a b c d t", real=True, positive=True)
202        function = a * (1 / (1 + b / (c + t))) * exp(-d * t)
203
204    elif model == "dijkstra":
205        # === DIJKSTRA ===
206        a, b, c, d, t = symbols("a b c d t", real=True, positive=True)
207        function = a * exp((b * (1 - exp(-c * t)) / c) - d * t)
208
209    elif model == "prasad":
210        # === PRASAD ===
211        a, b, c, d, t = symbols("a b c d t", real=True, positive=True)
212        function = a + b * t + c * t**2 + d / t
213
214    elif model == "milkbot":
215        # === MILKBOT ===
216        a, b, c, d, t = symbols("a b c d t", real=True, positive=True)
217        function = a * (1 - exp((c - t) / b) / 2) * exp(-d * t)
218
219    else:
220        raise Exception("Unknown model")
221
222    # solve derivative for when it is zero to find the function for time of peak
223    fdiff = diff(function, t)
224
225    tpeak = None  # default
226
227    valid_roots = []
228
229    try:
230        candidate_roots = solve(fdiff, t)
231
232        for r in candidate_roots:
233            r_simplified = simplify(r)
234
235            if r_simplified.is_real is not False:
236                valid_roots.append(r_simplified)
237
238    except Exception:
239        valid_roots = []
240
241    tpeak = valid_roots[0] if valid_roots else None
242
243    # override to prevent sympy error
244    if model == "dijkstra":
245        tpeak = simplify(log(b / d) / c)
246
247    # define the end of lactation
248    T = lactation_length  # days in milk
249
250    # Persistency = average slope after peak, does not work for all models so therefore try except
251    persistency = None
252    try:
253        if tpeak:
254            tmp = (function.subs(t, T) - function.subs(t, tpeak)) / (T - tpeak)
255            tmp = tmp.cancel()  # light simplification
256            if is_valid_sympy_expr(tmp):
257                persistency = tmp
258    except Exception:
259        persistency = None
260
261    if characteristic != "cumulative_milk_yield":
262        if tpeak:  # Check if the list is not empty
263            peak_expr = simplify(function.subs(t, tpeak))
264        else:
265            raise Exception("No real solution for time to peak and peak yield found")
266
267    # find function for cumulative milk yield of the lactation.
268    # A lactation length of 305 is the default, but this value can
269    # also be provided as a parameter in the characteristic function.
270    cum_my_expr = integrate(function, (t, 0, lactation_length), meijerg=False)
271
272    if cum_my_expr == 0:
273        cum_my_expr = Integral(function, (t, 0, lactation_length))
274
275    # Sorted parameter list (exclude t)
276    params = tuple(
277        sorted([s for s in function.free_symbols if s.name != "t"], key=lambda x: x.name)
278    )
279
280    # ----------------------------------------------------
281    # Select requested characteristic
282    # ----------------------------------------------------
283    if characteristic == "time_to_peak":
284        expr = tpeak
285    elif characteristic == "peak_yield":
286        expr = peak_expr
287    elif characteristic == "cumulative_milk_yield":
288        expr = cum_my_expr
289    elif characteristic == "persistency":
290        if persistency is None:
291            raise Exception("Persistency could not be computed symbolically")
292        expr = persistency
293    else:
294        # Return all four if None or 'all'
295        expr = {
296            "time_to_peak": tpeak,
297            "peak_yield": peak_expr,
298            "persistency": persistency,  # possibly None
299            "cumulative_milk_yield": cum_my_expr,
300        }
301
302    # ----------------------------------------------------
303    # Build fast numeric function with lambdify
304    # ----------------------------------------------------
305    if isinstance(expr, dict):
306        func = {
307            name: lambdify(params, ex, modules=["numpy", "scipy"])
308            for name, ex in expr.items()
309            if ex is not None
310        }
311    else:
312        func = lambdify(params, expr, modules=["numpy", "scipy"])
313
314    # ----------------------------------------------------
315    # Store in cache
316    # ----------------------------------------------------
317    _LCC_CACHE[storage] = {"expr": expr, "params": params, "func": func}
318
319    return expr, params, func
320
321
322def calculate_characteristic(
323    dim,
324    milkrecordings,
325    model="wood",
326    characteristic="cumulative_milk_yield",
327    fitting="frequentist",
328    key=None,
329    parity=3,
330    breed="H",
331    continent="USA",
332    custom_priors=None,
333    milk_unit="kg",
334    persistency_method="derived",
335    lactation_length=305,
336) -> float | None:
337    """Evaluate a lactation curve characteristic from observed test-day data.
338
339    This function fits the requested model (frequentist or Bayesian via MilkBot),
340    retrieves model parameters, and evaluates the requested characteristic using the
341    symbolic expression (if available), falling back to numeric methods when needed.
342
343    Args:
344        dim (Int): Days in milk (DIM).
345        milkrecordings (Float): Milk recordings (kg or lbs) for each DIM.
346        model (str): Model name. Supported for this function:
347            'milkbot', 'wood', 'wilmink', 'ali_schaeffer', 'fischer'.
348        characteristic (str): One of 'time_to_peak', 'peak_yield',
349            'cumulative_milk_yield', 'persistency'.
350        fitting (str): 'frequentist' (default) or 'bayesian'.
351        key (str | None): API key for MilkBot Bayesian fitting.
352        parity (Int): Parity of the cow; values above 3 are considered as 3 (Bayesian).
353        breed (str): 'H' (Holstein) or 'J' (Jersey) (Bayesian).
354        custom_priors: provide your own priors for Bayesian fitting.
355            provide as dictionary using build_prior() function
356            to create the priors in the right format.
357            Alternative use priors from the literature provided by the string command 'CHEN'
358        milk_unit: Unit of milk recordings ('kg' or 'lb') for Bayesian fitting.
359        continent (str): 'USA' or 'EU' (Bayesian).
360        persistency_method (str): 'derived' (average slope after peak; default) or 'literature'
361            (only for Wood and MilkBot).
362        lactation_length (Int | str): Horizon for persistency calculation:
363            305 (default), 'max' (use max DIM in data), or an integer.
364
365    Returns:
366        float | None: The requested characteristic value.
367        Or None if value is negative (except persistency)
368        or cannot be computed symbolically.
369
370    Raises:
371        Exception: If inputs are invalid, the model/characteristic is unsupported,
372            an API key is missing for Bayesian fitting, or the characteristic cannot be computed.
373    """
374    # check and prepare input
375    inputs = validate_and_prepare_inputs(
376        dim,
377        milkrecordings,
378        model=model,
379        fitting=fitting,
380        breed=breed,
381        parity=parity,
382        continent=continent,
383        custom_priors=custom_priors,
384        milk_unit=milk_unit,
385        persistency_method=persistency_method,
386        lactation_length=lactation_length,
387    )
388
389    dim = inputs.dim
390    milkrecordings = inputs.milkrecordings
391    model = inputs.model
392    fitting = inputs.fitting
393    breed = inputs.breed
394    parity = inputs.parity
395    continent = inputs.continent
396    custom_priors = inputs.custom_priors
397    milk_unit = inputs.milk_unit
398    persistency_method = inputs.persistency_method
399    lactation_length = inputs.lactation_length
400
401    # After validation, these fields are guaranteed non-None for the paths below
402    assert model is not None
403    assert fitting is not None
404
405    if model not in ["milkbot", "wood", "wilmink", "ali_schaeffer", "fischer"]:
406        msg = (
407            "this function currently only works for"
408            " the milkbot, wood, wilmink,"
409            " ali_schaeffer and fischer models"
410        )
411        raise Exception(msg)
412
413    characteristic_options: list[str] = [
414        "time_to_peak",
415        "peak_yield",
416        "cumulative_milk_yield",
417        "persistency",
418    ]
419
420    if characteristic in characteristic_options:
421        if fitting == "frequentist":
422            # Get fitted parameters from your fitting function
423            fitted_params = get_lc_parameters(dim, milkrecordings, model)
424            assert fitted_params is not None
425
426            if characteristic != "persistency":
427                # Try symbolic formula first
428                assert isinstance(lactation_length, int)
429                try:
430                    expr, params, fn = lactation_curve_characteristic_function(
431                        model, characteristic, lactation_length
432                    )
433                    with np.errstate(
434                        divide="ignore", invalid="ignore"
435                    ):  # get rid of warnings for invalid operations
436                        value = fn(*fitted_params)
437                except Exception:
438                    value = None
439
440                # If symbolic formula fails or is invalid (use numeric approach)
441                if (
442                    value is None
443                    or not np.isfinite(value)
444                    or (characteristic == "time_to_peak" and value <= 0)
445                ):
446                    assert parity is not None
447                    assert breed is not None
448                    assert continent is not None
449                    if characteristic == "time_to_peak":
450                        value = numeric_time_to_peak(
451                            dim,
452                            milkrecordings,
453                            model,
454                            fitting=fitting,
455                            key=key,
456                            parity=parity,
457                            breed=breed,
458                            continent=continent,
459                        )
460                    elif characteristic == "peak_yield":
461                        value = numeric_peak_yield(
462                            dim,
463                            milkrecordings,
464                            model,
465                            fitting=fitting,
466                            key=key,
467                            parity=parity,
468                            breed=breed,
469                            continent=continent,
470                        )
471                    elif characteristic == "cumulative_milk_yield":
472                        value = numeric_cumulative_yield(
473                            dim,
474                            milkrecordings,
475                            model,
476                            fitting=fitting,
477                            lactation_length=lactation_length,
478                            key=key,
479                            parity=parity,
480                            breed=breed,
481                            continent=continent,
482                        )
483
484            else:
485                if persistency_method == "derived":
486                    # find lactation length from data
487                    if lactation_length == "max":
488                        lactation_length = max(dim)
489                    elif isinstance(lactation_length, int):
490                        lactation_length = lactation_length
491                    else:
492                        lactation_length = 305
493                    value = persistency_fitted_curve(
494                        dim,
495                        milkrecordings,
496                        model,
497                        fitting="frequentist",
498                        lactation_length=lactation_length,
499                    )
500
501                else:
502                    if model == "wood":
503                        value = persistency_wood(fitted_params[1], fitted_params[2])
504
505                    elif model == "milkbot":
506                        assert len(fitted_params) >= 4
507                        value = persistency_milkbot(fitted_params[3])
508
509                    else:
510                        msg = (
511                            "Currently only the Wood model"
512                            " and MilkBot model have a"
513                            " separate model function from"
514                            " the literature integrated for"
515                            " persistency. if"
516                            ' persistency="derived" is'
517                            " selected, persistency can be"
518                            " calculated for every model as"
519                            " the average slope of the"
520                            " lactation after the peak."
521                        )
522                        raise Exception(msg)
523            try:
524                if value is None:
525                    return None
526                if characteristic != "persistency" and value < 0:
527                    return None
528                return float(value)
529            except ValueError:
530                raise Exception("Could not compute characteristic")
531
532        else:
533            if model == "milkbot":
534                if key is None:
535                    raise Exception("Key needed to use Bayesian fitting engine MilkBot")
536                else:
537                    assert parity is not None
538                    assert breed is not None
539                    assert continent is not None
540                    assert isinstance(milk_unit, str)
541                    fitted_params_bayes = bayesian_fit_milkbot_single_lactation(
542                        dim,
543                        milkrecordings,
544                        key=key,
545                        parity=parity,
546                        breed=breed,
547                        custom_priors=custom_priors,
548                        continent=continent,
549                        milk_unit=milk_unit,
550                    )
551                    fitted_params_bayes = (
552                        fitted_params_bayes["scale"],
553                        fitted_params_bayes["ramp"],
554                        fitted_params_bayes["offset"],
555                        fitted_params_bayes["decay"],
556                    )
557
558                    if characteristic != "persistency":
559                        # Get the symbolic expression and model parameters
560                        try:
561                            expr, param_symbols, fn = lactation_curve_characteristic_function(
562                                model, characteristic
563                            )
564                            with np.errstate(divide="ignore", invalid="ignore"):
565                                value = fn(*fitted_params_bayes)
566                        except Exception:
567                            value = None
568
569                        if (
570                            value is None
571                            or not np.isfinite(value)
572                            or (characteristic == "time_to_peak" and value <= 0)
573                        ):
574                            if characteristic == "time_to_peak":
575                                value = numeric_time_to_peak(
576                                    dim,
577                                    milkrecordings,
578                                    model,
579                                    fitting=fitting,
580                                    key=key,
581                                    parity=parity,
582                                    breed=breed,
583                                    custom_priors=custom_priors,
584                                    milk_unit=milk_unit,
585                                    continent=continent,
586                                )
587                            elif characteristic == "peak_yield":
588                                value = numeric_peak_yield(
589                                    dim,
590                                    milkrecordings,
591                                    model,
592                                    fitting=fitting,
593                                    key=key,
594                                    parity=parity,
595                                    breed=breed,
596                                    custom_priors=custom_priors,
597                                    milk_unit=milk_unit,
598                                    continent=continent,
599                                )
600                            elif characteristic == "cumulative_milk_yield":
601                                if lactation_length == "max":
602                                    numeric_lactation_length = max(dim)
603                                elif isinstance(lactation_length, int):
604                                    numeric_lactation_length = lactation_length
605                                else:
606                                    numeric_lactation_length = 305
607                                value = numeric_cumulative_yield(
608                                    dim,
609                                    milkrecordings,
610                                    model,
611                                    fitting=fitting,
612                                    lactation_length=numeric_lactation_length,
613                                    key=key,
614                                    parity=parity,
615                                    breed=breed,
616                                    custom_priors=custom_priors,
617                                    milk_unit=milk_unit,
618                                    continent=continent,
619                                )
620
621                    else:
622                        if persistency_method == "derived":
623                            # find lactation length from data
624                            if lactation_length == "max":
625                                lactation_length = max(dim)
626                            elif isinstance(lactation_length, int):
627                                lactation_length = lactation_length
628                            else:
629                                lactation_length = 305
630                            assert isinstance(milk_unit, str)
631                            value = persistency_fitted_curve(
632                                dim,
633                                milkrecordings,
634                                model,
635                                fitting="bayesian",
636                                key=key,
637                                parity=parity,
638                                breed=breed,
639                                custom_priors=custom_priors,
640                                milk_unit=milk_unit,
641                                continent=continent,
642                                lactation_length=lactation_length,
643                            )
644
645                        else:
646                            value = persistency_milkbot(fitted_params_bayes[3])
647
648                try:
649                    if value is None:
650                        return None
651                    if characteristic != "persistency" and value < 0:
652                        return None
653                    return float(value)
654                except ValueError:
655                    raise Exception("Could not compute characteristic")
656            else:
657                raise Exception("Bayesian fitting is currently only implemented for MilkBot models")
658
659    else:
660        raise Exception("Unknown characteristic")
661
662
663# also define numeric approaches as back up if symbolic functions fail or throw invalid results
664def numeric_time_to_peak(
665    dim,
666    milkrecordings,
667    model,
668    fitting="frequentist",
669    key=None,
670    parity=3,
671    breed="H",
672    custom_priors=None,
673    milk_unit="kg",
674    continent="USA",
675) -> int:
676    """Compute time to peak using a numeric approach.
677
678    Fits the curve (frequentist or Bayesian), evaluates the predicted yields,
679    and returns the DIM corresponding to the maximum predicted yield.
680
681    Args:
682        dim: DIM values.
683        milkrecordings: Milk recordings (kg).
684        model: Model name.
685        fitting: 'frequentist' or 'bayesian'.
686        key: API key for MilkBot (Bayesian).
687        parity: Parity for Bayesian fitting.
688        breed: Breed for Bayesian fitting ('H' or 'J').
689        custom_priors: provide your own priors for Bayesian fitting.
690            provide as dictionary using build_prior() function to set
691            the priors in the right format. Alternative use priors from
692            the literature provided by the string command 'CHEN'
693        milk_unit: Unit of milk recordings ('kg' or 'lb') for Bayesian fitting.
694        continent: Prior source for Bayesian ('USA', 'EU').
695
696
697    Returns:
698        int: DIM at which the curve attains its maximum (1-indexed).
699    """
700    # Fit the curve to get predicted milk yields
701    yields = fit_lactation_curve(
702        dim,
703        milkrecordings,
704        model,
705        fitting=fitting,
706        key=key,
707        parity=parity,
708        breed=breed,
709        custom_priors=custom_priors,
710        milk_unit=milk_unit,
711        continent=continent,
712    )
713    # Find the index of the peak yield
714    peak_idx = np.argmax(yields)
715    # Return the corresponding DIM
716    return int(peak_idx + 1)  # +1 because DIM starts at 1, not 0
717
718
719def numeric_cumulative_yield(
720    dim, milkrecordings, model, fitting="frequentist", lactation_length=305, **kwargs
721) -> float:
722    """Compute cumulative milk yield numerically over a given horizon.
723
724    Adds up the fitted milk yield for the first `lactation_length` days of the
725    predicted yield curve.
726
727    Args:
728        dim: DIM values.
729        milkrecordings: Milk recordings (kg).
730        model: Model name.
731        fitting: 'frequentist' or 'bayesian'.
732        lactation_length: Number of days to integrate (default 305).
733        **kwargs: Additional arguments forwarded to `fit_lactation_curve`.
734
735    Returns:
736        float: Cumulative milk yield over the specified horizon.
737    """
738    y = fit_lactation_curve(dim, milkrecordings, model, fitting=fitting, **kwargs)
739    return float(np.trapezoid(y[:lactation_length], dx=1))
740
741
742def numeric_peak_yield(
743    dim,
744    milkrecordings,
745    model,
746    fitting="frequentist",
747    key=None,
748    parity=3,
749    breed="H",
750    custom_priors=None,
751    milk_unit="kg",
752    continent="USA",
753) -> float:
754    """Compute peak yield numerically from the fitted curve.
755
756    Args:
757        dim: DIM values.
758        milkrecordings: Milk recordings (kg).
759        model: Model name.
760        fitting: 'frequentist' or 'bayesian'.
761        key: API key for MilkBot (Bayesian).
762        parity: Parity for Bayesian fitting.
763        breed: Breed for Bayesian fitting ('H' or 'J').
764        continent: Prior source for Bayesian ('USA', 'EU', 'CHEN').
765
766    Returns:
767        float: Maximum predicted milk yield.
768    """
769    # Fit the curve to get predicted milk yields
770    yields = fit_lactation_curve(
771        dim,
772        milkrecordings,
773        model,
774        fitting=fitting,
775        key=key,
776        parity=parity,
777        breed=breed,
778        custom_priors=custom_priors,
779        milk_unit=milk_unit,
780        continent=continent,
781    )
782    # Find the peak yield
783    peak_yield = np.max(yields)
784    return peak_yield
785
786
787def persistency_wood(b, c) -> float:
788    """Persistency from Wood et al. (1984): `Persistency = -(b + 1) * ln(c)`.
789
790    Args:
791        b (float): Parameter `b` of the Wood model.
792        c (float): Parameter `c` of the Wood model.
793
794    Returns:
795        float: Persistency value from the Wood formula.
796    """
797    return float(-(b + 1) * ln(c))
798
799
800def persistency_milkbot(d) -> float:
801    """Persistency from the MilkBot model (Ehrlich, 2013): `Persistency = 0.693 / d`.
802
803    Args:
804        d (float): Parameter `d` of the MilkBot model.
805
806    Returns:
807        float: Persistency value from the MilkBot formula.
808    """
809    return 0.693 / d
810
811
812def persistency_fitted_curve(
813    dim,
814    milkrecordings,
815    model,
816    fitting="frequentist",
817    key=None,
818    parity=3,
819    breed="H",
820    custom_priors=None,
821    milk_unit="kg",
822    continent="USA",
823    lactation_length=305,
824) -> float:
825    """Persistency as the average slope after peak until end of lactation (numeric).
826
827    This is the default approach because symbolic derivation is not feasible for
828    all models. It computes:
829        `(yield_at_end - yield_at_peak) / (lactation_length - time_to_peak)`
830
831    Args:
832        dim: DIM values.
833        milkrecordings: Milk recordings (kg).
834        model: Model name. Options include 'milkbot' (Bayesian or frequentist),
835            'wood', 'wilmink', 'ali_schaeffer', 'fischer'.
836        fitting: 'frequentist' or 'bayesian'.
837        key: API key (only for Bayesian fitting).
838        parity: Parity of the cow; values above 3 treated as 3 (Bayesian).
839        breed: 'H' or 'J' (Bayesian).
840        custom_priors: provide your own priors for Bayesian fitting.
841            provide as dictionary using build_prior() function to
842            create the priors in the right format. Alternative use
843            priors from the literature provided by the string
844            command 'CHEN'
845        milk_unit: Unit of milk recordings ('kg' or 'lb') for Bayesian fitting.
846        continent: 'USA', 'EU', or 'CHEN' (Bayesian).
847        lactation_length (int | str): 305 (default), 'max' (use max DIM), or integer.
848
849    Returns:
850        float: Average slope after the peak until end of lactation.
851    """
852    # calculate time to peak
853    t_peak = numeric_time_to_peak(
854        dim,
855        milkrecordings,
856        model,
857        fitting=fitting,
858        key=key,
859        parity=parity,
860        breed=breed,
861        custom_priors=custom_priors,
862        milk_unit=milk_unit,
863        continent=continent,
864    )
865
866    # determine lactation length
867    if lactation_length == "max":
868        lactation_length = max(dim)
869    elif isinstance(lactation_length, int):
870        lactation_length = lactation_length
871    else:
872        lactation_length = 305
873
874    # calculate milk yield at peak
875    yields = fit_lactation_curve(
876        dim,
877        milkrecordings,
878        model,
879        fitting=fitting,
880        key=key,
881        parity=parity,
882        breed=breed,
883        custom_priors=custom_priors,
884        milk_unit=milk_unit,
885        continent=continent,
886    )
887    peak_yield = yields[int(t_peak) - 1]  # -1 to prevent index error
888    # calculate milk yield at end of lactation
889    end_yield = yields[int(lactation_length) - 1]  # -1 to prevent index error
890
891    # calculate persistency
892    persistency = (end_yield - peak_yield) / (lactation_length - t_peak)
893    return persistency
894
895
896# add main to test
897def demo():
898    from sympy import pprint
899
900    models = [
901        "brody",
902        "sikka",
903        "fischer",
904        "nelder",
905        "wood",
906        "dhanoa",
907        "emmans",
908        "ali_schaeffer",
909        "wilmink",
910        "hayashi",
911        "rook",
912        "dijkstra",
913        "prasad",
914        "milkbot",
915    ]
916
917    characteristics = ["time_to_peak", "peak_yield", "cumulative_milk_yield"]
918
919    # simple test parameter values (biologically reasonable-ish)
920    test_values = {"a": 10, "b": 2, "c": 0.1, "d": 1, "k": 0.05, "k1": 0.1, "k2": 0.05}
921
922    for model in models:
923        print("\n" + "=" * 60)
924        print(f"MODEL: {model.upper()}")
925        print("=" * 60)
926
927        for char in characteristics:
928            print(f"\n--- {char} ---")
929
930            try:
931                expr, params, func = lactation_curve_characteristic_function(
932                    model=model, characteristic=char
933                )
934
935                # print symbolic expression
936                print("Symbolic:")
937                pprint(expr)
938
939                # build numeric input in correct order
940                try:
941                    param_values = [test_values[p.name] for p in params]
942
943                    numeric_value = func(*param_values)
944
945                    print("Numeric:")
946                    print(numeric_value)
947
948                except Exception as e:
949                    print("Numeric evaluation failed:", e)
950
951            except Exception as e:
952                print("Error:", e)
953
954
955if __name__ == "__main__":
956    demo()