[Animal] Maintaining herd size

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-1195%20errors-red)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml)
+[![Mypy](https://img.shields.io/badge/Mypy-1194%20errors-red)](https://github.com/RuminantFarmSystems/MASM/actions/workflows/combined_format_lint_test_mypy.yml)
 
 
 # RuFaS: Ruminant Farm Systems

RUFAS/biophysical/animal/animal.py

@@ -1748,14 +1748,6 @@ def animal_life_stage_update(self, time: RufasTime) -> tuple[AnimalStatus, NewBo
             self.cull_reason = animal_constants.DEATH_CULL
             animal_status = AnimalStatus.DEAD
 
-        if (
-            self.animal_type.is_cow
-            and self.reproduction.do_not_breed
-            and self.milk_production.daily_milk_produced < AnimalConfig.cull_milk_production
-        ):
-            self.cull_reason = animal_constants.LOW_PROD_CULL
-            self.sold_at_day = time.simulation_day
-            animal_status = AnimalStatus.SOLD
         return animal_status, newborn_calf_config
 
     def _evaluate_calf_for_heiferI(self) -> bool:

RUFAS/biophysical/animal/animal_config.py

@@ -37,8 +37,6 @@ class AnimalConfig:
         Maximum day at which a heifer is culled if not pregnant, (simulation days).
     do_not_breed_time : int
         The duration after which breeding is stopped, (simulation days).
-    cull_milk_production : int
-        The threshold milk production below which cows are culled, (simulation days).
     semen_type : str
         Types of semen used for reproduction, e.g., "conventional", (unitless).
     male_calf_rate_conventional_semen : float
@@ -193,7 +191,6 @@ class AnimalConfig:
     dry_off_day_of_pregnancy: int = 218
     heifer_reproduction_cull_day: int = 500
     do_not_breed_time: int = 185
-    cull_milk_production: int = 30
 
     semen_type: str = "conventional"
     male_calf_rate_conventional_semen: float = 0.53
@@ -392,7 +389,6 @@ def initialize_animal_config(cls) -> None:
         cls.dry_off_day_of_pregnancy = animal_config_data["management_decisions"]["days_in_preg_when_dry"]
         cls.heifer_reproduction_cull_day = animal_config_data["management_decisions"]["heifer_repro_cull_time"]
         cls.do_not_breed_time = animal_config_data["management_decisions"]["do_not_breed_time"]
-        cls.cull_milk_production = animal_config_data["management_decisions"]["cull_milk_production"]
 
         cls.semen_type = animal_config_data["management_decisions"]["semen_type"]
         cls.male_calf_rate_conventional_semen = animal_config_data["farm_level"]["calf"][

RUFAS/biophysical/animal/animal_constants.py

@@ -44,6 +44,12 @@
 # heifer repro
 INJECT_CIDR = "inject CIDR"
 
+# herd size management
+MIN_DIM_FOR_REMOVAL = 60
+"""Minimum days in milk required for a cow to be eligible for removal."""
+MAX_DAYS_IN_PREG_FOR_REMOVAL = 180
+"""Maximum pregnancy duration for a cow to be eligible for removal."""
+
 # presynch protocols
 PRESYNCH_PERIOD_START = "Presynch period started"
 PRESYNCH_PERIOD_END = "Presynch period ended"
@@ -86,7 +92,7 @@
 
 # culling
 HEIFER_REPRO_CULL = "culled for heifer reproductive problem"
-LOW_PROD_CULL = "culled for low production"
+OVERSUPPLY_CULL = "culled for herd resize"
 DEATH_CULL = "culled for death"
 LAMENESS_CULL = "culled for lameness"
 INJURY_CULL = "culled for injury"

RUFAS/biophysical/animal/animal_module_reporter.py

@@ -596,8 +596,8 @@ def report_herd_statistics_data(cls, herd_statistics: HerdStatistics, simulation
             "data_origin": [("HerdManager", "daily_update")],
         }
         om.add_variable(
-            "sold_heiferIII_oversupply_num",
-            herd_statistics.sold_heiferIII_oversupply_num,
+            "sold_cow_oversupply_num",
+            herd_statistics.sold_cow_oversupply_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.LOW_PROD_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,

RUFAS/biophysical/animal/data_types/herd_statistics.py

@@ -47,8 +47,8 @@ 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_heiferIII_oversupply_num : int
-        Number of surplus "Heifer III" animals sold, (unitless).
+    sold_cow_oversupply_num : int
+        Number of surplus cow sold, (unitless).
     bought_heifer_num : int
         Number of heifers purchased during a specific period, (unitless).
     sold_heiferII_num : int
@@ -178,7 +178,7 @@ class HerdStatistics:
 
     stillborn_calf_num = 0
     sold_calf_num = 0
-    sold_heiferIII_oversupply_num = 0
+    sold_cow_oversupply_num = 0
     bought_heifer_num = 0
     sold_heiferII_num = 0
     cow_herd_exit_num = 0
@@ -256,7 +256,7 @@ def __init__(self) -> None:
         }
         self.cull_reason_stats = {
             animal_constants.DEATH_CULL: 0,
-            animal_constants.LOW_PROD_CULL: 0,
+            animal_constants.OVERSUPPLY_CULL: 0,
             animal_constants.LAMENESS_CULL: 0,
             animal_constants.INJURY_CULL: 0,
             animal_constants.MASTITIS_CULL: 0,
@@ -270,7 +270,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.LOW_PROD_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,
@@ -313,7 +313,7 @@ def reset_daily_stats(self) -> None:
 
         self.stillborn_calf_num = 0
         self.sold_calf_num = 0
-        self.sold_heiferIII_oversupply_num = 0
+        self.sold_cow_oversupply_num = 0
         self.bought_heifer_num = 0
         self.sold_heiferII_num = 0
         self.cow_herd_exit_num = 0

RUFAS/biophysical/animal/herd_manager.py

@@ -131,6 +131,9 @@ def __init__(
 
         self.herd_statistics = HerdStatistics()
         self.herd_statistics.herd_num = animal_config_data["herd_information"]["herd_num"]
+        self.adjustment_period = animal_config_data["herd_information"]["herd_size_adjustment_period"]
+        self.selling_threshold = animal_config_data["herd_information"]["herd_size_sell_threshold"]
+        self.buying_threshold = animal_config_data["herd_information"]["herd_size_buy_threshold"]
         self.herd_reproduction_statistics = HerdReproductionStatistics()
 
         self.housing = animal_config_data["housing"]
@@ -274,19 +277,6 @@ def phosphorus_concentration_by_animal_class(self) -> dict[AnimalType, float]:
 
         return phosphorus_concentration_by_animal_class
 
-    @property
-    def current_herd_size(self) -> int:
-        """
-        Calculates the current size of the herd based on the number of heiferIIIs and cows.
-
-        Returns
-        -------
-        int
-            The current size of the herd.
-
-        """
-        return len(self.heiferIIIs) + len(self.cows)
-
     @property
     def heiferII_events_by_id(self) -> dict[str, AnimalEvents]:
         """
@@ -601,18 +591,35 @@ def daily_routines(
 
         self._update_stillborn_calf_statistics(stillborn_newborn_calves)
 
-        removed_animals += self._check_if_heifers_need_to_be_sold(simulation_day=time.simulation_day)
-        newly_added_animals = self._check_if_replacement_heifers_needed(time=time)
-        self._update_herd_structure(
-            graduated_animals=graduated_animals,
-            newborn_calves=newborn_calves,
-            newly_added_animals=newly_added_animals,
-            removed_animals=removed_animals,
-            available_feeds=available_feeds,
-            current_day_conditions=weather.get_current_day_conditions(time),
-            total_inventory=total_inventory,
-            simulation_day=time.simulation_day,
-        )
+        adjust_herd_size: bool = time.simulation_day > 0 and time.simulation_day % self.adjustment_period == 0
+        if adjust_herd_size:
+            removed_animals += self._check_if_cows_need_to_be_sold(
+                simulation_day=time.simulation_day, removed_animal=removed_animals
+            )
+            self._update_sold_and_died_cow_statistics(removed_animals)
+            newly_added_animals = self._check_if_replacement_heifers_needed(time=time)
+
+            self._update_herd_structure(
+                graduated_animals=graduated_animals,
+                newborn_calves=newborn_calves,
+                newly_added_animals=newly_added_animals,
+                removed_animals=removed_animals,
+                available_feeds=available_feeds,
+                current_day_conditions=weather.get_current_day_conditions(time),
+                total_inventory=total_inventory,
+                simulation_day=time.simulation_day,
+            )
+        else:
+            self._update_herd_structure(
+                graduated_animals=graduated_animals,
+                newborn_calves=newborn_calves,
+                newly_added_animals=[],
+                removed_animals=removed_animals,
+                available_feeds=available_feeds,
+                current_day_conditions=weather.get_current_day_conditions(time),
+                total_inventory=total_inventory,
+                simulation_day=time.simulation_day,
+            )
 
         self.record_pen_history(time.simulation_day)
         enteric_methane_emission_by_pen: dict[str, float] = {}
@@ -625,6 +632,25 @@ def daily_routines(
 
         self.update_herd_statistics()
 
+        no_milk_cow_num = len(
+            [
+                cow
+                for cow in self.cows
+                if cow.milk_production.daily_milk_produced == 0 and cow.is_milking and cow.days_in_milk > 1
+            ]
+        )
+
+        if no_milk_cow_num > 0:
+            self.om.add_warning(
+                "Warning: Lactating cows with no production.",
+                f"There are {no_milk_cow_num} lactating cows with no milking production on simulation"
+                f" day {time.simulation_day}.",
+                info_map={
+                    "class": self.__class__.__name__,
+                    "function": self.daily_routines.__name__,
+                    "simulation_day": time.simulation_day,
+                },
+            )
         AnimalModuleReporter.report_enteric_methane_emission(enteric_methane_emission_by_pen)
         AnimalModuleReporter.report_daily_animal_population(self.herd_statistics, time.simulation_day)
         AnimalModuleReporter.report_herd_statistics_data(self.herd_statistics, time.simulation_day)
@@ -657,8 +683,7 @@ def _report_ration(self, simulation_day: int) -> None:
 
     def _create_newborn_calf(self, newborn_calf_config: NewBornCalfValuesTypedDict, simulation_day: int) -> Animal:
         """
-        Creates a new newborn calf instance and records its entry event in the herd if it
-        is not sold.
+        Creates a new newborn calf instance and records its entry event in the herd if it is not sold.
 
         Parameters
         ----------
@@ -679,50 +704,49 @@ def _create_newborn_calf(self, newborn_calf_config: NewBornCalfValuesTypedDict,
             newborn_calf.events.add_event(newborn_calf.days_born, simulation_day, animal_constants.ENTER_HERD)
         return newborn_calf
 
-    def _check_if_heifers_need_to_be_sold(
-        self,
-        simulation_day: int,
-    ) -> list[Animal]:
-        """
-        Checks if surplus heifers need to be sold based on herd size.
+    def _get_cow_removal_index(self, removed_animal: list[Animal]) -> int | None:
+        """Finds the indices of cows with the lowest daily milk production among cows that meet the specified
+        days-in-milk and days-pregnant criteria."""
+        eligible_indices = []
 
-        This method evaluates if the current number of heifers and cows exceeds a
-        specified threshold (defined as 3% over the herd statistics' target
-        herd size). If the threshold is surpassed, heiferIIIs are removed from the
-        herd until the herd size falls within the acceptable range.
+        for index, cow in enumerate(self.cows):
+            if cow in removed_animal:
+                continue
+            eligible_for_removal = (
+                cow.days_in_milk > animal_constants.MIN_DIM_FOR_REMOVAL
+                and cow.days_in_pregnancy < animal_constants.MAX_DAYS_IN_PREG_FOR_REMOVAL
+            )
+            if eligible_for_removal:
+                eligible_indices.append(index)
 
-        Parameters
-        ----------
-        simulation_day : int
-            The simulation day on which the check and potential sale is conducted.
+        if not eligible_indices:
+            return None
 
-        Returns
-        -------
-        list[Animal]
-            A list of heiferIIIs to be sold.
+        return min(eligible_indices, key=lambda i: self.cows[i].milk_production.daily_milk_produced)
 
-        """
+    def _check_if_cows_need_to_be_sold(self, simulation_day: int, removed_animal: list[Animal]) -> list[Animal]:
+        """Checks if surplus cows need to be sold based on herd size."""
         animals_removed: list[Animal] = []
-        while (
-            self.current_herd_size > self.herd_statistics.herd_num * animal_constants.SELLING_THRESHOLD
-            and len(self.heiferIIIs) > 0
-        ):
-            removed_heiferIII = self.heiferIIIs.pop()
-            animals_removed.append(removed_heiferIII)
-            removed_heiferIII.sold_at_day = simulation_day
-            self.herd_statistics.sold_heiferIIIs_info.append(
-                SoldAnimalTypedDict(
-                    id=removed_heiferIII.id,
-                    animal_type=removed_heiferIII.animal_type.value,
-                    sold_at_day=removed_heiferIII.sold_at_day,
-                    body_weight=removed_heiferIII.body_weight,
-                    cull_reason="NA",
-                    days_in_milk="NA",
-                    parity="NA",
+
+        while len(self.cows) > self.selling_threshold and len(self.cows) > 0:
+            remove_index = self._get_cow_removal_index(removed_animal)
+
+            if remove_index is None:
+                info_map = {
+                    "class": self.__class__.__name__,
+                    "function": self._check_if_cows_need_to_be_sold.__name__,
+                    "simulation_day": simulation_day,
+                }
+                self.om.add_error(
+                    "Unable to adjust herd size", "There are no cow that's qualified to be sold.", info_map
                 )
-            )
-            self.herd_statistics.sold_heiferIII_oversupply_num += 1
-            self.herd_statistics.heiferIII_num -= 1
+                break
+
+            removed_cow = self.cows.pop(remove_index)
+            removed_cow.sold_at_day = simulation_day
+            removed_cow.cull_reason = "culled for herd resize"
+            animals_removed.append(removed_cow)
+
         return animals_removed
 
     def _check_if_replacement_heifers_needed(self, time: RufasTime) -> list[Animal]:
@@ -746,9 +770,7 @@ def _check_if_replacement_heifers_needed(self, time: RufasTime) -> list[Animal]:
         """
         animals_added: list[Animal] = []
         while (
-            self.current_herd_size + self.herd_statistics.bought_heifer_num
-            < self.herd_statistics.herd_num * animal_constants.BUYING_THRESHOLD
-            and time.simulation_day > 1
+            len(self.cows) + self.herd_statistics.bought_heifer_num < self.buying_threshold and time.simulation_day > 1
         ):
             if len(self.replacement_market) == 0:
                 break
@@ -1810,6 +1832,7 @@ def _update_sold_and_died_cow_statistics(self, sold_and_died_cows: list[Animal])
         sum_cow_culling_age = self.herd_statistics.avg_cow_culling_age * self.herd_statistics.cow_herd_exit_num + sum(
             [cow.days_born for cow in sold_and_died_cows]
         )
+        self.herd_statistics.cow_num -= len(sold_and_died_cows)
         self.herd_statistics.cow_herd_exit_num += len(sold_and_died_cows)
         self.herd_statistics.avg_cow_culling_age = (
             (sum_cow_culling_age / self.herd_statistics.cow_herd_exit_num)
@@ -1834,6 +1857,9 @@ 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
+
         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 += [
             SoldAnimalTypedDict(

RUFAS/biophysical/field/soil/infiltration.py

@@ -101,7 +101,7 @@ def infiltrate(self, rainfall: float) -> None:
 
     # --- static methods ---
     @staticmethod
-    def _determine_first_moisture_condition_parameter(second_moisture_condition: float):
+    def _determine_first_moisture_condition_parameter(second_moisture_condition: float) -> float:
         """
         Determine the curve number for dry (wilting point) conditions.
 
@@ -125,7 +125,7 @@ def _determine_first_moisture_condition_parameter(second_moisture_condition: flo
         return second_moisture_condition - (top / bottom)
 
     @staticmethod
-    def _determine_third_moisture_condition_parameter(second_moisture_condition: float):
+    def _determine_third_moisture_condition_parameter(second_moisture_condition: float) -> float:
         """
         Determine the curve number for wet (field capacity) conditions.