From 7b850b8ee43fb7f8a0b2a1565ed01102d40b0a14 Mon Sep 17 00:00:00 2001 From: Arnur Nigmetov Date: Mon, 5 Sep 2016 13:32:05 +0200 Subject: Technical changes for R integration Avoid including iostream (R complains about that). All output protected by preprocessor directive (R checker should not see an instance of std::cout << in your code). Also added getWassersteinCost to be in line with the Dionysus implementation used in TDA. --- .../wasserstein/include/auction_runner_gs.h | 6 ++-- .../wasserstein/include/auction_runner_jac.h | 28 +++++++++-------- geom_matching/wasserstein/include/basic_defs_ws.h | 5 ++- geom_matching/wasserstein/include/def_debug_ws.h | 8 +++-- .../include/dnn/geometry/euclidean-fixed.h | 2 +- .../wasserstein/include/dnn/local/kd-tree.hpp | 9 ++---- .../wasserstein/include/dnn/parallel/tbb.h | 2 +- geom_matching/wasserstein/include/wasserstein.h | 36 ++++++++++++++++++++++ 8 files changed, 70 insertions(+), 26 deletions(-) (limited to 'geom_matching/wasserstein/include') diff --git a/geom_matching/wasserstein/include/auction_runner_gs.h b/geom_matching/wasserstein/include/auction_runner_gs.h index 34a91e8..80aa9f0 100644 --- a/geom_matching/wasserstein/include/auction_runner_gs.h +++ b/geom_matching/wasserstein/include/auction_runner_gs.h @@ -76,8 +76,9 @@ public: const double _initialEpsilon, const double _epsFactor); void setEpsilon(double newVal) { assert(epsilon > 0.0); epsilon = newVal; }; - double getEpsilon(void) const { return epsilon; } - double getWassersteinDistance(void); + double getEpsilon() const { return epsilon; } + double getWassersteinDistance(); + double getWassersteinCost(); static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling private: // private data @@ -94,6 +95,7 @@ private: double epsilonCommonRatio; // next epsilon = current epsilon / epsilonCommonRatio double weightAdjConst; double wassersteinDistance; + double wassersteinCost; // to get the 2 best items std::unique_ptr 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 ae0cb56..22d42b0 100644 --- a/geom_matching/wasserstein/include/auction_runner_jac.h +++ b/geom_matching/wasserstein/include/auction_runner_jac.h @@ -45,12 +45,13 @@ using AuctionOracle = AuctionOracleKDTreeRestricted; // 1. epsilonCommonRatio // 2. maxIterNum -class AuctionRunnerJak { +class AuctionRunnerJac { public: - AuctionRunnerJak(const std::vector& A, const std::vector& B, const double q, const double _delta, const double _internal_p); + AuctionRunnerJac(const std::vector& A, const std::vector& B, const double q, const double _delta, const double _internal_p); void setEpsilon(double newVal) { assert(epsilon > 0.0); epsilon = newVal; }; - double getEpsilon(void) const { return epsilon; } - double getWassersteinDistance(void); + double getEpsilon() const { return epsilon; } + double getWassersteinDistance(); + double getWassersteinCost(); static constexpr double epsilonCommonRatio { 5 }; // next epsilon = current epsilon / epsilonCommonRatio static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling private: @@ -66,6 +67,7 @@ private: double internal_p; double weightAdjConst; double wassersteinDistance; + double wassersteinCost; std::vector bidTable; // to get the 2 best items std::unique_ptr oracle; @@ -76,18 +78,18 @@ private: // private methods void assignGoodToBidder(const IdxType bidderIdx, const IdxType itemsIdx); void assignToBestBidder(const IdxType itemsIdx); - void clearBidTable(void); - void runAuction(void); - void runAuctionPhase(void); + void clearBidTable(); + void runAuction(); + void runAuctionPhase(); void submitBid(IdxType bidderIdx, const IdxValPair& itemsBidValuePair); - void flushAssignment(void); + void flushAssignment(); // for debug only - void sanityCheck(void); - void printDebug(void); - int countUnhappy(void); - void printMatching(void); - double getDistanceToQthPowerInternal(void); + void sanityCheck(); + void printDebug(); + int countUnhappy(); + void printMatching(); + double getDistanceToQthPowerInternal(); }; diff --git a/geom_matching/wasserstein/include/basic_defs_ws.h b/geom_matching/wasserstein/include/basic_defs_ws.h index 9e9c4ec..365c3bd 100644 --- a/geom_matching/wasserstein/include/basic_defs_ws.h +++ b/geom_matching/wasserstein/include/basic_defs_ws.h @@ -33,7 +33,6 @@ derivative works thereof, in binary and source code form. #include #include #include -#include #include #include @@ -58,7 +57,9 @@ struct Point { bool operator!=(const Point& other) const; Point(double ax, double ay) : x(ax), y(ay) {} Point() : x(0.0), y(0.0) {} +#ifndef FOR_R_TDA friend std::ostream& operator<<(std::ostream& output, const Point p); +#endif }; struct DiagramPoint @@ -76,7 +77,9 @@ struct DiagramPoint bool isNormal(void) const { return type == NORMAL; } double getRealX() const; // return the x-coord double getRealY() const; // return the y-coord +#ifndef FOR_R_TDA friend std::ostream& operator<<(std::ostream& output, const DiagramPoint p); +#endif struct LexicographicCmp { diff --git a/geom_matching/wasserstein/include/def_debug_ws.h b/geom_matching/wasserstein/include/def_debug_ws.h index 7323c18..6751556 100644 --- a/geom_matching/wasserstein/include/def_debug_ws.h +++ b/geom_matching/wasserstein/include/def_debug_ws.h @@ -25,12 +25,16 @@ derivative works thereof, in binary and source code form. */ -#ifndef DEF_DEBUG_H -#define DEF_DEBUG_H +#ifndef DEF_DEBUG_WS_H +#define DEF_DEBUG_WS_H //#define DEBUG_BOUND_MATCH //#define DEBUG_NEIGHBOUR_ORACLE //#define DEBUG_MATCHING //#define DEBUG_AUCTION +// This symbol should be defined only in the version +// for R package TDA, to comply with some CRAN rules +// like no usage of cout, cerr, cin, exit, etc. +//#define FOR_R_TDA #endif diff --git a/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h index a6ccef7..e2c5b44 100644 --- a/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h +++ b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h @@ -7,7 +7,7 @@ #include #include -#include +//#include #include #include #include diff --git a/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp b/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp index 151a4ad..6b0852c 100644 --- a/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp +++ b/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp @@ -6,6 +6,7 @@ #include #include "../parallel/tbb.h" +#include "def_debug_ws.h" template dnn::KDTree:: @@ -127,13 +128,8 @@ search(PointHandle q, ResultsFunctor& rf) const // TODO: use tbb::scalable_allocator for the queue std::queue nodes; - - nodes.push(KDTreeNode(tree_.begin(), tree_.end(), 0)); - - //std::cout << "started kdtree::search" << std::endl; - while (!nodes.empty()) { HCIterator b, e; size_t i; @@ -163,7 +159,6 @@ search(PointHandle q, ResultsFunctor& rf) const nodes.push(KDTreeNode(b, m, i)); } } - //std::cout << "exited kdtree::search" << std::endl; } template @@ -290,6 +285,7 @@ void dnn::KDTree:: printWeights(void) { +#ifndef FOR_R_TDA std::cout << "weights_:" << std::endl; for(const auto ph : indices_) { std::cout << "idx = " << ph.second << ": (" << (ph.first)->at(0) << ", " << (ph.first)->at(1) << ") weight = " << weights_[ph.second] << std::endl; @@ -298,6 +294,7 @@ printWeights(void) for(size_t idx = 0; idx < subtree_weights_.size(); ++idx) { std::cout << idx << " : " << subtree_weights_[idx] << std::endl; } +#endif } diff --git a/geom_matching/wasserstein/include/dnn/parallel/tbb.h b/geom_matching/wasserstein/include/dnn/parallel/tbb.h index 4aa6805..64c59e0 100644 --- a/geom_matching/wasserstein/include/dnn/parallel/tbb.h +++ b/geom_matching/wasserstein/include/dnn/parallel/tbb.h @@ -1,7 +1,7 @@ #ifndef PARALLEL_H #define PARALLEL_H -#include +//#include #include #include diff --git a/geom_matching/wasserstein/include/wasserstein.h b/geom_matching/wasserstein/include/wasserstein.h index 38ac6bd..fcc9e74 100644 --- a/geom_matching/wasserstein/include/wasserstein.h +++ b/geom_matching/wasserstein/include/wasserstein.h @@ -49,6 +49,15 @@ double wassersteinDistVec(const std::vector& A, const double _initialEpsilon = 0.0, const double _epsFactor = 0.0); +// get Wasserstein cost (distance^q) between two persistence diagrams +double wassersteinCostVec(const std::vector& A, + const std::vector& B, + const double q, + const double delta, + const double _internal_p = std::numeric_limits::infinity(), + const double _initialEpsilon = 0.0, + const double _epsFactor = 0.0); + // compare as multisets template @@ -98,6 +107,33 @@ double wassersteinDist(PairContainer& A, PairContainer& B, const double q, const return wassersteinDistVec(dgmA, dgmB, q, delta, _internal_p, _initialEpsilon, _epsFactor); } +template +double wassersteinCost(PairContainer& A, PairContainer& B, const double q, const double delta, const double _internal_p = std::numeric_limits::infinity(), const double _initialEpsilon = 0.0, const double _epsFactor = 0.0) +{ + if (areEqual(A, B)) { + return 0.0; + } + + std::vector dgmA, dgmB; + // loop over A, add projections of A-points to corresponding positions + // in B-vector + for(auto& pairA : A) { + double x = pairA.first; + double y = pairA.second; + dgmA.push_back(DiagramPoint(x, y, DiagramPoint::NORMAL)); + dgmB.push_back(DiagramPoint(x, y, DiagramPoint::DIAG)); + } + // the same for B + for(auto& pairB : B) { + double x = pairB.first; + double y = pairB.second; + dgmA.push_back(DiagramPoint(x, y, DiagramPoint::DIAG)); + dgmB.push_back(DiagramPoint(x, y, DiagramPoint::NORMAL)); + } + + return wassersteinCostVec(dgmA, dgmB, q, delta, _internal_p, _initialEpsilon, _epsFactor); +} + // fill in result with points from file fname // return false if file can't be opened -- cgit v1.2.3