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)

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

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

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

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

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

782def build_prior(
783    scale_mean: float,
784    scale_sd: float,
785    ramp_mean: float,
786    ramp_sd: float,
787    decay_mean: float,
788    decay_sd: float,
789    offset_mean: float,
790    offset_sd: float,
791    se_milk: float = 4,
792) -> dict:
793    return {
794        "scale": {"mean": scale_mean, "sd": scale_sd},
795        "ramp": {"mean": ramp_mean, "sd": ramp_sd},
796        "decay": {"mean": decay_mean, "sd": decay_sd},
797        "offset": {"mean": offset_mean, "sd": offset_sd},
798        "seMilk": se_milk,
799    }