7894856040

Committed 14 Feb 2024 12:52AM UTC coverage: 46.404% (-5.7%) from 52.085%

Build # 7894856040

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #38 from js51/SplitP-rewrite

Re-organise modules

Run Details

403 of 880 new or added lines in 12 files covered. (45.8%)

413 of 890 relevant lines covered (46.4%)

1.39 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

24.53

/splitp/model.py

import numpy as np
from scipy.linalg import expm
from . import constants

class model:
    def __init__(self, name=None, state_space=None, init_dist=None, rate_matrix=None):
        """A model of evolution."""
        # Set chosen model properties
        self.name = name if name else "Generic model"
        self.state_space = state_space if state_space else None
        self.init_dist = init_dist if init_dist else None
        self.rate_matrix = rate_matrix if rate_matrix else None
    
    def transition_matrix(self, t):
        """Return the transition matrix for the model."""
        return expm(t * self.rate_matrix)


class GTR(model):
    """
    A general time reversible model of evolution.
    """

    def __init__(
        self,
        name="Generic GTR model",
        state_space=None,
        equilibrium_distribution=None,
        additional_parameters=None
    ):
        self.state_space = state_space
        self.init_dist = equilibrium_distribution
        self.name = name
        self.rate_matrix = None
        n = len(state_space)
        
        ### Error checking
        # Check that we have the correct number of parameters
        if len(additional_parameters) != n * (n - 1) / 2 or len(equilibrium_distribution) != n:
            raise ValueError("Incorrect number of parameters for GTR model.")
        # Check that the equilibrium distribution sums to 1
        if sum(equilibrium_distribution) != 1:
            raise ValueError("Equilibrium distribution must sum to 1.")
        # Check that the parameters are positive
        if any([x < 0 for x in equilibrium_distribution + additional_parameters]):
            raise ValueError("All parameters must be positive.")
        
        ### Construct rate matrix
        # Construct matrix from parameters
        Q = np.zeros((n, n))
        Q[np.triu_indices(n, 1)] = additional_parameters
        Q += Q.T
        # Equilibrium parameters
        E = np.tile(equilibrium_distribution, (n, 1))
        E = E - np.diag(equilibrium_distribution)
        # Construict rate matrix
        Q = np.multiply(Q, E)
        Q -=  np.diag(np.sum(Q, axis=1))
        # Normalise rate matrix
        scale = -np.dot(equilibrium_distribution, np.diag(Q))
        self.rate_matrix = Q / scale

    def __str__(self) -> str:
        return self.name

    @classmethod
    def JukesCantor(cls, rate=1, state_space=constants.DNA_state_space, name="Jukes-Cantor model"):
        """ A Jukes-Cantor model of evolution. """
        n = len(state_space)
        return cls(
            name, 
            state_space, 
            [1/n for _ in range(n)], 
            [rate for _ in range(int(n*(n-1)/2))]
        )

    @classmethod
    def Kimura(cls, transversion_rate=1, transition_rate=1, state_space=constants.DNA_state_space, name="Kimura model"):
        """ A Kimura model of evolution. """
        n = len(state_space)
        # Transitions A-G and C-T
        temp = np.ones((n, n)) * transversion_rate
        A = state_space.index('A')
        G = state_space.index('G')
        C = state_space.index('C')
        T = state_space.index('T')
        temp[A, G] = transition_rate
        temp[G, A] = transition_rate
        temp[C, T] = transition_rate
        temp[T, C] = transition_rate
        parameters = list(temp[np.triu_indices(n, 1)])
        return cls(
            name, 
            state_space, 
            [1/n for _ in range(n)], 
            parameters
        )

