Added data for simulations, output trait matrices

data/matrices/precision3x3Full.tsv

@@ -0,0 +1,2 @@
+c_00	c_11	c_22	c_01	c_02	c_12
+0.5	2.0	1.0	0.0	0.5	0.0

data/matrices/precision3x3Primates.tsv

@@ -0,0 +1,2 @@
+c_00	c_11	c_22	c_01	c_02	c_12
+0.5	2.0	1.0	0.1	-0.45	0.3

data/trees/placnr.lht.rootedtree.nhx

@@ -0,0 +1 @@
+((((Trichechus_manatus_latirostris:0.035513,(Procavia_capensis:0.084637,Loxodonta_africana:0.041516)1:0.003159)1:0.012832,(Orycteropus_afer_afer:0.061819,(Elephantulus_edwardii:0.113286,Echinops_telfairi:0.132412)1:0.002801)1:0.003326)1:0.030446,(Dasypus_novemcinctus:0.062179,Choloepus_hoffmanni:0.052537)1:0.028553)1:0.0128868,((((Sorex_araneus:0.151181,Erinaceus_europaeus:0.132436)1:0.010489,Condylura_cristata:0.099251)1:0.014601,((((Rousettus_aegyptiacus:0.02426,Pteropus_vampyrus:0.017977)1:0.041808,(Rhinolophus_sinicus:0.033606,Hipposideros_armiger:0.034245)1:0.022424)1:0.005466,(Miniopterus_natalensis:0.048559,(Myotis_lucifugus:0.02061,Eptesicus_fuscus:0.020745)1:0.032909)1:0.026648)1:0.014977,(((Equus_caballus:0.040208,Ceratotherium_simum_simum:0.031375)1:0.01757,((Vicugna_pacos:0.01123,Camelus_bactrianus:0.007435)1:0.048955,(Sus_scrofa:0.061861,((Odocoileus_virginianus_texanus:0.022142,(Bos_taurus:0.015258,(Ovis_aries:0.007108,Capra_hircus:0.005201)1:0.01242)1:0.005885)1:0.046478,((Physeter_catodon:0.012098,(Lipotes_vexillifer:0.010043,((Tursiops_truncatus:0.004507,Orcinus_orca:0.002912)1:0.003737,Delphinapterus_leucas:0.005948)1:0.00282)1:0.004819)1:0.001789,Balaenoptera_acutorostrata_scammoni:0.011414)1:0.021555)1:0.008733)1:0.00488)1:0.02127)1:0.00171,((Canis_familiaris:0.045402,((Ursus_maritimus:0.009802,Ailuropoda_melanoleuca:0.01071)1:0.019944,(Mustela_putorius:0.014168,Enhydra_lutris_kenyoni:0.010979)1:0.033406)1:0.008825)1:0.009892,(Panthera_pardus:0.004889,(Acinonyx_jubatus:0.005956,Felis_catus:0.004223)1:0.001612)1:0.040037)1:0.028291)1:0.001352)1:0.00366)1:0.011854,((Tupaia_belangeri:0.09974,((Oryctolagus_cuniculus:0.057311,Ochotona_princeps:0.105655)1:0.054285,((Marmota_marmota_marmota:0.007619,Ictidomys_tridecemlineatus:0.00925)1:0.071266,((((Octodon_degus:0.060538,Chinchilla_lanigera:0.038902)1:0.006926,Cavia_porcellus:0.06173)1:0.012433,Heterocephalus_glaber:0.049792)1:0.058724,((Jaculus_jaculus:0.103095,(((Rattus_norvegicus:0.042005,Mus_musculus:0.03846)1:0.026492,Meriones_unguiculatus:0.06177)1:0.007886,(Peromyscus_maniculatus:0.043119,(Microtus_ochrogaster:0.051557,Mesocricetus_auratus:0.048603)1:0.002904)1:0.012733)1:0.068248)1:0.022736,(Dipodomys_ordii:0.109334,Castor_canadensis:0.067389)1:0.010586)1:0.009351)1:0.002793)1:0.013804)1:0.005752)1:0.001773,((Otolemur_garnettii:0.068772,(Propithecus_coquereli:0.020483,Microcebus_murinus:0.028823)1:0.018976)1:0.01813,(((Nomascus_leucogenys:0.011348,(Pongo_abelii:0.009387,((Pan_troglodytes:0.003287,Homo_sapiens:0.003197)1:0.000928,Gorilla_gorilla:0.004417)1:0.004317)1:0.00142)1:0.005095,((Chlorocebus_sabaeus:0.005997,(Macaca_mulatta:0.003914,(Papio_anubis:0.003289,(Mandrillus_leucophaeus:0.00288,Cercocebus_atys:0.003073)1:0.000625)1:0.000757)1:0.001692)1:0.002705,Colobus_angolensis:0.009421)1:0.010106)1:0.009753,((Saimiri_boliviensis:0.014318,Cebus_capucinus:0.012274)1:0.001913,(Callithrix_jacchus:0.017752,Aotus_nancymaae:0.012014)1:0.000671)1:0.019849)1:0.037902)1:0.009114)1:0.01138)1:0.00067825);

