# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
# AQAL model: simulating the lifetime trajectory of a dairy cow
# Version 2.1.2 - 03/03/2021
# Major modifications from version 2.0:
# - Impact of gestation and lactation on allocation to growth
# - Energetics of gravid uterus
# Working version for health developments with NF
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
import numpy
# --- Import constitutive modelling elements (variables and solver) ---
# If a new element (like a probability or a flow) is added in the cow biology,
# Then it should be declare in cowElements
from AQAL_V2_cowElements import Solver, FlowAlloc, Acquisition, FlowEnergy, \
QtEnergy, Material, FlowMass, Proba, HealthDisturb, TheoreticalStructMass
###################################################
# ----------------------------------------------- #
# --------- THE COW BIOLOGICAL MODEL ------------ #
# ----------------------------------------------- #
###################################################
# An object of the class CowMod corresponds to a single cow model
# It has its own lifetime trajectory depending on its attributes and methods
# Energy elements are in kcal
class CowMod:
# -------------------------------------------------------------
# -- SOME "PLUG-IN" FUNCTIONS to connect with in/out files ----
# -------------------------------------------------------------
# ------ LOG FUNCTION: looking at the cow --------------
# This method is "empty" because everything is defined in the main file
# The log just need to know time and which cow to log
# Events are logged the same way, with log-ev_xxx, where xxx is the event name.
# The idea is to separate the biological part (= AQAL equations) from the simulation aspects (= observations)
@staticmethod
def log(time, cow):
pass
# ------ NUTRITIONAL ENVIRONMENT: reading values in input files --------
# The connexion with the environment is ensured by "empty" method of the Cow Class
# These methods are connected with functions in the MAIN script (reading txt files)
# Quantity of feed available in the environment (kg DM/cow/day)
@staticmethod
def dm_offer(now):
pass
# dm_offer is linked to env scenario file defined in the main script
# Gross energy content of feed (kcal/kg DM)
@staticmethod
def ge_res(time):
pass
# ge_res is linked to env scenario file defined in the main script
# Health pressure
@staticmethod
def health_pressure(time):
pass
# linked to env scenario file defined in the main script
# Proportion of concentrate in feed (not used in this version)
@staticmethod
def co_prop_dmi(now):
return 0.4
# Proportion of NDF in feed (not used in this version)
@staticmethod
def ndf_propdmi_ref(now):
return 0.384
# -----------------------------------------------------
# ---------- ATTRIBUTES of the class CowMod -----------
# -----------------------------------------------------
# These parameters don't change between individuals, they are common to all objects, ie all cows
# They can be considered as species level parameters
# The difference between species and individual parameters is not set in stone and can change
# -------- NUTRITIONAL REFERENCE VALUES defining a standard environment ---------
# --- Standard gross energy content of DM offer
# This value IS NOT used during cow lifetime simulation (ge_res is used, see above)
# NOT USED in this version (used before for estimating the structural mass asymptote)
ge_res_ref = 4424
# --- Proportion of metabolizable energy in gross energy
# This value is used together with ge_res (from external file) to define intake of ME
# So far, it is a fixed value but can be changed (for instance depends on NDF)
me_pctge_ref = 0.578
# --- Standard metabolizable energy content
# This value is only used for estimating structural mass asymptote
me_res_ref = ge_res_ref * me_pctge_ref
# --------- PARAMETERS for BIOLOGICAL FUNCTIONS --------------
# --- Allocation components ---
# Value to split genetic scaling for lactation allocation between AllocPc and Soma (S2)
# The total quantity of priority flowing between S2 and Pc has to be split in order to get allocation to lactation
# at the beginning of lactation (initialisation of AllocPc)
k_lac0 = 0.25
# Parameters defining the flow of priority between soma and future progeny (Hill function)
# Therefore, these parameters define shape of allocation to gestation
h_pf_0 = 0.000001 # 0.000001
# h_pf_1 is defined as an object attribute as its value is updated after 1st parturition
# we assume that for first gestation, more energy should be allocated to gestation
h_pf_2 = 6 # 6
h_pf_3 = 205
# Transfer rate defining the flow of priority between current progeny and survival
# It drives decrease in allocation to lactation
k_flow_pc2s = 0.0014 # V1 = 0.0024 - V0 = 0.0022
# Transfer rate defining the flow of priority between soma and current progeny (mass action law)
# It drives increase in allocation to lactation
k_flow_s2pc = 0.20 # 0.3
# Transfer rate for growth
k_flow_g2S = 0.0024
# Parameters for the depressive effect of allocation to gestation and lactation on allocation to growth
# This effect represent a decrease in allocation to growth when the female is gestating or lactating
k_Pf_on_g2s = 0.1 # was 0.1
k_Pc_on_g2s = 0.01255 # 0.0122 # was k_AllocG_lac
# Parameters of feedback from Growth allocation and Futur progeny allocation on flow from Pc to S
k_G_on_pc2s = 0.01 #0.001 #0.01 #0.05
k_Pf_on_pc2s = 0.1 #0.1
# --- Acquisition components ---
# Parameters defining the shape of acquisition during lactation
# NEW: lactation acquisition dynamic was changed to not start below zero.
# The dynamic is based on a 3 compartments system A-->B-->C
# with B being the lactation acquisition dynamic.
# kd_acq_lac_init represents the initial charge of B
# ToDo : these parameters will probably be in calibration proc
# --> should be object attributes not class attributes
kd_acq_lac_in = 0.05
kd_acq_lac_out = 0.002 # 0.0026 #0.0025 # was 0.0016
kd_acq_lac_init = 0 # set to zero otherwise no mobilisation the first 2 days of lactation
# Parameters defining the shape of acquisition during gestation
# The shape is a Hill function with 4 parameters but the minimum doesn't change (0.00001), so 3 parameters
kd_acq_gest_1 = 1.5 # was 1.5
kd_acq_gest_2 = 6
kd_acq_gest_3 = 190
# --- Energetics parameters (for energy-material conversion) ---
# Metabolic efficiencies, used to convert ME to NE (kcal/kcal)
eff_gain_lab = 0.67 # efficiency of labile mass gain
eff_lac = 0.64 # efficiency of lactation
eff_loss_lab = 0.82 # efficiency of labile mass loss
eff_mnt = 0.63 # efficiency of maintenance
eff_preg = 0.10 # efficiency of pregnancy
# Energy values of material components (kcal of Net energy/kg of material)
ev_gain_lab = 7500 # gain of 1 kg of labile mass
ev_loss_lab = 7500 # loss of 1 kg of labile mass - if need update, 7440
ev_milk = 739.2 # production of 1 kg of milk (INRAE 2018 ref, 0.42 UFl, 4% fat and 3.1% prot)
# Hill function parameters driving changes in maintenance cost with age
# These values lead to natural death around 20 years
h_ev_mnt_0 = 98 #was 93
h_ev_mnt_1 = 800
h_ev_mnt_2 = 5
h_ev_mnt_3 = 15000
# Parameters of the logistic function, linking energy cost for structural mass gain and age (allometry)
# This function for gain aggregates net energy content and efficiency ME/NE so unit = kcal ME/kg
max_cost_struct = 6480
min_cost_struct = 5138
k_cost_struct = 0.0082
# Hill parameters for change in net energy content of gravid uterus (kcal NE/kg)
h_ev_gu_0 = 300
h_ev_gu_1 = 1380
h_ev_gu_2 = 5
h_ev_gu_3 = 200
# --- Parameters used for body composition ---
# Theoretical structural mass is a breed characteristic (NOT USED here)
breed_mass_struct = 410
# Lipid content in the structural mass
lipcont_mass_struct = 0.025
# Proportion of labile mass at birth
k_lab_birth = 0.03
# Parameters defining change of protein content in structural mass through time
k_P_struct_1 = 0.22
k_P_struct_2 = 0.05 # initial content = k1-k2 = 0.17 at birth
k_P_struct_3 = 0.01
# --- Parameters (in days) related to reproduction physiology ---
# Gestation length
len_gest = 282 # was 282
# Oestrus length (Dennis et al., 2018)
len_oest = 21
# Waiting period length: time before 1st insemination
# len_wait = 60 (not used here)
# Postpartum anestrous
len_ppa = 121 # reproscope IVV = 422d with 280d of gestation --> conception @ 142d after calving - 1 cycle 21d
# Number of days before calving to dry off the cow
len_time_dry = 60
# Age at first conception
age_first_conc = 574 # was 455 (Dennis et al., 2018) # 574 for calving at 28months, 800 for calving 36 months
# Use a fixed calving interval if conception is not random and regular
calv_int = 422 # was 370
# Use a vector of age at conception if conception is not random and not regular
# (here 4 parturitions defined, after a default 365d interval is used)
# list_age_conc = [417, 787, 1157, 1527, 1897, 2267, 2637, 3007] # NZ cow
list_age_conc =[412, 902, 1248, 1620] # JB cow nb 1 (last value is just to have full carreer)
# list_age_dry = [977, 1347, 1717, 2087, 2457, 2827, 3197] # NZ cow
list_age_dry = [1061, 1481, 1845, 2190] # JB cow nb 1
# Time limit for being pregnant again (not used here)
# len_notPreg_multi = 140
# -----------------------------------------------------------
# --------- PARAMETERS for HANDLING SIMULATION ASPECTS ------
# -----------------------------------------------------------
# ------ Management parameters -------
# In a herd architecture, these parameters would be defined in a farm context
# To not multiply script when working on 1 cow, there are defined hereafter
# --- Calendar dates used as time limit for drying and culling for reproductive reason
# t=0 is the 01/01/2000
# Time in year where mating starts (relative to 01/01) --> 284, which is 11/10
date_start_mating = 284
# Time in year when mating season closed (here 70d after start, 10 weeks)
date_end_mating = 354
# Time in year when herd is dried-off --> 141 is mid may
date_herd_dry = 141
# --- Indexes for state variables names ------
# Ensure the consistency between name-number in state vector (ie compartments)
# This enables to write the name of the component instead of the index (avoiding errors)
allocG = 0
allocS1 = 1
allocS2 = 2
allocPc = 3
allocPf = 4
mass_struct = 5
mass_labile = 6
mass_uterus = 7
central_pool = 8
mass_uterus_pot = 9
alloc_Pf_pot = 10
# ----------------------------------------------
# ----------- PHYSIOLOGICAL TIMING -------------
# ----------------------------------------------
def age(self, now):
return now - self.birth_date
def days_in_gestation(self, now):
if self.gest_stat:
return now - self.last_conception_date
else:
return float('nan')
def days_in_lactation(self, now):
if self.lac_stat:
return now - self.last_parturition_date
else:
return float('nan')
# -----------------------------------------------------
# ---------- INITIALIZATION of a new cow object ------
# -----------------------------------------------------
# It is where an object of the class CowMod gets its own attributes and initial values of model components
def __init__(self, now, health_resist, recovery, acq_bas_pot, acq_bas_lac, init_G, s2pc_pot, cow_id):
# --- Initiate a new cow with an ID, date of birth, living status, dates for reproduction
self.cow_id = cow_id
self.birth_date = now
self.living = True
self.last_conception_date = float('nan')
self.last_parturition_date = float('nan')
self.next_oestrus_date = self.age_first_conc + now
print("{} {} {} {}".format(now, self.cow_id, "next_oestrus", self.next_oestrus_date))
# --- Creating elements of the single model (variables, flows...) using generic classes from cowElements
self.allocation = FlowAlloc()
self.acquisition = Acquisition()
self.flow_energy = FlowEnergy()
self.flow_mass = FlowMass()
self.mat_comp = Material()
self.qt_energy = QtEnergy()
self.proba = Proba()
# --- Defining genetic-scaling parameters of the cow (defined when creating the object)
# So far, 4 parameters are defined as "genetic-scaling parameters" but could be more (add in _init)
# Initial value of the growth allocation compartment, at birth
# It means the first day of life, 46% of energy intake will go for growth function
self.init_G = init_G
# Transfer rate of priority between allocation to soma and lactation (dimensionless)
self.s2pc_pot = s2pc_pot
# Maximum value (asymptote) for basal acquisition (kg DM/d)
self.acq_bas_pot = acq_bas_pot
# Maximum value (peak) for lactation acquisition (kg DM/d)
self.acq_lac_pot = acq_bas_lac
# --- Acquisition (kg DM/d) ---
# Basal acquisition at birth
self.acq_basal_init = 2 #1.5 # was 2.5
# Parameter of basal acquisition function (logistic equation)
self.k_acq_bas_mat = 0.012 # 0.0095
self.k_acq_bas_coef = (self.acq_bas_pot - self.acq_basal_init) / self.acq_basal_init
# Maximum of acquisition reached during a lactation (multiplier of the shape of acquisition)
self.acq_lac_max = 1
# List of values for increase in acquisition with parities (proportion of max)
# 1 by default
self.acq_lac_mat = [0.75, 0.83, 0.90, 1, 1, 1] #[0.61, 0.76, 0.87, 0.91, 0.975, 1]
# Proportion of the maximum
self.acq_lac_max_pct = 1
# Dry matter intake, depending on acquisition and feed offer
# Not defined as an instance of Acquisition Class in cowElements because it depends on external thing
self.intake = float('nan')
# Energy value to maintain protein
# Effective energy is 1.63 unit EE/kg P (and EE is somewhere between ME and NE)
# Conversion of 1.63 EE in ME gave 2.1 MJ/kg P
# In Emmans, there is an equation to convert EE into ME based on CP in diet. If diet has 11% Digestible CP in OM,
# which is close to 15% CP in DM, the conversion factor is 1.29 and 1.29 * 1.63 = 2.1 of ME
self.ev_mnt_struct_base = 2.3
# --- Activity costs (component of maintenance requirements)
# 1.44 is AFRC reference (kcal of ME / kg BW)
self.activ_cost = 6.5 # 5.75
# --- Basal conception rate (further modified by milk, BCS and energy balance)
self.p_conc_basal = 0.6
# Coefficient used to increase conception for virgin heifer (if > 1)
self.p_conc_heifer = 1
# --- Initial body mass at birth
self.mass_birth = 42.8 # JB cow#1 = 42.8
# --- Gestation allocation
# Maximum value of the Hill equation for gestation allocation (may not be reached)
self.h_pf_1 = 0.35 # was 0.22
# ----- Reproductive status (0: false, 1: true)
self.gest_stat = False
self.lac_stat = False
self.culling_stat = False
# ----- Reproductive counters
self.gest_nb = 0
self.lac_nb = 0
# ------ State vector (ie structure of compartments) containing allocation, mass components and mnt deficit
# Initialization of the state vector (zero vector and then initial values of compartments)
self.state = numpy.zeros((1, 11), dtype=float)
self.state[0, CowMod.allocG] = self.init_G
self.state[0, CowMod.allocS1] = 1 - self.init_G
self.state[0, CowMod.mass_struct] = self.mass_birth * (1 - CowMod.k_lab_birth)
self.state[0, CowMod.mass_labile] = self.mass_birth * CowMod.k_lab_birth
# ----------- Work in progress --------------
self.fetus_req = 0
self.maternal_mass = 0
self.fetus_mass = 0
#self.ge_res = 3721 # was 3721
# ------ NEW: HEALTH COMPONENT -------
# A health problem is the result of 2 processes:
# 1) a health event occurs, 2) this event may lead to a disturbance
# Occurrence of a health event depends on the disease pressure in the environment
# This pressure is defined between 0 and 1 (to be consistent with other env pressures (e.g. predation))
# This pressure is translated into a frequency (high pressure = high frequency of events)
# Here, it is 1/frequency of exposure
# In this version, the value is passed to the cow only at init so it won't vary with time
# TODO: check health pressure dynamically
self.health_press = CowMod.health_pressure(now)
# The frequency of health events is translated into a lag exposure, i.e. nb of days between 2 events
self.lag_expose_max = 40000 # was 400
self.lag_expose = min(1/self.health_press, self.lag_expose_max)
# After an event, next one will occurs lag_expose days after
self.date_next_health_expose = now + self.lag_expose
# Occurrence of an health event can lead to a disturbance depending on animal resistance
# This resistance is a probability of not being affected by the health event
# It is represented as a random drawn in a uniform distribution with the threshold being the resistance
# Resistance is a genetic-scaling parameter defined at object creation
self.healthEvent_thresh = health_resist
# If the random drawn result is positive:
# The health status changes (true=healthy, false=not healthy)
self.health_status = True
# The cow starts being disturbed, with disturbance represented by a PLM-like system
# Health disturbance is a new type of object (see cowElements), containing the characteristics of a disturbance
# That is to say: intensity, collapse, recovery, time of start (parameters of the PLM like approach)
self.p_illness_intensity = 0.6
self.p_illness_collapse = 0.5
self.p_illness_recovery = recovery
# Each time a cow has an health_event realized, an object "health disturbance" is created
# We need a list to store all health disturbances
self.list_health_disturb = list()
# Add a counter of disturbances
self.count_health_disturb = 0
# And finally all the disturbances achieve an degree of illness
self.illness_degree = 1
# -----------------------------------------------------
# ---- ESTIMATION of STRUCTURAL MASS at maturity ------
# -----------------------------------------------------
# A key point in AQAL is that mature BW is not an input, it is the result of acquisition x allocation dynamics
# As a result, it is not possible to know what will be mature weight and degree of maturity (BW(t) / BW mature)
# Mass at maturity is a breed parameter, used to compute maintenance requirement associated to protein mass
# To get an estimation, the functions hereafter are used, based on the following assumption:
# The mature structural BW will be the structural mass of an animal in a reference environment, not reproducing
# It is only defined by allocation to growth and basal acquisition
# As algebraic integration of equations is very complicated, we use numerical estimation
# Import the functions required to estimate structural mass
estimate_mass = TheoreticalStructMass()
self.mass_struct_th = round(estimate_mass.compute_asymptote(0, 3000, CowMod.me_res_ref, self.acq_bas_pot,
self.acq_basal_init, self.k_acq_bas_mat,
self.init_G, CowMod.k_flow_g2S, 0.22, 0.05, 0.01,
self.mass_birth, CowMod.k_lab_birth), 1)
# self.mass_struct_th = CowMod.breed_mass_struct
print("{} {} {}".format("Structural mass estimated:", self.mass_struct_th, " kg "))
# -----------------------------------------------------
# --------------- COMPUTE SPEED VECTOR ----------------
# -----------------------------------------------------
# These equations define flows between compartments
# Speed tells how the system will change between 2 time steps
def speed(self, now, state):
# Dead cow doesn't change (no speed defined)
if not self.living:
return numpy.zeros((1, 11), dtype=float)
self.compute_aux(now, state) # update auxiliary variables before computing speed vector
# Here it is the definition of the differential system (flows between compartments)
# Each variation (= speed) of a state variable (= compartment) is defined by in/out flows
speed_vec = numpy.zeros((1, 11), dtype=float)
speed_vec[0, CowMod.allocG] = - self.allocation.f_prio_g2s + self.allocation.f_prio_s2g
speed_vec[0, CowMod.allocS1] = self.allocation.f_prio_g2s - self.allocation.f_prio_s2pf_pot + \
self.allocation.f_prio_pc2s - self.allocation.f_prio_s2g
speed_vec[0, CowMod.allocS2] = - self.allocation.f_prio_s2pc
speed_vec[0, CowMod.allocPc] = self.allocation.f_prio_s2pc - self.allocation.f_prio_pc2s
speed_vec[0, CowMod.allocPf] = self.allocation.f_prio_s2pf_pot
speed_vec[0, CowMod.mass_struct] = self.flow_mass.f_mass_struct
speed_vec[0, CowMod.mass_labile] = self.flow_mass.f_mass_rep - self.flow_mass.f_mass_mob
speed_vec[0, CowMod.mass_uterus] = self.flow_mass.f_mass_gravid
speed_vec[0, CowMod.central_pool] = self.flow_energy.f_me_intake \
- self.flow_energy.f_me_gravid_uterus \
- self.flow_energy.f_me_growth \
- self.flow_energy.f_me_milk \
- self.flow_energy.f_me_soma
speed_vec[0, CowMod.mass_uterus_pot] = self.flow_mass.f_mass_gravid_pot
speed_vec[0, CowMod.alloc_Pf_pot] = self.allocation.f_prio_s2pf_pot
return speed_vec
# ---------------------------------------------------
# --------- COMPUTE ALLOCATION FLOWS ----------------
# ---------------------------------------------------
def compute_allocation(self, now, state):
# here, we suppose that the cow is alive
# --- Use the class FlowAlloc to create flows objects
self.allocation = FlowAlloc()
# --- Then define equations for flows
# Flow driving allocation to growth, from Growth to Soma
# k_flow_g2s_corr = CowMod.k_flow_g2S * (1.45 + (1 - 1.45)*((self.age(now)**2)/(self.age(now)**2 + 75**2)))
k_flow_g2s_corr = CowMod.k_flow_g2S
self.allocation.f_prio_g2s = state[0, CowMod.allocG] * (
k_flow_g2s_corr + CowMod.k_Pf_on_g2s * state[0, CowMod.allocPf] )
self.allocation.f_prio_s2g = state[0, CowMod.allocG] * (CowMod.k_Pc_on_g2s * state[0, CowMod.allocPc])
# Flow driving allocation to gestation, from Soma to Future Progeny
if self.gest_stat:
gest_time = self.days_in_gestation(now)
self.allocation.f_prio_s2pf = (self.h_pf_1 - CowMod.h_pf_0) * \
(numpy.exp(-CowMod.h_pf_2 * numpy.log(CowMod.h_pf_3)) *
CowMod.h_pf_2 * gest_time ** (CowMod.h_pf_2 - 1)) / \
(1 + (numpy.exp(-CowMod.h_pf_2 *
numpy.log(CowMod.h_pf_3)) * gest_time ** CowMod.h_pf_2)) ** 2
self.allocation.f_prio_s2pf_pot = (self.flow_energy.f_me_gravid_uterus_pot / \
self.flow_energy.f_me_intake) - state[0, CowMod.allocPf]
else:
self.allocation.f_prio_s2pf = 0.0
self.allocation.f_prio_s2pf_pot = 0
# Flows driving allocation to lactation, between Soma and Current Progeny
if self.lac_stat:
if self.gest_stat:
self.allocation.f_prio_pc2s = state[0, CowMod.allocPc] * \
(CowMod.k_flow_pc2s + state[0, CowMod.allocPf] * CowMod.k_Pf_on_pc2s +
state[0, CowMod.allocG] * CowMod.k_G_on_pc2s)
else:
self.allocation.f_prio_pc2s = state[0, CowMod.allocPc] * (CowMod.k_flow_pc2s +
state[0, CowMod.allocG] * CowMod.k_G_on_pc2s) # was 0.02
self.allocation.f_prio_s2pc = CowMod.k_flow_s2pc * state[0, CowMod.allocS2]
else:
self.allocation.f_prio_pc2s = 0.0
self.allocation.f_prio_s2pc = 0.0
# ---------------------------------------------------
# --------- COMPUTE ACQUISITION VARIABLES -----------
# ---------------------------------------------------
def compute_acquisition(self, now, state):
# Here, we suppose that the cow is alive
# --- Use the class Acquisition to create variables
self.acquisition = Acquisition()
# --- Define equations for acquisition
# Basal acquisition
self.acquisition.acq_bas = self.acq_bas_pot / (1 + self.k_acq_bas_coef *
numpy.exp(-self.k_acq_bas_mat * self.age(now)))
# Lactation acquisition
if self.lac_stat:
lac_time = self.days_in_lactation(now)
# First, definition of the shape of lactation acquisition
# Which is based on the algebraic form of the intermediary compartment of a 3-compartments network
coef_acq_lac = (numpy.log(CowMod.kd_acq_lac_in) - numpy.log(CowMod.kd_acq_lac_out)) /\
(CowMod.kd_acq_lac_in - CowMod.kd_acq_lac_out)
max_acq_lac_dyn = (1 - CowMod.kd_acq_lac_init) * (CowMod.kd_acq_lac_in / (CowMod.kd_acq_lac_out -
CowMod.kd_acq_lac_in)) * \
(numpy.exp(-CowMod.kd_acq_lac_in * coef_acq_lac) -
numpy.exp(-CowMod.kd_acq_lac_out * coef_acq_lac)) +\
CowMod.kd_acq_lac_init
self.acquisition.acq_lac_dyn = (((1 - CowMod.kd_acq_lac_init) * (CowMod.kd_acq_lac_in / (
CowMod.kd_acq_lac_out - CowMod.kd_acq_lac_in)) * (
numpy.exp(-lac_time * CowMod.kd_acq_lac_in) -
numpy.exp(- lac_time * CowMod.kd_acq_lac_out))) +
CowMod.kd_acq_lac_init)
# Then the shape of lactation acquisition is scaled with acq_lac_max
# This latter reflects the genetic scaling parameter + age effect (value updated at parturition)
self.acquisition.acq_lac = self.acquisition.acq_lac_dyn * (self.acq_lac_max / max_acq_lac_dyn)
else:
self.acquisition.acq_lac_dyn = 0.0
self.acquisition.acq_lac = 0.0
if self.gest_stat and not self.lac_stat:
# Acquisition has a gestation part, ie the cow is able to ingest a bit more during gestation
# This is represented with a Hill function
# The maximum DMI during gestation (Hill parameter 1) is higher for multiparous than primiparous cows
if self.gest_nb == 1:
max_acq_gest = CowMod.kd_acq_gest_1
else:
max_acq_gest = CowMod.kd_acq_gest_1 * 2 # was 1.25
self.acquisition.acq_gest = 0.00001 + (max_acq_gest - 0.00001) * \
((self.days_in_gestation(now)**CowMod.kd_acq_gest_2) /
((self.days_in_gestation(now)**CowMod.kd_acq_gest_2) +
CowMod.kd_acq_gest_3**CowMod.kd_acq_gest_2))
else:
self.acquisition.acq_gest = 0
# Total acquisition
# there is nothing that prevents acq_lac to make acq_tot negative
# this is why there is a max function, to prevent acq_tot being negative
# NEW: degree of illness impacts acquisition
self.acquisition.acq_tot = max(0, self.acquisition.acq_bas + self.acquisition.acq_lac +
self.acquisition.acq_gest) * self.illness_degree
# --- And finally defined DM intake, maximum between environment offer and genetic capacity (total acq)
self.intake = min(self.acquisition.acq_tot, CowMod.dm_offer(now))
# dm_offer will be replaced by possible intake
# --- Conversion of GE into ME depends on DM metabolisability = digestibility and loss in urine and CH4
# Metabolisability is affected by the level of intake and NDF
self.acquisition.ME_in_GE = (50.9 - 1.01 * 100 *
((self.intake/self.state[0, CowMod.mass_struct])/1.59) + 7.9)/100
#print(self.acquisition.ME_in_GE)
#self.acquisition.ME_in_GE = CowMod.me_pctge_ref
# ---------------------------------------------------
# --------- COMPUTE Metabolizable ENERGY FLOWS ------
# ---------------------------------------------------
# These flows are the result of acquisition (ME entering the system) and allocation (proportion for each function)
def compute_flowEnergy(self, now, state): # Requires acquisition
# here, we suppose that the cow is alive
# if self.lac_stat and self.lac_nb > 1:
# self.ge_res = 4750
# if self.lac_stat and self.lac_nb == 1:
# self.ge_res = 4400
# if not self.lac_stat and self.lac_nb >= 1:
# self.ge_res = 4400
self.flow_energy = FlowEnergy()
self.flow_energy.f_me_intake = self.intake * CowMod.ge_res(now) * self.acquisition.ME_in_GE
if self.age(now) < self.list_age_conc[0] + CowMod.len_gest:
self.flow_energy.f_me_intake = self.flow_energy.f_me_intake * self.MEreduction
self.flow_energy.f_me_growth = self.flow_energy.f_me_intake * state[0, CowMod.allocG]
self.flow_energy.f_me_gravid_uterus = self.flow_energy.f_me_intake * state[0, CowMod.allocPf]
self.flow_energy.f_me_milk = self.flow_energy.f_me_intake * state[0, CowMod.allocPc]
self.flow_energy.f_me_soma = self.flow_energy.f_me_intake * (state[0, CowMod.allocS1] + state[0, CowMod.allocS2])
if self.days_in_gestation(now) > 0:
self.flow_energy.f_me_gravid_uterus_pot = (self.flow_mass.f_mass_gravid_pot *
self.qt_energy.ev_gain_uterus) / CowMod.eff_preg
else:
self.flow_energy.f_me_gravid_uterus_pot = 0
# ----------------------------------------------------------------
# ---- COMPUTE MATERIAL Elements (body components and milk) ------
# ----------------------------------------------------------------
def compute_material(self, now, state): # Requires acquisition, flowEnergy
# here, we suppose that the cow is alive
self.mat_comp = Material()
self.mat_comp.empty_bw = state[0, CowMod.mass_struct] + state[0, CowMod.mass_labile]
# Digestive contents are a fixed part of BW (coef = 16.3%)
# with DC = coef*(EBW+DC) <-> DC = coef/(1-coef)*EBW
self.mat_comp.mass_digcont = (0.163 / (1 - 0.163)) * self.mat_comp.empty_bw
self.mat_comp.mass_tot = state[0, CowMod.mass_struct] + state[0, CowMod.mass_labile] + \
state[0, CowMod.mass_uterus] + self.mat_comp.mass_digcont
self.mat_comp.milk_std = (self.flow_energy.f_me_milk * CowMod.eff_lac) / CowMod.ev_milk
# Conversion of standard milk into raw milk based NZ data
if self.lac_stat:
if self.mat_comp.milk_std > 2:
# New Zealand equation (K. MacDonald data)
#self.mat_comp.milk_raw = self.mat_comp.milk_std * 1.0149 - 1.8979
# France equation (Defilait data)
self.mat_comp.milk_raw = self.mat_comp.milk_std * 1.044 - 0.2065
else:
self.mat_comp.milk_raw = self.mat_comp.milk_std
# condition on milk std to avoid negative milk_raw
else:
self.mat_comp.milk_raw = 0
# Empty body fat estimation (kg lipids)
# Estimation 1 is based on Martin and Sauvant 2010 (limit: degree maturity estimation)
self.mat_comp.mass_lip_1 = state[0, CowMod.mass_labile] * 0.75 + \
state[0, CowMod.mass_struct] * CowMod.lipcont_mass_struct * \
(state[0, CowMod.mass_struct] / self.mass_struct_th)
# Estimation 2 is based on fat in body reserves (made of 75% fat) (remove: "representing 90% of total body fat")
self.mat_comp.mass_lip_2 = (state[0, CowMod.mass_labile] * 0.75) # 0.9 division removed
# BCS estimation n°1, based on Martin et Sauvant, 2010
self.mat_comp.bcs_1 = ((self.mat_comp.mass_lip_1 / ((1 - numpy.exp(-0.0068 * (
self.mat_comp.empty_bw - self.mat_comp.mass_lip_1))) ** 3.7235)) - 17.9) / 39.8
# BCS estimation n°2 based on prop_fat_ebw = slope x ( BCS - intercept )
# Which is equivalent to BCS = 1/slope * prop_fat_ebw + intercept
self.mat_comp.slope = 0.11
self.mat_comp.intercept = 0.35
self.mat_comp.bcs_2 = (self.mat_comp.mass_lip_2/self.mat_comp.empty_bw) * (1/self.mat_comp.slope) +\
self.mat_comp.intercept
# BCS estimation n°3, testing the idea of proportion of mature weight in one unit of CS (from Zygoyiannis 1997)
unit_CS = ((self.mass_struct_th / 0.7) / 0.84) * 0.13
self.mat_comp.bcs_3 = 2.5 + ((state[0, CowMod.mass_labile] - 0.3 * self.mass_struct_th) / unit_CS)
# -------------------------------------------
# ---------- COMPUTE PROBABILITIES ----------
# -------------------------------------------
def compute_proba(self, now, state): # Requires Material, QtEnergy
# here, we suppose that the cow is alive
self.proba = Proba()
self.proba.p_conc_bcs = min(1, 2 / (1 + numpy.exp(-0.425 * (self.mat_comp.bcs_1 - 3))))
self.proba.p_conc_milk = min(1, 2 / (1 + numpy.exp(0.012 * self.mat_comp.milk_std)))
self.proba.p_conc_eb = min(1, 2 / (1 + numpy.exp(-0.027 * - self.qt_energy.me_lab_mob / 239)))
# Probability of conception is computed differently for virgin heifer versus cow
# For cow (ie when gestation nb > 1), the probability depends on basal value (which is a max) corrected for:
# BCS, milk production and energy balance (from Phuong et al. 2015)
if self.gest_nb > 0:
self.proba.p_conc = min(self.p_conc_basal,
self.p_conc_basal * self.proba.p_conc_bcs * self.proba.p_conc_milk * self.proba.p_conc_eb)
# For heifer, the probability is the basal value (that can be increase with a coefficient)
# For example, in NZ, virgin heifer are hyper fertile (10% more) so p_conc_heifer = 1.1
else:
self.proba.p_conc = self.p_conc_basal * self.p_conc_heifer
# Probability of death depends on extreme values of body reserves proportion
self.proba.p_death = 0
if state[0, CowMod.mass_labile] / self.mat_comp.empty_bw < 0.01:
self.proba.p_death = 1
if state[0, CowMod.mass_labile] / self.mat_comp.empty_bw > 0.6:
self.proba.p_death = 1
# -----------------------------------------------
# ------- COMPUTE QUANTITIES of ENERGY ----------
# -----------------------------------------------
def compute_qtEnergy(self, now, state): # Requires flowEnergy, material
self.qt_energy = QtEnergy()
# ----- COST of GROWTH for structural mass (kcal NE/kg structural mass)
# Compute proportion of protein in gain of 1 kg mass structural (changing with age)
self.qt_energy.P_in_struct = CowMod.k_P_struct_1 - CowMod.k_P_struct_2 * numpy.exp(
-CowMod.k_P_struct_3 * self.age(now))
# Compute NET energy content of 1 kg structural mass (with 2.5% fat)
# 5.7 Mcal/kg of protein and 9.5 Mcal/hg lipids (equivalent to 23.85 MJ and 39.75 MJ)
self.qt_energy.ev_gain_struct = self.qt_energy.P_in_struct * 5700 + CowMod.lipcont_mass_struct * 9500
# Efficiency ME/NE is computed in flow mass part
# Old Equation
# self.qt_energy.ev_gain_struct = (CowMod.max_cost_struct / (
# 1 + ((CowMod.max_cost_struct - CowMod.min_cost_struct) / CowMod.min_cost_struct) *
# numpy.exp(- CowMod.k_cost_struct * self.age(now))))
# ------ COST of GRAVID UTERUS (depending on gestation stage) in kcal NE/kg gravid uterus
self.qt_energy.ev_gain_uterus = CowMod.h_ev_gu_0 + (CowMod.h_ev_gu_1 - CowMod.h_ev_gu_0) * (
(self.days_in_gestation(now) ** CowMod.h_ev_gu_2) / (CowMod.h_ev_gu_3 ** CowMod.h_ev_gu_2 +
self.days_in_gestation(now) ** CowMod.h_ev_gu_2))
# ------ COST of MAINTENANCE (kcal ME/kg of BW)
# This cost depends on age to simulate senescence
# When the cow is old, the cost becomes high leading to no more reserves and death probability
# self.qt_energy.ev_mnt = CowMod.h_ev_mnt_0 + (CowMod.h_ev_mnt_1 - CowMod.h_ev_mnt_0) * (
# self.age(now) ** CowMod.h_ev_mnt_2) / (
# CowMod.h_ev_mnt_3 ** CowMod.h_ev_mnt_2 + self.age(now) ** CowMod.h_ev_mnt_2)
# Fixed cost of maintenance, in net energy
self.qt_energy.ev_mnt = 93 # equivalent to 143 kcal of ME with 0.65 efficiency ME/NE
# Maternal mass with the assumption that fetus mass represents 0.597 of gravid uterus
self.maternal_mass = self.state[0, CowMod.mass_struct] + self.state[0, CowMod.mass_labile] + \
self.mat_comp.mass_digcont + self.state[0, CowMod.mass_uterus] * (1 - 0.597)
self.fetus_mass = self.state[0, CowMod.mass_uterus] * 0.597
# --------------- WORK IN PROGRESS: NOT USED, JUST COMPUTED ----------------------
# Assumption: fetus requirements are twice maternal ones
# Need to add a condition on gestation stage
# With only gest_stat = TRUE, Python finds an invalid scalar (?)
if self.days_in_gestation(now) > 30:
# DEBUG : print("{} {} ".format(now, self.state[0, CowMod.mass_uterus]))
self.fetus_req = (self.fetus_mass ** 0.75) * self.qt_energy.ev_mnt * 2
else:
self.fetus_req = 0
# ----------------------------------------------------------------------------------
# ----- Compute maintenance requirements based on maternal mass
self.qt_energy.me_mnt_req_1 = (self.activ_cost * self.mat_comp.mass_tot) + (
(self.mat_comp.mass_tot ** 0.75) * self.qt_energy.ev_mnt) / CowMod.eff_mnt
# ---- Other option: compute maintenance requirement based on structural mass and protein
# this new maintenance requirement is directly in ME so no need to convert with an efficiency coefficient
# Further activity is added with an additional cost, not by multiplying maintenance requirement
# The value 1.0277 corresponds to conversion of AFRC value (0.0043 MJ/kg BW) into kcal
self.qt_energy.Pmat = self.mass_struct_th * (1 - CowMod.lipcont_mass_struct) * 0.2224
self.qt_energy.P = self.state[0, CowMod.mass_struct] * (1 - CowMod.lipcont_mass_struct) * \
self.qt_energy.P_in_struct
self.qt_energy.ev_mnt_struct_age = 2 * numpy.exp(-0.004 * self.age(now)) + self.ev_mnt_struct_base
self.qt_energy.me_mnt_Pstruct = (self.qt_energy.P * ((self.qt_energy.ev_mnt_struct_age /
(self.qt_energy.Pmat ** 0.27)) / 4.184 * 1000))
self.qt_energy.me_mnt_activity = self.activ_cost * self.mat_comp.mass_tot
self.qt_energy.me_mnt_req_2 = self.qt_energy.me_mnt_Pstruct + self.qt_energy.me_mnt_activity
# ---- Option retained to calculate maintenance
self.qt_energy.me_mnt_req = self.qt_energy.me_mnt_req_2
# ---- ENERGY BALANCE to see if the cow is mobilizing or reconstituting body reserves,
# based on the difference between maintenance requirements and allocation to survival
self.qt_energy.me_var_stock = self.flow_energy.f_me_soma - self.qt_energy.me_mnt_req
# Update the boolean variable that indicates if energy allocated is enough or not to cover maintenance
if self.qt_energy.me_var_stock >= 0:
me_var_stat = 1
else:
me_var_stat = 0
# Based on this status, compute energy that should be mobilised (max = size of labile mass)
self.qt_energy.me_lab_mob = min(((1 - me_var_stat) * (-self.qt_energy.me_var_stock)),
(state[0, CowMod.mass_labile] * CowMod.ev_loss_lab * CowMod.eff_loss_lab))
# Compute if energy can be stored
self.qt_energy.me_lab_rep = me_var_stat * self.qt_energy.me_var_stock
# Based on the balance requirements/reserves, compute energy for maintenance (equal to requirement or below)
self.qt_energy.me_mnt = self.flow_energy.f_me_soma - self.qt_energy.me_lab_rep + self.qt_energy.me_lab_mob
# -----------------------------------------------
# ------- COMPUTE FLOWS of MASS COMPONENTS ------
# -----------------------------------------------
def compute_flowMass(self, now, state): # Requires flowEnergy, qt_energy
# Here, we suppose that the cow is alive
self.flow_mass = FlowMass()
# ---- Mass flow to structural mass is equal to ME allocated divided by ME cost of 1 kg
# Net cost is computed in qt_energy
# Efficiency is based on partial efficiencies of protein and lipid (Williams and Jenkins 2003)
eff_gain_struct = 0.75/(1 + 2.75 * ((self.qt_energy.P_in_struct * 5700) / self.qt_energy.ev_gain_struct))
self.flow_mass.f_mass_struct = (self.flow_energy.f_me_growth * eff_gain_struct) / self.qt_energy.ev_gain_struct
# ---- Mass flow to gravid uterus is equal to ME allocated, multiplied by efficiency ME-->NE
# and divided by NE cost of 1 kg
if self.gest_stat:
self.flow_mass.f_mass_gravid = self.flow_energy.f_me_gravid_uterus * (
CowMod.eff_preg / self.qt_energy.ev_gain_uterus)
self.flow_mass.f_mass_gravid_pot = (2.31347143 * self.mass_birth * 0.024 * \
numpy.exp(0.024*(240+self.days_in_gestation(now)))) / \
(numpy.exp(0.024 * 240) +
numpy.exp(0.024 * self.days_in_gestation(now)))**2
else:
self.flow_mass.f_mass_gravid = 0.0
self.flow_mass.f_mass_gravid_pot = 0
# ---- Mass flow mobilized from labile mass is equal to ME allocated, multiplied by 1/efficiency ME/NE
# and divided by NE cost of 1 kg
self.flow_mass.f_mass_mob = self.qt_energy.me_lab_mob * (1 / CowMod.eff_loss_lab) * (1 / CowMod.ev_loss_lab)
# ---- Mass flow toward labile mass is equal to ME allocated, multiplied by efficiency ME/NE
# and divided by NE cost of 1 kg
self.flow_mass.f_mass_rep = self.qt_energy.me_lab_rep * CowMod.eff_gain_lab * (1 / CowMod.ev_gain_lab)
# ------------------------------------------
# ------- COMPUTE HEALTH DISTURBANCES ------
# ------------------------------------------
# ------ Compute dynamics of disturbances and total degree of illness
def compute_health_dist(self, now, state):
# create an empty list, to stock values of each disturbance
temp_stock = list()
# Loop to calculate value of each disturbance
for dist in self.list_health_disturb:
coef_temp = (dist.p_intensity * dist.p_collapse) / (dist.p_collapse - dist.p_recovery)
time_dist = now - dist.p_start
# The equation is the algebraic form of the 2d compartment in a 3 compartments system (PLM like)
dist_degree = 1 - (coef_temp * (numpy.exp(-dist.p_recovery * time_dist)
- numpy.exp(-dist.p_collapse * time_dist)))
temp_stock.append(dist_degree)
# The total quantity of disturbance is the product of each disturbance
self.illness_degree = numpy.prod(temp_stock)
# -------------------------------------------------------------------------
# -------- EXECUTE computation of auxiliary variables and flows -----------
# -------------------------------------------------------------------------
# This function will do the calculation of the different model elements, defined before
# Starting with allocation and acquisition to determine energy for each function,
# Then these energies are converted into material (BW components, milk) and probabilities
def compute_aux(self, now, state):
if not self.living:
return # warning, when cow is dead, values are undefined (so old values may be here)
# -- NEW: compute disturbances
self.compute_health_dist(now, state) # Need to be before acquisition and allocation
self.compute_allocation(now, state)
self.compute_acquisition(now, state)
self.compute_flowEnergy(now, state) # require acquisition
self.compute_material(now, state) # require acquisition, flow energy
self.compute_qtEnergy(now, state) # require material, flow_energy
self.compute_flowMass(now, state) # require flow_energy, qt_energy
self.compute_proba(now, state) # require material, qt_energy
# -----------------------------------------
# ------- DISCRETE EVENTS definition ------
# -----------------------------------------
# The logic in AQAL is to have the biological energetic processes (continuous) that are modified by some events,
# Reflecting time point at which something happen (like at parturition, gravid uterus is suddenly emptied)
# ----- Actions of DISCRETE EVENTS, related to reproduction and management -------
# ------- CONCEPTION -----------
@staticmethod
def evlog_conception(time, cow, msg):
pass
def event_conception(self, now, msg=""):
if not self.living:
return
CowMod.evlog_conception(now, self, msg) # ask Pierre about order of args in function
self.last_conception_date = now
self.next_oestrus_date = float('nan') # this enables to deactivate a potential conception event
self.gest_stat = True
self.gest_nb += 1
# modification: value of allocation just before parturition should be higher in 1st gestation
# This way heifer and multi have the same size of gravid uterus
if self.gest_nb == 1:
# update of the maximum value of allocation to gestation
self.h_pf_1 = 0.35
if self.gest_nb > 1:
self.h_pf_1 = 0.30
# ------- PARTURITION -------------
@staticmethod
def evlog_parturition(time, cow, msg):
pass
def event_parturition(self, now, msg=""):
if not self.living:
return
# NEW: empty uterus before log
# Emptying of uterus = mass equal zero
self.state[0, CowMod.mass_uterus] = 0
self.state[0, CowMod.mass_uterus_pot] = 0
CowMod.evlog_parturition(now, self, msg)
# update dates
self.last_parturition_date = now
# update physiological status
self.gest_stat = False
self.lac_stat = True
self.lac_nb += 1
# the next oestrus is set from calving + waiting period (anoestrus post-partum) + the first oestrus cycle
self.next_oestrus_date = now + self.len_ppa + self.len_oest
# Initialization of compartments involved in allocation to lactation
# S2 is the part of survival priority that will be invested in lactation.
# The charge of S2 conditions how much priority will flow through lactation compartment
# First of all, no more allocation to gestation
self.state[0, CowMod.allocPf] = 0
# Charge of alloc Pc, as a proportion of genetic potential:
self.state[0, CowMod.allocPc] = CowMod.k_lac0 * self.s2pc_pot
# Charge of allocS2: this initial charge defines how much priority will flow from Soma to lactation
# The initial charge is then adjusted for growth allocation
self.state[0, CowMod.allocS2] = (1 - CowMod.k_lac0) * self.s2pc_pot - self.state[0, CowMod.allocG]
# allocS1 is adjusted to ensure that the sum of state variables are equal to 1
self.state[0, CowMod.allocS1] = 1 - self.state[0, CowMod.allocPc] - self.state[0, CowMod.allocG] - \
self.state[0, CowMod.allocS2] - self.state[0, CowMod.allocPf]
# update parameters for lactation acquisition, depending on parity (scaling of lactation acquisition)
if self.gest_nb > len(self.acq_lac_mat):
self.acq_lac_max_pct = 1
else:
self.acq_lac_max_pct = min(1, self.acq_lac_mat[self.gest_nb - 1])
self.acq_lac_max = self.acq_lac_pot * self.acq_lac_max_pct
# --------- DRYING OFF --------------
@staticmethod
def evlog_dryoff(time, cow, msg):
pass
def event_dryoff(self, now, msg=""):
if not self.living:
return
CowMod.evlog_dryoff(now, self, msg)
self.lac_stat = False
self.state[0, CowMod.allocS1] = self.state[0, CowMod.allocS1] + self.state[0, CowMod.allocPc] + self.state[
0, CowMod.allocS2]
self.state[0, CowMod.allocPc] = 0
self.state[0, CowMod.allocS2] = 0
if self.culling_stat:
self.event_die(now, "culled for repro reason")
# --------- CULLING for reproductive reason -------
@staticmethod
def evlog_culling_repro(time, cow, msg=""):
pass
def event_culling_repro(self, now, msg=""):
if not self.living:
return
CowMod.evlog_culling_repro(now, self, msg="")
# if the cow is lactating and pointed out for culling, wait until dryoff
if self.lac_stat:
self.culling_stat = True
# if the cow is not lactating, remove her
else:
self.event_die(now, "culled for repro reason")
# ----------- DEATH -------------
@staticmethod
def evlog_die(time, cow, msg):
pass
def event_die(self, now, msg): # warning, here msg is required, and must describe the cause of death
if not self.living:
return
CowMod.evlog_die(now, self, msg)
print("{} {} {} {} {}".format(now, self.cow_id, "dead at", self.age(now), msg))
self.living = False
# -------- NEW: HEALTH event ----------
def new_health_event(self, now):
# To test the code, the next line is commented as the event trigger occurs at fixed time
# If random is needed, remove the # and add 4 spaces on each line after the "if"
if numpy.random.uniform(0,1) >= self.healthEvent_thresh:
self.health_status = False
self.count_health_disturb += 1
print("{} {}".format("number of perturb: ", self.count_health_disturb))
# create a new disturbance
new_dist = HealthDisturb(self.p_illness_intensity, self.p_illness_collapse, self.p_illness_recovery, now)
# add the disturbance to the cow's list of health disturbances
self.list_health_disturb.append(new_dist)
print("{} {}".format("Not healthy at: ", now))
# ---------- Triggers of DISCRETE EVENTS --------------
def compute_events(self, now):
# Events are only compute if cow is alive
if not self.living:
return
# Hereafter, 3 options of simulating conception:
# 1) with no alea + fixed calving interval
# 2) with no alea + used defined interval
# 3) with alea (probability defined by biological variables
# Add or remove # to comment or not the block and activate one or the other option
# ------- CONCEPTION with NO ALEA --------
# --- Option 1 : fixed calving interval
#if (self.age(now) - CowMod.age_first_conc) % CowMod.calv_int <= 0.1 and self.age(now) >= \
# CowMod.age_first_conc and not (self.gest_stat):
# print("{} {} {} {}".format(now, self.cow_id, "fixed conception at", self.age(now)))
# self.event_conception(now)
# --- Option 2: user defined calving intervals
if self.gest_nb + 1 <= len(self.list_age_conc):
if self.age(now) == self.list_age_conc[self.gest_nb] and not self.gest_stat:
print("{} {} {} {}".format(now, self.cow_id, "fixed conception (user) at", self.age(now)))
self.event_conception(now)
#else:
# if (self.age(now) - CowMod.age_first_conc) % CowMod.calv_int <= 0.1 and self.age(now) >= \
#CowMod.age_first_conc and not (self.gest_stat):
# print("{} {} {} {}".format(now, self.cow_id, "fixed conception (interval) at", self.age(now)))
# self.event_conception(now)
# --- Option 3: random conception based on cow probability
# Is the cow ready for conception --> not gestating and at oestrus date ?
# # --> YES
# if not self.gest_stat and self.next_oestrus_date == now:
# print("{} {} {} {}".format(now, self.cow_id, "oestrus at", self.age(now)))
# # #
# # # # Is the mating season opened ?
# # # # Compute time on a yearly basis to further check
# # # # Maiden heifer is mated 10d before cows, so +10 on year time to take that into account
# if self.lac_nb == 0:
# year_time = (now - 365 * int(now / 365)) + 10
# else:
# year_time = (now - 365 * int(now / 365))
# # # # --> YES mating opened: opportunity for conception
# if CowMod.date_start_mating <= year_time <= CowMod.date_end_mating:
# print("{} {} {}".format(now, year_time, "mating season"))
# # # # Is the fecondation a success? ie the result of the random draw + mating season opened
# # # # --> YES
# if numpy.random.uniform(0, 1) <= self.proba.p_conc:
# print("{} {} {} {}".format(now, self.cow_id, "random conception at", self.age(now)))
# self.event_conception(now)
# # # # --> NO
# else:
# self.next_oestrus_date = now + self.len_oest
# # # # --> NO mating closed: just compute next oestrus date
# else:
# self.next_oestrus_date = now + self.len_oest
# end of the condition "if not self.gest_stat and self.next_oestrus_date == now"
# --------- PARTURITION -----------
if self.gest_stat and (self.days_in_gestation(now) - CowMod.len_gest) >= 0.001:
print("{} {} {} {}".format(now, self.cow_id, "parturition", self.age(now)))
self.event_parturition(now)
# --------- DRYING OFF -----------
# --- First rule for drying off is a cow rule
# --> drying off "len_time_dry" days before calving (90d for NZ)
#if self.lac_stat and self.days_in_gestation(now) >= (CowMod.len_gest - CowMod.len_time_dry):
# print("{} {} {} {}".format(now, self.cow_id, "individual drying-off", self.age(now)))
# self.event_dryoff(now)
# --- Second rule for drying off is a calendar rule
# --> drying off at date_herd_dry (farmer time limit for milking)
# if self.lac_stat and (now - 365 * int(now / 365)) == CowMod.date_herd_dry :
# print("{} {} {} {}".format(now, self.cow_id, "herd drying-off", self.age(now)))
# self.event_dryoff(now)
# --- Third rule is a user defined date for drying off (useful for calibration with real data)
if self.lac_stat and len(self.list_age_dry) >= self.lac_nb and \
(self.age(now) - self.list_age_dry[self.lac_nb-1]) >= 0.001:
print("{} {} {} {}".format(now, self.cow_id, "user defined drying-off", self.age(now)))
self.event_dryoff(now)
#if self.lac_stat and len(self.list_age_dry) < self.lac_nb \
# and self.days_in_gestation(now) >= (CowMod.len_gest - CowMod.len_time_dry):
#print("{} {} {} {}".format(now, self.cow_id, "individual drying-off", self.age(now)))
#self.event_dryoff(now)
# --------- CULLING FOR REPRODUCTIVE REASON --------
# Rule for a culling depending on cow only (time after calving when not pregnant)
# if self.lac_nb > 1 and not self.gest_stat and now - self.last_parturition_date > CowMod.len_notPreg_multi\
# and not self.culling_stat:
# print("{} {} {} {}".format(now, self.cow_id, "detected for culling", self.age(now)))
# self.event_culling_repro(now)
# Rule for a culling that depends on a calendar timing: at the end of mating, non pregnant cows culled
# if not self.gest_stat and now - 365 * int(now/365) == CowMod.date_end_mating and not self.culling_stat\
# and self.age(now) > self.age_first_conc:
# print("{} {} {} {}".format(now, self.cow_id, "detected for culling", self.age(now)))
# self.event_culling_repro(now)
# --------- DEATH ----------
if numpy.random.uniform(0, 1) <= self.proba.p_death:
print("{} {} {} {}".format(now, self.cow_id, "natural death", self.age(now)))
self.event_die(now, "natural death")
# --------NEW: HEALTH EVENT --------
if now >= self.date_next_health_expose:
#to test the code, the time of occurrence can be fixed
#if now == 1180 or now == 1300:
print("{} {} ".format(now, "health event"))
#self.date_last_health_expose = now
self.date_next_health_expose = now + self.lag_expose
self.new_health_event(now)
# end def compute_events(self, now)
# ----------------------------------------------------
# ------- SIMULATION OF COW LIFETIME TRAJECTORY ------
# ----------------------------------------------------
# this function corresponds to a loop executing all previous functions at each time step
def simul_life(self, begin_time, end_time):
time = begin_time
dt = 1.0 # time step for integration in days
continue_loop = True
solver = Solver(dt) # instantiate an object from the solver class
while time < end_time and continue_loop: # check if time step is still within the range of simulation
if end_time - time > dt:
step = dt
else:
step = end_time - time # adjust the last integration step to finish at end_time
continue_loop = False
# Compute aux, events and then log
self.compute_aux(time, self.state)
self.compute_events(time)
CowMod.log(time, self)
# Compute the new state of the system (used at the beginning of next time step)
self.state = solver.solve(self.state, time, (time + step), self.speed)
# next time step
time += step