[Avida-SVN] r2231 - / development development/source/analyze extras/source/testsuites

kaben at myxo.css.msu.edu kaben at myxo.css.msu.edu
Wed Dec 19 12:29:16 PST 2007


Author: kaben
Date: 2007-12-19 15:29:15 -0500 (Wed, 19 Dec 2007)
New Revision: 2231

Added:
   development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.cc
   development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.h
   development/source/analyze/cAnalyzeTreeStats_Gamma.cc
   development/source/analyze/cAnalyzeTreeStats_Gamma.h
Removed:
   development/source/analyze/cAnalyzeGenotypeTreeStats.cc
   development/source/analyze/cAnalyzeGenotypeTreeStats.h
Modified:
   /
   development/CMakeLists.txt
   development/source/analyze/cAnalyze.cc
   development/source/analyze/cAnalyze.h
   extras/source/testsuites/nAnalyze.cpp
   extras/source/testsuites/nAnalyzeGenotype.cpp
   extras/source/testsuites/nMemTracker.cpp
   extras/source/testsuites/nString.cpp
Log:
 r2269 at vallista:  kaben | 2007-12-15 00:13:05 -0800
 Complete implementation of Gabe's gamma statistic.
 



Property changes on: 
___________________________________________________________________
Name: svk:merge
   - 079b078a-dbed-46b9-b3da-37668d4295ca:/avida/xcode-test:1653
c457ea80-0a68-11dc-9323-a45eea2efad5:/private:2228
   + 079b078a-dbed-46b9-b3da-37668d4295ca:/avida/xcode-test:1653
c457ea80-0a68-11dc-9323-a45eea2efad5:/private:2269

Modified: development/CMakeLists.txt
===================================================================
--- development/CMakeLists.txt	2007-12-19 20:28:31 UTC (rev 2230)
+++ development/CMakeLists.txt	2007-12-19 20:29:15 UTC (rev 2231)
@@ -121,7 +121,8 @@
 SET(ANALYZE_SOURCES
   ${ANALYZE_DIR}/cAnalyze.cc
   ${ANALYZE_DIR}/cAnalyzeGenotype.cc
-  ${ANALYZE_DIR}/cAnalyzeGenotypeTreeStats.cc
+  ${ANALYZE_DIR}/cAnalyzeTreeStats_CumulativeStemminess.cc
+  ${ANALYZE_DIR}/cAnalyzeTreeStats_Gamma.cc
   ${ANALYZE_DIR}/cAnalyzeJobQueue.cc
   ${ANALYZE_DIR}/cAnalyzeJobWorker.cc
   ${ANALYZE_DIR}/cMutationalNeighborhood.cc
@@ -566,7 +567,66 @@
   INSTALL_TARGETS(/work unit-tests)
 ENDIF(AVD_UNIT_TESTS)
 
+OPTION(AVD_TESTSUITES
+  "Enable the testsuites executable.  Running this target will test various low level functionality."
+  OFF
+)
+IF(AVD_TESTSUITES)
+  SET(TESTSUITES_DIR
+    "${PROJECT_SOURCE_DIR}/../extras/source"
+    CACHE PATH
+    "The path to testsuite files."
+  )
+  IF(EXISTS ${TESTSUITES_DIR})
+    INCLUDE_DIRECTORIES(${TESTSUITES_DIR}/tools)
+    SET(TESTSUITES_TOOL_SOURCES
+      ${TESTSUITES_DIR}/tools/cAnalyzeTestFixture.cpp
+      ${TESTSUITES_DIR}/tools/cConsoleCatcher.cpp
+      ${TESTSUITES_DIR}/tools/cMemTracker.cpp
+      ${TESTSUITES_DIR}/tools/cTestDriver.cpp
+      ${TESTSUITES_DIR}/tools/cTestLib.cpp
+      ${TESTSUITES_DIR}/tools/cTestSettings.cpp
+    )
+    SET(TESTSUITES_SOURCES
+      ${TESTSUITES_DIR}/testsuites/nAGenotype.cpp
+      ${TESTSUITES_DIR}/testsuites/nAnalyze.cpp
+      ${TESTSUITES_DIR}/testsuites/nAnalyzeGenotype.cpp
+      ${TESTSUITES_DIR}/testsuites/nAnalyzeTestFixture.cpp
+      ${TESTSUITES_DIR}/testsuites/nChangeList.cpp
+      ${TESTSUITES_DIR}/testsuites/nConsoleCatcher.cpp
+      ${TESTSUITES_DIR}/testsuites/nDataEntry.cpp
+      ${TESTSUITES_DIR}/testsuites/nDataFile.cpp
+      ${TESTSUITES_DIR}/testsuites/nFile.cpp
+      ${TESTSUITES_DIR}/testsuites/nFixedCoords.cpp
+      ${TESTSUITES_DIR}/testsuites/nInitFile.cpp
+      ${TESTSUITES_DIR}/testsuites/nMemTracker.cpp
+      ${TESTSUITES_DIR}/testsuites/nRandom.cpp
+      ${TESTSUITES_DIR}/testsuites/nString.cpp
+      ${TESTSUITES_DIR}/testsuites/nStringList.cpp
+      ${TESTSUITES_DIR}/testsuites/nTemplate.cpp
+      ${TESTSUITES_DIR}/testsuites/nTestDriver.cpp
+      ${TESTSUITES_DIR}/testsuites/nTestLib.cpp
+      ${TESTSUITES_DIR}/testsuites/nTestSettings.cpp
+      ${TESTSUITES_DIR}/testsuites/nTreeStats.cpp
+    )
+    SET(TESTSUITES_MAIN_SOURCES
+      ${TESTSUITES_DIR}/targets/TestAvida/TestAvida.cpp
+    )
+    ADD_EXECUTABLE(testsuites ${TESTSUITES_TOOL_SOURCES} ${TESTSUITES_SOURCES} ${TESTSUITES_MAIN_SOURCES})
+    INSTALL_TARGETS(/work testsuites)
 
+    SET(TESTSUITES_LIBS avidacore)  
+    IF(NOT MSVC)
+      LIST(APPEND TESTSUITES_LIBS pthread)
+    ENDIF(NOT MSVC)
+    IF(AVD_ENABLE_TCMALLOC)
+      LIST(APPEND TESTSUITES_LIBS tcmalloc)
+    ENDIF(AVD_ENABLE_TCMALLOC)  
+    TARGET_LINK_LIBRARIES(testsuites ${TESTSUITES_LIBS})
+
+  ENDIF(EXISTS ${TESTSUITES_DIR})
+ENDIF(AVD_TESTSUITES)
+
 # Default Configuration Files
 # - Installed into the work directory alongside selected targets
 # ------------------------------------------------------------------------------

Modified: development/source/analyze/cAnalyze.cc
===================================================================
--- development/source/analyze/cAnalyze.cc	2007-12-19 20:28:31 UTC (rev 2230)
+++ development/source/analyze/cAnalyze.cc	2007-12-19 20:29:15 UTC (rev 2231)
@@ -40,7 +40,8 @@
 #include "cAnalyzeFlowCommandDef.h"
 #include "cAnalyzeFunction.h"
 #include "cAnalyzeGenotype.h"
-#include "cAnalyzeGenotypeTreeStats.h"
+#include "cAnalyzeTreeStats_CumulativeStemminess.h"
+#include "cAnalyzeTreeStats_Gamma.h"
 #include "tAnalyzeJob.h"
 #include "cAvidaContext.h"
 #include "cDataFile.h"
@@ -3242,109 +3243,64 @@
   fp << "# 1: Average cumulative stemminess" << endl;
   fp << endl;
   
-  cAnalyzeGenotypeTreeStats agts(m_world);
+  cAnalyzeTreeStats_CumulativeStemminess agts(m_world);
   agts.AnalyzeBatchTree(batch[cur_batch].List());
 
   fp << agts.AverageStemminess();
   fp << endl;
+}
 
