In [13]:

    

import pyevodyn.utils as utils
import pyevodyn.games as games
import numpy as np
import math


    

In [14]:

    

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))


    

In [21]:

    

np.ones(2, dtype=np.float128)


    

    
Out[21]:





array([ 1.0,  1.0], dtype=float128)

In [90]:

    

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)


    

In [91]:

    

%%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

In [92]:

    

print a


    

    




8.46314233069e-15

In [93]:

    

%%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

In [94]:

    

print a


    

    




8.46314233069e-15

In [55]:

    

matriz  = np.array(range(16)).reshape(4,4)


    

In [56]:

    

matriz


    

    
Out[56]:





array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

In [69]:

    

matriz[0:2:1,0:2:1]


    

    
Out[69]:





array([[0, 1],
       [4, 5]])

In [66]:

    

matriz[0:3:2,0:3:2]


    

    
Out[66]:





array([[ 0,  2],
       [ 8, 10]])

In [68]:

    

matriz[0:4:3,0:4:3]


    

    
Out[68]:





array([[ 0,  3],
       [12, 15]])

In [80]:

    

matriz[1::2,1::2]


    

    
Out[80]:





array([[ 5,  7],
       [13, 15]])

In [84]:

    

matriz[[2,0],:][:,[2,0]]


    

    
Out[84]:





array([[10,  8],
       [ 2,  0]])

In [87]:

    

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]]])


    

In [88]:

    

subarray(matriz, 2,0)


    

    
Out[88]:





array([[10,  8],
       [ 2,  0]])

In [89]:

    

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))


    

In [ ]: