examples/fitness/avida-landscape.cpp
 author Dmitriy Morozov Fri, 06 Nov 2009 14:19:08 -0800 branch dev changeset 173 5fd3f43e6fbf parent 49 7558e122ba4f permissions -rw-r--r--
Added arbitrary coefficients in boundaries (to C++ and Python) + cleaned up CohomologyPersistence.add overloading in Python
```
#include <iostream>
#include <vector>
#include <algorithm>
#include <string>
#include "avida-population-detail.h"

#include <boost/program_options.hpp>

#include <topology/lowerstarfiltration.h>

// Lower-star filtration
typedef         AvidaPopulationDetail::OrganismIndex                OrganismIndex;
struct          OrganismVertexType;
typedef         std::vector<OrganismVertexType>                     VertexVector;
typedef         VertexVector::iterator                              VertexIndex;
typedef         LowerStarFiltration<VertexIndex>                    LSFiltration;
typedef         LSFiltration::Simplex                               Simplex;

struct          OrganismVertexType: public LSFiltration::VertexType<OrganismIndex>
{
typedef     LSFiltration::VertexType<OrganismIndex>             Parent;
OrganismVertexType(OrganismIndex i): Parent(i)      {}
};

struct          OrganismVertexComparison
{
public:
bool    operator()(VertexIndex i, VertexIndex j) const
{ return i->index()->fitness() > j->index()->fitness(); }
// > because of -fitness, so that maxima turn into minima
};

typedef         LSFiltration::Vineyard                              LSVineyard;
class           StaticEvaluator: public LSVineyard::Evaluator
{
public:
StaticEvaluator(float max_fitness):
max_fitness_(max_fitness)                       {}

virtual RealType
value(const Simplex& s) const
{ return s.get_attachment()->index()->fitness()/max_fitness_; }

private:
float   max_fitness_;
};

std::ostream& operator<<(std::ostream& out, VertexIndex i)
{ return (out << (i->index())); }

// Distance filtration
typedef         SimplexWithValue<VertexIndex>                       DistanceSimplex;
typedef         std::vector<DistanceSimplex>                        DistanceSimplexVector;
typedef         Filtration<DistanceSimplex>                         DistanceSimplexFiltration;

namespace po = boost::program_options;

int main(int argc, char** argv)
{
#ifdef LOGGING
stderrLog.subscribeTo(RLOG_CHANNEL("error"));
//stdoutLog.subscribeTo(RLOG_CHANNEL("info"));
#endif

typedef     AvidaOrganismDetail::DistanceType       DistanceType;

DistanceType connected_distance;
bool connect_mst = false;
std::string population_input_fn;

// Parse program options
po::options_description hidden("Hidden options");
("input-file",  po::value<std::string>(&population_input_fn),
"Avida population file");

po::options_description visible("Allowed options");
("help,h",      "produce help message")
("distance,d",  po::value<DistanceType>(&connected_distance)->default_value(0),
"set connected distance")
("mst,m",       "connect minimum spanning tree");
po::positional_options_description p;

po::variables_map vm;
po::store(po::command_line_parser(argc, argv).
options(all).positional(p).run(), vm);
po::notify(vm);

if (vm.count("mst"))            { connect_mst = true; }
if (vm.count("help") || !vm.count("input-file"))
{
std::cout << "Usage: " << argv[0] << " [options] POPULATION" << std::endl;
std::cout << visible << std::endl;
return 1;
}

AvidaPopulationDetail population(population_input_fn);
const AvidaPopulationDetail::OrganismVector& organisms = population.get_organisms();

rInfo("Number of organisms: %d", organisms.size());
float max_fitness = organisms[0].fitness();
for (int i = 0; i < population.get_organisms().size(); ++i)
{
max_fitness = std::max(max_fitness, organisms[i].fitness());
rInfo("%d (%s) %f %d %d", organisms[i].id(),
organisms[i].genome().c_str(),
organisms[i].fitness(),
organisms[i].length(),
organisms[i].genome().size());
}

// Order vertices
StaticEvaluator     evaluator(max_fitness);
LSVineyard          vineyard(&evaluator);
VertexVector        vertices;
for (OrganismIndex cur = organisms.begin(); cur != organisms.end(); ++cur)  vertices.push_back(cur);
LSFiltration        fitness_filtration(vertices.begin(), vertices.end(), OrganismVertexComparison(), &vineyard);

// Compute MST and insert its edges if asked
if (connect_mst)
{
DistanceSimplexFiltration filtration;
{   // Scope so that simplices is deleted once it's not needed
// Distance function filtration
DistanceSimplexVector simplices;

// Insert edges
for (VertexIndex i = vertices.begin(); i != vertices.end(); ++i)
{
simplices.push_back(DistanceSimplex());

for (VertexIndex j = boost::next(i); j != vertices.end(); ++j)
{
simplices.push_back(DistanceSimplex(i->index()->genome_distance(*(j->index()))));
}
}
std::sort(simplices.begin(), simplices.end(), DimensionValueComparison<DistanceSimplex>());

for (DistanceSimplexVector::const_iterator  cur = simplices.begin();
cur != simplices.end(); ++cur)
filtration.append(*cur);
}

filtration.fill_simplex_index_map();
filtration.pair_simplices(false);            // pair simplices without storing trails

for (DistanceSimplexFiltration::Index i = filtration.begin(); i != filtration.end(); ++i)
{
if (i->is_paired() && !i->sign())
{
Simplex s(*i);
if (i->get_value() > connected_distance)    // <= will be connected below
fitness_filtration.append(s);
}
}
}

for (VertexIndex cur = vertices.begin(); cur != vertices.end(); ++cur)
if (cur->index()->genome_distance(*(link->index())) <= connected_distance)
{
fitness_filtration.append(s);
}
rInfo("Number of simplices: %d", fitness_filtration.size());

// Pair simplices
fitness_filtration.fill_simplex_index_map();
fitness_filtration.pair_simplices(false);            // pair simplices without storing trails

//std::cout << "Outputting persistence pairs" << std::endl;
for (LSFiltration::Index i = fitness_filtration.begin(); i != fitness_filtration.end(); ++i)
{
if (i->is_paired())
{
if (i->sign())
{
AssertMsg(i->dimension() == 0, "Expecting only 0-dimensional diagram");
if (i->pair()->get_attachment() == i->vertices()[0]) continue;     // skip non-critical pairs
std::cout << i->dimension() << " "
<< evaluator.value(*i) << " "
<< evaluator.value(*(i->pair())) << std::endl;
}
}
else
{
if (i->dimension() != 0) continue;
std::cout << i->dimension() << " "
<< evaluator.value(*i) << " "
<< 0 << std::endl;
// The infinite pair does not make sense since we are interested in
// max of fitness, rather than min of -fitness. However min value for
// fitness is 0, so it's a natural choice of the answer.
}
}
}
```