Traktor/myenv/Lib/site-packages/sympy/stats/random_matrix_models.py

458 lines
15 KiB
Python
Raw Permalink Normal View History

2024-05-26 05:12:46 +02:00
from sympy.concrete.products import Product
from sympy.concrete.summations import Sum
from sympy.core.basic import Basic
from sympy.core.function import Lambda
from sympy.core.numbers import (I, pi)
from sympy.core.singleton import S
from sympy.core.symbol import Dummy
from sympy.functions.elementary.complexes import Abs
from sympy.functions.elementary.exponential import exp
from sympy.functions.special.gamma_functions import gamma
from sympy.integrals.integrals import Integral
from sympy.matrices.expressions.matexpr import MatrixSymbol
from sympy.matrices.expressions.trace import Trace
from sympy.tensor.indexed import IndexedBase
from sympy.core.sympify import _sympify
from sympy.stats.rv import _symbol_converter, Density, RandomMatrixSymbol, is_random
from sympy.stats.joint_rv_types import JointDistributionHandmade
from sympy.stats.random_matrix import RandomMatrixPSpace
from sympy.tensor.array import ArrayComprehension
__all__ = [
'CircularEnsemble',
'CircularUnitaryEnsemble',
'CircularOrthogonalEnsemble',
'CircularSymplecticEnsemble',
'GaussianEnsemble',
'GaussianUnitaryEnsemble',
'GaussianOrthogonalEnsemble',
'GaussianSymplecticEnsemble',
'joint_eigen_distribution',
'JointEigenDistribution',
'level_spacing_distribution'
]
@is_random.register(RandomMatrixSymbol)
def _(x):
return True
class RandomMatrixEnsembleModel(Basic):
"""
Base class for random matrix ensembles.
It acts as an umbrella and contains
the methods common to all the ensembles
defined in sympy.stats.random_matrix_models.
"""
def __new__(cls, sym, dim=None):
sym, dim = _symbol_converter(sym), _sympify(dim)
if dim.is_integer == False:
raise ValueError("Dimension of the random matrices must be "
"integers, received %s instead."%(dim))
return Basic.__new__(cls, sym, dim)
symbol = property(lambda self: self.args[0])
dimension = property(lambda self: self.args[1])
def density(self, expr):
return Density(expr)
def __call__(self, expr):
return self.density(expr)
class GaussianEnsembleModel(RandomMatrixEnsembleModel):
"""
Abstract class for Gaussian ensembles.
Contains the properties common to all the
gaussian ensembles.
References
==========
.. [1] https://en.wikipedia.org/wiki/Random_matrix#Gaussian_ensembles
.. [2] https://arxiv.org/pdf/1712.07903.pdf
"""
def _compute_normalization_constant(self, beta, n):
"""
Helper function for computing normalization
constant for joint probability density of eigen
values of Gaussian ensembles.
References
==========
.. [1] https://en.wikipedia.org/wiki/Selberg_integral#Mehta's_integral
"""
n = S(n)
prod_term = lambda j: gamma(1 + beta*S(j)/2)/gamma(S.One + beta/S(2))
j = Dummy('j', integer=True, positive=True)
term1 = Product(prod_term(j), (j, 1, n)).doit()
term2 = (2/(beta*n))**(beta*n*(n - 1)/4 + n/2)
term3 = (2*pi)**(n/2)
return term1 * term2 * term3
def _compute_joint_eigen_distribution(self, beta):
"""
Helper function for computing the joint
probability distribution of eigen values
of the random matrix.
"""
n = self.dimension
Zbn = self._compute_normalization_constant(beta, n)
l = IndexedBase('l')
i = Dummy('i', integer=True, positive=True)
j = Dummy('j', integer=True, positive=True)
k = Dummy('k', integer=True, positive=True)
term1 = exp((-S(n)/2) * Sum(l[k]**2, (k, 1, n)).doit())
sub_term = Lambda(i, Product(Abs(l[j] - l[i])**beta, (j, i + 1, n)))
term2 = Product(sub_term(i).doit(), (i, 1, n - 1)).doit()
syms = ArrayComprehension(l[k], (k, 1, n)).doit()
return Lambda(tuple(syms), (term1 * term2)/Zbn)
class GaussianUnitaryEnsembleModel(GaussianEnsembleModel):
@property
def normalization_constant(self):
n = self.dimension
return 2**(S(n)/2) * pi**(S(n**2)/2)
def density(self, expr):
n, ZGUE = self.dimension, self.normalization_constant
h_pspace = RandomMatrixPSpace('P', model=self)
H = RandomMatrixSymbol('H', n, n, pspace=h_pspace)
return Lambda(H, exp(-S(n)/2 * Trace(H**2))/ZGUE)(expr)
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S(2))
def level_spacing_distribution(self):
s = Dummy('s')
f = (32/pi**2)*(s**2)*exp((-4/pi)*s**2)
return Lambda(s, f)
class GaussianOrthogonalEnsembleModel(GaussianEnsembleModel):
@property
def normalization_constant(self):
n = self.dimension
_H = MatrixSymbol('_H', n, n)
return Integral(exp(-S(n)/4 * Trace(_H**2)))
def density(self, expr):
n, ZGOE = self.dimension, self.normalization_constant
h_pspace = RandomMatrixPSpace('P', model=self)
H = RandomMatrixSymbol('H', n, n, pspace=h_pspace)
return Lambda(H, exp(-S(n)/4 * Trace(H**2))/ZGOE)(expr)
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S.One)
def level_spacing_distribution(self):
s = Dummy('s')
f = (pi/2)*s*exp((-pi/4)*s**2)
return Lambda(s, f)
class GaussianSymplecticEnsembleModel(GaussianEnsembleModel):
@property
def normalization_constant(self):
n = self.dimension
_H = MatrixSymbol('_H', n, n)
return Integral(exp(-S(n) * Trace(_H**2)))
def density(self, expr):
n, ZGSE = self.dimension, self.normalization_constant
h_pspace = RandomMatrixPSpace('P', model=self)
H = RandomMatrixSymbol('H', n, n, pspace=h_pspace)
return Lambda(H, exp(-S(n) * Trace(H**2))/ZGSE)(expr)
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S(4))
def level_spacing_distribution(self):
s = Dummy('s')
f = ((S(2)**18)/((S(3)**6)*(pi**3)))*(s**4)*exp((-64/(9*pi))*s**2)
return Lambda(s, f)
def GaussianEnsemble(sym, dim):
sym, dim = _symbol_converter(sym), _sympify(dim)
model = GaussianEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def GaussianUnitaryEnsemble(sym, dim):
"""
Represents Gaussian Unitary Ensembles.
Examples
========
>>> from sympy.stats import GaussianUnitaryEnsemble as GUE, density
>>> from sympy import MatrixSymbol
>>> G = GUE('U', 2)
>>> X = MatrixSymbol('X', 2, 2)
>>> density(G)(X)
exp(-Trace(X**2))/(2*pi**2)
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = GaussianUnitaryEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def GaussianOrthogonalEnsemble(sym, dim):
"""
Represents Gaussian Orthogonal Ensembles.
Examples
========
>>> from sympy.stats import GaussianOrthogonalEnsemble as GOE, density
>>> from sympy import MatrixSymbol
>>> G = GOE('U', 2)
>>> X = MatrixSymbol('X', 2, 2)
>>> density(G)(X)
exp(-Trace(X**2)/2)/Integral(exp(-Trace(_H**2)/2), _H)
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = GaussianOrthogonalEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def GaussianSymplecticEnsemble(sym, dim):
"""
Represents Gaussian Symplectic Ensembles.
Examples
========
>>> from sympy.stats import GaussianSymplecticEnsemble as GSE, density
>>> from sympy import MatrixSymbol
>>> G = GSE('U', 2)
>>> X = MatrixSymbol('X', 2, 2)
>>> density(G)(X)
exp(-2*Trace(X**2))/Integral(exp(-2*Trace(_H**2)), _H)
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = GaussianSymplecticEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
class CircularEnsembleModel(RandomMatrixEnsembleModel):
"""
Abstract class for Circular ensembles.
Contains the properties and methods
common to all the circular ensembles.
References
==========
.. [1] https://en.wikipedia.org/wiki/Circular_ensemble
"""
def density(self, expr):
# TODO : Add support for Lie groups(as extensions of sympy.diffgeom)
# and define measures on them
raise NotImplementedError("Support for Haar measure hasn't been "
"implemented yet, therefore the density of "
"%s cannot be computed."%(self))
def _compute_joint_eigen_distribution(self, beta):
"""
Helper function to compute the joint distribution of phases
of the complex eigen values of matrices belonging to any
circular ensembles.
"""
n = self.dimension
Zbn = ((2*pi)**n)*(gamma(beta*n/2 + 1)/S(gamma(beta/2 + 1))**n)
t = IndexedBase('t')
i, j, k = (Dummy('i', integer=True), Dummy('j', integer=True),
Dummy('k', integer=True))
syms = ArrayComprehension(t[i], (i, 1, n)).doit()
f = Product(Product(Abs(exp(I*t[k]) - exp(I*t[j]))**beta, (j, k + 1, n)).doit(),
(k, 1, n - 1)).doit()
return Lambda(tuple(syms), f/Zbn)
class CircularUnitaryEnsembleModel(CircularEnsembleModel):
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S(2))
class CircularOrthogonalEnsembleModel(CircularEnsembleModel):
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S.One)
class CircularSymplecticEnsembleModel(CircularEnsembleModel):
def joint_eigen_distribution(self):
return self._compute_joint_eigen_distribution(S(4))
def CircularEnsemble(sym, dim):
sym, dim = _symbol_converter(sym), _sympify(dim)
model = CircularEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def CircularUnitaryEnsemble(sym, dim):
"""
Represents Circular Unitary Ensembles.
Examples
========
>>> from sympy.stats import CircularUnitaryEnsemble as CUE
>>> from sympy.stats import joint_eigen_distribution
>>> C = CUE('U', 1)
>>> joint_eigen_distribution(C)
Lambda(t[1], Product(Abs(exp(I*t[_j]) - exp(I*t[_k]))**2, (_j, _k + 1, 1), (_k, 1, 0))/(2*pi))
Note
====
As can be seen above in the example, density of CiruclarUnitaryEnsemble
is not evaluated because the exact definition is based on haar measure of
unitary group which is not unique.
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = CircularUnitaryEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def CircularOrthogonalEnsemble(sym, dim):
"""
Represents Circular Orthogonal Ensembles.
Examples
========
>>> from sympy.stats import CircularOrthogonalEnsemble as COE
>>> from sympy.stats import joint_eigen_distribution
>>> C = COE('O', 1)
>>> joint_eigen_distribution(C)
Lambda(t[1], Product(Abs(exp(I*t[_j]) - exp(I*t[_k])), (_j, _k + 1, 1), (_k, 1, 0))/(2*pi))
Note
====
As can be seen above in the example, density of CiruclarOrthogonalEnsemble
is not evaluated because the exact definition is based on haar measure of
unitary group which is not unique.
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = CircularOrthogonalEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def CircularSymplecticEnsemble(sym, dim):
"""
Represents Circular Symplectic Ensembles.
Examples
========
>>> from sympy.stats import CircularSymplecticEnsemble as CSE
>>> from sympy.stats import joint_eigen_distribution
>>> C = CSE('S', 1)
>>> joint_eigen_distribution(C)
Lambda(t[1], Product(Abs(exp(I*t[_j]) - exp(I*t[_k]))**4, (_j, _k + 1, 1), (_k, 1, 0))/(2*pi))
Note
====
As can be seen above in the example, density of CiruclarSymplecticEnsemble
is not evaluated because the exact definition is based on haar measure of
unitary group which is not unique.
"""
sym, dim = _symbol_converter(sym), _sympify(dim)
model = CircularSymplecticEnsembleModel(sym, dim)
rmp = RandomMatrixPSpace(sym, model=model)
return RandomMatrixSymbol(sym, dim, dim, pspace=rmp)
def joint_eigen_distribution(mat):
"""
For obtaining joint probability distribution
of eigen values of random matrix.
Parameters
==========
mat: RandomMatrixSymbol
The matrix symbol whose eigen values are to be considered.
Returns
=======
Lambda
Examples
========
>>> from sympy.stats import GaussianUnitaryEnsemble as GUE
>>> from sympy.stats import joint_eigen_distribution
>>> U = GUE('U', 2)
>>> joint_eigen_distribution(U)
Lambda((l[1], l[2]), exp(-l[1]**2 - l[2]**2)*Product(Abs(l[_i] - l[_j])**2, (_j, _i + 1, 2), (_i, 1, 1))/pi)
"""
if not isinstance(mat, RandomMatrixSymbol):
raise ValueError("%s is not of type, RandomMatrixSymbol."%(mat))
return mat.pspace.model.joint_eigen_distribution()
def JointEigenDistribution(mat):
"""
Creates joint distribution of eigen values of matrices with random
expressions.
Parameters
==========
mat: Matrix
The matrix under consideration.
Returns
=======
JointDistributionHandmade
Examples
========
>>> from sympy.stats import Normal, JointEigenDistribution
>>> from sympy import Matrix
>>> A = [[Normal('A00', 0, 1), Normal('A01', 0, 1)],
... [Normal('A10', 0, 1), Normal('A11', 0, 1)]]
>>> JointEigenDistribution(Matrix(A))
JointDistributionHandmade(-sqrt(A00**2 - 2*A00*A11 + 4*A01*A10 + A11**2)/2
+ A00/2 + A11/2, sqrt(A00**2 - 2*A00*A11 + 4*A01*A10 + A11**2)/2 + A00/2 + A11/2)
"""
eigenvals = mat.eigenvals(multiple=True)
if not all(is_random(eigenval) for eigenval in set(eigenvals)):
raise ValueError("Eigen values do not have any random expression, "
"joint distribution cannot be generated.")
return JointDistributionHandmade(*eigenvals)
def level_spacing_distribution(mat):
"""
For obtaining distribution of level spacings.
Parameters
==========
mat: RandomMatrixSymbol
The random matrix symbol whose eigen values are
to be considered for finding the level spacings.
Returns
=======
Lambda
Examples
========
>>> from sympy.stats import GaussianUnitaryEnsemble as GUE
>>> from sympy.stats import level_spacing_distribution
>>> U = GUE('U', 2)
>>> level_spacing_distribution(U)
Lambda(_s, 32*_s**2*exp(-4*_s**2/pi)/pi**2)
References
==========
.. [1] https://en.wikipedia.org/wiki/Random_matrix#Distribution_of_level_spacings
"""
return mat.pspace.model.level_spacing_distribution()