Prepare to output real relative error.

In Wasserstein distance computation AuctionRunner class now has a field
relativeError which contains the relative error which we can guarantee,
so if the user asked for 0.1 accuracy, but the result is 0.03 accurate,
this information can be retrieved.

6 files changed, 14 insertions, 4 deletions

diff --git a/geom_matching/wasserstein/example/wasserstein_dist.cpp b/geom_matching/wasserstein/example/wasserstein_dist.cpp
index e92ab54..c699282 100644
--- a/geom_matching/wasserstein/example/wasserstein_dist.cpp
+++ b/geom_matching/wasserstein/example/wasserstein_dist.cpp
@@ -46,7 +46,7 @@ int main(int argc, char* argv[])
     PairVector diagramA, diagramB;
 
     if (argc < 3 ) {
-        std::cerr << "Usage: " << argv[0] << " file1 file2 [wasserstein_degree] [relative_error] [internal norm]. By default power is 1.0, relative error is 0.01, internal norm is l_infinity." << std::endl;
+        std::cerr << "Usage: " << argv[0] << " file1 file2 [wasserstein_degree] [relative_error] [internal norm] [output_actual_error]. By default power is 1.0, relative error is 0.01, internal norm is l_infinity, actual relative error is not printed." << std::endl;
         return 1;
     }
 
diff --git a/geom_matching/wasserstein/include/auction_runner_gs.h b/geom_matching/wasserstein/include/auction_runner_gs.h
index 80aa9f0..7968fa9 100644
--- a/geom_matching/wasserstein/include/auction_runner_gs.h
+++ b/geom_matching/wasserstein/include/auction_runner_gs.h
@@ -79,6 +79,7 @@ public:
     double getEpsilon() const { return epsilon; }
     double getWassersteinDistance();
     double getWassersteinCost();
+    double getRelativeError() const { return relativeError; };
     static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling
 private:
     // private data
@@ -96,6 +97,7 @@ private:
     double weightAdjConst;
     double wassersteinDistance;
     double wassersteinCost;
+    double relativeError;
     // to get the 2 best items
     std::unique_ptr<AuctionOracle> oracle;
 #ifdef KEEP_UNASSIGNED_ORDERED
diff --git a/geom_matching/wasserstein/include/auction_runner_jac.h b/geom_matching/wasserstein/include/auction_runner_jac.h
index 22d42b0..524498a 100644
--- a/geom_matching/wasserstein/include/auction_runner_jac.h
+++ b/geom_matching/wasserstein/include/auction_runner_jac.h
@@ -52,6 +52,7 @@ public:
     double getEpsilon() const { return epsilon; }
     double getWassersteinDistance();
     double getWassersteinCost();
+    double getRelativeError() const { return relativeError; };
     static constexpr double epsilonCommonRatio { 5 }; // next epsilon = current epsilon / epsilonCommonRatio
     static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling
 private:
@@ -69,6 +70,7 @@ private:
     double wassersteinDistance;
     double wassersteinCost;
     std::vector<IdxValPair> bidTable;
+    double relativeError;
     // to get the 2 best items
     std::unique_ptr<AuctionOracle> oracle;
     std::list<size_t> unassignedBidders;
diff --git a/geom_matching/wasserstein/src/auction_runner_gs.cpp b/geom_matching/wasserstein/src/auction_runner_gs.cpp
index bd25442..5865325 100644
--- a/geom_matching/wasserstein/src/auction_runner_gs.cpp
+++ b/geom_matching/wasserstein/src/auction_runner_gs.cpp
@@ -122,6 +122,7 @@ void AuctionRunnerGS::flushAssignment(void)
 
 void AuctionRunnerGS::runAuction(void)
 {
+    relativeError = std::numeric_limits<double>::max();
 #ifdef PRINT_DETAILED_TIMING
     std::chrono::high_resolution_clock hrClock;
     std::chrono::time_point<std::chrono::high_resolution_clock> startMoment;
@@ -169,6 +170,7 @@ void AuctionRunnerGS::runAuction(void)
             std::cout << "; error bound: " << numerator / denominator << std::endl;
 #endif
 #endif
+            relativeError = numerator / denominator;
             // if relative error is greater than delta, continue
             notDone = ( numerator / denominator > delta );
         }
@@ -243,6 +245,7 @@ double AuctionRunnerGS::getDistanceToQthPowerInternal(void)
         auto pA = bidders[bIdx];
         assert( 0 <= biddersToItems[bIdx] and biddersToItems[bIdx] < static_cast<int>(items.size()) );
         auto pB = items[biddersToItems[bIdx]];
+        std::cout << "pA = " << pA << ", pB = " << pB << ", pow(distLp(pA, pB, internal_p), wassersteinPower) = " << pow(distLp(pA, pB, internal_p), wassersteinPower) << ", dist = " << distLp(pA, pB, internal_p) << std::endl;
         result += pow(distLp(pA, pB, internal_p), wassersteinPower);
     }
     wassersteinCost = result;
diff --git a/geom_matching/wasserstein/src/auction_runner_jac.cpp b/geom_matching/wasserstein/src/auction_runner_jac.cpp
index b892643..c807ec9 100644
--- a/geom_matching/wasserstein/src/auction_runner_jac.cpp
+++ b/geom_matching/wasserstein/src/auction_runner_jac.cpp
@@ -188,7 +188,7 @@ void AuctionRunnerJac::flushAssignment(void)
 
 void AuctionRunnerJac::runAuction(void)
 {
-    // relative error 
+    relativeError = std::numeric_limits<double>::max();
     // choose some initial epsilon
     oracle->setEpsilon(oracle->maxVal / 4.0);
     assert( oracle->getEpsilon() > 0 );
@@ -210,6 +210,7 @@ void AuctionRunnerJac::runAuction(void)
         } else {
             denominator = pow(denominator, 1.0 / wassersteinPower);
             double numerator = currentResult - denominator;
+            relativeError = numerator / denominator;
             //std::cout << " numerator: " << numerator << " denominator: " << denominator << std::endl;
             //std::cout << " error bound: " << numerator / denominator << std::endl;
             // if relative error is greater than delta, continue
diff --git a/geom_matching/wasserstein/src/wasserstein.cpp b/geom_matching/wasserstein/src/wasserstein.cpp
index fc1b662..2c536ed 100644
--- a/geom_matching/wasserstein/src/wasserstein.cpp
+++ b/geom_matching/wasserstein/src/wasserstein.cpp
@@ -79,7 +79,8 @@ double wassersteinDistVec(const std::vector<DiagramPoint>& A,
 #else
     AuctionRunnerJac auction(A, B, q,  delta, _internal_p);
 #endif
-    return auction.getWassersteinDistance();
+    double result = auction.getWassersteinDistance();
+    return result;
 }
 
 double wassersteinCostVec(const std::vector<DiagramPoint>& A, 
@@ -119,7 +120,8 @@ double wassersteinCostVec(const std::vector<DiagramPoint>& A,
 #else
     AuctionRunnerJac auction(A, B, q,  delta, _internal_p);
 #endif
-    return auction.getWassersteinCost();
+    double result = auction.getWassersteinCost();
+    return result;
 }
 
 bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result)