-  /*
-  Below is the original implementation by Ofria.
-  -- kgn
-  */
 
-  //cAnalyzeGenotype * genotype = NULL;
-  //tListIterator<cAnalyzeGenotype> batch_it(batch[cur_batch].List());
-  //const int num_gens = batch[cur_batch].List().GetSize();
-  //
-  //// Put all of the genotypes in an array for easy reference and collect
-  //// other information on them as we process them.
-  //tArray<cAnalyzeGenotype *> gen_array(num_gens);
-  //tHashTable<int, int> id_hash;  // Store array pos for each id.
-  //tArray<int> id_array(num_gens), pid_array(num_gens);
-  //tArray<int> depth_array(num_gens), birth_array(num_gens);
-  //int array_pos = 0;
-  //while ((genotype = batch_it.Next()) != NULL) {
-  //  // Put the genotype in an array.
-  //  gen_array[array_pos] = genotype;
-  //  id_hash.Add(genotype->GetID(), array_pos);
-  //  id_array[array_pos] = genotype->GetID();
-  //  pid_array[array_pos] = genotype->GetParentID();
-  //  depth_array[array_pos] = genotype->GetDepth();
-  //  birth_array[array_pos] = genotype->GetUpdateBorn();
-  //  array_pos++;
-  //}
-  //
-  //// Now collect information about the offspring of each individual.
-  //tArray<int> ppos_array(num_gens), offspring_count(num_gens);
-  //offspring_count.SetAll(0);
-  //for (int pos = 0; pos < num_gens; pos++) {
-  //  int parent_id = gen_array[pos]->GetParentID();
-  //  if (parent_id == -1) {  // Organism has no parent (i.e., ancestor)
-  //    ppos_array[pos] = -1;
-  //    continue;
-  //  }
-  //  int parent_pos = -1;
-  //  id_hash.Find(parent_id, parent_pos);
-  //  ppos_array[pos] = parent_pos;
-  //  offspring_count[parent_pos]++;
-  //}
-  //
-  //// For each genotype, figure out how far back you need to go to get to a
-  //// branch point.
-  //tArray<int> branch_dist_array(num_gens);
-  //tArray<int> branch_pos_array(num_gens);
-  //branch_dist_array.SetAll(-1);
-  //branch_pos_array.SetAll(-1);
-  //bool found = true;
-  //int loop_count = 0;
-  //while (found == true) {
-  //  found = false;
-  //  for (int pos = 0; pos < num_gens; pos++) {
-  //    if (branch_dist_array[pos] > -1) continue; // continue if its set.
-  //    found = true;
-  //    int parent_pos = ppos_array[pos];
-  //    if (parent_pos == -1) branch_dist_array[pos] = 0;  // Org is root.
-  //    else if (offspring_count[parent_pos] > 1) {        // Parent is branch.
-  //      branch_dist_array[pos] = 1;
-  //      branch_pos_array[pos] = parent_pos;
-  //    }
-  //    else if (branch_dist_array[parent_pos] > -1) {     // Parent calculated.
-  //      branch_dist_array[pos] = branch_dist_array[parent_pos]+1;
-  //      branch_pos_array[pos] = branch_pos_array[parent_pos];
-  //    }
-  //    // Otherwise, we are not yet ready to calculate this entry.
-  //  }
-  //  loop_count++;
-  //}
-  //
-  //
-  //// Cumulative Stemminess
-  //for (int pos = 0; pos < num_gens; pos++) {
-  //  // We're only interested in internal n-furcating nodes.
-  //  if (pid_array[pos] == -1) continue;  // Don't want root.
-  //  if (offspring_count[pos] <= 1) continue; // No leaves or nonfurcating nodes
-  //  
-  //  // @CAO Find distance to all children.
-  //  // @CAO Find distance to parent branch.
-  //  // @CAO DO math.
-  //}
-  //
-  //
-  //cout << "LOOP COUNT:" << loop_count << endl;
-  //for (int i = 0; i < num_gens; i++) {
-  //  int branch_pos = branch_pos_array[i];
-  //  int branch_id = (branch_pos == -1) ? -1 : id_array[branch_pos];
-  //  cout << i << " "
-  //    << id_array[i] << " "
-  //    << offspring_count[i] << " "
-  //    << branch_dist_array[i] << " "
-  //    << branch_id << " "
-  //    << endl;
-  //}
+// Calculate cumulative stemmines for trees in population.
+void cAnalyze::CommandPrintCumulativeStemminess(cString cur_string)
+{
+  if (m_world->GetVerbosity() >= VERBOSE_ON) cout << "Printing cumulative stemmines for batch "
+    << cur_batch << endl;
+  else cout << "Printing cumulative stemmines..." << endl;
+  
+  // Load in the variables...
+  cString filename("cumulative_stemminess.dat");
+  if (cur_string.GetSize() != 0) filename = cur_string.PopWord();
+  
+  ofstream& fp = m_world->GetDataFileOFStream(filename);
+  
+  fp << "# Legend:" << endl;
+  fp << "# 1: Average cumulative stemminess" << endl;
+  fp << endl;
+  
+  cAnalyzeTreeStats_CumulativeStemminess agts(m_world);
+  agts.AnalyzeBatchTree(batch[cur_batch].List());
+  
+  fp << agts.AverageStemminess();
+  fp << endl;
 }
 
 
