cleared redudant parts

brnlenMCMC.py

@@ -169,7 +169,7 @@ def replaedglen():
 def mcmcbrn(postorder, patterns, rates):
     nodes = readnewick(treenewick())
 
-    run = 2
+    run = 1
     mcmc = 0
     output = os.path.join('brnlenMCMC_results.txt')
     newf = open(output, 'w')
@@ -189,19 +189,21 @@ def mcmcbrn(postorder, patterns, rates):
     for r in range(run):
         b = 0.1
 
-        for i,nd in enumerate(postorder):
+#         for i,nd in enumerate(postorder):
+        for i in range(len(postorder)):
+
             print 'mcmc gen=', mcmc+1,
             print
 
-    #         print i, nd.number
 #             print 'current',
             print 'current_node.edgelen=', nodes[i].number, nodes[i].edgelen
 #             print 
             preedgelen = nodes[i].edgelen
 #             print 'nodes[i]=', nodes[i].edgelen
-    #         prelike = prepareTree(postorder, patterns, rates)
             currentlike = prepareTree(nodes, patterns, rates)
+#             currentlike = 0.0
             currentprior = (-nodes[i].edgelen/b)-(log(b))
+#             newf.write('%s\t'%(currentprior))
             current_ln_posterior = currentlike + currentprior
 
 
@@ -211,7 +213,10 @@ def mcmcbrn(postorder, patterns, rates):
             print 'proposed_node.edgelen=',nodes[i].number, nodes[i].edgelen 
 #             print 'new_edgelen=', nodes[i].edgelen
             proposedprior = (-nodes[i].edgelen/b)-(log(b))
+#             newf.write('%s\t'%(proposedprior))
+
             proposedtlike = prepareTree(nodes, patterns, rates)
+#             proposedtlike = 0.0
             proposed_ln_posterior = proposedtlike + proposedprior
             print 
             print 'current prior, currentlike=', currentprior, currentlike
@@ -248,7 +253,7 @@ def mcmcbrn(postorder, patterns, rates):
 
 def treenewick():
     script_dir = os.path.dirname(os.path.realpath(sys.argv[0]))
-    path = os.path.join(script_dir, 'tree2.tre')
+    path = os.path.join(script_dir, 'tree.tre')
     with open(path, 'r') as content:
         newick = content.read()
     return newick
@@ -462,9 +467,9 @@ def calcExpectedHeight(num_species, mu_over_s):
 
 
 if __name__ == '__main__':
-    random_seed = 146549 # 7632557, 12345
+    random_seed = 346549 # 7632557, 12345
     number_of_species = 5
-    mutation_speciation_rate_ratio = 0.4 # 0.689655172  # yields tree height 1 for 6 species
+    mutation_speciation_rate_ratio = 0.689655172 # 0.689655172  # yields tree height 1 for 6 species
     random.seed(random_seed)
     species_tree_root = yuleTree(number_of_species, mutation_speciation_rate_ratio)
 #     print '#########'
@@ -479,11 +484,12 @@ def calcExpectedHeight(num_species, mu_over_s):
 #     print '  expected height                =',expected_height
     actual_height = calcActualHeight(species_tree_root)
 #     print '  actual height                  =',actual_height
-#     print '  newick: ',newick[0]
+    print '  newick: ',newick[0]
 
     alpha = 0.5 ### gamma shape parameter rate categories
 #     yuletree = '(((1:0.54019,(5:0.40299,10:0.40299):0.1372):0.72686,(6:0.10576,4:0.10576):1.16129):0.42537,(2:0.58122,(9:0.21295,(7:0.16691,(8:0.14622,3:0.14622):0.02069):0.04604):0.36827):1.1112)'
     rates_list = gammaRates(alpha)
     postorder = readnewick(treenewick())
     result =  prepareTree(postorder, readSeq.patterns(), rates_list)
+#     try1 = readSeq.patterns()
     result2 = mcmcbrn(postorder, readSeq.patterns(), rates_list)

readSeq.py

@@ -24,7 +24,7 @@ def patterns():
 
         m = re.search('nchar\s*=\s*(\d+)', f, re.M | re.S)
         nchar = int(m.group(1))
-#         print 'nchar=', nchar
+        print 'nchar=', nchar
 
         m = re.search('Matrix\s+(.+?);', f, re.M | re.S)
         matrix = m.group(1).strip()

simSequences.py

@@ -165,10 +165,10 @@ def readnewick(tree):
 
     output_filename = os.path.join('simulated_output.nexus')
 
-    yuletree = '(((1:0.54019,(5:0.40299,10:0.40299):0.13720):0.72686,(6:0.10576,4:0.10576):1.16129):0.42537,(2:0.58122,(9:0.21295,(7:0.16691,(8:0.14622,3:0.14622):0.02069):0.04604):0.36827):1.11120)'
+    yuletree = '(5:1.86010,((3:0.47109,2:0.47109):0.492,(4:0.05805,1:0.05805):0.906):0.896)'
 
     preorder = readnewick(yuletree)
-    ntax = 10
-    num_sites = 10000
+    ntax = 5
+    num_sites = 100
     result = simulate(preorder, ntax, num_sites, output_filename)
 

tree.tre

@@ -1 +1 @@
-(((1:0.00039,5:0.03915):0.987,(4:0.42253,2:0.000042):0.84):0.118,3:0.1433)
+(5:1.8601,((3:0.47109,2:0.47109):0.492,(4:0.05805,1:0.05805):0.906):0.896)

treeLike_JCGAMMA.py

@@ -160,7 +160,7 @@ def allocatePartial(node, patterns, rates):
 
 def treenewick():
     script_dir = os.path.dirname(os.path.realpath(sys.argv[0]))
-    path = os.path.join(script_dir, 'tree2.tre')
+    path = os.path.join(script_dir, 'treecopy.tre')
     with open(path, 'r') as content:
         newick = content.read()
     return newick
@@ -308,6 +308,8 @@ def readnewick(tree):
 
     post = pre[:]
     post.reverse()
+    for nd in post:
+        print nd.number, nd.edgelen
     return post
 
 def Makenewick(pre):