--- a/.issues/8a14b4849071f910/new/1221008555.M667846P30017Q17.cole Thu Dec 25 13:09:00 2008 -0800
+++ b/.issues/8a14b4849071f910/new/1221008555.M667846P30017Q17.cole Sat Dec 27 14:51:38 2008 -0800
@@ -1,8 +1,9 @@
From: Dmitriy Morozov <morozov@cs.duke.edu>
Date: Thu, 10 Jan 2008 04:36:03 -0500
-State: new
+State: resolved
Subject: Remove maintenance of "lazy decomposition"
Message-Id: <8a14b4849071f910-0-artemis@metatron>
+resolution: fixed
The maintenance of "lazy decomposition" (added in [a0736dd3c671]) is
incorrect (due to original theoretical errors). Remove it completely.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/8a14b4849071f910/new/1229837212.M576947P31212Q2.rufus Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,9 @@
+From: Dmitriy Morozov <dmitriy@mrzv.org>
+Date: Sat, 20 Dec 2008 21:26:52
+Subject: properties changes (state, resolution)
+Message-Id: <8a14b4849071f910-c82a77e7b320c7e8-artemis@rufus>
+References: <8a14b4849071f910-0-artemis@metatron>
+In-Reply-To: <8a14b4849071f910-0-artemis@metatron>
+
+state=resolved
+resolution=fixed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/957a589d7c6c3fa8/new/1230356422.M300560P11955Q1.rufus Sat Dec 27 14:51:38 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/.issues/96b51b44f7764f5c/new/1221008595.M830007P30168Q1.cole Thu Dec 25 13:09:00 2008 -0800
+++ b/.issues/96b51b44f7764f5c/new/1221008595.M830007P30168Q1.cole Sat Dec 27 14:51:38 2008 -0800
@@ -1,7 +1,8 @@
From: Dmitriy Morozov <dmitriy@mrzv.org>
Date: Tue, 09 Sep 2008 18:03:00
-State: new
+State: resolved
Subject: Boost 1.36
Message-Id: <96b51b44f7764f5c-0-artemis@cole>
+resolution: fixed
Make the software compile with Boost 1.36
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/96b51b44f7764f5c/new/1229837196.M796227P31208Q2.rufus Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,9 @@
+From: Dmitriy Morozov <dmitriy@mrzv.org>
+Date: Sat, 20 Dec 2008 21:26:36
+Subject: properties changes (state, resolution)
+Message-Id: <96b51b44f7764f5c-85f856488e2d4827-artemis@rufus>
+References: <96b51b44f7764f5c-0-artemis@cole>
+In-Reply-To: <96b51b44f7764f5c-0-artemis@cole>
+
+state=resolved
+resolution=fixed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/de674a2ac5f6c18c/new/1229837097.M859247P31105Q1.rufus Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,10 @@
+From: Dmitriy Morozov <dmitriy@mrzv.org>
+Date: Sat, 20 Dec 2008 21:21:53
+State: new
+Subject: Specialize ChainWrapper<C>::add() for linked lists
+Message-Id: <de674a2ac5f6c18c-0-artemis@rufus>
+
+The current implementation of add() in ChainWrapper is container agnostic: it
+uses a temporary container, and then swaps it into place. There should be a
+specialized add() for linked lists (in particular, List from circular_list.h)
+since it is would likely be more efficient then the generic one.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/f6496b3e37275888/new/1229837161.M623187P31203Q1.rufus Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,10 @@
+From: Dmitriy Morozov <dmitriy@mrzv.org>
+Date: Sat, 20 Dec 2008 21:24:59
+State: new
+Subject: Add field arithmetic
+Message-Id: <f6496b3e37275888-0-artemis@rufus>
+
+Add support for field arithmetic (i.e., not only Z_2 field as it is now).
+However, make sure that there is a specialization of ChainWrapper<C>::add() for
+Z_2 field (since it can be made much more efficient both in terms of space and
+time).
--- a/CMakeLists.txt Thu Dec 25 13:09:00 2008 -0800
+++ b/CMakeLists.txt Sat Dec 27 14:51:38 2008 -0800
@@ -1,107 +1,108 @@
-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 python serialization signals)
-find_package (Doxygen)
+find_package (Boost REQUIRED COMPONENTS program_options python 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 (examples)
+add_subdirectory (tests)
+add_subdirectory (tools)
add_subdirectory (bindings)
--- a/bindings/python/simplex.cpp Thu Dec 25 13:09:00 2008 -0800
+++ b/bindings/python/simplex.cpp Sat Dec 27 14:51:38 2008 -0800
@@ -57,7 +57,7 @@
template<class V, class T>
int data_dimension_comparison(const Simplex<V,T>& a, const Simplex<V,T>& b)
{
- return ThreeOutcomeCompare<typename Simplex<V,T>::DataDimensionComparison>().compare(a,b);
+ return ThreeOutcomeCompare<DataDimensionComparison<Simplex<V,T> > >().compare(a,b);
}
#include "python-simplex.h" // defines SimplexVD, Vertex, and Data
--- a/examples/CMakeLists.txt Thu Dec 25 13:09:00 2008 -0800
+++ b/examples/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/alphashapes/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/cech-complex/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/cech-complex/cech-complex.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/fitness/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/fitness/avida-distance.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/fitness/avida-rips-distance.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/poincare/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/poincare/poincare.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/rips/CMakeLists.txt Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/rips/rips.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/examples/triangle/triangle.cpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/topology/chain.h Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/topology/chain.hpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/topology/conesimplex.h Sat Dec 27 14:51:38 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/filtration.h Thu Dec 25 13:09:00 2008 -0800
+++ b/include/topology/filtration.h Sat Dec 27 14:51:38 2008 -0800
@@ -10,6 +10,10 @@
// Class: Filtration
//
+// Filtration keeps track of the ordering of the simplices in a complex.
+// The most significant function it provides is <boundary()> which converts
+// the boundary of a simplex at a given index into a list of indices.
+//
// TODO: this is really specialized for an std::vector<> Complex; eventually generalize
// TODO: should we derive from Order?
template<class Complex_,
--- a/include/topology/filtrationcontainer.h Thu Dec 25 13:09:00 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 13:09:00 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 13:09:00 2008 -0800
+++ b/include/topology/lowerstarfiltration.h Sat Dec 27 14:51:38 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/lowerstarfiltration.hpp Thu Dec 25 13:09:00 2008 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,226 +0,0 @@
-#include "utilities/counter.h"
-
-/* Implementations */
-
-#ifdef LOGGING
-static rlog::RLogChannel* rlLowerStar = DEF_CHANNEL("lowerstar", rlog::Log_Debug);
-#endif // LOGGING
-
-#ifdef COUNTERS
-static Counter* cLowerStarTransposition = GetCounter("lowerstar");
-static Counter* cLowerStarChangedAttachment = GetCounter("lowerstar/ChangedAttachment");
-#endif // COUNTERS
-
-
-/**
- * Copies vertices [begin, end), orders them by cmp, and
- * records which vineyard to use in case of transpositions.
- */
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-template<class VertexCmp>
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-LowerStarFiltration(VertexIndex begin, VertexIndex end, const VertexCmp& cmp, Vineyard* vineyard):
- Parent(vineyard)
-{
- // Record VertexIndexes in a temporary list
- typedef std::list<VertexIndex> VertexIndexList;
- VertexIndexList tmp_list;
- while (begin != end)
- tmp_list.push_back(begin++);
-
- // Sort the temporary list
- tmp_list.sort(cmp);
-
- // Record vertex order
- for(typename VertexIndexList::const_iterator cur = tmp_list.begin(); cur != tmp_list.end(); ++cur)
- (*cur)->set_order(vertex_order.push_back(VertexDescriptor(*cur, Parent::append(Simplex(0, *cur)))));
-}
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-typename LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::Index
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-append(const Simplex& s)
-{
- AssertMsg(s.dimension() != 0, "All vertices must have been inserted in the constructor");
-
- // Find the highest vertex
- typename Simplex::VertexContainer::const_iterator cur = s.vertices().begin(), max = cur++;
- for (; cur != s.vertices().end(); ++cur)
- if (!vertex_cmp((*cur)->get_order(), (*max)->get_order()))
- max = cur;
-
- Index ms = (*max)->get_order()->simplex_index; Index prior;
- do { prior = ms++; } while (ms->dimension() <= s.dimension() && ms != Parent::end() && ms->get_attachment() == *max);
-
- Index i = Parent::insert(prior, s);
- i->set_attachment(*max);
-
- return i;
-}
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-bool
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::SimplexAttachmentComparison::
-operator()(const Simplex& first, const Simplex& second) const
-{
- int cmp = vertex_cmp.compare(first.get_attachment()->get_order(), second.get_attachment()->get_order());
- if (cmp == 0)
- return first.dimension() < second.dimension();
- else
- return cmp == -1;
-}
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-bool
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-transpose_vertices(const VertexOrderIndex& order)
-{
- Count(cLowerStarTransposition);
- rLog(rlLowerStar, "Transposing vertices (%s, %s)",
- tostring(order->vertex_index).c_str(),
- tostring(boost::next(order)->vertex_index).c_str());
-
- Index i = order->simplex_index;
- Index i_prev = boost::prior(i);
- Index i_next = boost::next(order)->simplex_index;
- Index i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
- Index j = boost::next(i_next);
-
- const VertexIndex& v_i = order->vertex_index;
- const VertexIndex& v_i_next = boost::next(order)->vertex_index;
- bool nbghrs = neighbors(v_i, v_i_next);
-
- bool result = false;
-
- // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
- while (i_next_prev != i_prev)
- {
- result |= transpose(i_next_prev);
- i_next_prev = boost::prior(i_next);
- }
- rLog(rlLowerStar, "Done moving the vertex");
-
- // Second, move the simplices attached to it
- rLog(rlLowerStar, "Moving attached simplices");
- while (j != Parent::end() && j->get_attachment() == v_i_next)
- {
- rLog(rlLowerStar, " Considering %s", tostring(*j).c_str());
- if (nbghrs && j->contains(v_i)) // short circuit
- {
- Count(cLowerStarChangedAttachment);
- rLog(rlLowerStar, " Attachment changed for %s", tostring(*j).c_str());
- j->set_attachment(v_i);
- ++j;
- continue;
- }
-
- Index j_prev = j; ++j;
- while ((--j_prev)->get_attachment() != v_i_next) // i.e., until we have reached v_i_next
- // (and the simplices that follow it) again
- {
- rLog(rlLowerStar, " Moving: %s, %s",
- tostring(*j_prev).c_str(), tostring(*boost::next(j_prev)).c_str());
- AssertMsg(j_prev->get_attachment() == v_i, "Simplex preceding the one being moved must be attached to v_i");
- result |= transpose(j_prev);
- --j_prev;
- }
- }
- rLog(rlLowerStar, "Done moving attached simplices");
- vertex_order.swap(order, boost::next(order));
-
- return result;
-}
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-bool
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-transpose(Index i)
-{
- Index j = boost::next(i);
-
- rLog(rlLowerStar, " Transposing (%s, %s) and (%s, %s)",
- tostring(*i).c_str(), tostring(*(i->pair())).c_str(),
- tostring(*j).c_str(), tostring(*(j->pair())).c_str());
-
- assert_pairing(i);
- assert_pairing(j);
-
- bool res = Parent::transpose(i, false);
- rLog(rlLowerStar, " %s: %s, %s: %s",
- tostring(*j).c_str(), tostring(*(j->pair())).c_str(),
- tostring(*i).c_str(), tostring(*(i->pair())).c_str());
-
- assert_pairing(j);
- assert_pairing(i);
-
- return res;
-}
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-void
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-assert_pairing(Index i)
-{
-#ifndef NDEBUG
- AssertMsg(i != Parent::end(), "Cannot assert pairing of end()");
- if (!i->sign())
- {
- if (i->pair() != i->cycle().top(Parent::get_cycles_cmp()))
- {
- rLog(rlLowerStar, "i (negative): %s", tostring(*i).c_str());
- rLog(rlLowerStar, "pair(i): %s", tostring(*(i->pair())).c_str());
- rLog(rlLowerStar, "i->cycle().top(): %s",
- tostring(*(i->cycle().top(Parent::get_cycles_cmp()))).c_str());
- AssertMsg(0, "Pairing not matching the matrix at %s", tostring(*i).c_str());
- }
- } else
- {
- if (i->pair() != i)
- {
- if (i->pair()->cycle().top(Parent::get_cycles_cmp()) != i)
- {
- rLog(rlLowerStar, "i (positive): %s", tostring(*i).c_str());
- rLog(rlLowerStar, "pair(i): %s", tostring(*(i->pair())).c_str());
- rLog(rlLowerStar, "pair(i)->cycle(): %s", tostring(i->pair()->cycle()).c_str());
- rLog(rlLowerStar, "pair->cycle().top(): %s", tostring(*(i->pair()->cycle().top(Parent::get_cycles_cmp()))).c_str());
- AssertMsg(0, "Pairing not matching the matrix at %s", tostring(*(i->pair())).c_str());
- }
- }
- }
-#endif
-}
-
-
-template<class VI, class Smplx, class FltrSmplx, class Vnrd>
-template<class Archive>
-void
-LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
-load(Archive& ar, version_type )
-{
-/*
- ar >> BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
-
- // Count the number of vertices
- VertexIndex num_vertices = 0;
- for (Index cur = begin(); cur != end(); ++cur)
- if (dimension(cur) == 0)
- num_vertices++;
-
- // Second pass to record vertex_order
- vertex_order.resize(num_vertices);
- inverse_vertex_order.resize(num_vertices);
- num_vertices = 0;
- for (Index cur = begin(); cur != end(); ++cur)
- {
- if (dimension(cur) == 0)
- {
- vertex_order[num_vertices].index = cur;
- vertex_order[num_vertices].vertex_index = *(cur->get_vertices().begin());
- inverse_vertex_order[vertex_order[num_vertices].vertex_index] = num_vertices;
- ++num_vertices;
- }
- }
-*/
-}
-
-
--- a/include/topology/persistence-diagram.h Thu Dec 25 13:09:00 2008 -0800
+++ b/include/topology/persistence-diagram.h Sat Dec 27 14:51:38 2008 -0800
@@ -15,12 +15,14 @@
*
* Stores birth-death pair plus any additional information provided by `Data` template parameter.
*/
-template<class Data_>
+template<class Data_ = Empty<> >
class PDPoint
{
public:
typedef Data_ Data;
+ PDPoint(const PDPoint& other):
+ point_(other.point_) {}
PDPoint(RealType x = 0, RealType y = 0, const Data& data = Data());
RealType x() const { return point_.first().first; }
@@ -28,7 +30,7 @@
const Data& data() const { return point_.second(); }
Data& data() { return point_.second(); }
- std::ostream& operator<<(std::ostream& out) const { return (out << x() << " " << y() << " " << data()); }
+ std::ostream& operator<<(std::ostream& out) const { return (out << x() << " " << y()); } // << " " << data()); }
struct Visitor
{
@@ -79,6 +81,9 @@
PersistenceDiagram() {}
+ template<class OtherData>
+ PersistenceDiagram(const PersistenceDiagram<OtherData>& other);
+
template<class Iterator, class Evaluator>
PersistenceDiagram(Iterator bg, Iterator end,
const Evaluator& eval = Evaluator());
--- a/include/topology/persistence-diagram.hpp Thu Dec 25 13:09:00 2008 -0800
+++ b/include/topology/persistence-diagram.hpp Sat Dec 27 14:51:38 2008 -0800
@@ -12,6 +12,18 @@
point_.second() = data;
}
+
+template<class D>
+template<class OtherData>
+PersistenceDiagram<D>::
+PersistenceDiagram(const PersistenceDiagram<OtherData>& other)
+{
+ points_.reserve(other.size());
+ for (typename PersistenceDiagram<OtherData>::PointVector::const_iterator cur = points_.begin();
+ cur != points_.end(); ++cur)
+ push_back(Point(cur->x(), cur->y()));
+}
+
template<class D>
template<class Iterator, class Evaluator>
PersistenceDiagram<D>::
--- a/include/topology/rips.h Thu Dec 25 13:09:00 2008 -0800
+++ b/include/topology/rips.h Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/topology/rips.hpp Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/topology/simplex.h Sat Dec 27 14:51:38 2008 -0800
@@ -33,7 +33,7 @@
template<class V, class T = Empty<> >
class Simplex
{
- public:
+ public:
/* Typedefs: Public Types
*
* Vertex - vertex type (template parameter V)
@@ -41,49 +41,44 @@
* Self -
* Boundary - type in which the boundary is stored
*/
- typedef V Vertex;
+ typedef V Vertex;
typedef T Data;
- typedef Simplex<Vertex, Data> Self;
- class BoundaryIterator;
+ typedef Simplex<Vertex, Data> Self;
+ class BoundaryIterator;
/* 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
+ /// @{
///
/// Functions: boundary_begin(), boundary_end()
/// Returns the iterators over the boundary of the simplex
@@ -91,25 +86,26 @@
BoundaryIterator boundary_end() 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()); }
- /// @}
+ /// @}
- 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
*
@@ -128,18 +124,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 );
};
@@ -163,20 +159,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 13:09:00 2008 -0800
+++ b/include/topology/simplex.hpp Sat Dec 27 14:51:38 2008 -0800
@@ -57,16 +57,27 @@
return BoundaryIterator(vertices().end(), vertices());
}
+#if 0
template<class V, class T>
bool
Simplex<V,T>::
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));
}
-
+#endif
+
+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>::
@@ -75,7 +86,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
@@ -87,27 +98,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/topology/static-persistence.hpp Thu Dec 25 13:09:00 2008 -0800
+++ b/include/topology/static-persistence.hpp Sat Dec 27 14:51:38 2008 -0800
@@ -23,7 +23,7 @@
ocmp_(ocmp)
{
OrderIndex ocur = begin();
- OffsetMap<size_t, OrderIndex> om(0, ocur); // TODO: this is customized for std::vector Order
+ OffsetMap<typename Filtration::IntermediateIndex, OrderIndex> om(0, ocur); // TODO: this is customized for std::vector Order
for (typename Filtration::Index cur = filtration.begin(); cur != filtration.end(); ++cur, ++ocur)
{
// Convert the Filtration::IndexBoundary into a Cycle, and
--- a/include/utilities/containers.h Thu Dec 25 13:09:00 2008 -0800
+++ b/include/utilities/containers.h Sat Dec 27 14:51:38 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 13:09:00 2008 -0800
+++ b/include/utilities/counter.h Sat Dec 27 14:51:38 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)
#include <map>
@@ -31,46 +31,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"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/CMakeLists.txt Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,11 @@
+find_package (Qt4 REQUIRED)
+set (QT_USE_QTOPENGL TRUE)
+set (QT_USE_QTXML TRUE)
+include (${QT_USE_FILE})
+
+#find_library (gle_LIBRARY NAMES gle)
+#find_library (QGLViewer_LIBRARY NAMES QGLViewer)
+#find_path (QGLViewer_INCLUDE_DIR QGLViewer/qglviewer.h)
+#include_directories (${QGLViewer_INCLUDE_DIR})
+
+add_subdirectory (diagram-viewer)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/diagram-viewer/CMakeLists.txt Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,18 @@
+set (diagram-viewerSources
+ diagram.cpp
+ diagram-viewer-main.cpp)
+
+set (diagram-viewerHeaders
+ diagram.h)
+
+qt4_wrap_cpp (diagram-viewerMocSources ${diagram-viewerHeaders})
+
+set (libraries ${libraries}
+ ${Boost_SERIALIZATION_LIBRARY}
+ ${Boost_PROGRAM_OPTIONS_LIBRARY}
+ ${QT_LIBRARIES})
+
+add_executable (diagram-viewer ${diagram-viewerSources}
+ ${diagram-viewerMocSources})
+
+target_link_libraries (diagram-viewer ${libraries})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/diagram-viewer/diagram-viewer-main.cpp Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,59 @@
+#include <qapplication.h>
+#include <QtGui>
+
+#include "diagram.h"
+
+#include <fstream>
+#include <map>
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/serialization/map.hpp>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+
+int main (int argc, char *argv[])
+{
+ std::string diagrams_filename;
+ int dimension;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("diagrams-file", po::value<std::string>(&diagrams_filename), "The collection of persistence diagrams")
+ ("dimension", po::value<int>(&dimension), "Dimension of the diagram to show");
+
+ po::positional_options_description p;
+ p.add("diagrams-file", 1);
+ p.add("dimension", 2);
+
+ po::options_description all; all.add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(p).run(), vm);
+ po::notify(vm);
+
+ if (!vm.count("diagrams-file") || !vm.count("dimension"))
+ {
+ std::cout << "Usage: " << argv[0] << " diagrams-file dimension" << std::endl;
+ std::cout << hidden << std::endl;
+ return 1;
+ }
+
+ std::map<Dimension, PDiagram> dgms;
+ std::ifstream ifs(diagrams_filename.c_str());
+ boost::archive::binary_iarchive ia(ifs);
+ ia >> dgms;
+
+
+ QApplication application(argc, argv);
+
+ std::cout << dimension << std::endl;
+ std::cout << dgms[dimension] << std::endl;
+
+ DgmViewer pd(dgms[dimension]);
+ pd.show();
+
+ // Run main loop.
+ return application.exec();
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/diagram-viewer/diagram.cpp Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,95 @@
+#include <iostream>
+
+#include <QtGui>
+#include <QRectF>
+
+#include "diagram.h"
+#include <cmath>
+
+//static const double ellipse_size = 0.035;
+static const double ellipse_size = 3;
+
+/* DgmViewer Implementation */
+DgmViewer::DgmViewer(const PDiagram& dgm):
+ min_x(0), min_y(0), max_x(0), max_y(0)
+{
+ points.reserve(dgm.size());
+
+ for (PDiagram::const_iterator cur = dgm.begin(); cur != dgm.end(); ++cur)
+ {
+ min_x = std::min(min_x, cur->x());
+ min_y = std::min(min_y, cur->y());
+ max_x = std::max(max_x, cur->x());
+ max_y = std::max(max_y, cur->y());
+
+ points.push_back(new DgmPoint(*cur, ellipse_size));
+ }
+
+ addDgmPoints();
+ setWindowTitle(QString("Persistence Diagram"));
+}
+DgmViewer::~DgmViewer()
+{
+ for (PointsVector::iterator cur = points.begin(); cur != points.end(); ++cur)
+ delete *cur;
+}
+
+
+void DgmViewer::addDgmPoints()
+{
+ RealType min = std::min(min_x, min_y);
+ RealType max = std::max(max_x, max_y);
+
+ QGraphicsLineItem* diagonal = new QGraphicsLineItem(QLineF(min, -min, max, -max));
+ QGraphicsLineItem* y_axis = new QGraphicsLineItem(QLineF(0, -min_y, 0, -max_y));
+ QGraphicsLineItem* x_axis = new QGraphicsLineItem(QLineF(min_x, 0, max_x, 0));
+
+ scene.addItem(diagonal);
+ scene.addItem(y_axis);
+ scene.addItem(x_axis);
+
+ for (PointsVector::const_iterator cur = points.begin(); cur != points.end(); ++cur)
+ scene.addItem(*cur);
+
+ //scale(100,100);
+ setScene(&scene);
+ setRenderHint(QPainter::Antialiasing);
+ ensureVisible(scene.itemsBoundingRect());
+ //setMinimumSize( (int)(maxX - minX)*100 + 100, (int) (maxY - minY)*100 + 100);
+}
+
+
+DgmPoint::DgmPoint(QGraphicsItem* parent):
+ QGraphicsItem(parent)
+{
+}
+
+DgmPoint::DgmPoint(const Parent& pt, qreal size, QGraphicsItem *parent):
+ Parent(pt), ellipse_size(size), QGraphicsItem(parent)
+{
+ setToolTip(QString("(%1, %2)").arg(getX()).arg(getY()));
+}
+
+DgmPoint::DgmPoint(RealType b, RealType d, qreal size, QGraphicsItem *parent):
+ Parent(b, d), ellipse_size(size), QGraphicsItem(parent)
+{
+ setToolTip(QString("(%1, %2)").arg(getX()).arg(getY()));
+}
+
+void DgmPoint::paint(QPainter *painter, const QStyleOptionGraphicsItem *option, QWidget *widget)
+{
+ Q_UNUSED(option);
+ Q_UNUSED(widget);
+
+ //QBrush solidFill(unselectColor);
+ //QBRush selectSolidFill(selectColor);
+ painter->setBrush(Qt::SolidPattern);
+ //painter->setPen(selectColor);
+ painter->drawEllipse(QRectF(getX() - ellipse_size, -getY() - ellipse_size, 2*ellipse_size, 2*ellipse_size));
+}
+
+
+QRectF DgmPoint::boundingRect() const
+{
+ return QRectF(getX() - ellipse_size, -getY() - ellipse_size, 2*ellipse_size, 2*ellipse_size);
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/diagram-viewer/diagram.h Sat Dec 27 14:51:38 2008 -0800
@@ -0,0 +1,60 @@
+#ifndef __DIAGRAM_H__
+#define __DIAGRAM_H__
+
+#include <QtGui>
+#include <QObject>
+#include <QColor>
+
+#include <map>
+
+#include <utilities/types.h>
+#include <topology/persistence-diagram.h>
+
+typedef PersistenceDiagram<> PDiagram;
+typedef std::map<Dimension, PDiagram> Diagrams;
+
+class DgmPoint;
+
+class DgmViewer: public QGraphicsView
+{
+ Q_OBJECT
+
+ public:
+ typedef std::vector<DgmPoint*> PointsVector;
+
+ DgmViewer(const PDiagram& dgm);
+ ~DgmViewer();
+
+ void addDgmPoints();
+
+ private:
+ PointsVector points;
+ QGraphicsScene scene;
+ RealType min_x, min_y, max_x, max_y;
+};
+
+
+class DgmPoint: public PDPoint<>, public QGraphicsItem
+{
+ public:
+ typedef PDPoint<> Parent;
+
+ DgmPoint(QGraphicsItem* parent = 0);
+ DgmPoint(const Parent& pt, qreal size, QGraphicsItem *parent = 0);
+ DgmPoint(RealType b, RealType d, qreal size, QGraphicsItem *parent = 0);
+
+ void paint(QPainter *painter, const QStyleOptionGraphicsItem *option, QWidget *widget = 0);
+ QRectF boundingRect() const;
+
+ qreal getX() const { return Parent::x(); }
+ qreal getY() const { return Parent::y(); }
+
+ int type() const { return QGraphicsItem::UserType + 1; }
+
+ private:
+ // size of rectangle containing ellipses
+ qreal ellipse_size;
+};
+
+
+#endif // __DIAGRAM_H__