data/trees/primates.nwk

@@ -0,0 +1 @@
+(((((((Perodicticus_potto:20.3881,Arctocebus_calabarensis:20.3881)PerodicArctoce:16.8815,(Loris_tardigradus:26.2208,Nycticebus_coucang:26.2208)Loris_tNyctice:11.0488)PeroArctLoriNyct:0.517898,(Otolemur_crassicaudatus:18.6,Galago_senegalensis:18.6)OtolemuGalago_:19.1875)PerArcLorNycOtoGal:21.5369,(Daubentonia_madagascariensis:54.6986,(((Avahi_laniger:25.329,Propithecus_verreauxi:25.329)Avahi_lPropith:12.4536,(Cheirogaleus_medius:35.3,(Microcebus_murinus:18.4471,Mirza_coquereli:18.4471)MicroceMirza_c:16.8529)CheirMicroMirza:2.48264)AvaProCheMicMir:0.581348,(Varecia_variegata:25.4576,((Hapalemur_griseus:13.7458,Lemur_catta:13.7458)HapalemLemur_c:5.76808,Eulemur_rufus:19.5138)HapalLemurEulem:5.94372)VareHapaLemuEule:12.9064)AvPrChMiMiVaHaLeEu:16.3346)DaAvPrChMiMiVaHaLeEu:4.62574)PALNOGDAPCMMVHLE:14.5125,(((((Pithecia_pithecia:11.3305,(Cacajao_calvus:7.31,Chiropotes_chiropotes:7.31)CacajaoChiropo:4.02055)PitheCacajChiro:6.66267,Plecturocebus_donacophilus:17.9932)PithCacaChirPlec:3.88632,((((Saguinus_fuscicollis:13.91,((Callimico_goeldii:10.9036,Callithrix_jacchus:10.9036)CallimiCallith:1.5529,Leontopithecus_rosalia:12.4565)CalliCalliLeont:1.45346)SaguCallCallLeon:4.47386,Aotus_azarai:18.3839)SagCalCalLeoAot:1.2972,(Saimiri_oerstedii:16.0705,Sapajus_apella:16.0705)SaimiriSapajus:3.61059)SaCaCaLeAoSaSa:0.98254,(Alouatta_palliata:14.8903,(Ateles_belzebuth:12.1315,(Lagothrix_cana:9.78613,Brachyteles_arachnoides:9.78613)LagothrBrachyt:2.3454)AteleLagotBrach:2.7588)AlouAtelLagoBrac:5.77326)SaCaCaLeAoSaSaAlAtLaBr:1.21594)PCCPSCCLASSAALB:21.2718,(((((Mandrillus_sphinx:4.59037,Cercocebus_torquatus:4.59037)MandrilCercoce:7.80963,Macaca_mulatta:12.4)MandrCercoMacac:1.34957,(((Erythrocebus_patas:5.88567,Chlorocebus_aethiops:5.88567)ErythroChloroc:6.47182,(Cercopithecus_albogularis:11.3,Miopithecus_ogouensis:11.3)CercopiMiopith:1.05748)ErytChloCercMiop:0.590152,Allenopithecus_nigroviridis:12.9476)EryChlCerMioAll:0.801935)MaCeMaErChCeMiAl:5.67273,(((((Rhinopithecus_brelichi:7.16132,Pygathrix_cinerea:7.16132)RhinopiPygathr:0.575314,Nasalis_larvatus:7.73663)RhinoPygatNasal:2.23843,Semnopithecus_entellus:9.97506)RhinPygaNasaSemn:0.137949,(Presbytis_melalophos:8.82307,Trachypithecus_francoisi:8.82307)PresbytTrachyp:1.28994)RhiPygNasSemPreTra:3.90969,(Piliocolobus_badius:12.8,Colobus_angolensis:12.8)PiliocoColobus:1.2227)RhPyNaSePrTrPiCo:5.39961)MCMECCMARPNSPTPC:10.0192,(((Hylobates_agilis:8.129,Symphalangus_syndactylus:8.129)HylobatSymphal:2.14664,Nomascus_siki:10.2756)HylobSymphNomas:9.91357,((Gorilla_gorilla:9.0631,(Pan_paniscus:6.6509,Homo_sapiens:6.6509)Pan_panHomo_sa:2.41219)GorilPan_pHomo_:6.69907,Pongo_pygmaeus:15.7622)GoriPan_HomoPong:4.42705)HySyNoGoPaHoPo:9.25233)MCMECCMARPNSPTPCHSNGPHP:13.7097)PCCPSCCLASSAALBMCMECCMARPNSPTPCHSNGPHP:23.9098,Carlito_syrichta:67.0611)PCCPSCCLASSAALBMCMECCMARPNSPTPCHSNGPHPC:6.77584)PALNOGDAPCMMVHLEPCCPSCCLASSAALBMCMECCMARPNSPTPCHSNGPHPC:1.86153,(Cynocephalus_volans:18.626,Galeopterus_variegatus:18.626)CynocepGaleopt:57.0725)PALNOGDAPCMMVHLEPCCPSCCLASSAALBMCMECCMARPNSPTPCHSNGPHPCCG:6.40378,Tupaia_minor:82.1022);

