Commit before merging in Python branch:
* converted more examples: poincare, rips, some fitness
* two generators for Rips complexes
* ChainWrapper uses stl algorithms for everything,
added CountingBackInserter and switched SizeStorage to use operators
* retabbing files along the way
* added #957a (namespace dionysus)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/957a589d7c6c3fa8/new/1230356422.M300560P11955Q1.rufus Sat Dec 27 14:33:25 2008 -0800
@@ -0,0 +1,10 @@
+From: Dmitriy Morozov <dmitriy@mrzv.org>
+Date: Fri, 26 Dec 2008 21:34:29
+State: new
+Subject: namespace dionysus
+Message-Id: <957a589d7c6c3fa8-0-artemis@rufus>
+
+Put everything in namespace dionysus. Besides making things more organized, it
+will also allow us to get rid of the problem of having to typedef Simplex<...>
+to Smplx. If it were in dionysus namespace, then
+`typedef dionysus::Simplex<...> Simplex;` would work.
--- a/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,107 +1,107 @@
-project (Dionysus)
+project (Dionysus)
cmake_minimum_required (VERSION 2.4)
-option (logging "Build Dionysus with logging on" OFF)
-option (counters "Build Dionysus with counters on" OFF)
-option (debug "Build Dionysus with debugging on" OFF)
-option (optimize "Build Dionysus with optimization" ON)
-option (use_dsrpdb "Build examples that use DSR-PDB" OFF)
-option (use_synaps "Build examples that use SYNAPS" OFF)
+option (logging "Build Dionysus with logging on" OFF)
+option (counters "Build Dionysus with counters on" OFF)
+option (debug "Build Dionysus with debugging on" OFF)
+option (optimize "Build Dionysus with optimization" ON)
+option (use_dsrpdb "Build examples that use DSR-PDB" OFF)
+option (use_synaps "Build examples that use SYNAPS" OFF)
# Find everything that's always required
-find_package (Boost REQUIRED COMPONENTS program_options serialization signals)
-find_package (Doxygen)
+find_package (Boost REQUIRED COMPONENTS program_options serialization signals)
+find_package (Doxygen)
if (use_dsrpdb)
- find_library (dsrpdb_LIBRARY NAMES dsrpdb)
- find_path (dsrpdb_INCLUDE_DIR dsrpdb/Protein.h)
+ find_library (dsrpdb_LIBRARY NAMES dsrpdb)
+ find_path (dsrpdb_INCLUDE_DIR dsrpdb/Protein.h)
endif (use_dsrpdb)
#CGAL
-execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile libpaths
- OUTPUT_VARIABLE cgal_libpaths)
-execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile ldflags
- OUTPUT_VARIABLE cgal_ldflags)
-execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile cxxflags
- OUTPUT_VARIABLE cgal_cxxflags)
-execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile libpath
- OUTPUT_VARIABLE cgal_libpath)
-#string (REPLACE "\n" "" cgal_libpaths ${cgal_libpaths})
-#string (REPLACE "\n" "" cgal_ldflags ${cgal_ldflags})
-string (REPLACE "\n" "" cgal_cxxflags ${cgal_cxxflags})
-string (REPLACE "\n" "" cgal_libpath ${cgal_libpath})
-find_library (cgal_LIBRARY NAMES CGAL
- PATHS ${cgal_libpath})
-find_library (core_LIBRARY NAMES CGALcore++
- PATHS ${cgal_libpath})
-find_library (mpfr_LIBRARY NAMES mpfr)
-find_library (gmp_LIBRARY NAMES gmp)
-find_library (gmpxx_LIBRARY NAMES gmpxx)
-find_library (m_LIBRARY NAMES m)
+execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile libpaths
+ OUTPUT_VARIABLE cgal_libpaths)
+execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile ldflags
+ OUTPUT_VARIABLE cgal_ldflags)
+execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile cxxflags
+ OUTPUT_VARIABLE cgal_cxxflags)
+execute_process (COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile libpath
+ OUTPUT_VARIABLE cgal_libpath)
+#string (REPLACE "\n" "" cgal_libpaths ${cgal_libpaths})
+#string (REPLACE "\n" "" cgal_ldflags ${cgal_ldflags})
+string (REPLACE "\n" "" cgal_cxxflags ${cgal_cxxflags})
+string (REPLACE "\n" "" cgal_libpath ${cgal_libpath})
+find_library (cgal_LIBRARY NAMES CGAL
+ PATHS ${cgal_libpath})
+find_library (core_LIBRARY NAMES CGALcore++
+ PATHS ${cgal_libpath})
+find_library (mpfr_LIBRARY NAMES mpfr)
+find_library (gmp_LIBRARY NAMES gmp)
+find_library (gmpxx_LIBRARY NAMES gmpxx)
+find_library (m_LIBRARY NAMES m)
-set (cgal_libraries ${cgal_LIBRARY}
- ${core_LIBRARY}
- ${mpfr_LIBRARY}
- ${gmp_LIBRARY}
- ${gmpxx_LIBRARY}
- ${m_LIBRARY})
+set (cgal_libraries ${cgal_LIBRARY}
+ ${core_LIBRARY}
+ ${mpfr_LIBRARY}
+ ${gmp_LIBRARY}
+ ${gmpxx_LIBRARY}
+ ${m_LIBRARY})
add_definitions (-DCGAL_NO_ASSERTIONS -DCGAL_NO_PRECONDITIONS)
# SYNAPS
if (use_synaps)
- add_definitions (-DBOOST_UBLAS_TYPE_CHECK=0)
- find_library (synaps_LIBRARY NAMES synaps)
- set (synaps_libraries ${synaps_LIBRARY}
- ${gmp_LIBRARY}
- ${gmpxx_LIBRARY})
+ add_definitions (-DBOOST_UBLAS_TYPE_CHECK=0)
+ find_library (synaps_LIBRARY NAMES synaps)
+ set (synaps_libraries ${synaps_LIBRARY}
+ ${gmp_LIBRARY}
+ ${gmpxx_LIBRARY})
endif (use_synaps)
# Debugging
-if (debug)
- if (optimize)
- set (cxx_flags ${CMAKE_CXX_FLAGS_RELWITHDEBINFO})
- else (optimize)
- set (cxx_flags ${CMAKE_CXX_FLAGS_DEBUG})
- endif (optimize)
-else (debug)
- if (optimize)
- set (cxx_flags ${CMAKE_CXX_FLAGS_RELEASE})
- else (optimize)
- set (cxx_flags ${CMAKE_CXX_FLAGS})
- endif (optimize)
-endif (debug)
-add_definitions (${cxx_flags})
+if (debug)
+ if (optimize)
+ set (cxx_flags ${CMAKE_CXX_FLAGS_RELWITHDEBINFO})
+ else (optimize)
+ set (cxx_flags ${CMAKE_CXX_FLAGS_DEBUG})
+ endif (optimize)
+else (debug)
+ if (optimize)
+ set (cxx_flags ${CMAKE_CXX_FLAGS_RELEASE})
+ else (optimize)
+ set (cxx_flags ${CMAKE_CXX_FLAGS})
+ endif (optimize)
+endif (debug)
+add_definitions (${cxx_flags})
# Logging
-if (logging)
- find_library (rlog_LIBRARY NAMES rlog)
- find_path (rlog_INCLUDE_DIR rlog/rlog.h)
- set (rlog_INCLUDE_DIR ${rlog_INCLUDE_DIR})
- add_definitions (-DLOGGING -DRLOG_COMPONENT=dionysus)
- set (libraries ${libraries} ${rlog_LIBRARY})
-endif (logging)
+if (logging)
+ find_library (rlog_LIBRARY NAMES rlog)
+ find_path (rlog_INCLUDE_DIR rlog/rlog.h)
+ set (rlog_INCLUDE_DIR ${rlog_INCLUDE_DIR})
+ add_definitions (-DLOGGING -DRLOG_COMPONENT=dionysus)
+ set (libraries ${libraries} ${rlog_LIBRARY})
+endif (logging)
# Counters
-if (counters)
- add_definitions (-DCOUNTERS)
-endif (counters)
+if (counters)
+ add_definitions (-DCOUNTERS)
+endif (counters)
# Set includes
-include_directories (${CMAKE_CURRENT_BINARY_DIR}
- ${CMAKE_CURRENT_SOURCE_DIR}/include
- ${Boost_INCLUDE_DIR}
- ${dsrpdb_INCLUDE_DIR}
- ${cwd_INCLUDE_DIR}
- ${rlog_INCLUDE_DIR})
+include_directories (${CMAKE_CURRENT_BINARY_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/include
+ ${Boost_INCLUDE_DIR}
+ ${dsrpdb_INCLUDE_DIR}
+ ${cwd_INCLUDE_DIR}
+ ${rlog_INCLUDE_DIR})
# Doxygen (FIXME)
-if (DOXYGEN_FOUND)
-# add_custom_target (docs ALL
-# ${DOXYGEN_EXECUTABLE} Doxyfile
-# DEPENDS Doxyfile)
-endif (DOXYGEN_FOUND)
+if (DOXYGEN_FOUND)
+# add_custom_target (docs ALL
+# ${DOXYGEN_EXECUTABLE} Doxyfile
+# DEPENDS Doxyfile)
+endif (DOXYGEN_FOUND)
# Process subdirectories
-add_subdirectory (examples)
-add_subdirectory (tests)
-add_subdirectory (tools)
+add_subdirectory (examples)
+add_subdirectory (tests)
+add_subdirectory (tools)
--- a/examples/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,8 +1,8 @@
add_subdirectory (alphashapes)
-add_subdirectory (ar-vineyard)
+#add_subdirectory (ar-vineyard)
add_subdirectory (cech-complex)
add_subdirectory (fitness)
-add_subdirectory (grid)
+#add_subdirectory (grid)
add_subdirectory (triangle)
add_subdirectory (poincare)
add_subdirectory (rips)
--- a/examples/alphashapes/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/alphashapes/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,15 +1,15 @@
-set (targets
- alphashapes3d
- alphashapes2d
- alpharadius
+set (targets
+ alphashapes3d
+ #alphashapes2d
+ #alpharadius
compare-diagrams)
-
+
set (libraries ${libraries}
${Boost_SERIALIZATION_LIBRARY}
${Boost_PROGRAM_OPTIONS_LIBRARY})
-add_definitions (${cgal_cxxflags})
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${cgal_libraries})
-endforeach (t ${targets})
+add_definitions (${cgal_cxxflags})
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${cgal_libraries})
+endforeach (t ${targets})
--- a/examples/cech-complex/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/cech-complex/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,7 +1,7 @@
-set (targets
- cech-complex)
-
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries})
-endforeach (t ${targets})
+set (targets
+ cech-complex)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries})
+endforeach (t ${targets})
--- a/examples/cech-complex/cech-complex.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/cech-complex/cech-complex.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -111,7 +111,7 @@
std::sort(sv.begin(), sv.end(), Smplx::VertexComparison());
// Set up the filtration
- CechFiltration cf(sv.begin(), sv.end(), Smplx::DataDimensionComparison());
+ CechFiltration cf(sv.begin(), sv.end(), DataDimensionComparison<Smplx>());
rInfo("Filtration initialized");
// Compute persistence
--- a/examples/fitness/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/fitness/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,9 +1,10 @@
-set (targets
- avida-distance
- avida-rips-distance
- avida-landscape)
-
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-endforeach (t ${targets})
+set (targets
+ avida-distance
+ avida-rips-distance
+ #avida-landscape
+ )
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endforeach (t ${targets})
--- a/examples/fitness/avida-distance.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/fitness/avida-distance.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -5,90 +5,90 @@
#include <topology/filtration.h>
#include <topology/simplex.h>
+#include <topology/static-persistence.h>
-typedef SimplexWithValue<AvidaOrganismDetail::IDType> Simplex;
-typedef std::vector<Simplex> SimplexVector;
-typedef Filtration<Simplex> SimplexFiltration;
+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"));
+ rlog::RLogInit(argc, argv);
+ //stdoutLog.subscribeTo(RLOG_CHANNEL("info"));
#endif
- if (argc < 2)
- {
- std::cout << "USAGE: avida FILENAME" << std::endl;
- return 0;
- }
+ if (argc < 2)
+ {
+ std::cout << "USAGE: avida FILENAME" << std::endl;
+ return 0;
+ }
- AvidaPopulationDetail population(argv[1]);
- const AvidaPopulationDetail::OrganismVector& organisms = population.get_organisms();
+ 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());
+ 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
- SimplexVector simplices;
+ // Distance function filtration
+ Complex simplices;
- // Insert edges
- AvidaOrganismDetail::DistanceType avg_distance = 0;
- for (AvidaOrganismDetail::CountType i = 0; i < organisms.size(); ++i)
- {
- simplices.push_back(0);
- simplices.back().add(organisms[i].id());
+ // 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(Simplex(organisms[i].genome_distance(organisms[j])));
- simplices.back().add(organisms[i].id());
- simplices.back().add(organisms[j].id());
- }
- }
- std::sort(simplices.begin(), simplices.end(), DimensionValueComparison<Simplex>());
- rInfo("Average distance: %f", float(avg_distance)/
- ((organisms.size()*organisms.size() - organisms.size())/2));
+ 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));
- SimplexFiltration filtration;
- for (SimplexVector::const_iterator cur = simplices.begin(); cur != simplices.end(); ++cur)
- {
- rInfo("Simplex: %s", tostring(*cur).c_str());
- filtration.append(*cur);
- }
-
- filtration.fill_simplex_index_map();
- filtration.pair_simplices(false); // pair simplices without storing trails
+ 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;
- 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(i->get_value() == 0, "Expecting only 0 birth values in 0-D diagram ");
- deaths.push_back(i->pair()->get_value());
- }
- }
+ 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
+ // 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
}
--- a/examples/fitness/avida-rips-distance.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/fitness/avida-rips-distance.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -4,84 +4,86 @@
#include "avida-population-detail.h"
#include <topology/filtration.h>
-#include <topology/simplex.h>
#include <topology/rips.h>
+#include <topology/static-persistence.h>
typedef ExplicitDistances<AvidaPopulationDetail> ExplicitDist;
-typedef Rips<ExplicitDist> RipsComplex;
-typedef RipsComplex::Simplex Simplex;
-typedef RipsComplex::SimplexVector SimplexVector;
-typedef Filtration<Simplex> SimplexFiltration;
+typedef RipsGenerator<ExplicitDist> RipsGen;
+typedef RipsGen::Simplex Smplx;
+typedef RipsGen::SimplexVector Complex;
+
+typedef Filtration<Complex, unsigned> Fltr;
+typedef StaticPersistence<> Persistence;
int main(int argc, char** argv)
{
#ifdef LOGGING
- rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo(RLOG_CHANNEL("info"));
- //stdoutLog.subscribeTo(RLOG_CHANNEL("rips/info"));
+ 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;
- }
+ if (argc < 2)
+ {
+ std::cout << "USAGE: avida FILENAME" << std::endl;
+ return 0;
+ }
- AvidaPopulationDetail population(argv[1]);
+ AvidaPopulationDetail population(argv[1]);
ExplicitDist distances(population);
- RipsComplex rips(distances);
- RipsComplex::Evaluator evaluator(rips.distances());
+ RipsGen rips(distances);
+ RipsGen::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());
+ 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());
+ 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);
+ Complex c;
+ rips.generate(c, 1, rips.max_distance()/2);
+ std::sort(c.begin(), c.end(), Smplx::VertexComparison());
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);
- }
+ Fltr f(c.begin(), c.end(), RipsGen::Comparison(rips.distances()));
- filtration.fill_simplex_index_map();
- filtration.pair_simplices(false); // pair simplices without storing trails
+ Persistence p(f);
+ p.pair_simplices();
- 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())));
- }
- }
+ std::cout << "Outputting histogram of death values" << std::endl;
+ typedef std::vector<RealType> DeathVector;
+ DeathVector deaths;
+ 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(evaluator(s) == 0, "Expecting only 0 birth values in 0-D diagram ");
+ deaths.push_back(evaluator(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
+ // 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
}
--- a/examples/poincare/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/poincare/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,7 +1,7 @@
-set (targets
- poincare)
-
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-endforeach (t ${targets})
+set (targets
+ poincare)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endforeach (t ${targets})
--- a/examples/poincare/poincare.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/poincare/poincare.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -1,13 +1,19 @@
+#include "topology/simplex.h"
#include "topology/filtration.h"
-#include "topology/simplex.h"
+#include "topology/static-persistence.h"
+#include "topology/persistence-diagram.h"
+
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <boost/program_options.hpp>
-typedef SimplexWithValue<int> Simplex;
-typedef Filtration<Simplex> SimplexFiltration;
+typedef Simplex<unsigned, unsigned> Smplx;
+typedef std::vector<Smplx> Complex;
+typedef Filtration<Complex> Fltr;
+typedef StaticPersistence<> Persistence;
+typedef PersistenceDiagram<> PDgm;
namespace po = boost::program_options;
@@ -42,9 +48,7 @@
}
- Evaluator<Simplex> e;
- SimplexFiltration::Vineyard v(&e);
- SimplexFiltration f(&v);
+ Complex c;
std::ifstream in(infilename.c_str());
unsigned int i = 0;
@@ -53,7 +57,7 @@
while(in)
{
std::istringstream linestream(s);
- Simplex simplex(float(i++));
+ Smplx simplex(i++);
unsigned int vertex;
linestream >> vertex;
while(linestream)
@@ -62,19 +66,25 @@
linestream >> vertex;
}
std::cout << simplex << std::endl;
- f.append(simplex);
+ c.push_back(simplex);
std::getline(in, s);
}
-
- f.fill_simplex_index_map();
- f.pair_simplices();
- v.start_vines(f.begin(), f.end());
-
- std::cout << "Filtration size: " << f.size() << std::endl;
- for (SimplexFiltration::Index cur = f.begin(); cur != f.end(); ++cur)
- if (cur->sign())
- std::cout << cur->dimension() << " "
- << cur->get_value() << " "
- << cur->pair()->get_value() << std::endl;
+
+ std::sort(c.begin(), c.end(), Smplx::VertexComparison());
+ Fltr f(c.begin(), c.end(), Smplx::DataComparison());
+ Persistence pers(f);
+ pers.pair_simplices();
+
+ std::map<Dimension, PDgm> dgms;
+ init_diagrams(dgms, pers.begin(), pers.end(),
+ evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(pers.begin(), f.begin()),
+ evaluate_through_filtration(f, Smplx::DataEvaluator())),
+ evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(pers.begin(), f.begin()),
+ evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+
+ std::cout << 0 << std::endl << dgms[0] << std::endl;
+ std::cout << 1 << std::endl << dgms[1] << std::endl;
+ std::cout << 2 << std::endl << dgms[2] << std::endl;
+ std::cout << 3 << std::endl << dgms[3] << std::endl;
}
--- a/examples/rips/CMakeLists.txt Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/rips/CMakeLists.txt Sat Dec 27 14:33:25 2008 -0800
@@ -1,7 +1,7 @@
-set (targets
- rips)
-
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries})
-endforeach (t ${targets})
+set (targets
+ rips)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries})
+endforeach (t ${targets})
--- a/examples/rips/rips.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/rips/rips.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -24,15 +24,20 @@
Distances distances;
#if 0
+ typedef RipsGenerator<ExplicitDistances<Distances> > RipsGenerator;
ExplicitDistances<Distances> explicit_distances(distances);
- Rips<ExplicitDistances<Distances> > rips(explicit_distances);
+ RipsGenerator rips(explicit_distances);
#else
- Rips<Distances> rips(distances);
+ typedef RipsGeneratorMemory<Distances> RipsGenerator;
+ RipsGenerator rips(distances);
#endif
- //rips.generate(3, distances.size());
- rips.generate(3, 50);
+ RipsGenerator::SimplexVector complex;
+ //rips.generate(complex, 3, distances.size());
+ rips.generate(complex, 3, 50);
//rips.print();
- std::cout << "Size: " << rips.size() << std::endl;
+ std::cout << "Size: " << complex.size() << std::endl;
+// for (RipsGenerator::SimplexVector::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+// std::cout << *cur << std::endl;
}
--- a/examples/triangle/triangle.cpp Thu Dec 25 14:24:02 2008 -0800
+++ b/examples/triangle/triangle.cpp Sat Dec 27 14:33:25 2008 -0800
@@ -78,7 +78,7 @@
}
#endif
- Fltr f(c.begin(), c.end(), Smplx::DataDimensionComparison());
+ Fltr f(c.begin(), c.end(), Smplx::DataComparison());
std::cout << "Filtration initialized" << std::endl;
std::cout << f << std::endl;
@@ -90,9 +90,9 @@
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(p.begin(), f.begin()),
+ evaluate_through_map(make_offset_map(p.begin(), f.begin()),
evaluate_through_filtration(f, Smplx::DataEvaluator())),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(p.begin(), f.begin()),
+ evaluate_through_map(make_offset_map(p.begin(), f.begin()),
evaluate_through_filtration(f, Smplx::DimensionExtractor())));
std::cout << 0 << std::endl << dgms[0] << std::endl;
--- a/include/topology/chain.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/chain.h Sat Dec 27 14:33:25 2008 -0800
@@ -31,87 +31,87 @@
class ChainWrapper: public Container_,
public SizeStorage<Container_>
{
- public:
- /// \name Template Parameters
- /// @{
+ public:
+ /// \name Template Parameters
+ /// @{
typedef Container_ Container;
typedef typename boost::iterator_value<typename Container::iterator>::type Item;
- /// @}
-
- typedef ChainWrapper Self;
- typedef Container ChainRepresentation;
+ /// @}
+
+ typedef ChainWrapper Self;
+ typedef Container ChainRepresentation;
- /// \name Accessor typedefs
- /// @{
- typedef typename ChainRepresentation::iterator iterator;
- typedef typename ChainRepresentation::const_iterator const_iterator;
- typedef typename ChainRepresentation::const_reference const_reference;
- typedef typename ChainRepresentation::reference reference;
- typedef typename ChainRepresentation::pointer pointer;
- typedef typename ChainRepresentation::const_pointer const_pointer;
+ /// \name Accessor typedefs
+ /// @{
+ typedef typename ChainRepresentation::iterator iterator;
+ typedef typename ChainRepresentation::const_iterator const_iterator;
+ typedef typename ChainRepresentation::const_reference const_reference;
+ typedef typename ChainRepresentation::reference reference;
+ typedef typename ChainRepresentation::pointer pointer;
+ typedef typename ChainRepresentation::const_pointer const_pointer;
typedef Item value_type;
- /// @}
-
- public:
- ChainWrapper();
- ChainWrapper(const ChainWrapper& c);
-
- /// \name Whole Chain operations
- /// @{
- /** Add c to *this assuming both Chains are sorted in increasing order according to cmp. */
+ /// @}
+
+ public:
+ ChainWrapper();
+ ChainWrapper(const ChainWrapper& c);
+
+ /// \name Whole Chain operations
+ /// @{
+ /** Add c to *this assuming both Chains are sorted in increasing order according to cmp. */
template<class ConsistencyComparison>
- Self& add(const Self& c, const ConsistencyComparison& cmp);
+ Self& add(const Self& c, const ConsistencyComparison& cmp);
- void swap(ChainWrapper& c); ///< Swaps the contents of c and *this (like STL's swap destroys c)
-
+ void swap(ChainWrapper& c); ///< Swaps the contents of c and *this (like STL's swap destroys c)
+
template<class ConsistencyComparison>
- void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
+ void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
size_t size() const { return SizeStorage<Container>::size(*this); }
- using ChainRepresentation::empty;
- using ChainRepresentation::clear;
+ using ChainRepresentation::empty;
+ using ChainRepresentation::clear;
/// @}
-
- /// \name Modifiers
- /// @{
- using ChainRepresentation::erase;
+
+ /// \name Modifiers
+ /// @{
+ using ChainRepresentation::erase;
template<class ConsistencyComparison>
- void append(const_reference x, const ConsistencyComparison& cmp);
- /// @}
-
- /// \name Accessors
- /// @{
- using ChainRepresentation::begin;
- using ChainRepresentation::end;
+ void append(const_reference x, const ConsistencyComparison& cmp);
+ /// @}
+
+ /// \name Accessors
+ /// @{
+ using ChainRepresentation::begin;
+ using ChainRepresentation::end;
template<class OrderComparison>
- const_reference top(const OrderComparison& cmp) const; ///< First element in cmp order
+ const_reference top(const OrderComparison& cmp) const; ///< First element in cmp order
- boost::optional<const_iterator> contains(const_reference x) const; ///< tests whether chain contains x
- boost::optional<iterator> contains(reference x); ///< tests whether chain contains x
- /// @}
-
- /// \name Debugging
- /// @{
+ boost::optional<const_iterator> contains(const_reference x) const; ///< tests whether chain contains x
+ boost::optional<iterator> contains(reference x); ///< tests whether chain contains x
+ /// @}
+
+ /// \name Debugging
+ /// @{
template<class OrderComparison>
- const_reference get_first(const OrderComparison& cmp) const; ///< First element in cmp order
+ const_reference get_first(const OrderComparison& cmp) const; ///< First element in cmp order
template<class OutputMap>
std::string tostring(const OutputMap& outmap = OutputMap()) const;
- /// @}
-
- private:
- using ChainRepresentation::push_back;
- using ChainRepresentation::insert;
-
- private:
- // Serialization
- typedef Container Parent;
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, boost::serialization::version_type );
+ /// @}
+
+ private:
+ using ChainRepresentation::push_back;
+ using ChainRepresentation::insert;
+
+ private:
+ // Serialization
+ typedef Container Parent;
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, boost::serialization::version_type );
};
#include "chain.hpp"
--- a/include/topology/chain.hpp Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/chain.hpp Sat Dec 27 14:33:25 2008 -0800
@@ -14,12 +14,12 @@
using boost::serialization::make_binary_object;
#ifdef LOGGING
-static rlog::RLogChannel* rlChain = DEF_CHANNEL( "topology/chain", rlog::Log_Debug);
+static rlog::RLogChannel* rlChain = DEF_CHANNEL( "topology/chain", rlog::Log_Debug);
#endif // LOGGING
#ifdef COUNTERS
-static Counter* cChainAddBasic = GetCounter("chain/add/basic");
-static Counter* cChainAddComparison = GetCounter("chain/add/comparison");
+static Counter* cChainAddBasic = GetCounter("chain/add/basic");
+static Counter* cChainAddComparison = GetCounter("chain/add/comparison");
#endif // COUNTERS
template<class C>
@@ -36,38 +36,28 @@
template<class ConsistencyCmp>
void
ChainWrapper<C>::
-append(const_reference x, const ConsistencyCmp& cmp)
+append(const_reference x, const ConsistencyCmp& cmp)
{
- SizeStorage<Container>::increment();
+ SizeStorage<Container>::operator++();
- // First try the special cases that x goes at the end
- const_reference last = ChainRepresentation::back();
- if (empty() || cmp(last, x))
- {
- push_back(x);
- return;
- }
+ // First try the special cases that x goes at the end
+ if (empty() || cmp(ChainRepresentation::back(), x))
+ {
+ push_back(x);
+ return;
+ }
- for (iterator cur = begin(); cur != end(); ++cur)
- if (cmp(x, *cur))
- {
- insert(cur, x);
- return;
- }
+ insert(std::upper_bound(begin(), end(), x, cmp), x);
}
-
+
template<class C>
template<class OrderComparison>
-typename ChainWrapper<C>::const_reference
+typename ChainWrapper<C>::const_reference
ChainWrapper<C>::
top(const OrderComparison& cmp) const
{
- AssertMsg(!empty(), "Chain must not be empty for low()");
- const_iterator min = begin();
- for (const_iterator cur = ++begin(); cur != end(); ++cur)
- if (cmp(*cur, *min))
- min = cur;
- return *min;
+ AssertMsg(!empty(), "Chain must not be empty for low()");
+ return *std::min_element(begin(), end(), cmp);
}
template<class C>
@@ -75,7 +65,7 @@
ChainWrapper<C>::
swap(ChainWrapper& c)
{
- ChainRepresentation::swap(c);
+ ChainRepresentation::swap(c);
SizeStorage<Container>::swap(c);
}
@@ -93,13 +83,6 @@
ChainWrapper<C>::
contains(const_reference x) const
{
-#if 0
- for (const_iterator cur = begin(); cur != end(); ++cur)
- if (!cmp(*cur, x) && !cmp(x, *cur))
- return make_optional(cur);
-
- return make_optional(false);
-#endif
const_iterator res = std::find(begin(), end(), x);
return make_optional(res != end(), res);
}
@@ -120,13 +103,13 @@
tostring(const OutputMap& outmap) const
{
std::string str;
- for (const_iterator cur = begin(); cur != end(); ++cur)
- {
+ for (const_iterator cur = begin(); cur != end(); ++cur)
+ {
if (cur != begin()) str += ", ";
- str += outmap(*cur);
- }
- // str += "(last: " + *last + ")"; // For debugging only
- return str;
+ str += outmap(*cur);
+ }
+ // str += "(last: " + *last + ")"; // For debugging only
+ return str;
}
template<class C>
@@ -139,69 +122,14 @@
// however, I believe it creates costly overhead for Containers that are lists.
// Need to put some thought into this.
ChainRepresentation tmp;
- SizeStorage<Container> size;
- // FIXME: need to do something about the output
- rLog(rlChain, "Adding chains"); //: ": %s + %s", tostring(*this).c_str(), tostring(c).c_str());
-
- iterator cur1 = begin();
- const_iterator cur2 = c.begin();
-
- while (cur2 != c.end())
- {
- if (cur1 == end())
- {
- while (cur2 != c.end())
- {
- tmp.push_back(*cur2++);
- size.increment();
- Count(cChainAddBasic);
- }
- rLog(rlChain, "After addition"); //: %s", tostring(*this).c_str());
-
- static_cast<ChainRepresentation*>(this)->swap(tmp);
- static_cast<SizeStorage<Container>*>(this)->swap(size);
- return *this;
- }
-
- // mod 2
- int res = cmp.compare(*cur1, *cur2);
- Count(cChainAddComparison);
- rLog(rlChain, "Comparison result: %i", res);
- if (res == 0) // *cur1 == *cur2
- {
- rLog(rlChain, "Equality");
- //cur1 = erase(cur1); // erase cur1 --- as a result cur1 will be pointing at old_cur1++
- ++cur1;
- ++cur2;
- } else if (res < 0) // *cur1 < *cur2
- {
- rLog(rlChain, "Less than");
- //cur1++;
- tmp.push_back(*cur1++);
- size.increment();
- } else if (res > 0) // *cur1 > *cur2
- {
- rLog(rlChain, "Greater than");
- //insert(cur1, *cur2);
- //++cur2;
- tmp.push_back(*cur2++);
- size.increment();
- }
- Count(cChainAddBasic);
- }
- while (cur1 != end())
- {
- tmp.push_back(*cur1++);
- size.increment();
- Count(cChainAddBasic);
- }
-
- rLog(rlChain, "After addition"); //: %s", tostring(*this).c_str());
+ CountingBackInserter<ChainRepresentation> bi(tmp);
+ std::set_symmetric_difference(begin(), end(), c.begin(), c.end(), bi, cmp);
static_cast<ChainRepresentation*>(this)->swap(tmp);
- static_cast<SizeStorage<Container>*>(this)->swap(size);
- return *this;
+ static_cast<SizeStorage<Container>*>(this)->swap(bi);
+
+ return *this;
}
template<class C>
@@ -211,12 +139,12 @@
get_first(const OrderComparison& cmp) const
{ return top(cmp); }
-
+
template<class C>
template<class Archive>
-void
+void
ChainWrapper<C>::
serialize(Archive& ar, boost::serialization::version_type )
{
- ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
+ ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
}
--- a/include/topology/conesimplex.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/conesimplex.h Sat Dec 27 14:33:25 2008 -0800
@@ -10,38 +10,53 @@
#include <iostream>
#include "utilities/types.h"
+#include "simplex.h"
-template<class S>
-class ConeSimplex: public S
+template<class V, class T = Empty>
+class ConeSimplex: public Simplex<V,T>
{
- public:
- typedef S Parent;
- typedef ConeSimplex<S> Self;
- typedef std::list<Self> Cycle;
+ public:
+ typedef Simplex<V,T> Parent;
+ typedef ConeSimplex<V,T> Self;
public:
- ConeSimplex(const Self& s):
+ ConeSimplex(const Self& s):
Parent(s), coned_(s.coned_) {}
- ConeSimplex(const Parent& parent,
- bool coned = false):
- Parent(parent), coned_(coned) {}
-
- Cycle boundary() const;
- bool coned() const { return coned_; }
+ ConeSimplex(const Parent& parent,
+ bool coned = false):
+ Parent(parent), coned_(coned) {}
+
+ Cycle boundary() const;
+ bool coned() const { return coned_; }
Dimension dimension() const { return coned_ ? (Parent::dimension() + 1) : Parent::dimension(); }
- bool operator<(const Self& other) const { if (coned_ ^ other.coned_) return !coned_; else return Parent::operator<(other); }
bool operator==(const Self& other) const { return !(coned_ ^ other.coned_) && Parent::operator==(other); }
- std::ostream& operator<<(std::ostream& out) const;
-
- private:
- bool coned_;
+ std::ostream& operator<<(std::ostream& out) const;
+
+ struct ConedVertexComparison;
+
+ private:
+ bool coned_;
};
-template<class S>
-std::ostream& operator<<(std::ostream& out, const ConeSimplex<S>& s) { return s.operator<<(out); }
+template<class V, class T>
+struct ConeSimplex<V,T>::ConedVertexComparison: public typename Simplex<V,T>::VertexComparison
+{
+ typedef typename Simplex<V,T>::VertexComparison Parent;
+
+ bool operator()(const Self& a, const Self& b) const
+ {
+ if (a.coned() ^ b.coned())
+ return b.coned(); // coned simplices shall come after non-coned ones
+ else
+ return Parent::operator()(a,b); // within coned/non-coned order by vertices
+ }
+};
+
+template<class V, class T>
+std::ostream& operator<<(std::ostream& out, const ConeSimplex<V,T>& s) { return s.operator<<(out); }
#include "conesimplex.hpp"
--- a/include/topology/filtrationcontainer.h Thu Dec 25 14:24:02 2008 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2006
- */
-
-#ifndef __FILTRATIONCONTAINER_H__
-#define __FILTRATIONCONTAINER_H__
-
-#include "utilities/consistencylist.h"
-#include "cycle.h"
-
-/**
- * FiltrationContainer class. Serves as a parent of Filtration that
- * describes the container functionality. Used by FiltrationSimplex
- * to get Cycle representation.
- *
- * \ingroup topology
- */
-template<class FltrSmplx>
-class FiltrationContainer: public ConsistencyList<FltrSmplx>
-{
- public:
- typedef FltrSmplx FiltrationSimplex;
- typedef ConsistencyList<FiltrationSimplex> ConsistencyLst;
-
- /// \name Cycles and Trails
- /// @{
- /// Index is and therfore acts like an iterator. The name is preserved for historical reasons.
- typedef typename ConsistencyLst::iterator Index;
- /// const_Index is a const_iterator
- typedef typename ConsistencyLst::const_iterator const_Index;
- /// @}
-
- /// \name Cycles and Trails
- /// @{
- typedef typename ConsistencyLst::GreaterThanComparison CyclesComparator;
- typedef typename ConsistencyLst::LessThanComparison TrailsComparator;
- typedef typename ConsistencyLst::ConsistencyComparison ConsistencyComparator;
- typedef ::Cycle<Index, CyclesComparator, ConsistencyComparator> Cycle;
- typedef ::Cycle<Index, TrailsComparator, ConsistencyComparator> Trail;
- /// @}
-
- template<class U>
- struct rebind { typedef FiltrationContainer<U> other; };
-};
-
-#endif // __FILTRATIONCONTAINER_H__
--- a/include/topology/filtrationsimplex.h Thu Dec 25 14:24:02 2008 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,92 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2006
- */
-
-#ifndef __FILTRATIONSIMPLEX_H__
-#define __FILTRATIONSIMPLEX_H__
-
-#include "filtrationcontainer.h"
-#include "vineyard.h"
-#include "utilities/types.h"
-
-#include <list>
-
-#if 0
-#include <boost/serialization/access.hpp>
-#include <boost/serialization/nvp.hpp>
-#include <boost/serialization/list.hpp>
-#endif
-
-/**
- * Evaluator is a base class for the structures that are able to return a value
- * given a simplex.
- *
- * \ingroup topology
- */
-template<class Smplx>
-class Evaluator
-{
- public:
- typedef Smplx Simplex;
-
- virtual RealType time() const { return 0; }
- virtual RealType value(const Simplex& s) const { return 0; }
-
- virtual ~Evaluator() {}
-};
-
-/**
- * FiltrationSimplex stores information needed for the RU-decomposition:
- * cycle (column of R), trail (row of U), and pair.
- *
- * \ingroup topology
- */
-template<class Smplx>
-class FiltrationSimplex: public Smplx
-{
- public:
- typedef Smplx Simplex;
- typedef FiltrationSimplex<Simplex> Self;
- typedef FiltrationContainer<Self> Container;
- typedef Simplex Parent;
-
- typedef Vine<Simplex> VineS;
- typedef typename Container::Cycle Cycle;
- typedef typename Container::Trail Trail;
- typedef typename Container::Index Index;
-
- typedef Evaluator<Simplex> EvaluatorS;
-
- FiltrationSimplex(const Simplex& s):
- Simplex(s), vine_(0) {}
-
-
- /// \name Core functionality
- /// @{
- void set_pair(Index pair) { pair_ = pair; }
- bool sign() const { return cycles_column.empty(); }
- bool is_paired() const { return pair() != pair()->pair(); }
- void set_vine(VineS* v) { vine_ = v; }
- using Parent::dimension;
- /// @}
-
-
- /// \name Accessors
- /// @{
- Cycle& cycle() { return cycles_column; }
- Trail& trail() { return trails_row; }
- const Cycle& cycle() const { return cycles_column; }
- const Trail& trail() const { return trails_row; }
- Index pair() const { return pair_; }
- VineS* vine() const { return vine_; }
- /// @}
-
- private:
- Cycle cycles_column;
- Trail trails_row;
- Index pair_;
- VineS* vine_;
-};
-
-#endif // __FILTRATIONSIMPLEX_H__
--- a/include/topology/lowerstarfiltration.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/lowerstarfiltration.h Sat Dec 27 14:33:25 2008 -0800
@@ -6,138 +6,55 @@
#ifndef __LOWERSTARFILTRATION_H__
#define __LOWERSTARFILTRATION_H__
-#include "utilities/log.h"
-
-#include "filtration.h"
-#include "simplex.h"
-#include "utilities/consistencylist.h"
-#include <boost/utility.hpp>
-#include <list>
-#include "utilities/types.h"
-
#include <boost/serialization/access.hpp>
#include <boost/serialization/vector.hpp>
-#include <boost/serialization/map.hpp>
-#include <boost/serialization/base_object.hpp>
#include <boost/serialization/nvp.hpp>
+/**
+ * Struct: MaxVertexComparison
+ *
+ * Functor that determines which simplex has a higher vertex with respect to VertexComparison_
+ */
+template<class Simplex_, class VertexComparison_>
+struct MaxVertexComparison
+{
+ typedef VertexComparison_ VertexComparison;
+ typedef Simplex_ Simplex;
+
+ MaxVertexComparison(const VertexComparison& vcmp):
+ vcmp_(vcmp) {}
+
+ bool operator()(const Simplex& s1, const Simplex& s2) const
+ {
+ return std::max_element(s1.vertices().begin(), s1.vertices().end(), vcmp) <
+ std::max_element(s2.vertices().begin(), s2.vertices().end(), vcmp);
+ }
+
+ VertexComparison vcmp_;
+};
+
-template<class VI,
- class Smplx = SimplexWithAttachment<VI>,
- class FltrSmplx = FiltrationSimplex<Smplx>,
- class Vnrd = Vineyard<FltrSmplx> >
-class LowerStarFiltration: public Filtration<Smplx, FltrSmplx, Vnrd>
+/**
+ * Map from i-th vertex to its index in the filtration.
+ */
+template<class Index_>
+class VertexSimplexMap
{
- public:
- // Treat VertexIndex as an iterator
- typedef VI VertexIndex;
- typedef Smplx Simplex;
- typedef Filtration<Simplex> Parent;
- typedef typename Parent::Vineyard Vineyard;
-
- typedef typename Parent::Index Index;
- typedef typename Parent::const_Index const_Index;
- typedef typename Parent::Cycle Cycle;
- typedef typename Parent::Trail Trail;
- typedef typename Simplex::Cycle SimplexBoundaryCycle;
-
- template<class IndexType>
- class VertexType;
-
- struct VertexDescriptor;
- typedef ConsistencyList<VertexDescriptor> VertexOrder;
- typedef typename VertexOrder::iterator VertexOrderIndex;
- typedef typename VertexOrder::const_iterator const_VertexOrderIndex;
- typedef typename VertexOrder::LessThanComparison VertexOrderComparison;
- struct SimplexAttachmentComparison;
-
- public:
- template<class VertexCmp>
- LowerStarFiltration(VertexIndex begin, VertexIndex end, const VertexCmp& cmp, Vineyard* vineyard);
+ public:
+ typedef Index_ Index;
+ typedef std::vector<FiltrationIndex> VertexVector;
+
+ VertexSimplexMap(Index begin, Index end, const Map& m)
+ {
+ for (FiltrationIndex cur = begin; cur != end; ++cur)
+ if (m[cur].dimension() == 0)
+ vertices_.push_back(cur);
+ }
- using Parent::size;
- using Parent::begin;
- using Parent::end;
- VertexIndex num_vertices() const { return vertex_order.size(); }
- const VertexOrderComparison&
- get_vertex_cmp() const { return vertex_cmp; }
-
- Index append(const Simplex& s);
- bool transpose_vertices(const VertexOrderIndex& voi);
-
- protected:
- /// Hint function: if unsure, should return true
- virtual bool neighbors(VertexIndex v1, VertexIndex v2) const { return true; }
-
- private:
- bool transpose(Index i);
- void assert_pairing(Index i);
-
- private:
- VertexOrder vertex_order;
- VertexOrderComparison vertex_cmp;
-
- /* Serialization */
- protected:
- LowerStarFiltration() {}
-
- private:
- friend class boost::serialization::access;
-
- template<class Archive>
- void save(Archive& ar, version_type ) const { ar << BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent); }
-
- template<class Archive>
- void load(Archive& ar, version_type );
-
- BOOST_SERIALIZATION_SPLIT_MEMBER()
+ private:
+ VertexVector vertices_;
};
-/**
- * Helper class that describes lower star vertex type. Defines essential
- * methods that LowerStarFiltration expects from its vertex type. Actual
- * vertex types should inherit from it.
- */
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-template<class IndexType_>
-class LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-VertexType
-{
- public:
- typedef IndexType_ IndexType;
-
- VertexType(IndexType ii = 0): i_(ii) {}
-
- IndexType index() const { return i_; }
- void set_index(IndexType i) { i_ = i; }
- VertexOrderIndex get_order() const { return order_; }
- void set_order(const VertexOrderIndex& o) { order_ = o; }
-
- private:
- IndexType i_;
- VertexOrderIndex order_;
-};
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-struct LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-VertexDescriptor
-{
- VertexDescriptor(VertexIndex vi, Index si):
- vertex_index(vi), simplex_index(si)
- {}
-
- VertexIndex vertex_index;
- Index simplex_index;
-};
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-struct LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-SimplexAttachmentComparison
-{
- bool operator()(const Simplex& first, const Simplex& second) const;
- VertexOrderComparison vertex_cmp;
-};
-
-#include "lowerstarfiltration.hpp"
+// TODO: transpose_vertices(Index, Filtration, Persistence, Visitor);
#endif // __LOWERSTARFILTRATION_H__
--- a/include/topology/rips.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/rips.h Sat Dec 27 14:33:25 2008 -0800
@@ -4,18 +4,19 @@
#include <vector>
#include <string>
#include "simplex.h"
-#include "filtrationsimplex.h"
/**
- * Rips complex class
+ * RipsBase class
+ *
+ * Base class for the generator of Rips complexes.
*
* Distances_ is expected to define types IndexType and DistanceType as well as
* provide operator()(...) which given two IndexTypes should return
* the distance between them. There should be methods begin() and end()
* for iterating over IndexTypes as well as a method size().
*/
-template<class Distances_, class Simplex_ = SimplexWithVertices<typename Distances_::IndexType> >
-class Rips
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
+class RipsBase
{
public:
typedef Distances_ Distances;
@@ -26,30 +27,69 @@
typedef std::vector<Simplex> SimplexVector;
class Evaluator;
+ class Comparison;
+ struct ComparePair;
public:
- Rips(const Distances& distances):
+ RipsBase(const Distances& distances):
distances_(distances) {}
- void generate(Dimension k, DistanceType max); /// generate k-skeleton of the Rips complex
- const SimplexVector&
- simplices() const { return simplices_; }
- size_t size() const { return simplices_.size(); }
-
const Distances& distances() const { return distances_; }
DistanceType max_distance() const;
-
- void print() const;
+
+ DistanceType distance(const Simplex& s1, const Simplex& s2) const;
private:
- struct ComparePair;
+ const Distances& distances_;
+};
+
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
+class RipsGenerator: public RipsBase<Distances_, Simplex_>
+{
+ public:
+ typedef RipsBase<Distances_, Simplex_> Parent;
+ typedef typename Parent::Distances Distances;
+ typedef typename Parent::Simplex Simplex;
+ typedef typename Parent::SimplexVector SimplexVector;
+ typedef typename Parent::DistanceType DistanceType;
+ typedef typename Parent::IndexType IndexType;
+ typedef typename Parent::ComparePair ComparePair;
- const Distances& distances_;
- SimplexVector simplices_;
+ RipsGenerator(const Distances& distances):
+ Parent(distances) {}
+
+ using Parent::distances;
+
+ /// generate k-skeleton of the Rips complex
+ void generate(SimplexVector& v, Dimension k, DistanceType max) const;
};
+// Much more memory efficient, but also much slower
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
+class RipsGeneratorMemory: public RipsBase<Distances_, Simplex_>
+{
+ public:
+ typedef RipsBase<Distances_, Simplex_> Parent;
+ typedef typename Parent::Distances Distances;
+ typedef typename Parent::Simplex Simplex;
+ typedef typename Parent::SimplexVector SimplexVector;
+ typedef typename Parent::DistanceType DistanceType;
+ typedef typename Parent::IndexType IndexType;
+ typedef typename Parent::ComparePair ComparePair;
+
+ RipsGeneratorMemory(const Distances& distances):
+ Parent(distances) {}
+
+ using Parent::distances;
+ using Parent::distance;
+
+ /// generate k-skeleton of the Rips complex
+ void generate(SimplexVector& v, Dimension k, DistanceType max) const;
+};
+
+
template<class Distances_, class Simplex_>
-class Rips<Distances_, Simplex_>::Evaluator: public ::Evaluator<Simplex_>
+class RipsBase<Distances_, Simplex_>::Evaluator
{
public:
typedef Simplex_ Simplex;
@@ -57,13 +97,27 @@
Evaluator(const Distances& distances):
distances_(distances) {}
- virtual DistanceType
- value(const Simplex& s) const;
+ DistanceType operator()(const Simplex& s) const;
private:
const Distances& distances_;
};
+template<class Distances_, class Simplex_>
+class RipsBase<Distances_, Simplex_>::Comparison
+{
+ public:
+ typedef Simplex_ Simplex;
+
+ Comparison(const Distances& distances):
+ eval_(distances) {}
+
+ bool operator()(const Simplex& s1, const Simplex& s2) const { return eval_(s1) < eval_(s2); }
+
+ private:
+ Evaluator eval_;
+};
+
/**
* ExplicitDistances stores the pairwise distances of Distances_ instance passed at construction.
* It's a protypical Distances template argument for the Rips complex.
--- a/include/topology/rips.hpp Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/rips.hpp Sat Dec 27 14:33:25 2008 -0800
@@ -3,9 +3,11 @@
#include <boost/utility.hpp>
#include <iostream>
#include <utilities/log.h>
+#include <utilities/counter.h>
#ifdef LOGGING
static rlog::RLogChannel* rlRips = DEF_CHANNEL("rips/info", rlog::Log_Debug);
+static rlog::RLogChannel* rlRipsDebug = DEF_CHANNEL("rips/debug", rlog::Log_Debug);
#endif // LOGGING
#ifdef COUNTERS
@@ -13,7 +15,7 @@
#endif // COUNTERS
template<class Distances_, class Simplex_>
-struct Rips<Distances_, Simplex_>::ComparePair
+struct RipsBase<Distances_, Simplex_>::ComparePair
{
ComparePair(const Distances& distances):
distances_(distances) {}
@@ -27,51 +29,52 @@
template<class DistanceType_, class Simplex_>
void
-Rips<DistanceType_, Simplex_>::
-generate(Dimension k, DistanceType max)
+RipsGenerator<DistanceType_, Simplex_>::
+generate(SimplexVector& simplices, Dimension k, DistanceType max) const
{
// Order all the edges
typedef std::vector< std::pair<IndexType, IndexType> > EdgeVector;
EdgeVector edges;
- for (IndexType a = distances_.begin(); a != distances_.end(); ++a)
+ for (IndexType a = distances().begin(); a != distances().end(); ++a)
{
Simplex ssx; ssx.add(a);
- simplices_.push_back(ssx);
- for (IndexType b = boost::next(a); b != distances_.end(); ++b)
+ simplices.push_back(ssx);
+ for (IndexType b = boost::next(a); b != distances().end(); ++b)
{
- if (distances_(a,b) <= max)
+ if (distances()(a,b) <= max)
edges.push_back(std::make_pair(a,b));
}
}
- std::sort(edges.begin(), edges.end(), ComparePair(distances_));
+ std::sort(edges.begin(), edges.end(), ComparePair(distances()));
// Generate simplices
- std::vector<std::vector<size_t> > vertex_star(distances_.size());
+ std::vector<std::vector<size_t> > vertex_star(distances().size());
for(typename EdgeVector::const_iterator cur = edges.begin(); cur != edges.end(); ++cur)
{
- rLog(rlRips, "Current edge: %d %d", cur->first, cur->second);
+ rLog(rlRipsDebug, "Current edge: %d %d", cur->first, cur->second);
// Create the edge
Simplex edge; edge.add(cur->first); edge.add(cur->second);
- simplices_.push_back(edge);
+ simplices.push_back(edge);
if (k <= 1) continue;
- vertex_star[cur->first].push_back(simplices_.size() - 1);
- vertex_star[cur->second].push_back(simplices_.size() - 1);
+ vertex_star[cur->first].push_back(simplices.size() - 1);
+ vertex_star[cur->second].push_back(simplices.size() - 1);
// Go through a star
size_t sz = vertex_star[cur->first].size() - 1;
for (size_t i = 0; i < sz; ++i)
{
- const Simplex& ssx = simplices_[vertex_star[cur->first][i]];
- rLog(rlRips, " %s", tostring(ssx).c_str());
+ const Simplex& ssx = simplices[vertex_star[cur->first][i]];
+ // FIXME: eventually can uncomment, missing Empty::operator<<()
+ // rLog(rlRipsDebug, " %s", tostring(ssx).c_str());
bool accept = true;
for (typename Simplex::VertexContainer::const_iterator v = ssx.vertices().begin(); v != ssx.vertices().end(); ++v)
{
if (*v == cur->first) continue;
- if ( distances_(*v, cur->second) > distances_(cur->first, cur->second) ||
- ((distances_(*v, cur->second) == distances_(cur->first, cur->second)) &&
+ if ( distances()(*v, cur->second) > distances()(cur->first, cur->second) ||
+ ((distances()(*v, cur->second) == distances()(cur->first, cur->second)) &&
(*v > cur->first)))
{
accept = false;
@@ -81,30 +84,71 @@
if (accept)
{
Simplex tsx(ssx); tsx.add(cur->second);
- simplices_.push_back(tsx);
- rLog(rlRips, " Accepting: %s", tostring(tsx).c_str());
+ simplices.push_back(tsx);
+ // rLog(rlRipsDebug, " Accepting: %s", tostring(tsx).c_str());
// Update stars
if (tsx.dimension() < k - 1)
- for (typename Simplex::VertexContainer::const_iterator v = tsx.vertices().begin(); v != tsx.vertices().end(); ++v)
- vertex_star[*v].push_back(simplices_.size() - 1);
+ for (typename Simplex::VertexContainer::const_iterator v = static_cast<const Simplex&>(tsx).vertices().begin();
+ v != static_cast<const Simplex&>(tsx).vertices().end();
+ ++v)
+ vertex_star[*v].push_back(simplices.size() - 1);
}
}
}
}
-template<class Distances_, class Simplex_>
+template<class DistanceType_, class Simplex_>
void
-Rips<Distances_, Simplex_>::
-print() const
+RipsGeneratorMemory<DistanceType_, Simplex_>::
+generate(SimplexVector& simplices, Dimension k, DistanceType max) const
{
- for (typename SimplexVector::const_iterator cur = simplices_.begin(); cur != simplices_.end(); ++cur)
- std::cout << *cur << std::endl;
+ for (IndexType v = distances().begin(); v != distances().end(); ++v)
+ {
+ simplices.push_back(Simplex());
+ simplices.back().add(v);
+ }
+ size_t last_vertex = simplices.size() - 1;
+ size_t begin_previous_dimension = 0;
+ size_t end_previous_dimension = simplices.size() - 1;
+ typename Simplex::VertexComparison vcmp;
+
+ for (Dimension d = 1; d < k; ++d)
+ {
+ //rLog(rlRips, "Generating dimension %d", d);
+ //rLog(rlRips, " Begin previous dimension: %d", begin_previous_dimension);
+ //rLog(rlRips, " End previous dimension: %d", end_previous_dimension);
+ for (size_t i = 0; i <= last_vertex; ++i)
+ {
+ for (size_t j = begin_previous_dimension; j <= end_previous_dimension; ++j)
+ if (!simplices[j].contains(simplices[i]) &&
+ vcmp(simplices[i], simplices[j]) &&
+ distance(simplices[i], simplices[j]) <= max)
+ {
+ simplices.push_back(Simplex(simplices[j]));
+ simplices.back().join(simplices[i]);
+ }
+ }
+ begin_previous_dimension = end_previous_dimension + 1;
+ end_previous_dimension = simplices.size() - 1;
+ }
}
template<class Distances_, class Simplex_>
-typename Rips<Distances_, Simplex_>::DistanceType
-Rips<Distances_, Simplex_>::
+typename RipsBase<Distances_, Simplex_>::DistanceType
+RipsBase<Distances_, Simplex_>::
+distance(const Simplex& s1, const Simplex& s2) const
+{
+ DistanceType mx = 0;
+ for (typename Simplex::VertexContainer::const_iterator a = s1.vertices().begin(); a != s1.vertices().end(); ++a)
+ for (typename Simplex::VertexContainer::const_iterator b = s2.vertices().begin(); b != s2.vertices().end(); ++b)
+ mx = std::max(mx, distances_(*a,*b));
+ return mx;
+}
+
+template<class Distances_, class Simplex_>
+typename RipsBase<Distances_, Simplex_>::DistanceType
+RipsBase<Distances_, Simplex_>::
max_distance() const
{
DistanceType mx = 0;
@@ -115,9 +159,9 @@
}
template<class Distances_, class Simplex_>
-typename Rips<Distances_, Simplex_>::DistanceType
-Rips<Distances_, Simplex_>::Evaluator::
-value(const Simplex& s) const
+typename RipsBase<Distances_, Simplex_>::DistanceType
+RipsBase<Distances_, Simplex_>::Evaluator::
+operator()(const Simplex& s) const
{
DistanceType mx = 0;
for (typename Simplex::VertexContainer::const_iterator a = s.vertices().begin(); a != s.vertices().end(); ++a)
--- a/include/topology/simplex.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/simplex.h Sat Dec 27 14:33:25 2008 -0800
@@ -32,7 +32,7 @@
template<class V, class T = Empty>
class Simplex
{
- public:
+ public:
/* Typedefs: Public Types
*
* Vertex - vertex type (template parameter V)
@@ -40,74 +40,70 @@
* Self -
* Boundary - type in which the boundary is stored
*/
- typedef V Vertex;
+ typedef V Vertex;
typedef T Data;
- typedef Simplex<Vertex, Data> Self;
- typedef std::list<Self> Boundary;
-
+ typedef Simplex<Vertex, Data> Self;
+ typedef std::list<Self> Boundary;
+
/* Typedefs: Internal representation
*
* VertexContainer - internal representation of the vertices
* VerticesDataPair - `compressed_pair` of VertexContainer and Data
*/
- typedef std::vector<Vertex> VertexContainer;
+ typedef std::vector<Vertex> VertexContainer;
typedef boost::compressed_pair<VertexContainer, Data> VerticesDataPair;
-
- // TODO
-
- /// \name Constructors
- /// @{
+
+ /// \name Constructors
+ /// @{
//
/// Constructor: Simplex()
/// Default constructor
- Simplex() {}
+ Simplex() {}
/// Constructor: Simplex(Self other)
/// Copy constructor
- Simplex(const Self& other):
- vdpair_(other.vertices(), other.data()) {}
+ Simplex(const Self& other):
+ vdpair_(other.vertices(), other.data()) {}
+ /// Constructor: Simplex(Data d)
+ /// Initialize simplex with data
+ Simplex(const Data& d) { data() = d; }
/// Constructor: Simplex(Iterator bg, Iterator end)
/// Initialize simplex by copying vertices in range [bg, end)
- template<class Iterator>
- Simplex(Iterator bg, Iterator end, const Data& d = Data()) { join(bg, end); data() = d; }
+ template<class Iterator>
+ Simplex(Iterator bg, Iterator end, const Data& d = Data()) { join(bg, end); data() = d; }
/// Constructor: Simplex(VertexContainer v)
/// Initialize simplex by copying the given VertexContainer
- Simplex(const VertexContainer& v, const Data& d = Data()):
- vdpair_(v, d) { std::sort(vertices().begin(), vertices().end()); }
- /// Constructor: Simplex(Dimension d, Vertex v)
- /// Initialize simplex of dimension d, and set its first vertex to v
- Simplex(Dimension d, Vertex v) { vertices().reserve(d+1); add(v); }
- /// Constructor: Simplex(Dimension d)
- /// Initialize simplex of dimension d
- Simplex(Dimension d) { vertices().resize(d+1); }
- /// @}
-
- /// \name Core
- /// @{
+ Simplex(const VertexContainer& v, const Data& d = Data()):
+ vdpair_(v, d) { std::sort(vertices().begin(), vertices().end()); }
+ /// @}
+
+ /// \name Core
+ /// @{
///
/// Function: boundary()
/// Returns the boundary of the simplex (of type Boundary)
- Boundary boundary() const;
+ Boundary boundary() const;
/// Function: dimension()
/// Returns the dimension of the simplex
- Dimension dimension() const { return vertices().size() - 1; }
- /// @}
-
- const Data& data() const { return vdpair_.second(); }
+ Dimension dimension() const { return vertices().size() - 1; }
+ /// @}
+
+ const Data& data() const { return vdpair_.second(); }
Data& data() { return vdpair_.second(); }
- const VertexContainer& vertices() const { return vdpair_.first(); }
-
- /// \name Vertex manipulation
- /// @{
- bool contains(const Vertex& v) const;
- void add(const Vertex& v);
+ const VertexContainer& vertices() const { return vdpair_.first(); }
+
+ /// \name Vertex manipulation
+ /// @{
+ bool contains(const Vertex& v) const;
+ bool contains(const Self& s) const;
+ void add(const Vertex& v);
template<class Iterator>
void join(Iterator bg, Iterator end);
- void join(const Self& other) { join(other.vertices.begin(), other.vertices().end()); }
- /// @}
+ void join(const Self& other) { join(other.vertices().begin(), other.vertices().end()); }
+ /// @}
- const Self& operator=(const Self& s) { vdpair_ = s.vdpair_; return *this; }
+ const Self& operator=(const Self& s) { vdpair_ = s.vdpair_; return *this; }
- std::ostream& operator<<(std::ostream& out) const;
+ std::ostream& operator<<(std::ostream& out) const;
/* Classes: Comparisons
*
@@ -126,18 +122,18 @@
struct DataEvaluator;
struct DimensionExtractor;
-
- private:
- VertexContainer& vertices() { return vdpair_.first(); }
+
+ private:
+ VertexContainer& vertices() { return vdpair_.first(); }
VerticesDataPair vdpair_;
- private:
- /* Serialization */
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ private:
+ /* Serialization */
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
@@ -161,20 +157,21 @@
bool operator()(const Self& a, const Self& b) const { return a.data() < b.data(); }
};
-template<class V, class T>
-struct Simplex<V,T>::DataDimensionComparison
+template<class S>
+struct DataDimensionComparison
{
- typedef Self first_argument_type;
- typedef Self second_argument_type;
+ typedef S Simplex;
+ typedef Simplex first_argument_type;
+ typedef Simplex second_argument_type;
typedef bool result_type;
- bool operator()(const Self& a, const Self& b) const
- {
- if (a.dimension() == b.dimension())
- return a.data() < b.data();
- else
- return a.dimension() < b.dimension();
- }
+ bool operator()(const Simplex& a, const Simplex& b) const
+ {
+ if (a.dimension() == b.dimension())
+ return a.data() < b.data();
+ else
+ return a.dimension() < b.dimension();
+ }
};
template<class V, class T>
--- a/include/topology/simplex.hpp Thu Dec 25 14:24:02 2008 -0800
+++ b/include/topology/simplex.hpp Sat Dec 27 14:33:25 2008 -0800
@@ -9,19 +9,19 @@
/* Simplex */
template<class V, class T>
-typename Simplex<V,T>::Boundary
+typename Simplex<V,T>::Boundary
Simplex<V,T>::
boundary() const
{
typedef std::not_equal_to<Vertex> NotEqualVertex;
- Boundary bdry;
+ Boundary bdry;
if (dimension() == 0) return bdry;
-
+
for (typename VertexContainer::const_iterator cur = vertices().begin(); cur != vertices().end(); ++cur)
bdry.push_back(Self(boost::make_filter_iterator(std::bind2nd(NotEqualVertex(), *cur), vertices().begin(), vertices().end()),
boost::make_filter_iterator(std::bind2nd(NotEqualVertex(), *cur), vertices().end(), vertices().end())));
-
+
return bdry;
}
@@ -31,10 +31,19 @@
contains(const Vertex& v) const
{
// TODO: would std::find() be faster? (since most simplices we deal with are low dimensional)
- typename VertexContainer::const_iterator location = std::lower_bound(vertices().begin(), vertices().end(), v);
- return ((location != vertices().end()) && (*location == v));
+ typename VertexContainer::const_iterator location = std::lower_bound(vertices().begin(), vertices().end(), v);
+ return ((location != vertices().end()) && (*location == v));
}
-
+
+template<class V, class T>
+bool
+Simplex<V,T>::
+contains(const Self& s) const
+{
+ return std::includes( vertices().begin(), vertices().end(),
+ s.vertices().begin(), s.vertices().end());
+}
+
template<class V, class T>
void
Simplex<V,T>::
@@ -43,7 +52,7 @@
// TODO: would find() or lower_bound() followed by insert be faster?
vertices().push_back(v); std::sort(vertices().begin(), vertices().end());
}
-
+
template<class V, class T>
template<class Iterator>
void
@@ -55,27 +64,31 @@
}
template<class V, class T>
-std::ostream&
+std::ostream&
Simplex<V,T>::
operator<<(std::ostream& out) const
{
- for (typename VertexContainer::const_iterator cur = vertices().begin(); cur != vertices().end(); ++cur)
- out << *cur;
- out << " [" << data() << "] ";
+ typename VertexContainer::const_iterator cur = vertices().begin();
+ out << *cur;
+ for (++cur; cur != vertices().end(); ++cur)
+ {
+ out << ", " << *cur;
+ }
+ out << " [" << data() << "] ";
- return out;
+ return out;
}
-
+
template<class V, class T>
template<class Archive>
void
Simplex<V,T>::
-serialize(Archive& ar, version_type )
+serialize(Archive& ar, version_type )
{
ar & make_nvp("vertices", vertices());
ar & make_nvp("data", data());
}
template<class V, class T>
-std::ostream& operator<<(std::ostream& out, const Simplex<V,T>& s)
+std::ostream& operator<<(std::ostream& out, const Simplex<V,T>& s)
{ return s.operator<<(out); }
--- a/include/utilities/containers.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/utilities/containers.h Sat Dec 27 14:33:25 2008 -0800
@@ -30,11 +30,16 @@
{
public:
typedef Container_ Container;
+ typedef SizeStorage<Container> Self;
- SizeStorage(): size_(0) {}
+ SizeStorage(size_t size = 0): size_(size) {}
- void increment(size_t inc = 1) { size_ += inc; }
- void decrement(size_t dec = 1) { size_ -= dec; }
+ Self& operator+=(size_t inc) { size_ += inc; return *this; }
+ Self& operator-=(size_t dec) { size_ -= dec; return *this; }
+ Self& operator++() { ++size_; return *this; }
+ Self operator++(int) { Self tmp = *this; size_++; return tmp; }
+ Self& operator--() { --size_; return *this; }
+ Self operator--(int) { Self tmp = *this; size_--; return tmp; }
size_t size(const Container& c) const { return size_; }
void swap(SizeStorage& other) { std::swap(size_, other.size_); }
@@ -42,6 +47,26 @@
size_t size_;
};
+/**
+ * Class: CountingBackInserter<C>
+ *
+ * Derives from std::back_insert_iterator<C> and SizeStorage<C>,
+ * and keeps track of the number of inserted elements.
+ */
+template<class C>
+struct CountingBackInserter: public std::back_insert_iterator<C>,
+ public SizeStorage<C>
+{
+ typedef CountingBackInserter Self;
+ typedef std::back_insert_iterator<C> ParentIterator;
+ typedef SizeStorage<C> ParentSize;
+
+ CountingBackInserter(C& c):
+ ParentIterator(c) {}
+
+ Self& operator++() { ParentSize::operator++(); ParentIterator::operator++(); return *this; }
+ Self operator++(int i) { Self tmp = *this; ParentSize::operator++(i); ParentIterator::operator++(i); return tmp; }
+};
/* Specializations */
@@ -66,9 +91,9 @@
static void reserve(Container& c, size_t sz) { }
static void sort(Container& c, const Comparison& cmp = Comparison())
{
- std::vector<Item> tmp(c.begin(), c.end());
- std::sort(tmp.begin(), tmp.end(), cmp);
- std::copy(tmp.begin(), tmp.end(), c.begin());
+ std::vector<Item> tmp(c.begin(), c.end());
+ std::sort(tmp.begin(), tmp.end(), cmp);
+ std::copy(tmp.begin(), tmp.end(), c.begin());
}
};
@@ -79,9 +104,16 @@
{
public:
typedef std::vector<T> Container;
+ typedef SizeStorage<Container> Self;
+
+ SizeStorage(size_t size = 0) {}
- void increment(size_t inc = 1) {}
- void decrement(size_t dec = 1) {}
+ Self& operator+=(size_t inc) { return *this; }
+ Self& operator-=(size_t dec) { return *this; }
+ Self& operator++() { return *this; }
+ Self operator++(int) { return *this; }
+ Self& operator--() { return *this; }
+ Self operator--(int) { return *this; }
size_t size(const Container& c) const { return c.size(); }
void swap(SizeStorage& other) {}
};
--- a/include/utilities/counter.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/utilities/counter.h Sat Dec 27 14:33:25 2008 -0800
@@ -8,19 +8,19 @@
#ifndef COUNTERS
- #define GetCounter(path) 0
- #define Count(x)
- #define CountNum(x,y)
- #define CountBy(x,y)
- #define SetFrequency(x, freq)
- #define SetTrigger(x, y)
+ #define GetCounter(path) 0
+ #define Count(x)
+ #define CountNum(x,y)
+ #define CountBy(x,y)
+ #define SetFrequency(x, freq)
+ #define SetTrigger(x, y)
#else // COUNTERS
- #define GetCounter(path) get_counter(path)
- #define Count(x) do { x->count++; if ((x->count % x->frequency == 0)) x->trigger->print(); } while (0)
- #define CountNum(x,y) do { x->subcount[y]++; } while (0)
- #define CountBy(x,y) do { x->count += y; } while (0)
- #define SetFrequency(x, freq) do { x->frequency = freq; } while (0)
- #define SetTrigger(x, y) do { x->trigger = y; } while(0)
+ #define GetCounter(path) get_counter(path)
+ #define Count(x) do { x->count++; if ((x->count % x->frequency == 0)) x->trigger->print(); } while (0)
+ #define CountNum(x,y) do { x->subcount[y]++; } while (0)
+ #define CountBy(x,y) do { x->count += y; } while (0)
+ #define SetFrequency(x, freq) do { x->frequency = freq; } while (0)
+ #define SetTrigger(x, y) do { x->trigger = y; } while(0)
#endif // COUNTERS
@@ -32,46 +32,46 @@
class Counter
{
- public:
- typedef unsigned long CounterType;
- typedef std::map<std::string, Counter*> SubCounterMap;
- typedef std::map<int, CounterType> SubCountMap;
+ public:
+ typedef unsigned long CounterType;
+ typedef std::map<std::string, Counter*> SubCounterMap;
+ typedef std::map<int, CounterType> SubCountMap;
- public:
- CounterType count;
- CounterType frequency;
- SubCountMap subcount;
- Counter* trigger;
+ public:
+ CounterType count;
+ CounterType frequency;
+ SubCountMap subcount;
+ Counter* trigger;
- public:
- Counter(const std::string& full_name = "",
- CounterType freq = std::numeric_limits<CounterType>::max());
- ~Counter();
+ public:
+ Counter(const std::string& full_name = "",
+ CounterType freq = std::numeric_limits<CounterType>::max());
+ ~Counter();
- Counter* get_child(const std::string& path, std::string::size_type pos);
- void print();
+ Counter* get_child(const std::string& path, std::string::size_type pos);
+ void print();
- private:
- SubCounterMap subcounters_;
- std::string full_name_;
-
- static const char* start_color;
- static const char* finish_color;
- static const char green_color[];
- static const char normal_color[];
- static const char empty_string[];
+ private:
+ SubCounterMap subcounters_;
+ std::string full_name_;
+
+ static const char* start_color;
+ static const char* finish_color;
+ static const char green_color[];
+ static const char normal_color[];
+ static const char empty_string[];
};
-const char Counter::green_color[] = "\033[32m";
-const char Counter::normal_color[] = "\033[0m";
-const char Counter::empty_string[] = "";
-const char* Counter::start_color = 0;
-const char* Counter::finish_color = 0;
+const char Counter::green_color[] = "\033[32m";
+const char Counter::normal_color[] = "\033[0m";
+const char Counter::empty_string[] = "";
+const char* Counter::start_color = 0;
+const char* Counter::finish_color = 0;
-static Counter rootCounter;
+static Counter rootCounter;
-Counter* get_counter(const char* path)
+Counter* get_counter(const char* path)
{
- return rootCounter.get_child(path, 0);
+ return rootCounter.get_child(path, 0);
}
#include "counter.hpp"
--- a/include/utilities/types.h Thu Dec 25 14:24:02 2008 -0800
+++ b/include/utilities/types.h Sat Dec 27 14:33:25 2008 -0800
@@ -17,7 +17,7 @@
typedef const unsigned int& version_type;
struct Empty {};
-//std::ostream& operator<<(std::ostream& out, Empty e) { return out; }
+std::ostream& operator<<(std::ostream& out, Empty e) { return out; }
enum SwitchType
{