# Copyright (c) 2019-2021 Elias Fernandez
#
# This file is part of EGTtools.
#
# EGTtools is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# EGTtools is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with EGTtools. If not, see <http://www.gnu.org/licenses/>
"""
This python module contains some utility functions
to find saddle points and plot gradients in 2 player, 2 strategy games.
"""
from warnings import warn
import numpy as np
from scipy.linalg import schur, eigvals
from scipy.sparse import csr_matrix, csc_matrix, issparse
from scipy.sparse.linalg import eigs as sparse_eigs
from typing import Optional, List, Generator, Union, Callable, Dict
from networkx import Graph
from egttools.games import AbstractGame
[docs]
def find_saddle_type_and_gradient_direction(gradient: Union[List[float], np.ndarray],
saddle_points_idx: Union[List[int], np.ndarray],
offset: float = 0.01):
"""
Finds whether a saddle point is stable or not. And defines the direction of the
gradient among stable and unstable points.
Parameters
----------
gradient : Union[List[float], numpy.ndarray[numpy.float64[1,m]]]
array containing the gradient of selection for all states of the population
saddle_points_idx : Union[List[int], numpy.ndarray[numpy.int64[1,m]]]
array containing the saddle points indices
offset : float
offset for the gradient_directions, so that arrows don't overlap with point
Returns
-------
Tuple[List[bool], List[float]]
Tuple containing an array that indicates the type of saddle points and another array indicating
the direction of the gradient between unstable and stable points
"""
saddle_type = []
gradient_direction = []
nb_points = len(gradient)
real_offset = offset * (nb_points - 1)
for i, point in enumerate(saddle_points_idx):
if point < nb_points - 1:
if point > 0:
if gradient[point + 1] > 0:
saddle_type.append(False)
gradient_direction.append((point, saddle_points_idx[i + 1] - real_offset))
elif gradient[point - 1] < 0:
saddle_type.append(False)
gradient_direction.append((point, saddle_points_idx[i - 1] + real_offset))
else:
saddle_type.append(True)
else:
if gradient[point + 1] > 0:
saddle_type.append(False)
gradient_direction.append((point, saddle_points_idx[i + 1] - real_offset))
else:
saddle_type.append(True)
else:
if gradient[point - 1] < 0:
saddle_type.append(False)
gradient_direction.append((point, saddle_points_idx[i - 1] + real_offset))
else:
saddle_type.append(True)
saddle_type = np.asarray(saddle_type)
gradient_direction = np.asarray(gradient_direction) / (nb_points - 1)
return saddle_type, gradient_direction
[docs]
def get_payoff_function(strategy_i: int,
strategy_j: int,
nb_strategies: int,
game: AbstractGame) -> Callable[[int, int, Optional[List]], float]:
"""
Returns a function which gives the payoff of strategy i against strategy j.
The returned function will return the payoff of strategy i
given k individuals of strategy i and group_size - k j strategists.
Parameters
----------
strategy_i : int
index of strategy i
strategy_j : int
index of strategy j
nb_strategies : int
Total number of strategies in the population.
game: egttools.games.AbstractGame
A game object which contains the method `egttools.games.AbstractGame.payoff` which returns the payoff of
a strategy given a group composition.
Returns
-------
Callable[[int, int, Optional[List]], float]
A function which will return the payoff of strategy i,
given k individuals of strategy i and group_size - k j strategists.
"""
def get_payoff(k: int, group_size: int, *args: Optional[List]) -> float:
"""
Returns the payoff given k individuals of strategy i and group_size - k j strategists.
This function is useful when we can assume that there will only be two
strategies in a group at any moment in time.
Parameters
----------
k : int
Number of players adopting strategy i in the group.
group_size: int
Total size of the group.
args: Optional[List]
Extra arguments which may be required ot calculate the payoff
Returns
-------
float
The payoff of strategy i in a group with k i strategists and group_size - k
j strategists.
"""
if k > group_size:
raise Exception("You have indicated a wrong group composition. k must be smaller or equal to group_size.")
group_composition = np.zeros(shape=(nb_strategies,), dtype=int)
if strategy_i != strategy_j:
group_composition[strategy_i] = k
group_composition[strategy_j] = group_size - k
else:
group_composition[strategy_i] = group_size
return game.payoff(strategy_i, group_composition.tolist())
return get_payoff
[docs]
def calculate_stationary_distribution(transition_matrix: Union[np.ndarray, csr_matrix, csc_matrix]) -> np.ndarray:
"""
Calculates stationary distribution from a transition matrix of Markov chain.
The use of this function is not recommended if the matrix is non-Hermitian. Please use
calculate_stationary_distribution_non_hermitian instead in this case.
The stationary distribution is the normalized eigenvector associated with the eigenvalue 1
Parameters
----------
transition_matrix : Union[numpy.ndarray, scipy.sparse.csr_matrix, scipy.sparse.csc_matrix]
A 2 dimensional transition matrix
Returns
-------
numpy.ndarray
A 1-dimensional vector containing the stationary distribution
See Also
--------
egttools.utils.calculate_stationary_distribution_non_hermitian
"""
if (type(transition_matrix) == csr_matrix) or (type(transition_matrix) == csc_matrix):
tmp = transition_matrix.toarray()
else:
tmp = transition_matrix
# Check if there is any transition with value 1 - this would mean that the game is degenerate
if np.isclose(tmp, 1., atol=1e-11).any():
warn(
"Some of the entries in the transition matrix are close to 1 (with a tolerance of 1e-11). "
"This could result in more than one eigenvalue of magnitute 1 "
"(the Markov Chain is degenerate), so please be careful when analysing the results.", RuntimeWarning)
# calculate stationary distributions using eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(tmp)
index_stationary = np.argmin(abs(eigenvalues - 1.0)) # look for the element closest to 1 in the list of eigenvalues
sd = abs(eigenvectors[:, index_stationary].T.real) # it is essential to access the matrix by column
return sd / sd.sum() # normalize
[docs]
def calculate_stationary_distribution_non_hermitian(
transition_matrix: Union[np.ndarray, csr_matrix, csc_matrix]) -> np.ndarray:
"""
Calculates stationary distribution from a transition matrix of Markov chain which is not hermitian.
The stationary distribution is the normalized eigenvector associated with the eigenvalue 1
Parameters
----------
transition_matrix : Union[numpy.ndarray, scipy.sparse.csr_matrix, scipy.sparse.csc_matrix]
A 2 dimensional transition matrix
Returns
-------
numpy.ndarray
A 1-dimensional vector containing the stationary distribution
See Also
--------
egttools.utils.calculate_stationary_distribution
"""
if (type(transition_matrix) == csr_matrix) or (type(transition_matrix) == csc_matrix):
tmp = transition_matrix.toarray()
else:
tmp = transition_matrix
# Check if there is any transition with value 1 - this would mean that the game is degenerate
if np.isclose(tmp, 1., atol=1e-11).any():
warn(
"Some of the entries in the transition matrix are close to 1 (with a tolerance of 1e-11). "
"This could result in more than one eigenvalue of magnitute 1 "
"(the Markov Chain is degenerate), so please be careful when analysing the results.", RuntimeWarning)
# calculate stationary distributions using eigenvalues and eigenvectors
# noinspection PyTupleAssignmentBalance
schur_form, eigenvectors = schur(tmp)
eigenvalues = eigvals(schur_form)
index_stationary = np.argmin(abs(eigenvalues - 1.0)) # look for the element closest to 1 in the list of eigenvalues
sd = abs(eigenvectors[:, index_stationary].T.real) # it is essential to access the matrix by column
return sd / sd.sum() # normalize
[docs]
def calculate_spectral_gap(
transition_matrix: Union[np.ndarray, csr_matrix, csc_matrix],
tol: float = 1e-10,
) -> float:
"""Compute the spectral gap of a Markov chain transition matrix.
The spectral gap is ``1 - |lambda_2|`` where ``|lambda_2|`` is the second-largest
eigenvalue modulus (the largest eigenvalue of an ergodic chain is always 1).
A larger gap means faster mixing. The chain needs approximately
``ceil(log(epsilon) / log(1 - gap))`` steps from the worst initial state to be
epsilon-close to stationarity, providing a principled lower bound for the
``transitory`` parameter in Monte Carlo estimation.
Uses ``scipy.sparse.linalg.eigs`` (ARPACK) to avoid dense conversion; falls back
to a 2x2 dense solve only when the matrix is too small for ARPACK (n <= 3).
Parameters
----------
transition_matrix : numpy.ndarray, csr_matrix, or csc_matrix
Square row- or column-stochastic transition matrix of shape (n, n).
Eigenvalues are invariant to transposition so either convention works.
tol : float, optional
ARPACK convergence tolerance (default 1e-10).
Returns
-------
float
Spectral gap in (0, 1]. Returns 1.0 for a single-state chain (n <= 1).
"""
if issparse(transition_matrix):
P = transition_matrix.tocsr()
else:
P = csr_matrix(transition_matrix)
n = P.shape[0]
if n <= 1:
return 1.0
# ARPACK requires k < n; for n=2 or n=3 we can only ask for k=1 eigenvalue
# which is the dominant one (≈1). Return gap = 0 as a conservative bound.
if n <= 3:
return 0.0
try:
# Request 2 eigenvalues of largest magnitude from P (right eigenproblem).
# Eigenvalues of P and P^T coincide, so transposition convention doesn't matter.
vals, _ = sparse_eigs(P, k=2, which="LM", tol=tol, maxiter=10000)
except Exception:
return 0.0 # conservative: unknown gap
mods = np.sort(np.abs(vals))[::-1] # descending by magnitude
lambda2 = float(mods[1]) if len(mods) > 1 else 0.0
return float(1.0 - lambda2)
[docs]
def combine(values: List[Union[int, str]], length: int) -> Generator:
"""
Outputs a generator that will generate an ordered list
with the possible combinations of values with length.
Each time the generator is called it will output a list
of length :param length which contains a combinations of the elements in
the list of values.
Parameters
----------
values : List
elements to combine
length : int
size of the output
Returns
-------
Generator
A generator which outputs ordered combinations of value as a list of size
length
Examples
--------
>>> for value in combine([1, 2], 2):
... print(value)
[1, 1]
[2, 1]
[1, 2]
[2, 2]
"""
output = [values[0]] * length
max_combinations = len(values) ** length
yield output
for i in range(1, max_combinations):
for j in range(length):
output[j] = values[(i // len(values) ** j) % len(values)]
yield output
[docs]
def calculate_nb_unique_combinations(slots_per_bin: Union[List[int], np.ndarray]) -> int:
"""
Calculates the number of unique combinations given the required number
of elements of each group, which should be given in List format in the
slots_per_bin parameter.
Parameters
----------
slots_per_bin: Union[List[int], numpy.ndarray]
The list should contain the required number of elements of each group
that should be combined in a tuple of length sum(slots_per_bin).
Returns
-------
int
The total number of unique combinations of the groups in the available slots.
"""
total_slots = np.sum(slots_per_bin)
divider = np.prod([np.math.factorial(x) for x in slots_per_bin])
return np.math.factorial(total_slots) // divider