lactationcurve.fitting

150def ali_schaeffer_model(t, a, b, c, d, k) -> np.floating | np.ndarray:
151    """Ali & Schaeffer lactation curve model.
152
153    Args:
154        t: Time since calving in days (DIM). Use `t >= 1` to avoid `log(0)`.
155        a: Fitted coefficient (numerical).
156        b: Fitted coefficient (numerical).
157        c: Fitted coefficient (numerical).
158        d: Fitted coefficient (numerical).
159        k: Fitted coefficient (numerical).
160
161    Returns:
162        Predicted milk yield at `t`.
163
164    Notes:
165        Formula: `t_scaled = t / 305` and `log_term = ln(305 / t)`
166        then `y(t) = a + b * t_scaled + c * t_scaled^2 + d * log_term + k * log_term^2`.
167    """
168    t_scaled = t / 305
169    log_term = np.log(305 / t)
170    return a + b * t_scaled + c * (t_scaled**2) + d * log_term + k * (log_term**2)

790def bayesian_fit_milkbot_single_lactation(
791    dim,
792    milkrecordings,
793    key: str,
794    parity=3,
795    breed="H",
796    custom_priors: MilkBotPriors | str | None = None,
797    continent="USA",
798    milk_unit="kg",
799) -> dict:
800    """
801    Fit a single lactation using the MilkBot API.
802
803    Args:
804        dim: List/array of DIM values.
805        milkrecordings: List/array of milk recordings (kg).
806        key: API key for MilkBot.
807        parity: Lactation number; values >= 3 are treated as one group in priors.
808        breed: "H" (Holstein) or "J" (Jersey).
809        custom_priors:
810            - "CHEN"  → Chen et al. published priors
811            - dict    → Custom priors in MilkBot format (overrides `continent`)
812        continent: priors used by MilkBot API for fitting; options:
813            - "USA"   → MilkBot USA priors
814            - "EU"    → MilkBot EU priors > estimates lower milk production
815
816
817    Returns:
818        Dictionary with fitted parameters and metadata:
819            {
820                "scale": float,
821                "ramp": float,
822                "decay": float,
823                "offset": float,
824                "nPoints": int
825            }
826
827    Raises:
828        requests.HTTPError: For unsuccessful HTTP response codes.
829        RuntimeError: If the response format is unexpected.
830
831    Notes:
832        - When `continent == "CHEN"`, Chen et al. priors are included in the request payload.
833        - EU calls use the GCP EU endpoint; others use `milkbot.com`.
834    """
835    # check and prepare input
836    inputs = validate_and_prepare_inputs(
837        dim,
838        milkrecordings,
839        breed=breed,
840        parity=parity,
841        custom_priors=custom_priors,
842        continent=continent,
843        milk_unit=milk_unit,
844    )
845
846    dim = inputs.dim
847    milkrecordings = inputs.milkrecordings
848    breed = inputs.breed
849    parity = inputs.parity
850    continent = inputs.continent
851    custom_priors = inputs.custom_priors
852    milk_unit = inputs.milk_unit
853
854    # -----------------------------
855    # Select server (USA vs EU)
856    # -----------------------------
857    if continent == "EU":
858        base_url = "https://europe-west1-numeric-analogy-337601.cloudfunctions.net/milkBot-fitter"
859    else:
860        base_url = "https://milkbot.com"
861
862    # -----------------------------
863    # Prepare headers
864    # -----------------------------
865    headers = {"Content-Type": "application/json", "X-API-KEY": key}
866
867    # -----------------------------
868    # Prepare milk points
869    # -----------------------------
870    points = sorted(
871        ({"dim": int(d), "milk": float(m)} for d, m in zip(dim, milkrecordings)),
872        key=lambda p: p["dim"],
873    )
874
875    # -----------------------------
876    # Lactation metadata
877    # -----------------------------
878    payload = {
879        "lactation": {
880            "lacKey": "single_lactation_fit",
881            "breed": breed,
882            "parity": parity,
883            "points": points,
884        },
885        "options": {
886            "returnInputData": False,
887            "returnPath": False,
888            "returnDiscriminatorPath": False,
889            # "fitEngine": "AnnealingFitter@2.0", #comment out to use the default fitter
890            # "fitObjective": "MB2@2.0",
891            "preferredMilkUnit": milk_unit,
892        },
893    }
894
895    # -----------------------------
896    # Add priors if provided or when using Chen et al. priors
897    # -----------------------------
898    if custom_priors == "CHEN":
899        assert parity is not None, "parity is required for Chen priors"
900        payload["priors"] = get_chen_priors(parity)
901
902    elif isinstance(custom_priors, dict):
903        payload["priors"] = dict(custom_priors)
904
905    # -----------------------------
906    # Call API
907    # -----------------------------
908    response = requests.post(f"{base_url}/fitLactation", headers=headers, json=payload)
909    response.raise_for_status()
910    res = response.json()
911
912    # -----------------------------
913    # Normalize response (USA vs EU)
914    # -----------------------------
915    if "fittedParams" in res:
916        # USA-style response
917        fitted = res["fittedParams"]
918    elif "params" in res:
919        fitted = res["params"]
920    else:
921        raise RuntimeError(f"Unexpected MilkBot response format: {res}")
922
923    # -----------------------------
924    # Parse result
925    # -----------------------------
926
927    return {
928        "scale": fitted["scale"],
929        "ramp": fitted["ramp"],
930        "decay": fitted["decay"],
931        "offset": fitted["offset"],
932        "nPoints": len(points),
933    }

191def brody_model(t, a, k) -> float:
192    """First Brody lactation curve model.
193
194    Args:
195        t: Time since calving in days (DIM).
196        a: Scale parameter (numerical).
197        k: Decay parameter (numerical).
198
199    Returns:
200        Predicted milk yield at `t`.
201
202    Notes:
203        Formula: `y(t) = a * exp(-k * t)`.
204
205        This was the first lactation curve model ever developed in 1923.
206    """
207    return a * np.exp(-k * t)

246def dhanoa_model(t, a, b, c) -> float:
247    """Dhanoa lactation curve model.
248
249    Args:
250        t: Time since calving in days (DIM).
251        a: Scale parameter (numerical).
252        b: Shape parameter (numerical).
253        c: Decay parameter (numerical).
254
255    Returns:
256        Predicted milk yield at `t`.
257
258    Notes:
259        Formula: `y(t) = a * t ** (b * c) * exp(-c * t)`.
260    """
261    return a * t ** (b * c) * np.exp(-c * t)

321def dijkstra_model(t, a, b, c, d) -> float:
322    """Dijkstra lactation curve model.
323
324    Args:
325        t: Time since calving in days (DIM).
326        a: Scale parameter (numerical).
327        b: Growth parameter (numerical).
328        c: Saturation rate parameter (numerical).
329        d: Decay parameter (numerical).
330
331    Returns:
332        Predicted milk yield at `t`.
333
334    Notes:
335        Formula: `y(t) = a * exp((b * (1 - exp(-c * t)) / c) - d * t)`.
336    """
337    return a * np.exp((b * (1 - np.exp(-c * t)) / c) - d * t)

264def emmans_model(t, a, b, c, d) -> float:
265    """Emmans lactation curve model.
266
267    Args:
268        t: Time since calving in days (DIM).
269        a: Scale parameter (numerical).
270        b: Growth parameter (numerical).
271        c: Decay parameter (numerical).
272        d: Location parameter in nested exponential (numerical).
273
274    Returns:
275        Predicted milk yield at `t`.
276
277    Notes:
278        Formula: `y(t) = a * exp(-exp(d - b*t)) * exp(-c*t)`.
279    """
280    return a * np.exp(-np.exp(d - b * t)) * np.exp(-c * t)

173def fischer_model(t, a, b, c) -> np.floating | np.ndarray:
174    """Fischer lactation curve model.
175
176    Args:
177        t: Time since calving in days (DIM).
178        a: Fitted coefficient (numerical).
179        b: Fitted coefficient (numerical).
180        c: Fitted coefficient (numerical).
181
182    Returns:
183        Predicted milk yield at `t`.
184
185    Notes:
186        Formula: `y(t) = a - b * t - a * exp(-c * t)`.
187    """
188    return a - b * t - a * np.exp(-c * t)

