Skip to content

Commit 93a6bd7

Browse files
4ment4ment
committed
simultron can simulate relaxed local clocks
1 parent bd61849 commit 93a6bd7

File tree

3 files changed

+102
-22
lines changed

3 files changed

+102
-22
lines changed

‎src/phyc/treeio.c

+23-8
Original file line numberDiff line numberDiff line change
@@ -60,18 +60,33 @@ void Tree_print_nexus_header_figtree( FILE *pf, Tree *tree ){
6060
}
6161

6262
// Taxlabels are in preorder
63+
void Tree_print_nexus_taxa_block( FILE *pf, Tree *tree ){
64+
65+
fprintf(pf, "Begin taxa;\n");
66+
fprintf(pf, "\tDimensions ntax=%d;\n", Tree_tip_count(tree) );
67+
fprintf(pf, "\tTaxlabels\n");
68+
69+
Node **nodes = get_tips( tree, PREORDER );
70+
71+
for ( int i = 0; i < Tree_tip_count(tree); i++) {
72+
if( Node_isleaf(nodes[i]) ) fprintf(pf, "\t\t%s\n", nodes[i]->name);
73+
}
74+
fprintf(pf, "\t\t;\nEnd;\n\n");
75+
free(nodes);
76+
}
77+
6378
void Tree_print_nexus_header_figtree_Taxa( FILE *pf, Tree *tree ){
64-
65-
fprintf(pf, "#NEXUS\n\nBegin taxa;\n");
66-
fprintf(pf, "\tDimensions ntax=%d;\n", Tree_tip_count(tree) );
67-
fprintf(pf, "\tTaxlabels\n");
68-
69-
Node **nodes = get_tips( tree, PREORDER );
70-
79+
80+
fprintf(pf, "#NEXUS\n\nBegin taxa;\n");
81+
fprintf(pf, "\tDimensions ntax=%d;\n", Tree_tip_count(tree) );
82+
fprintf(pf, "\tTaxlabels\n");
83+
84+
Node **nodes = get_tips( tree, PREORDER );
85+
7186
for ( int i = 0; i < Tree_tip_count(tree); i++) {
7287
if( Node_isleaf(nodes[i]) ) fprintf(pf, "\t\t%s\n", nodes[i]->name);
7388
}
74-
fprintf(pf, "\t\t;\nEnd;\n\n");
89+
fprintf(pf, "\t\t;\nEnd;\n\n");
7590
free(nodes);
7691
}
7792

‎src/phyc/treeio.h

+1
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,7 @@ void print_tree_extended_nexus_double( FILE *pf, Tree *tree, double *info );
9595
// Should not be used with trees with different topologies
9696
void Tree_print_nexus_header_figtree( FILE *pf, Tree *tree );
9797
void Tree_print_nexus_header_figtree_Taxa( FILE *pf, Tree *tree );
98+
void Tree_print_nexus_taxa_block( FILE *pf, Tree *tree );
9899
void Tree_print_nexus_header_figtree_BeginTrees( FILE *pf, Tree *tree );
99100
void Tree_print_nexus_with_annotation2( FILE *pf, Tree *tree, bool time );
100101
void Tree_print_nexus_with_annotation( FILE *pf, Tree *tree );

‎src/simultron.c

+78-14
Original file line numberDiff line numberDiff line change
@@ -264,6 +264,7 @@ int main( int argc, char *argv[] ){
264264
}
265265

266266
model_string = String_clone(model_type);
267+
free(model_type);
267268