data/trees/primates.short.nwk

@@ -0,0 +1 @@
+(((((Perodicticus_potto:20.3881,Arctocebus_calabarensis:20.3881)PerodicArctoce:16.8815,(Loris_tardigradus:26.2208,Nycticebus_coucang:26.2208)Loris_tNyctice:11.0488)PeroArctLoriNyct:0.517898,(Otolemur_crassicaudatus:18.6,Galago_senegalensis:18.6)OtolemuGalago_:19.1875)PerArcLorNycOtoGal:21.5369,(Daubentonia_madagascariensis:54.6986,(((Cheirogaleus_medius:35.3,Microcebus_murinus:35.3)CheirogMicroce:2.48264,Propithecus_verreauxi:37.7826)CheirMicroPropi:0.581348,(Varecia_variegata:25.4576,(Hapalemur_griseus:13.7458,Lemur_catta:13.7458)HapalemLemur_c:11.7118)VarecHapalLemur:12.9064)CheMicProVarHapLem:16.3346)DaChMiPrVaHaLe:4.62574)PALNOGDCMPVHL:14.5125,(((((Saguinus_fuscicollis:13.91,((Callimico_goeldii:10.9036,Callithrix_jacchus:10.9036)CallimiCallith:1.5529,Leontopithecus_rosalia:12.4565)CalliCalliLeont:1.45346)SaguCallCallLeon:5.77106,Sapajus_apella:19.6811)SagCalCalLeoSap:0.98254,(Alouatta_palliata:14.8903,Ateles_belzebuth:14.8903)AlouattAteles_:5.77326)SaCaCaLeSaAlAt:1.21594,(Pithecia_pithecia:11.3305,Cacajao_calvus:11.3305)PitheciCacajao:10.549)SaCaCaLeSaAlAtPiCa:21.2718,(((((Mandrillus_sphinx:4.59037,Cercocebus_torquatus:4.59037)MandrilCercoce:7.80963,Macaca_mulatta:12.4)MandrCercoMacac:1.34957,(Erythrocebus_patas:5.88567,Chlorocebus_aethiops:5.88567)ErythroChloroc:7.86391)ManCerMacEryChl:5.67273,((Semnopithecus_entellus:9.97506,Nasalis_larvatus:9.97506)SemnopiNasalis:0.137949,Trachypithecus_francoisi:10.113)SemnoNasalTrach:9.3093)MaCeMaErChSeNaTr:10.0192,((Gorilla_gorilla:9.0631,(Pan_paniscus:6.6509,Homo_sapiens:6.6509)Pan_panHomo_sa:2.41219)GorilPan_pHomo_:6.69907,Pongo_pygmaeus:15.7622)GoriPan_HomoPong:13.6794)MaCeMaErChSeNaTrGoPaHoPo:13.7097)SCCLSAAPCMCMECSNTGPHP:30.6856);

