import pyevodyn.utils as utils
import pyevodyn.games as games
import numpy as np
import math
def fixation_probability(mutant_index, resident_index, intensity_of_selection, population_size, payoff_function=None, game_matrix = None, number_of_strategies=None ,mapping='EXP', **kwargs):
suma = []
for k in xrange(1, population_size):
mult = []
for j in xrange(1, k + 1):
if (payoff_function is not None and game_matrix is None):
if (number_of_strategies==None):
raise ValueError('When using a custom payoff_function you must specify number_of_strategies.')
strategies = np.zeros(number_of_strategies, dtype=int)
strategies[mutant_index] = j
strategies[resident_index] = population_size - j
payoffMutant = payoff_function(mutant_index, population_composition=strategies, **kwargs)
payoffResident = payoff_function(resident_index, population_composition=strategies, **kwargs)
elif (game_matrix is not None and payoff_function is None):
(payoffMutant,payoffResident) = payoff_from_matrix(mutant_index,resident_index, game_matrix, j, population_size)
else:
raise ValueError('No valid payoff structure given, please specify a game_matrix or a payoff_function.')
if (mapping=='EXP'):
fitnessMutant = math.e ** (intensity_of_selection * payoffMutant)
fitnessResident = math.e ** (intensity_of_selection * payoffResident)
elif (mapping =='LIN'):
fitnessMutant = 1 - intensity_of_selection + intensity_of_selection*payoffMutant
fitnessResident = 1 - intensity_of_selection + intensity_of_selection*payoffResident
else:
raise ValueError('No valid mapping given. Use EXP or LIN for exponential and linear respectively.')
mult.append(fitnessResident/fitnessMutant)
suma.append(utils.kahan_product(mult))
if any(np.isinf(suma)):
return 0.0
try:
complex_expression = utils.kahan_sum(suma)
except OverflowError:
return 0.0
if np.isinf(complex_expression):
return 0.0
return 1.0 / (1.0 + complex_expression)
def payoff_from_matrix(focal_index, other_index, game_matrix, number_of_focal_individuals, population_size):
"""
Computes a vector of payoffs from a game_matrix. The first element is the payoff of the strategy with index focal_index, the second element is
the payoff of the strategy with index other_index. The game is given by game_matrix, and assumes a population cmposed of number_of_focal_individuals of strategy focal_index
and population_size - number_of_focal_individuals copies of strategy other_index
"""
sub_matrix = np.array([[game_matrix[focal_index, focal_index],game_matrix[focal_index, other_index]],[game_matrix[other_index, focal_index],game_matrix[other_index, other_index]]])
return (1.0/(population_size-1.0))*(np.dot(sub_matrix,np.array([number_of_focal_individuals,population_size-number_of_focal_individuals]))-np.diagonal(sub_matrix))
np.ones(2, dtype=np.float128)
array([ 1.0, 1.0], dtype=float128)
def fixation_probability_fast(mutant_index, resident_index, intensity_of_selection, population_size, payoff_function=None, game_matrix = None, number_of_strategies=None ,mapping='EXP', **kwargs):
suma = np.ones(population_size, dtype=np.float128)
for k in xrange(1, population_size):
mult = np.ones(k, dtype=np.float128)
for j in xrange(1, k + 1):
if (payoff_function is not None and game_matrix is None):
if (number_of_strategies==None):
raise ValueError('When using a custom payoff_function you must specify number_of_strategies.')
strategies = np.zeros(number_of_strategies, dtype=int)
strategies[mutant_index] = j
strategies[resident_index] = population_size - j
payoffMutant = payoff_function(mutant_index, population_composition=strategies, **kwargs)
payoffResident = payoff_function(resident_index, population_composition=strategies, **kwargs)
elif (game_matrix is not None and payoff_function is None):
(payoffMutant,payoffResident) = payoff_from_matrix_faster(mutant_index,resident_index, game_matrix, j, population_size)
else:
raise ValueError('No valid payoff structure given, please specify a game_matrix or a payoff_function.')
if (mapping=='EXP'):
fitnessMutant = math.e ** (intensity_of_selection * payoffMutant)
fitnessResident = math.e ** (intensity_of_selection * payoffResident)
elif (mapping =='LIN'):
fitnessMutant = 1 - intensity_of_selection + intensity_of_selection*payoffMutant
fitnessResident = 1 - intensity_of_selection + intensity_of_selection*payoffResident
else:
raise ValueError('No valid mapping given. Use EXP or LIN for exponential and linear respectively.')
mult[j-1]=(fitnessResident/fitnessMutant)
suma[k-1]=(np.prod(mult))
if np.any(np.isinf(suma)):
return 0.0
try:
complex_expression = np.sum(suma)
except OverflowError:
return 0.0
if np.isinf(complex_expression):
return 0.0
return 1.0 / (1.0 + complex_expression)
%%timeit
a = fixation_probability(0,2, 1.0, 500, game_matrix=games.allc_tft_alld())
1 loops, best of 3: 7.83 s per loop
print a
8.46314233069e-15
%%timeit
a = fixation_probability_fast(0,2, 1.0, 500, game_matrix=games.allc_tft_alld())
1 loops, best of 3: 8.32 s per loop
print a
8.46314233069e-15
matriz = np.array(range(16)).reshape(4,4)
matriz
array([[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11], [12, 13, 14, 15]])
matriz[0:2:1,0:2:1]
array([[0, 1], [4, 5]])
matriz[0:3:2,0:3:2]
array([[ 0, 2], [ 8, 10]])
matriz[0:4:3,0:4:3]
array([[ 0, 3], [12, 15]])
matriz[1::2,1::2]
array([[ 5, 7], [13, 15]])
matriz[[2,0],:][:,[2,0]]
array([[10, 8], [ 2, 0]])
def subarray(game_matrix, focal_index, other_index):
return np.array([[game_matrix[focal_index, focal_index],game_matrix[focal_index, other_index]],[game_matrix[other_index, focal_index],game_matrix[other_index, other_index]]])
subarray(matriz, 2,0)
array([[10, 8], [ 2, 0]])
def payoff_from_matrix_faster(focal_index, other_index, game_matrix, number_of_focal_individuals, population_size):
"""
Computes a vector of payoffs from a game_matrix. The first element is the payoff of the strategy with index focal_index, the second element is
the payoff of the strategy with index other_index. The game is given by game_matrix, and assumes a population cmposed of number_of_focal_individuals of strategy focal_index
and population_size - number_of_focal_individuals copies of strategy other_index
"""
sub_matrix = matriz[[focal_index,other_index],:][:,[focal_index,other_index]]
return (1.0/(population_size-1.0))*(np.dot(sub_matrix,np.array([number_of_focal_individuals,population_size-number_of_focal_individuals]))-np.diagonal(sub_matrix))