268269
if( array_of_string_contains(model_string, nucleotide_models, sizeof(nucleotide_models) / sizeof(nucleotide_models[0])) ){
269270
dataType = new_NucleotideDataType();
@@ -334,6 +335,7 @@ int main( int argc, char *argv[] ){
334335
int rateCount = 0;
335336
double *rates = args_get_double_array(argc, argv, "-r", ',', &rateCount);
336337
SubstitutionModel_set_rates(m,rates);
338+
free(rates);
337339
}
338340

339341
}
@@ -519,19 +521,79 @@ int main( int argc, char *argv[] ){
519521
StringBuffer_append_format(buffer, "\nRates from tree file\n");
520522
}
521523
else {
522-
if ( strcasecmp(rate_dist, "logn") == 0 ) {
523-
bool success = false;
524-
double mean = args_get_double(argc, argv, "-M", &success );
525-
if( !success || mean <= 0 ) error("Could not read the mean of the lognormal distribution of categories [-M]");
524+
if ( strcasecmp(rate_dist, "logn") == 0 || strcasecmp(rate_dist, "exp") == 0 ) {
525+
526+
char* mean_str = args_get_string(argc, argv, "-M");
527+
if( mean_str == NULL ) error("Could not read the mean of the lognormal distribution of categories [-M]");
526528

527-
double logsigma = args_get_double(argc, argv, "-S", &success );
528-
if( !success || logsigma <= 0 ) error("Could not read the standard deviation of the distribution of categories [-S]");
529+
unsigned l = 0;
530+
double* means = String_to_double_array(mean_str, ',', &l);
531+
for (int i = 0; i < l; i++) {
532+
if( means[i] <= 0 ) error("Mean cannot be negative [-M]");
533+
}
534+
free(mean_str);
529535

530-
lognormal_discretize(log(mean), logsigma, rates, Tree_node_count(tree)-1);
531-
print_dvector(rates, Tree_node_count(tree)-1);
532-
randomize_dvector( rates, Tree_node_count(tree)-1 );
536+
double* logsigmas = NULL;
537+
if(strcasecmp(rate_dist, "logn") == 0){
538+
char* logsigma_str = args_get_string(argc, argv, "-S");
539+
if( logsigma_str == NULL ) error("Could not read the standard deviation of the distribution of categories [-S]");
540+
unsigned l2 = 0;
541+
logsigmas = String_to_double_array(logsigma_str, ',', &l2);
542+
for (int i = 0; i < l2; i++) {
543+
if( logsigmas[i] <= 0 ) error("logsigma cannot be negative [-M]");
544+
}
545+
free(logsigma_str);
546+
547+
if(l!=l2){
548+
error("-M != -S");
549+
}
550+
}
533551

534-
StringBuffer_append_format(buffer, "\nLognormal distribution mean = %f log(stdev) = %f\n", mean, logsigma);
552+
if(l > 1){
553+
int* counts = calloc(l, sizeof(int));
554+
for ( int i = 0; i < Tree_node_count(tree); i++ ) {
555+
if( !Node_isroot(nodes[i]) ){
556+
counts[Node_get_int_from_info(nodes[i], "class=")]++;
557+
}
558+
}
559+
double* temp_rates = dvector(Tree_node_count(tree)-1);
560+
561+
for (int i = 0; i < l; i++) {
562+
if(strcasecmp(rate_dist, "logn") == 0){
563+
lognormal_discretize(log(means[i]), logsigmas[i], temp_rates, counts[i]);
564+
StringBuffer_append_format(buffer, "\n%d Lognormal distribution mean = %f log(stdev) = %f", i, means[i], logsigmas[i]);
565+
}
566+
else{
567+
exponential_discretize(means[i], temp_rates, counts[i]);
568+
StringBuffer_append_format(buffer, "\n%d Exponential distribution mean = %f", i, means[i]);
569+
}
570+
randomize_dvector( temp_rates, counts[i]);
571+
int k = 0;
572+
for ( int j = 0; j < Tree_node_count(tree); j++ ) {
573+
int class = Node_get_int_from_info(nodes[j], "class=");
574+
if( !Node_isroot(nodes[j]) && class == i){
575+
rates[j] = temp_rates[k];
576+
k++;
577+
}
578+
}
579+
}
580+
StringBuffer_append_string(buffer, "\n");
581+
free(temp_rates);
582+
free(counts);
583+
}
584+
else{
585+
if(strcasecmp(rate_dist, "logn") == 0){
586+
lognormal_discretize(log(means[0]), logsigmas[0], rates, Tree_node_count(tree)-1);
587+
}
588+
else{
589+
exponential_discretize(means[0], rates, Tree_node_count(tree)-1);
590+
}
591+
randomize_dvector( rates, Tree_node_count(tree)-1);
592+
}
593+
free(means);
594+
if(logsigmas != NULL) {
595+
free(logsigmas);
596+
}
535597
}
536598
else if ( strcasecmp(rate_dist, "exp") == 0 ) {
537599
bool success = false;
@@ -546,7 +608,6 @@ int main( int argc, char *argv[] ){
546608
else {
547609
error("Coud not read the type of distribution -d exp or -d logn\n");
548610
}
549-
550611
StringBuffer *info = new_StringBuffer(100);
551612
for ( int i = 0; i < Tree_node_count(tree); i++ ) {
552613
Node_empty_annotation(nodes[i]);
@@ -556,7 +617,8 @@ int main( int argc, char *argv[] ){
556617
Node_set_annotation(nodes[i], "rate", info->c);
557618

558619
StringBuffer_empty(info);
559-
StringBuffer_append_format(info, "%d", i);
620+
int class = Node_get_int_from_info(nodes[i], "class=");
621+
StringBuffer_append_format(info, "%d", class);
560622
Node_set_annotation(nodes[i], "class", info->c);
561623
}
562624
}
@@ -658,8 +720,9 @@ int main( int argc, char *argv[] ){
658720

659721
if( outputtree != NULL ){
660722
FILE *testFile = fopen(outputtree,"w");
661-
662-
Tree_print_nexus_header_figtree_Taxa(testFile, tree);
723+
fprintf(testFile, "#NEXUS\n");
724+
fprintf(testFile, "%s\n\n", buffer->c);
725+
Tree_print_nexus_taxa_block(testFile, tree);
663726
Tree_print_nexus_header_figtree_BeginTrees(testFile, tree);
664727

665728
fprintf(testFile, "tree TREE = [&R] ");
@@ -693,6 +756,7 @@ int main( int argc, char *argv[] ){
693756
if ( freq_codon != NULL ) {
694757
free(freq_codon);
695758
}
759+
free(model_string);
696760
free_Sequences(sim);
697761
free(output);
698762
free_SiteModel(sm);

0 commit comments

Comments
 (0)