+// Calculate Pybus-Harvey gamma statistic for trees in population.
+void cAnalyze::CommandPrintGamma(cString cur_string)
+{
+  if (m_world->GetVerbosity() >= VERBOSE_ON) cout << "Printing Pybus-Harvey gamma statistic for batch "
+    << cur_batch << endl;
+  else cout << "Printing Pybus-Harvey gamma statistic..." << endl;
+  
+  // Load in the variables...
+  cString filename("gamma.dat");
+  if (cur_string.GetSize() != 0) filename = cur_string.PopWord();
+  
+  ofstream& fp = m_world->GetDataFileOFStream(filename);
+  
+  fp << "# Legend:" << endl;
+  fp << "# 1: Pybus-Harvey gamma statistic" << endl;
+  fp << endl;
+  
+  cAnalyzeTreeStats_Gamma agts(m_world);
+  agts.AnalyzeBatch(batch[cur_batch].List());
+  
+  fp << agts.Gamma();
+  fp << endl;
+}
+
+
 void cAnalyze::AnalyzeCommunityComplexity(cString cur_string)
 {
   /////////////////////////////////////////////////////////////////////
@@ -8618,6 +8574,8 @@
   AddLibraryDef("PRINT_PHENOTYPES", &cAnalyze::CommandPrintPhenotypes);
   AddLibraryDef("PRINT_DIVERSITY", &cAnalyze::CommandPrintDiversity);
   AddLibraryDef("PRINT_TREE_STATS", &cAnalyze::CommandPrintTreeStats);
+  AddLibraryDef("PRINT_CUMULATIVE_STEMMINESS", &cAnalyze::CommandPrintCumulativeStemminess);
+  AddLibraryDef("PRINT_GAMMA", &cAnalyze::CommandPrintGamma);
   AddLibraryDef("COMMUNITY_COMPLEXITY", &cAnalyze::AnalyzeCommunityComplexity);
   AddLibraryDef("PRINT_RESOURCE_FITNESS_MAP", &cAnalyze::CommandPrintResourceFitnessMap);
   

Modified: development/source/analyze/cAnalyze.h
===================================================================
--- development/source/analyze/cAnalyze.h	2007-12-19 20:28:31 UTC (rev 2230)
+++ development/source/analyze/cAnalyze.h	2007-12-19 20:29:15 UTC (rev 2231)
@@ -258,6 +258,8 @@
   void CommandPrintPhenotypes(cString cur_string);
   void CommandPrintDiversity(cString cur_string);
   void CommandPrintTreeStats(cString cur_string);
+  void CommandPrintCumulativeStemminess(cString cur_string);
+  void CommandPrintGamma(cString cur_string);
   void PhyloCommunityComplexity(cString cur_string);
   void AnalyzeCommunityComplexity(cString cur_string);
   void CommandPrintResourceFitnessMap(cString cur_string);

Deleted: development/source/analyze/cAnalyzeGenotypeTreeStats.cc
===================================================================
--- development/source/analyze/cAnalyzeGenotypeTreeStats.cc	2007-12-19 20:28:31 UTC (rev 2230)
+++ development/source/analyze/cAnalyzeGenotypeTreeStats.cc	2007-12-19 20:29:15 UTC (rev 2231)
@@ -1,393 +0,0 @@
-
-#include "cAnalyzeGenotypeTreeStats.h"
-
-#include "cAnalyzeGenotype.h"
-#include "tHashTable.h"
-#include "cWorld.h"
-
-
-void cAnalyzeGenotypeTreeStats::PrintAGLData(tArray<cAGLData> &agl){
-  for(int i=0; i < agl.GetSize(); i++){
-    cout << i << ":";
-    cout << " " << agl[i].id;
-    cout << " " << agl[i].pid;
-    cout << " " << agl[i].depth;
-    cout << " " << agl[i].birth;
-    cout << " " << agl[i].ppos;
-    cout << " " << agl[i].offspring_count;
-    cout << " " << agl[i].anc_branch_dist;
-    cout << " " << agl[i].anc_branch_id;
-    cout << " " << agl[i].anc_branch_pos;
-    cout << " " << agl[i].off_branch_dist_acc;
-    cout << " " << agl[i].cumulative_stemminess;
-    cout << " " << agl[i].traversal_visited;
-    for(int j=0; j < agl[i].offspring_positions.GetSize(); j++){
-      cout << " " << agl[i].offspring_positions[j];
-    }
-    cout << endl;
-  }
-}
-
-void cAnalyzeGenotypeTreeStats::AnalyzeBatchTree(tList<cAnalyzeGenotype> &genotype_list){
-  cAnalyzeGenotype * genotype = NULL;
-  tListIterator<cAnalyzeGenotype> batch_it(genotype_list);
-  const int num_gens = genotype_list.GetSize();
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Number of genotypes: " << num_gens << endl;
-  }
-
-
-  int array_pos = 0;
-
-  /*
-  Put all of the genotypes in an array for easy reference and collect other
-  information on them as we process them. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Scanning genotypes..." << endl;
-  }
-  tArray<cAnalyzeGenotype *> gen_array(num_gens);
-  tHashTable<int, int> id_hash;  // Store array pos for each id.
-  tArray<int> id_array(num_gens), pid_array(num_gens);
-  tArray<int> depth_array(num_gens), birth_array(num_gens);
-
-  array_pos = 0;
-  batch_it.Reset();
-  while ((genotype = batch_it.Next()) != NULL) {
-    id_hash.Add(genotype->GetID(), array_pos);
-    array_pos++;
-  }
-
-  m_agl.Resize(num_gens);
-  array_pos = 0;
-  batch_it.Reset();
-  while ((genotype = batch_it.Next()) != NULL) {
-    // Put the genotype in an array.
-    m_agl[array_pos].genotype = genotype;
-    m_agl[array_pos].id = genotype->GetID();
-    m_agl[array_pos].pid = genotype->GetParentID();
-    m_agl[array_pos].depth = genotype->GetDepth();
-    m_agl[array_pos].birth = genotype->GetUpdateBorn();
-    array_pos++;
-  }
-
-  //// Now collect information about the offspring of each individual. {{{4
-  tArray<int> ppos_array(num_gens), offspring_count(num_gens);
-  offspring_count.SetAll(0);
-
-  // For each genotype, figure out how far back you need to go to get to a branch point. {{{4
-  tArray<int> anc_branch_dist_array(num_gens);
-  tArray<int> anc_branch_pos_array(num_gens);
-  anc_branch_dist_array.SetAll(-1);
-  anc_branch_pos_array.SetAll(-1);
-  bool found = true;
-  int loop_count = 0;
-
-  /*
-  Link each offspring to its parent. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Assembling tree..." << endl;
-  }
-  for (int pos = 0; pos < num_gens; pos++) {
-    //cAnalyzeGenotype * genotype = gen_array[pos];
-    cAnalyzeGenotype * genotype = m_agl[pos].genotype;
-    int parent_id = genotype->GetParentID();
-    if (-1 != parent_id){
-      bool found_parent = id_hash.Find(parent_id, m_agl[pos].ppos);
-      if (found_parent){
-        int parent_position = m_agl[pos].ppos;
-        m_agl[parent_position].offspring_positions.Push(pos);
-        ///* XXX I think I'll be able to remove this. */
-        //cAnalyzeGenotype * parent_genotype = m_agl[parent_position].genotype;
-        //genotype->LinkParent(parent_genotype);
-      } else {
-        if (m_world->GetVerbosity() >= VERBOSE_ON) {
-          cerr << "Error: the parent of a non-root tree node is missing - " << endl;
-        }
-        return;
-      }
-    }
-  }
-
-  /*
-  Count offspring of each parent. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Loading offspring counts..." << endl;
-  }
-  for (int pos = 0; pos < num_gens; pos++) {
-    m_agl[pos].offspring_count = m_agl[pos].offspring_positions.GetSize();
-  }
-
-
-  /*
-  For each genotype, figure out how far back you need to go to get to a branch point. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Finding branch points..." << endl;
-  }
-  found = true;
-  loop_count = 0;
-  while (found == true) {
-    if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
-      cout << endl << "  wave " << loop_count << "...";
-    }
-    found = false;
-
-    int genotypes_already_set = 0;
-    int genotypes_identified_as_root = 0;
-    int genotypes_with_branching_parent = 0;
-    int genotypes_with_calced_parent = 0;
-    int genotypes_not_ready = 0;
-
-    for (int pos = 0; pos < num_gens; pos++) {
-      if (m_agl[pos].anc_branch_dist > -1){
-        genotypes_already_set++;
-        continue; // continue if its set.
-      }
-      found = true;
-      int parent_pos = m_agl[pos].ppos;
-      if (parent_pos == -1) {
-        genotypes_identified_as_root++;
-        m_agl[pos].anc_branch_dist = 0;  // Org is root.
-      } else if (m_agl[parent_pos].offspring_count > 1) {        // Parent is branch.
-        genotypes_with_branching_parent++;
-        m_agl[pos].anc_branch_dist = 1;
-        m_agl[pos].anc_branch_id = m_agl[parent_pos].id;
-        m_agl[pos].anc_branch_pos = parent_pos;
-      } else if (m_agl[parent_pos].anc_branch_dist > -1) {     // Parent calculated.
-        genotypes_with_calced_parent++;
-        m_agl[pos].anc_branch_dist = m_agl[parent_pos].anc_branch_dist + 1;
-        m_agl[pos].anc_branch_id = m_agl[parent_pos].anc_branch_id;
-        m_agl[pos].anc_branch_pos = m_agl[parent_pos].anc_branch_pos;
-      } else {
-        genotypes_not_ready++;
-        // Otherwise, we are not yet ready to calculate this entry.
-      }
-    }
-    if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
-      cout << " (" << genotypes_already_set << "," << genotypes_identified_as_root << "," << genotypes_with_branching_parent << "," << genotypes_with_calced_parent << "," << genotypes_not_ready << ")" << endl;
-    }
-    loop_count++;
-  }
-
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Building n-furcating subtree..." << endl;
-  }
-  // compute number of subtree nodes.
-  int branch_tree_size = 0;
-  for (int pos = 0; pos < num_gens; pos++) {
-    //if(m_agl[pos].offspring_count != 1){
-    if(m_agl[pos].offspring_count > 1){
-      branch_tree_size++;
-    }
-  }
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Number of n-furcating nodes: " << branch_tree_size << endl;
-  }
-  if (branch_tree_size <= 0){
-    if (m_world->GetVerbosity() >= VERBOSE_ON) {
-      cerr << "Error: no branches found in tree - " << endl;
-    }
-    return;
-  }
-
-  m_agl2.Resize(branch_tree_size);  // Store agl data for each id.
-  tHashTable<int, int> id_hash_2;
-  int array_pos_2 = 0;
-  if (true) for (int pos = 0; pos < num_gens; pos++) {
-    int offs_count = m_agl[pos].offspring_count;
-    //if (offs_count != 1){
-    if (offs_count > 1){
-      m_agl2[array_pos_2].id = m_agl[pos].id;
-      m_agl2[array_pos_2].pid = m_agl[pos].pid;
-      m_agl2[array_pos_2].depth = m_agl[pos].depth;
-      m_agl2[array_pos_2].birth = m_agl[pos].birth;
-      m_agl2[array_pos_2].anc_branch_dist = m_agl[pos].anc_branch_dist;
-      m_agl2[array_pos_2].anc_branch_id = m_agl[pos].anc_branch_id;
-      /*
-      missing still are ppos (skip this), offspring_count (redundant),
-      anc_branch_pos, off_branch_dist_acc (to be calculated),
-      offspring_positions
-      */
-      id_hash_2.Add(m_agl2[array_pos_2].id, array_pos_2);
-      array_pos_2++;
-    }
-  }
-
-  // find branch ancestor positions. {{{4
-  for (int pos = 0; pos < branch_tree_size; pos++){
-    int anc_branch_id = m_agl2[pos].anc_branch_id;
-    id_hash_2.Find(anc_branch_id, m_agl2[pos].anc_branch_pos);
-    int anc_branch_pos = m_agl2[pos].anc_branch_pos;
-    if(0 <= anc_branch_pos){
-      m_agl2[anc_branch_pos].offspring_positions.Push(pos);
-    }
-  }
-  
-  /*
-  For DFS of branch tree, locate root. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Finding root of n-furcating subtree..." << endl;
-  }
-  cAGLData *root = 0;
-  for (int pos = 0; pos < branch_tree_size; pos++){
-    m_agl2[pos].traversal_visited = false;
-    /*
-    root : anc_branch_dist: 0 anc_branch_id: -1 anc_branch_pos: -1 
-
-    Only one of these conditions should be needed. I'm mixing-in a
-    sanity check here. I ought to move it...
-    */
-    if( (m_agl2[pos].anc_branch_dist == 0)
-      ||(m_agl2[pos].anc_branch_id == -1)
-      ||(m_agl2[pos].anc_branch_pos == -1)
-    ){
-      root = &(m_agl2[pos]);
-      /* Sanity check. */
-      if(
-        !(
-          //(m_agl2[pos].anc_branch_dist == 0) &&
-          (m_agl2[pos].anc_branch_id == -1) &&
-          (m_agl2[pos].anc_branch_pos == -1)
-        )
-      ){
-        if (m_world->GetVerbosity() >= VERBOSE_ON) {
-          cerr << "Error: while looking for root of subtree, found inconsistencies - " << endl;
-          cerr << " root ancestor-branch-distance: " << m_agl2[pos].anc_branch_dist << endl;
-          cerr << " root ancestor-branch-id: " << m_agl2[pos].anc_branch_id << endl;
-          cerr << " root ancestor-branch-pos: " << m_agl2[pos].anc_branch_id << endl;
-        }
-        return;
-      }
-    }
-  }
-
-  /*
-  DFS of branch tree, to accumulate branch distances. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Accumulating branch distances..." << endl;
-  }
-  tList<cAGLData> dfs_stack;
-  if(0 != root){
-    /* DFS. */
-    dfs_stack.Push(root);
-    cAGLData *node = 0;
-    while (0 < dfs_stack.GetSize()){
-      node = dfs_stack.Pop();
-      if (0 != node){
-        if (! node->traversal_visited){
-          dfs_stack.Push(node);
-          node->off_branch_dist_acc = 0;
-          for (int i = 0; i < node->offspring_positions.GetSize(); i++){
-            int pos = node->offspring_positions[i];
-            if (! m_agl2[pos].traversal_visited){
-              dfs_stack.Push(&(m_agl2[pos]));
-            }
-          }
-          node->traversal_visited = true;
-        } else {
-          /*
-          Child nodes, if any, have been visited and have added their
-          off_branch_dist_acc to this node.
-          */
-          if(0 <= node->anc_branch_pos){
-            /*
-            Only accumulate to parent if there is a parent (i.e., this
-            is not the root.)
-            */
-            m_agl2[node->anc_branch_pos].off_branch_dist_acc += node->anc_branch_dist;
-            m_agl2[node->anc_branch_pos].off_branch_dist_acc += node->off_branch_dist_acc;
-          }
-        }
-      }
-    }
-  } else {
-    if (m_world->GetVerbosity() >= VERBOSE_ON) {
-      cerr << "Error: couldn't find root of subtree - " << endl;
-    }
-    return;
-  }
-
-  /*
-  Compute cumulative stemminesses. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Computing cumulative stemminesses..." << endl;
-  }
-  for (int pos = 0; pos < branch_tree_size; pos++){
-    if (0 == m_agl2[pos].anc_branch_dist){
-      /* Correct stemminess for root... */
-      /*
-      Let me rephrase that. If anc_branch_dist is zero and all is well,
-      we're dealing with the root.
-      */
-      m_agl2[pos].cumulative_stemminess = 0.0;
-    } else {
-      m_agl2[pos].cumulative_stemminess =
-      (
-        (double)(m_agl2[pos].anc_branch_dist)
-        /
-        (
-          (double)(m_agl2[pos].off_branch_dist_acc) + (double)(m_agl2[pos].anc_branch_dist)
-        )
-      );
-    }
-  }
-
-  /*
-  Compute average cumulative stemminess. {{{4
-  */
-  if (m_world->GetVerbosity() >= VERBOSE_ON) {
-    cout << "Computing average cumulative stemminess..." << endl;
-  }
-  m_stemminess_sum = 0.0;
-  m_average_stemminess = 0.0;
-  m_inner_nodes = 0;
-  if (1 < branch_tree_size) {
-    for (int pos = 0; pos < branch_tree_size; pos++){
-      bool not_leaf = true;
-      if (m_should_exclude_leaves) {
-        not_leaf = (0 < m_agl2[pos].off_branch_dist_acc);
-      }
-      bool not_root = (0 < m_agl2[pos].anc_branch_id);
-      if (not_leaf && not_root){
-        m_stemminess_sum += m_agl2[pos].cumulative_stemminess;
-        m_inner_nodes++;
-      }
-    }
-  }
-  if(0 < m_inner_nodes){
-    m_average_stemminess = m_stemminess_sum / (m_inner_nodes);
-  }
-
-  /*
-  Print branch tree. {{{4
-  Use to get expected data for making following test.
-  */
-
-  if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
-    for (int pos = 0; pos < branch_tree_size; pos++){
-      cout << "id: " << m_agl2[pos].id;
-      cout << " offspring_count: " << m_agl2[pos].offspring_positions.GetSize();
-      cout << " anc_branch_id: " << m_agl2[pos].anc_branch_id;
-      cout << " anc_branch_dist: " << m_agl2[pos].anc_branch_dist;
-      cout << " off_branch_dist_acc: " << m_agl2[pos].off_branch_dist_acc;
-      cout << " cumulative_stemminess: " << m_agl2[pos].cumulative_stemminess;
-      cout << endl;
-      if (0 < m_agl2[pos].offspring_positions.GetSize()) {
-        cout << "  offspring ids:";
-        for (int i = 0; i < m_agl2[pos].offspring_positions.GetSize(); i++){
-          cout << " " << m_agl2[pos].offspring_positions[i];
-        }
-        cout << endl;
-      }
-    }
-    cout << "m_stemminess_sum: " << m_stemminess_sum << endl;
-    cout << "m_average_stemminess: " << m_average_stemminess << endl;
-    cout << "m_inner_nodes: " << m_inner_nodes << endl;
-  }
-}

