Fix bug in point cloud version.

internal_p was not respected by DynamicPointTraits, l_2 distance was always used.

7 files changed, 84 insertions, 30 deletions

diff --git a/geom_matching/wasserstein/CMakeLists.txt b/geom_matching/wasserstein/CMakeLists.txt
index b59378c..c6fba2c 100644
--- a/geom_matching/wasserstein/CMakeLists.txt
+++ b/geom_matching/wasserstein/CMakeLists.txt
@@ -46,17 +46,17 @@ file(GLOB WS_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/include/*.h ${CMAKE_CURRENT_SOU
 find_package (Threads)
 set (libraries ${libraries}     ${CMAKE_THREAD_LIBS_INIT})
 
-add_executable(wasserstein_dist ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist.cpp ${WS_HEADERS})
+add_executable(wasserstein_dist ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist.cpp ${WS_HEADERS} include/hera_infinity.h)
 target_link_libraries(wasserstein_dist PUBLIC ${libraries})
 
-add_executable(wasserstein_dist_dipha ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist_dipha.cpp ${WS_HEADERS})
+add_executable(wasserstein_dist_dipha ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist_dipha.cpp ${WS_HEADERS} include/hera_infinity.h)
 target_link_libraries(wasserstein_dist_dipha PUBLIC ${libraries})
 
 # pure geometric version, arbitrary dimension
-add_executable(wasserstein_dist_point_cloud ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist_point_cloud.cpp ${WS_HEADERS})
+add_executable(wasserstein_dist_point_cloud ${CMAKE_CURRENT_SOURCE_DIR}/example/wasserstein_dist_point_cloud.cpp ${WS_HEADERS} include/hera_infinity.h)
 target_link_libraries(wasserstein_dist_point_cloud PUBLIC ${libraries})
 
 # Tests
-add_executable(wasserstein_test  ${CMAKE_CURRENT_SOURCE_DIR}/tests/tests_main.cpp ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_hera_wasserstein.cpp)
+add_executable(wasserstein_test  ${CMAKE_CURRENT_SOURCE_DIR}/tests/tests_main.cpp ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_hera_wasserstein.cpp include/hera_infinity.h)
 #add_executable(wasserstein_test EXCLUDE_FROM_ALL ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_hera_wasserstein.cpp)
 target_link_libraries(wasserstein_test PUBLIC ${libraries})
diff --git a/geom_matching/wasserstein/include/auction_oracle_kdtree_pure_geom.hpp b/geom_matching/wasserstein/include/auction_oracle_kdtree_pure_geom.hpp
index a6bdf10..eaf54cf 100644
--- a/geom_matching/wasserstein/include/auction_oracle_kdtree_pure_geom.hpp
+++ b/geom_matching/wasserstein/include/auction_oracle_kdtree_pure_geom.hpp
@@ -111,11 +111,9 @@ AuctionOracleKDTreePureGeom<Real_, PointContainer_>::get_optimal_bid_debug(IdxTy
     Real best_item_value = std::numeric_limits<Real>::max();
     Real second_best_item_value = std::numeric_limits<Real>::max();
 
-    for(IdxType item_idx = 0; item_idx < this->items.size(); ++item_idx) {
+    for(size_t item_idx = 0; item_idx < this->items.size(); ++item_idx) {
         auto item = this->items[item_idx];
-        if (item.type != bidder.type and item_idx != bidder_idx)
-            continue;
-        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p), this->wasserstein_power, this->dim) + this->prices[item_idx];
+        auto item_value = std::pow(traits.distance(bidder, item), this->wasserstein_power) + this->prices[item_idx];
         if (item_value < best_item_value) {
             best_item_value = item_value;
             best_item_idx = item_idx;
@@ -126,11 +124,10 @@ AuctionOracleKDTreePureGeom<Real_, PointContainer_>::get_optimal_bid_debug(IdxTy
 
     for(size_t item_idx = 0; item_idx < this->items.size(); ++item_idx) {
         auto item = this->items[item_idx];
-        if (item.type != bidder.type and item_idx != bidder_idx)
-            continue;
         if (item_idx == best_item_idx)
             continue;
-        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p), this->wasserstein_power, this->dim) + this->prices[item_idx];
+
+        auto item_value = std::pow(traits.distance(bidder, item), this->wasserstein_power) + this->prices[item_idx];
         if (item_value < second_best_item_value) {
             second_best_item_value = item_value;
             second_best_item_idx = item_idx;
@@ -166,6 +163,12 @@ IdxValPair<Real_> AuctionOracleKDTreePureGeom<Real_, PointContainer_>::get_optim
     result.first = best_item_idx;
     result.second = ( second_best_item_value - best_item_value ) + this->prices[best_item_idx] + this->epsilon;
 
+#ifdef DEBUG_KDTREE_RESTR_ORACLE
+    auto bid_debug = get_optimal_bid_debug(bidder_idx);
+    assert(fabs(bid_debug.best_item_value - best_item_value) < 0.000000001);
+    assert(fabs(bid_debug.second_best_item_value - second_best_item_value) < 0.000000001);
+#endif
+
     return result;
 }
 
diff --git a/geom_matching/wasserstein/include/auction_oracle_kdtree_restricted.hpp b/geom_matching/wasserstein/include/auction_oracle_kdtree_restricted.hpp
index 0e6f780..7817bf3 100644
--- a/geom_matching/wasserstein/include/auction_oracle_kdtree_restricted.hpp
+++ b/geom_matching/wasserstein/include/auction_oracle_kdtree_restricted.hpp
@@ -247,7 +247,7 @@ AuctionOracleKDTreeRestricted<Real_, PointContainer_>::get_optimal_bid_debug(Idx
         auto item = this->items[item_idx];
         if (item.type != bidder.type and item_idx != bidder_idx)
             continue;
-        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p), this->wasserstein_power) + this->prices[item_idx];
+        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p, 2), this->wasserstein_power) + this->prices[item_idx];
         if (item_value < best_item_value) {
             best_item_value = item_value;
             best_item_idx = item_idx;
@@ -262,7 +262,7 @@ AuctionOracleKDTreeRestricted<Real_, PointContainer_>::get_optimal_bid_debug(Idx
             continue;
         if (item_idx == best_item_idx)
             continue;
-        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p), this->wasserstein_power) + this->prices[item_idx];
+        auto item_value = std::pow(dist_lp(bidder, item, this->internal_p, 2), this->wasserstein_power) + this->prices[item_idx];
         if (item_value < second_best_item_value) {
             second_best_item_value = item_value;
             second_best_item_idx = item_idx;
diff --git a/geom_matching/wasserstein/include/auction_runner_gs.hpp b/geom_matching/wasserstein/include/auction_runner_gs.hpp
index d9f419d..960c707 100644
--- a/geom_matching/wasserstein/include/auction_runner_gs.hpp
+++ b/geom_matching/wasserstein/include/auction_runner_gs.hpp
@@ -244,6 +244,12 @@ void AuctionRunnerGS<R, AO, PC>::run_auction_phases(const int max_num_phases, co
         flush_assignment();
         run_auction_phase();
         Real current_result = getDistanceToQthPowerInternal();
+//        Real current_result_1 = 0.0;
+//        for(size_t i = 0; i < num_bidders; ++i) {
+//            current_result_1 += oracle.traits.distance(bidders[i], items[bidders_to_items[i]]);
+//        }
+//        current_result = current_result_1;
+//        assert(fabs(current_result - current_result_1) < 0.001);
         Real denominator = current_result - num_bidders * oracle.get_epsilon();
         current_result = pow(current_result, 1.0 / wasserstein_power);
 #ifdef LOG_AUCTION
@@ -259,6 +265,7 @@ void AuctionRunnerGS<R, AO, PC>::run_auction_phases(const int max_num_phases, co
             denominator = pow(denominator, 1.0 / wasserstein_power);
             Real numerator = current_result - denominator;
             relative_error = numerator / denominator;
+            // spdlog::get("console")->info("relative error = {} / {} = {}, result = {}", numerator, denominator, relative_error, current_result);
 #ifdef LOG_AUCTION
             console_logger->info("error = {0} / {1} = {2}",
                     numerator, denominator, relative_error);
@@ -319,7 +326,7 @@ void AuctionRunnerGS<R, AO, PC>::run_auction_phase()
 
 #ifdef DEBUG_AUCTION
     for(size_t bidder_idx = 0; bidder_idx < num_bidders; ++bidder_idx) {
-        if ( bidders_to_items[bidder_idx] < 0 or bidders_to_items[bidder_idx] >= num_bidders) {
+        if ( bidders_to_items[bidder_idx] < 0 or bidders_to_items[bidder_idx] >= (IdxType)num_bidders) {
             std::cerr << "After auction terminated bidder " << bidder_idx;
             std::cerr << " has no items assigned" << std::endl;
             throw std::runtime_error("Auction did not give a perfect matching");
@@ -333,8 +340,7 @@ template<class R, class AO, class PC>
 R AuctionRunnerGS<R, AO, PC>::get_item_bidder_cost(const size_t item_idx, const size_t bidder_idx, const bool tolerate_invalid_idx) const
 {
     if (item_idx != k_invalid_index and bidder_idx != k_invalid_index) {
-        return std::pow(dist_lp(bidders[bidder_idx], items[item_idx], internal_p, dimension),
-                        wasserstein_power);
+        return std::pow(dist_lp(bidders[bidder_idx], items[item_idx], internal_p, dimension), wasserstein_power);
     } else {
         if (tolerate_invalid_idx)
             return R(0.0);
@@ -416,7 +422,7 @@ void AuctionRunnerGS<R, AO, PC>::sanity_check()
     }
 
     for(size_t bidder_idx = 0; bidder_idx < num_bidders; ++bidder_idx) {
-        assert( bidders_to_items[bidder_idx] == k_invalid_index or ( bidders_to_items[bidder_idx] < num_items and bidders_to_items[bidder_idx] >= 0));
+        assert( bidders_to_items[bidder_idx] == k_invalid_index or ( bidders_to_items[bidder_idx] < (IdxType)num_items and bidders_to_items[bidder_idx] >= 0));
 
         if ( bidders_to_items[bidder_idx] != k_invalid_index) {
 
@@ -440,7 +446,7 @@ void AuctionRunnerGS<R, AO, PC>::sanity_check()
     }
 
     for(IdxType item_idx = 0; item_idx < static_cast<IdxType>(num_bidders); ++item_idx) {
-        assert( items_to_bidders[item_idx] == k_invalid_index or ( items_to_bidders[item_idx] < num_items and items_to_bidders[item_idx] >= 0));
+        assert( items_to_bidders[item_idx] == k_invalid_index or ( items_to_bidders[item_idx] < static_cast<IdxType>(num_items) and items_to_bidders[item_idx] >= 0));
         if ( items_to_bidders.at(item_idx) != k_invalid_index) {
 
             // check for uniqueness
diff --git a/geom_matching/wasserstein/include/basic_defs_ws.h b/geom_matching/wasserstein/include/basic_defs_ws.h
index 28f7452..1c5928f 100644
--- a/geom_matching/wasserstein/include/basic_defs_ws.h
+++ b/geom_matching/wasserstein/include/basic_defs_ws.h
@@ -51,6 +51,7 @@ derivative works thereof, in binary and source code form.
 #include "spdlog/fmt/ostr.h"
 #endif
 
+#include "hera_infinity.h"
 #include "dnn/geometry/euclidean-dynamic.h"
 #include "def_debug_ws.h"
 
@@ -59,17 +60,17 @@ derivative works thereof, in binary and source code form.
 namespace hera
 {
 
-template<class Real = double>
-inline bool is_infinity(const Real& x)
-{
-    return x == Real(-1);
-};
-
-template<class Real = double>
-inline Real get_infinity()
-{
-    return Real( -1 );
-}
+//template<class Real = double>
+//inline bool is_infinity(const Real& x)
+//{
+//    return x == Real(-1);
+//};
+//
+//template<class Real = double>
+//inline Real get_infinity()
+//{
+//    return Real( -1 );
+//}
 
 template<class Real = double>
 inline bool is_p_valid_norm(const Real& p)
diff --git a/geom_matching/wasserstein/include/dnn/geometry/euclidean-dynamic.h b/geom_matching/wasserstein/include/dnn/geometry/euclidean-dynamic.h
index 4b98309..b003906 100644
--- a/geom_matching/wasserstein/include/dnn/geometry/euclidean-dynamic.h
+++ b/geom_matching/wasserstein/include/dnn/geometry/euclidean-dynamic.h
@@ -8,6 +8,8 @@
 #include <boost/serialization/vector.hpp>
 #include <cmath>
 
+#include "hera_infinity.h"
+
 namespace hera
 {
 namespace ws
@@ -89,7 +91,27 @@ struct DynamicPointTraits
                     DynamicPointTraits(unsigned dim = 0):
                         dim_(dim)                                       {}
 
-    DistanceType    distance(PointType p1, PointType p2) const          { return sqrt(sq_distance(p1,p2)); }
+    DistanceType    distance(PointType p1, PointType p2) const
+        {
+            Real result = 0.0;
+            if (hera::is_infinity(internal_p)) {
+                // max norm
+                for (unsigned i = 0; i < dimension(); ++i)
+                    result = std::max(result, fabs(coordinate(p1,i) - coordinate(p2,i)));
+            } else if (internal_p == Real(1.0)) {
+                // l1-norm
+                for (unsigned i = 0; i < dimension(); ++i)
+                    result += fabs(coordinate(p1,i) - coordinate(p2,i));
+            } else if (internal_p == Real(2.0)) {
+                result = sqrt(sq_distance(p1,p2));
+            } else {
+                assert(internal_p > 1.0);
+                for (unsigned i = 0; i < dimension(); ++i)
+                    result += std::pow(fabs(coordinate(p1,i) - coordinate(p2,i)), internal_p);
+                result = std::pow(result, Real(1.0) / internal_p);
+            }
+            return result;
+        }
     DistanceType    distance(PointHandle p1, PointHandle p2) const      { return distance(PointType({p1.p}), PointType({p2.p})); }
     DistanceType    sq_distance(PointType p1, PointType p2) const       { Real res = 0; for (unsigned i = 0; i < dimension(); ++i) { Real c1 = coordinate(p1,i), c2 = coordinate(p2,i); res += (c1 - c2)*(c1 - c2); } return res; }
     DistanceType    sq_distance(PointHandle p1, PointHandle p2) const   { return sq_distance(PointType({p1.p}), PointType({p2.p})); }
diff --git a/geom_matching/wasserstein/include/hera_infinity.h b/geom_matching/wasserstein/include/hera_infinity.h
new file mode 100644
index 0000000..5a446e7
--- /dev/null
+++ b/geom_matching/wasserstein/include/hera_infinity.h
@@ -0,0 +1,22 @@
+#ifndef WASSERSTEIN_HERA_INFINITY_H
+#define WASSERSTEIN_HERA_INFINITY_H
+
+// we cannot assume that template parameter Real will always provide infinity() value,
+// so value -1.0 is used to encode infinity (l_inf norm is used by default)
+
+namespace hera {
+
+    template<class Real = double>
+    inline bool is_infinity(const Real& x)
+    {
+        return x == Real(-1);
+    };
+
+    template<class Real = double>
+    inline Real get_infinity()
+    {
+        return Real(-1);
+    }
+}
+
+#endif //WASSERSTEIN_HERA_INFINITY_H