402def fit_lactation_curve(
403    dim,
404    milkrecordings,
405    model="wood",
406    fitting="frequentist",
407    breed="H",
408    parity=3,
409    continent="USA",
410    custom_priors=None,
411    key=None,
412    milk_unit="kg",
413) -> np.ndarray:
414    """Fit lactation data to a lactation curve model and return predictions.
415
416    Depending on `fitting`:
417    - **frequentist**: Fits parameters using `minimize` and/or `curve_fit`
418      for the specified `model`, then predicts over DIM 1–305 (or up to `max(dim)` if greater).
419    - **bayesian**: (MilkBot only) Calls the MilkBot Bayesian fitting API and
420      returns predictions using the fitted parameters.
421
422    Args:
423        dim (Int): List/array of days in milk (DIM).
424        milkrecordings (Float): List/array of milk recordings (kg).
425        model (Str): Model name (lowercase), default "wood".
426            Supported for frequentist: "wood", "wilmink", "ali_schaeffer", "fischer", "milkbot".
427        fitting (Str): "frequentist" (default) or "bayesian".
428            Bayesian fitting is currently implemented only for "milkbot".
429        breed (Str): "H" (Holstein, default) or "J" (Jersey). Only used for Bayesian.
430        parity (Int): Lactation number; all parities >= 3 considered one group in priors.
431            Only used for Bayesian.
432        continent (Str): priors chosen by MilkBot API based on continent averages.
433            Only used for Bayesian, options: "USA" (default) and "EU".
434        custom_priors (Dict | str | None): Custom prior
435            distributions for Bayesian fitting.
436            If a dict is provided, it must be a dictionary
437            of prior distributions for each parameter
438            in the model.
439            Set the correct dictionary using the `build_prior` helper function.
440            If the string "CHEN" is provided, the default Chen et al. priors are used.
441            Only used for Bayesian.
442        key = Str: API key for MilkBot API (required for Bayesian fitting).
443            Only used for Bayesian.
444        milk_unit (Str): Unit of milk yield measurements. Must be either "kg" or "lbs".
445            Default is "kg".
446            Only used for Bayesian.
447
448
449    Returns:
450        List/array of predicted milk yield for DIM 1–305 (or up to the maximum DIM if > 305).
451
452    Raises:
453        Exception: If an unknown model is requested (frequentist),
454            or Bayesian is requested for a non-MilkBot model,
455            or `key` is missing when Bayesian fitting is requested.
456
457    Notes:
458        Uses `validate_and_prepare_inputs` for input checking and normalization.
459    """
460    # check and prepare input
461    inputs = validate_and_prepare_inputs(
462        dim,
463        milkrecordings,
464        model=model,
465        fitting=fitting,
466        breed=breed,
467        parity=parity,
468        continent=continent,
469        custom_priors=custom_priors,
470        milk_unit=milk_unit,
471    )
472
473    dim = inputs.dim
474    milkrecordings = inputs.milkrecordings
475    model = inputs.model
476    fitting = inputs.fitting
477    breed = inputs.breed
478    parity = inputs.parity
479    continent = inputs.continent
480    custom_priors = inputs.custom_priors
481    milk_unit = inputs.milk_unit
482
483    if fitting == "frequentist":
484        if model == "wood":
485            params = get_lc_parameters(dim, milkrecordings, model)
486            assert params is not None, "Failed to fit Wood model parameters"
487            a_w, b_w, c_w = params[0], params[1], params[2]
488            if max(dim) > 305:
489                t_range = np.arange(1, (max(dim) + 1))
490                y_w = wood_model(t_range, a_w, b_w, c_w)
491            else:
492                t_range = np.arange(1, 306)
493                y_w = wood_model(t_range, a_w, b_w, c_w)
494            return np.asarray(y_w)
495
496        elif model == "wilmink":
497            params = get_lc_parameters(dim, milkrecordings, model)
498            assert params is not None, "Failed to fit Wilmink model parameters"
499            a_wil, b_wil, c_wil, k_wil = (params[0], params[1], params[2], params[3])
500            if max(dim) > 305:
501                t_range = np.arange(1, (max(dim) + 1))
502                y_wil = wilmink_model(t_range, a_wil, b_wil, c_wil, k_wil)
503            else:
504                t_range = np.arange(1, 306)
505                y_wil = wilmink_model(t_range, a_wil, b_wil, c_wil, k_wil)
506            return np.asarray(y_wil)
507
508        elif model == "ali_schaeffer":
509            params = get_lc_parameters(dim, milkrecordings, model)
510            assert params is not None, "Failed to fit Ali & Schaeffer model parameters"
511            a_as, b_as, c_as, d_as, k_as = (params[0], params[1], params[2], params[3], params[4])
512            if max(dim) > 305:
513                t_range = np.arange(1, (max(dim) + 1))
514                y_as = ali_schaeffer_model(t_range, a_as, b_as, c_as, d_as, k_as)
515            else:
516                t_range = np.arange(1, 306)
517                y_as = ali_schaeffer_model(t_range, a_as, b_as, c_as, d_as, k_as)
518            return np.asarray(y_as)
519
520        elif model == "fischer":
521            params = get_lc_parameters(dim, milkrecordings, model)
522            assert params is not None, "Failed to fit Fischer model parameters"
523            a_f, b_f, c_f = params[0], params[1], params[2]
524            if max(dim) > 305:
525                t_range = np.arange(1, (max(dim) + 1))
526                y_f = fischer_model(t_range, a_f, b_f, c_f)
527            else:
528                t_range = np.arange(1, 306)
529                y_f = fischer_model(t_range, a_f, b_f, c_f)
530            return np.asarray(y_f)
531
532        elif model == "milkbot":
533            params = get_lc_parameters(dim, milkrecordings, model)
534            assert params is not None, "Failed to fit MilkBot model parameters"
535            a_mb, b_mb, c_mb, d_mb = (params[0], params[1], params[2], params[3])
536            if max(dim) > 305:
537                t_range = np.arange(1, (max(dim) + 1))
538            else:
539                t_range = np.arange(1, 306)
540
541            y_mb = milkbot_model(t_range, a_mb, b_mb, c_mb, d_mb)
542            return np.asarray(y_mb)
543
544        else:
545            raise Exception("Unknown model")
546    else:
547        if model == "milkbot":
548            if key is None:
549                raise Exception("Key needed to use Bayesian fitting engine milkbot")
550            else:
551                assert parity is not None, "parity is required for Bayesian fitting"
552                assert breed is not None, "breed is required for Bayesian fitting"
553                assert continent is not None, "continent is required for Bayesian fitting"
554                parameters = bayesian_fit_milkbot_single_lactation(
555                    dim,
556                    milkrecordings,
557                    key,
558                    parity,
559                    breed,
560                    custom_priors,
561                    continent,
562                    milk_unit or "kg",
563                )
564                if max(dim) > 305:
565                    t_range = np.arange(1, (max(dim) + 1))
566                    y_mb_bay = milkbot_model(
567                        t_range,
568                        parameters["scale"],
569                        parameters["ramp"],
570                        parameters["offset"],
571                        parameters["decay"],
572                    )
573                else:
574                    t_range = np.arange(1, 306)
575                    y_mb_bay = milkbot_model(
576                        t_range,
577                        parameters["scale"],
578                        parameters["ramp"],
579                        parameters["offset"],
580                        parameters["decay"],
581                    )
582                return np.asarray(y_mb_bay)
583        else:
584            raise Exception("Bayesian fitting is currently only implemented for milkbot models")

