[Avida-SVN] r2595 - development/source/analyze
kaben at myxo.css.msu.edu
kaben at myxo.css.msu.edu
Sun May 25 21:23:14 PDT 2008
Author: kaben
Date: 2008-05-26 00:23:14 -0400 (Mon, 26 May 2008)
New Revision: 2595
Modified:
development/source/analyze/cAnalyze.cc
development/source/analyze/cAnalyzeTreeStats_Gamma.cc
development/source/analyze/cAnalyzeTreeStats_Gamma.h
Log:
Added optional argument to PRINT_GAMMA analyze command, specifying
filename to store data for a lineage-through-time plot. PRINT_GAMMA has the form:
PRINT_GAMMA end-of-run gamma-filename lineage-through-time-datafile
The analyze.cfg file using this command should read something like this:
# Analyze config file
SET_BATCH 0
PURGE_BATCH
LOAD data/historic-1000.pop
LOAD data/detail-1000.pop
PRINT_GAMMA 1000 gamma.dat gamma-ltt.dat
Modified: development/source/analyze/cAnalyze.cc
===================================================================
--- development/source/analyze/cAnalyze.cc 2008-05-26 00:57:09 UTC (rev 2594)
+++ development/source/analyze/cAnalyze.cc 2008-05-26 04:23:14 UTC (rev 2595)
@@ -3299,7 +3299,24 @@
cString filename("gamma.dat");
if (cur_string.GetSize() != 0) filename = cur_string.PopWord();
- int furcation_time_convention = (cur_string.GetSize()) ? cur_string.PopWord().AsInt() : 1; // #1
+ cString lineage_thru_time_fname("");
+ if (cur_string.GetSize() != 0) lineage_thru_time_fname = cur_string.PopWord();
+
+ /*
+ I've hardwired the option 'furcation_time_convention' to '1'.
+
+ 'furcation_time_convention' refers to the time at which a 'speciation' event
+ occurs (I'm not sure 'speciation' is the right word for this). If a parent
+ genotype produces two distinct surviving lineages, then the time of
+ speciation could be:
+ - 1: The parent genotype's birth time
+ - 2: The elder child genotype's birth time
+ - 3: The younger child genotype's birth time
+
+ @kgn
+ */
+ // int furcation_time_convention = (cur_string.GetSize()) ? cur_string.PopWord().AsInt() : 1;
+ int furcation_time_convention = 1;
ofstream& fp = m_world->GetDataFileOFStream(filename);
@@ -3307,11 +3324,25 @@
fp << "# 1: Pybus-Harvey gamma statistic" << endl;
fp << endl;
- cAnalyzeTreeStats_Gamma agts(m_world);
- agts.AnalyzeBatch(batch[cur_batch].List(), end_time, furcation_time_convention);
+ cAnalyzeTreeStats_Gamma atsg(m_world);
+ atsg.AnalyzeBatch(batch[cur_batch].List(), end_time, furcation_time_convention);
- fp << agts.Gamma();
+ fp << atsg.Gamma();
fp << endl;
+
+ if(lineage_thru_time_fname != ""){
+ ofstream& ltt_fp = m_world->GetDataFileOFStream(lineage_thru_time_fname);
+
+ ltt_fp << "# Legend:" << endl;
+ ltt_fp << "# 1: num_lineages" << endl;
+ ltt_fp << "# 2: furcation_time" << endl;
+ ltt_fp << endl;
+
+ int size = atsg.FurcationTimes().GetSize();
+ for(int i = 0; i < size; i++){
+ ltt_fp << i+2 << " " << atsg.FurcationTimes()[i] << endl;
+ }
+ }
}
Modified: development/source/analyze/cAnalyzeTreeStats_Gamma.cc
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_Gamma.cc 2008-05-26 00:57:09 UTC (rev 2594)
+++ development/source/analyze/cAnalyzeTreeStats_Gamma.cc 2008-05-26 04:23:14 UTC (rev 2595)
@@ -194,16 +194,15 @@
int FurcationTimePolicy_FirstChildBirth(cAnalyzeLineageFurcation &furcation);
int FurcationTimePolicy_SecondChildBirth(cAnalyzeLineageFurcation &furcation);
*/
- tArray<cAnalyzeLineageFurcation> furcations;
- FindFurcations(lineage, furcations);
+ FindFurcations(lineage, m_furcations);
- int size = furcations.GetSize();
+ int size = m_furcations.GetSize();
out_furcation_times.Resize(size, 0);
for(int i = 0; i < size; i++){
- out_furcation_times[i] = furcation_time_policy(furcations[i]);
+ out_furcation_times[i] = furcation_time_policy(m_furcations[i]);
if (m_world->GetVerbosity() >= VERBOSE_DETAILS){
cout << "FurcationTime("
- << "parent_id=" << furcations[i].m_parent->GetID() << ") = "
+ << "parent_id=" << m_furcations[i].m_parent->GetID() << ") = "
<< out_furcation_times[i] << endl;
}
}
@@ -251,8 +250,9 @@
}
/*
- For convenience, make a copy of m_g with the indices offset.
- This permits use of the same indices as in the formula for gamma.
+ For convenience, make a copy of inode_dists (internode distances) with the
+ indices offset. This permits use of the same indices as in the formula for
+ gamma.
*/
tArray<int> g = tArray<int>(2, 0) + inode_dists;
@@ -286,18 +286,30 @@
return m_gamma;
}
+
+
+// Accessors.
+const tArray<int> &cAnalyzeTreeStats_Gamma::FurcationTimes(void) const {
+ return m_furcation_times;
+}
+
+const tArray<int> &cAnalyzeTreeStats_Gamma::InternodeDistances(void) const {
+ return m_internode_distances;
+}
+
double cAnalyzeTreeStats_Gamma::Gamma(void){
return m_gamma;
}
+
+
+// Commands.
void cAnalyzeTreeStats_Gamma::AnalyzeBatch(
tList<cAnalyzeGenotype> &genotype_list,
int end_time,
int furcation_time_convention
){
tHashTable<int, int> mapping;
- tArray<int> furcation_times;
- tArray<int> internode_distances;
int (*furcation_time_policy)(cAnalyzeLineageFurcation &furcation);
furcation_time_policy = 0;
@@ -320,16 +332,16 @@
LoadGenotypes(genotype_list);
EstablishLinks(m_gen_array, mapping);
- FindFurcationTimes(m_gen_array, furcation_time_policy, furcation_times);
+ FindFurcationTimes(m_gen_array, furcation_time_policy, m_furcation_times);
- if (end_time < furcation_times[furcation_times.GetSize() - 1]){
+ if (end_time < m_furcation_times[m_furcation_times.GetSize() - 1]){
/* Bad furcation time convention specified. */
if(m_world->GetVerbosity() >= VERBOSE_ON) {
cerr << "Error: Lineage end time comes before last furcation." << endl;
}
}
- FindInternodeDistances(furcation_times, end_time, internode_distances);
- CalculateGamma(internode_distances);
+ FindInternodeDistances(m_furcation_times, end_time, m_internode_distances);
+ CalculateGamma(m_internode_distances);
}
Modified: development/source/analyze/cAnalyzeTreeStats_Gamma.h
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_Gamma.h 2008-05-26 00:57:09 UTC (rev 2594)
+++ development/source/analyze/cAnalyzeTreeStats_Gamma.h 2008-05-26 04:23:14 UTC (rev 2595)
@@ -66,7 +66,9 @@
public:
cWorld* m_world;
tArray<cAnalyzeGenotype *> m_gen_array;
- tArray<int> m_g;
+ tArray<cAnalyzeLineageFurcation> m_furcations;
+ tArray<int> m_furcation_times;
+ tArray<int> m_internode_distances;
double m_gamma;
public:
cAnalyzeTreeStats_Gamma(cWorld* world);
@@ -96,6 +98,10 @@
tArray<int> &out_internode_distances
);
double CalculateGamma(tArray<int> &inode_dists);
+
+ // Accessors.
+ const tArray<int> &FurcationTimes(void) const;
+ const tArray<int> &InternodeDistances(void) const;
double Gamma(void);
// Commands.
More information about the Avida-cvs
mailing list