本文整理汇总了Python中MatrixUtil.assert_transition_matrix方法的典型用法代码示例。如果您正苦于以下问题:Python MatrixUtil.assert_transition_matrix方法的具体用法?Python MatrixUtil.assert_transition_matrix怎么用?Python MatrixUtil.assert_transition_matrix使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类MatrixUtil
的用法示例。
在下文中一共展示了MatrixUtil.assert_transition_matrix方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: get_response_content
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_response_content(fs):
# precompute some transition matrices
P_drift_selection = pgmsinglesite.create_drift_selection_transition_matrix(
fs.npop, fs.selection_ratio)
MatrixUtil.assert_transition_matrix(P_drift_selection)
P_mutation = pgmsinglesite.create_mutation_transition_matrix(
fs.npop, fs.mutation_ab, fs.mutation_ba)
MatrixUtil.assert_transition_matrix(P_mutation)
# define the R table headers
headers = ['generation', 'number.of.mutants']
# compute the path samples
P = np.dot(P_drift_selection, P_mutation)
mypath = PathSampler.sample_endpoint_conditioned_path(
fs.nmutants_initial, fs.nmutants_final, fs.ngenerations, P)
arr = [[i, nmutants] for i, nmutants in enumerate(mypath)]
# create the R table string and scripts
# get the R table
table_string = RUtil.get_table_string(arr, headers)
# get the R script
script = get_ggplot()
# create the R plot image
device_name = Form.g_imageformat_to_r_function[fs.imageformat]
retcode, r_out, r_err, image_data = RUtil.run_plotter(
table_string, script, device_name)
if retcode:
raise RUtil.RError(r_err)
return image_data
示例2: get_absorption_variance
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_absorption_variance(P, plain, absorbing):
"""
Get expected times to absorption.
Note that if an index is indicated as absorbing by its presence
in the sequence of absorbing state indices,
then it will be treated as absorbing
even if the transition matrix P indicates otherwise.
@param P: transition matrix
@param plain: sequence of plain state indices
@param absorbing: sequence of absorbing state indices
@return: variance of times to absorption or 0 from absorbing states
"""
# check that P is really a transition matrix
MatrixUtil.assert_transition_matrix(P)
# define some state lists
states = np.hstack((plain, absorbing))
# check that the index sequences match the size of P
if sorted(states) != range(len(P)):
raise ValueError('P is not conformant with the index sequences')
# compute the time to absorption
Q = P[plain, :][:, plain]
c = np.ones_like(plain)
I = np.eye(len(plain))
t = linalg.solve(I - Q, c)
# compute the variance
vplain = 2*linalg.solve(I - Q, t) - t*(t+1)
v = np.hstack((vplain, np.zeros_like(absorbing)))
return v[inverse_permutation(states)]
示例3: test_row_sums
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def test_row_sums(self):
N = 20
k = 4
mutation, fitness = get_test_mutation_fitness()
P = get_transition_matrix(N, k, mutation, fitness)
MatrixUtil.assert_transition_matrix(mutation)
MatrixUtil.assert_transition_matrix(P)
示例4: test_invariant_selection_transition
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def test_invariant_selection_transition(self):
selection = 1.1
nchromosomes = 3
npositions = 2
P = get_selection_transition_matrix(
selection, nchromosomes, npositions)
MatrixUtil.assert_transition_matrix(P)
示例5: test_mutation
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def test_mutation(self):
npop = 10
mutation_ab = 0.1
mutation_ba = 0.2
P = create_mutation_transition_matrix(
npop, mutation_ab, mutation_ba)
MatrixUtil.assert_transition_matrix(P)
示例6: get_type_2_info
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_type_2_info(P):
"""
The expected time for a type 2 event is computed as follows.
It is the expected number of steps from AB to ab
conditional on not entering the states AB, Ab, or aB.
It should also include a bit of exponential delay that it takes
to leave the final fixed AB state before embark.
@param P: a huge transition matrix which is not modified
@return: expectation and variance of compensatory substitution time
"""
MatrixUtil.assert_transition_matrix(P)
nstates = len(P)
# define index sequences
plain = range(4, nstates)
forbidden = [0, 1, 2]
target = [3]
#
H = hittingtime.get_conditional_transition_matrix(
P, plain, forbidden, target)
t = hittingtime.get_absorption_time(
H, plain+forbidden, target)
v = hittingtime.get_absorption_variance(
H, plain+forbidden, target)
#
t0 = t[0]
v0 = v[0]
# add a geometric rv that depends on probability of leaving fixed AB
p = 1 - P[0, 0]
t0 += (1 - p) / p
v0 += (1 - p) / (p*p)
#
return t0, v0
示例7: get_response_content
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_response_content(fs):
# precompute some transition matrices
P_drift_selection = pgmsinglesite.create_drift_selection_transition_matrix(
fs.npop, fs.selection_ratio)
MatrixUtil.assert_transition_matrix(P_drift_selection)
P_mutation = pgmsinglesite.create_mutation_transition_matrix(
fs.npop, fs.mutation_ab, fs.mutation_ba)
MatrixUtil.assert_transition_matrix(P_mutation)
# define the R table headers
headers = [
'generation',
'number.of.mutants',
'probability',
'log.prob',
]
# compute the transition matrix
P = np.dot(P_drift_selection, P_mutation)
# Compute the endpoint conditional probabilities for various states
# along the unobserved path.
nstates = fs.npop + 1
M = np.zeros((nstates, fs.ngenerations))
M[fs.nmutants_initial, 0] = 1.0
M[fs.nmutants_final, fs.ngenerations-1] = 1.0
for i in range(fs.ngenerations-2):
A_exponent = i + 1
B_exponent = fs.ngenerations - 1 - A_exponent
A = np.linalg.matrix_power(P, A_exponent)
B = np.linalg.matrix_power(P, B_exponent)
weights = np.zeros(nstates)
for k in range(nstates):
weights[k] = A[fs.nmutants_initial, k] * B[k, fs.nmutants_final]
weights /= np.sum(weights)
for k, p in enumerate(weights):
M[k, i+1] = p
arr = []
for g in range(fs.ngenerations):
for k in range(nstates):
p = M[k, g]
if p:
logp = math.log(p)
else:
logp = float('-inf')
row = [g, k, p, logp]
arr.append(row)
# create the R table string and scripts
# get the R table
table_string = RUtil.get_table_string(arr, headers)
# get the R script
script = get_ggplot()
# create the R plot image
device_name = Form.g_imageformat_to_r_function[fs.imageformat]
retcode, r_out, r_err, image_data = RUtil.run_plotter(
table_string, script, device_name)
if retcode:
raise RUtil.RError(r_err)
return image_data
示例8: test_invariant_mutation_transition_s
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def test_invariant_mutation_transition_s(self):
mutation = 0.01
nchromosomes = 3
npositions = 2
ci_to_short, short_to_count, sorted_chrom_lists = get_state_space_info(
nchromosomes, npositions)
P = get_mutation_transition_matrix_s(
ci_to_short, short_to_count, sorted_chrom_lists,
mutation, nchromosomes, npositions)
MatrixUtil.assert_transition_matrix(P)
示例9: sample_endpoint_conditioned_path
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def sample_endpoint_conditioned_path(
initial_state, final_state, path_length, P):
"""
Return a sequence of states.
The returned sequence starts at the initial state
and ends at the final state.
Consecutive states may be the same
if the transition matrix has positive diagonal elements.
This function is derived from an earlier function
in an old SamplePath module.
@param initial_state: the first state as a python integer
@param final_state: the last state as a python integer
@param path_length: the number of states in the returned path
@param P: ndarray state transition matrix
@return: list of integer states
"""
# get the size of the state space and do some input validation
MatrixUtil.assert_transition_matrix(P)
nstates = len(P)
if not (0 <= initial_state < nstates):
raise ValueError('invalid initial state')
if not (0 <= final_state < nstates):
raise ValueError('invalid final state')
# take care of edge cases
if path_length == 0:
return []
elif path_length == 1:
if initial_state != final_state:
raise ValueError('unequal states for a path of length one')
return [initial_state]
elif path_length == 2:
return [initial_state, final_state]
# create transition matrices raised to various powers
max_power = path_length - 2
matrix_powers = [1]
matrix_powers.append(P)
for i in range(2, max_power + 1):
matrix_powers.append(np.dot(matrix_powers[i-1], P))
# sample the path
path = [initial_state]
for i in range(1, path_length-1):
previous_state = path[i-1]
weight_state_pairs = []
for state_index in range(nstates):
weight = 1
weight *= P[previous_state, state_index]
weight *= matrix_powers[path_length - 1 - i][state_index, final_state]
weight_state_pairs.append((weight, state_index))
next_state = Util.weighted_choice(weight_state_pairs)
path.append(next_state)
path.append(final_state)
return path
示例10: get_sample_distn
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_sample_distn(M_large, M_small, mut_ratio, fit_ratio, h):
"""
@param M_large: an array of microstate counts
@param M_small: an array of microstate counts
@param mut_ratio: positive mutation rate ratio
@param fit_ratio: positive fitness ratio
@param h: recessiveness or dominance between 0 and 1
@return: a distribution over polymorphic samples
"""
# N is the population haploid count
# n is the sample count
N = sum(M_large[0])
n = sum(M_small[0])
nstates_large = len(M_large)
nstates_small = len(M_small)
if N-1 != nstates_large:
raise Exception('internal error unexpected population state space')
if N % 2:
raise Exception('expected an even number of haploids in the population')
if n-1 != nstates_small:
raise Exception('internal error unexpected sample state space')
# Get the value of s to construct the diallelic transition matrix.
# It is positive when the lowercase allele is less fit.
s = 1.0 - fit_ratio
# Get the diallelic transition matrix.
# The first state is allele A fixation
# and the last state is allele a fixation.
N_diploid = N // 2
P = np.exp(wfengine.create_diallelic_recessive(N_diploid, s, h))
MatrixUtil.assert_transition_matrix(P)
# Get the expected number of visits
# to polymorphic states, starting from (1, N-1) and from (N-1, 1),
# assuming that fixed states are absorbing.
I = np.eye(N - 1)
Q = P[1:-1, 1:-1]
B = np.zeros((N - 1, 2))
B[0,0] = 1
B[-1,-1] = 1
X = linalg.solve((I - Q).T, B)
# At this point X has two columns, each giving an array of expectations.
# The first column gives expectations starting from a single A allele.
# The second column gives expectations starting from a single a allele.
# To get the expectations defined by the mixture,
# take the weighted sum of the columns.
e_large = X[:,0] * mut_ratio + X[:,1]
# Compute a multivariate hypergeometric reduction
# from the population state space expectations
# to the sample state space expectations.
e_small = wfengine.reduce_hypergeometric(M_large, M_small, e_large)
# Compute the distribution.
return e_small / e_small.sum()
示例11: get_type_1_info
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_type_1_info(P):
"""
The expected time for a type 1 event is computed as follows.
It is the sum of two expected times.
The first time is the expected number of steps from AB to either Ab or aB
conditional on not entering the states AB or ab.
The second time is the expected number of steps from Ab to ab
conditional on not entering AB.
Note that this formulation depends on the assumption
that the process associated with the first
step is equally likely to end up in Ab as in aB.
It should also include a bit of exponential delay that it takes
to leave the final fixed AB state before embark.
@param P: a huge transition matrix which is not modified
@return: expectation and variance of compensatory substitution time
"""
MatrixUtil.assert_transition_matrix(P)
nstates = len(P)
# get the expected time for the first stage
plain = range(4, nstates)
forbidden = [0, 3]
target = [1, 2]
H = hittingtime.get_conditional_transition_matrix(
P, plain, forbidden, target)
t = hittingtime.get_absorption_time(
H, plain+forbidden, target)
v = hittingtime.get_absorption_variance(
H, plain+forbidden, target)
t_first = t[0]
v_first = v[0]
# get the expected time for the second stage
plain = [1, 2] + range(4, nstates)
forbidden = [0]
target = [3]
H = hittingtime.get_conditional_transition_matrix(
P, plain, forbidden, target)
t = hittingtime.get_absorption_time(
H, plain+forbidden, target)
v = hittingtime.get_absorption_variance(
H, plain+forbidden, target)
t_second = t[1]
v_second = v[1]
# add a geometric rv that depends on probability of leaving fixed AB
p = 1 - P[0, 0]
t_third = (1 - p) / p
v_third = (1 - p) / (p*p)
# return the moments of the distribution accounting for both stages
return t_first + t_second + t_third, v_first + v_second + v_third
示例12: get_fixation_conditioned_matrix
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_fixation_conditioned_matrix(P, x):
"""
UNFINISHED
@param P: unconditioned transition matrix
@param x: fixation probabilities
@return: a smaller transition matrix
"""
MatrixUtil.assert_transition_matrix(P)
N_haploid = len(P) - 1
B = A - np.eye(N_haploid)
b = -np.ones(N_haploid)
B[-1] = np.zeros(N_haploid)
B[-1, -1] = 1
b[-1] = 0
y = linalg.solve(B, b)
print y[0]
示例13: test_wright_fisher
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def test_wright_fisher(self):
N_diallelic = 3
s = 0.03
# compute the wright fisher transition matrix directly
Ma = np.exp(wfengine.create_genic_diallelic(N_diallelic, s))
MatrixUtil.assert_transition_matrix(Ma)
# compute the wright fisher transition matrix less directly
N = 2*N_diallelic
fi = 1.0
fj = 1.0 - s
log_distns = np.zeros((N+1, 2))
for h in range(0, N+1):
probs = genic_diallelic(fi, fj, h, (N-h))
log_distns[h] = np.log(np.array(probs))
Mb = np.exp(wfengine.expand_multinomials(N, log_distns))
MatrixUtil.assert_transition_matrix(Mb)
# compare the transition matrices
self.assertTrue(np.allclose(Ma, Mb))
示例14: get_conditional_transition_matrix
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_conditional_transition_matrix(P, plain, forbid, target):
"""
Get a conditional transition matrix.
The new transition matrix is conditioned on reaching a state
in the sequence of target states before reaching a state
in the sequence of forbidden states.
@param P: transition matrix
@param plain: sequence of plain state indices
@param forbid: sequence of forbidden state indices
@param target: sequence of target state indices
@return: a conditional transition matrix conformant with P
"""
# check that P is really a transition matrix
MatrixUtil.assert_transition_matrix(P)
# define some state lists
special = np.hstack((forbid, target))
states = np.hstack((plain, special))
# check that the index sequences match the size of P
if sorted(states) != range(len(P)):
raise ValueError('P is not conformant with the index sequences')
# Q is the part of the transition matrix that gives
# transition probabilities from transient to transient states.
# In general it is not row stochastic.
Q = P[plain, :][:, plain]
# R is the part of the transition matrix that gives
# transition probabilities from transient to absorbing states.
# In general it is not row stochastic.
R = P[plain, :][:, special]
# Note that Q and R completely define an absorbing process,
# because there are no transitions away from absorbing states.
# Each row of B is a probability distribution.
B = linalg.solve(np.eye(len(plain)) - Q, R)
# Use the distributions in B to compute weights
# to apply to the transitions.
c = np.hstack((np.zeros_like(forbid), np.ones_like(target)))
w = np.hstack((np.dot(B, c), c))
# Make a new transition matrix that is conditioned
# on hitting a target state before hitting a forbidden state.
# Permute the weights to conform to the original matrix indices.
# Rescale rows to restore row-stochasticity.
H = P * w[inverse_permutation(states)]
v = np.sum(H, axis=1)
H /= v[:, np.newaxis]
return H
示例15: get_sample_distn
# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import assert_transition_matrix [as 别名]
def get_sample_distn(N, n, mutation_rate, fitness_ratio, h):
"""
@param N: haploid pop size
@param n: allele sample size
@param mutation_rate: mutation rate
@param fitness_ratio: fitness ratio
@param h: dominance parameter 0.5 when additive
@return: a distribution over n+1 mono-/di-morphic sample states
"""
s = 1.0 - fitness_ratio
P = np.exp(wfengine.create_diallelic_recessive(N // 2, s, h))
MatrixUtil.assert_transition_matrix(P)
# allow mutation out of the fixed states
P[0, 0] = 1.0 - mutation_rate
P[0, 1] = mutation_rate
P[-1, -1] = 1.0 - mutation_rate
P[-1, -2] = mutation_rate
MatrixUtil.assert_transition_matrix(P)
# get the population stationary distribution
v_large = MatrixUtil.get_stationary_distribution(P)
# get the allele distribution
v_small = StatsUtil.subsample_pmf_without_replacement(v_large, n)
return v_small