723def get_chen_priors(parity: int) -> dict:
724    """
725    Return Chen et al. priors in MilkBot format.
726
727    Args:
728        parity: Lactation number (1, 2, or >= 3).
729
730    Returns:
731        Dictionary with parameter priors:
732        - "scale": {"mean", "sd"}
733        - "ramp": {"mean", "sd"}
734        - "decay": {"mean", "sd"}
735        - "offset": {"mean", "sd"}
736        - "seMilk": Standard error of milk measurement.
737        - "milkUnit": Unit string (e.g., "kg").
738    """
739    if parity == 1:
740        return {
741            "scale": {"mean": 34.11, "sd": 7},
742            "ramp": {"mean": 29.96, "sd": 3},
743            "decay": {"mean": 0.001835, "sd": 0.000738},
744            "offset": {"mean": -0.5, "sd": 0.02},
745            "seMilk": 4,
746            "milkUnit": "kg",
747        }
748
749    if parity == 2:
750        return {
751            "scale": {"mean": 44.26, "sd": 9.57},
752            "ramp": {"mean": 22.52, "sd": 3},
753            "decay": {"mean": 0.002745, "sd": 0.000979},
754            "offset": {"mean": -0.78, "sd": 0.07},
755            "seMilk": 4,
756            "milkUnit": "kg",
757        }
758
759    # parity >= 3
760    return {
761        "scale": {"mean": 48.41, "sd": 10.66},
762        "ramp": {"mean": 22.54, "sd": 8.724},
763        "decay": {"mean": 0.002997, "sd": 0.000972},
764        "offset": {"mean": 0.0, "sd": 0.03},
765        "seMilk": 4,
766        "milkUnit": "kg",
767    }