Deleted: development/source/analyze/cAnalyzeGenotypeTreeStats.h
===================================================================
--- development/source/analyze/cAnalyzeGenotypeTreeStats.h	2007-12-19 20:28:31 UTC (rev 2230)
+++ development/source/analyze/cAnalyzeGenotypeTreeStats.h	2007-12-19 20:29:15 UTC (rev 2231)
@@ -1,98 +0,0 @@
-#ifndef cAnalyzeGenotypeTreeStats_h
-#define cAnalyzeGenotypeTreeStats_h
-
-#ifndef tArray_h
-#include "tArray.h"
-#endif
-#ifndef tList_h
-#include "tList.h"
-#endif
-
-
-class cAnalyzeGenotype;
-class cWorld;
-
-
-/*
-Convenience class cAGLData serves to collect various data I found I needed
-while building the class cAnalyzeGenotypeTreeStats.
-*/
-class cAGLData {
-public:
-  cAnalyzeGenotype *genotype;
-  int id;
-  int pid;
-  int depth;
-  int birth;
-  int ppos;
-  int offspring_count;
-  int anc_branch_dist;
-  int anc_branch_id;
-  int anc_branch_pos;
-  int off_branch_dist_acc;
-  double cumulative_stemminess;
-  bool traversal_visited;
-  tArray<int> offspring_positions;
-
-  cAGLData()
-  : genotype(0)
-  , id(-1)
-  , pid(-1)
-  , depth(-1)
-  , birth(-1)
-  , ppos(-1)
-  , offspring_count(-1)
-  , anc_branch_dist(-1)
-  , anc_branch_id(-1)
-  , anc_branch_pos(-1)
-  , off_branch_dist_acc(-1)
-  , cumulative_stemminess(-1.)
-  , traversal_visited(false)
-  , offspring_positions(0)
-  {}
-};
-
-/*
-class cAnalyzeGenotypeTreeStats
-
-This class will eventually collect various kinds of statistics about various
-kinds of trees (e.g., phylogenetic).
-
-For now it only figures out "average cumulative stemminess".
-*/
-class cAnalyzeGenotypeTreeStats {
-public:
-  tArray<cAGLData> m_agl;
-  tArray<cAGLData> m_agl2;
-
-  double m_stemminess_sum;
-  double m_average_stemminess;
-  int m_inner_nodes;
-  bool m_should_exclude_leaves;
-
-  cWorld* m_world;
-
-public:
-  cAnalyzeGenotypeTreeStats(cWorld* world)
-  : m_agl(0)
-  , m_agl2(0)
-  , m_stemminess_sum(0.0)
-  , m_average_stemminess(0.0)
-  , m_inner_nodes(0)
-  , m_should_exclude_leaves(true)
-  , m_world(world)
-  {}
-
-  tArray<cAGLData> &AGL(){ return m_agl; }
-  tArray<cAGLData> &AGL2(){ return m_agl2; }
-  double StemminessSum(){ return m_stemminess_sum; }
-  double AverageStemminess(){ return m_average_stemminess; }
-  int InnerNodes(){ return m_inner_nodes; }
-  bool ExcludesLeaves(){ return m_should_exclude_leaves; }
-
-  void PrintAGLData(tArray<cAGLData> &agl);
-
-  void AnalyzeBatchTree(tList<cAnalyzeGenotype> &genotype_list);
-};
-
-#endif

Copied: development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.cc (from rev 2226, development/source/analyze/cAnalyzeGenotypeTreeStats.cc)
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.cc	                        (rev 0)
+++ development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.cc	2007-12-19 20:29:15 UTC (rev 2231)
@@ -0,0 +1,443 @@
+/*
+ *  cAnalyzeTreeStats_CumulativeStemminess.cc
+ *  Avida at vallista
+ *
+ *  Created by Kaben Nanlohy on 2007.12.03.
+ *  Copyright 1999-2007 Michigan State University. All rights reserved.
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; version 2
+ *  of the License.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ */
+
+#include "cAnalyzeTreeStats_CumulativeStemminess.h"
+
+#include "cAnalyzeGenotype.h"
+#include "tHashTable.h"
+#include "cWorld.h"
+
+
+cAGLData::cAGLData()
+: genotype(0)
+, id(-1)
+, pid(-1)
+, depth(-1)
+, birth(-1)
+, ppos(-1)
+, offspring_count(-1)
+, anc_branch_dist(-1)
+, anc_branch_id(-1)
+, anc_branch_pos(-1)
+, off_branch_dist_acc(-1)
+, cumulative_stemminess(-1.)
+, traversal_visited(false)
+, offspring_positions(0)
+{}
+
+cAnalyzeTreeStats_CumulativeStemminess::cAnalyzeTreeStats_CumulativeStemminess(cWorld* world)
+: m_agl(0)
+, m_agl2(0)
+, m_stemminess_sum(0.0)
+, m_average_stemminess(0.0)
+, m_inner_nodes(0)
+, m_should_exclude_leaves(true)
+, m_world(world)
+{}
+
+
+void cAnalyzeTreeStats_CumulativeStemminess::PrintAGLData(tArray<cAGLData> &agl){
+  for(int i=0; i < agl.GetSize(); i++){
+    cout << i << ":";
+    cout << " " << agl[i].id;
+    cout << " " << agl[i].pid;
+    cout << " " << agl[i].depth;
+    cout << " " << agl[i].birth;
+    cout << " " << agl[i].ppos;
+    cout << " " << agl[i].offspring_count;
+    cout << " " << agl[i].anc_branch_dist;
+    cout << " " << agl[i].anc_branch_id;
+    cout << " " << agl[i].anc_branch_pos;
+    cout << " " << agl[i].off_branch_dist_acc;
+    cout << " " << agl[i].cumulative_stemminess;
+    cout << " " << agl[i].traversal_visited;
+    for(int j=0; j < agl[i].offspring_positions.GetSize(); j++){
+      cout << " " << agl[i].offspring_positions[j];
+    }
+    cout << endl;
+  }
+}
+
+void cAnalyzeTreeStats_CumulativeStemminess::AnalyzeBatchTree(tList<cAnalyzeGenotype> &genotype_list){
+  cAnalyzeGenotype * genotype = NULL;
+  tListIterator<cAnalyzeGenotype> batch_it(genotype_list);
+  const int num_gens = genotype_list.GetSize();
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Number of genotypes: " << num_gens << endl;
+  }
+
+
+  int array_pos = 0;
+
+  /*
+  Put all of the genotypes in an array for easy reference and collect other
+  information on them as we process them. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Scanning genotypes..." << endl;
+  }
+  tArray<cAnalyzeGenotype *> gen_array(num_gens);
+  tHashTable<int, int> id_hash;  // Store array pos for each id.
+  tArray<int> id_array(num_gens), pid_array(num_gens);
+  tArray<int> depth_array(num_gens), birth_array(num_gens);
+
+  array_pos = 0;
+  batch_it.Reset();
+  while ((genotype = batch_it.Next()) != NULL) {
+    id_hash.Add(genotype->GetID(), array_pos);
+    array_pos++;
+  }
+
+  m_agl.Resize(num_gens);
+  array_pos = 0;
+  batch_it.Reset();
+  while ((genotype = batch_it.Next()) != NULL) {
+    // Put the genotype in an array.
+    m_agl[array_pos].genotype = genotype;
+    m_agl[array_pos].id = genotype->GetID();
+    m_agl[array_pos].pid = genotype->GetParentID();
+    m_agl[array_pos].depth = genotype->GetDepth();
+    m_agl[array_pos].birth = genotype->GetUpdateBorn();
+    array_pos++;
+  }
+
+  //// Now collect information about the offspring of each individual. {{{4
+  tArray<int> ppos_array(num_gens), offspring_count(num_gens);
+  offspring_count.SetAll(0);
+
+  // For each genotype, figure out how far back you need to go to get to a branch point. {{{4
+  tArray<int> anc_branch_dist_array(num_gens);
+  tArray<int> anc_branch_pos_array(num_gens);
+  anc_branch_dist_array.SetAll(-1);
+  anc_branch_pos_array.SetAll(-1);
+  bool found = true;
+  int loop_count = 0;
+
+  /*
+  Link each offspring to its parent. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Assembling tree..." << endl;
+  }
+  for (int pos = 0; pos < num_gens; pos++) {
+    //cAnalyzeGenotype * genotype = gen_array[pos];
+    cAnalyzeGenotype * genotype = m_agl[pos].genotype;
+    int parent_id = genotype->GetParentID();
+    if (-1 != parent_id){
+      bool found_parent = id_hash.Find(parent_id, m_agl[pos].ppos);
+      if (found_parent){
+        int parent_position = m_agl[pos].ppos;
+        m_agl[parent_position].offspring_positions.Push(pos);
+        ///* XXX I think I'll be able to remove this. */
+        //cAnalyzeGenotype * parent_genotype = m_agl[parent_position].genotype;
+        //genotype->LinkParent(parent_genotype);
+      } else {
+        if (m_world->GetVerbosity() >= VERBOSE_ON) {
+          cerr << "Error: the parent of a non-root tree node is missing - " << endl;
+        }
+        return;
+      }
+    }
+  }
+
+  /*
+  Count offspring of each parent. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Loading offspring counts..." << endl;
+  }
+  for (int pos = 0; pos < num_gens; pos++) {
+    m_agl[pos].offspring_count = m_agl[pos].offspring_positions.GetSize();
+  }
+
+
+  /*
+  For each genotype, figure out how far back you need to go to get to a branch point. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Finding branch points..." << endl;
+  }
+  found = true;
+  loop_count = 0;
+  while (found == true) {
+    if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
+      cout << endl << "  wave " << loop_count << "...";
+    }
+    found = false;
+
+    int genotypes_already_set = 0;
+    int genotypes_identified_as_root = 0;
+    int genotypes_with_branching_parent = 0;
+    int genotypes_with_calced_parent = 0;
+    int genotypes_not_ready = 0;
+
+    for (int pos = 0; pos < num_gens; pos++) {
+      if (m_agl[pos].anc_branch_dist > -1){
+        genotypes_already_set++;
+        continue; // continue if its set.
+      }
+      found = true;
+      int parent_pos = m_agl[pos].ppos;
+      if (parent_pos == -1) {
+        genotypes_identified_as_root++;
+        m_agl[pos].anc_branch_dist = 0;  // Org is root.
+      } else if (m_agl[parent_pos].offspring_count > 1) {        // Parent is branch.
+        genotypes_with_branching_parent++;
+        m_agl[pos].anc_branch_dist = 1;
+        m_agl[pos].anc_branch_id = m_agl[parent_pos].id;
+        m_agl[pos].anc_branch_pos = parent_pos;
+      } else if (m_agl[parent_pos].anc_branch_dist > -1) {     // Parent calculated.
+        genotypes_with_calced_parent++;
+        m_agl[pos].anc_branch_dist = m_agl[parent_pos].anc_branch_dist + 1;
+        m_agl[pos].anc_branch_id = m_agl[parent_pos].anc_branch_id;
+        m_agl[pos].anc_branch_pos = m_agl[parent_pos].anc_branch_pos;
+      } else {
+        genotypes_not_ready++;
+        // Otherwise, we are not yet ready to calculate this entry.
+      }
+    }
+    if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
+      cout << " (" << genotypes_already_set << "," << genotypes_identified_as_root << "," << genotypes_with_branching_parent << "," << genotypes_with_calced_parent << "," << genotypes_not_ready << ")" << endl;
+    }
+    loop_count++;
+  }
+
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Building n-furcating subtree..." << endl;
+  }
+  // compute number of subtree nodes.
+  int branch_tree_size = 0;
+  for (int pos = 0; pos < num_gens; pos++) {
+    //if(m_agl[pos].offspring_count != 1){
+    if(m_agl[pos].offspring_count > 1){
+      branch_tree_size++;
+    }
+  }
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Number of n-furcating nodes: " << branch_tree_size << endl;
+  }
+  if (branch_tree_size <= 0){
+    if (m_world->GetVerbosity() >= VERBOSE_ON) {
+      cerr << "Error: no branches found in tree - " << endl;
+    }
+    return;
+  }
+
+  m_agl2.Resize(branch_tree_size);  // Store agl data for each id.
+  tHashTable<int, int> id_hash_2;
+  int array_pos_2 = 0;
+  if (true) for (int pos = 0; pos < num_gens; pos++) {
+    int offs_count = m_agl[pos].offspring_count;
+    //if (offs_count != 1){
+    if (offs_count > 1){
+      m_agl2[array_pos_2].id = m_agl[pos].id;
+      m_agl2[array_pos_2].pid = m_agl[pos].pid;
+      m_agl2[array_pos_2].depth = m_agl[pos].depth;
+      m_agl2[array_pos_2].birth = m_agl[pos].birth;
+      m_agl2[array_pos_2].anc_branch_dist = m_agl[pos].anc_branch_dist;
+      m_agl2[array_pos_2].anc_branch_id = m_agl[pos].anc_branch_id;
+      /*
+      missing still are ppos (skip this), offspring_count (redundant),
+      anc_branch_pos, off_branch_dist_acc (to be calculated),
+      offspring_positions
+      */
+      id_hash_2.Add(m_agl2[array_pos_2].id, array_pos_2);
+      array_pos_2++;
+    }
+  }
+
+  // find branch ancestor positions. {{{4
+  for (int pos = 0; pos < branch_tree_size; pos++){
+    int anc_branch_id = m_agl2[pos].anc_branch_id;
+    id_hash_2.Find(anc_branch_id, m_agl2[pos].anc_branch_pos);
+    int anc_branch_pos = m_agl2[pos].anc_branch_pos;
+    if(0 <= anc_branch_pos){
+      m_agl2[anc_branch_pos].offspring_positions.Push(pos);
+    }
+  }
+  
+  /*
+  For DFS of branch tree, locate root. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Finding root of n-furcating subtree..." << endl;
+  }
+  cAGLData *root = 0;
+  for (int pos = 0; pos < branch_tree_size; pos++){
+    m_agl2[pos].traversal_visited = false;
+    /*
+    root : anc_branch_dist: 0 anc_branch_id: -1 anc_branch_pos: -1 
+
+    Only one of these conditions should be needed. I'm mixing-in a
+    sanity check here. I ought to move it...
+    */
+    if( (m_agl2[pos].anc_branch_dist == 0)
+      ||(m_agl2[pos].anc_branch_id == -1)
+      ||(m_agl2[pos].anc_branch_pos == -1)
+    ){
+      root = &(m_agl2[pos]);
+      /* Sanity check. */
+      if(
+        !(
+          //(m_agl2[pos].anc_branch_dist == 0) &&
+          (m_agl2[pos].anc_branch_id == -1) &&
+          (m_agl2[pos].anc_branch_pos == -1)
+        )
+      ){
+        if (m_world->GetVerbosity() >= VERBOSE_ON) {
+          cerr << "Error: while looking for root of subtree, found inconsistencies - " << endl;
+          cerr << " root ancestor-branch-distance: " << m_agl2[pos].anc_branch_dist << endl;
+          cerr << " root ancestor-branch-id: " << m_agl2[pos].anc_branch_id << endl;
+          cerr << " root ancestor-branch-pos: " << m_agl2[pos].anc_branch_id << endl;
+        }
+        return;
+      }
+    }
+  }
+
+  /*
+  DFS of branch tree, to accumulate branch distances. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Accumulating branch distances..." << endl;
+  }
+  tList<cAGLData> dfs_stack;
+  if(0 != root){
+    /* DFS. */
+    dfs_stack.Push(root);
+    cAGLData *node = 0;
+    while (0 < dfs_stack.GetSize()){
+      node = dfs_stack.Pop();
+      if (0 != node){
+        if (! node->traversal_visited){
+          dfs_stack.Push(node);
+          node->off_branch_dist_acc = 0;
+          for (int i = 0; i < node->offspring_positions.GetSize(); i++){
+            int pos = node->offspring_positions[i];
+            if (! m_agl2[pos].traversal_visited){
+              dfs_stack.Push(&(m_agl2[pos]));
+            }
+          }
+          node->traversal_visited = true;
+        } else {
+          /*
+          Child nodes, if any, have been visited and have added their
+          off_branch_dist_acc to this node.
+          */
+          if(0 <= node->anc_branch_pos){
+            /*
+            Only accumulate to parent if there is a parent (i.e., this
+            is not the root.)
+            */
+            m_agl2[node->anc_branch_pos].off_branch_dist_acc += node->anc_branch_dist;
+            m_agl2[node->anc_branch_pos].off_branch_dist_acc += node->off_branch_dist_acc;
+          }
+        }
+      }
+    }
+  } else {
+    if (m_world->GetVerbosity() >= VERBOSE_ON) {
+      cerr << "Error: couldn't find root of subtree - " << endl;
+    }
+    return;
+  }
+
+  /*
+  Compute cumulative stemminesses. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Computing cumulative stemminesses..." << endl;
+  }
+  for (int pos = 0; pos < branch_tree_size; pos++){
+    if (0 == m_agl2[pos].anc_branch_dist){
+      /* Correct stemminess for root... */
+      /*
+      Let me rephrase that. If anc_branch_dist is zero and all is well,
+      we're dealing with the root.
+      */
+      m_agl2[pos].cumulative_stemminess = 0.0;
+    } else {
+      m_agl2[pos].cumulative_stemminess =
+      (
+        (double)(m_agl2[pos].anc_branch_dist)
+        /
+        (
+          (double)(m_agl2[pos].off_branch_dist_acc) + (double)(m_agl2[pos].anc_branch_dist)
+        )
+      );
+    }
+  }
+
+  /*
+  Compute average cumulative stemminess. {{{4
+  */
+  if (m_world->GetVerbosity() >= VERBOSE_ON) {
+    cout << "Computing average cumulative stemminess..." << endl;
+  }
+  m_stemminess_sum = 0.0;
+  m_average_stemminess = 0.0;
+  m_inner_nodes = 0;
+  if (1 < branch_tree_size) {
+    for (int pos = 0; pos < branch_tree_size; pos++){
+      bool not_leaf = true;
+      if (m_should_exclude_leaves) {
+        not_leaf = (0 < m_agl2[pos].off_branch_dist_acc);
+      }
+      bool not_root = (0 < m_agl2[pos].anc_branch_id);
+      if (not_leaf && not_root){
+        m_stemminess_sum += m_agl2[pos].cumulative_stemminess;
+        m_inner_nodes++;
+      }
+    }
+  }
+  if(0 < m_inner_nodes){
+    m_average_stemminess = m_stemminess_sum / (m_inner_nodes);
+  }
+
+  /*
+  Print branch tree. {{{4
+  Use to get expected data for making following test.
+  */
+
+  if (false && m_world->GetVerbosity() >= VERBOSE_ON) {
+    for (int pos = 0; pos < branch_tree_size; pos++){
+      cout << "id: " << m_agl2[pos].id;
+      cout << " offspring_count: " << m_agl2[pos].offspring_positions.GetSize();
+      cout << " anc_branch_id: " << m_agl2[pos].anc_branch_id;
+      cout << " anc_branch_dist: " << m_agl2[pos].anc_branch_dist;
+      cout << " off_branch_dist_acc: " << m_agl2[pos].off_branch_dist_acc;
+      cout << " cumulative_stemminess: " << m_agl2[pos].cumulative_stemminess;
+      cout << endl;
+      if (0 < m_agl2[pos].offspring_positions.GetSize()) {
+        cout << "  offspring ids:";
+        for (int i = 0; i < m_agl2[pos].offspring_positions.GetSize(); i++){
+          cout << " " << m_agl2[pos].offspring_positions[i];
+        }
+        cout << endl;
+      }
+    }
+    cout << "m_stemminess_sum: " << m_stemminess_sum << endl;
+    cout << "m_average_stemminess: " << m_average_stemminess << endl;
+    cout << "m_inner_nodes: " << m_inner_nodes << endl;
+  }
+}

