# 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)
                    gene-wide `omega`.

            Returns:
                `wr` (list)
                    list of `omega_r` values of length `nsites`

        Following
        `Spielman and Wilke,
        MBE, 32:1097-1108 <https://doi.org/10.1093/molbev/msv003>`_,
        we can predict the `dN/dS` value for each site `r`,
        :math:`\\rm{spielman}\\omega_r`, from the `ExpCM`.

        When `norm` is `False`, the `omega_r` values are defined as
        :math:`\\rm{spielman}\\omega_r = \\frac{\\sum_x\
        \\sum_{y \\in N_x}p_{r,x}P_{r,xy}}\
        {\\sum_x \\sum_{y \\in Nx}p_{r,x}Q_{xy}}`,
        where `r,x,y`, :math:`p_{r,x}`, :math:`P_{r,xy}`, and :math:`Q_{x,y}`
        have the same definitions as in the main `ExpCM` doc string
        and :math:`N_{x}`is the set of codons which are non-synonymous to
        codon `x` and differ from `x` by one nucleotide.

        When `norm` is `True`, the `omega_r` values above are divided by the
        ExpCM `omega` value.
        """
        wr = []
        for r in range(self.nsites):
            num = 0
            den = 0
            for i in range(N_CODON):
                # indices where single mut and non syn true

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

                    gene-wide `omega`.

            Returns:
                `wr` (list)
                    list of `omega_r` values of length `nsites`

        Following
        `Spielman and Wilke,
        MBE, 32:1097-1108 <https://doi.org/10.1093/molbev/msv003>`_,
        we can predict the `dN/dS` value for each site `r`,
        :math:`\\rm{spielman}\\omega_r`, from the `ExpCM`.

        When `norm` is `False`, the `omega_r` values are defined as
        :math:`\\rm{spielman}\\omega_r = \\frac{\\sum_x\
        \\sum_{y \\in N_x}p_{r,x}P_{r,xy}}\
        {\\sum_x \\sum_{y \\in Nx}p_{r,x}Q_{xy}}`,
        where `r,x,y`, :math:`p_{r,x}`, :math:`P_{r,xy}`, and :math:`Q_{x,y}`
        have the same definitions as in the main `ExpCM` doc string
        and :math:`N_{x}`is the set of codons which are non-synonymous to
        codon `x` and differ from `x` by one nucleotide.

        When `norm` is `True`, the `omega_r` values above are divided by the
        ExpCM `omega` value.
        """
        wr = []
        for r in range(self.nsites):
            num = 0
            den = 0
            for i in range(N_CODON):
                # indices where single mut and non syn true

                j = numpy.where((CODON_SINGLEMUT[i] * CODON_NONSYN[i]) == 1)[0]