# This file was automatically generated by SWIG (http://www.swig.org).
# Version 2.0.11
#
# Do not make changes to this file unless you know what you are doing--modify
# the SWIG interface file instead.
from sys import version_info
if version_info >= (2,6,0):
def swig_import_helper():
from os.path import dirname
import imp
fp = None
try:
fp, pathname, description = imp.find_module('_likelihood', [dirname(__file__)])
except ImportError:
import _likelihood
return _likelihood
if fp is not None:
try:
_mod = imp.load_module('_likelihood', fp, pathname, description)
finally:
fp.close()
return _mod
_likelihood = swig_import_helper()
del swig_import_helper
else:
import _likelihood
del version_info
try:
_swig_property = property
except NameError:
pass # Python < 2.2 doesn't have 'property'.
def _swig_setattr_nondynamic(self,class_type,name,value,static=1):
if (name == "thisown"): return self.this.own(value)
if (name == "this"):
if type(value).__name__ == 'SwigPyObject':
self.__dict__[name] = value
return
method = class_type.__swig_setmethods__.get(name,None)
if method: return method(self,value)
if (not static):
self.__dict__[name] = value
else:
raise AttributeError("You cannot add attributes to %s" % self)
def _swig_setattr(self,class_type,name,value):
return _swig_setattr_nondynamic(self,class_type,name,value,0)
def _swig_getattr(self,class_type,name):
if (name == "thisown"): return self.this.own()
method = class_type.__swig_getmethods__.get(name,None)
if method: return method(self)
raise AttributeError(name)
def _swig_repr(self):
try: strthis = "proxy of " + self.this.__repr__()
except: strthis = ""
return "<%s.%s; %s >" % (self.__class__.__module__, self.__class__.__name__, strthis,)
try:
_object = object
_newclass = 1
except AttributeError:
class _object : pass
_newclass = 0
def _dcprogs_dtype():
return _likelihood._dcprogs_dtype()
_dcprogs_dtype = _likelihood._dcprogs_dtype
[docs]def eig(*args):
"""Computes eigenvalues of a matrix"""
return _likelihood.eig(*args)
[docs]def inv(*args):
"""Computes inverse of a non-singular matrix"""
return _likelihood.inv(*args)
[docs]def det(*args):
"""Computes determinant of a matrix"""
return _likelihood.det(*args)
def svd(*args):
return _likelihood.svd(*args)
svd = _likelihood.svd
[docs]def expm(*args):
"""Computes singular value decomposition of a matrix"""
return _likelihood.expm(*args)
[docs]class QMatrix(_object):
"""
Transition rate matrix.
Open-open transitions should be located in the top-left corner.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, QMatrix, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, QMatrix, name)
__swig_setmethods__["nopen"] = _likelihood.QMatrix_nopen_set
__swig_getmethods__["nopen"] = _likelihood.QMatrix_nopen_get
if _newclass:nopen = _swig_property(_likelihood.QMatrix_nopen_get, _likelihood.QMatrix_nopen_set)
__swig_setmethods__["matrix"] = _likelihood.QMatrix_matrix_set
__swig_getmethods__["matrix"] = _likelihood.QMatrix_matrix_get
if _newclass:matrix = _swig_property(_likelihood.QMatrix_matrix_get, _likelihood.QMatrix_matrix_set)
def __init__(self, *args):
this = _likelihood.new_QMatrix(*args)
try: self.this.append(this)
except: self.this = this
def __str__(self): return _likelihood.QMatrix___str__(self)
def __repr__(self): return _likelihood.QMatrix___repr__(self)
[docs] def transpose(self):
"""
Swaps AA block with FF, AF block with FA, nopen with nshut.
"""
return _likelihood.QMatrix_transpose(self)
def __getitem__(self, *args): return self.matrix.__getitem__(*args)
def __setitem__(self, *args): return self.matrix.__setitem__(*args)
__swig_destroy__ = _likelihood.delete_QMatrix
__del__ = lambda self : None;
QMatrix_swigregister = _likelihood.QMatrix_swigregister
QMatrix_swigregister(QMatrix)
QMatrix.aa = property(lambda self: self.matrix[:self.nopen, :self.nopen],
doc=""" Open to open transitions. """)
QMatrix.af = property(lambda self: self.matrix[:self.nopen, self.nopen:],
doc=""" Open to close transitions. """)
QMatrix.fa = property(lambda self: self.matrix[self.nopen:, :self.nopen],
doc=""" Open to close transitions. """)
QMatrix.ff = property(lambda self: self.matrix[self.nopen:, self.nopen:],
doc=""" Open to close transitions. """)
QMatrix.nshut = property(lambda self: self.matrix.shape[0] - self.nopen,
lambda self, value: \
setattr(self, "nopen", self.matrix.shape[0] - value),
doc=""" Number of shut states. """);
[docs]class IdealG(_object):
"""
Ideal Likelihood.
This object can be instantiated one of several way:
- With a matrix and an integer
>>> idealg = IdealG(array([...]), 2)
- With a QMatrix
>>> matrix = QMatrix(array([...]), 2)
>>> idealg = IdealG(matrix)
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, IdealG, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, IdealG, name)
def __init__(self, *args):
this = _likelihood.new_IdealG(*args)
try: self.this.append(this)
except: self.this = this
__swig_getmethods__["initial_occupancies"] = _likelihood.IdealG_initial_occupancies_get
if _newclass:initial_occupancies = _swig_property(_likelihood.IdealG_initial_occupancies_get)
__swig_getmethods__["final_occupancies"] = _likelihood.IdealG_final_occupancies_get
if _newclass:final_occupancies = _swig_property(_likelihood.IdealG_final_occupancies_get)
__swig_getmethods__["nopen"] = _likelihood.IdealG_nopen_get
if _newclass:nopen = _swig_property(_likelihood.IdealG_nopen_get)
__swig_getmethods__["nshut"] = _likelihood.IdealG_nshut_get
if _newclass:nshut = _swig_property(_likelihood.IdealG_nshut_get)
[docs] def af(self, *args):
"""
af(self, t) -> DCProgs::t_rmatrix
AF transitions with respect to time
Implements the ideal likelihood:
.. math::
e^{t\mathcal{Q}_{FF}}\mathcal{Q}_{FA}.
:param t:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.IdealG_af(self, *args)
[docs] def fa(self, *args):
"""
fa(self, t) -> DCProgs::t_rmatrix
FA transitions with respect to time
Implements the ideal likelihood:
.. math::
e^{t\mathcal{Q}_{AA}}\mathcal{Q}_{AF}.
:param t:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.IdealG_fa(self, *args)
[docs] def laplace_af(self, *args):
"""
laplace_af(self, s) -> DCProgs::t_rmatrix
AF transitions with respect to scale
Implements the laplace transform of the likelihood:
.. math::
(sI - \mathcal{Q}_{AA})^{-1}\mathcal{Q}_{AF}.
:param s:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.IdealG_laplace_af(self, *args)
[docs] def laplace_fa(self, *args):
"""
laplace_fa(self, s) -> DCProgs::t_rmatrix
FA transitions with respect to scale
Implements the laplace transform of the likelihood:
.. math::
(sI - \mathcal{Q}_{FF})^{-1}\mathcal{Q}_{FA}.
:param s:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.IdealG_laplace_fa(self, *args)
def __str__(self): return _likelihood.IdealG___str__(self)
def __repr__(self): return _likelihood.IdealG___repr__(self)
__swig_destroy__ = _likelihood.delete_IdealG
__del__ = lambda self : None;
IdealG_swigregister = _likelihood.IdealG_swigregister
IdealG_swigregister(IdealG)
[docs]class DeterminantEq(_object):
"""
Compute determinant W needed to approximate missed event G
This object can be instantiated from a square matrix, the number of open states, and the
resolution, or maximum length of missed events :math:`\\tau`:
>>> DeterminantEq(matrix, nopen, tau)
or, it can be instantiated from a :py:class:`~dcprogs.likelihood.QMatrix` instance
and :math:`\\tau`
>>> DeterminantEq(qmatrix, tau)
:param matrix:
Transition rate matrix, where the upper left corner contain open-open transitions
:param integer nopen:
Number of open states in the transition matrix.
:param qmatrix:
:py:class:`QMatrix` instance
:param number tau:
Max length of missed events.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, DeterminantEq, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, DeterminantEq, name)
def __init__(self, *args):
this = _likelihood.new_DeterminantEq(*args)
try: self.this.append(this)
except: self.this = this
def transpose(self): return _likelihood.DeterminantEq_transpose(self)
def _get_tau(self): return _likelihood.DeterminantEq__get_tau(self)
def _set_tau(self, *args): return _likelihood.DeterminantEq__set_tau(self, *args)
__swig_getmethods__["tau"] = _get_tau
__swig_setmethods__["tau"] = _set_tau
if _newclass: tau = property(_get_tau, _set_tau, doc="Max length of mixed events.")
[docs] def __call__(self, *args):
"""
__call__(self, _s) -> DCProgs::t_real
__call__(self, _s, _tau) -> DCProgs::t_real
Computes the determinant of :math:`W(s)`.
:param number s:
The laplace scale. It can be a scalar or a numpy array of any shape.
In the latter case, the return is a array of the same shape where
each element is a matrix corresponding to the element in the input
array.
:param number tau:
*Optional*. If present, it is the max length of missed events.
:returns: a numpy array.
"""
return _likelihood.DeterminantEq___call__(self, *args)
__call__.__doc__ = "Computes determinant W\n\n" \
"Parameters:\n" \
" s: scalar, tuple, list, array\n" \
" The laplace scale.\n" \
" tau: *optional* number\n" \
" If present, it is the max length of missed events.\n\n" \
"Returns: If a scalar, returns a scalar. " \
"Otherwise returns a numpy array."
[docs] def H(self, *args):
"""
H(self, s) -> DCProgs::t_rmatrix
H(self, s, tau) -> DCProgs::t_rmatrix
Computes the matrix H
H is defined as :math:`\mathcal{Q}_{AA} + \mathcal{Q}_{AF}\ \int_0^\\tau e^{-st}e^{\mathcal{Q}_{FF}t}\partial\,t\ \mathcal{Q}_{FA}`.
:param number s:
The laplace scale. It can be a scalar or a numpy array of any shape.
In the latter case, the return is a array of the same shape where
each element is a matrix corresponding to the element in the input
array.
:param number tau:
*Optional*. If present, it is the max length of missed events.
:returns: a numpy array.
"""
return _likelihood.DeterminantEq_H(self, *args)
[docs] def s_derivative(self, *args):
"""
s_derivative(self, s) -> DCProgs::t_rmatrix
s_derivative(self, s, tau) -> DCProgs::t_rmatrix
Computes the derivative of H versus s
H is defined as :math:`\mathcal{Q}_{AA} + \mathcal{Q}_{AF}\ \int_0^\\tau e^{-st}e^{\mathcal{Q}_{FF}t}\partial\,t\ \mathcal{Q}_{FA}`.
:param number s:
The laplace scale. It can be a scalar or a numpy array of any shape.
In the latter case, the return is a array of the same shape where
each element is a matrix corresponding to the element in the input
array.
:param number tau:
*Optional*. If present, it is the max length of missed events.
:returns: a numpy array.
"""
return _likelihood.DeterminantEq_s_derivative(self, *args)
def __str__(self): return _likelihood.DeterminantEq___str__(self)
def __repr__(self): return _likelihood.DeterminantEq___repr__(self)
__swig_destroy__ = _likelihood.delete_DeterminantEq
__del__ = lambda self : None;
DeterminantEq_swigregister = _likelihood.DeterminantEq_swigregister
DeterminantEq_swigregister(DeterminantEq)
class Asymptotes(_object):
"""
Computes asymptotic values of missed-event G.
This function tkaes on input a value or an array of values in laplace space and returns the asymptotic likelihood.
:param determinant:
A :class:`DeterminantEq` instance from which to compute asymptotic values.
:param roots:
A sequence of `(root, multiplicity)` tuples.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, Asymptotes, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, Asymptotes, name)
__repr__ = _swig_repr
def __init__(self, *args):
this = _likelihood.new_Asymptotes(*args)
try: self.this.append(this)
except: self.this = this
def __call__(self, *args): return _likelihood.Asymptotes___call__(self, *args)
__swig_destroy__ = _likelihood.delete_Asymptotes
__del__ = lambda self : None;
Asymptotes_swigregister = _likelihood.Asymptotes_swigregister
Asymptotes_swigregister(Asymptotes)
def _find_lower_bound_for_roots(*args):
return _likelihood._find_lower_bound_for_roots(*args)
_find_lower_bound_for_roots = _likelihood._find_lower_bound_for_roots
[docs]def find_lower_bound_for_roots(determinant, start=0e0, alpha=2e0, itermax=100):
""" Figures out lower bound for root.
The lower bound is obtained by iteratively checking the lowest eigenvalue
:math:`\\epsilon_i^s` of :math:`H(s_i)`, where :math:`s_i` is the guess at iteration
:math:`i`. If the lowest eigenvalue is smaller than :math:`s_i`, then
:math:`s_{i+1} = \\epsilon_i^s + \\alpha(\\epsilon_i^s - s_i)` is created.
:param DeterminantEq det:
Function for which to guess bound for lower root.
:param Number start:
Starting point.
:param Number alpha:
Factor from which to determine next best guess.
:param Integer itermax:
Maximum number of iterations before bailing out.
"""
return _likelihood._find_lower_bound_for_roots(determinant, start, alpha, itermax)
def _find_upper_bound_for_roots(*args):
return _likelihood._find_upper_bound_for_roots(*args)
_find_upper_bound_for_roots = _likelihood._find_upper_bound_for_roots
[docs]def find_upper_bound_for_roots(determinant, start=0e0, alpha=2e0, itermax=100):
""" Figures out upper bound for root.
The upper bound is obtained by iteratively checking the largest eigenvalue
:math:`\\epsilon_i^s` of :math:`H(s_i)`, where :math:`s_i` is the guess at iteration
:math:`i`. If the largest eigenvalue is larger than :math:`s_i`, then
:math:`s_{i+1} = \\epsilon_i^s + \\alpha(\\epsilon_i^s - s_i)` is created.
:param DeterminantEq det:
Function for which to guess bound for lower root.
:param Number start:
Starting point.
:param Number alpha:
Factor from which to determine next best guess.
:param Integer itermax:
Maximum number of iterations before bailing out.
"""
return _likelihood._find_upper_bound_for_roots(determinant, start, alpha, itermax)
def _find_root_intervals(*args):
return _likelihood._find_root_intervals(*args)
_find_root_intervals = _likelihood._find_root_intervals
[docs]def find_root_intervals(determinant, lower_bound=None, upper_bound=None, tolerance=1e-8):
""" Returns intervals for searching roots.
Finds a set of intervals between lower_bound and upper_bound, each containing a single root.
:param determinant:
:class:`DeterminantEq` instance for which to guess bound for lower root.
:param lower_bound: Number
Lower bound of the interval within which to search for roots. If None, the lower bound
is determined using :py:func:`find_lower_bound_for_roots`.
:param Number upper_bound:
Upper bound of the interval within which to search for roots. If None, the upper bound
is determined using :py:func:`find_upper_bound_for_roots`.
:param Number tolerance:
Size of the smallest possible intervals. Intervals smaller than this are likely to
have multiple roots.
:raises: ArithmeticError when NaN is encountered or when eigenvalue solver fails.
:returns:
A list `[(a, b), multiplicity]`, where `(a, b)` denotes an interval, and
`multiplicity` the multiplicity of the root. All roots with `multiplicity`:math:`> 1`
will have a size of `tolerance` or smaller.
"""
from numpy import NaN
if lower_bound is None: lower_bound = NaN
if upper_bound is None: upper_bound = NaN
return _likelihood._find_root_intervals(determinant, lower_bound, upper_bound, tolerance)
def _find_root_intervals_brute_force(*args):
return _likelihood._find_root_intervals_brute_force(*args)
_find_root_intervals_brute_force = _likelihood._find_root_intervals_brute_force
[docs]def find_root_intervals_brute_force( determinant, resolution=1e-1, lower_bound=None, upper_bound=None,
tolerance=1e-8 ):
""" Find intervals for roots using brute force.
Computes all values between lower_bound and upper_bound, for a given resolution. If the
determinant changes sign between two values, or if it comes to within tolerance of zero,
then the eigenvalues of H are used to determine possible multiplicity.
:param determinant:
Instance of :py:class:`DeterminantEq` for which to find root intervals
:param Number resolution:
Resolution at which to sample interval.
:param Number lower_bound:
Lower bound of interval. If None, then :py:func:`find_lower_bound_for_roots` is called.
:param Number upper_bound:
Upper bound of interval. If None, then :py:func:`find_upper_bound_for_roots` is called.
:param Number tolerance:
Tolerance below which the value of the determinant is considered *close* *to* *zero*.
:raises: ArithmeticError when NaN is encountered or when eigenvalue solver fails.\n\n"
:returns:
A list `[(a, b), multiplicity]`, where `(a, b)` denotes an interval, and `multiplicity`
the multiplicity of the root. All roots with `multiplicity > 1` will have a size of
`tolerance` or smaller.
"""
from numpy import NaN
if lower_bound is None: lower_bound = NaN
if upper_bound is None: upper_bound = NaN
return _likelihood._find_root_intervals_brute_force(
determinant, resolution, lower_bound, upper_bound, tolerance)
[docs]class ExactSurvivor(_object):
"""
Computes exact missed-event survivor function.
:param qmatrix:
The Q-matrix for which to compute the exact likelihood.
:param float double:
The maximum length of missed events.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, ExactSurvivor, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, ExactSurvivor, name)
def __init__(self, *args):
this = _likelihood.new_ExactSurvivor(*args)
try: self.this.append(this)
except: self.this = this
[docs] def af(self, *args):
"""
Open to close transitions
Open to close transitions
"""
return _likelihood.ExactSurvivor_af(self, *args)
[docs] def fa(self, *args):
"""
Closed to open transitions
Closed to open transitions
"""
return _likelihood.ExactSurvivor_fa(self, *args)
[docs] def recursion_af(self, *args):
""":math:`C_{iml}` recursion matrices for af."""
return _likelihood.ExactSurvivor_recursion_af(self, *args)
[docs] def recursion_fa(self, *args):
""":math:`C_{iml}` recursion matrices for fa."""
return _likelihood.ExactSurvivor_recursion_fa(self, *args)
[docs] def D_af(self, *args):
""":math:`D_i` matrices for af."""
return _likelihood.ExactSurvivor_D_af(self, *args)
[docs] def D_fa(self, *args):
""":math:`D_i` matrices for fa."""
return _likelihood.ExactSurvivor_D_fa(self, *args)
__swig_getmethods__["tau"] = _likelihood.ExactSurvivor_tau_get
if _newclass:tau = _swig_property(_likelihood.ExactSurvivor_tau_get)
__swig_getmethods__["eigenvalues_af"] = _likelihood.ExactSurvivor_eigenvalues_af_get
if _newclass:eigenvalues_af = _swig_property(_likelihood.ExactSurvivor_eigenvalues_af_get)
__swig_getmethods__["eigenvalues_fa"] = _likelihood.ExactSurvivor_eigenvalues_fa_get
if _newclass:eigenvalues_fa = _swig_property(_likelihood.ExactSurvivor_eigenvalues_fa_get)
def __str__(self): return _likelihood.ExactSurvivor___str__(self)
def __repr__(self): return _likelihood.ExactSurvivor___repr__(self)
__swig_destroy__ = _likelihood.delete_ExactSurvivor
__del__ = lambda self : None;
ExactSurvivor_swigregister = _likelihood.ExactSurvivor_swigregister
ExactSurvivor_swigregister(ExactSurvivor)
[docs]class ApproxSurvivor(_object):
"""
Computes approximate missed-event survivor function.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, ApproxSurvivor, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, ApproxSurvivor, name)
def __init__(self, *args, **kwargs):
""" Initializes the approximate survivor function.
There are three possible ways to instantiate this object:
>>> ApproxSurvivor(determinant_af, roots_af, determinant_fa, roots_fa)
>>> ApproxSurvivor(matrix, nopen, tau, **kwargs)
>>> ApproxSurvivor(qmatrix, tau, **kwargs)
The parameters between brackets are optional. The last two versions will try and calculate
the roots of the determinant equations automatically. A number of parameters can be given
to control this process. The keyword arguments are listed below.
:param determinant_af:
A :class:`DeterminantEq` instance, specifically for the af block.
:param roots_af:
The roots of the af determinant equation. The should come in the format
`[(root, multiplicity), (root, multiplicity), ...]`.
:param determinant_fa:
A :class:`DeterminantEq` instance, specifically for the fa block.
It should the transpose of `determinant_af`. It is required so
it need not be recomputed, since it most likely already exists.
:param roots_fa:
The roots of the fa determinant equation. The should come in the format
`[(root, multiplicity), (root, multiplicity), ...]`.
:param matrix:
An object convertible to a square matrix.
:param integer nopen:
Number of open states.
:param qmatrix:
A :class:`QMatrix` instance.
:param float xtol:
Tolerance criteria when computing roots using :func:`~dcprogs.likelihood.brentq`.
Defaults to 1e-12.
:param float rtol:
Tolerance criteria when computing roots using :func:`~dcprogs.likelihood.brentq`.
Defaults to 1e-12.
:param float itermax:
Maximum number of iterations when computing roots using
:func:`~dcprogs.likelihood.brentq`.
Defaults to 100.
:param float lower_bound:
Lower bound for *all* roots. Defaults to None, in which the case the lower bound is
computed from :py:func:`find_lower_bound_for_roots`.
:param float upper_bound:
Upper bound for *all* roots. Defaults to None, in which the case the upper bound is
computed from :py:func:`find_upper_bound_for_roots`.
"""
from numpy import NaN, array
from .. import internal_dtype
from . import QMatrix
# Try first initalization mode
det_af = args[0] if len(args) > 0 else kwargs.get("det_af", None)
if isinstance(det_af, DeterminantEq):
roots_af = args[1] if len(args) > 1 else kwargs.get("roots_af", None)
if roots_af is None: raise TypeError('Missing roots_af argument')
det_fa = args[2] if len(args) > 2 else kwargs.get("det_fa", None)
if det_fa is None: raise TypeError('Missing det_fa argument')
roots_fa = args[3] if len(args) > 3 else kwargs.get("roots_fa", None)
if roots_fa is None: raise TypeError('Missing roots_fa argument')
this = _likelihood.new_ApproxSurvivor(det_af, roots_af, det_fa, roots_fa)
try: self.this.append(this)
except: self.this = this
return
# Try second and third instanciation
matrix = args[0] if len(args) > 0 else kwargs.get('matrix', kwargs.get('qmatrix', None))
if matrix is None:
raise TypeError("Expected the q-matrix or a matrix + number of open states on input.")
index = 1
if not isinstance(matrix, QMatrix):
matrix = array(matrix, internal_dtype)
nopen = int(args[index]) if len(args) > index else kwargs.get('nopen', None)
if nopen is None:
raise TypeError("A matrix was given on input, but without a nopen parameter.")
matrix = QMatrix(matrix, nopen)
index += 1
tau = args[index] if len(args) > index else kwargs.get('tau', None)
if tau is None: raise TypeError("No tau given on input.")
xtol = args[index+1] if len(args) > index+1 else kwargs.get('xtol', 1e-12)
rtol = args[index+3] if len(args) > index+2 else kwargs.get('rtol', 1e-12)
itermax = int(args[index+3] if len(args) > index+3 else kwargs.get('itermax', 100))
lower_bound = args[index+4] if len(args) > index+4 else kwargs.get('lower_bound', None)
upper_bound = args[index+5] if len(args) > index+5 else kwargs.get('upper_bound', None)
if lower_bound is None: lower_bound = NaN
if upper_bound is None: upper_bound = NaN
this = _likelihood.new_ApproxSurvivor(matrix, tau, xtol, rtol, itermax, lower_bound, upper_bound)
try: self.this.append(this)
except: self.this = this
[docs] def af(self, *args):
"""
Open to close transitions
Open to close transitions
"""
return _likelihood.ApproxSurvivor_af(self, *args)
[docs] def fa(self, *args):
"""
Closed to open transitions
Closed to open transitions
"""
return _likelihood.ApproxSurvivor_fa(self, *args)
__swig_getmethods__["af_components"] = _likelihood.ApproxSurvivor_af_components_get
if _newclass:af_components = _swig_property(_likelihood.ApproxSurvivor_af_components_get)
__swig_getmethods__["fa_components"] = _likelihood.ApproxSurvivor_fa_components_get
if _newclass:fa_components = _swig_property(_likelihood.ApproxSurvivor_fa_components_get)
def __str__(self): return _likelihood.ApproxSurvivor___str__(self)
def __repr__(self): return _likelihood.ApproxSurvivor___repr__(self)
__swig_destroy__ = _likelihood.delete_ApproxSurvivor
__del__ = lambda self : None;
ApproxSurvivor_swigregister = _likelihood.ApproxSurvivor_swigregister
ApproxSurvivor_swigregister(ApproxSurvivor)
[docs]class MissedEventsG(_object):
"""
Computes missed-events likelihood.
Exact calculations take place for times smaller than :math:`n_{\mathrm{max}}\\tau`. Asymptotic calculations take over for larger times.
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, MissedEventsG, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, MissedEventsG, name)
def __init__(self, *args, **kwargs):
""" Initializes the missed-events likelihood.
Exact calculations take place for times smaller than `nmax * tau`.
Asymptotic calculations take over for larger times.
There are three possible ways to instantiate this object:
>>> MissedEvents(determinant_af, roots_af, determinant_fa, roots_fa[, nmax=2])
>>> MissedEvents(matrix, nopen, tau, **kwargs)
>>> MissedEvents(qmatrix, tau, **kwargs)
The parameters between brackets are optional. The last two versions will try
and calculate the roots of the determinant equations automatically. A number of parameters
can be given to control this process.
:param determinant_af:
A :class:`DeterminantEq` instance, specifically for the af block.
:param roots_af:
The roots of the af determinant equation. The should come in the format
`[(root, multiplicity), (root, multiplicity), ...]`.
:param determinant_fa:
A :class:`DeterminantEq` instance, specifically for the fa block.
It should the transpose of `determinant_af`. It is required so
it need not be recomputed, since it most likely already exists.
:param roots_fa:
The roots of the fa determinant equation. The should come in the format
`[(root, multiplicity), (root, multiplicity), ...]`.
:param matrix:
An object convertible to a square matrix.
:param integer nopen:
Number of open states.
:param qmatrix:
A :class:`QMatrix` instance.
:param int nmax:
The exact missed event likelihood will be computed for times
:math:`t \in [0, n_{\mathrm{max}} \tau`. Defaults to 3.
:param float xtol:
Tolerance criteria when computing roots using :func:`~dcprogs.likelihood.brentq`.
Defaults to 1e-12.
:param float rtol:
Tolerance criteria when computing roots using :func:`~dcprogs.likelihood.brentq`.
Defaults to 1e-12.
:param float itermax:
Maximum number of iterations when computing roots using
:func:`~dcprogs.likelihood.brentq`.
Defaults to 100.
:param float lower_bound:
Lower bound for *all* roots. Defaults to None, in which the case the lower bound is
computed from :py:func:`find_lower_bound_for_roots`.
:param float upper_bound:
Upper bound for *all* roots. Defaults to None, in which the case the upper bound is
computed from :py:func:`find_upper_bound_for_roots`.
"""
from numpy import NaN, array
from .. import internal_dtype
from . import QMatrix
# Try first initalization mode
det_af = args[0] if len(args) > 0 else kwargs.get("det_af", None)
if isinstance(det_af, DeterminantEq):
roots_af = args[1] if len(args) > 1 else kwargs.get("roots_af", None)
if roots_af is None: raise TypeError('Missing roots_af argument')
det_fa = args[2] if len(args) > 2 else kwargs.get("det_fa", None)
if det_fa is None: raise TypeError('Missing det_fa argument')
roots_fa = args[3] if len(args) > 3 else kwargs.get("roots_fa", None)
if roots_fa is None: raise TypeError('Missing roots_fa argument')
nmax = int(args[4] if len(args) > 4 else kwargs.get('nmax', 3))
this = _likelihood.new_MissedEventsG(det_af, roots_af, det_fa, roots_fa, nmax)
try: self.this.append(this)
except: self.this = this
return
# Try second and third instanciation
matrix = args[0] if len(args) > 0 else kwargs.get('matrix', kwargs.get('qmatrix', None))
if matrix is None:
raise ValueError("Expected the q-matrix or a matrix + number of open states on input.")
index = 1
if not isinstance(matrix, QMatrix):
matrix = array(matrix, internal_dtype)
nopen = int(args[index]) if len(args) > index else kwargs.get('nopen', None)
if nopen is None:
raise ValueError("A matrix was given on input, but without a nopen parameter.")
matrix = QMatrix(matrix, nopen)
index += 1
tau = args[index] if len(args) > index else kwargs.get('tau', None)
if tau is None: raise ValueError("No tau given on input.")
nmax = int(args[index+1] if len(args) > index+1 else kwargs.get('nmax', 3))
xtol = args[index+2] if len(args) > index+2 else kwargs.get('xtol', 1e-12)
rtol = args[index+3] if len(args) > index+3 else kwargs.get('rtol', 1e-12)
itermax = int(args[index+4] if len(args) > index+4 else kwargs.get('itermax', 100))
lower_bound = args[index+5] if len(args) > index+5 else kwargs.get('lower_bound', None)
upper_bound = args[index+6] if len(args) > index+6 else kwargs.get('upper_bound', None)
if lower_bound is None: lower_bound = NaN
if upper_bound is None: upper_bound = NaN
this = _likelihood.new_MissedEventsG(matrix, tau, nmax, xtol, rtol, itermax, lower_bound, upper_bound)
try: self.this.append(this)
except: self.this = this
[docs] def af(self, *args):
"""
af(self, t) -> DCProgs::t_rmatrix
Likelihood of an observed open time of length ``t``
:param t:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.MissedEventsG_af(self, *args)
[docs] def fa(self, *args):
"""
fa(self, t) -> DCProgs::t_rmatrix
Likelihood of a shut time of length ``t``
:param t:
A scalar or something to a numpy array. In the latter case, the return is an array
of matrices.
"""
return _likelihood.MissedEventsG_fa(self, *args)
[docs] def laplace_af(self, *args):
"""
laplace_af(self, s) -> DCProgs::t_rmatrix
Exact missed-events G in Laplace space.
The exact expression is :math:`^{e}\mathcal{G}_{AF}(s) = {}^AR(s) e^{-s\\tau}\mathcal{Q}_{AF}e^{\mathcal{Q}_{FF}\\tau}`, with :math:`{}^AR(s) = [sI - \mathcal{Q}_{AA} - \mathcal{Q}_{AF} \\int_0^\\tau e^{-st}e^{\mathcal{Q}_{FF}t}\\partial t \mathcal{Q}_{FA}]^{-1}`.
:param s:
The laplace scale. A real scalar or something convertible to a numpy array.
:returns:
A matrix if the input is scalar, an array of matrices otherwise,
with the shape of the input.
"""
return _likelihood.MissedEventsG_laplace_af(self, *args)
[docs] def laplace_fa(self, *args):
"""
laplace_fa(self, s) -> DCProgs::t_rmatrix
Exact missed-events G in Laplace space.
The exact expression is :math:`^{e}\mathcal{G}_{FA}(s) = {}^FR(s) e^{-s\\tau}\mathcal{Q}_{FA}e^{\mathcal{Q}_{AA}\\tau}`, with :math:`{}^FR(s) = [sI - \mathcal{Q}_{FF} - \mathcal{Q}_{FA} \\int_0^\\tau e^{-st}e^{\mathcal{Q}_{AA}t}\\partial t \mathcal{Q}_{AF}]^{-1}`.
:param s:
The laplace scale. A real scalar or something convertible to a numpy array.
:returns:
A matrix if the input is scalar, an array of matrices otherwise,
with the shape of the input.
"""
return _likelihood.MissedEventsG_laplace_fa(self, *args)
[docs] def initial_CHS_occupancies(self, *args):
"""Computes initial CHS occupancies."""
return _likelihood.MissedEventsG_initial_CHS_occupancies(self, *args)
[docs] def final_CHS_occupancies(self, *args):
"""Computes final CHS occupancies."""
return _likelihood.MissedEventsG_final_CHS_occupancies(self, *args)
__swig_getmethods__["tau"] = _likelihood.MissedEventsG_tau_get
if _newclass:tau = _swig_property(_likelihood.MissedEventsG_tau_get)
__swig_getmethods__["af_factor"] = _likelihood.MissedEventsG_af_factor_get
if _newclass:af_factor = _swig_property(_likelihood.MissedEventsG_af_factor_get)
__swig_getmethods__["fa_factor"] = _likelihood.MissedEventsG_fa_factor_get
if _newclass:fa_factor = _swig_property(_likelihood.MissedEventsG_fa_factor_get)
__swig_getmethods__["initial_occupancies"] = _likelihood.MissedEventsG_initial_occupancies_get
if _newclass:initial_occupancies = _swig_property(_likelihood.MissedEventsG_initial_occupancies_get)
__swig_getmethods__["final_occupancies"] = _likelihood.MissedEventsG_final_occupancies_get
if _newclass:final_occupancies = _swig_property(_likelihood.MissedEventsG_final_occupancies_get)
def __str__(self): return _likelihood.MissedEventsG___str__(self)
def __repr__(self): return _likelihood.MissedEventsG___repr__(self)
__swig_getmethods__["tmax"] = _likelihood.MissedEventsG_tmax_get
if _newclass:tmax = _swig_property(_likelihood.MissedEventsG_tmax_get)
__swig_setmethods__["nmax"] = _likelihood.MissedEventsG_nmax_set
__swig_getmethods__["nmax"] = _likelihood.MissedEventsG_nmax_get
if _newclass:nmax = _swig_property(_likelihood.MissedEventsG_nmax_get, _likelihood.MissedEventsG_nmax_set)
__swig_getmethods__["nopen"] = _likelihood.MissedEventsG_nopen_get
if _newclass:nopen = _swig_property(_likelihood.MissedEventsG_nopen_get)
__swig_getmethods__["nshut"] = _likelihood.MissedEventsG_nshut_get
if _newclass:nshut = _swig_property(_likelihood.MissedEventsG_nshut_get)
__swig_destroy__ = _likelihood.delete_MissedEventsG
__del__ = lambda self : None;
MissedEventsG_swigregister = _likelihood.MissedEventsG_swigregister
MissedEventsG_swigregister(MissedEventsG)
[docs]class Log10Likelihood(_object):
"""
Creates a functor for optimizing the likelihood.
This functor can as input a :math:`\mathcal{Q}`-matrix and returns the resulting
likelihood. The root finding process is performed by :py:func:`find_roots`. Most of the
arguments to initialize this functor are parameters for that process.
.. attribute:: nopen
Number of open states. Should be strictly positive.
.. attribute:: tau
Maximum length of missed events. Should be positive.
.. attribute:: tcritical
:math:`t_{\mathrm{crit}}` for CHS vectors. If negative, ``None``, or ``NaN``,
then the initial and final states are the equilibrium occupancies.
.. attribute:: nmax
The exact missed-events likelihood will be computed for
:math:`t < n_{\mathrm{max}} au`.
.. attribute:: xtol
Tolerance criteria when computing roots via :py:func:`likelihood.brentq`.
.. attribute:: rtol
Tolerance criteria when computing roots via :py:func:`likelihood.brentq`.
.. attribute:: itermax
Maximum number of iterations when computing roots via :py:func:`likelihood.brentq`.
.. attribute:: lower_bound
Lower bound bracketing all roots. If None or NaN, then the lower bound will be
computed during calls by :py:func:`~dcprogs.likelihood.find_lower_bound_for_root`
.. attribute:: upper_bound
Upper bound bracketing all roots. If None or NaN, then the upper bound will be
computed during calls by :py:func:`~dcprogs.likelihood.find_upper_bound_for_root`
"""
__swig_setmethods__ = {}
__setattr__ = lambda self, name, value: _swig_setattr(self, Log10Likelihood, name, value)
__swig_getmethods__ = {}
__getattr__ = lambda self, name: _swig_getattr(self, Log10Likelihood, name)
__swig_setmethods__["nopen"] = _likelihood.Log10Likelihood_nopen_set
__swig_getmethods__["nopen"] = _likelihood.Log10Likelihood_nopen_get
if _newclass:nopen = _swig_property(_likelihood.Log10Likelihood_nopen_get, _likelihood.Log10Likelihood_nopen_set)
__swig_setmethods__["tau"] = _likelihood.Log10Likelihood_tau_set
__swig_getmethods__["tau"] = _likelihood.Log10Likelihood_tau_get
if _newclass:tau = _swig_property(_likelihood.Log10Likelihood_tau_get, _likelihood.Log10Likelihood_tau_set)
__swig_setmethods__["tcritical"] = _likelihood.Log10Likelihood_tcritical_set
__swig_getmethods__["tcritical"] = _likelihood.Log10Likelihood_tcritical_get
if _newclass:tcritical = _swig_property(_likelihood.Log10Likelihood_tcritical_get, _likelihood.Log10Likelihood_tcritical_set)
__swig_setmethods__["nmax"] = _likelihood.Log10Likelihood_nmax_set
__swig_getmethods__["nmax"] = _likelihood.Log10Likelihood_nmax_get
if _newclass:nmax = _swig_property(_likelihood.Log10Likelihood_nmax_get, _likelihood.Log10Likelihood_nmax_set)
__swig_setmethods__["xtol"] = _likelihood.Log10Likelihood_xtol_set
__swig_getmethods__["xtol"] = _likelihood.Log10Likelihood_xtol_get
if _newclass:xtol = _swig_property(_likelihood.Log10Likelihood_xtol_get, _likelihood.Log10Likelihood_xtol_set)
__swig_setmethods__["rtol"] = _likelihood.Log10Likelihood_rtol_set
__swig_getmethods__["rtol"] = _likelihood.Log10Likelihood_rtol_get
if _newclass:rtol = _swig_property(_likelihood.Log10Likelihood_rtol_get, _likelihood.Log10Likelihood_rtol_set)
__swig_setmethods__["itermax"] = _likelihood.Log10Likelihood_itermax_set
__swig_getmethods__["itermax"] = _likelihood.Log10Likelihood_itermax_get
if _newclass:itermax = _swig_property(_likelihood.Log10Likelihood_itermax_get, _likelihood.Log10Likelihood_itermax_set)
__swig_setmethods__["lower_bound"] = _likelihood.Log10Likelihood_lower_bound_set
__swig_getmethods__["lower_bound"] = _likelihood.Log10Likelihood_lower_bound_get
if _newclass:lower_bound = _swig_property(_likelihood.Log10Likelihood_lower_bound_get, _likelihood.Log10Likelihood_lower_bound_set)
__swig_setmethods__["upper_bound"] = _likelihood.Log10Likelihood_upper_bound_set
__swig_getmethods__["upper_bound"] = _likelihood.Log10Likelihood_upper_bound_get
if _newclass:upper_bound = _swig_property(_likelihood.Log10Likelihood_upper_bound_get, _likelihood.Log10Likelihood_upper_bound_set)
def __init__( self, bursts, nopen, tau, tcritical=None, nmax=3, xtol=1e-10, rtol=1e-10,
itermax=100, lower_bound=None, upper_bound=None ):
""" Initializes the missed-events likelihood.
:param bursts:
A sequence of bursts. Each burst is a sequence of open and shut intervals.
There should always be an odd number of intervals. And the first interval
is open by default.
:param integer nopen:
Number of open states. Should be strictly positive.
:param float tau:
Maximum length of missed events. Should be positive.
:param float tcritical:
:math:`t_{\mathrm{crit}}` for CHS vectors. If None, then the initial and final states
are the equilibrium occupancies.
:param integer nmax:
The exact missed-events likelihood will be computed for
:math:`t < n_{\matrm{max}}\tau`.
:param float xtol:
Tolerance criteria when computing roots via :meth:`likelihood.brentq`
:param float rtol:
Tolerance criteria when computing roots via :meth:`likelihood.brentq`.
:param interger itermax:
Maximum number of iterations when computing roots via :meth:`likelihood.brentq`.
:param float lower_bound:
Lower bound for *all* roots. Defaults to None, in which the case the lower bound is
computed from :py:func:`find_lower_bound_for_roots`.
:param float upper_bound:
Upper bound for *all* roots. Defaults to None, in which the case the upper bound is
computed from :py:func:`find_upper_bound_for_roots`.
"""
from numpy import NaN
if lower_bound is None: lower_bound = NaN
if upper_bound is None: upper_bound = NaN
if tcritical is None: tcritical = -1e0
this = _likelihood.new_Log10Likelihood( bursts, nopen, tau, tcritical, nmax, xtol, rtol, itermax, lower_bound,
upper_bound )
try: self.this.append(this)
except: self.this = this
def __getitem__(self, *args):
"""
Returns the nth burst.
"""
return _likelihood.Log10Likelihood___getitem__(self, *args)
[docs] def insert(self, *args):
"""
Inserts a burst at given index.
"""
return _likelihood.Log10Likelihood_insert(self, *args)
def __setitem__(self, *args):
"""
Sets the nth burst.
:param index:
Index to the burst that should be changed.
:param value:
A list of time intervals with an odd number of components, in seconds.
"""
return _likelihood.Log10Likelihood___setitem__(self, *args)
def __delitem__(self, *args):
"""
Deletes the nth burst.
"""
return _likelihood.Log10Likelihood___delitem__(self, *args)
def __len__(self):
"""
Number of bursts.
"""
return _likelihood.Log10Likelihood___len__(self)
def append(self, *args): return _likelihood.Log10Likelihood_append(self, *args)
def __str__(self): return _likelihood.Log10Likelihood___str__(self)
def __repr__(self): return _likelihood.Log10Likelihood___repr__(self)
[docs] def vector(self, *args):
"""
vector(self, _Q) -> DCProgs::t_rvector
Computes the :math:`log_{10}`-likelihood for each individual burst
The sum over the return of this function is the return of :py:meth:`__call__`.
This is merely a convenience function if one wants the likelihood associate with
each burst
:param matrix:
Can be either a square numpy array or a :py:class:`QMatrix` object. In the first
case, the number of open states is set to attribute `nopen`.
:return:
A numpy array where each component contains the :math:`log_{10}` likelihood of a
single burst.
"""
return _likelihood.Log10Likelihood_vector(self, *args)
[docs] def __call__(self, *args):
"""
__call__(self, _Q) -> DCProgs::t_real
Computes the :math:`log_{10}`-likelihood for the input Q-matrix
:param matrix:
Can be either a square numpy array or a :py:class:`QMatrix` object. In the first
case, the number of open states is set to attribute `nopen`.
"""
return _likelihood.Log10Likelihood___call__(self, *args)
__swig_destroy__ = _likelihood.delete_Log10Likelihood
__del__ = lambda self : None;
Log10Likelihood_swigregister = _likelihood.Log10Likelihood_swigregister
Log10Likelihood_swigregister(Log10Likelihood)
def _brentq(*args):
return _likelihood._brentq(*args)
_brentq = _likelihood._brentq
[docs]def brentq(function, a, b, xtol=1e-12, rtol=4.440892098e-16, itermax=100):
""" Computes zeros via brentq.
This is the exact same algorithm as scipy.optimize.brentq. Only, the errors and c types
have been changed to accomodate DCProgs and protect the innocents.
:param callable function:
Function for which to solve :math:`f(x) = 0`.
:param float xstart:
Beginning of the interval
:param float xend:
End of the interval
:param float xtol:
Tolerance for interval size.
:param float rtol:
Tolerance for interval size. The convergence criteria is an affine function of the
root:
:math:`x_{\\mathrm{tol}}+r_{\\mathrm{tol}} x_{\\mathrm{current}} = \\frac{|x_a-x_b|}{2}`.
:param int itermax:
maximum number of iterations.
:returns: the tuple (x, iterations, times the function was called)
"""
if xtol <= 0: raise ValueError("xtol must be strictly positive");
if rtol <= 0: raise ValueError("rtol must be strictly positive");
if itermax <= 0: raise ValueError("itermax must be strictly positive");
if b <= a: raise ValueError("Interval is not valid: {0} >= {1}.".format(a, b))
return _likelihood._brentq(function, a, b, xtol, rtol, itermax)
def time_filter(*args):
"""
Filters out time series for intervals smaller that :math:`\tau`
This procedure erases small intervals from the time series.
:param times:
A timer series to filter. It can also arrays of time series. In that case, The filters are applied independently to all but the innermost sequence.
:returns: A numpy array, or a list of numpy arrays.
"""
return _likelihood.time_filter(*args)
def interval_filter(*args):
"""
Filters out intervals smaller that :math:`\tau`
This procedure erases small intervals from the time series that can be created from the input.
:param times:
A timer series to filter. It can also arrays of time series. In that case, The filters are applied independently to all but the innermost sequence.
:returns: A numpy array, or a list of numpy arrays.
"""
return _likelihood.interval_filter(*args)
def chained_likelihood(*args):
return _likelihood.chained_likelihood(*args)
chained_likelihood = _likelihood.chained_likelihood
# This file is compatible with both classic and new-style classes.