Made storing negative simplices in StaticPersistence parameter-based, so that DynamicPersistence works again
#include <iostream>
#include <vector>
#include <algorithm>
#include "avida-population-detail.h"
#include <topology/filtration.h>
#include <topology/simplex.h>
#include <topology/static-persistence.h>
typedef Simplex<AvidaOrganismDetail::IDType, double> Smplx;
typedef std::vector<Smplx> Complex;
typedef Filtration<Complex> Fltr;
typedef StaticPersistence<> Persistence;
int main(int argc, char** argv)
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
//stdoutLog.subscribeTo(RLOG_CHANNEL("info"));
#endif
if (argc < 2)
{
std::cout << "USAGE: avida FILENAME" << std::endl;
return 0;
}
AvidaPopulationDetail population(argv[1]);
const AvidaPopulationDetail::OrganismVector& organisms = population.get_organisms();
rInfo("Number of organisms: %d", organisms.size());
for (int i = 0; i < population.get_organisms().size(); ++i)
rInfo("%d (%s) %f %d %d", organisms[i].id(),
organisms[i].genome().c_str(),
organisms[i].fitness(),
organisms[i].length(),
organisms[i].genome().size());
// Distance function filtration
Complex simplices;
// Insert edges between all the organisms
AvidaOrganismDetail::DistanceType avg_distance = 0;
for (AvidaOrganismDetail::CountType i = 0; i < organisms.size(); ++i)
{
simplices.push_back(0); // all vertices have 0 value
simplices.back().add(organisms[i].id());
for (AvidaOrganismDetail::CountType j = i+1; j < organisms.size(); ++j)
{
avg_distance += organisms[i].genome_distance(organisms[j]);
simplices.push_back(Smplx(organisms[i].genome_distance(organisms[j])));
simplices.back().add(organisms[i].id());
simplices.back().add(organisms[j].id());
}
}
std::sort(simplices.begin(), simplices.end(), Smplx::VertexComparison());
rInfo("Average distance: %f", float(avg_distance)/
((organisms.size()*organisms.size() - organisms.size())/2));
Fltr f(simplices.begin(), simplices.end(), DataDimensionComparison<Smplx>());
Persistence p(f);
p.pair_simplices();
std::cout << "Outputting histogram of death values" << std::endl;
typedef std::vector<RealType> DeathVector;
DeathVector deaths;
Smplx::DataEvaluator eval;
OffsetMap<Persistence::OrderIndex, Fltr::Index> m(p.begin(), f.begin());
for (Persistence::OrderIndex i = p.begin(); i != p.end(); ++i)
{
if (i == i->pair) continue;
if (i->sign())
{
const Smplx& s = f.simplex(m[i]);
const Smplx& t = f.simplex(m[i->pair]);
AssertMsg(s.dimension() == 0, "Expecting only 0-dimensional diagram");
AssertMsg(eval(s) == 0, "Expecting only 0 birth values in 0-D diagram ");
deaths.push_back(eval(t));
}
}
// Produce histogram
std::sort(deaths.begin(), deaths.end());
for (DeathVector::iterator cur = deaths.begin(); cur != deaths.end(); )
{
DeathVector::iterator nw = std::find_if(cur, deaths.end(),
std::bind2nd(std::greater<RealType>(), *cur));
std::cout << *cur << "\t" << (nw - cur) << std::endl;
cur = nw;
}
std::cout << "Total: " << deaths.size() + 1; // +1 for the unpaired
}