650def get_lc_parameters(dim, milkrecordings, model="wood") -> tuple[float, ...]:
651    """Fit lactation data to a model and return fitted parameters (frequentist).
652
653    Depending on `model`, this uses `scipy.optimize.minimize` and/or
654    `scipy.optimize.curve_fit` with model-specific starting values and bounds.
655
656    Args:
657        dim (int): List/array of DIM values.
658        milkrecordings (float): List/array of milk recordings (kg).
659        model (str): One of "wood", "wilmink", "ali_schaeffer", "fischer", "milkbot".
660
661    Returns:
662        Fitted parameters as floats, in alphabetical order by parameter name:
663            - wood: (a, b, c)
664            - wilmink: (a, b, c, k) with k fixed at -0.05
665            - ali_schaeffer: (a, b, c, d, k)
666            - fischer: (a, b, c)
667            - milkbot: (a, b, c, d)
668    """
669    # check and prepare input
670    inputs = validate_and_prepare_inputs(dim, milkrecordings, model=model)
671
672    dim = inputs.dim
673    milkrecordings = inputs.milkrecordings
674    model = inputs.model
675
676    if model == "wood":
677        wood_guess = [30, 0.2, 0.01]
678        wood_bounds = [(1, 100), (0.01, 1.5), (0.0001, 0.1)]
679        wood_res = minimize(
680            wood_objective, wood_guess, args=(dim, milkrecordings), bounds=wood_bounds
681        )
682        a_w, b_w, c_w = wood_res.x
683        return a_w, b_w, c_w
684
685    elif model == "wilmink":
686        wil_guess = [10, 0.1, 30]
687        wil_params, _ = curve_fit(wilmink_model, dim, milkrecordings, p0=wil_guess)
688        a_wil, b_wil, c_wil = wil_params
689        k_wil = -0.05  # set fixed
690        return a_wil, b_wil, c_wil, k_wil
691
692    elif model == "ali_schaeffer":
693        ali_schaeffer_guess = [10, 10, -5, 1, 1]
694        ali_schaeffer_params, _ = curve_fit(
695            ali_schaeffer_model, dim, milkrecordings, p0=ali_schaeffer_guess
696        )
697        a_as, b_as, c_as, d_as, k_as = ali_schaeffer_params
698        return a_as, b_as, c_as, d_as, k_as
699
700    elif model == "fischer":
701        fischer_guess = [max(milkrecordings), 0.01, 0.01]
702        fischer_bounds = [(0, 100), (0, 1), (0.0001, 1)]
703        fischer_params, _ = curve_fit(
704            fischer_model,
705            dim,
706            milkrecordings,
707            p0=fischer_guess,
708            bounds=np.transpose(fischer_bounds),
709        )
710        a_f, b_f, c_f = fischer_params
711        return a_f, b_f, c_f
712
713    elif model == "milkbot":
714        mb_guess = [max(milkrecordings), 20.0, -0.7, 0.022]
715        mb_bounds = [(1, 100), (1, 100), (-600, 300), (0.0001, 0.1)]
716        mb_res = minimize(milkbot_objective, mb_guess, args=(dim, milkrecordings), bounds=mb_bounds)
717        a_mb, b_mb, c_mb, d_mb = mb_res.x
718        return a_mb, b_mb, c_mb, d_mb
719
720    raise ValueError(f"Unknown model: {model}")

