# Path: phydmslib/constants.py
# N_CODON = len(CODON_TO_INDEX)
# 
# CODON_NONSYN = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_TRANSITION = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_SINGLEMUT = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# ALMOST_ZERO = 1.0e-6
# 
# CODON_NT_INDEX = numpy.full((3, N_CODON), -1, dtype='int')
# 
# STOP_CODON_TO_NT_INDICES = []
# 
# N_STOP = 0
# 
# N_NT = len(INDEX_TO_NT)
# 
# N_AA = len(INDEX_TO_AA)
# 
# INDEX_TO_NT = dict(enumerate('ACGT'))
# 
# CODON_NT_MUT = numpy.full((N_NT, N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_NT_COUNT = numpy.zeros((N_NT, N_CODON), dtype='int')
# 
# CODON_TO_AA = []
# 
# INDEX_TO_AA = dict(enumerate('ACDEFGHIKLMNPQRSTVWY'))
# 
# AA_TO_INDEX = {aa: i for (i, aa) in INDEX_TO_AA.items()}
# 
# CODON_NT = numpy.full((3, N_NT, N_CODON), False, dtype='bool')


import copy
import six
import abc
import warnings
import numpy
import scipy.misc
import scipy.optimize
import scipy.linalg
import scipy.stats
import scipy.special
    import doctest
from phydmslib.numutils import (broadcastMatrixMultiply,
                                broadcastMatrixVectorMultiply,
                                broadcastGetCols)
