From bb6b8d7e1d3593a6f43a9806af560d45b4ec8326 Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Mon, 24 Aug 2026 15:55:59 -0400 Subject: [PATCH 01/13] draft --- RUFAS/biophysical/animal/animal.py | 225 +++++++++++------- RUFAS/biophysical/animal/animal_config.py | 172 +------------ RUFAS/biophysical/animal/animal_constants.py | 50 +++- .../animal/animal_module_reporter.py | 7 +- .../animal/data_types/herd_statistics.py | 14 +- RUFAS/input/metadata/properties/default.json | 143 +---------- .../data/animal/example_freestall_animal.json | 28 +-- .../data/animal/example_open_lot_animal.json | 28 +-- .../animal_cross_validation.json | 176 -------------- .../test_animal/test_animal/test_animal.py | 104 ++++---- .../test_animal/test_animal_config.py | 56 +---- .../test_animal_module_reporter.py | 4 +- .../test_data_types/test_herd_statistics.py | 2 +- .../test_herd_manager/pytest_fixtures.py | 28 +-- .../test_herd_manager_herd_statistics.py | 98 ++------ 15 files changed, 269 insertions(+), 866 deletions(-) diff --git a/RUFAS/biophysical/animal/animal.py b/RUFAS/biophysical/animal/animal.py index b60edcc87e..955ca7a1d1 100644 --- a/RUFAS/biophysical/animal/animal.py +++ b/RUFAS/biophysical/animal/animal.py @@ -1,6 +1,6 @@ import sys from datetime import timedelta -from random import random, randint +from random import random, randint, uniform from typing import Callable, cast from scipy.stats import truncnorm @@ -1767,9 +1767,12 @@ def daily_reproduction_update( self.milk_production.set_wood_parameters( wood_parameters["l"], wood_parameters["m"], wood_parameters["n"] ) - self.future_death_date = self.determine_future_death_date() - self._future_death_reason = animal_constants.DEATH_CULL - self.future_cull_date, self.cull_reason = self.determine_future_cull_date() + # A heifer becomes a cow at her first calving; give her an initial mortality / + # acute-sale assessment now so she is not risk-free until the next Jan 1. Later + # lactations are not re-rolled here -- removal risk is reassessed annually + # (see ``_assess_annual_removal_risk``), decoupling it from reproduction. + if self.calves == 1: + self._assess_annual_removal_risk() self.events += reproduction_outputs.events @@ -1807,6 +1810,13 @@ def daily_routines(self, time: RufasTime) -> DailyRoutinesOutput: newborn_calf_config, daily_routines_output.herd_reproduction_statistics = self.daily_reproduction_update(time) + # Reassess mortality / acute-sale risk once a year (Jan 1, Julian day 1) for every cow, so + # risk accrues with time in the herd rather than only at calving (issue #2694). The guard + # inside the assessment skips risk types that already have a pending event, so a cow that + # first calved earlier today is not double-rolled. + if self.animal_type.is_cow and time.current_julian_day == 1: + self._assess_annual_removal_risk() + daily_routines_output.animal_status, daily_routines_output.newborn_calf_config = self.animal_life_stage_update( time ) @@ -2422,110 +2432,161 @@ def _get_cow_values(self) -> CowValuesTypedDict: parity=self.calves, ) + def _assess_annual_removal_risk(self) -> None: + """ + Roll a cow's annual mortality and acute-sale risk and schedule any resulting removal. + + Called on each Jan 1 for every cow, and once when a heifer first calves, so removal risk + accrues with time spent in the herd rather than only at calving (issue #2694). Death and + acute sale are rolled independently. A risk type that already has a pending future event is + left untouched: this prevents overwriting or double-scheduling an event a prior roll placed + (the days-in-milk timing can land a scheduled event more than a year out, across a Jan 1). + + Notes + ------- + [AN.ANM.1], [AN.ANM.2] + + """ + if self.future_death_date == sys.maxsize: + death_date = self.determine_future_death_date() + if death_date != sys.maxsize: + self.future_death_date = death_date + self._future_death_reason = animal_constants.DEATH_CULL + + if self.future_cull_date == sys.maxsize: + cull_date, cull_reason = self.determine_future_cull_date() + if cull_date != sys.maxsize: + self.future_cull_date = cull_date + self.cull_reason = cull_reason + + def _parity_index(self) -> int: + """Return the 0-based index into a by-parity array, capping parity 4+ at the last entry.""" + return 3 if self.calves >= 4 else self.calves - 1 + + @staticmethod + def _interpolate_cdf_at_day(day: float, cdf: list[float], breakpoints: list[int]) -> float: + """ + Linearly interpolate a cumulative-distribution value at a given day in milk. + + Parameters + ---------- + day : float + Day in milk at which to evaluate the CDF. + cdf : list[float] + Cumulative-distribution values at each breakpoint (non-decreasing, 0.0 to 1.0). + breakpoints : list[int] + Days-in-milk breakpoints that partition the CDF (same length as ``cdf``). + + Returns + ------- + float + The interpolated CDF value, clamped to ``cdf[0]`` below the first breakpoint and + ``cdf[-1]`` at or beyond the last. + + """ + if day <= breakpoints[0]: + return cdf[0] + for i in range(len(breakpoints) - 1): + if breakpoints[i] <= day < breakpoints[i + 1]: + slope = (cdf[i + 1] - cdf[i]) / (breakpoints[i + 1] - breakpoints[i]) + return cdf[i] + slope * (day - breakpoints[i]) + return cdf[-1] + + def _sample_removal_date(self, timing_cdf: list[float]) -> int: + """ + Sample an absolute simulation day for a scheduled removal from a days-in-milk timing CDF. + + The day is shaped by ``timing_cdf`` (the distribution of removal timing across a lactation) + but conditioned to fall after the cow's current days in milk, so a cow selected for removal + is always scheduled to leave in the future rather than "escaping" the event. If the cow is + already past the CDF's last breakpoint (an extended lactation), the event is instead placed + uniformly within ``REMOVAL_FALLBACK_WINDOW_DAYS`` of today. + + Parameters + ---------- + timing_cdf : list[float] + Cumulative-distribution values of removal timing at each of + ``REMOVAL_TIMING_DAY_BREAKPOINTS``. + + Returns + ------- + int + The absolute simulation day (in ``days_born`` terms) on which the removal occurs. + + Notes + ------- + [AN.ANM.1], [AN.ANM.2] + + """ + breakpoints = animal_constants.REMOVAL_TIMING_DAY_BREAKPOINTS + current_days_in_milk = self.days_in_milk + lactation_start = self.days_born - current_days_in_milk + + if current_days_in_milk >= breakpoints[-1]: + return self.days_born + randint(1, animal_constants.REMOVAL_FALLBACK_WINDOW_DAYS) + + # Draw uniformly on the CDF mass that remains after the current day in milk, then invert + # back to a day in milk so the timing keeps its lactation-stage shape. + lower_cdf_value = self._interpolate_cdf_at_day(current_days_in_milk, timing_cdf, breakpoints) + removal_cdf_value = uniform(lower_cdf_value, 1.0) + for i in range(len(timing_cdf) - 1): + if timing_cdf[i] <= removal_cdf_value < timing_cdf[i + 1]: + slope = (breakpoints[i + 1] - breakpoints[i]) / (timing_cdf[i + 1] - timing_cdf[i]) + day_in_milk = breakpoints[i] + slope * (removal_cdf_value - timing_cdf[i]) + return round(lactation_start + day_in_milk) + # ``uniform`` can return its upper bound (1.0), which no half-open segment contains. + return round(lactation_start + breakpoints[-1]) + def determine_future_death_date(self) -> int: """ - Determine the future death date of the animal based on its parity. + Roll the cow's annual mortality risk and, if selected, schedule a death day. + + The parity-indexed :attr:`AnimalConfig.parity_death_probability` is now an *annual* + probability (issue #2694), rolled by the annual removal assessment rather than once per + lactation. When the cow is selected to die, the day is drawn from the days-in-milk death + timing CDF via :meth:`_sample_removal_date`. Returns ------- int - Calculated future death date in simulation days. + Calculated future death date in simulation days, or ``sys.maxsize`` if not selected. Notes ------- [AN.ANM.1] """ - if self.calves >= 4: - death_rate = AnimalConfig.parity_death_probability[3] - else: - death_rate = AnimalConfig.parity_death_probability[self.calves - 1] - death_rand = random() - if death_rand <= death_rate: - death_probability_upper_limit = death_probability_lower_limit = 0.0 - death_time_upper_limit = death_time_lower_limit = 0.0 - death_date_random = random() - for i in range(len(AnimalConfig.death_day_probability) - 1): - if ( - AnimalConfig.death_day_probability[i] - <= death_date_random - < AnimalConfig.death_day_probability[i + 1] - ): - death_probability_lower_limit = AnimalConfig.death_day_probability[i] - death_probability_upper_limit = AnimalConfig.death_day_probability[i + 1] - death_time_lower_limit = AnimalConfig.cull_day_count[i] - death_time_upper_limit = AnimalConfig.cull_day_count[i + 1] - n = (death_time_upper_limit - death_time_lower_limit) / ( - death_probability_upper_limit - death_probability_lower_limit - ) - return round( - death_time_lower_limit + n * (death_date_random - death_probability_lower_limit) + self.days_born - ) + death_rate = AnimalConfig.parity_death_probability[self._parity_index()] + if random() <= death_rate: + return self._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) return sys.maxsize def determine_future_cull_date(self) -> tuple[int, str]: """ - Determine the future cull date and reason for the animal based on parity-specific probabilities. + Roll the cow's annual acute-sale risk and, if selected, schedule an acute-sale day. + + An acute sale ("forced" / "involuntary" / "spontaneous" removal) is a cow that must leave + the herd immediately regardless of whether a replacement is available. The parity-indexed + :attr:`AnimalConfig.parity_acute_sale_probability` is now an *annual* probability + (issue #2694); the former six disease-specific reasons are collapsed into the single + :data:`animal_constants.ACUTE_SALE_CULL` reason with one timing CDF. Returns ------- tuple[int, str] - - Future cull date in simulation days. - - Reason for culling. + - Future acute-sale date in simulation days (``sys.maxsize`` if not selected). + - Reason for removal (empty string if not selected). Notes ------- [AN.ANM.2] """ - cull_reason = "" - future_cull_date = sys.maxsize - if self.calves >= 4: - inv_cull_rate = AnimalConfig.parity_cull_probability[3] - else: - inv_cull_rate = AnimalConfig.parity_cull_probability[self.calves - 1] - cull_rand = random() - if cull_rand <= inv_cull_rate: - cull_reason_rand = random() - cull_prob = 0.0 - if cull_reason_rand <= (cull_prob := cull_prob + AnimalConfig.feet_leg_cull_probability): - cull_reason_cull_prob = AnimalConfig.feet_leg_cull_day_probability - cull_reason = animal_constants.LAMENESS_CULL - - elif cull_reason_rand <= (cull_prob := cull_prob + AnimalConfig.injury_cull_probability): - cull_reason_cull_prob = AnimalConfig.injury_cull_day_probability - cull_reason = animal_constants.INJURY_CULL - - elif cull_reason_rand <= (cull_prob := cull_prob + AnimalConfig.mastitis_cull_probability): - cull_reason_cull_prob = AnimalConfig.mastitis_cull_day_probability - cull_reason = animal_constants.MASTITIS_CULL - - elif cull_reason_rand <= (cull_prob := cull_prob + AnimalConfig.disease_cull_probability): - cull_reason_cull_prob = AnimalConfig.disease_cull_day_probability - cull_reason = animal_constants.DISEASE_CULL - - elif cull_reason_rand <= (cull_prob + AnimalConfig.udder_cull_probability): - cull_reason_cull_prob = AnimalConfig.udder_cull_day_probability - cull_reason = animal_constants.UDDER_CULL - - else: - cull_reason_cull_prob = AnimalConfig.unknown_cull_day_probability - cull_reason = animal_constants.UNKNOWN_CULL - - cull_time_rand = random() - cull_reason_upper_limit = cull_reason_lower_limit = cull_time_upper_limit = cull_time_lower_limit = 0.0 - for i in range(len(cull_reason_cull_prob) - 1): - if cull_reason_cull_prob[i] <= cull_time_rand < cull_reason_cull_prob[i + 1]: - cull_reason_lower_limit = cull_reason_cull_prob[i] - cull_reason_upper_limit = cull_reason_cull_prob[i + 1] - cull_time_lower_limit = AnimalConfig.cull_day_count[i] - cull_time_upper_limit = AnimalConfig.cull_day_count[i + 1] - x = (cull_time_upper_limit - cull_time_lower_limit) / (cull_reason_upper_limit - cull_reason_lower_limit) - future_cull_date = round( - cull_time_lower_limit + x * (cull_time_rand - cull_reason_lower_limit) + self.days_born - ) - - return future_cull_date, cull_reason + acute_sale_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] + if random() <= acute_sale_rate: + future_cull_date = self._sample_removal_date(animal_constants.ACUTE_SALE_TIMING_DAY_PROBABILITY) + return future_cull_date, animal_constants.ACUTE_SALE_CULL + return sys.maxsize, "" def update_pen_history(self, current_pen: int, current_day: int, animal_types_in_pen: set[AnimalType]) -> None: """ diff --git a/RUFAS/biophysical/animal/animal_config.py b/RUFAS/biophysical/animal/animal_config.py index a11a7cde68..fb89aaf780 100644 --- a/RUFAS/biophysical/animal/animal_config.py +++ b/RUFAS/biophysical/animal/animal_config.py @@ -146,37 +146,10 @@ class AnimalConfig: third_pregnancy_check_loss_rate : float Pregnancy loss probability during the third pregnancy check, (unitless). parity_death_probability : list[float] - List of probabilities of death based on parity number, (unitless). - death_day_probability : list[float] - Cumulative probability of cow death as a function of days in production, (unitless). - parity_cull_probability : list[float] - List of culling probabilities based on parity number, (unitless). - cull_day_count : list[int] - List of day intervals for culling analysis, (simulation day). - feet_leg_cull_probability : float - Probability of feet and leg-related culling, (unitless). - feet_leg_cull_day_probability : list[float] - Feet and leg-related culling probability over time, (unitless). - injury_cull_probability : float - Probability of culling due to injuries, (unitless). - injury_cull_day_probability : list[float] - Cumulative distribution for injury-related culling over time, (unitless). - mastitis_cull_probability : float - Probability of culling due to mastitis, (unitless). - mastitis_cull_day_probability : list[float] - Cumulative distribution for mastitis-related culling over time, (unitless). - disease_cull_probability : float - Probability of culling due to diseases, (unitless). - disease_cull_day_probability : list[float] - Cumulative distribution for disease-related culling over time, (unitless). - udder_cull_probability : float - Probability of culling due to udder-related issues, (unitless). - udder_cull_day_probability : list[float] - Cumulative distribution for udder-related culling over time, (unitless). - unknown_cull_probability : float - Probability of culling for unknown reasons, (unitless). - unknown_cull_day_probability : list[float] - Cumulative distribution for unknown reasons of culling over time, (unitless). + Annual, parity-indexed probability that a cow dies during a given year, (unitless). + parity_acute_sale_probability : list[float] + Annual, parity-indexed probability that a cow is sold for an acute / involuntary + reason during a given year, (unitless). methane_mitigation_method : str The mitigation method applied for methane reduction, e.g., "None", (unitless). methane_mitigation_additive_amount : float @@ -267,113 +240,11 @@ class AnimalConfig: third_pregnancy_check_day: int = 200 third_pregnancy_check_loss_rate: float = 0.017 + # Annual, parity-indexed probabilities that a cow dies (parity_death_probability) or is sold + # for an acute / involuntary reason (parity_acute_sale_probability) during a given year. See + # issue #2694: these were formerly evaluated once per lactation and are now evaluated annually. parity_death_probability: list[float] = [0.039, 0.056, 0.085, 0.117] - death_day_probability: list[float] = [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1] - - parity_cull_probability: list[float] = [0.169, 0.233, 0.301, 0.408] - cull_day_count: list[int] = [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530] - feet_leg_cull_probability: float = 0.1633 - feet_leg_cull_day_probability: list[float] = [ - 0, - 0.03, - 0.08, - 0.16, - 0.25, - 0.36, - 0.48, - 0.59, - 0.69, - 0.78, - 0.85, - 0.90, - 0.95, - 1, - ] - injury_cull_probability: float = 0.2883 - injury_cull_day_probability: list[float] = [ - 0, - 0.08, - 0.18, - 0.28, - 0.38, - 0.47, - 0.56, - 0.64, - 0.71, - 0.78, - 0.85, - 0.90, - 0.95, - 1, - ] - mastitis_cull_probability: float = 0.2439 - mastitis_cull_day_probability: list[float] = [ - 0, - 0.06, - 0.12, - 0.19, - 0.30, - 0.43, - 0.56, - 0.68, - 0.78, - 0.85, - 0.90, - 0.94, - 0.97, - 1, - ] - disease_cull_probability: float = 0.1391 - disease_cull_day_probability: list[float] = [ - 0, - 0.04, - 0.12, - 0.24, - 0.34, - 0.42, - 0.50, - 0.57, - 0.64, - 0.72, - 0.81, - 0.89, - 0.95, - 1, - ] - udder_cull_probability: float = 0.0645 - udder_cull_day_probability: list[float] = [ - 0, - 0.12, - 0.24, - 0.33, - 0.41, - 0.48, - 0.55, - 0.62, - 0.68, - 0.76, - 0.82, - 0.89, - 0.95, - 1, - ] - unknown_cull_probability: float = 0.1009 - unknown_cull_day_probability: list[float] = [ - 0, - 0.05, - 0.11, - 0.18, - 0.27, - 0.37, - 0.45, - 0.54, - 0.62, - 0.70, - 0.77, - 0.84, - 0.92, - 1, - ] + parity_acute_sale_probability: list[float] = [0.169, 0.233, 0.301, 0.408] methane_model: dict[str, Any] = { "calves": "Pattanaik", @@ -515,32 +386,7 @@ def initialize_animal_config(cls) -> None: cls.third_pregnancy_check_loss_rate = animal_config_data["from_literature"]["repro"]["preg_loss_rate_3"] cls.parity_death_probability = animal_config_data["from_literature"]["culling"]["parity_death_prob"] - cls.death_day_probability = animal_config_data["from_literature"]["culling"]["death_day_prob"] - - cls.parity_cull_probability = animal_config_data["from_literature"]["culling"]["parity_cull_prob"] - cls.cull_day_count = animal_config_data["from_literature"]["culling"]["cull_day_count"] - cls.feet_leg_cull_probability = animal_config_data["from_literature"]["culling"]["feet_leg_cull"]["probability"] - cls.feet_leg_cull_day_probability = animal_config_data["from_literature"]["culling"]["feet_leg_cull"][ - "cull_day_prob" - ] - cls.injury_cull_probability = animal_config_data["from_literature"]["culling"]["injury_cull"]["probability"] - cls.injury_cull_day_probability = animal_config_data["from_literature"]["culling"]["injury_cull"][ - "cull_day_prob" - ] - cls.mastitis_cull_probability = animal_config_data["from_literature"]["culling"]["mastitis_cull"]["probability"] - cls.mastitis_cull_day_probability = animal_config_data["from_literature"]["culling"]["mastitis_cull"][ - "cull_day_prob" - ] - cls.disease_cull_probability = animal_config_data["from_literature"]["culling"]["disease_cull"]["probability"] - cls.disease_cull_day_probability = animal_config_data["from_literature"]["culling"]["disease_cull"][ - "cull_day_prob" - ] - cls.udder_cull_probability = animal_config_data["from_literature"]["culling"]["udder_cull"]["probability"] - cls.udder_cull_day_probability = animal_config_data["from_literature"]["culling"]["udder_cull"]["cull_day_prob"] - cls.unknown_cull_probability = animal_config_data["from_literature"]["culling"]["unknown_cull"]["probability"] - cls.unknown_cull_day_probability = animal_config_data["from_literature"]["culling"]["unknown_cull"][ - "cull_day_prob" - ] + cls.parity_acute_sale_probability = animal_config_data["from_literature"]["culling"]["parity_acute_sale_prob"] cls.methane_model = animal_data["methane_model"] methane_mitigation_data = animal_data["methane_mitigation"] diff --git a/RUFAS/biophysical/animal/animal_constants.py b/RUFAS/biophysical/animal/animal_constants.py index 3515ff06e5..6672112c0c 100644 --- a/RUFAS/biophysical/animal/animal_constants.py +++ b/RUFAS/biophysical/animal/animal_constants.py @@ -96,17 +96,55 @@ HEIFER_REPRO_CULL = "culled for heifer reproductive problem" OVERSUPPLY_CULL = "culled for herd resize" DEATH_CULL = "culled for death" -LAMENESS_CULL = "culled for lameness" -INJURY_CULL = "culled for injury" -MASTITIS_CULL = "culled for mastitis" -DISEASE_CULL = "culled for disease" -UDDER_CULL = "culled for udder" -UNKNOWN_CULL = "culled for unknown" +# A single "acute sale" (a.k.a. forced / involuntary / spontaneous) removal reason replaces the +# former six disease-specific cull reasons (feet-and-leg, injury, mastitis, disease, udder, +# unknown). These represent animals that must leave the herd immediately regardless of whether a +# replacement is available; the disease-level detail is not recoverable from farm records. +ACUTE_SALE_CULL = "sold for acute (involuntary) reason" # youngstock mortality (a loss from death, not a cull) CALF_MORTALITY_LOSS = "died from pre-wean mortality" HEIFER_MORTALITY_LOSS = "died from post-wean mortality" +# Removal-event timing distributions (moved out of user inputs; see issue #2694). +# +# These describe *when*, within a cow's lactation (days in milk), a scheduled death or acute +# sale occurs. They are model constants -- not recoverable from farm records -- so they live +# here rather than in animal.json. REMOVAL_TIMING_DAY_BREAKPOINTS are the days-in-milk +# breakpoints that partition each cumulative distribution function (CDF); the two probability +# arrays give the CDF value at each breakpoint (they must start at 0.0, end at 1.0, be +# non-decreasing, and match the breakpoints in length). +REMOVAL_TIMING_DAY_BREAKPOINTS = [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530] + +# CDF of death timing by days in milk. +DEATH_TIMING_DAY_PROBABILITY = [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1] + +# CDF of acute-sale timing by days in milk. This single curve replaces the six former +# reason-specific curves. The placeholder values below are the mean of those six curves; the +# final values should be confirmed via literature review / expert input (issue #2694). +ACUTE_SALE_TIMING_DAY_PROBABILITY = [ + 0, + 0.063, + 0.142, + 0.230, + 0.325, + 0.422, + 0.517, + 0.607, + 0.687, + 0.765, + 0.833, + 0.893, + 0.948, + 1, +] + +# When a cow selected for removal is already past the timing CDF's last breakpoint (e.g. an +# extended lactation with days in milk >= REMOVAL_TIMING_DAY_BREAKPOINTS[-1]), the day-in-milk +# curve has no remaining support, so the event is instead scheduled uniformly within this many +# days from the current day. +REMOVAL_FALLBACK_WINDOW_DAYS = 60 + # STATS STDI = 2 diff --git a/RUFAS/biophysical/animal/animal_module_reporter.py b/RUFAS/biophysical/animal/animal_module_reporter.py index 474fd6f8ab..995844703d 100644 --- a/RUFAS/biophysical/animal/animal_module_reporter.py +++ b/RUFAS/biophysical/animal/animal_module_reporter.py @@ -842,12 +842,7 @@ def report_herd_statistics_data(cls, herd_statistics: HerdStatistics, simulation cull_reason_stats_units = { animal_constants.DEATH_CULL: MeasurementUnits.UNITLESS, animal_constants.OVERSUPPLY_CULL: MeasurementUnits.UNITLESS, - animal_constants.LAMENESS_CULL: MeasurementUnits.UNITLESS, - animal_constants.INJURY_CULL: MeasurementUnits.UNITLESS, - animal_constants.MASTITIS_CULL: MeasurementUnits.UNITLESS, - animal_constants.DISEASE_CULL: MeasurementUnits.UNITLESS, - animal_constants.UDDER_CULL: MeasurementUnits.UNITLESS, - animal_constants.UNKNOWN_CULL: MeasurementUnits.UNITLESS, + animal_constants.ACUTE_SALE_CULL: MeasurementUnits.UNITLESS, } om.add_variable( "cull_reason_stats", diff --git a/RUFAS/biophysical/animal/data_types/herd_statistics.py b/RUFAS/biophysical/animal/data_types/herd_statistics.py index 718b2ee349..0859744fcd 100644 --- a/RUFAS/biophysical/animal/data_types/herd_statistics.py +++ b/RUFAS/biophysical/animal/data_types/herd_statistics.py @@ -267,12 +267,7 @@ def __init__(self) -> None: self.cull_reason_stats = { animal_constants.DEATH_CULL: 0, animal_constants.OVERSUPPLY_CULL: 0, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 0, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 0, } self.parity_culling_stats_range = {"1": 0, "2": 0, "3": 0, "4": 0, "5": 0, "greater_than_5": 0} self.num_cow_for_parity = {"1": 0, "2": 0, "3": 0, "4": 0, "5": 0, "greater_than_5": 0} @@ -281,12 +276,7 @@ def __init__(self) -> None: self.cull_reason_stats_percent = { animal_constants.DEATH_CULL: 0.0, animal_constants.OVERSUPPLY_CULL: 0.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 0.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 0.0, } self.percent_cow_for_parity = { "1": 0.0, diff --git a/RUFAS/input/metadata/properties/default.json b/RUFAS/input/metadata/properties/default.json index 3a4d2208b6..202b30b60a 100644 --- a/RUFAS/input/metadata/properties/default.json +++ b/RUFAS/input/metadata/properties/default.json @@ -621,154 +621,23 @@ }, "culling": { "type": "object", - "description": "Defines probabilities and distributions for death and six health-related reasons (feet-and-leg, injury, mastitis, disease, udder, unknown) an animal might be removed (culled) from the herd.", - "cull_day_count": { - "type": "array", - "description": "Defines breakpoints that partition the cumulative distribution function (CDF) for culling probabilities into segments. These values correspond to the 'cull_day_prob' array, allowing for more accurate definition of the CDF.The numbers in the array represent days into the lactation (days in milk).", - "properties": { - "type": "number", - "minimum": 0 - } - }, - "feet_leg_cull": { - "type": "object", - "description": "Cull probabilities due to feet-and-leg-related health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to feet-and-leg-related issues. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of feet-and-leg-related culling over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, - "injury_cull": { - "type": "object", - "description": "Cull probabilities due to injury-related health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to injury-related issues. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of injury-related culling over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, - "mastitis_cull": { - "type": "object", - "description": "Cull probabilities due to mastitis-related health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to mastitis-related issues. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of mastitis-related culling over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, - "disease_cull": { - "type": "object", - "description": "Cull probabilities due to disease-related health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to general disease-related issues. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of general disease-related culling over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, - "udder_cull": { - "type": "object", - "description": "Cull probabilities due to udder-related health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to udder-related issues. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of udder-related culling over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, - "unknown_cull": { - "type": "object", - "description": "Cull probabilities due to other health issues.", - "probability": { - "type": "number", - "description": "Conditional probability that a culled (sold) animal is removed due to unknown reasons. This probability is used to determine the reason for culling after it has been decided that the animal will be culled during the current lactation. The sum of probabilities for the six culling reasons equals 1.", - "minimum": 0, - "maximum": 1 - }, - "cull_day_prob": { - "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of culling for unknown reasons over time. The 'cull_day_count' array defines the segments of the CDF.", - "properties": { - "type": "number", - "minimum": 0, - "maximum": 1 - } - } - }, + "description": "Defines the annual, by-parity probabilities that a cow dies or is sold for an acute (involuntary) reason. The timing of each event within a lactation, previously configured here, is now a model constant (see animal_constants.py, issue #2694).", "parity_death_prob": { "type": "array", - "description": "Death Probability, by Parity", - "properties": { - "type": "number", - "description": "Death Probability, by Parity Group -- The probability of death for cows of a single parity group (first lactation, second lactation, etc.); a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total)", - "minimum": 0, - "maximum": 1 - } - }, - "parity_cull_prob": { - "type": "array", - "description": "Cull Probability, by Parity", + "description": "Annual Death Probability, by Parity", "properties": { "type": "number", - "description": "Cull Probability, by Parity Group -- The probability of culling for cows of a single parity group (first lactation, second lactation, etc.); a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total). Culling refers to removal from the herd while alive.", + "description": "Annual Death Probability, by Parity Group -- The probability that a cow of a single parity group (first lactation, second lactation, etc.) dies during a given year; a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total).", "minimum": 0, "maximum": 1 } }, - "death_day_prob": { + "parity_acute_sale_prob": { "type": "array", - "description": "Cumulative distribution function (CDF) values associated with the likelihood of death over time. The 'cull_day_count' array defines the segments of the CDF.", + "description": "Annual Acute-Sale Probability, by Parity", "properties": { "type": "number", + "description": "Annual Acute-Sale Probability, by Parity Group -- The probability that a cow of a single parity group is sold for an acute (forced / involuntary) reason during a given year, regardless of whether a replacement is available; a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total).", "minimum": 0, "maximum": 1 } diff --git a/input/data/animal/example_freestall_animal.json b/input/data/animal/example_freestall_animal.json index b62c3d0914..bc56673da5 100644 --- a/input/data/animal/example_freestall_animal.json +++ b/input/data/animal/example_freestall_animal.json @@ -112,34 +112,8 @@ "std_estrus_cycle_after_pgf": 2 }, "culling": { - "cull_day_count": [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530], - "feet_leg_cull": { - "probability": 0.1633, - "cull_day_prob": [0, 0.03, 0.08, 0.16, 0.25, 0.36, 0.48, 0.59, 0.69, 0.78, 0.85, 0.90, 0.95, 1] - }, - "injury_cull": { - "probability": 0.2883, - "cull_day_prob": [0, 0.08, 0.18, 0.28, 0.38, 0.47, 0.56, 0.64, 0.71, 0.78, 0.85, 0.90, 0.95, 1] - }, - "mastitis_cull": { - "probability": 0.2439, - "cull_day_prob": [0, 0.06, 0.12, 0.19, 0.30, 0.43, 0.56, 0.68, 0.78, 0.85, 0.90, 0.94, 0.97, 1] - }, - "disease_cull": { - "probability": 0.1391, - "cull_day_prob": [0, 0.04, 0.12, 0.24, 0.34, 0.42, 0.50, 0.57, 0.64, 0.72, 0.81, 0.89, 0.95, 1] - }, - "udder_cull": { - "probability": 0.0645, - "cull_day_prob": [0, 0.12, 0.24, 0.33, 0.41, 0.48, 0.55, 0.62, 0.68, 0.76, 0.82, 0.89, 0.95, 1] - }, - "unknown_cull": { - "probability": 0.1009, - "cull_day_prob": [0, 0.05, 0.11, 0.18, 0.27, 0.37, 0.45, 0.54, 0.62, 0.70, 0.77, 0.84, 0.92, 1] - }, "parity_death_prob": [0.039,0.056,0.085,0.117], - "parity_cull_prob": [0.169, 0.233, 0.301, 0.408], - "death_day_prob": [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1] + "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408] }, "life_cycle": { "still_birth_rate": 0.065 diff --git a/input/data/animal/example_open_lot_animal.json b/input/data/animal/example_open_lot_animal.json index 66c0aaf167..7d515431b5 100644 --- a/input/data/animal/example_open_lot_animal.json +++ b/input/data/animal/example_open_lot_animal.json @@ -112,34 +112,8 @@ "std_estrus_cycle_after_pgf": 2 }, "culling": { - "cull_day_count": [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530], - "feet_leg_cull": { - "probability": 0.1633, - "cull_day_prob": [0, 0.03, 0.08, 0.16, 0.25, 0.36, 0.48, 0.59, 0.69, 0.78, 0.85, 0.90, 0.95, 1] - }, - "injury_cull": { - "probability": 0.2883, - "cull_day_prob": [0, 0.08, 0.18, 0.28, 0.38, 0.47, 0.56, 0.64, 0.71, 0.78, 0.85, 0.90, 0.95, 1] - }, - "mastitis_cull": { - "probability": 0.2439, - "cull_day_prob": [0, 0.06, 0.12, 0.19, 0.30, 0.43, 0.56, 0.68, 0.78, 0.85, 0.90, 0.94, 0.97, 1] - }, - "disease_cull": { - "probability": 0.1391, - "cull_day_prob": [0, 0.04, 0.12, 0.24, 0.34, 0.42, 0.50, 0.57, 0.64, 0.72, 0.81, 0.89, 0.95, 1] - }, - "udder_cull": { - "probability": 0.0645, - "cull_day_prob": [0, 0.12, 0.24, 0.33, 0.41, 0.48, 0.55, 0.62, 0.68, 0.76, 0.82, 0.89, 0.95, 1] - }, - "unknown_cull": { - "probability": 0.1009, - "cull_day_prob": [0, 0.05, 0.11, 0.18, 0.27, 0.37, 0.45, 0.54, 0.62, 0.70, 0.77, 0.84, 0.92, 1] - }, "parity_death_prob": [0.039,0.056,0.085,0.117], - "parity_cull_prob": [0.169, 0.233, 0.301, 0.408], - "death_day_prob": [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1] + "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408] }, "life_cycle": { "still_birth_rate": 0.065 diff --git a/input/metadata/cross_validation/animal_cross_validation.json b/input/metadata/cross_validation/animal_cross_validation.json index d592fcdb1b..66cbb7e2db 100644 --- a/input/metadata/cross_validation/animal_cross_validation.json +++ b/input/metadata/cross_validation/animal_cross_validation.json @@ -1053,182 +1053,6 @@ } ] }, - { - "description": "Sum of all cull reason probabilities must equal 1.0", - "aliases": { - "variables": { - "feet_leg_prob": "animal.animal_config.from_literature.culling.feet_leg_cull.probability", - "injury_prob": "animal.animal_config.from_literature.culling.injury_cull.probability", - "mastitis_prob": "animal.animal_config.from_literature.culling.mastitis_cull.probability", - "disease_prob": "animal.animal_config.from_literature.culling.disease_cull.probability", - "udder_prob": "animal.animal_config.from_literature.culling.udder_cull.probability", - "unknown_prob": "animal.animal_config.from_literature.culling.unknown_cull.probability" - }, - "constants": { - "one": 1.0 - } - }, - "rules": [ - { - "left_hand": { - "aggregation": { - "operation": "sum", - "operands": [ - "feet_leg_prob", - "injury_prob", - "mastitis_prob", - "disease_prob", - "udder_prob", - "unknown_prob" - ] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "operands": ["one"] - } - }, - "relationship": "equal" - } - ] - }, - { - "description": "All cull_day_prob and death_day_prob must have the same length as cull_day_count", - "aliases": { - "variables": { - "cull_day_count": "animal.animal_config.from_literature.culling.cull_day_count", - "feet_leg_cull_day_prob": "animal.animal_config.from_literature.culling.feet_leg_cull.cull_day_prob", - "injury_cull_day_prob": "animal.animal_config.from_literature.culling.injury_cull.cull_day_prob", - "mastitis_cull_day_prob": "animal.animal_config.from_literature.culling.mastitis_cull.cull_day_prob", - "disease_cull_day_prob": "animal.animal_config.from_literature.culling.disease_cull.cull_day_prob", - "udder_cull_day_prob": "animal.animal_config.from_literature.culling.udder_cull.cull_day_prob", - "unknown_cull_day_prob": "animal.animal_config.from_literature.culling.unknown_cull.cull_day_prob", - "death_day_prob": "animal.animal_config.from_literature.culling.death_day_prob" - } - }, - "rules": [ - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["feet_leg_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["injury_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["mastitis_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["disease_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["udder_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["unknown_cull_day_prob"] - } - }, - "relationship": "is_equal_length" - }, - { - "left_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["cull_day_count"] - } - }, - "right_hand": { - "aggregation": { - "operation": "no_op", - "mode": "element_wise", - "operands": ["death_day_prob"] - } - }, - "relationship": "is_equal_length" - } - ] - }, { "description": "The maximum days carried calf cannot be more than gestation length", "aliases": { diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal.py b/tests/test_biophysical/test_animal/test_animal/test_animal.py index 0d2da131b0..bc3849533f 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal.py @@ -2384,7 +2384,8 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix ), ) mocker.patch.object(AnimalType, "is_cow", new_callable=PropertyMock, return_value=True) - mocker.patch.object(Animal, "calves", new_callable=PropertyMock, return_value=100) + # First calving (parity 1) triggers the initial annual removal-risk assessment. + mocker.patch.object(Animal, "calves", new_callable=PropertyMock, return_value=1) mocker.patch.object(Animal, "calving_interval_history", new_callable=PropertyMock, return_value=[100]) mocker.patch.object(AnimalEvents, "get_most_recent_date", return_value=2) result, _ = animal.daily_reproduction_update(MagicMock(RufasTime)) @@ -3049,12 +3050,15 @@ def test_determine_future_death_date_no_death(mock_lactating_cow: Animal, mocker def test_determine_future_death_date_with_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """When the annual death roll selects the cow, the day is delegated to _sample_removal_date.""" animal = mock_lactating_cow animal.calves = 5 - animal.days_born = 12 + # random() <= parity_death_probability[3] (0.117) selects the cow for death this year. mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0005) + mock_sample = mocker.patch.object(animal, "_sample_removal_date", return_value=42) result = animal.determine_future_death_date() - assert result == 12 + assert result == 42 + mock_sample.assert_called_once_with(animal_constants.DEATH_TIMING_DAY_PROBABILITY) def test_setup_calf_mortality_disabled_when_rate_zero(mock_calf: Animal, mocker: MockerFixture) -> None: @@ -3217,73 +3221,49 @@ def test_setup_heifer_mortality_not_committed_when_day_already_passed( assert animal._future_death_date is None -def patch_random_first_call(mocker: MockerFixture, first_value: float, second_value: float) -> None: - called = False - - def side_effect() -> float: - nonlocal called - if not called: - called = True - return first_value - return second_value - - mocker.patch("RUFAS.biophysical.animal.animal.random", side_effect=side_effect) - - -def test_determine_future_cull_date_feet_leg(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.05) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (159, animal_constants.LAMENESS_CULL) - - -def test_determine_future_cull_date_injury(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 6 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.25) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (186, animal_constants.INJURY_CULL) - - -def test_determine_future_cull_date_mastitis(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.46) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (295, animal_constants.MASTITIS_CULL) - - -def test_determine_future_cull_date_disease(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.7) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (465, animal_constants.DISEASE_CULL) +def test_determine_future_cull_date_with_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """When the annual acute-sale roll selects the cow, the day is delegated to _sample_removal_date.""" + animal = mock_lactating_cow + animal.calves = 1 + # random() <= parity_acute_sale_probability[0] (0.169) selects the cow for an acute sale. + mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.05) + mock_sample = mocker.patch.object(animal, "_sample_removal_date", return_value=159) + result = animal.determine_future_cull_date() + assert result == (159, animal_constants.ACUTE_SALE_CULL) + mock_sample.assert_called_once_with(animal_constants.ACUTE_SALE_TIMING_DAY_PROBABILITY) -def test_determine_future_cull_date_udder(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: +def test_determine_future_cull_date_no_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.85) + mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) result = mock_lactating_cow.determine_future_cull_date() - assert result == (551, animal_constants.UDDER_CULL) + assert result == (sys.maxsize, "") -def test_determine_future_cull_date_unknown(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - patch_random_first_call(mocker, 0.05, 0.9) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (618, animal_constants.UNKNOWN_CULL) +def test_sample_removal_date_conditions_on_current_days_in_milk( + mock_lactating_cow: Animal, mocker: MockerFixture +) -> None: + """The event is placed at the sampled day in milk, anchored to the current lactation start.""" + animal = mock_lactating_cow + animal.days_born = 150 + animal.days_in_milk = 10 + # A cow at 10 DIM has already passed CDF value 0.25 on the death curve; uniform draws the + # remaining mass. Force the draw to the very start of that remaining mass (0.25) so the + # inverted day in milk is exactly the current DIM (10), landing the event at lactation_start + 10. + mocker.patch("RUFAS.biophysical.animal.animal.uniform", return_value=0.25) + result = animal._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) + # lactation_start (150 - 10) + sampled day in milk (10) == 150. + assert result == 150 -def test_determine_future_cull_date_no_cull(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mock_lactating_cow.days_born = 150 - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (sys.maxsize, "") +def test_sample_removal_date_falls_back_past_last_breakpoint(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A cow past the timing curve's last breakpoint gets a short fixed fallback window.""" + animal = mock_lactating_cow + animal.days_born = 900 + animal.days_in_milk = animal_constants.REMOVAL_TIMING_DAY_BREAKPOINTS[-1] + 5 + mocker.patch("RUFAS.biophysical.animal.animal.randint", return_value=30) + result = animal._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) + assert result == 930 def test_set_nutrient_standard() -> None: diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal_config.py b/tests/test_biophysical/test_animal/test_animal/test_animal_config.py index 633ac25e0d..9e0a702e01 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal_config.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal_config.py @@ -130,34 +130,8 @@ def _make_base_animal_config(repro_sub_protocol: str, heifer_repro_method: str) "std_estrus_cycle_after_pgf": 2, }, "culling": { - "cull_day_count": [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530], - "feet_leg_cull": { - "probability": 0.1633, - "cull_day_prob": [0, 0.03, 0.08, 0.16, 0.25, 0.36, 0.48, 0.59, 0.69, 0.78, 0.85, 0.90, 0.95, 1], - }, - "injury_cull": { - "probability": 0.2883, - "cull_day_prob": [0, 0.08, 0.18, 0.28, 0.38, 0.47, 0.56, 0.64, 0.71, 0.78, 0.85, 0.90, 0.95, 1], - }, - "mastitis_cull": { - "probability": 0.2439, - "cull_day_prob": [0, 0.06, 0.12, 0.19, 0.30, 0.43, 0.56, 0.68, 0.78, 0.85, 0.90, 0.94, 0.97, 1], - }, - "disease_cull": { - "probability": 0.1391, - "cull_day_prob": [0, 0.04, 0.12, 0.24, 0.34, 0.42, 0.50, 0.57, 0.64, 0.72, 0.81, 0.89, 0.95, 1], - }, - "udder_cull": { - "probability": 0.0645, - "cull_day_prob": [0, 0.12, 0.24, 0.33, 0.41, 0.48, 0.55, 0.62, 0.68, 0.76, 0.82, 0.89, 0.95, 1], - }, - "unknown_cull": { - "probability": 0.1009, - "cull_day_prob": [0, 0.05, 0.11, 0.18, 0.27, 0.37, 0.45, 0.54, 0.62, 0.70, 0.77, 0.84, 0.92, 1], - }, "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_cull_prob": [0.169, 0.233, 0.301, 0.408], - "death_day_prob": [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1], + "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], }, "life_cycle": {"still_birth_rate": 0.065}, }, @@ -411,34 +385,8 @@ def test_initialize_animal_config_adds_warning_when_third_check_after_or_on_dryo "std_estrus_cycle_after_pgf": 2, }, "culling": { - "cull_day_count": [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530], - "feet_leg_cull": { - "probability": 0.1633, - "cull_day_prob": [0, 0.03, 0.08, 0.16, 0.25, 0.36, 0.48, 0.59, 0.69, 0.78, 0.85, 0.90, 0.95, 1], - }, - "injury_cull": { - "probability": 0.2883, - "cull_day_prob": [0, 0.08, 0.18, 0.28, 0.38, 0.47, 0.56, 0.64, 0.71, 0.78, 0.85, 0.90, 0.95, 1], - }, - "mastitis_cull": { - "probability": 0.2439, - "cull_day_prob": [0, 0.06, 0.12, 0.19, 0.30, 0.43, 0.56, 0.68, 0.78, 0.85, 0.90, 0.94, 0.97, 1], - }, - "disease_cull": { - "probability": 0.1391, - "cull_day_prob": [0, 0.04, 0.12, 0.24, 0.34, 0.42, 0.50, 0.57, 0.64, 0.72, 0.81, 0.89, 0.95, 1], - }, - "udder_cull": { - "probability": 0.0645, - "cull_day_prob": [0, 0.12, 0.24, 0.33, 0.41, 0.48, 0.55, 0.62, 0.68, 0.76, 0.82, 0.89, 0.95, 1], - }, - "unknown_cull": { - "probability": 0.1009, - "cull_day_prob": [0, 0.05, 0.11, 0.18, 0.27, 0.37, 0.45, 0.54, 0.62, 0.70, 0.77, 0.84, 0.92, 1], - }, "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_cull_prob": [0.169, 0.233, 0.301, 0.408], - "death_day_prob": [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1], + "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], }, "life_cycle": {"still_birth_rate": 0.065}, }, diff --git a/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py b/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py index 8f6630b041..675da98c3c 100644 --- a/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py +++ b/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py @@ -1106,7 +1106,7 @@ def test_report_sold_animal_information(mocker: MockerFixture) -> None: animal_type="LacCow", sold_at_day=123, body_weight=456.78, - cull_reason=animal_constants.UDDER_CULL, + cull_reason=animal_constants.ACUTE_SALE_CULL, days_in_milk=18, parity=2, genetic_history="", @@ -1126,7 +1126,7 @@ def test_report_sold_animal_information(mocker: MockerFixture) -> None: animal_type="DryCow", sold_at_day=123, body_weight=456.78, - cull_reason=animal_constants.LAMENESS_CULL, + cull_reason=animal_constants.ACUTE_SALE_CULL, days_in_milk=0, parity=3, genetic_history="", diff --git a/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py b/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py index 78f71ba185..7c68a568b6 100644 --- a/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py +++ b/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py @@ -235,7 +235,7 @@ def test_reset_cull_reason_stats(herd_statistics: HerdStatistics) -> None: """Test that reset_cull_reason_stats resets cull reason-based attributes correctly.""" # Set non-zero values herd_statistics.cull_reason_stats[animal_constants.DEATH_CULL] = 3 - herd_statistics.cull_reason_stats_percent[animal_constants.LAMENESS_CULL] = 40.5 + herd_statistics.cull_reason_stats_percent[animal_constants.ACUTE_SALE_CULL] = 40.5 herd_statistics.reset_cull_reason_stats() diff --git a/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py b/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py index 3604985a44..30c6f03f32 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py @@ -141,34 +141,8 @@ def animal_json() -> dict[str, Any]: "std_estrus_cycle_after_pgf": 2, }, "culling": { - "cull_day_count": [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530], - "feet_leg_cull": { - "probability": 0.1633, - "cull_day_prob": [0, 0.03, 0.08, 0.16, 0.25, 0.36, 0.48, 0.59, 0.69, 0.78, 0.85, 0.90, 0.95, 1], - }, - "injury_cull": { - "probability": 0.2883, - "cull_day_prob": [0, 0.08, 0.18, 0.28, 0.38, 0.47, 0.56, 0.64, 0.71, 0.78, 0.85, 0.90, 0.95, 1], - }, - "mastitis_cull": { - "probability": 0.2439, - "cull_day_prob": [0, 0.06, 0.12, 0.19, 0.30, 0.43, 0.56, 0.68, 0.78, 0.85, 0.90, 0.94, 0.97, 1], - }, - "disease_cull": { - "probability": 0.1391, - "cull_day_prob": [0, 0.04, 0.12, 0.24, 0.34, 0.42, 0.50, 0.57, 0.64, 0.72, 0.81, 0.89, 0.95, 1], - }, - "udder_cull": { - "probability": 0.0645, - "cull_day_prob": [0, 0.12, 0.24, 0.33, 0.41, 0.48, 0.55, 0.62, 0.68, 0.76, 0.82, 0.89, 0.95, 1], - }, - "unknown_cull": { - "probability": 0.1009, - "cull_day_prob": [0, 0.05, 0.11, 0.18, 0.27, 0.37, 0.45, 0.54, 0.62, 0.70, 0.77, 0.84, 0.92, 1], - }, "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_cull_prob": [0.169, 0.233, 0.301, 0.408], - "death_day_prob": [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1], + "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], }, "life_cycle": {"still_birth_rate": 0.065}, }, diff --git a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py index 830a4b1768..9ea03e6d39 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py @@ -176,23 +176,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st { animal_constants.DEATH_CULL: 0, animal_constants.OVERSUPPLY_CULL: 0, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 0, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 0, }, 0, { animal_constants.DEATH_CULL: 0.0, animal_constants.OVERSUPPLY_CULL: 0.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 0.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 0.0, }, ), # 2. One reason has all culls, matches exit_num -> 100% that reason @@ -200,23 +190,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st { animal_constants.DEATH_CULL: 5, animal_constants.OVERSUPPLY_CULL: 0, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 0, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 0, }, 5, { animal_constants.DEATH_CULL: 100.0, animal_constants.OVERSUPPLY_CULL: 0.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 0.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 0.0, }, ), # 3. Multiple reasons evenly split @@ -225,23 +205,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st { animal_constants.DEATH_CULL: 5, animal_constants.OVERSUPPLY_CULL: 5, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 0, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 0, }, 10, { animal_constants.DEATH_CULL: 50.0, animal_constants.OVERSUPPLY_CULL: 50.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 0.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 0.0, }, ), # 4. Partial distribution @@ -250,23 +220,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st { animal_constants.DEATH_CULL: 3, animal_constants.OVERSUPPLY_CULL: 2, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 0, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 0, }, 10, { animal_constants.DEATH_CULL: 30.0, animal_constants.OVERSUPPLY_CULL: 20.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 0.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 0.0, }, ), # 5. Non-zero exit, some reasons zero @@ -276,23 +236,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st { animal_constants.DEATH_CULL: 2, animal_constants.OVERSUPPLY_CULL: 0, - animal_constants.LAMENESS_CULL: 0, - animal_constants.INJURY_CULL: 0, - animal_constants.MASTITIS_CULL: 0, - animal_constants.DISEASE_CULL: 8, - animal_constants.UDDER_CULL: 0, - animal_constants.UNKNOWN_CULL: 0, + animal_constants.ACUTE_SALE_CULL: 8, }, 10, { animal_constants.DEATH_CULL: 20.0, animal_constants.OVERSUPPLY_CULL: 0.0, - animal_constants.LAMENESS_CULL: 0.0, - animal_constants.INJURY_CULL: 0.0, - animal_constants.MASTITIS_CULL: 0.0, - animal_constants.DISEASE_CULL: 80.0, - animal_constants.UDDER_CULL: 0.0, - animal_constants.UNKNOWN_CULL: 0.0, + animal_constants.ACUTE_SALE_CULL: 80.0, }, ), ], @@ -540,12 +490,7 @@ def test_update_sold_and_died_cow_statistics( """Unit test for _update_sold_and_died_cow_statistics()""" cull_reasons = [ animal_constants.OVERSUPPLY_CULL, - animal_constants.LAMENESS_CULL, - animal_constants.INJURY_CULL, - animal_constants.MASTITIS_CULL, - animal_constants.DISEASE_CULL, - animal_constants.UDDER_CULL, - animal_constants.UNKNOWN_CULL, + animal_constants.ACUTE_SALE_CULL, ] num_sold_cows, num_dead_cows = randint(1, 100), randint(1, 100) @@ -612,12 +557,7 @@ def test_update_sold_and_died_cow_statistics( current_cull_reason_stats = { animal_constants.DEATH_CULL: randint(0, num_total_sold_and_died_cows), animal_constants.OVERSUPPLY_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.LAMENESS_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.INJURY_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.MASTITIS_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.DISEASE_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.UDDER_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.UNKNOWN_CULL: randint(0, num_total_sold_and_died_cows), + animal_constants.ACUTE_SALE_CULL: randint(0, num_total_sold_and_died_cows), } herd_manager.herd_statistics.cull_reason_stats = current_cull_reason_stats expected_cull_reason_stats = { @@ -625,18 +565,8 @@ def test_update_sold_and_died_cow_statistics( + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.DEATH_CULL]), animal_constants.OVERSUPPLY_CULL: current_cull_reason_stats[animal_constants.OVERSUPPLY_CULL] + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.OVERSUPPLY_CULL]), - animal_constants.LAMENESS_CULL: current_cull_reason_stats[animal_constants.LAMENESS_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.LAMENESS_CULL]), - animal_constants.INJURY_CULL: current_cull_reason_stats[animal_constants.INJURY_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.INJURY_CULL]), - animal_constants.MASTITIS_CULL: current_cull_reason_stats[animal_constants.MASTITIS_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.MASTITIS_CULL]), - animal_constants.DISEASE_CULL: current_cull_reason_stats[animal_constants.DISEASE_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.DISEASE_CULL]), - animal_constants.UDDER_CULL: current_cull_reason_stats[animal_constants.UDDER_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.UDDER_CULL]), - animal_constants.UNKNOWN_CULL: current_cull_reason_stats[animal_constants.UNKNOWN_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.UNKNOWN_CULL]), + animal_constants.ACUTE_SALE_CULL: current_cull_reason_stats[animal_constants.ACUTE_SALE_CULL] + + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.ACUTE_SALE_CULL]), } current_sold_cow_num = randint(0, current_cow_herd_exit_num) From dcee6fe1e79106a880c9b4afa0c9114760b29d1e Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 27 Aug 2026 18:56:57 +0000 Subject: [PATCH 02/13] Apply Black Formatting From ab2f54a740a08a8ebf3707e3183e07d696a1ae15 Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Mon, 7 Sep 2026 23:02:45 -0400 Subject: [PATCH 03/13] implement feedback --- RUFAS/biophysical/animal/animal.py | 155 ++++-------------- RUFAS/biophysical/animal/animal_constants.py | 45 +---- .../test_animal/test_animal/test_animal.py | 74 ++------- 3 files changed, 48 insertions(+), 226 deletions(-) diff --git a/RUFAS/biophysical/animal/animal.py b/RUFAS/biophysical/animal/animal.py index 955ca7a1d1..dbe31b6923 100644 --- a/RUFAS/biophysical/animal/animal.py +++ b/RUFAS/biophysical/animal/animal.py @@ -1,6 +1,6 @@ import sys from datetime import timedelta -from random import random, randint, uniform +from random import random, randint from typing import Callable, cast from scipy.stats import truncnorm @@ -1767,12 +1767,8 @@ def daily_reproduction_update( self.milk_production.set_wood_parameters( wood_parameters["l"], wood_parameters["m"], wood_parameters["n"] ) - # A heifer becomes a cow at her first calving; give her an initial mortality / - # acute-sale assessment now so she is not risk-free until the next Jan 1. Later - # lactations are not re-rolled here -- removal risk is reassessed annually - # (see ``_assess_annual_removal_risk``), decoupling it from reproduction. if self.calves == 1: - self._assess_annual_removal_risk() + self._assess_removal_risk() self.events += reproduction_outputs.events @@ -1810,12 +1806,8 @@ def daily_routines(self, time: RufasTime) -> DailyRoutinesOutput: newborn_calf_config, daily_routines_output.herd_reproduction_statistics = self.daily_reproduction_update(time) - # Reassess mortality / acute-sale risk once a year (Jan 1, Julian day 1) for every cow, so - # risk accrues with time in the herd rather than only at calving (issue #2694). The guard - # inside the assessment skips risk types that already have a pending event, so a cow that - # first calved earlier today is not double-rolled. - if self.animal_type.is_cow and time.current_julian_day == 1: - self._assess_annual_removal_risk() + if self.animal_type.is_cow: + self._assess_removal_risk() daily_routines_output.animal_status, daily_routines_output.newborn_calf_config = self.animal_life_stage_update( time @@ -2432,15 +2424,15 @@ def _get_cow_values(self) -> CowValuesTypedDict: parity=self.calves, ) - def _assess_annual_removal_risk(self) -> None: + def _assess_removal_risk(self) -> None: """ - Roll a cow's annual mortality and acute-sale risk and schedule any resulting removal. + Roll a cow's daily mortality and acute-sale risk and schedule any resulting removal. - Called on each Jan 1 for every cow, and once when a heifer first calves, so removal risk + Called each day for every cow, and once when a heifer first calves, so removal risk accrues with time spent in the herd rather than only at calving (issue #2694). Death and acute sale are rolled independently. A risk type that already has a pending future event is - left untouched: this prevents overwriting or double-scheduling an event a prior roll placed - (the days-in-milk timing can land a scheduled event more than a year out, across a Jan 1). + left untouched, preventing a later roll from overwriting a removal a prior roll scheduled. + A cow selected for removal is scheduled to leave the herd immediately (on the current day). Notes ------- @@ -2448,145 +2440,58 @@ def _assess_annual_removal_risk(self) -> None: """ if self.future_death_date == sys.maxsize: - death_date = self.determine_future_death_date() - if death_date != sys.maxsize: - self.future_death_date = death_date + if self.will_die_tomorrow(): + self.future_death_date = self.days_born self._future_death_reason = animal_constants.DEATH_CULL if self.future_cull_date == sys.maxsize: - cull_date, cull_reason = self.determine_future_cull_date() - if cull_date != sys.maxsize: - self.future_cull_date = cull_date - self.cull_reason = cull_reason + if self.will_be_sold_tomorrow(): + self.future_cull_date = self.days_born + self.cull_reason = animal_constants.ACUTE_SALE_CULL def _parity_index(self) -> int: """Return the 0-based index into a by-parity array, capping parity 4+ at the last entry.""" return 3 if self.calves >= 4 else self.calves - 1 - @staticmethod - def _interpolate_cdf_at_day(day: float, cdf: list[float], breakpoints: list[int]) -> float: - """ - Linearly interpolate a cumulative-distribution value at a given day in milk. - - Parameters - ---------- - day : float - Day in milk at which to evaluate the CDF. - cdf : list[float] - Cumulative-distribution values at each breakpoint (non-decreasing, 0.0 to 1.0). - breakpoints : list[int] - Days-in-milk breakpoints that partition the CDF (same length as ``cdf``). - - Returns - ------- - float - The interpolated CDF value, clamped to ``cdf[0]`` below the first breakpoint and - ``cdf[-1]`` at or beyond the last. - - """ - if day <= breakpoints[0]: - return cdf[0] - for i in range(len(breakpoints) - 1): - if breakpoints[i] <= day < breakpoints[i + 1]: - slope = (cdf[i + 1] - cdf[i]) / (breakpoints[i + 1] - breakpoints[i]) - return cdf[i] + slope * (day - breakpoints[i]) - return cdf[-1] - - def _sample_removal_date(self, timing_cdf: list[float]) -> int: - """ - Sample an absolute simulation day for a scheduled removal from a days-in-milk timing CDF. - - The day is shaped by ``timing_cdf`` (the distribution of removal timing across a lactation) - but conditioned to fall after the cow's current days in milk, so a cow selected for removal - is always scheduled to leave in the future rather than "escaping" the event. If the cow is - already past the CDF's last breakpoint (an extended lactation), the event is instead placed - uniformly within ``REMOVAL_FALLBACK_WINDOW_DAYS`` of today. - - Parameters - ---------- - timing_cdf : list[float] - Cumulative-distribution values of removal timing at each of - ``REMOVAL_TIMING_DAY_BREAKPOINTS``. - - Returns - ------- - int - The absolute simulation day (in ``days_born`` terms) on which the removal occurs. - - Notes - ------- - [AN.ANM.1], [AN.ANM.2] - - """ - breakpoints = animal_constants.REMOVAL_TIMING_DAY_BREAKPOINTS - current_days_in_milk = self.days_in_milk - lactation_start = self.days_born - current_days_in_milk - - if current_days_in_milk >= breakpoints[-1]: - return self.days_born + randint(1, animal_constants.REMOVAL_FALLBACK_WINDOW_DAYS) - - # Draw uniformly on the CDF mass that remains after the current day in milk, then invert - # back to a day in milk so the timing keeps its lactation-stage shape. - lower_cdf_value = self._interpolate_cdf_at_day(current_days_in_milk, timing_cdf, breakpoints) - removal_cdf_value = uniform(lower_cdf_value, 1.0) - for i in range(len(timing_cdf) - 1): - if timing_cdf[i] <= removal_cdf_value < timing_cdf[i + 1]: - slope = (breakpoints[i + 1] - breakpoints[i]) / (timing_cdf[i + 1] - timing_cdf[i]) - day_in_milk = breakpoints[i] + slope * (removal_cdf_value - timing_cdf[i]) - return round(lactation_start + day_in_milk) - # ``uniform`` can return its upper bound (1.0), which no half-open segment contains. - return round(lactation_start + breakpoints[-1]) - - def determine_future_death_date(self) -> int: + def will_die_tomorrow(self) -> bool: """ - Roll the cow's annual mortality risk and, if selected, schedule a death day. + Roll the cow's daily mortality risk. - The parity-indexed :attr:`AnimalConfig.parity_death_probability` is now an *annual* - probability (issue #2694), rolled by the annual removal assessment rather than once per - lactation. When the cow is selected to die, the day is drawn from the days-in-milk death - timing CDF via :meth:`_sample_removal_date`. + The parity-indexed annual :attr:`AnimalConfig.parity_death_probability` is converted to a + daily rate by dividing by 365 and rolled once per day (issue #2694). Returns ------- - int - Calculated future death date in simulation days, or ``sys.maxsize`` if not selected. + bool + ``True`` if the cow is selected to die, ``False`` otherwise. Notes ------- [AN.ANM.1] """ - death_rate = AnimalConfig.parity_death_probability[self._parity_index()] - if random() <= death_rate: - return self._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) - return sys.maxsize + death_rate = AnimalConfig.parity_death_probability[self._parity_index()] / 365 + return random() <= death_rate - def determine_future_cull_date(self) -> tuple[int, str]: + def will_be_sold_tomorrow(self) -> bool: """ - Roll the cow's annual acute-sale risk and, if selected, schedule an acute-sale day. + Roll the cow's daily acute-sale (forced / involuntary) risk. - An acute sale ("forced" / "involuntary" / "spontaneous" removal) is a cow that must leave - the herd immediately regardless of whether a replacement is available. The parity-indexed - :attr:`AnimalConfig.parity_acute_sale_probability` is now an *annual* probability - (issue #2694); the former six disease-specific reasons are collapsed into the single - :data:`animal_constants.ACUTE_SALE_CULL` reason with one timing CDF. + The parity-indexed annual :attr:`AnimalConfig.parity_acute_sale_probability` is converted to + a daily rate by dividing by 365 and rolled once per day (issue #2694). Returns ------- - tuple[int, str] - - Future acute-sale date in simulation days (``sys.maxsize`` if not selected). - - Reason for removal (empty string if not selected). + bool + ``True`` if the cow is selected for an acute sale, ``False`` otherwise. Notes ------- [AN.ANM.2] """ - acute_sale_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] - if random() <= acute_sale_rate: - future_cull_date = self._sample_removal_date(animal_constants.ACUTE_SALE_TIMING_DAY_PROBABILITY) - return future_cull_date, animal_constants.ACUTE_SALE_CULL - return sys.maxsize, "" + acute_sale_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] / 365 + return random() <= acute_sale_rate def update_pen_history(self, current_pen: int, current_day: int, animal_types_in_pen: set[AnimalType]) -> None: """ diff --git a/RUFAS/biophysical/animal/animal_constants.py b/RUFAS/biophysical/animal/animal_constants.py index 6672112c0c..cf905ffc77 100644 --- a/RUFAS/biophysical/animal/animal_constants.py +++ b/RUFAS/biophysical/animal/animal_constants.py @@ -96,55 +96,12 @@ HEIFER_REPRO_CULL = "culled for heifer reproductive problem" OVERSUPPLY_CULL = "culled for herd resize" DEATH_CULL = "culled for death" -# A single "acute sale" (a.k.a. forced / involuntary / spontaneous) removal reason replaces the -# former six disease-specific cull reasons (feet-and-leg, injury, mastitis, disease, udder, -# unknown). These represent animals that must leave the herd immediately regardless of whether a -# replacement is available; the disease-level detail is not recoverable from farm records. -ACUTE_SALE_CULL = "sold for acute (involuntary) reason" +ACUTE_SALE_CULL = "culled for acute sale" # youngstock mortality (a loss from death, not a cull) CALF_MORTALITY_LOSS = "died from pre-wean mortality" HEIFER_MORTALITY_LOSS = "died from post-wean mortality" -# Removal-event timing distributions (moved out of user inputs; see issue #2694). -# -# These describe *when*, within a cow's lactation (days in milk), a scheduled death or acute -# sale occurs. They are model constants -- not recoverable from farm records -- so they live -# here rather than in animal.json. REMOVAL_TIMING_DAY_BREAKPOINTS are the days-in-milk -# breakpoints that partition each cumulative distribution function (CDF); the two probability -# arrays give the CDF value at each breakpoint (they must start at 0.0, end at 1.0, be -# non-decreasing, and match the breakpoints in length). -REMOVAL_TIMING_DAY_BREAKPOINTS = [0, 5, 15, 45, 90, 135, 180, 225, 270, 330, 380, 430, 480, 530] - -# CDF of death timing by days in milk. -DEATH_TIMING_DAY_PROBABILITY = [0, 0.18, 0.32, 0.42, 0.48, 0.54, 0.60, 0.65, 0.70, 0.77, 0.83, 0.89, 0.95, 1] - -# CDF of acute-sale timing by days in milk. This single curve replaces the six former -# reason-specific curves. The placeholder values below are the mean of those six curves; the -# final values should be confirmed via literature review / expert input (issue #2694). -ACUTE_SALE_TIMING_DAY_PROBABILITY = [ - 0, - 0.063, - 0.142, - 0.230, - 0.325, - 0.422, - 0.517, - 0.607, - 0.687, - 0.765, - 0.833, - 0.893, - 0.948, - 1, -] - -# When a cow selected for removal is already past the timing CDF's last breakpoint (e.g. an -# extended lactation with days in milk >= REMOVAL_TIMING_DAY_BREAKPOINTS[-1]), the day-in-milk -# curve has no remaining support, so the event is instead scheduled uniformly within this many -# days from the current day. -REMOVAL_FALLBACK_WINDOW_DAYS = 60 - # STATS STDI = 2 diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal.py b/tests/test_biophysical/test_animal/test_animal/test_animal.py index bc3849533f..f00746b428 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal.py @@ -2356,10 +2356,7 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix animal.animal_type = AnimalType.HEIFER_II mock_determine_days_in_milk = mocker.patch.object(animal, "_determine_days_in_milk", return_value=3) mock_set_wood_parameters = mocker.patch.object(MilkProduction, "set_wood_parameters") - mock_determine_future_death_date = mocker.patch.object(animal, "determine_future_death_date", return_value=3) - mock_determine_future_cull_date = mocker.patch.object( - animal, "determine_future_cull_date", return_value=(3, "test") - ) + mock_assess_removal_risk = mocker.patch.object(animal, "_assess_removal_risk") mock_get_wood_parameters = mocker.patch.object( LactationCurve, "get_wood_parameters", return_value={"l": 10.2, "m": 41.2, "n": 41.8} ) @@ -2384,7 +2381,7 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix ), ) mocker.patch.object(AnimalType, "is_cow", new_callable=PropertyMock, return_value=True) - # First calving (parity 1) triggers the initial annual removal-risk assessment. + # First calving (parity 1) triggers the initial removal-risk assessment. mocker.patch.object(Animal, "calves", new_callable=PropertyMock, return_value=1) mocker.patch.object(Animal, "calving_interval_history", new_callable=PropertyMock, return_value=[100]) mocker.patch.object(AnimalEvents, "get_most_recent_date", return_value=2) @@ -2394,11 +2391,8 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix mock_get_wood_parameters.assert_called_once() mock_set_wood_parameters.assert_called_once() mock_determine_days_in_milk.assert_called_once() - mock_determine_future_cull_date.assert_called_once() - mock_determine_future_death_date.assert_called_once() + mock_assess_removal_risk.assert_called_once() - assert animal.future_cull_date == 3 - assert animal.cull_reason == "test" assert animal.days_in_milk == 3 assert animal.body_weight == 10 assert animal.days_in_pregnancy == 12 @@ -3040,25 +3034,21 @@ def test_get_cow_values(mock_lactating_cow: Animal) -> None: assert mock_lactating_cow._get_cow_values() == expected -def test_determine_future_death_date_no_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: +def test_will_die_tomorrow_no_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: animal = mock_lactating_cow animal.calves = 1 animal.days_born = 150 mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - result = animal.determine_future_death_date() - assert result == sys.maxsize + assert animal.will_die_tomorrow() is False -def test_determine_future_death_date_with_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - """When the annual death roll selects the cow, the day is delegated to _sample_removal_date.""" +def test_will_die_tomorrow_with_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A roll below the daily death rate (annual parity probability / 365) selects the cow.""" animal = mock_lactating_cow animal.calves = 5 - # random() <= parity_death_probability[3] (0.117) selects the cow for death this year. - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0005) - mock_sample = mocker.patch.object(animal, "_sample_removal_date", return_value=42) - result = animal.determine_future_death_date() - assert result == 42 - mock_sample.assert_called_once_with(animal_constants.DEATH_TIMING_DAY_PROBABILITY) + # daily death rate = parity_death_probability[3] (0.117) / 365 ~= 0.00032. + mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) + assert animal.will_die_tomorrow() is True def test_setup_calf_mortality_disabled_when_rate_zero(mock_calf: Animal, mocker: MockerFixture) -> None: @@ -3221,49 +3211,19 @@ def test_setup_heifer_mortality_not_committed_when_day_already_passed( assert animal._future_death_date is None -def test_determine_future_cull_date_with_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - """When the annual acute-sale roll selects the cow, the day is delegated to _sample_removal_date.""" +def test_will_be_sold_tomorrow_with_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A roll below the daily acute-sale rate (annual parity probability / 365) selects the cow.""" animal = mock_lactating_cow animal.calves = 1 - # random() <= parity_acute_sale_probability[0] (0.169) selects the cow for an acute sale. - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.05) - mock_sample = mocker.patch.object(animal, "_sample_removal_date", return_value=159) - result = animal.determine_future_cull_date() - assert result == (159, animal_constants.ACUTE_SALE_CULL) - mock_sample.assert_called_once_with(animal_constants.ACUTE_SALE_TIMING_DAY_PROBABILITY) + # daily acute-sale rate = parity_acute_sale_probability[0] (0.169) / 365 ~= 0.00046. + mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) + assert animal.will_be_sold_tomorrow() is True -def test_determine_future_cull_date_no_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: +def test_will_be_sold_tomorrow_no_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: mock_lactating_cow.calves = 1 mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - result = mock_lactating_cow.determine_future_cull_date() - assert result == (sys.maxsize, "") - - -def test_sample_removal_date_conditions_on_current_days_in_milk( - mock_lactating_cow: Animal, mocker: MockerFixture -) -> None: - """The event is placed at the sampled day in milk, anchored to the current lactation start.""" - animal = mock_lactating_cow - animal.days_born = 150 - animal.days_in_milk = 10 - # A cow at 10 DIM has already passed CDF value 0.25 on the death curve; uniform draws the - # remaining mass. Force the draw to the very start of that remaining mass (0.25) so the - # inverted day in milk is exactly the current DIM (10), landing the event at lactation_start + 10. - mocker.patch("RUFAS.biophysical.animal.animal.uniform", return_value=0.25) - result = animal._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) - # lactation_start (150 - 10) + sampled day in milk (10) == 150. - assert result == 150 - - -def test_sample_removal_date_falls_back_past_last_breakpoint(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - """A cow past the timing curve's last breakpoint gets a short fixed fallback window.""" - animal = mock_lactating_cow - animal.days_born = 900 - animal.days_in_milk = animal_constants.REMOVAL_TIMING_DAY_BREAKPOINTS[-1] + 5 - mocker.patch("RUFAS.biophysical.animal.animal.randint", return_value=30) - result = animal._sample_removal_date(animal_constants.DEATH_TIMING_DAY_PROBABILITY) - assert result == 930 + assert mock_lactating_cow.will_be_sold_tomorrow() is False def test_set_nutrient_standard() -> None: From c259a2221060a8fb9561d9fbb3bac7f12fdf1835 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 8 Sep 2026 03:04:39 +0000 Subject: [PATCH 04/13] Apply Black Formatting From 2f4535f1d0a7b7caada77d6e0e23a53591483f39 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 8 Sep 2026 03:12:36 +0000 Subject: [PATCH 05/13] Apply Black Formatting From 38bc8dc305ffec327971ef317d26f4693d0cdb79 Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Mon, 7 Sep 2026 23:18:18 -0400 Subject: [PATCH 06/13] Update changelog_WIP.md --- changelog_WIP.md | 1 + 1 file changed, 1 insertion(+) diff --git a/changelog_WIP.md b/changelog_WIP.md index dcb8212dfb..6c583fc7b6 100644 --- a/changelog_WIP.md +++ b/changelog_WIP.md @@ -119,3 +119,4 @@ This **WIP Changelog** records development changes in progress and not yet inclu - [3161](https://github.com/RuminantFarmSystems/RuFaS/pull/3161) - [minor change] [Docs] [Governance] [NoInputChange] [NoOutputChange] Adds Strategy.md defining the RuFaS vision, mission, strategic pillars, core values, and operating principles, and aligns README.md, CONTRIBUTING.md, and FORKING.md with the approved strategic direction. - [3223](https://github.com/RuminantFarmSystems/RuFaS/pull/3223) - [minor change] [PostProcessing] [OutputManager] [NoInputChange] [NoOutputChange] Establishes new overhauled version of OutputManager and the subclasses it oversees. - [3235](https://github.com/RuminantFarmSystems/RuFaS/pull/3235) - [minor change] [Dependabot] [NoInputChange] [NoOutputChange] Updates file-target of dependabot-change PRs for tagging dev-team members for review. +- [3241](https://github.com/RuminantFarmSystems/RuFaS/pull/3241) - [minor change] [Animal] [InputChange] [OutputChange] Replaces the per-lactation culling model with annual, parity-indexed death and acute-sale probabilities rolled daily. From f089395a44fa6d3c0f7156572bc333e99deb5c1e Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 8 Sep 2026 03:20:10 +0000 Subject: [PATCH 07/13] Apply Black Formatting From d9b243ee14bf621a30e0d42d499c99887529484e Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Wed, 16 Sep 2026 13:22:23 -0400 Subject: [PATCH 08/13] implement % fresh cow as culling/death assessment factor --- RUFAS/biophysical/animal/animal.py | 68 +++++++++---------- RUFAS/biophysical/animal/animal_config.py | 3 - RUFAS/biophysical/animal/herd_manager.py | 51 +++++++++++++- .../test_animal/test_animal/test_animal.py | 10 ++- .../test_herd_manager_daily_routines.py | 2 +- 5 files changed, 84 insertions(+), 50 deletions(-) diff --git a/RUFAS/biophysical/animal/animal.py b/RUFAS/biophysical/animal/animal.py index dbe31b6923..4075a04760 100644 --- a/RUFAS/biophysical/animal/animal.py +++ b/RUFAS/biophysical/animal/animal.py @@ -1767,8 +1767,6 @@ def daily_reproduction_update( self.milk_production.set_wood_parameters( wood_parameters["l"], wood_parameters["m"], wood_parameters["n"] ) - if self.calves == 1: - self._assess_removal_risk() self.events += reproduction_outputs.events @@ -1806,9 +1804,6 @@ def daily_routines(self, time: RufasTime) -> DailyRoutinesOutput: newborn_calf_config, daily_routines_output.herd_reproduction_statistics = self.daily_reproduction_update(time) - if self.animal_type.is_cow: - self._assess_removal_risk() - daily_routines_output.animal_status, daily_routines_output.newborn_calf_config = self.animal_life_stage_update( time ) @@ -2424,74 +2419,73 @@ def _get_cow_values(self) -> CowValuesTypedDict: parity=self.calves, ) - def _assess_removal_risk(self) -> None: + def assess_removal_risk(self, percent_fresh: float, time: RufasTime) -> None: """ Roll a cow's daily mortality and acute-sale risk and schedule any resulting removal. - Called each day for every cow, and once when a heifer first calves, so removal risk - accrues with time spent in the herd rather than only at calving (issue #2694). Death and - acute sale are rolled independently. A risk type that already has a pending future event is - left untouched, preventing a later roll from overwriting a removal a prior roll scheduled. - A cow selected for removal is scheduled to leave the herd immediately (on the current day). + Parameters + ---------- + percent_fresh : float + Fraction of the herd currently fresh, used to scale the daily death and + acute-sale selection probabilities. + time : RufasTime + Current simulation time, used to record the day of death or sale. - Notes + Returns ------- - [AN.ANM.1], [AN.ANM.2] - + None """ if self.future_death_date == sys.maxsize: - if self.will_die_tomorrow(): + if self.is_selected_for_death(percent_fresh): self.future_death_date = self.days_born - self._future_death_reason = animal_constants.DEATH_CULL + self.cull_reason = self._future_death_reason = animal_constants.DEATH_CULL + self.dead_at_day = time.simulation_day if self.future_cull_date == sys.maxsize: - if self.will_be_sold_tomorrow(): + if self.is_selected_for_acute_sale(percent_fresh): self.future_cull_date = self.days_born self.cull_reason = animal_constants.ACUTE_SALE_CULL + self.sold_at_day = time.simulation_day def _parity_index(self) -> int: """Return the 0-based index into a by-parity array, capping parity 4+ at the last entry.""" return 3 if self.calves >= 4 else self.calves - 1 - def will_die_tomorrow(self) -> bool: + def is_selected_for_death(self, percent_fresh: float) -> bool: """ Roll the cow's daily mortality risk. - The parity-indexed annual :attr:`AnimalConfig.parity_death_probability` is converted to a - daily rate by dividing by 365 and rolled once per day (issue #2694). - Returns ------- bool ``True`` if the cow is selected to die, ``False`` otherwise. + """ + average_daily_death_rate = AnimalConfig.parity_death_probability[self._parity_index()] / 365 + percent_other = 1 - percent_fresh - Notes - ------- - [AN.ANM.1] + daily_death_risk_other = average_daily_death_rate / (2.55 * percent_fresh + percent_other) + daily_death_risk_fresh = 2.55 * daily_death_risk_other - """ - death_rate = AnimalConfig.parity_death_probability[self._parity_index()] / 365 - return random() <= death_rate + daily_death_risk = daily_death_risk_fresh if self.days_in_milk < 50 else daily_death_risk_other + return random() <= daily_death_risk - def will_be_sold_tomorrow(self) -> bool: + def is_selected_for_acute_sale(self, percent_fresh: float) -> bool: """ Roll the cow's daily acute-sale (forced / involuntary) risk. - The parity-indexed annual :attr:`AnimalConfig.parity_acute_sale_probability` is converted to - a daily rate by dividing by 365 and rolled once per day (issue #2694). - Returns ------- bool ``True`` if the cow is selected for an acute sale, ``False`` otherwise. + """ + average_daily_removal_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] / 365 + percent_other = 1 - percent_fresh - Notes - ------- - [AN.ANM.2] + daily_removal_risk_other = average_daily_removal_rate / (2.55 * percent_fresh + percent_other) + daily_removal_risk_fresh = 2.55 * daily_removal_risk_other - """ - acute_sale_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] / 365 - return random() <= acute_sale_rate + daily_removal_risk = daily_removal_risk_fresh if self.days_in_milk < 50 else daily_removal_risk_other + return random() <= daily_removal_risk def update_pen_history(self, current_pen: int, current_day: int, animal_types_in_pen: set[AnimalType]) -> None: """ diff --git a/RUFAS/biophysical/animal/animal_config.py b/RUFAS/biophysical/animal/animal_config.py index fb89aaf780..359f9211b7 100644 --- a/RUFAS/biophysical/animal/animal_config.py +++ b/RUFAS/biophysical/animal/animal_config.py @@ -240,9 +240,6 @@ class AnimalConfig: third_pregnancy_check_day: int = 200 third_pregnancy_check_loss_rate: float = 0.017 - # Annual, parity-indexed probabilities that a cow dies (parity_death_probability) or is sold - # for an acute / involuntary reason (parity_acute_sale_probability) during a given year. See - # issue #2694: these were formerly evaluated once per lactation and are now evaluated annually. parity_death_probability: list[float] = [0.039, 0.056, 0.085, 0.117] parity_acute_sale_probability: list[float] = [0.169, 0.233, 0.301, 0.408] diff --git a/RUFAS/biophysical/animal/herd_manager.py b/RUFAS/biophysical/animal/herd_manager.py index 18c91519f9..699a976920 100644 --- a/RUFAS/biophysical/animal/herd_manager.py +++ b/RUFAS/biophysical/animal/herd_manager.py @@ -569,6 +569,42 @@ def _perform_daily_routines_for_animals( animal.update_genetic_history(simulation_day=time.simulation_day) return (graduated_animals, sold_animals, stillborn_newborn_calves, newborn_calves, sold_newborn_calves) + def _assess_removal_risk(self, animals: list[Animal], time: RufasTime) -> tuple[list[Animal], list[Animal]]: + """ + Assess daily removal risk for each cow and collect those removed. + + Computes the fresh fraction across all cows in the herd, then rolls death and + acute-sale risk for every cow in ``animals`` via + :meth:`Animal.assess_removal_risk`. Non-cow animals are skipped. + + Parameters + ---------- + animals : list of Animal + Animals to assess on the current day. + time : RufasTime + Current simulation time, passed through to the per-animal risk assessment. + + Returns + ------- + tuple[list[Animal], list[Animal]] + A ``(sold_cows, dead_cows)`` pair listing the cows selected for acute sale + and the cows selected to die, respectively. + """ + sold_cows: list[Animal] = [] + dead_cows: list[Animal] = [] + + all_cows = [animal for animal in self.all_animals if animal.animal_type.is_cow] + fresh_cows: list[Animal] = [cow for cow in all_cows if cow.days_in_milk < 50] + percent_fresh_cows = len(fresh_cows) / len(all_cows) if len(all_cows) > 0 else 0 + for animal in animals: + if animal.animal_type.is_cow: + animal.assess_removal_risk(percent_fresh_cows, time) + if animal.sold: + sold_cows.append(animal) + if animal.dead: + dead_cows.append(animal) + return (sold_cows, dead_cows) + def _update_genetic_values_at_lactation_start(self, animal: Animal, time: RufasTime) -> None: """ Updates the genetic values of an animal at the start of a new lactation. @@ -642,17 +678,26 @@ def _process_daily_herd_updates(self, time: RufasTime) -> DailyHerdUpdates: group_sold_newborn_calves, ) = self._perform_daily_routines_for_animals(time, animals) collect_birth_results = animal_group_name in ["heiferIIIs", "cows"] - daily_herd_updates.graduated_animals += group_graduated_animals - daily_herd_updates.removed_animals += sold_animals if collect_birth_results: daily_herd_updates.stillborn_newborn_calves += group_stillborn_newborn_calves daily_herd_updates.newborn_calves += group_newborn_calves daily_herd_updates.sold_newborn_calves += group_sold_newborn_calves if animal_group_name == "heiferIIs": daily_herd_updates.sold_heiferIIs = sold_animals + elif animal_group_name == "heiferIIIs": + sold_cows, dead_cows = self._assess_removal_risk(group_graduated_animals, time) + sold_animals.extend(sold_cows) + sold_animals.extend(dead_cows) + daily_herd_updates.sold_and_died_cows.extend(sold_cows) + daily_herd_updates.sold_and_died_cows.extend(dead_cows) elif animal_group_name == "cows": - daily_herd_updates.sold_and_died_cows = sold_animals + sold_cows, dead_cows = self._assess_removal_risk(animals, time) + sold_animals.extend(sold_cows) + sold_animals.extend(dead_cows) + daily_herd_updates.sold_and_died_cows.extend(sold_animals) + daily_herd_updates.graduated_animals += group_graduated_animals + daily_herd_updates.removed_animals += sold_animals return daily_herd_updates def _apply_daily_herd_structure_updates( diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal.py b/tests/test_biophysical/test_animal/test_animal/test_animal.py index f00746b428..463c69eab6 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal.py @@ -2356,7 +2356,6 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix animal.animal_type = AnimalType.HEIFER_II mock_determine_days_in_milk = mocker.patch.object(animal, "_determine_days_in_milk", return_value=3) mock_set_wood_parameters = mocker.patch.object(MilkProduction, "set_wood_parameters") - mock_assess_removal_risk = mocker.patch.object(animal, "_assess_removal_risk") mock_get_wood_parameters = mocker.patch.object( LactationCurve, "get_wood_parameters", return_value={"l": 10.2, "m": 41.2, "n": 41.8} ) @@ -2391,7 +2390,6 @@ def test_daily_reproduction_update(mock_lactating_cow: Animal, mocker: MockerFix mock_get_wood_parameters.assert_called_once() mock_set_wood_parameters.assert_called_once() mock_determine_days_in_milk.assert_called_once() - mock_assess_removal_risk.assert_called_once() assert animal.days_in_milk == 3 assert animal.body_weight == 10 @@ -3039,7 +3037,7 @@ def test_will_die_tomorrow_no_death(mock_lactating_cow: Animal, mocker: MockerFi animal.calves = 1 animal.days_born = 150 mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - assert animal.will_die_tomorrow() is False + assert animal.is_selected_for_death(percent_fresh=0.15) is False def test_will_die_tomorrow_with_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: @@ -3048,7 +3046,7 @@ def test_will_die_tomorrow_with_death(mock_lactating_cow: Animal, mocker: Mocker animal.calves = 5 # daily death rate = parity_death_probability[3] (0.117) / 365 ~= 0.00032. mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) - assert animal.will_die_tomorrow() is True + assert animal.is_selected_for_death(percent_fresh=0.15) is True def test_setup_calf_mortality_disabled_when_rate_zero(mock_calf: Animal, mocker: MockerFixture) -> None: @@ -3217,13 +3215,13 @@ def test_will_be_sold_tomorrow_with_acute_sale(mock_lactating_cow: Animal, mocke animal.calves = 1 # daily acute-sale rate = parity_acute_sale_probability[0] (0.169) / 365 ~= 0.00046. mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) - assert animal.will_be_sold_tomorrow() is True + assert animal.is_selected_for_acute_sale(percent_fresh=0.15) is True def test_will_be_sold_tomorrow_no_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: mock_lactating_cow.calves = 1 mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - assert mock_lactating_cow.will_be_sold_tomorrow() is False + assert mock_lactating_cow.is_selected_for_acute_sale(percent_fresh=0.15) is False def test_set_nutrient_standard() -> None: diff --git a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py index ad29de3569..e9c95a75a0 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py @@ -748,7 +748,7 @@ def test_daily_routines(herd_manager: HerdManager, mock_herd: dict[str, list[Ani mock_update_sold_animal_statistics.assert_called_once_with( sold_newborn_calves=[], sold_heiferIIs=sold_heiferIIs, - sold_and_died_cows=sold_and_died_cows, + sold_and_died_cows=graduated_heiferIIIs + sold_and_died_cows, ) assert mock_check_if_cows_need_to_be_sold.call_count == 0 assert mock_check_if_replacement_heifers_needed.call_count == 0 From 7eec5cdfb8586332d3d900e4932dd3ce0c62fdfe Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Thu, 17 Sep 2026 11:14:20 -0400 Subject: [PATCH 09/13] Update animal.py --- RUFAS/biophysical/animal/animal.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/RUFAS/biophysical/animal/animal.py b/RUFAS/biophysical/animal/animal.py index 4075a04760..a56f6c5441 100644 --- a/RUFAS/biophysical/animal/animal.py +++ b/RUFAS/biophysical/animal/animal.py @@ -2435,17 +2435,17 @@ def assess_removal_risk(self, percent_fresh: float, time: RufasTime) -> None: ------- None """ - if self.future_death_date == sys.maxsize: - if self.is_selected_for_death(percent_fresh): - self.future_death_date = self.days_born - self.cull_reason = self._future_death_reason = animal_constants.DEATH_CULL - self.dead_at_day = time.simulation_day - if self.future_cull_date == sys.maxsize: if self.is_selected_for_acute_sale(percent_fresh): self.future_cull_date = self.days_born self.cull_reason = animal_constants.ACUTE_SALE_CULL self.sold_at_day = time.simulation_day + return + if self.future_death_date == sys.maxsize: + if self.is_selected_for_death(percent_fresh): + self.future_death_date = self.days_born + self.cull_reason = self._future_death_reason = animal_constants.DEATH_CULL + self.dead_at_day = time.simulation_day def _parity_index(self) -> int: """Return the 0-based index into a by-parity array, capping parity 4+ at the last entry.""" From 8309ba8605e936d5fa427f5b9b42825aaae1c825 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 17 Sep 2026 15:17:29 +0000 Subject: [PATCH 10/13] Apply Black Formatting From 69ad9b34194fbc7d743271b761faa5b2b2bb2579 Mon Sep 17 00:00:00 2001 From: Allister Liu Date: Mon, 28 Sep 2026 11:16:14 -0400 Subject: [PATCH 11/13] Julie's comment --- RUFAS/biophysical/animal/animal.py | 73 +++++++---- RUFAS/biophysical/animal/animal_config.py | 39 ++++-- RUFAS/biophysical/animal/animal_constants.py | 6 +- .../animal/animal_module_reporter.py | 6 +- .../animal/data_types/animal_enums.py | 9 ++ .../animal/data_types/herd_statistics.py | 10 +- RUFAS/biophysical/animal/herd_manager.py | 96 ++++++++++++-- RUFAS/input/metadata/properties/default.json | 38 ++++-- .../data/animal/example_freestall_animal.json | 6 +- .../data/animal/example_open_lot_animal.json | 6 +- .../animal_cross_validation.json | 58 +++++++++ .../test_animal/test_animal/test_animal.py | 123 ++++++++++++++---- .../test_animal/test_animal_config.py | 26 +++- .../test_animal_module_reporter.py | 4 +- .../test_data_types/test_herd_statistics.py | 4 +- .../test_herd_manager/pytest_fixtures.py | 11 +- .../test_herd_manager_daily_routines.py | 82 +++++++++++- .../test_herd_manager_herd_statistics.py | 28 ++-- 18 files changed, 493 insertions(+), 132 deletions(-) diff --git a/RUFAS/biophysical/animal/animal.py b/RUFAS/biophysical/animal/animal.py index a56f6c5441..66683b852f 100644 --- a/RUFAS/biophysical/animal/animal.py +++ b/RUFAS/biophysical/animal/animal.py @@ -10,7 +10,7 @@ from RUFAS.biophysical.animal.animal_config import AnimalConfig from RUFAS.biophysical.animal.animal_genetics.animal_genetics import Genetics from RUFAS.biophysical.animal.animal_module_constants import AnimalModuleConstants -from RUFAS.biophysical.animal.data_types.animal_enums import Breed, Sex, AnimalStatus +from RUFAS.biophysical.animal.data_types.animal_enums import Breed, Sex, AnimalStatus, CowParity from RUFAS.biophysical.animal.data_types.animal_events import AnimalEvents from RUFAS.biophysical.animal.data_types.body_weight_history import BodyWeightHistory from RUFAS.biophysical.animal.data_types.daily_routines_output import DailyRoutinesOutput @@ -2419,7 +2419,13 @@ def _get_cow_values(self) -> CowValuesTypedDict: parity=self.calves, ) - def assess_removal_risk(self, percent_fresh: float, time: RufasTime) -> None: + def assess_removal_risk( + self, + percent_fresh: float, + annual_death_risk_by_parity: dict[CowParity, float], + annual_acute_sale_risk_by_parity: dict[CowParity, float], + time: RufasTime, + ) -> None: """ Roll a cow's daily mortality and acute-sale risk and schedule any resulting removal. @@ -2428,6 +2434,10 @@ def assess_removal_risk(self, percent_fresh: float, time: RufasTime) -> None: percent_fresh : float Fraction of the herd currently fresh, used to scale the daily death and acute-sale selection probabilities. + annual_death_risk_by_parity : dict[CowParity, float] + Annual death risk for each parity group, keyed by :class:`CowParity`, (unitless). + annual_acute_sale_risk_by_parity : dict[CowParity, float] + Annual acute-sale risk for each parity group, keyed by :class:`CowParity`, (unitless). time : RufasTime Current simulation time, used to record the day of death or sale. @@ -2436,56 +2446,63 @@ def assess_removal_risk(self, percent_fresh: float, time: RufasTime) -> None: None """ if self.future_cull_date == sys.maxsize: - if self.is_selected_for_acute_sale(percent_fresh): + if self.is_selected_for_removal(annual_acute_sale_risk_by_parity[self.parity_index], percent_fresh): self.future_cull_date = self.days_born self.cull_reason = animal_constants.ACUTE_SALE_CULL self.sold_at_day = time.simulation_day return if self.future_death_date == sys.maxsize: - if self.is_selected_for_death(percent_fresh): + if self.is_selected_for_removal(annual_death_risk_by_parity[self.parity_index], percent_fresh): self.future_death_date = self.days_born self.cull_reason = self._future_death_reason = animal_constants.DEATH_CULL self.dead_at_day = time.simulation_day - def _parity_index(self) -> int: - """Return the 0-based index into a by-parity array, capping parity 4+ at the last entry.""" - return 3 if self.calves >= 4 else self.calves - 1 - - def is_selected_for_death(self, percent_fresh: float) -> bool: + @property + def parity_index(self) -> CowParity: """ - Roll the cow's daily mortality risk. + The parity group of the animal. Returns ------- - bool - ``True`` if the cow is selected to die, ``False`` otherwise. + CowParity + ``CowParity.ONE`` or ``CowParity.TWO`` for cows in their 1st or 2nd lactation, + ``CowParity.THREE_PLUS`` for cows in any later lactation, and ``CowParity.NONE`` for + animals that are not cows. """ - average_daily_death_rate = AnimalConfig.parity_death_probability[self._parity_index()] / 365 - percent_other = 1 - percent_fresh - - daily_death_risk_other = average_daily_death_rate / (2.55 * percent_fresh + percent_other) - daily_death_risk_fresh = 2.55 * daily_death_risk_other - - daily_death_risk = daily_death_risk_fresh if self.days_in_milk < 50 else daily_death_risk_other - return random() <= daily_death_risk + if self.animal_type.is_cow: + if self.calves == 1: + return CowParity.ONE + elif self.calves == 2: + return CowParity.TWO + else: + return CowParity.THREE_PLUS + return CowParity.NONE - def is_selected_for_acute_sale(self, percent_fresh: float) -> bool: + def is_selected_for_removal(self, annual_risk: float, percent_fresh: float) -> bool: """ - Roll the cow's daily acute-sale (forced / involuntary) risk. + Roll the cow's daily removal (death or acute-sale) risk. + + Parameters + ---------- + annual_risk : float + Annual removal risk for the cow's parity group, (unitless). + percent_fresh : float + Fraction of the herd currently fresh, used to scale the daily risk so that fresh + cows carry a higher risk than the rest of the herd. Returns ------- bool - ``True`` if the cow is selected for an acute sale, ``False`` otherwise. + ``True`` if the cow is selected for removal, ``False`` otherwise. """ - average_daily_removal_rate = AnimalConfig.parity_acute_sale_probability[self._parity_index()] / 365 + average_daily_risk = annual_risk / 365 percent_other = 1 - percent_fresh - daily_removal_risk_other = average_daily_removal_rate / (2.55 * percent_fresh + percent_other) - daily_removal_risk_fresh = 2.55 * daily_removal_risk_other + daily_risk_other = average_daily_risk / (2.55 * percent_fresh + percent_other) + daily_risk_fresh = 2.55 * daily_risk_other - daily_removal_risk = daily_removal_risk_fresh if self.days_in_milk < 50 else daily_removal_risk_other - return random() <= daily_removal_risk + daily_risk = daily_risk_fresh if self.days_in_milk < 50 else daily_risk_other + return random() <= daily_risk def update_pen_history(self, current_pen: int, current_day: int, animal_types_in_pen: set[AnimalType]) -> None: """ diff --git a/RUFAS/biophysical/animal/animal_config.py b/RUFAS/biophysical/animal/animal_config.py index 359f9211b7..7c9db3443a 100644 --- a/RUFAS/biophysical/animal/animal_config.py +++ b/RUFAS/biophysical/animal/animal_config.py @@ -1,5 +1,6 @@ from typing import Any +from RUFAS.biophysical.animal.data_types.animal_enums import CowParity from RUFAS.biophysical.animal.data_types.repro_protocol_enums import ( HeiferReproductionProtocol, CowReproductionProtocol, @@ -145,11 +146,14 @@ class AnimalConfig: Third pregnancy check day post-breeding, (simulation day). third_pregnancy_check_loss_rate : float Pregnancy loss probability during the third pregnancy check, (unitless). - parity_death_probability : list[float] - Annual, parity-indexed probability that a cow dies during a given year, (unitless). - parity_acute_sale_probability : list[float] - Annual, parity-indexed probability that a cow is sold for an acute / involuntary - reason during a given year, (unitless). + annual_death_probability : float + Annual probability that a cow in the whole cow herd dies during a given year, (unitless). + parity_death_distribution : dict[CowParity, float] + Fractions of all cow deaths that come from the 1st, 2nd, and 3rd+ parity groups, (unitless). + annual_sale_probability : float + Annual probability that a cow in the whole cow herd is sold during a given year, (unitless). + parity_sale_distribution : dict[CowParity, float] + Fractions of all cow sales that come from the 1st, 2nd, and 3rd+ parity groups, (unitless). methane_mitigation_method : str The mitigation method applied for methane reduction, e.g., "None", (unitless). methane_mitigation_additive_amount : float @@ -240,8 +244,18 @@ class AnimalConfig: third_pregnancy_check_day: int = 200 third_pregnancy_check_loss_rate: float = 0.017 - parity_death_probability: list[float] = [0.039, 0.056, 0.085, 0.117] - parity_acute_sale_probability: list[float] = [0.169, 0.233, 0.301, 0.408] + annual_death_probability: float = 0.05 + parity_death_distribution: dict[CowParity, float] = { + CowParity.ONE: 0.2, + CowParity.TWO: 0.3, + CowParity.THREE_PLUS: 0.5, + } + annual_sale_probability: float = 0.32 + parity_sale_distribution: dict[CowParity, float] = { + CowParity.ONE: 0.2, + CowParity.TWO: 0.3, + CowParity.THREE_PLUS: 0.5, + } methane_model: dict[str, Any] = { "calves": "Pattanaik", @@ -382,8 +396,15 @@ def initialize_animal_config(cls) -> None: cls.third_pregnancy_check_day = animal_config_data["from_literature"]["repro"]["preg_check_day_3"] cls.third_pregnancy_check_loss_rate = animal_config_data["from_literature"]["repro"]["preg_loss_rate_3"] - cls.parity_death_probability = animal_config_data["from_literature"]["culling"]["parity_death_prob"] - cls.parity_acute_sale_probability = animal_config_data["from_literature"]["culling"]["parity_acute_sale_prob"] + culling_data = animal_config_data["from_literature"]["culling"] + cls.annual_death_probability = culling_data["annual_death_prob"] + cls.parity_death_distribution[CowParity.ONE] = culling_data["parity_death_dist"][0] + cls.parity_death_distribution[CowParity.TWO] = culling_data["parity_death_dist"][1] + cls.parity_death_distribution[CowParity.THREE_PLUS] = culling_data["parity_death_dist"][2] + cls.annual_sale_probability = culling_data["annual_sale_prob"] + cls.parity_sale_distribution[CowParity.ONE] = culling_data["parity_sale_dist"][0] + cls.parity_sale_distribution[CowParity.TWO] = culling_data["parity_sale_dist"][1] + cls.parity_sale_distribution[CowParity.THREE_PLUS] = culling_data["parity_sale_dist"][2] cls.methane_model = animal_data["methane_model"] methane_mitigation_data = animal_data["methane_mitigation"] diff --git a/RUFAS/biophysical/animal/animal_constants.py b/RUFAS/biophysical/animal/animal_constants.py index cf905ffc77..3410065574 100644 --- a/RUFAS/biophysical/animal/animal_constants.py +++ b/RUFAS/biophysical/animal/animal_constants.py @@ -44,7 +44,7 @@ # heifer repro INJECT_CIDR = "inject CIDR" -# oversupply cull ranking criteria +# low production cull ranking criteria CULL_RANKING_CRITERIA_MILK = "milk" """Cull ranking criterion that ranks eligible cows by their daily milk production.""" CULL_RANKING_CRITERIA_305_DAY_MILK = "305_day_milk" @@ -94,9 +94,11 @@ # culling HEIFER_REPRO_CULL = "culled for heifer reproductive problem" -OVERSUPPLY_CULL = "culled for herd resize" +LOW_PRODUCTION_CULL = "culled for low production" DEATH_CULL = "culled for death" ACUTE_SALE_CULL = "culled for acute sale" +ACUTE_SALE_FRACTION = 0.3 +"""Fraction of all annual cow sales that are acute (forced / involuntary) sales, (unitless).""" # youngstock mortality (a loss from death, not a cull) CALF_MORTALITY_LOSS = "died from pre-wean mortality" diff --git a/RUFAS/biophysical/animal/animal_module_reporter.py b/RUFAS/biophysical/animal/animal_module_reporter.py index 995844703d..7b2cf12ba0 100644 --- a/RUFAS/biophysical/animal/animal_module_reporter.py +++ b/RUFAS/biophysical/animal/animal_module_reporter.py @@ -596,8 +596,8 @@ def report_herd_statistics_data(cls, herd_statistics: HerdStatistics, simulation "is_daily_variable": True, } om.add_variable( - "sold_cow_oversupply_num", - herd_statistics.sold_cow_oversupply_num, + "sold_cow_low_production_num", + herd_statistics.sold_cow_low_production_num, dict(info_map, **{"units": MeasurementUnits.ANIMALS}), ) om.add_variable( @@ -841,7 +841,7 @@ def report_herd_statistics_data(cls, herd_statistics: HerdStatistics, simulation ) cull_reason_stats_units = { animal_constants.DEATH_CULL: MeasurementUnits.UNITLESS, - animal_constants.OVERSUPPLY_CULL: MeasurementUnits.UNITLESS, + animal_constants.LOW_PRODUCTION_CULL: MeasurementUnits.UNITLESS, animal_constants.ACUTE_SALE_CULL: MeasurementUnits.UNITLESS, } om.add_variable( diff --git a/RUFAS/biophysical/animal/data_types/animal_enums.py b/RUFAS/biophysical/animal/data_types/animal_enums.py index ee87620cf2..61ab3aae70 100644 --- a/RUFAS/biophysical/animal/data_types/animal_enums.py +++ b/RUFAS/biophysical/animal/data_types/animal_enums.py @@ -23,3 +23,12 @@ class AnimalStatus(Enum): DEAD = "dead" SOLD = "sold" STILLBORN = "stillborn" + + +class CowParity(Enum): + """Enum indicating the parity group of the animal, used to look up parity-group-specific removal risks.""" + + NONE = "none" + ONE = "1" + TWO = "2" + THREE_PLUS = "3+" diff --git a/RUFAS/biophysical/animal/data_types/herd_statistics.py b/RUFAS/biophysical/animal/data_types/herd_statistics.py index 0859744fcd..e162bec6f0 100644 --- a/RUFAS/biophysical/animal/data_types/herd_statistics.py +++ b/RUFAS/biophysical/animal/data_types/herd_statistics.py @@ -49,7 +49,7 @@ class HerdStatistics: Number of stillborn calves during a specific period, (unitless). sold_calf_num : int Number of calves sold during a specific period, (unitless). - sold_cow_oversupply_num : int + sold_cow_low_production_num : int Number of surplus cow sold, (unitless). bought_heifer_num : int Number of heifers purchased during a specific period, (unitless). @@ -184,7 +184,7 @@ class HerdStatistics: stillborn_calf_num = 0 sold_calf_num = 0 - sold_cow_oversupply_num = 0 + sold_cow_low_production_num = 0 bought_heifer_num = 0 sold_heiferII_num = 0 cow_herd_exit_num = 0 @@ -266,7 +266,7 @@ def __init__(self) -> None: } self.cull_reason_stats = { animal_constants.DEATH_CULL: 0, - animal_constants.OVERSUPPLY_CULL: 0, + animal_constants.LOW_PRODUCTION_CULL: 0, animal_constants.ACUTE_SALE_CULL: 0, } self.parity_culling_stats_range = {"1": 0, "2": 0, "3": 0, "4": 0, "5": 0, "greater_than_5": 0} @@ -275,7 +275,7 @@ def __init__(self) -> None: self.avg_age_for_parity = {"1": 0, "2": 0, "3": 0, "4": 0, "5": 0, "greater_than_5": 0} self.cull_reason_stats_percent = { animal_constants.DEATH_CULL: 0.0, - animal_constants.OVERSUPPLY_CULL: 0.0, + animal_constants.LOW_PRODUCTION_CULL: 0.0, animal_constants.ACUTE_SALE_CULL: 0.0, } self.percent_cow_for_parity = { @@ -320,7 +320,7 @@ def reset_daily_stats(self) -> None: self.stillborn_calf_num = 0 self.sold_calf_num = 0 - self.sold_cow_oversupply_num = 0 + self.sold_cow_low_production_num = 0 self.bought_heifer_num = 0 self.sold_heiferII_num = 0 self.cow_herd_exit_num = 0 diff --git a/RUFAS/biophysical/animal/herd_manager.py b/RUFAS/biophysical/animal/herd_manager.py index 699a976920..f3e416b68e 100644 --- a/RUFAS/biophysical/animal/herd_manager.py +++ b/RUFAS/biophysical/animal/herd_manager.py @@ -11,7 +11,7 @@ from RUFAS.biophysical.animal.animal_module_constants import AnimalModuleConstants from RUFAS.biophysical.animal.animal_module_reporter import AnimalModuleReporter from RUFAS.biophysical.animal.calf_retention_policy import CalfRetentionPolicy -from RUFAS.biophysical.animal.data_types.animal_enums import AnimalStatus +from RUFAS.biophysical.animal.data_types.animal_enums import AnimalStatus, CowParity from RUFAS.biophysical.animal.data_types.animal_events import AnimalEvents from RUFAS.biophysical.animal.data_types.animal_population import AnimalPopulation from RUFAS.biophysical.animal.data_types.animal_typed_dicts import ( @@ -107,13 +107,13 @@ class HerdManager: buying_threshold : int | float Herd size threshold below which replacement animals may be purchased. cull_eligibility_minimum_days_in_milk : int | float - Minimum days in milk a cow must have reached to be eligible for an oversupply cull, + Minimum days in milk a cow must have reached to be eligible for a low production cull, (simulation days). Protects fresh cows from being sold. cull_eligibility_maximum_days_carried_calf : int | float Maximum days carrying a calf (days in pregnancy) a cow may have to remain eligible for an - oversupply cull, (simulation days). Protects late-pregnant cows from being sold. + low production cull, (simulation days). Protects late-pregnant cows from being sold. cull_ranking_criteria : str - Attribute used to rank eligible cows when selecting which to sell for an oversupply cull. + Attribute used to rank eligible cows when selecting which to sell for a low production cull. One of ``"milk"`` (daily milk production) or ``"305_day_milk"`` (305-day milk yield); the lowest-ranked eligible cows are sold first. housing : dict[str, Any] @@ -569,13 +569,59 @@ def _perform_daily_routines_for_animals( animal.update_genetic_history(simulation_day=time.simulation_day) return (graduated_animals, sold_animals, stillborn_newborn_calves, newborn_calves, sold_newborn_calves) + def _calculate_annual_risk_by_parity( + self, + herd_annual_risk: float, + parity_distribution: dict[CowParity, float], + parity_group_fractions: dict[CowParity, float], + ) -> dict[CowParity, float]: + """ + Converts a whole-herd annual risk into an annual risk for each parity group. + + Parameters + ---------- + herd_annual_risk : float + Annual risk of the event for the whole cow herd, (unitless). + parity_distribution : dict[CowParity, float] + Fractions of all events that come from each parity group, keyed by :class:`CowParity`, + (unitless). + parity_group_fractions : dict[CowParity, float] + Fractions of the cow herd currently in each parity group, keyed by :class:`CowParity`, + (unitless). + + Returns + ------- + dict[CowParity, float] + Annual risk of the event for a cow in each of the ``ONE``, ``TWO``, and ``THREE_PLUS`` + parity groups, (unitless). A parity group with no cows is assigned a risk of 0. + + Notes + ----- + The expected number of events coming from parity group ``p`` is + ``herd_annual_risk * parity_distribution[p] * N``, spread over the + ``parity_group_fractions[p] * N`` cows in that group, so the per-cow risk for the group is + ``herd_annual_risk * parity_distribution[p] / parity_group_fractions[p]``. + """ + annual_risk_by_parity: dict[CowParity, float] = { + CowParity.ONE: 0.0, + CowParity.TWO: 0.0, + CowParity.THREE_PLUS: 0.0, + } + for parity in annual_risk_by_parity.keys(): + if parity_group_fractions[parity] > 0: + annual_risk_by_parity[parity] = ( + herd_annual_risk * parity_distribution[parity] / parity_group_fractions[parity] + ) + return annual_risk_by_parity + def _assess_removal_risk(self, animals: list[Animal], time: RufasTime) -> tuple[list[Animal], list[Animal]]: """ Assess daily removal risk for each cow and collect those removed. - Computes the fresh fraction across all cows in the herd, then rolls death and - acute-sale risk for every cow in ``animals`` via - :meth:`Animal.assess_removal_risk`. Non-cow animals are skipped. + Computes the fresh fraction and the 1st, 2nd, and 3rd+ parity group fractions across all + cows in the herd, converts the whole-herd annual death and acute-sale risks into annual + risks for each parity group, then rolls death and acute-sale risk for every cow in + ``animals`` via :meth:`Animal.assess_removal_risk`. Non-cow animals are skipped. Parameters ---------- @@ -594,15 +640,36 @@ def _assess_removal_risk(self, animals: list[Animal], time: RufasTime) -> tuple[ dead_cows: list[Animal] = [] all_cows = [animal for animal in self.all_animals if animal.animal_type.is_cow] + num_cows = len(all_cows) fresh_cows: list[Animal] = [cow for cow in all_cows if cow.days_in_milk < 50] - percent_fresh_cows = len(fresh_cows) / len(all_cows) if len(all_cows) > 0 else 0 + percent_fresh_cows = len(fresh_cows) / num_cows if num_cows > 0 else 0 + + parity_group_counts = {CowParity.ONE: 0, CowParity.TWO: 0, CowParity.THREE_PLUS: 0} + for cow in all_cows: + parity_group_counts[cow.parity_index] += 1 + parity_group_fractions: dict[CowParity, float] = { + parity: count / num_cows if num_cows > 0 else 0.0 for parity, count in parity_group_counts.items() + } + + annual_death_risk_by_parity = self._calculate_annual_risk_by_parity( + AnimalConfig.annual_death_probability, AnimalConfig.parity_death_distribution, parity_group_fractions + ) + annual_acute_sale_risk_by_parity = self._calculate_annual_risk_by_parity( + animal_constants.ACUTE_SALE_FRACTION * AnimalConfig.annual_sale_probability, + AnimalConfig.parity_sale_distribution, + parity_group_fractions, + ) + for animal in animals: if animal.animal_type.is_cow: - animal.assess_removal_risk(percent_fresh_cows, time) + animal.assess_removal_risk( + percent_fresh_cows, annual_death_risk_by_parity, annual_acute_sale_risk_by_parity, time + ) if animal.sold: sold_cows.append(animal) if animal.dead: dead_cows.append(animal) + self.herd_statistics.animals_deaths_by_stage[animal.animal_type] += 1 return (sold_cows, dead_cows) def _update_genetic_values_at_lactation_start(self, animal: Animal, time: RufasTime) -> None: @@ -1020,7 +1087,8 @@ def _create_newborn_calf(self, newborn_calf_config: NewBornCalfValuesTypedDict, return newborn_calf def _cull_ranking_value(self, cow: Animal) -> float: - """Returns the value used to rank ``cow`` for an oversupply cull, based on the user-defined + """ + Returns the value used to rank ``cow`` for a low production cull, based on the user-defined ``cull_ranking_criteria``. Parameters @@ -1051,7 +1119,7 @@ def _cull_ranking_value(self, cow: Animal) -> float: def _get_cow_removal_index(self, removed_animal: list[Animal]) -> int | None: """ - Finds the index of the lowest-ranked cow that is eligible for an oversupply cull. + Finds the index of the lowest-ranked cow that is eligible for a low production cull. Eligibility is governed by the ``cull_eligibility_minimum_days_in_milk`` and ``cull_eligibility_maximum_days_carried_calf`` user inputs (protecting fresh and @@ -2260,8 +2328,10 @@ def _update_sold_and_died_cow_statistics(self, sold_and_died_cows: list[Animal]) [cow for cow in sold_and_died_cows if cow.cull_reason == cull_reason] ) - oversupply_cows_num = sum(cow.cull_reason == animal_constants.OVERSUPPLY_CULL for cow in sold_and_died_cows) - self.herd_statistics.sold_cow_oversupply_num += oversupply_cows_num + low_production_cows_num = sum( + cow.cull_reason == animal_constants.LOW_PRODUCTION_CULL for cow in sold_and_died_cows + ) + self.herd_statistics.sold_cow_low_production_num += low_production_cows_num sold_cows: list[Animal] = [cow for cow in sold_and_died_cows if cow.cull_reason != animal_constants.DEATH_CULL] self.herd_statistics.sold_cows_info += [ diff --git a/RUFAS/input/metadata/properties/default.json b/RUFAS/input/metadata/properties/default.json index abdca2691a..0b2a4e69bf 100644 --- a/RUFAS/input/metadata/properties/default.json +++ b/RUFAS/input/metadata/properties/default.json @@ -240,19 +240,19 @@ }, "cull_eligibility_minimum_days_in_milk": { "type": "number", - "description": "Minimum days in milk a cow must have reached to be eligible for an oversupply cull (protects fresh cows).", + "description": "Minimum days in milk a cow must have reached to be eligible for a low production cull (protects fresh cows).", "minimum": 0, "default": 60 }, "cull_eligibility_maximum_days_carried_calf": { "type": "number", - "description": "Maximum days carrying a calf (days pregnant) a cow may have to remain eligible for an oversupply cull (protects late-pregnant cows).", + "description": "Maximum days carrying a calf (days pregnant) a cow may have to remain eligible for a low production cull (protects late-pregnant cows).", "minimum": 0, "default": 180 }, "cull_ranking_criteria": { "type": "string", - "description": "Attribute used to rank eligible cows when selecting which to sell for an oversupply cull. One of 'milk' (daily milk production) or '305_day_milk' (305-day milk yield).", + "description": "Attribute used to rank eligible cows when selecting which to sell for a low production cull. One of 'milk' (daily milk production) or '305_day_milk' (305-day milk yield).", "pattern": "^(milk|305_day_milk)$", "default": "milk" }, @@ -621,23 +621,41 @@ }, "culling": { "type": "object", - "description": "Defines the annual, by-parity probabilities that a cow dies or is sold for an acute (involuntary) reason. The timing of each event within a lactation, previously configured here, is now a model constant (see animal_constants.py, issue #2694).", - "parity_death_prob": { + "description": "Defines the whole-herd annual death and sale probabilities for cows, and how deaths and sales are distributed across the 1st, 2nd, and 3rd+ parity groups. The fraction of sales that are acute (involuntary) is a model constant (see animal_constants.py).", + "annual_death_prob": { + "type": "number", + "description": "Annual Death Probability, Whole Cow Herd -- The probability that a cow in the herd dies during a given year.", + "minimum": 0, + "maximum": 1, + "default": 0.05 + }, + "parity_death_dist": { "type": "array", - "description": "Annual Death Probability, by Parity", + "description": "Distribution of Deaths, by Parity Group", + "minimum_length": 3, + "maximum_length": 3, "properties": { "type": "number", - "description": "Annual Death Probability, by Parity Group -- The probability that a cow of a single parity group (first lactation, second lactation, etc.) dies during a given year; a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total).", + "description": "Distribution of Deaths, by Parity Group -- The fraction of all cow deaths that come from a single parity group; a separate entry should be included for 1st, 2nd, and 3rd+ parities (3 entries total, summing to 1).", "minimum": 0, "maximum": 1 } }, - "parity_acute_sale_prob": { + "annual_sale_prob": { + "type": "number", + "description": "Annual Sale Probability, Whole Cow Herd -- The probability that a cow in the herd is sold during a given year (total sales).", + "minimum": 0, + "maximum": 1, + "default": 0.32 + }, + "parity_sale_dist": { "type": "array", - "description": "Annual Acute-Sale Probability, by Parity", + "description": "Distribution of Sales, by Parity Group", + "minimum_length": 3, + "maximum_length": 3, "properties": { "type": "number", - "description": "Annual Acute-Sale Probability, by Parity Group -- The probability that a cow of a single parity group is sold for an acute (forced / involuntary) reason during a given year, regardless of whether a replacement is available; a separate entry should be included for 1st, 2nd, 3rd, and 4th+ parities (4 entries total).", + "description": "Distribution of Sales, by Parity Group -- The fraction of all cow sales that come from a single parity group; a separate entry should be included for 1st, 2nd, and 3rd+ parities (3 entries total, summing to 1).", "minimum": 0, "maximum": 1 } diff --git a/input/data/animal/example_freestall_animal.json b/input/data/animal/example_freestall_animal.json index bc56673da5..36b2f1998f 100644 --- a/input/data/animal/example_freestall_animal.json +++ b/input/data/animal/example_freestall_animal.json @@ -112,8 +112,10 @@ "std_estrus_cycle_after_pgf": 2 }, "culling": { - "parity_death_prob": [0.039,0.056,0.085,0.117], - "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408] + "annual_death_prob": 0.05, + "parity_death_dist": [0.2, 0.3, 0.5], + "annual_sale_prob": 0.32, + "parity_sale_dist": [0.2, 0.3, 0.5] }, "life_cycle": { "still_birth_rate": 0.065 diff --git a/input/data/animal/example_open_lot_animal.json b/input/data/animal/example_open_lot_animal.json index 7d515431b5..128f96fd37 100644 --- a/input/data/animal/example_open_lot_animal.json +++ b/input/data/animal/example_open_lot_animal.json @@ -112,8 +112,10 @@ "std_estrus_cycle_after_pgf": 2 }, "culling": { - "parity_death_prob": [0.039,0.056,0.085,0.117], - "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408] + "annual_death_prob": 0.05, + "parity_death_dist": [0.2, 0.3, 0.5], + "annual_sale_prob": 0.32, + "parity_sale_dist": [0.2, 0.3, 0.5] }, "life_cycle": { "still_birth_rate": 0.065 diff --git a/input/metadata/cross_validation/animal_cross_validation.json b/input/metadata/cross_validation/animal_cross_validation.json index 66cbb7e2db..caab5d0eda 100644 --- a/input/metadata/cross_validation/animal_cross_validation.json +++ b/input/metadata/cross_validation/animal_cross_validation.json @@ -1053,6 +1053,64 @@ } ] }, + { + "description": "Sum of parity death distribution fractions must equal 1.0", + "aliases": { + "variables": { + "parity_death_dist": "animal.animal_config.from_literature.culling.parity_death_dist" + }, + "constants": { + "one": 1.0 + } + }, + "rules": [ + { + "left_hand": { + "aggregation": { + "operation": "sum", + "mode": "aggregate", + "operands": ["parity_death_dist"] + } + }, + "right_hand": { + "aggregation": { + "operation": "no_op", + "operands": ["one"] + } + }, + "relationship": "equal" + } + ] + }, + { + "description": "Sum of parity sale distribution fractions must equal 1.0", + "aliases": { + "variables": { + "parity_sale_dist": "animal.animal_config.from_literature.culling.parity_sale_dist" + }, + "constants": { + "one": 1.0 + } + }, + "rules": [ + { + "left_hand": { + "aggregation": { + "operation": "sum", + "mode": "aggregate", + "operands": ["parity_sale_dist"] + } + }, + "right_hand": { + "aggregation": { + "operation": "no_op", + "operands": ["one"] + } + }, + "relationship": "equal" + } + ] + }, { "description": "The maximum days carried calf cannot be more than gestation length", "aliases": { diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal.py b/tests/test_biophysical/test_animal/test_animal/test_animal.py index 463c69eab6..a7cd78692d 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal.py @@ -11,7 +11,7 @@ from RUFAS.biophysical.animal.animal_config import AnimalConfig from RUFAS.biophysical.animal.animal_genetics.animal_genetics import Genetics from RUFAS.biophysical.animal.animal_module_constants import AnimalModuleConstants -from RUFAS.biophysical.animal.data_types.animal_enums import Breed, Sex, AnimalStatus +from RUFAS.biophysical.animal.data_types.animal_enums import Breed, Sex, AnimalStatus, CowParity from RUFAS.biophysical.animal.data_types.animal_events import AnimalEvents from RUFAS.biophysical.animal.data_types.animal_typed_dicts import ( NewBornCalfValuesTypedDict, @@ -3032,21 +3032,103 @@ def test_get_cow_values(mock_lactating_cow: Animal) -> None: assert mock_lactating_cow._get_cow_values() == expected -def test_will_die_tomorrow_no_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - animal = mock_lactating_cow - animal.calves = 1 - animal.days_born = 150 +def test_is_selected_for_removal_not_selected(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A roll above the daily removal risk does not select the cow.""" + mock_lactating_cow.days_in_milk = 150 mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - assert animal.is_selected_for_death(percent_fresh=0.15) is False + assert mock_lactating_cow.is_selected_for_removal(annual_risk=0.05, percent_fresh=0.15) is False -def test_will_die_tomorrow_with_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - """A roll below the daily death rate (annual parity probability / 365) selects the cow.""" - animal = mock_lactating_cow - animal.calves = 5 - # daily death rate = parity_death_probability[3] (0.117) / 365 ~= 0.00032. - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) - assert animal.is_selected_for_death(percent_fresh=0.15) is True +@pytest.mark.parametrize( + "days_in_milk,roll,expected", + [ + # other: 0.365 / 365 / (2.55 * 0.2 + 0.8) ~= 0.000763; fresh: 2.55x ~= 0.001947. + (150, 0.0007, True), + (150, 0.0008, False), + (10, 0.0019, True), + (10, 0.0020, False), + ], +) +def test_is_selected_for_removal_fresh_scaling( + mock_lactating_cow: Animal, mocker: MockerFixture, days_in_milk: int, roll: float, expected: bool +) -> None: + """The daily removal risk is the annual risk / 365, scaled up for fresh cows and down for the rest.""" + mock_lactating_cow.days_in_milk = days_in_milk + mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=roll) + assert mock_lactating_cow.is_selected_for_removal(annual_risk=0.365, percent_fresh=0.2) is expected + + +@pytest.mark.parametrize( + "calves,expected", + [ + (1, CowParity.ONE), + (2, CowParity.TWO), + (3, CowParity.THREE_PLUS), + (4, CowParity.THREE_PLUS), + (7, CowParity.THREE_PLUS), + ], +) +def test_parity_index_cow(mock_lactating_cow: Animal, calves: int, expected: CowParity) -> None: + """A cow's parity maps to the 1st, 2nd, or 3rd+ parity group.""" + mock_lactating_cow.calves = calves + assert mock_lactating_cow.parity_index == expected + + +def test_parity_index_non_cow(mock_heiferI: Animal) -> None: + """Animals that are not cows have no parity group.""" + assert mock_heiferI.parity_index == CowParity.NONE + + +DEATH_RISK_BY_PARITY = {CowParity.ONE: 0.01, CowParity.TWO: 0.02, CowParity.THREE_PLUS: 0.03} +ACUTE_SALE_RISK_BY_PARITY = {CowParity.ONE: 0.04, CowParity.TWO: 0.05, CowParity.THREE_PLUS: 0.06} + + +def test_assess_removal_risk_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A cow selected for acute sale is sold today and death risk is not rolled.""" + mock_lactating_cow.calves = 2 + mock_lactating_cow.days_born = 900 + mock_lactating_cow._future_cull_date = None + mock_lactating_cow._future_death_date = None + mock_select = mocker.patch.object(mock_lactating_cow, "is_selected_for_removal", return_value=True) + time = MagicMock(simulation_day=42) + + mock_lactating_cow.assess_removal_risk(0.1, DEATH_RISK_BY_PARITY, ACUTE_SALE_RISK_BY_PARITY, time) + + mock_select.assert_called_once_with(0.05, 0.1) + assert mock_lactating_cow.future_cull_date == 900 + assert mock_lactating_cow.cull_reason == animal_constants.ACUTE_SALE_CULL + assert mock_lactating_cow.sold_at_day == 42 + assert mock_lactating_cow.future_death_date == sys.maxsize + + +def test_assess_removal_risk_death(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """A cow not selected for acute sale but selected for death dies today.""" + mock_lactating_cow.calves = 5 + mock_lactating_cow.days_born = 1500 + mock_lactating_cow._future_cull_date = None + mock_lactating_cow._future_death_date = None + mock_select = mocker.patch.object(mock_lactating_cow, "is_selected_for_removal", side_effect=[False, True]) + time = MagicMock(simulation_day=7) + + mock_lactating_cow.assess_removal_risk(0.1, DEATH_RISK_BY_PARITY, ACUTE_SALE_RISK_BY_PARITY, time) + + assert mock_select.call_args_list == [call(0.06, 0.1), call(0.03, 0.1)] + assert mock_lactating_cow.future_death_date == 1500 + assert mock_lactating_cow.cull_reason == animal_constants.DEATH_CULL + assert mock_lactating_cow.dead_at_day == 7 + assert mock_lactating_cow.future_cull_date == sys.maxsize + + +def test_assess_removal_risk_skips_pending_removals(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: + """Risks are not rolled when a removal of that type is already scheduled.""" + mock_lactating_cow.calves = 1 + mock_lactating_cow._future_cull_date = 10 + mock_lactating_cow._future_death_date = 20 + mock_select = mocker.patch.object(mock_lactating_cow, "is_selected_for_removal") + + mock_lactating_cow.assess_removal_risk(0.1, DEATH_RISK_BY_PARITY, ACUTE_SALE_RISK_BY_PARITY, MagicMock()) + + mock_select.assert_not_called() def test_setup_calf_mortality_disabled_when_rate_zero(mock_calf: Animal, mocker: MockerFixture) -> None: @@ -3209,21 +3291,6 @@ def test_setup_heifer_mortality_not_committed_when_day_already_passed( assert animal._future_death_date is None -def test_will_be_sold_tomorrow_with_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - """A roll below the daily acute-sale rate (annual parity probability / 365) selects the cow.""" - animal = mock_lactating_cow - animal.calves = 1 - # daily acute-sale rate = parity_acute_sale_probability[0] (0.169) / 365 ~= 0.00046. - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.0001) - assert animal.is_selected_for_acute_sale(percent_fresh=0.15) is True - - -def test_will_be_sold_tomorrow_no_acute_sale(mock_lactating_cow: Animal, mocker: MockerFixture) -> None: - mock_lactating_cow.calves = 1 - mocker.patch("RUFAS.biophysical.animal.animal.random", return_value=0.95) - assert mock_lactating_cow.is_selected_for_acute_sale(percent_fresh=0.15) is False - - def test_set_nutrient_standard() -> None: test_standard = NutrientStandard.NASEM Animal.set_nutrient_standard(test_standard) diff --git a/tests/test_biophysical/test_animal/test_animal/test_animal_config.py b/tests/test_biophysical/test_animal/test_animal/test_animal_config.py index 9e0a702e01..3aa4d53fd2 100644 --- a/tests/test_biophysical/test_animal/test_animal/test_animal_config.py +++ b/tests/test_biophysical/test_animal/test_animal/test_animal_config.py @@ -5,6 +5,7 @@ import pytest_mock from RUFAS.biophysical.animal.animal_config import AnimalConfig +from RUFAS.biophysical.animal.data_types.animal_enums import CowParity from RUFAS.biophysical.animal.data_types.repro_protocol_enums import ( HeiferReproductionProtocol, HeiferTAISubProtocol, @@ -130,8 +131,10 @@ def _make_base_animal_config(repro_sub_protocol: str, heifer_repro_method: str) "std_estrus_cycle_after_pgf": 2, }, "culling": { - "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], + "annual_death_prob": 0.05, + "parity_death_dist": [0.2, 0.3, 0.5], + "annual_sale_prob": 0.32, + "parity_sale_dist": [0.2, 0.3, 0.5], }, "life_cycle": {"still_birth_rate": 0.065}, }, @@ -200,6 +203,19 @@ def get_data_side_effect(key: str) -> Any: assert AnimalConfig.cow_ovsynch_method == CowTAISubProtocol("OvSynch 56") assert AnimalConfig.cow_resynch_method == CowReSynchSubProtocol("TAIafterPD") + assert AnimalConfig.annual_death_probability == 0.05 + assert AnimalConfig.parity_death_distribution == { + CowParity.ONE: 0.2, + CowParity.TWO: 0.3, + CowParity.THREE_PLUS: 0.5, + } + assert AnimalConfig.annual_sale_probability == 0.32 + assert AnimalConfig.parity_sale_distribution == { + CowParity.ONE: 0.2, + CowParity.TWO: 0.3, + CowParity.THREE_PLUS: 0.5, + } + mock_om.add_warning.assert_not_called() @@ -385,8 +401,10 @@ def test_initialize_animal_config_adds_warning_when_third_check_after_or_on_dryo "std_estrus_cycle_after_pgf": 2, }, "culling": { - "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], + "annual_death_prob": 0.05, + "parity_death_dist": [0.2, 0.3, 0.5], + "annual_sale_prob": 0.32, + "parity_sale_dist": [0.2, 0.3, 0.5], }, "life_cycle": {"still_birth_rate": 0.065}, }, diff --git a/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py b/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py index 675da98c3c..c0874b27dc 100644 --- a/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py +++ b/tests/test_biophysical/test_animal/test_animal_module_reporter/test_animal_module_reporter.py @@ -801,7 +801,7 @@ def test_report_herd_statistics_data(mocker: MockerFixture) -> None: hs = HerdStatistics() # Set distinct non-zero values so each assertion catches a wrong-field mapping. - hs.sold_cow_oversupply_num = 1 + hs.sold_cow_low_production_num = 1 hs.bought_heifer_num = 2 hs.sold_heiferII_num = 3 hs.cow_herd_exit_num = 4 @@ -859,7 +859,7 @@ def test_report_herd_statistics_data(mocker: MockerFixture) -> None: reported = {c.args[0]: c.args[1] for c in mock_om_add_variable.call_args_list} # --- event counts --- - assert reported["sold_cow_oversupply_num"] == 1 + assert reported["sold_cow_low_production_num"] == 1 assert reported["bought_heifer_num"] == 2 assert reported["sold_heiferII_num"] == 3 assert reported["cow_herd_exit_num"] == 4 diff --git a/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py b/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py index 7c68a568b6..0333d55ec0 100644 --- a/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py +++ b/tests/test_biophysical/test_animal/test_data_types/test_herd_statistics.py @@ -41,7 +41,7 @@ def test_reset_daily_stats(herd_statistics: HerdStatistics) -> None: # --- event counts --- herd_statistics.stillborn_calf_num = 1 herd_statistics.sold_calf_num = 2 - herd_statistics.sold_cow_oversupply_num = 3 + herd_statistics.sold_cow_low_production_num = 3 herd_statistics.bought_heifer_num = 4 herd_statistics.sold_heiferII_num = 5 herd_statistics.cow_herd_exit_num = 6 @@ -131,7 +131,7 @@ def test_reset_daily_stats(herd_statistics: HerdStatistics) -> None: # --- event counts --- assert herd_statistics.stillborn_calf_num == 0 assert herd_statistics.sold_calf_num == 0 - assert herd_statistics.sold_cow_oversupply_num == 0 + assert herd_statistics.sold_cow_low_production_num == 0 assert herd_statistics.bought_heifer_num == 0 assert herd_statistics.sold_heiferII_num == 0 assert herd_statistics.cow_herd_exit_num == 0 diff --git a/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py b/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py index 30c6f03f32..efcfce1ebc 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py @@ -6,6 +6,7 @@ from RUFAS.biophysical.animal import animal_constants from RUFAS.biophysical.animal.animal import Animal +from RUFAS.biophysical.animal.data_types.animal_enums import CowParity from RUFAS.biophysical.animal.data_types.animal_events import AnimalEvents from RUFAS.biophysical.animal.data_types.animal_population import AnimalPopulation from RUFAS.biophysical.animal.data_types.animal_typed_dicts import SoldAnimalTypedDict, StillbornCalfTypedDict @@ -141,8 +142,10 @@ def animal_json() -> dict[str, Any]: "std_estrus_cycle_after_pgf": 2, }, "culling": { - "parity_death_prob": [0.039, 0.056, 0.085, 0.117], - "parity_acute_sale_prob": [0.169, 0.233, 0.301, 0.408], + "annual_death_prob": 0.05, + "parity_death_dist": [0.2, 0.3, 0.5], + "annual_sale_prob": 0.32, + "parity_sale_dist": [0.2, 0.3, 0.5], }, "life_cycle": {"still_birth_rate": 0.065}, }, @@ -564,6 +567,10 @@ def mock_animal( animal.sold = sold animal.stillborn = stillborn animal.calves = calves + if animal_type.is_cow: + animal.parity_index = {1: CowParity.ONE, 2: CowParity.TWO}.get(calves, CowParity.THREE_PLUS) + else: + animal.parity_index = CowParity.NONE animal.calving_interval = calving_interval animal.sold_at_day = sold_at_day animal.stillborn_day = stillborn_day diff --git a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py index e9c95a75a0..7b34cba90b 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py @@ -10,7 +10,7 @@ from RUFAS.biophysical.animal.animal_config import AnimalConfig from RUFAS.biophysical.animal.animal_genetics.animal_genetics import Genetics from RUFAS.biophysical.animal.bedding.bedding import Bedding -from RUFAS.biophysical.animal.data_types.animal_enums import AnimalStatus, Breed +from RUFAS.biophysical.animal.data_types.animal_enums import AnimalStatus, Breed, CowParity from RUFAS.biophysical.animal.data_types.animal_events import AnimalEvents from RUFAS.biophysical.animal.data_types.animal_population import AnimalPopulation from RUFAS.biophysical.animal.data_types.animal_typed_dicts import NewBornCalfValuesTypedDict @@ -473,7 +473,7 @@ def test_apply_daily_herd_structure_updates( graduated_animals = mock_herd["heiferIs"] newborn_calves = mock_herd["calves"] removed_animals = [mock_animal(AnimalType.LAC_COW, sold=True)] - sold_oversupply_cows = [mock_animal(AnimalType.LAC_COW, sold=True)] + sold_low_production_cows = [mock_animal(AnimalType.LAC_COW, sold=True)] replacement_heifers = [mock_animal(AnimalType.HEIFER_III)] mock_available_feeds: list[Feed] = [MagicMock(auto_spec=Feed)] mock_current_day_conditions = MagicMock(auto_spec=CurrentDayConditions) @@ -484,7 +484,7 @@ def test_apply_daily_herd_structure_updates( herd_manager.adjustment_period = 15 mock_check_if_cows_need_to_be_sold = mocker.patch.object( - herd_manager, "_check_if_cows_need_to_be_sold", return_value=sold_oversupply_cows + herd_manager, "_check_if_cows_need_to_be_sold", return_value=sold_low_production_cows ) mock_update_sold_and_died_cow_statistics = mocker.patch.object(herd_manager, "_update_sold_and_died_cow_statistics") mock_check_if_replacement_heifers_needed = mocker.patch.object( @@ -671,7 +671,7 @@ def test_daily_routines(herd_manager: HerdManager, mock_herd: dict[str, list[Ani [mock_animal(AnimalType.CALF, sold=False) for _ in range(2)], [mock_animal(AnimalType.CALF, sold=True) for _ in range(2)], ) - sold_oversupply_heiferIIIs = [mock_animal(AnimalType.HEIFER_III, sold=True) for _ in range(5)] + sold_low_production_heiferIIIs = [mock_animal(AnimalType.HEIFER_III, sold=True) for _ in range(5)] bought_replacement_heiferIIIs = [mock_animal(AnimalType.HEIFER_III, sold=False) for _ in range(5)] mock_perform_daily_routines_for_animals_side_effect: list[ @@ -692,7 +692,7 @@ def test_daily_routines(herd_manager: HerdManager, mock_herd: dict[str, list[Ani ) mock_update_sold_animal_statistics = mocker.patch.object(herd_manager, "_update_sold_animal_statistics") mock_check_if_cows_need_to_be_sold = mocker.patch.object( - herd_manager, "_check_if_cows_need_to_be_sold", return_value=sold_oversupply_heiferIIIs + herd_manager, "_check_if_cows_need_to_be_sold", return_value=sold_low_production_heiferIIIs ) mock_check_if_replacement_heifers_needed = mocker.patch.object( herd_manager, "_check_if_replacement_heifers_needed", return_value=bought_replacement_heiferIIIs @@ -846,7 +846,7 @@ def test_check_if_cows_need_to_be_sold_comprehensive(herd_manager: HerdManager, herd_manager.herd_statistics.herd_num = HERD_TARGET herd_manager.selling_threshold = SELLING_THRESHOLD herd_manager.herd_statistics.cow_num = 15 - herd_manager.herd_statistics.sold_cow_oversupply_num = 0 + herd_manager.herd_statistics.sold_cow_low_production_num = 0 herd_manager.herd_statistics.sold_cow_num = 0 herd_manager.herd_statistics.cow_herd_exit_num = 10 @@ -1343,3 +1343,73 @@ def test_get_cow_removal_index_invalid_criteria_raises(herd_manager: HerdManager with pytest.raises(ValueError, match="Invalid cull_ranking_criteria"): herd_manager._get_cow_removal_index([]) + + +def _by_parity(one: float, two: float, three_plus: float) -> dict[CowParity, float]: + """Builds a dict keyed by the 1st, 2nd, and 3rd+ parity groups.""" + return {CowParity.ONE: one, CowParity.TWO: two, CowParity.THREE_PLUS: three_plus} + + +@pytest.mark.parametrize( + "herd_annual_risk,parity_distribution,parity_group_fractions,expected", + [ + (0.05, _by_parity(0.2, 0.3, 0.5), _by_parity(0.4, 0.3, 0.3), _by_parity(0.025, 0.05, 0.05 * 0.5 / 0.3)), + (0.096, _by_parity(0.2, 0.3, 0.5), _by_parity(0.2, 0.3, 0.5), _by_parity(0.096, 0.096, 0.096)), + (0.05, _by_parity(0.2, 0.3, 0.5), _by_parity(0.0, 0.5, 0.5), _by_parity(0.0, 0.03, 0.05)), + (0.05, _by_parity(0.2, 0.3, 0.5), _by_parity(0.0, 0.0, 0.0), _by_parity(0.0, 0.0, 0.0)), + ], +) +def test_calculate_annual_risk_by_parity( + herd_manager: HerdManager, + herd_annual_risk: float, + parity_distribution: dict[CowParity, float], + parity_group_fractions: dict[CowParity, float], + expected: dict[CowParity, float], +) -> None: + """Whole-herd annual risk is split across parity groups by the event distribution and group size.""" + actual = herd_manager._calculate_annual_risk_by_parity( + herd_annual_risk, parity_distribution, parity_group_fractions + ) + + assert actual == pytest.approx(expected) + + +def test_assess_removal_risk(herd_manager: HerdManager, mocker: MockerFixture) -> None: + """Cows are assessed with the herd fresh fraction and parity-group risks; removed cows are collected.""" + fresh_first_parity_cow = mock_animal(AnimalType.LAC_COW, days_in_milk=10, calves=1) + sold_cow = mock_animal(AnimalType.LAC_COW, days_in_milk=100, calves=2) + dead_cow = mock_animal(AnimalType.DRY_COW, days_in_milk=300, calves=5) + heifer = mock_animal(AnimalType.HEIFER_III) + for cow in [fresh_first_parity_cow, sold_cow, dead_cow]: + cow.sold = False + cow.dead = False + sold_cow.assess_removal_risk.side_effect = lambda *_: setattr(sold_cow, "sold", True) + dead_cow.assess_removal_risk.side_effect = lambda *_: setattr(dead_cow, "dead", True) + cows = [fresh_first_parity_cow, sold_cow, dead_cow] + mocker.patch.object(HerdManager, "all_animals", new_callable=mocker.PropertyMock, return_value=cows + [heifer]) + mocker.patch.object(AnimalConfig, "annual_death_probability", 0.06) + mocker.patch.object(AnimalConfig, "parity_death_distribution", _by_parity(0.2, 0.3, 0.5)) + mocker.patch.object(AnimalConfig, "annual_sale_probability", 0.3) + mocker.patch.object(AnimalConfig, "parity_sale_distribution", _by_parity(0.1, 0.4, 0.5)) + herd_manager.herd_statistics.animals_deaths_by_stage[AnimalType.DRY_COW] = 0 + time = MagicMock(auto_spec=RufasTime) + + sold_cows, dead_cows = herd_manager._assess_removal_risk(cows + [heifer], time) + + assert sold_cows == [sold_cow] + assert dead_cows == [dead_cow] + assert herd_manager.herd_statistics.animals_deaths_by_stage[AnimalType.DRY_COW] == 1 + heifer.assess_removal_risk.assert_not_called() + expected_death_risk = pytest.approx(_by_parity(0.06 * 0.2 * 3, 0.06 * 0.3 * 3, 0.06 * 0.5 * 3)) + expected_acute_sale_risk = pytest.approx(_by_parity(0.09 * 0.1 * 3, 0.09 * 0.4 * 3, 0.09 * 0.5 * 3)) + for cow in cows: + cow.assess_removal_risk.assert_called_once_with( + pytest.approx(1 / 3), expected_death_risk, expected_acute_sale_risk, time + ) + + +def test_assess_removal_risk_no_cows(herd_manager: HerdManager, mocker: MockerFixture) -> None: + """With no cows in the herd, nothing is assessed or removed.""" + mocker.patch.object(HerdManager, "all_animals", new_callable=mocker.PropertyMock, return_value=[]) + + assert herd_manager._assess_removal_risk([], MagicMock(auto_spec=RufasTime)) == ([], []) diff --git a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py index 9ea03e6d39..a2b14d6567 100644 --- a/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py +++ b/tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_herd_statistics.py @@ -175,13 +175,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st ( { animal_constants.DEATH_CULL: 0, - animal_constants.OVERSUPPLY_CULL: 0, + animal_constants.LOW_PRODUCTION_CULL: 0, animal_constants.ACUTE_SALE_CULL: 0, }, 0, { animal_constants.DEATH_CULL: 0.0, - animal_constants.OVERSUPPLY_CULL: 0.0, + animal_constants.LOW_PRODUCTION_CULL: 0.0, animal_constants.ACUTE_SALE_CULL: 0.0, }, ), @@ -189,13 +189,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st ( { animal_constants.DEATH_CULL: 5, - animal_constants.OVERSUPPLY_CULL: 0, + animal_constants.LOW_PRODUCTION_CULL: 0, animal_constants.ACUTE_SALE_CULL: 0, }, 5, { animal_constants.DEATH_CULL: 100.0, - animal_constants.OVERSUPPLY_CULL: 0.0, + animal_constants.LOW_PRODUCTION_CULL: 0.0, animal_constants.ACUTE_SALE_CULL: 0.0, }, ), @@ -204,13 +204,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st ( { animal_constants.DEATH_CULL: 5, - animal_constants.OVERSUPPLY_CULL: 5, + animal_constants.LOW_PRODUCTION_CULL: 5, animal_constants.ACUTE_SALE_CULL: 0, }, 10, { animal_constants.DEATH_CULL: 50.0, - animal_constants.OVERSUPPLY_CULL: 50.0, + animal_constants.LOW_PRODUCTION_CULL: 50.0, animal_constants.ACUTE_SALE_CULL: 0.0, }, ), @@ -219,13 +219,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st ( { animal_constants.DEATH_CULL: 3, - animal_constants.OVERSUPPLY_CULL: 2, + animal_constants.LOW_PRODUCTION_CULL: 2, animal_constants.ACUTE_SALE_CULL: 0, }, 10, { animal_constants.DEATH_CULL: 30.0, - animal_constants.OVERSUPPLY_CULL: 20.0, + animal_constants.LOW_PRODUCTION_CULL: 20.0, animal_constants.ACUTE_SALE_CULL: 0.0, }, ), @@ -235,13 +235,13 @@ def test_calculate_cow_percentages(herd_manager: HerdManager, mock_herd: dict[st ( { animal_constants.DEATH_CULL: 2, - animal_constants.OVERSUPPLY_CULL: 0, + animal_constants.LOW_PRODUCTION_CULL: 0, animal_constants.ACUTE_SALE_CULL: 8, }, 10, { animal_constants.DEATH_CULL: 20.0, - animal_constants.OVERSUPPLY_CULL: 0.0, + animal_constants.LOW_PRODUCTION_CULL: 0.0, animal_constants.ACUTE_SALE_CULL: 80.0, }, ), @@ -489,7 +489,7 @@ def test_update_sold_and_died_cow_statistics( ) -> None: """Unit test for _update_sold_and_died_cow_statistics()""" cull_reasons = [ - animal_constants.OVERSUPPLY_CULL, + animal_constants.LOW_PRODUCTION_CULL, animal_constants.ACUTE_SALE_CULL, ] @@ -556,15 +556,15 @@ def test_update_sold_and_died_cow_statistics( current_cull_reason_stats = { animal_constants.DEATH_CULL: randint(0, num_total_sold_and_died_cows), - animal_constants.OVERSUPPLY_CULL: randint(0, num_total_sold_and_died_cows), + animal_constants.LOW_PRODUCTION_CULL: randint(0, num_total_sold_and_died_cows), animal_constants.ACUTE_SALE_CULL: randint(0, num_total_sold_and_died_cows), } herd_manager.herd_statistics.cull_reason_stats = current_cull_reason_stats expected_cull_reason_stats = { animal_constants.DEATH_CULL: current_cull_reason_stats[animal_constants.DEATH_CULL] + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.DEATH_CULL]), - animal_constants.OVERSUPPLY_CULL: current_cull_reason_stats[animal_constants.OVERSUPPLY_CULL] - + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.OVERSUPPLY_CULL]), + animal_constants.LOW_PRODUCTION_CULL: current_cull_reason_stats[animal_constants.LOW_PRODUCTION_CULL] + + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.LOW_PRODUCTION_CULL]), animal_constants.ACUTE_SALE_CULL: current_cull_reason_stats[animal_constants.ACUTE_SALE_CULL] + len([cow for cow in sold_and_died_cows if cow.cull_reason == animal_constants.ACUTE_SALE_CULL]), } From c66bd88a00b8b4a586092688f03060321d5b50d7 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Mon, 28 Sep 2026 15:23:40 +0000 Subject: [PATCH 12/13] Apply Black Formatting From ab155357363c586da5f114a685624ee9d76d6261 Mon Sep 17 00:00:00 2001 From: allisterakun Date: Mon, 28 Sep 2026 15:28:35 +0000 Subject: [PATCH 13/13] Update badges on README --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index ed66b0c6ba..3c09b15f78 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ [![Flake8](https://img.shields.io/badge/Flake8-passed-brightgreen)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml) [![Pytest](https://img.shields.io/badge/Pytest-passed-brightgreen)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml) [![Coverage](https://img.shields.io/badge/Coverage-99%25-brightgreen)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml) -[![Mypy](https://img.shields.io/badge/Mypy-1164%20errors-red)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml) +[![Mypy](https://img.shields.io/badge/Mypy-1170%20errors-red)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml) # RuFaS: Ruminant Farm Systems