587def get_lc_parameters_least_squares(
588    dim, milkrecordings, model="milkbot"
589) -> tuple[float, float, float, float]:
590    """Fit lactation data and return model parameters (least squares; frequentist).
591
592    This helper uses `scipy.optimize.least_squares` to fit the MilkBot model with bounds,
593    and returns the fitted parameters.
594    Currently implemented only for the MilkBot model, as it is
595    more complex and benefits from the robust optimization approach.
596    Other models can be fitted using `get_lc_parameters` with
597    numerical optimisation, which is generally faster for simpler
598    models.
599
600    Args:
601        dim (int): List/array of DIM values.
602        milkrecordings (float): List/array of milk recordings (kg).
603        model (str): Pre-defined model name; currently used with "milkbot".
604
605    Returns:
606        Parameters `(a, b, c, d)` as `np.float` in alphabetic order.
607
608    """
609    # check and prepare input
610    inputs = validate_and_prepare_inputs(dim, milkrecordings, model=model)
611
612    dim = inputs.dim
613    milkrecordings = inputs.milkrecordings
614    model = inputs.model
615
616    # ------------------------------
617    # Initial guess
618    # ------------------------------
619    a0 = np.max(milkrecordings)
620    b0 = 50.0
621    c0 = 30.0
622    d0 = 0.01
623    p0 = [a0, b0, c0, d0]
624
625    # ------------------------------
626    # Parameter bounds
627    # ------------------------------
628    lower = [np.max(milkrecordings) * 0.5, 1.0, -300.0, 1e-6]
629    upper = [np.max(milkrecordings) * 8.0, 400.0, 300.0, 1.0]
630
631    # ------------------------------
632    # Fit using least-squares
633    # ------------------------------
634    res = least_squares(
635        residuals_milkbot,
636        p0,
637        args=(dim, milkrecordings),
638        bounds=(lower, upper),
639        method="trf",  # trust region reflective, works well with bounds
640    )
641
642    # ------------------------------
643    # Extract parameters
644    # ------------------------------
645    a_mb, b_mb, c_mb, d_mb = res.x
646
647    return a_mb, b_mb, c_mb, d_mb

283def hayashi_model(t, a, b, c, d) -> float:
284    """Hayashi lactation curve model.
285
286    Args:
287        t: Time since calving in days (DIM).
288        a: Ratio parameter (> 0) (numerical).
289        b: Scale parameter (numerical).
290        c: Time constant for the first exponential term (numerical).
291        d: Parameter retained for compatibility with literature (unused in this expression).
292
293    Returns:
294        Predicted milk yield at `t`.
295
296    Notes:
297        Formula: `y(t) = b * (exp(-t / c) - exp(-t / (a * c)))`.
298    """
299    return b * (np.exp(-t / c) - np.exp(-t / (a * c)))

 93def milkbot_model(t, a, b, c, d) -> np.floating | np.ndarray:
 94    """MilkBot lactation curve model.
 95
 96    Args:
 97        t: Time since calving in days (DIM), scalar or array-like.
 98        a: Scale; overall level of milk production.
 99        b: Ramp; governs the rate of rise in early lactation.
100        c: Offset; small (usually minor) correction around the theoretical start of lactation.
101        d: Decay; exponential decline rate, evident in late lactation.
102
103    Returns:
104        Predicted milk yield at `t` (same shape as `t`).
105
106    Notes:
107        Formula: `y(t) = a * (1 - exp((c - t) / b) / 2) * exp(-d * t)`.
108    """
109    return a * (1 - np.exp((c - t) / b) / 2) * np.exp(-d * t)

