# 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)
        assert numpy.allclose(self.dphi_dbeta, dphi_dbeta_halfdx, atol=1e-5,
                              rtol=1e-4), ("The numerical derivative "
                                           "dphi_dbeta differs considerably "
                                           "in value for step dbeta = {0} and "
                                           "a step half that size, giving "
                                           "values of {1} and {2}."
                                           .format(dbeta, self.dphi_dbeta,
                                                   dphi_dbeta_halfdx))

    def _compute_empirical_phi(self, beta):
        """Returns empirical `phi` at the given value of `beta`.

        Does **not** set `phi` attribute, simply returns what
        should be value of `phi` given the current `g` and
        `pi_codon` attributes, plus the passed value of `beta`.
        Note that it uses the passed value of `beta`, **not**
        the current `beta` attribute.

        Initial guess is current value of `phi` attribute.
        """
        def F(phishort):
            """Difference between `g` and expected `g` given `phishort`."""
            phifull = numpy.append(phishort, 1 - phishort.sum())
            phiprod = numpy.ones(N_CODON, dtype='float')
            for w in range(N_NT):
                phiprod *= phifull[w]**CODON_NT_COUNT[w]
            frx_phiprod = frx * phiprod
            frx_phiprod_codonsum = frx_phiprod.sum(axis=1)
            gexpect = []
            for w in range(N_NT - 1):

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

        assert numpy.allclose(self.dphi_dbeta, dphi_dbeta_halfdx, atol=1e-5,
                              rtol=1e-4), ("The numerical derivative "
                                           "dphi_dbeta differs considerably "
                                           "in value for step dbeta = {0} and "
                                           "a step half that size, giving "
                                           "values of {1} and {2}."
                                           .format(dbeta, self.dphi_dbeta,
                                                   dphi_dbeta_halfdx))

    def _compute_empirical_phi(self, beta):
        """Returns empirical `phi` at the given value of `beta`.

        Does **not** set `phi` attribute, simply returns what
        should be value of `phi` given the current `g` and
        `pi_codon` attributes, plus the passed value of `beta`.
        Note that it uses the passed value of `beta`, **not**
        the current `beta` attribute.

        Initial guess is current value of `phi` attribute.
        """
        def F(phishort):
            """Difference between `g` and expected `g` given `phishort`."""
            phifull = numpy.append(phishort, 1 - phishort.sum())
            phiprod = numpy.ones(N_CODON, dtype='float')
            for w in range(N_NT):
                phiprod *= phifull[w]**CODON_NT_COUNT[w]
            frx_phiprod = frx * phiprod
            frx_phiprod_codonsum = frx_phiprod.sum(axis=1)
            gexpect = []
            for w in range(N_NT - 1):

                gexpect.append(((CODON_NT_COUNT[w] * frx_phiprod).sum(axis=1) /