# HG changeset patch
# User Dmitriy Morozov <morozov@cs.duke.edu>
# Date 1222145596 25200
# Node ID 3f94bf5ba9893032b94ebc8b7c09bad4b28ab33c
# Parent  900f42b643dec1b9394f9d5a3e5faed501e56224
Added Rips complex-based distance computation to fitness

diff -r 900f42b643de -r 3f94bf5ba989 examples/fitness/CMakeLists.txt
--- a/examples/fitness/CMakeLists.txt	Mon Sep 22 21:49:22 2008 -0700
+++ b/examples/fitness/CMakeLists.txt	Mon Sep 22 21:53:16 2008 -0700
@@ -1,5 +1,6 @@
 set							(targets						
 							 avida-distance
+							 avida-rips-distance
                              avida-landscape)
 	
 foreach 					(t ${targets})
diff -r 900f42b643de -r 3f94bf5ba989 examples/fitness/avida-population-detail.h
--- a/examples/fitness/avida-population-detail.h	Mon Sep 22 21:49:22 2008 -0700
+++ b/examples/fitness/avida-population-detail.h	Mon Sep 22 21:53:16 2008 -0700
@@ -53,7 +53,20 @@
 
 						AvidaPopulationDetail(std::string filename);
 
-		const OrganismVector& get_organisms() const						{ return organisms_; }
+		const OrganismVector& 
+                        get_organisms() const						    { return organisms_; }
+
+        /// \name Rips
+        /// @{
+        typedef         int                                             IndexType;
+        typedef         double                                          DistanceType;
+
+        DistanceType    operator()(IndexType a, IndexType b) const      { return organisms_[a].genome_distance(organisms_[b]); }
+
+        size_t          size() const                                    { return organisms_.size(); }
+        IndexType       begin() const                                   { return 0; }
+        IndexType       end() const                                     { return size(); }
+        /// @}
 
 	private:
 		OrganismVector	organisms_;
diff -r 900f42b643de -r 3f94bf5ba989 examples/fitness/avida-rips-distance.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/fitness/avida-rips-distance.cpp	Mon Sep 22 21:53:16 2008 -0700
@@ -0,0 +1,87 @@
+#include <iostream>
+#include <vector>
+#include <algorithm>
+#include "avida-population-detail.h"
+
+#include <topology/filtration.h>
+#include <topology/simplex.h>
+#include <topology/rips.h>
+
+
+typedef         ExplicitDistances<AvidaPopulationDetail>            ExplicitDist;
+typedef         Rips<ExplicitDist>                                  RipsComplex; 
+typedef         RipsComplex::Simplex                                Simplex;
+typedef         RipsComplex::SimplexVector                          SimplexVector;
+typedef			Filtration<Simplex>									SimplexFiltration;
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
+	stdoutLog.subscribeTo(RLOG_CHANNEL("info"));
+	//stdoutLog.subscribeTo(RLOG_CHANNEL("rips/info"));
+#endif
+
+	if (argc < 2)
+	{
+		std::cout << "USAGE: avida FILENAME" << std::endl;
+		return 0;
+	}
+
+	AvidaPopulationDetail population(argv[1]);
+    ExplicitDist distances(population);
+
+    RipsComplex rips(distances);
+    RipsComplex::Evaluator evaluator(rips.distances());
+    rInfo("Max distance: %f", rips.max_distance());
+
+	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());
+   */
+
+    rInfo("Starting to generate rips complex");
+    rips.generate(1, rips.max_distance()/2);
+    
+    rInfo("Generated Rips complex, filling filtration");
+	SimplexFiltration filtration;
+	for (SimplexVector::const_iterator cur = rips.simplices().begin(); cur != rips.simplices().end(); ++cur)
+	{
+		//rInfo("Simplex: %s %f", tostring(*cur).c_str(), evaluator.value(*cur));
+		filtration.append(*cur);
+	}
+
+	filtration.fill_simplex_index_map();
+	filtration.pair_simplices(false);			// pair simplices without storing trails
+
+	std::cout << "Outputting histogram of death values" << std::endl;
+	typedef std::vector<RealType> DeathVector;
+	DeathVector deaths;
+	for (SimplexFiltration::Index i = filtration.begin(); i != filtration.end(); ++i)
+	{
+		if (i->is_paired())
+			if (i->sign())
+			{
+				AssertMsg(i->dimension() == 0, "Expecting only 0-dimensional diagram");
+				AssertMsg(evaluator.value(*i) == 0, "Expecting only 0 birth values in 0-D diagram ");
+				deaths.push_back(evaluator.value(*(i->pair())));
+			}
+	}
+
+	// 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
+}