228def nelder_model(t, a, b, c) -> float:
229    """Nelder lactation curve model.
230
231    Args:
232        t: Time since calving in days (DIM).
233        a: Denominator intercept (numerical).
234        b: Denominator linear coefficient (numerical).
235        c: Denominator quadratic coefficient (numerical).
236
237    Returns:
238        Predicted milk yield at `t`.
239
240    Notes:
241        Formula: `y(t) = t / (a + b*t + c*t^2)`.
242    """
243    return t / (a + b * t + c * t**2)

340def prasad_model(t, a, b, c, d) -> float:
341    """Prasad lactation curve model.
342
343    Args:
344        t: Time since calving in days (DIM).
345        a: Intercept-like parameter (numerical).
346        b: Linear coefficient (numerical).
347        c: Quadratic coefficient (numerical).
348        d: Inverse-time coefficient (numerical).
349
350    Returns:
351        Predicted milk yield at `t`.
352
353    Notes:
354        Formula: `y(t) = a + b*t + c*t^2 + d/t`.
355    """
356    return a + b * t + c * t**2 + d / t

302def rook_model(t, a, b, c, d) -> float:
303    """Rook lactation curve model.
304
305    Args:
306        t: Time since calving in days (DIM).
307        a: Scale parameter (numerical).
308        b: Shape parameter in rational term (numerical).
309        c: Offset parameter in rational term (numerical).
310        d: Exponential decay parameter (numerical).
311
312    Returns:
313        Predicted milk yield at `t`.
314
315    Notes:
316        Formula: `y(t) = a * (1 / (1 + b / (c + t))) * exp(-d * t)`.
317    """
318    return a * (1 / (1 + b / (c + t))) * np.exp(-d * t)

210def sikka_model(t, a, b, c) -> float:
211    """Sikka lactation curve model.
212
213    Args:
214        t: Time since calving in days (DIM).
215        a: Scale parameter (numerical).
216        b: Growth parameter (numerical).
217        c: Quadratic decay parameter (numerical).
218
219    Returns:
220        Predicted milk yield at `t`.
221
222    Notes:
223        Formula: `y(t) = a * exp(b * t - c * t^2)`.
224    """
225    return a * np.exp(b * t - c * t**2)

130def wilmink_model(t, a, b, c, k=-0.05) -> np.floating | np.ndarray:
131    """Wilmink lactation curve model.
132
133    Args:
134        t: Time since calving in days (DIM), scalar or array-like.
135        a: Fitted coefficient (numerical).
136        b: Fitted coefficient (numerical).
137        c: Fitted coefficient (numerical).
138        k: Fixed exponential coefficient (numerical), default -0.05.
139
140    Returns:
141        Predicted milk yield at `t`.
142
143    Notes:
144        Formula: `y(t) = a + b * t + c * exp(k * t)`.
145    """
146    t = np.asarray(t)
147    return a + b * t + c * np.exp(k * t)

112def wood_model(t, a, b, c) -> np.floating | np.ndarray:
113    """Wood lactation curve model.
114
115    Args:
116        t: Time since calving in days (DIM), scalar or array-like.
117        a: Fitted coefficient (numerical).
118        b: Fitted coefficient (numerical).
119        c: Fitted coefficient (numerical).
120
121    Returns:
122        Predicted milk yield at `t`.
123
124    Notes:
125        Formula: `y(t) = a * t^b * exp(-c * t)`.
126    """
127    return a * (t**b) * np.exp(-c * t)

770def build_prior(
771    scale_mean: float,
772    scale_sd: float,
773    ramp_mean: float,
774    ramp_sd: float,
775    decay_mean: float,
776    decay_sd: float,
777    offset_mean: float,
778    offset_sd: float,
779    se_milk: float = 4,
780) -> dict:
781    return {
782        "scale": {"mean": scale_mean, "sd": scale_sd},
783        "ramp": {"mean": ramp_mean, "sd": ramp_sd},
784        "decay": {"mean": decay_mean, "sd": decay_sd},
785        "offset": {"mean": offset_mean, "sd": offset_sd},
786        "seMilk": se_milk,
787    }