# FEMALE FORWARD ITERATION, FOLLOWING A COHORT. # Parameter definitions and values. # offspring_loss = probability of foetus/infant loss given survival of mother at period d since conception. offspring_loss = (0.02, 0.02, 0.02, 0.02, 0.04, 0.04, 0.04, 0.04, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01) # offspring_energy = energetic requirements in kcal of foetus/infant at period d since conception. offspring_energy = (4500.0, 4500.0, 4500.0, 4500.0, 4500.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0) # max_time & min_time are the time points between which simulations are run. max_time = 250 min_time = 0 # min_size & max_size are the minimum and maximum sizes an adult female can attain. min_size = 8 max_size = 20 # min_condition & max_condition are the minimum and maximum conditions (biological ages) an adult female can attain. min_condition = 50 max_condition = 300 # min_depoff_age & max_depoff_age are the minimum and maximum ages of offspring dependent upon their mother. min_depoff = 0 max_depoff = 14 # size_interval = increment size used in the discretisation of body weights. size_interval = 0.5 # condition_interval = increment size used in the discretisation of condition. condition_interval = 1 # depoff_interval = increment size used in the discretisation of dependent offspring age. depoff_interval = 1 # mass_to_energy = energy (in kcal) obtained from the mobilisation of 1kg of body weight. mass_to_energy = 5600.0 # energy_to_mass = energy (in kcal) needed to put on 1kg of body weight. energy_to_mass = 8750.0 # ageing_rate_exponent = theta in the equation relating change in condition (ageing rate) to investment in maintenance. ageing_rate_exponent = 1.0 # IMR = initial mortality risk upon reaching adulthood. IMR = 0.00285 # RoA = rate of ageing: the increase in mortality rate with age. RoA = 0.123 # Parameters used in calculation of basal metabolic rate. BMR_multiplicand = 2520.0 BMR_exponent = 0.75 # dim_size, dim_condition & dim_depoff correspond to the size of each dimension in the discretised state space. dim_size = (max_size - min_size + 1)/size_interval dim_condition = (max_condition - min_condition + 1)/condition_interval dim_depoff = (max_depoff - min_depoff + 1)/depoff_interval import random import math # Find the best implementation available on this platform. try: from cStringIO import StringIO except: from StringIO import StringIO output_file = open('./output_female_forward_cohort.txt', 'wb') # The list P holds the states of individuals in a cohort. Each individual is represented by a tuple. # The first three entries in each tuple correspond to size, condition and age of dependent offspring. # At the start all individuals have just reached adulthood and are therefore of size 8.5kg, in condition 50 periods and do not have any dependent offspring. # The fourth entry in each tuple is set to 1 at time periods at the end of which the individual has given birth. At other times this is 0. # The fifth entry in each tuple records the times when an individual gives birth. At time periods other than when the individual gives birth this is set to 0. # The sixth entry in each tuple is set to 1 at time periods at the end of which the individual has successfully reared an offspring to independence. At other times this is 0. # The seventh entry in each tuple records the time period when an individual dies. At time periods other than when the individual dies this is set to 0. P = [] for i in range(0, 100, 1): P.append((8.5, 50, 0, 0, 0, 0, 0)) # Retrieve time-independent stable strategy from the file female_forward_input.txt. strategy_input = open('female_forward_input.txt', 'r') allocation_decisions = list(strategy_input) strategy_input.close() # Set time for first iteration. time = min_time # Write to a buffer. output_bf = StringIO() # Print the time and states of all organisms to the buffer. output_bf.write(str(time) + '\t',) output_file.write(str(time) + '\t',) for i in range(0, len(P)): for j in range(0, len(P[0]), 1): output_bf.write(str(P[i][j]) + '\t', ) output_file.write(str(P[i][j]) + '\t', ) # Retrieve data from the buffer. print(output_bf.getvalue()) output_bf.close() # Cycle over time points. while time <= max_time: # Write to a buffer. output_bf = StringIO() # Cycle over individuals in population. for p in range(0, len(P), 1): # If size is zero, the individual is dead and the algorithm moves onto the next individual. if P[p][0] == 0: P[p] = (0, 0, 0, 0, 0, 0, 0) else: # Calculate probability of surviving the current period. survival_probability = 1 - IMR*math.exp(RoA*P[p][1]/10) # Randomly select a number between zero and one and assign this value to the variable 'random1'. # If random1 if less than or equal to the survival probability, the individual survives. If not, she dies. # If the indiviual dies, her state is set to (0, 0, 0, 0, 0, 0, time) in P. random1 = random.random() if random1 > survival_probability: P[p] = (0, 0, 0, 0, 0, 0, time) else: # If the individual survives, read optimal strategy from backwards iteration output. current_strategy = allocation_decisions[int((dim_condition - 1)*dim_depoff*(P[p][0] - min_size - size_interval)/size_interval \ + dim_depoff*(P[p][1] - min_condition)/condition_interval + P[p][2])] depoff_investment = float(current_strategy.partition('\t')[0]) growth_and_maintenance = current_strategy.partition('\t')[2] size_investment = float(growth_and_maintenance.partition('\t')[0]) condition_investment = float(growth_and_maintenance.partition('\t')[2]) # Compute the expected state the individual would be in at the start of the next time point following this allocation strategy. if size_investment >= 0: size_next = min(P[p][0] + size_investment/energy_to_mass, max_size) else: size_next = max(P[p][0] + size_investment/mass_to_energy, min_size) half_BMR = 0.5*BMR_multiplicand*(P[p][0]**BMR_exponent) condition_next = min(P[p][1] + (half_BMR/(condition_investment))**ageing_rate_exponent, max_condition) depoff_next = ((P[p][2] + 1)%dim_depoff)*(1 - offspring_loss[int(P[p][2])])*(depoff_investment/(offspring_energy[int(P[p][2])])) # This exact state the organism would be in at the start of the next period may not fall on the discretised values of states used. # The organism is therefore assigned a state which will be one of the eight discretised values which form a cube around the exact value. # The probability of being assigned a given state is weighted according to how close the exact value is to the discretised one. size_lower = size_interval*(size_next//size_interval) size_upper = min(size_lower + size_interval, max_size) condition_lower = condition_interval*(condition_next//condition_interval) condition_upper = min(condition_lower + condition_interval, max_condition) depoff_lower = 0 depoff_upper = (P[p][2] + depoff_interval)%dim_depoff size_weighting = (size_next - size_lower)/size_interval condition_weighting = (condition_next - condition_lower)/condition_interval depoff_weighting = (1 - offspring_loss[int(P[p][2])])*(depoff_investment/(offspring_energy[int(P[p][2])])) # To keep track of the number of offspring raised to independence, define a variable offspring_potentially_reared. # This is set to one if an offspring is 14 periods old at the start of the current period and the mother invests in it. if P[p][2] == 14: offspring_potentially_reared = depoff_investment/offspring_energy[max_depoff] else: offspring_potentially_reared = 0 # Probabilities of being in each of eight possible states in the next period. prob_slo_clo_dlo = (1 - size_weighting)*(1 - condition_weighting)*(1 - depoff_weighting) prob_slo_clo_dup = (1 - size_weighting)*(1 - condition_weighting)*depoff_weighting prob_slo_cup_dlo = (1 - size_weighting)*condition_weighting*(1 - depoff_weighting) prob_slo_cup_dup = (1 - size_weighting)*condition_weighting*depoff_weighting prob_sup_clo_dlo = size_weighting*(1 - condition_weighting)*(1 - depoff_weighting) prob_sup_clo_dup = size_weighting*(1 - condition_weighting)*depoff_weighting prob_sup_cup_dlo = size_weighting*condition_weighting*(1 - depoff_weighting) prob_sup_cup_dup = size_weighting*condition_weighting*depoff_weighting # Randomly select a number between zero and one and assign this value to the variable 'random2'. # If random2 is less than the probability of taking the first of the eight discretised states, the organism will take this state in the next period. # If random2 is greater than the probability of taking the first discretised state abut less than the combined probability of taking the first or second # states, the individual will take the second state in the next period. And so on... # If the individual has raised an offspring to independence during the current time period, the sixth tuple element is set to one. # An offspring reaches independence if offspring_potentially_reared = one and the offspring survives the current period, the latter being the probability # of taking depoff_upper rather than depoff_lower. random2 = random.random() if random2 < (prob_slo_clo_dlo): P[p] = (size_lower, condition_lower, depoff_lower, 0, 0, 0, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup): P[p] = (size_lower, condition_lower, depoff_upper, 0, 0, offspring_potentially_reared, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup + prob_slo_cup_dlo): P[p] = (size_lower, condition_upper, depoff_lower, 0, 0, 0, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup + prob_slo_cup_dlo + prob_slo_cup_dup): P[p] = (size_lower, condition_upper, depoff_upper, 0, 0, offspring_potentially_reared, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup + prob_slo_cup_dlo + prob_slo_cup_dup + prob_sup_clo_dlo): P[p] = (size_upper, condition_lower, depoff_lower, 0, 0, 0, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup + prob_slo_cup_dlo + prob_slo_cup_dup + prob_sup_clo_dlo + prob_sup_clo_dup): P[p] = (size_upper, condition_lower, depoff_upper, 0, 0, offspring_potentially_reared, 0) elif random2 < (prob_slo_clo_dlo + prob_slo_clo_dup + prob_slo_cup_dlo + prob_slo_cup_dup + prob_sup_clo_dlo + prob_sup_clo_dup + prob_sup_cup_dlo): P[p] = (size_upper, condition_upper, depoff_lower, 0, 0, 0, 0) else: P[p] = (size_upper, condition_upper, depoff_upper, 0, 0, offspring_potentially_reared, 0) # If the individual has given birth at the end of the current time period, set the fourth tuple element to one and the fifth tuple element to time. # Otherwise set these to zero. if P[p][2] == 5: P[p] = (P[p][0], P[p][1], P[p][2], 1, time, P[p][5], P[p][6]) else: P[p] = (P[p][0], P[p][1], P[p][2], 0, 0, P[p][5], P[p][6]) # If size = 8 and/or condition = 300, the individual is dead and state is set to (0, P[p][1], P[p][2], P[p][3], P[p][4], P[p][5], time). if P[p][0] == 8: P[p] = (0, P[p][1], P[p][2], P[p][3], P[p][4], P[p][5], time) if P[p][1] == 300: P[p] = (0, P[p][1], P[p][2], P[p][3], P[p][4], P[p][5], time) # Print the time and states of all organisms to the buffer. output_bf.write(str(time) + '\t',) output_file.write(str(time) + '\t',) for i in range(0, len(P)): for j in range(0, len(P[0]), 1): output_bf.write(str(P[i][j]) + '\t', ) output_file.write(str(P[i][j]) + '\t', ) # Retrieve data from the buffer. print(output_bf.getvalue()) output_bf.close() # Replace t with (t + 1). time = time + 1 print('-- Computation End --\n');