few changes

brnlenMCMC.py

@@ -169,7 +169,7 @@ def replaedglen():
 def mcmcbrn(postorder, patterns, rates):
     nodes = readnewick(treenewick())
 
-    run = 5
+    run = 2
     mcmc = 0
     output = os.path.join('brnlenMCMC_results.txt')
     newf = open(output, 'w')
@@ -194,34 +194,32 @@ def mcmcbrn(postorder, patterns, rates):
             print
 
     #         print i, nd.number
-            print 'current',
-            for j in nodes:
-                print j.edgelen,
+#             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)
-            prelike = prepareTree(nodes, patterns, rates)
-            preprior = (-nodes[i].edgelen/b)-(log(b))
-            pre_ln_posterior = prelike + preprior
+            currentlike = prepareTree(nodes, patterns, rates)
+            currentprior = (-nodes[i].edgelen/b)-(log(b))
+            current_ln_posterior = currentlike + currentprior
 
 
             u = random.random()
             m = exp(0.2*(u-0.5))
             nodes[i].edgelen = preedgelen*m
+            print 'proposed_node.edgelen=',nodes[i].number, nodes[i].edgelen 
 #             print 'new_edgelen=', nodes[i].edgelen
-            postprior = (-nodes[i].edgelen/b)-(log(b))
-            postlike = prepareTree(nodes, patterns, rates)
-            post_ln_posterior = postlike + preprior
+            proposedprior = (-nodes[i].edgelen/b)-(log(b))
+            proposedtlike = prepareTree(nodes, patterns, rates)
+            proposed_ln_posterior = proposedtlike + proposedprior
             print 
-#             print 'preprior=',  preprior
-            print 'post_ln_posterior=',   post_ln_posterior 
-            print 'pre_ln_posterior=', pre_ln_posterior
-#             print 'postlike=',  postlike
+            print 'current prior, currentlike=', currentprior, currentlike
+            print 'proposed prior, proposedlike=', proposedprior, proposedtlike
             print
             hastings_ratio = log(m)
             print 'hastings_ratio', hastings_ratio
-            logR = post_ln_posterior - pre_ln_posterior + hastings_ratio
+            logR = proposed_ln_posterior - current_ln_posterior + hastings_ratio
             print 'logR=', logR
             u2 = random.random()
             print 'log-u2=', log(u2)
@@ -231,11 +229,7 @@ def mcmcbrn(postorder, patterns, rates):
             else:
                 nodes[i].edgelen = preedgelen
                 print 'failed accept new proposal..'
-    #             print 'nodes[i]-b=============================', nodes[i].edgelen
-    #         for i in nodes:
-    #             print 'new', i.edgelen
-
-    #         print 'pre, post=', preprior, postprior 
+
             newf.write('%s\t'%(mcmc+1))
             print
             print 'proposed=',
@@ -250,28 +244,6 @@ def mcmcbrn(postorder, patterns, rates):
         print '------------------------------------------------------------------------------'
 
 
-#         b = 0.1
-#         preedgelen = node.edgelen
-
-
-#         preprior = (-(node.edgelen)/b)-(log(b))
-#         pre_ln_posterior = preloglike + log(preprior)
-#     
-#         u = random.random()
-#         m = exp(0.2*(u-0.5))
-#         newedgelen = node.edgelen*m
-#         postprior = (-(node.newedgelen)/b)-(log(b))
-#         post_ln_posterior = postloglike + log(postprior)
-#     
-#         hastings_ratio = log(m)
-#         logR = post_ln_posterior - pre_ln_posterior + hastings_ratio
-#     
-#         u2 = random.random()
-#         if log(u2) < logR:
-#             node.edgelen = newedgelen
-#         else:
-#             node.edgelen = node.edgelen
-# 
 
 
 def treenewick():
@@ -490,10 +462,10 @@ def calcExpectedHeight(num_species, mu_over_s):
 
 
 if __name__ == '__main__':
-#     random_seed = 146549 # 7632557, 12345
+    random_seed = 146549 # 7632557, 12345
     number_of_species = 5
     mutation_speciation_rate_ratio = 0.4 # 0.689655172  # yields tree height 1 for 6 species
-#     random.seed(random_seed)
+    random.seed(random_seed)
     species_tree_root = yuleTree(number_of_species, mutation_speciation_rate_ratio)
 #     print '#########'
 #     print species_tree_root

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, 'tree.tre')
+    path = os.path.join(script_dir, 'tree2.tre')
     with open(path, 'r') as content:
         newick = content.read()
     return newick