Adds auction algorithm computations and distance matrix filter.

.gitignore

@@ -1,4 +1,10 @@
-ComputeDistanceMatrix
 .DS_Store
 *.o
-
+/build
+tests/data/
+tests/distance*
+tests/subsample*
+tests/sample.json
+tests/filter*
+bin/ComputeSubsample
+bin/ComputeDistances

CMakeLists.txt

@@ -21,6 +21,7 @@ endif("${isSystemDir}" STREQUAL "-1")
 #########
 
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -mtune=native -ffast-math")
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -fprofile-arcs -ftest-coverage")
 set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 -DNDEBUG" )
 
@@ -37,6 +38,9 @@ find_package(MPI REQUIRED)
 
 include_directories ( 
   ./include 
+  ./include/persistence/approximatedistances/geom_bottleneck/bottleneck/include/ANN
+  ./include/persistence/approximatedistances/geom_bottleneck/bottleneck/include
+  ./include/persistence/approximatedistances/geom_matching/wasserstein/include
   /usr/local/include 
   /opt/local/include 
   ${Boost_INCLUDE_DIRS} 
@@ -64,6 +68,7 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 # Recurse #
 ###########
 
+add_subdirectory (include/persistence/approximatedistances)
 add_subdirectory (source)
 
 #########

doc/index.raw

@@ -62,8 +62,8 @@ The output of the subsample program will be stored in the supplied filename `/pa
 
 The input to the distance program is the output from the subsample program. The arguments are
 ```bash
-/path/to/subsample.json /path/to/distance.txt
+error /path/to/subsample.json /path/to/distance.txt /path/to/distance_filter.txt
 ```
-where the first is a path to the subsample (which contains a path to the original sample), and the second path is the location the distance matrix is to be stored.
+where the first is a nonnegative number giving the relative error for computing the approximate distance between the persistence diagrams (supply 0.0 for exact distance computation), the second a path to the subsample (which contains a path to the original sample), the third path is the location the distance matrix is to be stored, and the fourth is an optional path to a filter for the distance matrix computations. If supplied, this distance matrix filter should be a sequence in {0,1} that is the same size as the expected output of /path/to/distance.txt. If not supplied, all pairwise distances will be computed.
 
 

include/persistence/BottleneckApproximateDistance.h

@@ -0,0 +1,79 @@
+/// BottleneckApproximateDistance.h
+/// Author: Rachel Levanger
+/// Date: Nov 6, 2015
+/// Edit history: 
+///   - Created from BottleneckDistance.h to interface approximate algorithms with Subsample code (Rache Levanger, Nov 6, 2015)
+
+
+#ifndef BOTTLENECKAPPROXIMATEDISTANCE_H
+#define BOTTLENECKAPPROXIMATEDISTANCE_H
+
+#include <vector>
+#include "persistence/PersistenceDiagram.h"
+#include "persistence/approximatedistances/geom_bottleneck/bottleneck/include/bottleneck.h"
+
+double 
+BottleneckApproximateDistance( subsample::PersistenceDiagram const& diagram_1, 
+                    subsample::PersistenceDiagram const& diagram_2,
+                    const double epsilon );
+
+
+namespace BottleneckApproximateDistance_wrapper {
+
+struct BottleneckApproximationWrapper {
+
+    /* Generators for two persistence diagrams which are going to be compared */
+    std::vector<subsample::Generator> Generators1;
+    std::vector<subsample::Generator> Generators2;
+
+    /* Read the generators into the approximate algorithm class. */
+    bool populateDiagramPointSets(std::vector<std::pair<double, double>>& A,
+                              std::vector<std::pair<double, double>>& B)
+    {
+        A.clear();
+        B.clear();
+
+        /* Read in generators from Generators1 to A */
+        for ( std::vector<subsample::Generator>::const_iterator cur = Generators1.begin(); 
+          cur != Generators1.end(); ++cur ) {
+            A.push_back(std::make_pair(cur->birth, cur->death));
+        }
+        /* Read in generators from Generators2 to B */
+        for ( std::vector<subsample::Generator>::const_iterator cur = Generators2.begin(); 
+          cur != Generators2.end(); ++cur ) {
+            B.push_back(std::make_pair(cur->birth, cur->death));
+        }
+        return true;
+    }  
+
+};
+
+} //namespace
+
+inline double 
+BottleneckApproximateDistance( subsample::PersistenceDiagram const& diagram_1, 
+                    subsample::PersistenceDiagram const& diagram_2,
+                    const double epsilon ) {
+
+    using namespace BottleneckApproximateDistance_wrapper;
+
+    BottleneckApproximationWrapper bw;
+
+    std::vector<subsample::Generator> & Generators1 = bw . Generators1;
+    std::vector<subsample::Generator> & Generators2 = bw . Generators2;
+    Generators1 . assign ( diagram_1 . begin (), diagram_1 . end () );
+    Generators2 . assign ( diagram_2 . begin (), diagram_2 . end () );
+
+    std::vector<std::pair<double, double>> A, B;
+    if (!bw.populateDiagramPointSets(A, B)) {
+        std::cout << "Could not convert PersistenceDiagrams to DiagramPointSets.\n";
+        return -1;
+    }
+
+    double distance;
+    distance = geom_bt::bottleneckDistApprox(A, B, epsilon);
+    return distance;
+
+}
+
+#endif

include/persistence/BottleneckDistance.h

@@ -13,8 +13,8 @@
 #include "persistence/PersistenceDiagram.h"
 
 double 
-BottleneckDistance( PersistenceDiagram const& diagram_1, 
-                    PersistenceDiagram const& diagram_2 );
+BottleneckDistance( subsample::PersistenceDiagram const& diagram_1, 
+                    subsample::PersistenceDiagram const& diagram_2 );
 
 
 namespace BottleneckDistance_detail {
@@ -31,8 +31,8 @@ class Edge{
 struct BottleneckProblem {
 
   /* Generators for two persistence diagrams which are going to be compared */
-  std::vector<Generator> Generators1;
-  std::vector<Generator> Generators2;
+  std::vector<subsample::Generator> Generators1;
+  std::vector<subsample::Generator> Generators2;
   /* Number of generators in Generator1 and Generator2 */
   unsigned int  Max_Size;
   /* Edges between all the nodes given by generators of the persistence diagrams */
@@ -50,7 +50,7 @@ struct BottleneckProblem {
   std::vector<int> Layers;
 
   void PrepareEdges( void ){
-    Generator::Distance distance;
+    subsample::Generator::Distance distance;
     /*Set the number of generators */
     Max_Size = Generators1.size() + Generators2.size();
     /*Clear edges */
@@ -61,22 +61,22 @@ struct BottleneckProblem {
         Edges.push_back(Edge(i, j, 0));
     /* Edges between real points */
     unsigned int i = 0;
-    for ( std::vector<Generator>::const_iterator cur1 = Generators1.begin(); 
+    for ( std::vector<subsample::Generator>::const_iterator cur1 = Generators1.begin(); 
           cur1 != Generators1.end(); ++cur1 ) {
       unsigned int j = Max_Size;
-      for ( std::vector<Generator>::const_iterator cur2 = Generators2.begin(); 
+      for ( std::vector<subsample::Generator>::const_iterator cur2 = Generators2.begin(); 
             cur2 != Generators2.end(); ++cur2 ) {
         Edges.push_back(Edge(i,j++, distance(*cur1, *cur2)));
       }
       ++i;
     }
     /* Edges between real points and their corresponding diagonal points */
     i = 0;
-    for ( std::vector<Generator>::const_iterator cur1 = Generators1.begin(); 
+    for ( std::vector<subsample::Generator>::const_iterator cur1 = Generators1.begin(); 
           cur1 != Generators1.end(); ++cur1, ++i)
       Edges.push_back( Edge( i, Max_Size + Generators2.size() + i, distance . diagonal ( *cur1 ) ) );
     i = Max_Size;
-    for ( std::vector<Generator>::const_iterator cur2 = Generators2.begin(); 
+    for ( std::vector<subsample::Generator>::const_iterator cur2 = Generators2.begin(); 
           cur2 != Generators2.end(); ++cur2, ++i)
       Edges.push_back( Edge( Generators1.size() + (i - Max_Size), i, distance . diagonal ( *cur2 ) ) );
     std::sort(Edges.begin(), Edges.end());
@@ -156,23 +156,24 @@ struct BottleneckProblem {
 } //namespace
 
 inline double 
-BottleneckDistance( PersistenceDiagram const& diagram_1, 
-                    PersistenceDiagram const& diagram_2 ) {
+BottleneckDistance( subsample::PersistenceDiagram const& diagram_1, 
+                    subsample::PersistenceDiagram const& diagram_2 ) {
   using namespace BottleneckDistance_detail;
   BottleneckProblem bp;
   std::vector<int> & Pair = bp.Pair;
   std::vector<int> & Layers = bp.Layers;
-  std::vector< std::vector < int > > & Connections = bp.Connections;
-  std::vector<Generator> & Generators1 = bp . Generators1;
-  std::vector<Generator> & Generators2 = bp . Generators2;
+  std::vector< std::vector < int > > & Connections = bp . Connections;
+  std::vector<subsample::Generator> & Generators1 = bp . Generators1;
+  std::vector<subsample::Generator> & Generators2 = bp . Generators2;
   unsigned int & Max_Size = bp . Max_Size;
   std::vector<Edge> & Edges = bp . Edges;
   Generators1 . assign ( diagram_1 . begin (), diagram_1 . end () );
   Generators2 . assign ( diagram_2 . begin (), diagram_2 . end () );
   /* If both diagrams are empty the distance is 0 */
   if( Generators1.size() == 0 &&  Generators2.size() == 0 ) return 0;
-  bp.PrepareEdges ();
 
+  bp.PrepareEdges ();
+
   /* Clear the pairing */
   Pair.clear( );
   Pair.assign( 2*Max_Size, -1);

include/persistence/PersistenceDiagram.h

@@ -20,6 +20,9 @@
 #include <boost/serialization/base_object.hpp>
 #include <boost/serialization/vector.hpp>
 
+namespace subsample
+{
+
 int64_t generator_distance_count = 0;
 struct Generator {
   double birth;
@@ -80,10 +83,12 @@ class PersistenceDiagram : private std::vector<Generator> {
   using Base::const_iterator;
   using Base::operator[];
   void load ( std::string const& filename ) {
+
     clear ();
     std::ifstream infile ( filename );
-    if ( not infile . good () ) {
-      throw std::runtime_error("PersistenceDiagram::load. File not found: " + filename );
+    if (not infile . good ()) { 
+      std::cout << "PersistenceDiagram::load. File not found. \n";
+      throw std::runtime_error("PersistenceDiagram::load. File not found.");
     }
     std::string line;
     while ( std::getline ( infile, line ) ) {
@@ -98,6 +103,7 @@ class PersistenceDiagram : private std::vector<Generator> {
     }
     infile . close ();
     // Replace -1's with something sensible
+
     double max_entry = 0;
     for ( int i = 0; i < size (); ++ i ) {
       max_entry = std::max ( (*this)[i].birth, max_entry );
@@ -121,4 +127,6 @@ class PersistenceDiagram : private std::vector<Generator> {
   }
 };
 
+}
+
 #endif

include/persistence/WassersteinApproximateDistance.h

@@ -0,0 +1,76 @@
+#ifndef WASSERSTEINAPPROXIMATEDISTANCE_H
+#define WASSERSTEINAPPROXIMATEDISTANCE_H
+#include <cmath>
+#include <cstring>
+#include <vector>
+#include <algorithm>
+
+#include "persistence/PersistenceDiagram.h"
+#include "persistence/approximatedistances/geom_matching/wasserstein/include/wasserstein.h"
+
+double 
+WassersteinApproximateDistance( subsample::PersistenceDiagram const& diagram_1, 
+                     subsample::PersistenceDiagram const& diagram_2,
+                     const double p,
+                     const double epsilon );
+
+
+namespace WassersteinApproximateDistance_wrapper {
+
+struct WassersteinApproximationWrapper {
+    /* Generators for two persistence diagrams which are going to be compared */
+    std::vector<subsample::Generator> Generators1;
+    std::vector<subsample::Generator> Generators2;
+
+    /* Read the generators into the approximate algorithm class. */
+        bool populateDiagramPointSets(std::vector<std::pair<double, double>>& A,
+                              std::vector<std::pair<double, double>>& B)
+    {
+        A.clear();
+        B.clear();
+        /* Read in generators from Generators1 to A */
+        for ( std::vector<subsample::Generator>::const_iterator cur = Generators1.begin(); 
+          cur != Generators1.end(); ++cur ) {
+            A.push_back(std::make_pair(cur->birth, cur->death));
+        }
+        /* Read in generators from Generators2 to B */
+        for ( std::vector<subsample::Generator>::const_iterator cur = Generators2.begin(); 
+          cur != Generators2.end(); ++cur ) {
+            B.push_back(std::make_pair(cur->birth, cur->death));
+        }
+        return true;
+    }  
+
+};
+
+} //namespace
+
+
+inline double 
+WassersteinApproximateDistance( subsample::PersistenceDiagram const& diagram_1, 
+                    subsample::PersistenceDiagram const& diagram_2,
+                    const double p,
+                    const double epsilon ) {
+
+    using namespace WassersteinApproximateDistance_wrapper;
+
+    WassersteinApproximationWrapper ww;
+
+    std::vector<subsample::Generator> & Generators1 = ww . Generators1;
+    std::vector<subsample::Generator> & Generators2 = ww . Generators2;
+    Generators1 . assign ( diagram_1 . begin (), diagram_1 . end () );
+    Generators2 . assign ( diagram_2 . begin (), diagram_2 . end () );
+
+    std::vector<std::pair<double, double>> A, B;
+    if (!ww.populateDiagramPointSets(A, B)) {
+        std::cout << "Could not convert PersistenceDiagrams to DiagramPointSets.\n";
+        return -1;
+    }
+
+    double distance;
+    distance = geom_ws::wassersteinDist(A, B, p, epsilon);
+    return distance;
+
+}
+
+#endif

include/persistence/WassersteinDistance.h

@@ -7,8 +7,8 @@
 #include "PersistenceDiagram.h"
 
 double 
-WassersteinDistance( PersistenceDiagram const& diagram_1, 
-                     PersistenceDiagram const& diagram_2,
+WassersteinDistance( subsample::PersistenceDiagram const& diagram_1, 
+                     subsample::PersistenceDiagram const& diagram_2,
                      double p );
 
 
@@ -142,11 +142,11 @@ namespace WassersteinDistance_detail {
 }
 
 inline double 
-WassersteinDistance( PersistenceDiagram const& diagram_1, 
-                     PersistenceDiagram const& diagram_2,
+WassersteinDistance( subsample::PersistenceDiagram const& diagram_1, 
+                     subsample::PersistenceDiagram const& diagram_2,
                      double p ) {
   using namespace WassersteinDistance_detail;
-  Generator::Distance distance; 
+  subsample::Generator::Distance distance; 
   // Distance matrix the generators in diagram_1 and diagram_2 
   double* distanceMatrix;
   // Distance matrix has (diagram_1.size +  diagram_2.size())^2 elements 

include/persistence/approximatedistances/CMakeLists.txt

@@ -0,0 +1,7 @@
+###########
+# Recurse #
+###########
+
+add_subdirectory (geom_bottleneck)
+add_subdirectory (geom_matching)
+

include/persistence/approximatedistances/geom_bottleneck/CMakeLists.txt

@@ -0,0 +1,6 @@
+###########
+# Recurse #
+###########
+
+add_subdirectory (bottleneck)
+