Copied: development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.h (from rev 2226, development/source/analyze/cAnalyzeGenotypeTreeStats.h)
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.h	                        (rev 0)
+++ development/source/analyze/cAnalyzeTreeStats_CumulativeStemminess.h	2007-12-19 20:29:15 UTC (rev 2231)
@@ -0,0 +1,101 @@
+/*
+ *  cAnalyzeTreeStats_CumulativeStemminess.h
+ *  Avida at vallista
+ *
+ *  Created by Kaben Nanlohy on 2007.12.03.
+ *  Copyright 1999-2007 Michigan State University. All rights reserved.
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; version 2
+ *  of the License.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ */
+
+#ifndef cAnalyzeTreeStats_CumulativeStemminess_h
+#define cAnalyzeTreeStats_CumulativeStemminess_h
+
+#ifndef tArray_h
+#include "tArray.h"
+#endif
+#ifndef tList_h
+#include "tList.h"
+#endif
+
+
+class cAnalyzeGenotype;
+class cWorld;
+
+
+/*
+Convenience class cAGLData collects data I needed while
+building the class cAnalyzeTreeStats_CumulativeStemminess.
+*/
+class cAGLData {
+public:
+  cAnalyzeGenotype *genotype;
+  int id;
+  int pid;
+  int depth;
+  int birth;
+  int ppos;
+  int offspring_count;
+  int anc_branch_dist;
+  int anc_branch_id;
+  int anc_branch_pos;
+  int off_branch_dist_acc;
+  double cumulative_stemminess;
+  bool traversal_visited;
+  tArray<int> offspring_positions;
+
+  cAGLData();
+};
+
+/*
+class cAnalyzeTreeStats_CumulativeStemminess
+
+This class will eventually collect various kinds of statistics about various
+kinds of trees (e.g., phylogenetic).
+
+For now it only figures out "average cumulative stemminess".
+*/
+class cAnalyzeTreeStats_CumulativeStemminess {
+public:
+  tArray<cAGLData> m_agl;
+  tArray<cAGLData> m_agl2;
+
+  double m_stemminess_sum;
+  double m_average_stemminess;
+  int m_inner_nodes;
+  bool m_should_exclude_leaves;
+
+  cWorld* m_world;
+
+public:
+  cAnalyzeTreeStats_CumulativeStemminess(cWorld* world);
+
+  // Accessors.
+  tArray<cAGLData> &AGL(){ return m_agl; }
+  tArray<cAGLData> &AGL2(){ return m_agl2; }
+  // Getters.
+  double StemminessSum(){ return m_stemminess_sum; }
+  double AverageStemminess(){ return m_average_stemminess; }
+  int InnerNodes(){ return m_inner_nodes; }
+  bool ExcludesLeaves(){ return m_should_exclude_leaves; }
+
+  void PrintAGLData(tArray<cAGLData> &agl);
+
+  // Commands.
+  void AnalyzeBatchTree(tList<cAnalyzeGenotype> &genotype_list);
+};
+
+#endif

