[Animal] ME305 lactation curve implementation

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

RUFAS/biophysical/animal/animal.py

@@ -206,6 +206,8 @@ def __init__(
         self.nutrition_supply: NutritionSupply = NutritionSupply.make_empty_nutrition_supply()
         self.nutrition_supply.dry_matter = AnimalModuleConstants.DEFAULT_DRY_MATTER_INTAKE
         self.previous_nutrition_supply: NutritionSupply | None = None
+        self.milk_yield_305_day: float = 0.0
+        self.mature_equivalent_milking_prediction_305_day: float = 0.0
 
         self._days_in_milk: int = 0
         self._milk_production_output_days_in_milk: int = 0
@@ -2464,6 +2466,26 @@ def calculate_nutrition_requirements(
 
         return requirements
 
+    def update_mature_equivalent_305_days_milk_production(self) -> None:
+        if self.days_in_milk < 305:
+            self.mature_equivalent_milking_prediction_305_day = self.milk_production.calculate_305_day_milk_yield(
+                self.milk_production.wood_l,
+                self.milk_production.wood_m,
+                self.milk_production.wood_n,
+                self.milk_production.milk_production_history,
+                self.days_in_milk,
+            )
+        else:
+            self.mature_equivalent_milking_prediction_305_day = (
+                self.milk_production.current_lactation_305_day_milk_produced
+            )
+
+        parity_factor = {1: 1.25, 2: 1.18}.get(self.calves, 1.0)
+
+        self.mature_equivalent_milking_prediction_305_day = (
+            self.mature_equivalent_milking_prediction_305_day * parity_factor
+        )
+
     def update_genetic_history(self, simulation_day: int) -> None:
         """
         Updates the genetic history record for the animal on the given simulation day.

RUFAS/biophysical/animal/herd_manager.py

@@ -812,7 +812,7 @@ def _get_cow_removal_index(self, removed_animal: list[Animal]) -> int | None:
         if not eligible_indices:
             return None
 
-        return min(eligible_indices, key=lambda i: self.cows[i].milk_production.daily_milk_produced)
+        return min(eligible_indices, key=lambda i: self.cows[i].mature_equivalent_milking_prediction_305_day)
 
     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."""
@@ -2175,3 +2175,7 @@ def _update_total_enteric_methane(self, digestive_outputs: list[dict[AnimalType,
                     self.herd_statistics.total_enteric_methane[animal_type] = {
                         k: float(current_totals.get(k, 0) + new_emissions.get(k, 0)) for k in all_keys
                     }
+
+    def update_milk_305_day_yield_predictions(self) -> None:
+        for cow in self.cows:
+            cow.update_mature_equivalent_305_days_milk_production()

RUFAS/biophysical/animal/milk/lactation_curve.py

@@ -377,7 +377,7 @@ def _estimate_305_day_milk_yield_by_parity(
     @staticmethod
     def _calculate_305_day_milk_yield_error(l_param: float, m_param: float, n_param: float, milk_yield: float) -> float:
         """Calculates absolute difference between an estimated 305 day milk yield and a predetermined one."""
-        return abs(MilkProduction.calc_305_day_milk_yield(l_param, m_param, n_param) - milk_yield)
+        return abs(MilkProduction.calculate_305_day_milk_yield(l_param, m_param, n_param, None) - milk_yield)
 
     @classmethod
     def _fit_wood_l_param_to_milk_yield(

RUFAS/biophysical/animal/milk/milk_production.py

@@ -260,7 +260,13 @@ def calculate_daily_milk_production(days_in_milk: int, l_param: float, m_param:
         return l_param * np.power(days_in_milk, m_param) * np.exp(-1 * n_param * days_in_milk)
 
     @staticmethod
-    def calc_305_day_milk_yield(l_param: float, m_param: float, n_param: float) -> float:
+    def calculate_305_day_milk_yield(
+        l_param: float,
+        m_param: float,
+        n_param: float,
+        milking_history: list[MilkProductionRecord] | None,
+        days_in_milk: int = 0,
+    ) -> float:
         """
         Calculates the total milk yield from day 1 to day 305 of the lactation.
 
@@ -276,16 +282,25 @@ def calc_305_day_milk_yield(l_param: float, m_param: float, n_param: float) -> f
             Wood's lactation curve parameter m.
         n_param: float
             Wood's lactation curve parameter n.
+        days_in_milk : int
+            Days in milk.
+        milking_history :
+            The milk production history if the animal.
 
         Returns
         -------
         float
             305 day milk yield for a cow with the given lactation curve (kg).
 
         """
+        production_history_sum = 0
+        if 305 > days_in_milk > 0 and milking_history is not None:
+            production_history_sum = sum(history["milk_production"] for history in milking_history[:days_in_milk])
 
-        result, _ = quad(MilkProduction.calculate_daily_milk_production, 1, 305, args=(l_param, m_param, n_param))
-        return result
+        result, _ = quad(
+            MilkProduction.calculate_daily_milk_production, days_in_milk + 1, 305, args=(l_param, m_param, n_param)
+        )
+        return result + production_history_sum
 
     def _get_milk_production_adjustment(self) -> float:
         """

RUFAS/simulation_engine.py

@@ -352,6 +352,7 @@ def _execute_ration_planning(self) -> None:
         """Checks if it's time to reformulate the ration and executes ration formulation if needed."""
         if self._is_time_to_reformulate_ration:
             self._formulate_ration()
+            self.herd_manager.update_milk_305_day_yield_predictions()
 
     @property
     def _is_time_to_reformulate_ration(self) -> bool:

tests/test_biophysical/test_animal/test_animal/test_animal.py

@@ -61,6 +61,53 @@
 from RUFAS.rufas_time import RufasTime
 
 
+def _make_cow_for_305_day_milk_prediction(
+    days_in_milk: int, calves: int, current_lactation_305_day_milk_produced: float = 0.0
+) -> Animal:
+    cow = Animal.__new__(Animal)
+    cow.animal_type = AnimalType.LAC_COW
+    cow.days_in_milk = days_in_milk
+    cow.reproduction = MagicMock()
+    cow.reproduction.calves = calves
+    cow.milk_production = MagicMock()
+    cow.milk_production.wood_l = 1.0
+    cow.milk_production.wood_m = 2.0
+    cow.milk_production.wood_n = 3.0
+    cow.milk_production.milk_production_history = []
+    cow.milk_production.current_lactation_305_day_milk_produced = current_lactation_305_day_milk_produced
+    cow.mature_equivalent_milking_prediction_305_day = 0.0
+    return cow
+
+
+@pytest.mark.parametrize("calves, parity_factor", [(1, 1.25), (2, 1.18), (3, 1.0)])
+def test_update_mature_equivalent_305_days_milk_production_for_partial_lactation(
+    calves: int, parity_factor: float
+) -> None:
+    """Test mature-equivalent 305-day milk prediction for cows before day 305."""
+    cow = _make_cow_for_305_day_milk_prediction(days_in_milk=120, calves=calves)
+    cow.milk_production.calculate_305_day_milk_yield.return_value = 8000.0
+
+    cow.update_mature_equivalent_305_days_milk_production()
+
+    cow.milk_production.calculate_305_day_milk_yield.assert_called_once_with(1.0, 2.0, 3.0, [], 120)
+    assert cow.mature_equivalent_milking_prediction_305_day == pytest.approx(8000.0 * parity_factor)
+
+
+@pytest.mark.parametrize("calves, parity_factor", [(1, 1.25), (2, 1.18), (3, 1.0)])
+def test_update_mature_equivalent_305_days_milk_production_for_completed_lactation(
+    calves: int, parity_factor: float
+) -> None:
+    """Test mature-equivalent 305-day milk prediction for cows at or after day 305."""
+    cow = _make_cow_for_305_day_milk_prediction(
+        days_in_milk=305, calves=calves, current_lactation_305_day_milk_produced=9000.0
+    )
+
+    cow.update_mature_equivalent_305_days_milk_production()
+
+    cow.milk_production.calculate_305_day_milk_yield.assert_not_called()
+    assert cow.mature_equivalent_milking_prediction_305_day == pytest.approx(9000.0 * parity_factor)
+
+
 @pytest.fixture
 def mock_time() -> RufasTime:
     mock_time = MagicMock(spec=RufasTime)

tests/test_biophysical/test_animal/test_herd_manager/pytest_fixtures.py

@@ -527,6 +527,7 @@ def mock_animal(
     ED_days: int = 0,
     breeding_to_preg_time: int = 0,
     daily_milk_produced: float = 0.0,
+    mature_equivalent_milking_prediction_305_day: float | None = None,
     milk_fat_content: float = 0.0,
     milk_protein_content: float = 0.0,
     sold_at_day: int | None = None,
@@ -578,6 +579,11 @@ def mock_animal(
     animal.milk_production.daily_milk_produced = daily_milk_produced
     animal.milk_production.fat_content = milk_fat_content
     animal.milk_production.true_protein_content = milk_protein_content
+    animal.mature_equivalent_milking_prediction_305_day = (
+        daily_milk_produced
+        if mature_equivalent_milking_prediction_305_day is None
+        else mature_equivalent_milking_prediction_305_day
+    )
 
     return animal
 

tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_daily_routines.py

@@ -514,6 +514,7 @@ def _create_sortable_mock_cow(
 
     cow.milk_production = MagicMock()
     cow.milk_production.daily_milk_produced = daily_milk
+    cow.mature_equivalent_milking_prediction_305_day = daily_milk
 
     cow.days_in_milk = days_in_milk
     cow.days_in_pregnancy = days_in_pregnancy

tests/test_biophysical/test_animal/test_herd_manager/test_herd_manager_misc.py

@@ -66,6 +66,21 @@ def test_collect_daily_feed_request(herd_manager: HerdManager) -> None:
     assert result == expected_total_requested_feed
 
 
+def test_update_milk_305_day_yield_predictions(herd_manager: HerdManager) -> None:
+    """Unit test for update_milk_305_day_yield_predictions()."""
+    cows = [
+        mock_animal(AnimalType.LAC_COW, id=1),
+        mock_animal(AnimalType.LAC_COW, id=2),
+        mock_animal(AnimalType.DRY_COW, id=3),
+    ]
+    herd_manager.cows = cows
+
+    herd_manager.update_milk_305_day_yield_predictions()
+
+    for cow in cows:
+        cow.update_mature_equivalent_305_days_milk_production.assert_called_once_with()
+
+
 def test_print_herd_snapshot(herd_manager: HerdManager, mocker: MockerFixture) -> None:
     """Unit test for print_herd_snapshot()"""
     mock_print = mocker.patch("builtins.print")

tests/test_biophysical/test_animal/test_milk_production/test_milk_production.py

@@ -358,9 +358,34 @@ def test_calculate_daily_milk_production(
         (0.3, 0.4, 0.5, quad(MilkProduction.calculate_daily_milk_production, 1, 305, args=(0.3, 0.4, 0.5))[0]),
     ],
 )
-def test_calc_305_day_milk_yield(l_param: float, m_param: float, n_param: float, expected: float) -> None:
-    """Test the calc_305_day_milk_yield method of the MilkProduction class."""
-    assert MilkProduction.calc_305_day_milk_yield(l_param, m_param, n_param) == pytest.approx(expected, rel=1e-6)
+def test_calculate_305_day_milk_yield(l_param: float, m_param: float, n_param: float, expected: float) -> None:
+    """Test the calculate_305_day_milk_yield method of the MilkProduction class."""
+    assert MilkProduction.calculate_305_day_milk_yield(l_param, m_param, n_param, None) == pytest.approx(
+        expected, rel=1e-6
+    )
+
+
+def test_calculate_305_day_milk_yield_uses_available_history() -> None:
+    """Test that partial 305-day milk yield projections include recorded production history."""
+    milk_production_history = [
+        {"simulation_day": 1, "days_in_milk": 1, "milk_production": 100.0, "days_born": 1000},
+        {"simulation_day": 2, "days_in_milk": 2, "milk_production": 100.0, "days_born": 1001},
+        {"simulation_day": 3, "days_in_milk": 0, "milk_production": 0.0, "days_born": 1002},
+        {"simulation_day": 4, "days_in_milk": 1, "milk_production": 10.0, "days_born": 1003},
+        {"simulation_day": 5, "days_in_milk": 2, "milk_production": 20.0, "days_born": 1004},
+    ]
+    l_param = 0.1
+    m_param = 0.2
+    n_param = 0.3
+    expected_projected_yield = quad(
+        MilkProduction.calculate_daily_milk_production, 3, 305, args=(l_param, m_param, n_param)
+    )[0]
+
+    result = MilkProduction.calculate_305_day_milk_yield(
+        l_param, m_param, n_param, milk_production_history, days_in_milk=2
+    )
+
+    assert result == pytest.approx(200.0 + expected_projected_yield, rel=1e-6)
 
 
 @pytest.mark.parametrize(