From 8ffc1482c88fd1c503524f6e328f69fe3fdede6c Mon Sep 17 00:00:00 2001 From: Arnur Nigmetov Date: Thu, 31 May 2018 20:44:59 +0200 Subject: Fix bug in point cloud version. internal_p was not respected by DynamicPointTraits, l_2 distance was always used. --- geom_matching/wasserstein/CMakeLists.txt | 8 ++++---- .../include/auction_oracle_kdtree_pure_geom.hpp | 17 ++++++++------- .../include/auction_oracle_kdtree_restricted.hpp | 4 ++-- .../wasserstein/include/auction_runner_gs.hpp | 16 ++++++++++----- geom_matching/wasserstein/include/basic_defs_ws.h | 23 +++++++++++---------- .../include/dnn/geometry/euclidean-dynamic.h | 24 +++++++++++++++++++++- geom_matching/wasserstein/include/hera_infinity.h | 22 ++++++++++++++++++++ 7 files changed, 84 insertions(+), 30 deletions(-) create mode 100644 geom_matching/wasserstein/include/hera_infinity.h 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::get_optimal_bid_debug(IdxTy Real best_item_value = std::numeric_limits::max(); Real second_best_item_value = std::numeric_limits::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::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 AuctionOracleKDTreePureGeom::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::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::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::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::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::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 R AuctionRunnerGS::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::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::sanity_check() } for(IdxType item_idx = 0; item_idx < static_cast(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(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 -inline bool is_infinity(const Real& x) -{ - return x == Real(-1); -}; - -template -inline Real get_infinity() -{ - return Real( -1 ); -} +//template +//inline bool is_infinity(const Real& x) +//{ +// return x == Real(-1); +//}; +// +//template +//inline Real get_infinity() +//{ +// return Real( -1 ); +//} template 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 #include +#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 + inline bool is_infinity(const Real& x) + { + return x == Real(-1); + }; + + template + inline Real get_infinity() + { + return Real(-1); + } +} + +#endif //WASSERSTEIN_HERA_INFINITY_H -- cgit v1.2.3