from phydmslib.constants import (N_CODON, CODON_NONSYN, CODON_TRANSITION,
                                 CODON_SINGLEMUT, ALMOST_ZERO, CODON_NT_INDEX,
                                 STOP_CODON_TO_NT_INDICES, N_STOP, N_NT, N_AA,
                                 INDEX_TO_NT, CODON_NT_MUT, CODON_NT_COUNT,
                                 CODON_TO_AA, INDEX_TO_AA, AA_TO_INDEX,
                                 CODON_NT)
        `dprx` (dict)
            Keyed by strings in `freeparams`, each value is `numpy.ndarray`
            of floats giving derivative of `prx` with respect to that param,
            or 0 if if `prx` does not depend on parameter. The shape of each
            array is `(nsites, N_CODON)` except for *eta*, for which it is
            `(N_NT - 1, nsites, N_CODON)` with the first index ranging over
            each element in `eta`. Equivalent to `dstationarystate[param]`.
        `D` (`numpy.ndarray` floats, shape `(nsites, N_CODON)`)
            Element `r` is diagonal of :math:`\mathbf{D_r}` diagonal matrix,
            :math:`\mathbf{P_r} = \mathbf{A_r}^{-1} \mathbf{D_r} \mathbf{A_r}`
        `A` (`numpy.ndarray` floats, shape `(nsites, N_CODON, N_CODON)`)
            Element r is :math:`\mathbf{A_r}` matrix.
        `Ainv` (`numpy.ndarray` floats, shape `(nsites, N_CODON, N_CODON)`)
            Element r is :math:`\mathbf{A_r}^{-1}` matrix.
        `B` (dict)
            Keyed by strings in `freeparams`, each value is `numpy.ndarray`
            of floats giving `B` matrix for that parameter. The shape
            is `(nsites, N_CODON, N_CODON)` except for *eta*, for which
            it is `(N_NT - 1, nsites, N_CODON, N_CODON)` with the first
            index ranging over each element in `eta`.
    """

    # class variables
    ALLOWEDPARAMS = ['kappa', 'omega', 'beta', 'eta', 'mu']
    _REPORTPARAMS = ['kappa', 'omega', 'beta', 'phi']
    _PARAMLIMITS = {'kappa': (0.01, 100.0),
                    'omega': (1.0e-5, 100.0),
                    'beta': (1.0e-5, 10),
                    'eta': (0.01, 0.99),
                    'phi': (0.001, 0.999),

# Path: phydmslib/constants.py
# N_CODON = len(CODON_TO_INDEX)
# 
# CODON_NONSYN = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_TRANSITION = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_SINGLEMUT = numpy.full((N_CODON, N_CODON), False, dtype='bool')
# 
# ALMOST_ZERO = 1.0e-6
# 
# CODON_NT_INDEX = numpy.full((3, N_CODON), -1, dtype='int')
# 
# STOP_CODON_TO_NT_INDICES = []
# 
# N_STOP = 0
# 
# N_NT = len(INDEX_TO_NT)
# 
# N_AA = len(INDEX_TO_AA)
# 
# INDEX_TO_NT = dict(enumerate('ACGT'))
# 
# CODON_NT_MUT = numpy.full((N_NT, N_CODON, N_CODON), False, dtype='bool')
# 
# CODON_NT_COUNT = numpy.zeros((N_NT, N_CODON), dtype='int')
# 
# CODON_TO_AA = []
# 
# INDEX_TO_AA = dict(enumerate('ACDEFGHIKLMNPQRSTVWY'))
# 
# AA_TO_INDEX = {aa: i for (i, aa) in INDEX_TO_AA.items()}
# 
# CODON_NT = numpy.full((3, N_NT, N_CODON), False, dtype='bool')


# Path: phydmslib/models.py
import copy
import six
import abc
import warnings
import numpy
import scipy.misc
import scipy.optimize
import scipy.linalg
import scipy.stats
import scipy.special
    import doctest
from phydmslib.numutils import (broadcastMatrixMultiply,
                                broadcastMatrixVectorMultiply,
                                broadcastGetCols)
from phydmslib.constants import (N_CODON, CODON_NONSYN, CODON_TRANSITION,
                                 CODON_SINGLEMUT, ALMOST_ZERO, CODON_NT_INDEX,
                                 STOP_CODON_TO_NT_INDICES, N_STOP, N_NT, N_AA,
                                 INDEX_TO_NT, CODON_NT_MUT, CODON_NT_COUNT,
                                 CODON_TO_AA, INDEX_TO_AA, AA_TO_INDEX,
                                 CODON_NT)

        `dprx` (dict)
            Keyed by strings in `freeparams`, each value is `numpy.ndarray`
            of floats giving derivative of `prx` with respect to that param,
            or 0 if if `prx` does not depend on parameter. The shape of each
            array is `(nsites, N_CODON)` except for *eta*, for which it is
            `(N_NT - 1, nsites, N_CODON)` with the first index ranging over
            each element in `eta`. Equivalent to `dstationarystate[param]`.
        `D` (`numpy.ndarray` floats, shape `(nsites, N_CODON)`)
            Element `r` is diagonal of :math:`\mathbf{D_r}` diagonal matrix,
            :math:`\mathbf{P_r} = \mathbf{A_r}^{-1} \mathbf{D_r} \mathbf{A_r}`
        `A` (`numpy.ndarray` floats, shape `(nsites, N_CODON, N_CODON)`)
            Element r is :math:`\mathbf{A_r}` matrix.
        `Ainv` (`numpy.ndarray` floats, shape `(nsites, N_CODON, N_CODON)`)
            Element r is :math:`\mathbf{A_r}^{-1}` matrix.
        `B` (dict)
            Keyed by strings in `freeparams`, each value is `numpy.ndarray`
            of floats giving `B` matrix for that parameter. The shape
            is `(nsites, N_CODON, N_CODON)` except for *eta*, for which
            it is `(N_NT - 1, nsites, N_CODON, N_CODON)` with the first
            index ranging over each element in `eta`.
    """

    # class variables
    ALLOWEDPARAMS = ['kappa', 'omega', 'beta', 'eta', 'mu']
    _REPORTPARAMS = ['kappa', 'omega', 'beta', 'phi']
    _PARAMLIMITS = {'kappa': (0.01, 100.0),
                    'omega': (1.0e-5, 100.0),
                    'beta': (1.0e-5, 10),
                    'eta': (0.01, 0.99),
                    'phi': (0.001, 0.999),

                    'pi': (ALMOST_ZERO, 1),