1	import numpy as np	3✔
2	from scipy.linalg import expm	3✔
3	from . import constants	3✔
4
5	class model:	3✔
6	def __init__(self, name=None, state_space=None, init_dist=None, rate_matrix=None):	3✔
7	"""A model of evolution."""
8	# Set chosen model properties
NEW 9	self.name = name if name else "Generic model"	×
NEW 10	self.state_space = state_space if state_space else None	×
NEW 11	self.init_dist = init_dist if init_dist else None	×
NEW 12	self.rate_matrix = rate_matrix if rate_matrix else None	×
13
14	def transition_matrix(self, t):	3✔
15	"""Return the transition matrix for the model."""
NEW 16	return expm(t * self.rate_matrix)	×
17
18
19	class GTR(model):	3✔
20	"""
21	A general time reversible model of evolution.
22	"""
23
24	def __init__(	3✔
25	self,
26	name="Generic GTR model",
27	state_space=None,
28	equilibrium_distribution=None,
29	additional_parameters=None
30	):
NEW 31	self.state_space = state_space	×
NEW 32	self.init_dist = equilibrium_distribution	×
NEW 33	self.name = name	×
NEW 34	self.rate_matrix = None	×
NEW 35	n = len(state_space)	×
36
37	### Error checking
38	# Check that we have the correct number of parameters
NEW 39	if len(additional_parameters) != n * (n - 1) / 2 or len(equilibrium_distribution) != n:	×
NEW 40	raise ValueError("Incorrect number of parameters for GTR model.")	×
41	# Check that the equilibrium distribution sums to 1
NEW 42	if sum(equilibrium_distribution) != 1:	×
NEW 43	raise ValueError("Equilibrium distribution must sum to 1.")	×
44	# Check that the parameters are positive
NEW 45	if any([x < 0 for x in equilibrium_distribution + additional_parameters]):	×
NEW 46	raise ValueError("All parameters must be positive.")	×
47
48	### Construct rate matrix
49	# Construct matrix from parameters
NEW 50	Q = np.zeros((n, n))	×
NEW 51	Q[np.triu_indices(n, 1)] = additional_parameters	×
NEW 52	Q += Q.T	×
53	# Equilibrium parameters
NEW 54	E = np.tile(equilibrium_distribution, (n, 1))	×
NEW 55	E = E - np.diag(equilibrium_distribution)	×
56	# Construict rate matrix
NEW 57	Q = np.multiply(Q, E)	×
NEW 58	Q -= np.diag(np.sum(Q, axis=1))	×
59	# Normalise rate matrix
NEW 60	scale = -np.dot(equilibrium_distribution, np.diag(Q))	×
NEW 61	self.rate_matrix = Q / scale	×
62
63	def __str__(self) -> str:	3✔
NEW 64	return self.name	×
65
66	@classmethod	3✔
67	def JukesCantor(cls, rate=1, state_space=constants.DNA_state_space, name="Jukes-Cantor model"):	3✔
68	""" A Jukes-Cantor model of evolution. """
NEW 69	n = len(state_space)	×
NEW 70	return cls(	×
71	name,
72	state_space,
73	[1/n for _ in range(n)],
74	[rate for _ in range(int(n*(n-1)/2))]
75	)
76
77	@classmethod	3✔
78	def Kimura(cls, transversion_rate=1, transition_rate=1, state_space=constants.DNA_state_space, name="Kimura model"):	3✔
79	""" A Kimura model of evolution. """
NEW 80	n = len(state_space)	×
81	# Transitions A-G and C-T
NEW 82	temp = np.ones((n, n)) * transversion_rate	×
NEW 83	A = state_space.index('A')	×
NEW 84	G = state_space.index('G')	×
NEW 85	C = state_space.index('C')	×
NEW 86	T = state_space.index('T')	×
NEW 87	temp[A, G] = transition_rate	×
NEW 88	temp[G, A] = transition_rate	×
NEW 89	temp[C, T] = transition_rate	×
NEW 90	temp[T, C] = transition_rate	×
NEW 91	parameters = list(temp[np.triu_indices(n, 1)])	×
NEW 92	return cls(	×
93	name,
94	state_space,
95	[1/n for _ in range(n)],
96	parameters
97	)

js51 / SplitP / 7894856040

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous