Python MatrixUtil.dict_to_row_major方法代码示例

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
 def __call__(self, X_logs):
     """
     The vth entry of X corresponds to the log rate of the branch above v.
     Return the quantity to be minimized (the neg log likelihood).
     @param X: vector of branch rate logs
     @return: negative log likelihood
     """
     X = [math.exp(x) for x in X_logs]
     B_subs = {}
     for v_parent, v_child in self.R:
         edge = frozenset([v_parent, v_child])
         r = X[v_child]
         t = self.B[edge]
         B_subs[edge] = r * t
     newick_string = FtreeIO.RBN_to_newick(self.R, B_subs, self.N_leaves)
     tree = Newick.parse(newick_string, Newick.NewickTree)
     # define the rate matrix object; horrible
     dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix() 
     ordered_states = list('ACGT') 
     row_major_rate_matrix = MatrixUtil.dict_to_row_major(
             dictionary_rate_matrix, ordered_states, ordered_states)
     rate_matrix_object = RateMatrix.RateMatrix(
             row_major_rate_matrix, ordered_states) 
     # get the log likelihood
     ll = PhyLikelihood.get_log_likelihood(
             tree, self.alignment, rate_matrix_object)
     return -ll

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def demo_rate_matrix_stationary_distribution():
    # get a sample codon rate matrix that is in the form of a dictionary
    rate_matrix = get_sample_codon_rate_matrix()
    states = list(sorted(Codon.g_non_stop_codons))
    # convert the rate matrix to row major form
    row_major_matrix = MatrixUtil.dict_to_row_major(rate_matrix, states, states)
    # show the eigensystem of the row major matrix
    print 'row sums:'
    for row in row_major_matrix:
        print sum(row)
    # convert the transition matrix to numpy array form
    numpy_array = np.array(row_major_matrix)
    # get the eigensystem
    left = True
    right = False
    w, vl = linalg.eig(numpy_array, None, left, right)
    # find the eigenvalue with the smallest absolute value
    best_eigenvalue_size, best_index = min((abs(value), i) for i, value in enumerate(w))
    print 'smallest absolute value of an eigenvalue:', best_eigenvalue_size
    # get the stationary distribution
    scaled_stationary_distribution = list(vl.T[best_index])
    scaled_total = sum(scaled_stationary_distribution)
    stationary_distribution = [value / scaled_total for value in scaled_stationary_distribution]
    # assert the validity of the stationary distribution
    if min(stationary_distribution) < 0:
        raise RateMatrixError('found a negative element of the stationary eigenvector: %s' % min(stationary_distribution))
    print 'stationary distribution:'
    for value in stationary_distribution:
        print value

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_unscaled_hky85_rate_matrix(distribution, kappa):
    """
    @param distribution: a dictionary mapping a nucleotide to its stationary frequency
    @param kappa: the transition / transversion substitution rate ratio
    @return: a nucleotide rate matrix object
    """
    assert len(distribution) == 4
    assert set(distribution) == set('ACGT')
    assert abs(sum(distribution.values()) - 1.0) < .0000001
    # Create the off-diagonal elements of the unscaled rate matrix.
    rate_matrix = {}
    for na in distribution:
        for nb, probability in distribution.items():
            if na != nb:
                rate = probability
                if na+nb in ('AG', 'GA', 'CT', 'TC'):
                    rate *= kappa
                rate_matrix[(na, nb)] = rate
    # Create the diagonal elements such that each row in the rate matrix sums to zero.
    for na in distribution:
        rate = sum(rate_matrix[(na, nb)] for nb in distribution if nb != na)
        rate_matrix[(na, na)] = -rate
    # Convert the dictionary rate matrix to a row major rate matrix
    ordered_states = list('ACGT')
    row_major_rate_matrix = MatrixUtil.dict_to_row_major(rate_matrix, ordered_states, ordered_states)
    rate_matrix_object = RateMatrix(row_major_rate_matrix, ordered_states)
    return rate_matrix_object

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_response_content(fs):
    # get the tree
    tree = Newick.parse(fs.tree, Newick.NewickTree)
    tree.assert_valid()
    # get the alignment
    try:
        alignment = Fasta.Alignment(fs.fasta.splitlines())
        alignment.force_nucleotide()
    except Fasta.AlignmentError as e:
        raise HandlingError(e)
    # define the jukes cantor rate matrix
    dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix()
    ordered_states = list('ACGT')
    row_major_rate_matrix = MatrixUtil.dict_to_row_major(
            dictionary_rate_matrix, ordered_states, ordered_states)
    rate_matrix_object = RateMatrix.RateMatrix(
            row_major_rate_matrix, ordered_states)
    # simulate the ancestral alignment
    try:
        alignment = PhyLikelihood.simulate_ancestral_alignment(
                tree, alignment, rate_matrix_object)
    except PhyLikelihood.SimulationError as e:
        raise HandlingError(e)
    # get the alignment string using an ordering defined by the tree
    arr = []
    for node in tree.preorder():
        arr.append(alignment.get_fasta_sequence(node.name))
    # return the response
    return '\n'.join(arr) + '\n'

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def create_rate_matrix(distribution, kappa, f):
    """
    The parameter f does not affect the stationary distribution.
    @param distribution: a dictionary mapping a nucleotide to its frequency
    @param kappa: the transition / transversion substitution rate ratio
    @param f: a WAG-like parameter between zero and one
    @return: a nucleotide rate matrix object
    """
    assert len(distribution) == 4
    assert set(distribution) == set('ACGT')
    assert abs(sum(distribution.values()) - 1.0) < .0000001
    # Create the off-diagonal elements of the unscaled rate matrix.
    rate_matrix = {}
    for na, pa in distribution.items():
        for nb, pb in distribution.items():
            if na != nb:
                if f == 1:
                    rate = pb
                else:
                    rate = (pb**f) / (pa**(1-f))
                if na+nb in ('AG', 'GA', 'CT', 'TC'):
                    rate *= kappa
                rate_matrix[(na, nb)] = rate
    # Create the diagonal elements 
    # such that each row in the rate matrix sums to zero.
    for na in distribution:
        rate = sum(rate_matrix[(na, nb)] for nb in distribution if nb != na)
        rate_matrix[(na, na)] = -rate
    # Convert the dictionary rate matrix to a row major rate matrix
    ordered_states = list('ACGT')
    row_major_rate_matrix = MatrixUtil.dict_to_row_major(
            rate_matrix, ordered_states, ordered_states)
    rate_matrix_object = RateMatrix.RateMatrix(
            row_major_rate_matrix, ordered_states)
    return rate_matrix_object

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_form():
    """
    @return: the body of a form
    """
    # define the default transition matrix
    dictionary_rate_matrix = EnglishModel.get_transition_matrix()
    labels = list(sorted(set(a for a, b in dictionary_rate_matrix)))
    T = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, labels, labels)
    T = np.array(T)
    form_objects = [
            Form.Matrix('matrix', 'transition matrix',
                T, MatrixUtil.assert_transition_matrix)]
    return form_objects

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_form():
    """
    @return: the body of a form
    """
    # define the default rate matrix
    dictionary_rate_matrix = RateMatrix.get_sample_codon_rate_matrix()
    labels = list(sorted(set(a for a, b in dictionary_rate_matrix)))
    R = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, labels, labels)
    R = np.array(R)
    form_objects = [
            Form.Matrix('matrix', 'rate matrix',
                R, MatrixUtil.assert_rate_matrix)]
    return form_objects

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_response_content(fs):
    # get a properly formatted newick tree with branch lengths
    tree = Newick.parse(fs.tree, SpatialTree.SpatialTree)
    tree.assert_valid()
    if tree.has_negative_branch_lengths():
        msg = 'drawing a tree with negative branch lengths is not implemented'
        raise HandlingError(msg)
    tree.add_branch_lengths()
    # get the dictionary mapping the branch name to the nucleotide
    name_to_nucleotide = {}
    # parse the column string
    for line in iterutils.stripped_lines(fs.column.splitlines()):
        name_string, nucleotide_string = SnippetUtil.get_state_value_pair(line)
        if nucleotide_string not in list('acgtACGT'):
            msg = '"%s" is not a valid nucleotide' % nucleotide_string
            raise HandlingError(msg)
        nucleotide_string = nucleotide_string.upper()
        if name_string in name_to_nucleotide:
            raise HandlingError('the name "%s" was duplicated' % name_string)
        name_to_nucleotide[name_string] = nucleotide_string
    # augment the tips with the nucleotide letters
    for name, nucleotide in name_to_nucleotide.items():
        try:
            node = tree.get_unique_node(name)
        except Newick.NewickSearchError as e:
            raise HandlingError(e)
        if node.children:
            msg = 'constraints on internal nodes are not implemented'
            raise HandlingError(msg)
        node.state = nucleotide
    # get the Jukes-Cantor rate matrix object
    dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix()
    ordered_states = list('ACGT')
    row_major_rate_matrix = MatrixUtil.dict_to_row_major(
            dictionary_rate_matrix, ordered_states, ordered_states)
    rate_matrix_object = RateMatrix.RateMatrix(
            row_major_rate_matrix, ordered_states)
    # simulate the ancestral nucleotides
    rate_matrix_object.simulate_ancestral_states(tree)
    # simulate a path on each branch
    # this breaks up the branch into a linear sequence of nodes and adds color
    for node in tree.gen_non_root_nodes():
        simulate_branch_path(tree, node)
    # do the layout
    EqualArcLayout.do_layout(tree)
    # draw the image
    try:
        ext = Form.g_imageformat_to_ext[fs.imageformat]
        return DrawTreeImage.get_tree_image(tree, (640, 480), ext)
    except CairoUtil.CairoUtilError as e:
        raise HandlingError(e)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
 def test_simulation(self):
     tree_string = '(((Human:0.1, Chimpanzee:0.2)to-chimp:0.8, Gorilla:0.3)to-gorilla:0.7, Orangutan:0.4, Gibbon:0.5)all;'
     # Parse the example tree.
     tree = Newick.parse(tree_string, Newick.NewickTree)
     tree.assert_valid()
     # Get header and sequence pairs.
     alignment = Fasta.Alignment(StringIO(Fasta.brown_example_alignment))
     # Get the Jukes-Cantor rate matrix object.
     dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix()
     ordered_states = list('ACGT')
     row_major_rate_matrix = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, ordered_states, ordered_states)
     rate_matrix_object = RateMatrix.RateMatrix(row_major_rate_matrix, ordered_states)
     # Simulate ancestral states.
     simulated_alignment = simulate_ancestral_alignment(tree, alignment, rate_matrix_object)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_form():
    """
    @return: the body of a form
    """
    # define the default rate matrix
    dictionary_rate_matrix = RateMatrix.get_sample_codon_rate_matrix()
    labels = Codon.g_sorted_non_stop_codons
    R = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, labels, labels)
    # define the form objects
    form_objects = [
            Form.Matrix('matrix', 'codon rate matrix',
                R, MatrixUtil.assert_rate_matrix),
            Form.Integer('maxcategories', 'maximum number of categories',
                5, low=2)]
    return form_objects

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_form():
    """
    @return: the body of a form
    """
    # define the default matrix lines
    dictionary_rate_matrix = EnglishModel.get_transition_matrix()
    labels = list(sorted(set(a for a, b in dictionary_rate_matrix)))
    T = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, labels, labels)
    T = np.array(T)
    # define the form objects
    form_objects = [
            Form.Matrix('matrix', 'matrix', T),
            Form.Integer('maxcategories', 'maximum number of categories',
                5, low=2)]
    return form_objects

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
 def test_likelihood(self):
     # Parse the example tree.
     tree_string = Newick.brown_example_tree
     tree = Newick.parse(tree_string, Newick.NewickTree)
     tree.assert_valid()
     # Get header and sequence pairs.
     alignment = Fasta.Alignment(StringIO(Fasta.brown_example_alignment))
     # Get the Jukes-Cantor rate matrix object.
     dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix()
     ordered_states = list('ACGT')
     row_major_rate_matrix = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, ordered_states, ordered_states)
     rate_matrix_object = RateMatrix.RateMatrix(row_major_rate_matrix, ordered_states)
     # Calculate the log likelihood.
     log_likelihood = get_log_likelihood(tree, alignment, rate_matrix_object)
     self.assertAlmostEqual(log_likelihood, -4146.26547208)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
 def gen_distance_matrices(self, count, max_steps):
     """
     Yield (ordered sequence list, distance matrix) pairs .
     The generator will stop if it sees that it cannot meet its goal
     in the allotted number of steps.
     @param count: the requested number of distance matrices
     @param max_steps: an upper bound on the allowed number of steps
     """
     # define the jukes cantor rate matrix
     dictionary_rate_matrix = RateMatrix.get_jukes_cantor_rate_matrix()
     ordered_states = list('ACGT')
     row_major_rate_matrix = MatrixUtil.dict_to_row_major(
             dictionary_rate_matrix, ordered_states, ordered_states)
     model = RateMatrix.RateMatrix(row_major_rate_matrix, ordered_states)
     # record the requested number of samples
     self.requested_matrix_count = count
     # do some rejection sampling
     while True:
         if self.get_complexity() >= max_steps:
             break
         if self.accepted_sample_count >= count:
             break
         # simulate an alignment from the tree
         alignment = PhyLikelihood.simulate_alignment(
                 self.tree, model, self.sequence_length)
         # extract the ordered list of sequences from the alignment object
         name_to_sequence = dict(zip(alignment.headers, alignment.sequences))
         sequence_list = [name_to_sequence[name]
                 for name in self.ordered_names]
         # get the estimated distance matrix
         distance_matrix = JC69.get_ML_distance_matrix(sequence_list)
         # look for degeneracies
         has_zero_off_diagonal = False
         has_inf_off_diagonal = False
         for i, row in enumerate(distance_matrix):
             for j, value in enumerate(row):
                 if i != j:
                     if value == 0.0:
                         has_zero_off_diagonal = True
                     if value == float('inf'):
                         has_inf_off_diagonal = True
         if has_zero_off_diagonal:
             self.rejected_zero_sample_count += 1
         elif has_inf_off_diagonal:
             self.rejected_inf_sample_count += 1
         else:
             self.accepted_sample_count += 1
             yield sequence_list, distance_matrix

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
def get_form():
    """
    @return: the body of a form
    """
    # define the default transition matrix
    dictionary_rate_matrix = EnglishModel.get_transition_matrix()
    labels = list(sorted(set(a for a, b in dictionary_rate_matrix)))
    T = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, labels, labels)
    T = np.array(T)
    # define the form objects
    form_objects = [
        Form.Matrix("matrix", "transition matrix", T, MatrixUtil.assert_transition_matrix),
        Form.SingleLine("first", "first letter", "a"),
        Form.SingleLine("last", "last letter", "b"),
        Form.Integer("count", "character count", 10, low=1, high=80),
    ]
    return form_objects

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import dict_to_row_major [as 别名]
 def test_rate_matrix(self):
     """
     Compare the transition matrix found using a RateMatrix object to the transition matrix found analytically.
     """
     # get the row major Jukes-Cantor rate matrix
     dictionary_rate_matrix = get_jukes_cantor_rate_matrix()
     ordered_states = list('ACGT')
     row_major_rate_matrix = MatrixUtil.dict_to_row_major(dictionary_rate_matrix, ordered_states, ordered_states)
     # get the rate matrix object
     rate_matrix_object = RateMatrix(row_major_rate_matrix, ordered_states)
     # get the transition matrix
     t = 2.0
     calculated_transition_matrix = rate_matrix_object.get_dictionary_transition_matrix(t)
     analytical_transition_matrix = get_jukes_cantor_transition_matrix(t)
     self.assertEqual(len(calculated_transition_matrix), 16)
     self.assertEqual(len(analytical_transition_matrix), 16)
     for key, value in calculated_transition_matrix.items():
         self.assertAlmostEqual(analytical_transition_matrix[key], value)