Added: development/source/analyze/cAnalyzeTreeStats_Gamma.cc
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_Gamma.cc	                        (rev 0)
+++ development/source/analyze/cAnalyzeTreeStats_Gamma.cc	2007-12-19 20:29:15 UTC (rev 2231)
@@ -0,0 +1,137 @@
+/*
+ *  cAnalyzeTreeStats_Gamma.cc
+ *  Avida at vallista
+ *
+ *  Created by Kaben Nanlohy on 2007.12.03.
+ *  Copyright 1999-2007 Michigan State University. All rights reserved.
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; version 2
+ *  of the License.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ */
+
+#include "cAnalyzeTreeStats_Gamma.h"
+
+#include "cAnalyze.h"
+#include "cAnalyzeGenotype.h"
+#include "tArray.h"
+#include "cWorld.h"
+
+#include <math.h>
+
+cAnalyzeTreeStats_Gamma::cAnalyzeTreeStats_Gamma(cWorld* world)
+: m_world(world)
+, m_gen_array(0)
+, m_g(0)
+, m_gamma(0.)
+{}
+
+void cAnalyzeTreeStats_Gamma::LoadGenotypes(tList<cAnalyzeGenotype> &genotype_list){
+  cAnalyzeGenotype * genotype = NULL;
+  tListIterator<cAnalyzeGenotype> batch_it(genotype_list);
+  
+  m_gen_array.Resize(genotype_list.GetSize());
+  int array_pos = 0;
+  batch_it.Reset();
+  while ((genotype = batch_it.Next()) != NULL) {
+    m_gen_array[array_pos] = genotype;
+    array_pos++;
+  }
+}
+int cAnalyzeTreeStats_Gamma::HeapSortGenotypes(void){
+  return HeapSortAGPhyloDepth(m_gen_array);
+}
+void cAnalyzeTreeStats_Gamma::CalculateInternodeDistances(void){
+  m_g.Resize(1 + m_gen_array.GetSize());
+  m_g[0] = 0;
+  m_g[1] = m_gen_array[0]->GetDepth() - 0;
+  for(int i = 1; i < m_gen_array.GetSize(); i++) {
+    m_g[i+1] = m_gen_array[i]->GetDepth() - m_gen_array[i-1]->GetDepth();
+  }  
+}
+void cAnalyzeTreeStats_Gamma::FixupInternodeDistances(void){
+  bool in_redundant_subsequence = false;
+  int saved_g = -1;
+  for(int i = 1; i < m_gen_array.GetSize(); i++) {
+    /* if we are entering a redundant subsequence, save and then clear first redundant g. */
+    if((m_gen_array[i]->GetDepth() == m_gen_array[i-1]->GetDepth()) && (!in_redundant_subsequence)) {
+      saved_g = m_g[i];
+      m_g[i] = 0;
+      in_redundant_subsequence = true;
+    }
+    
+    /* if we are exiting a redundant subsequence, restore the saved redundant g. */
+    if((m_gen_array[i]->GetDepth() == m_gen_array[i-1]->GetDepth()) && (in_redundant_subsequence)) {
+      m_g[i+1] = saved_g;
+      saved_g = -1;
+      in_redundant_subsequence = false;
+    }
+  }
+  /* if we fell off the end of a redundant subsequence, restore the saved redundant g. */
+  if(in_redundant_subsequence) {
+    m_g[m_gen_array.GetSize()] = saved_g;
+    // saved_g = -1;
+    // in_redundant_subsequence = false;
+  }
+}
+void cAnalyzeTreeStats_Gamma::CalculateGamma(void){
+  // n: number of leaves, constant for a given tree.
+  int n = m_gen_array.GetSize();
+  
+  int T = 0;
+  for(int j = 2; j <= n; j++) { T += j*m_g[j]; }
+  
+  // so: exterior summation
+  int so = 0;
+  for(int i = 2; i <= n-1; i++) {
+    // si: interior summation
+    int si = 0;
+    for(int k = 2; k <= i; k++) { si += k*m_g[k]; }
+    so += si;
+  }
+  
+  m_gamma = ( ( (1./(n-2.)) * so ) - (T/2.) ) / ( T*sqrt( 1. / (12.*(n-2.)) ) );
+}
+double cAnalyzeTreeStats_Gamma::Gamma(void){
+  return m_gamma;
+}
+
+void cAnalyzeTreeStats_Gamma::AnalyzeBatch(tList<cAnalyzeGenotype> &genotype_list){
+  LoadGenotypes(genotype_list);
+  HeapSortGenotypes();
+  CalculateInternodeDistances();
+  CalculateGamma();
+}
+
+
+
+// Comparison function for heapsort.
+int CompareAGPhyloDepth(const void * _a, const void * _b){
+  cAnalyzeGenotype a(**((cAnalyzeGenotype**)_a));
+  cAnalyzeGenotype b(**((cAnalyzeGenotype**)_b));
+  if(a.GetDepth() > b.GetDepth()){ return 1; }
+  if(a.GetDepth() < b.GetDepth()){ return -1; }
+  return 0;
+}
+int HeapSortAGPhyloDepth(tArray<cAnalyzeGenotype *> &gen_array){
+  const int size = gen_array.GetSize();
+  cAnalyzeGenotype *c_gen_array[size];
+  /* Copy unsorted array from gen_array into c_gen_array. */
+  for(int i = 0; i < size; i++){ c_gen_array[i] = (gen_array[i]); }
+  /* Heapsort c_gen_array. */
+  int result = heapsort(c_gen_array, size, sizeof(cAnalyzeGenotype*), CompareAGPhyloDepth);
+  /* If heapsort returned successfully, copy sorted array from c_gen_array into gen_array. */
+  if(result == 0){ for(int i = 0; i < size; i++){ gen_array[i] = (c_gen_array[i]); } }
+  return result;
+}  