src/lib/matrices.hpp

@@ -342,9 +342,8 @@ class CorrelationMatrix : public Matrix3x3 {
                 assert(v >= 0.0);
             }
         }
-        std::cout << "The input precision matrix is:\n" << precision_matrix << std::endl;
-        Matrix3x3 correlation_matrix = precision_matrix.inverse();
-        assert(correlation_matrix.transpose() == correlation_matrix);
+        Matrix3x3 covariance_matrix = precision_matrix.inverse();
+        assert(covariance_matrix.transpose() == covariance_matrix);
         for (int i = 0; i < dimensions; i++) {
             if (fix_pop_size and i == dim_population_size) { continue; }
             if (fix_mut_rate and i == dim_mutation_rate_per_generation) { continue; }
@@ -353,11 +352,29 @@ class CorrelationMatrix : public Matrix3x3 {
                 if (fix_pop_size and j == dim_population_size) { continue; }
                 if (fix_mut_rate and j == dim_mutation_rate_per_generation) { continue; }
                 if (fix_gen_time and j == dim_generation_time) { continue; }
-                (*this)(i, j) = correlation_matrix(i, j);
+                (*this)(i, j) = covariance_matrix(i, j);
             }
         }
-        std::cout << "The correlation (taking fixed effect into account) matrix is:\n"
-                  << *this << std::endl;
+        auto part_corr_matrix = precision_matrix;
+        auto corr_matrix = covariance_matrix;
+        for (int i = 0; i < dimensions; i++) {
+            corr_matrix.row(i) /= std::sqrt(covariance_matrix(i, i));
+            corr_matrix.col(i) /= std::sqrt(covariance_matrix(i, i));
+            part_corr_matrix.row(i) /= -std::sqrt(precision_matrix(i, i));
+            part_corr_matrix.col(i) /= std::sqrt(precision_matrix(i, i));
+        }
+        std::cout << "Covariances (taking fixed effect into account):\n" << *this << std::endl;
+        std::cout << "Correlation coefficients (taking fixed effect into account):\n"
+                  << corr_matrix << std::endl;
+        if (std::abs(corr_matrix.maxCoeff()) > 1.01) {
+            std::cerr
+                << "The precision matrix is not properly defined, specifically the "
+                   "correlation matrix has at least 1 entry greater than one (in absolute value)."
+                << std::endl;
+            exit(1);
+        }
+        std::cout << "Precisions:\n" << precision_matrix << std::endl;
+        std::cout << "Partial correlation coefficients:\n" << part_corr_matrix << std::endl;
     }
 
     void add_to_trace(Trace &trace) const {

src/lib/process.hpp

@@ -653,6 +653,14 @@ class Sequence {
         tree.set_tag(node, "mutation_rate_per_generation",
             d_to_string(log_multivariate.mutation_rate_per_generation()));
 
+        if (!tree.is_leaf(node)) {
+            tracer_fossils.add("NodeName", node_name);
+            double age = tree.max_distance_to_root() - time_from_root;
+            tracer_fossils.add("Age", age);
+            tracer_fossils.add("LowerBound", age);
+            tracer_fossils.add("UpperBound", age);
+        }
+
         if (tree.is_root(node)) { return; }
         assert(parent != nullptr);
         double geom_pop_size = piecewise_multivariate.GeometricPopSize();
@@ -711,22 +719,15 @@ class Sequence {
         tracer_sequences.add("CodonSequence", get_dna_str());
         tracer_sequences.add("AASequence", get_aa_str());
 
-        if (!tree.is_leaf(node)) {
-            tracer_fossils.add("NodeName", node_name);
-            double age = tree.max_distance_to_root() - time_from_root;
-            tracer_fossils.add("Age", age);
-            tracer_fossils.add("LowerBound", age * 0.9);
-            tracer_fossils.add("UpperBound", age * 1.1);
-            return;
-        }
+        if (!tree.is_leaf(node)) { return; }
         // If the node is a leaf, output the DNA sequence and name.
         write_sequence(output_filename, node_name, this->get_dna_str());
 
         tracer_traits.add("TaxonName", node_name);
-        tracer_traits.add("LogGenerationTime", log_multivariate.log_generation_time());
         tracer_traits.add("LogPopulationSize", log_multivariate.log_population_size());
         tracer_traits.add(
             "LogMutationRatePerGeneration", log_multivariate.log_mutation_rate_per_generation());
+        tracer_traits.add("LogGenerationTime", log_multivariate.log_generation_time());
     }
 
     // Method returning the DNA std::string corresponding to the codon sequence.
@@ -794,7 +795,10 @@ class Process {
         run_recursive(tree.root(), output_filename);
         std::ofstream nhx;
         nhx.open(output_filename + ".nhx");
-        nhx << tree.as_string() << std::endl;
+        nhx << tree.as_string(true) << std::endl;
+        nhx.close();
+        nhx.open(output_filename + ".tree");
+        nhx << tree.as_string(false) << std::endl;
         nhx.close();
     }
 

src/lib/tree.hpp

@@ -121,28 +121,32 @@ class Tree {
 
     void set_tag(NodeIndex node, TagName tag, TagValue value) { tags_.at(node)[tag] = value; }
 
-    std::string recursive_string(NodeIndex node) const {
+    std::string recursive_string(NodeIndex node, bool output_tags) const {
         std::string newick;
 
         if (not children(node).empty()) {
             // It's an internal node
             newick += "(";
-            for (auto const child : children(node)) { newick += recursive_string(child) + ","; };
+            for (auto const child : children(node)) {
+                newick += recursive_string(child, output_tags) + ",";
+            };
             newick.pop_back();
             newick += ")";
         }
         newick += name_.at(node);
         newick += ":" + std::to_string(length_.at(node));
 
-        if (not tags_.at(node).empty()) {
+        if (output_tags and not tags_.at(node).empty()) {
             newick += "[&&NHX";
             for (auto& it : tags_.at(node)) { newick += ":" + it.first + "=" + it.second; }
             newick += "]";
         }
         return newick;
     }
 
-    std::string as_string() const { return recursive_string(root()) + "; "; }
+    std::string as_string(bool output_tags = true) const {
+        return recursive_string(root(), output_tags) + "; ";
+    }
 
     void add_to_trace(Trace& trace) {
         trace.add("#tree_nodes", nb_nodes());