Added: development/source/analyze/cAnalyzeTreeStats_Gamma.h
===================================================================
--- development/source/analyze/cAnalyzeTreeStats_Gamma.h	                        (rev 0)
+++ development/source/analyze/cAnalyzeTreeStats_Gamma.h	2007-12-19 20:29:15 UTC (rev 2231)
@@ -0,0 +1,62 @@
+/*
+ *  cAnalyzeTreeStats_Gamma.h
+ *  Avida at vallista
+ *
+ *  Created by Kaben Nanlohy on 2007.12.03.
+ *  Copyright 1999-2007 Michigan State University. All rights reserved.
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; version 2
+ *  of the License.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ */
+
+#ifndef cAnalyzeTreeStats_Gamma_h
+#define cAnalyzeTreeStats_Gamma_h
+
+#ifndef tArray_h
+#include "tArray.h"
+#endif
+#ifndef tList_h
+#include "tList.h"
+#endif
+
+class cAnalyzeGenotype;
+class cWorld;
+
+class cAnalyzeTreeStats_Gamma {
+public:
+  cWorld* m_world;
+  tArray<cAnalyzeGenotype *> m_gen_array;
+  tArray<int> m_g;
+  double m_gamma;
+public:
+  cAnalyzeTreeStats_Gamma(cWorld* world);
+  
+  void LoadGenotypes(tList<cAnalyzeGenotype> &genotype_list);
+  int HeapSortGenotypes(void);
+  void CalculateInternodeDistances(void);
+  void FixupInternodeDistances(void);
+  void CalculateGamma(void);
+  
+  double Gamma(void);
+  
+  // Commands.
+  void AnalyzeBatch(tList<cAnalyzeGenotype> &genotype_list);
+};
+
+// Comparison function for heapsort.
+int CompareAGPhyloDepth(const void * _a, const void * _b);
+int HeapSortAGPhyloDepth(tArray<cAnalyzeGenotype *> &gen_array);
+
+#endif

Modified: extras/source/testsuites/nAnalyze.cpp
===================================================================
--- extras/source/testsuites/nAnalyze.cpp	2007-12-19 20:28:31 UTC (rev 2230)
+++ extras/source/testsuites/nAnalyze.cpp	2007-12-19 20:29:15 UTC (rev 2231)
@@ -3842,6 +3842,8 @@
 }
   cAddTestSuite t("cAnalyze_UnitTest_cAvidaConfigDestruction", testsuite);
 }
+  cAddTestSuite t("cAnalyze_UnitTest_cAvidaConfigDestruction", testsuite);
+}
 
     char * argv[] = {
       "cAnalyze_UnitTest_Destructor",

Modified: extras/source/testsuites/nAnalyzeGenotype.cpp
===================================================================

Modified: extras/source/testsuites/nMemTracker.cpp
===================================================================

Modified: extras/source/testsuites/nString.cpp
===================================================================




More information about the Avida-cvs mailing list