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. --- geom_bottleneck/bottleneck/include/ANN/ANN.h | 13 ++++- geom_bottleneck/bottleneck/include/ANN/bd_tree.h | 3 + geom_bottleneck/bottleneck/include/ANN/kd_tree.h | 7 +++ geom_bottleneck/bottleneck/include/basic_defs_bt.h | 12 +++- geom_bottleneck/bottleneck/include/bottleneck.h | 2 +- geom_bottleneck/bottleneck/include/def_debug_bt.h | 8 ++- geom_bottleneck/bottleneck/src/ann/ANN.cpp | 11 +++- geom_bottleneck/bottleneck/src/ann/bd_tree.cpp | 3 + geom_bottleneck/bottleneck/src/ann/kd_dump.cpp | 11 ++++ geom_bottleneck/bottleneck/src/ann/kd_tree.cpp | 8 ++- geom_bottleneck/bottleneck/src/basic_defs.cpp | 13 ++--- geom_bottleneck/bottleneck/src/bottleneck.cpp | 38 ++++++------- geom_bottleneck/bottleneck/src/bound_match.cpp | 22 +++++++- geom_bottleneck/bottleneck/src/neighb_oracle.cpp | 2 + .../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 ++++++++++++ geom_matching/wasserstein/src/auction_oracle.cpp | 57 ++++--------------- .../wasserstein/src/auction_runner_gs.cpp | 61 +++++++++++++++++--- .../wasserstein/src/auction_runner_jac.cpp | 52 ++++++++++++----- geom_matching/wasserstein/src/basic_defs.cpp | 10 +++- geom_matching/wasserstein/src/wasserstein.cpp | 65 ++++++++++++++++++++-- 27 files changed, 354 insertions(+), 140 deletions(-) diff --git a/geom_bottleneck/bottleneck/include/ANN/ANN.h b/geom_bottleneck/bottleneck/include/ANN/ANN.h index cd48d8e..004dfe2 100644 --- a/geom_bottleneck/bottleneck/include/ANN/ANN.h +++ b/geom_bottleneck/bottleneck/include/ANN/ANN.h @@ -96,7 +96,6 @@ #include // standard lib includes #include // math includes -#include // I/O streams #include // C-style strings #include #include @@ -130,6 +129,12 @@ #define ANNcopyright "David M. Mount and Sunil Arya" #define ANNlatestRev "Jan 27, 2010" +#include "def_debug_bt.h" + +#ifndef FOR_R_TDA +#include // I/O streams +#endif + namespace geom_bt { //---------------------------------------------------------------------- // ANNbool @@ -798,8 +803,10 @@ public: int bs = 1, // bucket size ANNsplitRule split = ANN_KD_SUGGEST); // splitting method +#ifndef FOR_R_TDA ANNkd_tree( // build from dump file std::istream& in); // input stream for dump file +#endif ~ANNkd_tree(); // tree destructor @@ -834,6 +841,7 @@ public: ANNpointArray thePoints() // return pointer to points { return pts; } +#ifndef FOR_R_TDA virtual void Print( // print the tree (for debugging) ANNbool with_pts, // print points as well? std::ostream& out); // output stream @@ -841,6 +849,7 @@ public: virtual void Dump( // dump entire tree ANNbool with_pts, // print points as well? std::ostream& out); // output stream +#endif virtual void getStats( // compute tree statistics ANNkdStats& st); // the statistics (modified) @@ -885,8 +894,10 @@ public: ANNsplitRule split = ANN_KD_SUGGEST, // splitting rule ANNshrinkRule shrink = ANN_BD_SUGGEST); // shrinking rule +#ifndef FOR_R_TDA ANNbd_tree( // build from dump file std::istream& in); // input stream for dump file +#endif }; //---------------------------------------------------------------------- diff --git a/geom_bottleneck/bottleneck/include/ANN/bd_tree.h b/geom_bottleneck/bottleneck/include/ANN/bd_tree.h index 0791429..38cecb7 100644 --- a/geom_bottleneck/bottleneck/include/ANN/bd_tree.h +++ b/geom_bottleneck/bottleneck/include/ANN/bd_tree.h @@ -29,6 +29,7 @@ #include // all ANN includes #include "kd_tree.h" // kd-tree includes +#include "def_debug_bt.h" namespace geom_bt { //---------------------------------------------------------------------- @@ -91,7 +92,9 @@ public: ANNkdStats &st, // statistics ANNorthRect &bnd_box); // bounding box virtual void print(int level, ostream &out);// print node +#ifndef FOR_R_TDA virtual void dump(ostream &out); // dump node +#endif virtual void ann_search(ANNdist); // standard search virtual void ann_pri_search(ANNdist); // priority search diff --git a/geom_bottleneck/bottleneck/include/ANN/kd_tree.h b/geom_bottleneck/bottleneck/include/ANN/kd_tree.h index 5fb362d..a1e53e5 100644 --- a/geom_bottleneck/bottleneck/include/ANN/kd_tree.h +++ b/geom_bottleneck/bottleneck/include/ANN/kd_tree.h @@ -31,6 +31,7 @@ #include // for std::pair #include // all ANN includes +#include "def_debug_bt.h" using namespace std; // make std:: available @@ -73,7 +74,9 @@ public: ANNorthRect &bnd_box) = 0; // bounding box // print node virtual void print(int level, ostream &out) = 0; +#ifndef FOR_R_TDA virtual void dump(ostream &out) = 0; // dump node +#endif friend class ANNkd_tree; // allow kd-tree to access us @@ -139,7 +142,9 @@ public: ANNkdStats &st, // statistics ANNorthRect &bnd_box); // bounding box virtual void print(int level, ostream &out);// print node +#ifndef FOR_R_TDA virtual void dump(ostream &out); // dump node +#endif virtual void ann_search(ANNdist); // standard search virtual void ann_pri_search(ANNdist); // priority search @@ -217,7 +222,9 @@ public: ANNkdStats &st, // statistics ANNorthRect &bnd_box); // bounding box virtual void print(int level, ostream &out);// print node +#ifndef FOR_R_TDA virtual void dump(ostream &out); // dump node +#endif virtual void ann_search(ANNdist); // standard search virtual void ann_pri_search(ANNdist); // priority search diff --git a/geom_bottleneck/bottleneck/include/basic_defs_bt.h b/geom_bottleneck/bottleneck/include/basic_defs_bt.h index 5334b6d..759f18e 100644 --- a/geom_bottleneck/bottleneck/include/basic_defs_bt.h +++ b/geom_bottleneck/bottleneck/include/basic_defs_bt.h @@ -27,16 +27,20 @@ #endif #include +#include #include #include #include #include -#include #include #include #include "def_debug_bt.h" +#ifndef FOR_R_TDA +#include +#endif + namespace geom_bt { @@ -50,7 +54,9 @@ struct Point { bool operator!=(const Point& other) const; Point(CoordinateType ax, CoordinateType 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 @@ -90,7 +96,9 @@ public: //return 0.5 * ( x + y); } +#ifndef FOR_R_TDA friend std::ostream& operator<<(std::ostream& output, const DiagramPoint p); +#endif }; struct PointHash { @@ -130,7 +138,9 @@ public: std::unordered_set::iterator end() { return points.end(); } std::unordered_set::const_iterator cbegin() const { return points.cbegin(); } std::unordered_set::const_iterator cend() const { return points.cend(); } +#ifndef FOR_R_TDA friend std::ostream& operator<<(std::ostream& output, const DiagramPointSet& ps); +#endif friend void addProjections(DiagramPointSet& A, DiagramPointSet& B); template DiagramPointSet(PairIterator first, PairIterator last); template void fillIn(PairIterator first, PairIterator last); diff --git a/geom_bottleneck/bottleneck/include/bottleneck.h b/geom_bottleneck/bottleneck/include/bottleneck.h index 19ae89a..3267c5d 100644 --- a/geom_bottleneck/bottleneck/include/bottleneck.h +++ b/geom_bottleneck/bottleneck/include/bottleneck.h @@ -22,7 +22,7 @@ #define BOTTLENECK_H -#include +//#include #include #include #include diff --git a/geom_bottleneck/bottleneck/include/def_debug_bt.h b/geom_bottleneck/bottleneck/include/def_debug_bt.h index eaf356d..60f834e 100644 --- a/geom_bottleneck/bottleneck/include/def_debug_bt.h +++ b/geom_bottleneck/bottleneck/include/def_debug_bt.h @@ -18,12 +18,16 @@ */ -#ifndef DEF_DEBUG_H -#define DEF_DEBUG_H +#ifndef DEF_DEBUG_BT_H +#define DEF_DEBUG_BT_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_bottleneck/bottleneck/src/ann/ANN.cpp b/geom_bottleneck/bottleneck/src/ann/ANN.cpp index 7bae577..83c7ef6 100644 --- a/geom_bottleneck/bottleneck/src/ann/ANN.cpp +++ b/geom_bottleneck/bottleneck/src/ann/ANN.cpp @@ -30,9 +30,11 @@ #include // make VS more conformal #endif +#include #include // C standard lib defs #include // all ANN includes #include // ANN performance +#include "def_debug_bt.h" @@ -80,10 +82,12 @@ void annPrintPt( // print a point int dim, // the dimension std::ostream &out) // output stream { +#ifndef FOR_R_TDA for (int j = 0; j < dim; j++) { out << pt[j]; if (j < dim-1) out << " "; } +#endif } //---------------------------------------------------------------------- @@ -197,11 +201,16 @@ bool ANNorthRect::intersects(const int dim, const ANNorthRect& r) const void annError(const char* msg, ANNerr level) { if (level == ANNabort) { +#ifndef FOR_R_TDA cerr << "ANN: ERROR------->" << msg << "<-------------ERROR\n"; - exit(1); +#endif + throw std::runtime_error(std::string("ANN: Error: ") + std::string(msg)); + //exit(1); } else { +#ifndef FOR_R_TDA cerr << "ANN: WARNING----->" << msg << "<-------------WARNING\n"; +#endif } } diff --git a/geom_bottleneck/bottleneck/src/ann/bd_tree.cpp b/geom_bottleneck/bottleneck/src/ann/bd_tree.cpp index 8c1ef6d..a5dd69c 100644 --- a/geom_bottleneck/bottleneck/src/ann/bd_tree.cpp +++ b/geom_bottleneck/bottleneck/src/ann/bd_tree.cpp @@ -31,6 +31,7 @@ #include "kd_split.h" // kd-tree splitting rules #include // performance evaluation +#include "def_debug_bt.h" namespace geom_bt { //---------------------------------------------------------------------- @@ -43,6 +44,7 @@ void ANNbd_shrink::print( // print shrinking node int level, // depth of node in tree ostream &out) // output stream { +#ifndef FOR_R_TDA child[ANN_OUT]->print(level+1, out); // print out-child out << " "; @@ -61,6 +63,7 @@ void ANNbd_shrink::print( // print shrinking node out << "\n"; child[ANN_IN]->print(level+1, out); // print in-child +#endif } //---------------------------------------------------------------------- diff --git a/geom_bottleneck/bottleneck/src/ann/kd_dump.cpp b/geom_bottleneck/bottleneck/src/ann/kd_dump.cpp index 64db9a7..ecaf7ea 100644 --- a/geom_bottleneck/bottleneck/src/ann/kd_dump.cpp +++ b/geom_bottleneck/bottleneck/src/ann/kd_dump.cpp @@ -33,6 +33,7 @@ #include "kd_tree.h" // kd-tree declarations #include "bd_tree.h" // bd-tree declarations +#include "def_debug_bt.h" using namespace std; // make std:: available @@ -101,6 +102,7 @@ namespace geom_bt { // ... (repeated n_bnds times) //---------------------------------------------------------------------- +#ifndef FOR_R_TDA void ANNkd_tree::Dump( // dump entire tree ANNbool with_pts, // print points as well? ostream &out) // output stream @@ -132,20 +134,24 @@ namespace geom_bt { } out.precision(0); // restore default precision } +#endif void ANNkd_split::dump( // dump a splitting node ostream &out) // output stream { +#ifndef FOR_R_TDA out << "split " << cut_dim << " " << cut_val << " "; out << cd_bnds[ANN_LO] << " " << cd_bnds[ANN_HI] << "\n"; child[ANN_LO]->dump(out); // print low child child[ANN_HI]->dump(out); // print high child +#endif } void ANNkd_leaf::dump( // dump a leaf node ostream &out) // output stream { +#ifndef FOR_R_TDA if (this == KD_TRIVIAL) { // canonical trivial leaf node out << "leaf 0\n"; // leaf no points } @@ -156,17 +162,20 @@ namespace geom_bt { } out << "\n"; } +#endif } void ANNbd_shrink::dump( // dump a shrinking node ostream &out) // output stream { +#ifndef FOR_R_TDA out << "shrink " << n_bnds << "\n"; for (int j = 0; j < n_bnds; j++) { out << bnds[j].cd << " " << bnds[j].cv << " " << bnds[j].sd << "\n"; } child[ANN_IN]->dump(out); // print in-child child[ANN_OUT]->dump(out); // print out-child +#endif } //---------------------------------------------------------------------- @@ -441,7 +450,9 @@ namespace geom_bt { } else { annError("Illegal node type in dump file", ANNabort); +#ifndef FOR_R_TDA exit(0); // to keep the compiler happy +#endif } } } diff --git a/geom_bottleneck/bottleneck/src/ann/kd_tree.cpp b/geom_bottleneck/bottleneck/src/ann/kd_tree.cpp index ad3a82d..e8f7f63 100644 --- a/geom_bottleneck/bottleneck/src/ann/kd_tree.cpp +++ b/geom_bottleneck/bottleneck/src/ann/kd_tree.cpp @@ -38,6 +38,7 @@ #include "kd_split.h" // kd-tree splitting rules #include "kd_util.h" // kd-tree utilities #include // performance evaluation +#include "def_debug_bt.h" namespace geom_bt { //---------------------------------------------------------------------- @@ -75,6 +76,7 @@ void ANNkd_split::print( // print splitting node int level, // depth of node in tree ostream &out) // output stream { +#ifndef FOR_R_TDA child[ANN_HI]->print(level+1, out); // print high child out << " "; for (int i = 0; i < level; i++) // print indentation @@ -85,13 +87,14 @@ void ANNkd_split::print( // print splitting node out << " np=" << actual_num_points; out << "\n"; child[ANN_LO]->print(level+1, out); // print low child +#endif } void ANNkd_leaf::print( // print leaf node int level, // depth of node in tree ostream &out) // output stream { - +#ifndef FOR_R_TDA out << " "; for (int i = 0; i < level; i++) // print indentation out << ".."; @@ -107,8 +110,10 @@ void ANNkd_leaf::print( // print leaf node } out << ">\n"; } +#endif } +#ifndef FOR_R_TDA void ANNkd_tree::Print( // print entire tree ANNbool with_pts, // print points as well? ostream &out) // output stream @@ -128,6 +133,7 @@ void ANNkd_tree::Print( // print entire tree root->print(0, out); // invoke printing at root } } +#endif //---------------------------------------------------------------------- // kd_tree statistics (for performance evaluation) diff --git a/geom_bottleneck/bottleneck/src/basic_defs.cpp b/geom_bottleneck/bottleneck/src/basic_defs.cpp index e09b119..76e6cc5 100644 --- a/geom_bottleneck/bottleneck/src/basic_defs.cpp +++ b/geom_bottleneck/bottleneck/src/basic_defs.cpp @@ -20,6 +20,7 @@ #include #include +#include "def_debug_bt.h" #include "basic_defs_bt.h" namespace geom_bt { @@ -36,6 +37,7 @@ bool Point::operator!=(const Point& other) const return !(*this == other); } +#ifndef FOR_R_TDA std::ostream& operator<<(std::ostream& output, const Point p) { output << "(" << p.x << ", " << p.y << ")"; @@ -51,6 +53,7 @@ std::ostream& operator<<(std::ostream& output, const PointSet& ps) output << "\b\b }"; return output; } +#endif double sqrDist(const Point& a, const Point& b) { @@ -90,6 +93,7 @@ bool DiagramPoint::operator!=(const DiagramPoint& other) const return !(*this == other); } +#ifndef FOR_R_TDA std::ostream& operator<<(std::ostream& output, const DiagramPoint p) { if ( p.type == DiagramPoint::DIAG ) { @@ -109,6 +113,7 @@ std::ostream& operator<<(std::ostream& output, const DiagramPointSet& ps) output << "\b\b }"; return output; } +#endif DiagramPoint::DiagramPoint(double xx, double yy, Type ttype, IdType uid) : x(xx), @@ -116,14 +121,8 @@ DiagramPoint::DiagramPoint(double xx, double yy, Type ttype, IdType uid) : type(ttype), id(uid) { - //if ( xx < 0 ) - //throw "Negative x coordinate"; - //if ( yy < 0) - //throw "Negative y coordinate"; - //if ( yy < xx ) - //throw "Point is below the diagonal"; if ( yy == xx and ttype != DiagramPoint::DIAG) - throw "Point on the main diagonal must have DIAG type"; + throw std::runtime_error("Point on the main diagonal must have DIAG type"); } diff --git a/geom_bottleneck/bottleneck/src/bottleneck.cpp b/geom_bottleneck/bottleneck/src/bottleneck.cpp index a5009c5..365995a 100644 --- a/geom_bottleneck/bottleneck/src/bottleneck.cpp +++ b/geom_bottleneck/bottleneck/src/bottleneck.cpp @@ -151,7 +151,7 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) const double approxDist = 0.5 * ( interval.first + interval.second); const double minDist = interval.first; const double maxDist = interval.second; - //std::cerr << std::setprecision(15) << "minDist = " << minDist << ", maxDist = " << maxDist << std::endl; + //std::cout << std::setprecision(15) << "minDist = " << minDist << ", maxDist = " << maxDist << std::endl; if ( delta == 0 ) { return interval.first; } @@ -163,17 +163,17 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) pointsA.push_back(ptA); } - //std::vector killDist; - //for(auto ptA : A) { - //for(auto ptB : B) { - //if ( distLInf(ptA, ptB) > minDist and distLInf(ptA, ptB) < maxDist) { - //killDist.push_back(distLInf(ptA, ptB)); - //std::cout << ptA << ", " << ptB << std::endl; + //std::vector killdist; + //for(auto pta : a) { + //for(auto ptb : b) { + //if ( distlinf(pta, ptb) > mindist and distlinf(pta, ptb) < maxdist) { + //killdist.push_back(distlinf(pta, ptb)); + //std::cout << pta << ", " << ptb << std::endl; //} //} //} - //std::sort(killDist.begin(), killDist.end()); - //for(auto d : killDist) { + //std::sort(killdist.begin(), killdist.end()); + //for(auto d : killdist) { //std::cout << d << std::endl; //} //std::cout << "*************" << std::endl; @@ -257,7 +257,7 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) std::sort(yCentersVec.begin(), yCentersVec.end(), compLambda); - // std::cout << "Sorted vector of y-centers:" << std::endl; + // std::cout << "Sorted vector of y-centers:" << std::endl; //for(auto coordPtPair : yCentersVec) { //std::cout << coordPtPair.first << ", id = " << coordPtPair.second.id << std::endl; //} @@ -271,13 +271,6 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) std::make_pair(ptB.getRealY() - delta, ptB), compLambda); - //if (ptB.id == 316) { - //std::cout << itStart - yCentersVec.begin() << " " << distLInf(itStart->second, ptB) << std::endl; - //std::cout << "maxDist = " << maxDist << std::endl; - //std::cout << "minDist = " << minDist << std::endl; - //double pwDistDebug = distLInf(itStart->second, ptB); - //std::cout << ( pwDistDebug >= minDist and pwDistDebug <= maxDist) << std::endl; - //} for(auto iterA = itStart; iterA < yCentersVec.end(); ++iterA) { if ( ptB.getRealY() < iterA->first - delta) { @@ -286,19 +279,16 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) double pwDist = distLInf(iterA->second, ptB); //std::cout << 1000*minDist << " <= " << 1000*pwDist << " <= " << 1000*maxDist << std::endl; if (pwDist >= minDist and pwDist <= maxDist) { - //if (ptB.id == 316) { - //std::cout << "adding " << pwDist << std::endl; - //} pairwiseDist.push_back(pwDist); } } } } - //std::cerr << "pairwiseDist.size = " << pairwiseDist.size() << " out of " << A.size() * A.size() << std::endl; + //std::cout << "pairwiseDist.size = " << pairwiseDist.size() << " out of " << A.size() * A.size() << std::endl; std::sort(pairwiseDist.begin(), pairwiseDist.end()); //for(auto ddd : pairwiseDist) { - //std::cerr << std::setprecision(15) << ddd << std::endl; + //std::cout << std::setprecision(15) << ddd << std::endl; //} return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist); @@ -516,7 +506,9 @@ bool readDiagramPointSet(const char* fname, std::vector> x >> y)) { +#ifndef FOR_R_TDA std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl; +#endif return false; } result.push_back(std::make_pair(x,y)); diff --git a/geom_bottleneck/bottleneck/src/bound_match.cpp b/geom_bottleneck/bottleneck/src/bound_match.cpp index 06d3b67..a9ec93a 100644 --- a/geom_bottleneck/bottleneck/src/bound_match.cpp +++ b/geom_bottleneck/bottleneck/src/bound_match.cpp @@ -18,9 +18,12 @@ along with GeomBottleneck. If not, see . */ -#include #include +#include "def_debug_bt.h" #include "bound_match.h" +#ifndef FOR_R_TDA +#include +#endif namespace geom_bt { /*static void printDebug(//bool isDebug, std::string s)*/ @@ -82,6 +85,7 @@ namespace geom_bt { //#endif /*}*/ +#ifndef FOR_R_TDA std::ostream& operator<<(std::ostream& output, const Matching& m) { output << "Matching: " << m.AToB.size() << " pairs ("; @@ -94,6 +98,7 @@ std::ostream& operator<<(std::ostream& output, const Matching& m) } return output; } +#endif void Matching::sanityCheck() const { @@ -103,8 +108,11 @@ void Matching::sanityCheck() const auto bToAPair = BToA.find(aToBPair.second); assert(bToAPair != BToA.end()); if (aToBPair.first != bToAPair->second) { +#ifndef FOR_R_TDA std::cerr << "failed assertion, in aToB " << aToBPair.first; std::cerr << ", in bToA " << bToAPair->second << std::endl; +#endif + assert(false); } assert( aToBPair.first == bToAPair->second); } @@ -151,7 +159,9 @@ void Matching::checkAugPath(const Path& augPath) const for(size_t idx = 0; idx < augPath.size(); ++idx) { bool mustBeExposed { idx == 0 or idx == augPath.size() - 1 }; if (isExposed(augPath[idx]) != mustBeExposed) { +#ifndef FOR_R_TDA std::cerr << "mustBeExposed = " << mustBeExposed << ", idx = " << idx << ", point " << augPath[idx] << std::endl; +#endif } assert( isExposed(augPath[idx]) == mustBeExposed ); DiagramPoint matchedVertex {0.0, 0.0, DiagramPoint::DIAG, 1}; @@ -275,10 +285,12 @@ bool BoundMatchOracle::buildAugmentingPath(const DiagramPoint startVertex, Path& // layer DiagramPoint nextVertexA {0.0, 0.0, DiagramPoint::DIAG, 1}; if (!M.getMatchedVertex(nextVertexB, nextVertexA)) { +#ifndef FOR_R_TDA std::cerr << "Vertices in even layers must be matched! Unmatched: "; std::cerr << nextVertexB << std::endl; std::cerr << evenLayerIdx << "; " << layerGraph.size() << std::endl; - throw "Unmatched vertex in even layer"; +#endif + throw std::runtime_error("Unmatched vertex in even layer"); } else { assert( ! (nextVertexA.getRealX() == 0 and nextVertexA.getRealY() == 0) ); result.push_back(nextVertexA); @@ -383,8 +395,10 @@ bool BoundMatchOracle::buildMatchingForThreshold(const double r) } if (augmentingPaths.empty()) { //printDebug(isDebug,"augmenting paths must exist, but were not found!", M); +#ifndef FOR_R_TDA std::cerr << "augmenting paths must exist, but were not found!" << std::endl; - throw "bad epsilon?"; +#endif + throw std::runtime_error("bad epsilon?"); } // swap all augmenting paths with matching to increase it //printDebug(isDebug,"before increase with augmenting paths:", M); @@ -403,6 +417,7 @@ bool BoundMatchOracle::buildMatchingForThreshold(const double r) void BoundMatchOracle::printLayerGraph(void) { #ifdef DEBUG_BOUND_MATCH +#ifndef FOR_R_TDA for(auto& layer : layerGraph) { std::cout << "{ "; for(auto& p : layer) { @@ -411,6 +426,7 @@ void BoundMatchOracle::printLayerGraph(void) std::cout << "\b\b }" << std::endl; } #endif +#endif } void BoundMatchOracle::buildLayerGraph(double r) diff --git a/geom_bottleneck/bottleneck/src/neighb_oracle.cpp b/geom_bottleneck/bottleneck/src/neighb_oracle.cpp index 40d043f..5005bc4 100644 --- a/geom_bottleneck/bottleneck/src/neighb_oracle.cpp +++ b/geom_bottleneck/bottleneck/src/neighb_oracle.cpp @@ -184,8 +184,10 @@ void NeighbOracleAnn::deletePoint(const DiagramPoint& p) diagonalPoints.erase(p, false); kdTree->delete_point(pointIdx); #ifdef DEBUG_NEIGHBOUR_ORACLE +#ifndef FOR_R_TDA kdTree->Print(ANNtrue, std::cout); #endif +#endif } bool NeighbOracleAnn::getNeighbour(const DiagramPoint& q, DiagramPoint& result) const 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 diff --git a/geom_matching/wasserstein/src/auction_oracle.cpp b/geom_matching/wasserstein/src/auction_oracle.cpp index 088936b..9c764da 100644 --- a/geom_matching/wasserstein/src/auction_oracle.cpp +++ b/geom_matching/wasserstein/src/auction_oracle.cpp @@ -68,8 +68,6 @@ AuctionOracleLazyHeap::AuctionOracleLazyHeap(const std::vector& b, weightMatrix.reserve(b.size()); //const double maxDistUpperBound = 3 * getFurthestDistance3Approx(b, g); //weightAdjConst = pow(maxDistUpperBound, wassersteinPower); - //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(b, g) << std::endl; - //std::cout << "in AuctionOracleLazyHeap weightAdjConst = " << weightAdjConst << std::endl; // init weight matrix for(const auto& pointA : bidders) { std::vector weightVec; @@ -134,7 +132,9 @@ void AuctionOracleLazyHeap::setPrice(IdxType itemIdx, double newPrice) { assert( prices.at(itemIdx) < newPrice ); #ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA std::cout << "price incremented by " << prices.at(itemIdx) - newPrice << std::endl; +#endif #endif prices[itemIdx] = newPrice; // lazy: record the moment we updated the price of the items, @@ -169,8 +169,10 @@ DebugOptimalBid AuctionOracleLazyHeap::getOptimalBidDebug(IdxType bidderIdx) result.secondBestItemValue = topIter->second; result.secondBestItemIdx = topIter->first; +#ifndef FOR_R_TDA //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; +#endif return result; } @@ -190,8 +192,10 @@ IdxValPair AuctionOracleLazyHeap::getOptimalBid(const IdxType bidderIdx) ++topIter; // now points to the second-best items double secondBestItemValue = topIter->second; +#ifndef FOR_R_TDA //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; +#endif // bid value: price + value difference + epsilon return std::make_pair(bestItemIdx, @@ -221,8 +225,6 @@ AuctionOracleLazyHeapRestricted::AuctionOracleLazyHeapRestricted(const std::vect weightMatrix.reserve(b.size()); //const double maxDistUpperBound = 3 * getFurthestDistance3Approx(b, g); //weightAdjConst = pow(maxDistUpperBound, wassersteinPower); - //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(b, g) << std::endl; - //std::cout << "in AuctionOracleLazyHeapRestricted weightAdjConst = " << weightAdjConst << std::endl; // init weight matrix for(const auto& pointA : bidders) { std::vector weightVec; @@ -319,7 +321,9 @@ void AuctionOracleLazyHeapRestricted::setPrice(IdxType itemIdx, double newPrice) { assert( prices.at(itemIdx) < newPrice ); #ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA std::cout << "price incremented by " << prices.at(itemIdx) - newPrice << std::endl; +#endif #endif prices[itemIdx] = newPrice; if (items[itemIdx].isNormal() ) { @@ -426,7 +430,6 @@ DebugOptimalBid AuctionOracleLazyHeapRestricted::getOptimalBidDebug(IdxType bidd //debugNaiveResult.secondBestItemIdx = itemIdx; //} //} - ////std::cout << "got naive result" << std::endl; //if ( fabs( debugMyResult.bestItemValue - debugNaiveResult.bestItemValue ) > 1e-6 or //fabs( debugNaiveResult.secondBestItemValue - debugMyResult.secondBestItemValue) > 1e-6 ) { @@ -648,7 +651,6 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector& _bidde } DnnTraits traits; traits.internal_p = internal_p; - //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl; kdtree = new dnn::KDTree(traits, dnnPointHandles, wassersteinPower); size_t dnnItemIdxAll { 0 }; @@ -658,7 +660,6 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector& _bidde DnnPoint p(dnnItemIdxAll++); p[0] = g.getRealX(); p[1] = g.getRealY(); - //std::cout << "to all tree: " << p[0] << ", " << p[1] << std::endl; dnnPointsAll.push_back(p); assert(dnnItemIdxAll == dnnPointsAll.size()); } @@ -666,7 +667,6 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector& _bidde for(size_t i = 0; i < dnnPointsAll.size(); ++i) { dnnPointHandlesAll.push_back(&dnnPointsAll[i]); } - //std::cout << "kdtreeAll: " << dnnPointHandlesAll.size() << " points" << std::endl; kdtreeAll = new dnn::KDTree(traits, dnnPointHandlesAll, wassersteinPower); size_t handleIdx {0}; @@ -677,13 +677,9 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector& _bidde } } //to-do: remove maxVal from - //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(_bidders, _items) << std::endl; maxVal = 3*getFurthestDistance3Approx(_bidders, _items); maxVal = pow(maxVal, wassersteinPower); weightAdjConst = maxVal; - //std::cout << "AuctionOracleKDTree: weightAdjConst = " << weightAdjConst << std::endl; - //std::cout << "AuctionOracleKDTree constructor done" << std::endl; - // for debug } DebugOptimalBid AuctionOracleKDTree::getOptimalBidDebug(IdxType bidderIdx) @@ -694,16 +690,12 @@ DebugOptimalBid AuctionOracleKDTree::getOptimalBidDebug(IdxType bidderIdx) bidderDnn[0] = bidder.getRealX(); bidderDnn[1] = bidder.getRealY(); - //std::cout << "bidder.x = " << bidderDnn[0] << std::endl; - //std::cout << "bidder.y = " << bidderDnn[1] << std::endl; - std::vector candItems; if ( bidder.isDiagonal() ) { // auto twoBestItems = kdtree->findK(bidderDnn, 2); - //std::cout << "twoBestItems for non-diag: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; candItems.push_back( std::make_pair(twoBestItems[0].p->id(), twoBestItems[0].d) ); candItems.push_back( std::make_pair(twoBestItems[1].p->id(), twoBestItems[1].d) ); assert(diagItemsHeap.size() > 1); @@ -719,7 +711,6 @@ DebugOptimalBid AuctionOracleKDTree::getOptimalBidDebug(IdxType bidderIdx) assert(candItems[1].second >= candItems[0].second); } else { auto twoBestItems = kdtreeAll->findK(bidderDnn, 2); - //std::cout << "twoBestItems for all: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; candItems.push_back( std::make_pair(twoBestItems[0].p->id(), twoBestItems[0].d) ); candItems.push_back( std::make_pair(twoBestItems[1].p->id(), twoBestItems[1].d) ); //size_t projItemIdx { biddersProjIndices.at(bidderIdx) }; @@ -734,7 +725,6 @@ DebugOptimalBid AuctionOracleKDTree::getOptimalBidDebug(IdxType bidderIdx) result.secondBestItemIdx = candItems[1].first; result.bestItemValue = candItems[0].second; result.secondBestItemValue = candItems[1].second; - //std::cout << "got result: " << result << std::endl; //double bestItemsPrice = prices[bestItemIdx]; //if (items[result.bestItemIdx].type == DiagramPoint::DIAG) { //double bestItemValue1 = pow( distLp(bidder, items[result.bestItemIdx]), q) + prices[result.bestItemIdx]; @@ -820,18 +810,12 @@ void AuctionOracleKDTree::setPrice(IdxType itemIdx, double newPrice) assert(newPrice > prices.at(itemIdx)); prices[itemIdx] = newPrice; if ( items[itemIdx].isNormal() ) { - //std::cout << "before increasing weight in kdtree " << std::endl; - //std::cout << kdtreeItems.at(itemIdx) << std::endl; assert(0 <= itemIdx and itemIdx < kdtreeItems.size()); assert(0 <= kdtreeItems[itemIdx] and kdtreeItems[itemIdx] < dnnPointHandles.size()); kdtree->increase_weight( dnnPointHandles[kdtreeItems[itemIdx]], newPrice); kdtreeAll->increase_weight( dnnPointHandlesAll[itemIdx], newPrice); - //std::cout << "after increasing weight in kdtree" << std::endl; } else { - //std::cout << "before decreasing weight in heap" << std::endl; - //std::cout << "diagHeapHandles.size = " << diagHeapHandles.size() << std::endl; kdtreeAll->increase_weight( dnnPointHandlesAll[itemIdx], newPrice); - //std::cout << "after increasing weight in kdtree" << std::endl; assert(diagHeapHandles.size() > heapHandlesIndices.at(itemIdx)); diagItemsHeap.decrease(diagHeapHandles[heapHandlesIndices[itemIdx]], std::make_pair(itemIdx, newPrice)); } @@ -986,7 +970,6 @@ AuctionOracleKDTreeRestricted::AuctionOracleKDTreeRestricted(const std::vector(traits, dnnPointHandles, wassersteinPower); size_t handleIdx {0}; @@ -997,12 +980,9 @@ AuctionOracleKDTreeRestricted::AuctionOracleKDTreeRestricted(const std::vectorid() }; double bestNormalItemValue { twoBestItems[0].d }; size_t secondBestNormalItemIdx { twoBestItems[1].p->id() }; @@ -1162,9 +1138,6 @@ IdxValPair AuctionOracleKDTreeRestricted::getOptimalBid(IdxType bidderIdx) DiagramPoint bidder = bidders[bidderIdx]; - //std::cout << "bidder.x = " << bidderDnn[0] << std::endl; - //std::cout << "bidder.y = " << bidderDnn[1] << std::endl; - // corresponding point is always considered as a candidate // if bidder is a diagonal point, projItem is a normal point, // and vice versa. @@ -1226,7 +1199,6 @@ IdxValPair AuctionOracleKDTreeRestricted::getOptimalBid(IdxType bidderIdx) bidderDnn[0] = bidder.getRealX(); bidderDnn[1] = bidder.getRealY(); auto twoBestItems = kdtree->findK(bidderDnn, 2); - //std::cout << "twoBestItems for all: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; size_t bestNormalItemIdx { twoBestItems[0].p->id() }; double bestNormalItemValue { twoBestItems[0].d }; // if there is only one off-diagonal point in the second diagram, @@ -1268,15 +1240,10 @@ void AuctionOracleKDTreeRestricted::setPrice(IdxType itemIdx, double newPrice) assert(newPrice > prices.at(itemIdx)); prices[itemIdx] = newPrice; if ( items[itemIdx].isNormal() ) { - //std::cout << "before increasing weight in kdtree " << std::endl; - //std::cout << kdtreeItems.at(itemIdx) << std::endl; assert(0 <= itemIdx and itemIdx < kdtreeItems.size()); assert(0 <= kdtreeItems[itemIdx] and kdtreeItems[itemIdx] < dnnPointHandles.size()); kdtree->increase_weight( dnnPointHandles[kdtreeItems[itemIdx]], newPrice); - //std::cout << "after increasing weight in kdtree" << std::endl; } else { - //std::cout << "before decreasing weight in heap" << std::endl; - //std::cout << "diagHeapHandles.size = " << diagHeapHandles.size() << std::endl; assert(diagHeapHandles.size() > heapHandlesIndices.at(itemIdx)); diagItemsHeap.decrease(diagHeapHandles[heapHandlesIndices[itemIdx]], std::make_pair(itemIdx, newPrice)); bestDiagonalItemsComputed = false; @@ -1300,10 +1267,10 @@ void AuctionOracleKDTreeRestricted::setEpsilon(double newVal) std::ostream& operator<< (std::ostream& output, const DebugOptimalBid& db) { - std::cout << "bestItemValue = " << db.bestItemValue; - std::cout << "; bestItemIdx = " << db.bestItemIdx; - std::cout << "; secondBestItemValue = " << db.secondBestItemValue; - std::cout << "; secondBestItemIdx = " << db.secondBestItemIdx; + output << "bestItemValue = " << db.bestItemValue; + output << "; bestItemIdx = " << db.bestItemIdx; + output << "; secondBestItemValue = " << db.secondBestItemValue; + output << "; secondBestItemIdx = " << db.secondBestItemIdx; return output; } diff --git a/geom_matching/wasserstein/src/auction_runner_gs.cpp b/geom_matching/wasserstein/src/auction_runner_gs.cpp index 4ff495f..bd25442 100644 --- a/geom_matching/wasserstein/src/auction_runner_gs.cpp +++ b/geom_matching/wasserstein/src/auction_runner_gs.cpp @@ -27,6 +27,7 @@ derivative works thereof, in binary and source code form. #include +#include #include #include #include @@ -36,6 +37,10 @@ derivative works thereof, in binary and source code form. #include "auction_runner_gs.h" #include "wasserstein.h" +#ifdef FOR_R_TDA +#include "Rcpp.h" +#endif + //#define PRINT_DETAILED_TIMING namespace geom_ws { @@ -144,11 +149,13 @@ void AuctionRunnerGS::runAuction(void) double denominator = currentResult - numBidders * oracle->getEpsilon(); currentResult = pow(currentResult, 1.0 / wassersteinPower); #ifdef PRINT_DETAILED_TIMING +#ifndef FOR_R_TDA iterResults.push_back(currentResult); iterTimes.push_back(hrClock.now()); std::cout << "Iteration " << iterNum << " finished. "; std::cout << "Current result is " << currentResult << ", epsilon = " << oracle->getEpsilon() << std::endl; std::cout << "Number of rounds (cumulative): " << numRounds << std::endl; +#endif #endif if ( denominator <= 0 ) { //std::cout << "Epsilon is too big." << std::endl; @@ -157,8 +164,10 @@ void AuctionRunnerGS::runAuction(void) denominator = pow(denominator, 1.0 / wassersteinPower); double numerator = currentResult - denominator; #ifdef PRINT_DETAILED_TIMING +#ifndef FOR_R_TDA std::cout << " numerator: " << numerator << " denominator: " << denominator; std::cout << "; error bound: " << numerator / denominator << std::endl; +#endif #endif // if relative error is greater than delta, continue notDone = ( numerator / denominator > delta ); @@ -166,13 +175,16 @@ void AuctionRunnerGS::runAuction(void) // decrease epsilon for the next iteration oracle->setEpsilon( oracle->getEpsilon() / epsilonCommonRatio ); if (iterNum > maxIterNum) { +#ifndef FOR_R_TDA std::cerr << "Maximum iteration number exceeded, exiting. Current result is:"; std::cerr << wassersteinDistance << std::endl; - std::exit(1); +#endif + throw std::runtime_error("Maximum iteration number exceeded"); } } while ( notDone ); //printMatching(); #ifdef PRINT_DETAILED_TIMING +#ifndef FOR_R_TDA for(size_t iterIdx = 0; iterIdx < iterResults.size(); ++iterIdx) { double trueRelError = ( iterResults.at(iterIdx) - currentResult ) / currentResult; auto iterCumulativeTime = iterTimes.at(iterIdx) - startMoment; @@ -183,6 +195,7 @@ void AuctionRunnerGS::runAuction(void) ", iteration time " << iterTime.count() << std::endl; } #endif +#endif } void AuctionRunnerGS::runAuctionPhase(void) @@ -200,15 +213,22 @@ void AuctionRunnerGS::runAuctionPhase(void) assignItemToBidder(optimalBid.first, bidderIdx); oracle->setPrice(optimalItemIdx, bidValue); //printDebug(); +#ifdef FOR_R_TDA + if ( numRounds % 10000 == 0 ) { + Rcpp::checkUserInterrupt(); + } +#endif } while (not unassignedBidders.empty()); //std::cout << "runAuctionPhase finished" << std::endl; #ifdef DEBUG_AUCTION for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { if ( biddersToItems[bidderIdx] < 0) { +#ifndef FOR_R_TDA std::cerr << "After auction terminated bidder " << bidderIdx; std::cerr << " has no items assigned" << std::endl; - throw "Auction did not give a perfect matching"; +#endif + throw std::runtime_error("Auction did not give a perfect matching"); } } #endif @@ -225,6 +245,7 @@ double AuctionRunnerGS::getDistanceToQthPowerInternal(void) auto pB = items[biddersToItems[bIdx]]; result += pow(distLp(pA, pB, internal_p), wassersteinPower); } + wassersteinCost = result; wassersteinDistance = pow(result, 1.0 / wassersteinPower); return result; } @@ -235,11 +256,20 @@ double AuctionRunnerGS::getWassersteinDistance(void) return wassersteinDistance; } +double AuctionRunnerGS::getWassersteinCost(void) +{ + runAuction(); + return wassersteinCost; +} + + + // Debug routines void AuctionRunnerGS::printDebug(void) { #ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA sanityCheck(); std::cout << "**********************" << std::endl; std::cout << "Current assignment:" << std::endl; @@ -259,6 +289,7 @@ void AuctionRunnerGS::printDebug(void) } std::cout << "**********************" << std::endl; #endif +#endif } @@ -266,13 +297,17 @@ void AuctionRunnerGS::sanityCheck(void) { #ifdef DEBUG_AUCTION if (biddersToItems.size() != numBidders) { +#ifndef FOR_R_TDA std::cerr << "Wrong size of biddersToItems, must be " << numBidders << ", is " << biddersToItems.size() << std::endl; - throw "Wrong size of biddersToItems"; +#endif + throw std::runtime_error("Wrong size of biddersToItems"); } if (itemsToBidders.size() != numBidders) { +#ifndef FOR_R_TDA std::cerr << "Wrong size of itemsToBidders, must be " << numBidders << ", is " << itemsToBidders.size() << std::endl; - throw "Wrong size of itemsToBidders"; +#endif + throw std::runtime_error("Wrong size of itemsToBidders"); } for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { @@ -281,18 +316,22 @@ void AuctionRunnerGS::sanityCheck(void) if ( std::count(biddersToItems.begin(), biddersToItems.end(), biddersToItems[bidderIdx]) > 1 ) { +#ifndef FOR_R_TDA std::cerr << "Item " << biddersToItems[bidderIdx]; std::cerr << " appears in biddersToItems more than once" << std::endl; - throw "Duplicate in biddersToItems"; +#endif + throw std::runtime_error("Duplicate in biddersToItems"); } if (itemsToBidders.at(biddersToItems[bidderIdx]) != static_cast(bidderIdx)) { +#ifndef FOR_R_TDA std::cerr << "Inconsitency: bidderIdx = " << bidderIdx; std::cerr << ", itemIdx in biddersToItems = "; std::cerr << biddersToItems[bidderIdx]; std::cerr << ", bidderIdx in itemsToBidders = "; std::cerr << itemsToBidders[biddersToItems[bidderIdx]] << std::endl; - throw "inconsistent mapping"; +#endif + throw std::runtime_error("inconsistent mapping"); } } } @@ -304,18 +343,22 @@ void AuctionRunnerGS::sanityCheck(void) if ( std::count(itemsToBidders.begin(), itemsToBidders.end(), itemsToBidders[itemIdx]) > 1 ) { +#ifndef FOR_R_TDA std::cerr << "Bidder " << itemsToBidders[itemIdx]; std::cerr << " appears in itemsToBidders more than once" << std::endl; - throw "Duplicate in itemsToBidders"; +#endif + throw std::runtime_error("Duplicate in itemsToBidders"); } // check for consistency if (biddersToItems.at(itemsToBidders[itemIdx]) != static_cast(itemIdx)) { +#ifndef FOR_R_TDA std::cerr << "Inconsitency: itemIdx = " << itemIdx; std::cerr << ", bidderIdx in itemsToBidders = "; std::cerr << itemsToBidders[itemIdx]; std::cerr << ", itemIdx in biddersToItems= "; std::cerr << biddersToItems[itemsToBidders[itemIdx]] << std::endl; - throw "inconsistent mapping"; +#endif + throw std::runtime_error("inconsistent mapping"); } } } @@ -325,6 +368,7 @@ void AuctionRunnerGS::sanityCheck(void) void AuctionRunnerGS::printMatching(void) { //#ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA sanityCheck(); for(size_t bIdx = 0; bIdx < biddersToItems.size(); ++bIdx) { if (biddersToItems[bIdx] >= 0) { @@ -335,6 +379,7 @@ void AuctionRunnerGS::printMatching(void) assert(false); } } +#endif //#endif } diff --git a/geom_matching/wasserstein/src/auction_runner_jac.cpp b/geom_matching/wasserstein/src/auction_runner_jac.cpp index 08ed202..b892643 100644 --- a/geom_matching/wasserstein/src/auction_runner_jac.cpp +++ b/geom_matching/wasserstein/src/auction_runner_jac.cpp @@ -37,10 +37,10 @@ derivative works thereof, in binary and source code form. namespace geom_ws { // ***************************** -// AuctionRunnerJak +// AuctionRunnerJac // ***************************** -AuctionRunnerJak::AuctionRunnerJak(const std::vector& A, const std::vector& B, const double q, const double _delta, const double _internal_p) : +AuctionRunnerJac::AuctionRunnerJac(const std::vector& A, const std::vector& B, const double q, const double _delta, const double _internal_p) : bidders(A), items(B), numBidders(A.size()), @@ -57,7 +57,7 @@ AuctionRunnerJak::AuctionRunnerJak(const std::vector& A, const std oracle = std::unique_ptr(new AuctionOracle(bidders, items, wassersteinPower, internal_p)); } -void AuctionRunnerJak::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) +void AuctionRunnerJac::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) { //sanityCheck(); IdxType myOldItem = biddersToItems[bidderIdx]; @@ -79,7 +79,9 @@ void AuctionRunnerJak::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) // RE: this cannot be necessary. I submitted the best bid, hence I was // an unassigned bidder. if (myOldItem != -1) { +#ifndef FOR_R_TDA std::cout << "This is not happening" << std::endl; +#endif assert(false); itemsToBidders[myOldItem] = -1; } @@ -88,7 +90,9 @@ void AuctionRunnerJak::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) biddersToItems[currItemOwner] = myOldItem; // add the old owner of bids to the list of if ( -1 != myOldItem ) { +#ifndef FOR_R_TDA std::cout << "This is not happening" << std::endl; +#endif assert(false); // if I had something, update itemsToBidders, too // RE: nonsense: if I had something, I am not unassigned and did not @@ -103,18 +107,17 @@ void AuctionRunnerJak::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) } -void AuctionRunnerJak::assignToBestBidder(IdxType itemIdx) +void AuctionRunnerJac::assignToBestBidder(IdxType itemIdx) { assert( itemIdx >= 0 and itemIdx < static_cast(numItems) ); assert( bidTable[itemIdx].first != -1); IdxValPair bestBid { bidTable[itemIdx] }; assignGoodToBidder(itemIdx, bestBid.first); - //std::cout << "About to call setPrice" << std::endl; oracle->setPrice(itemIdx, bestBid.second); //dynamic_cast(oracle)->setNai } -void AuctionRunnerJak::clearBidTable(void) +void AuctionRunnerJac::clearBidTable(void) { for(auto& itemWithBidIdx : itemsWithBids) { itemReceivedBidVec[itemWithBidIdx] = 0; @@ -124,7 +127,7 @@ void AuctionRunnerJak::clearBidTable(void) itemsWithBids.clear(); } -void AuctionRunnerJak::submitBid(IdxType bidderIdx, const IdxValPair& itemsBidValuePair) +void AuctionRunnerJac::submitBid(IdxType bidderIdx, const IdxValPair& itemsBidValuePair) { IdxType itemIdx = itemsBidValuePair.first; double bidValue = itemsBidValuePair.second; @@ -139,9 +142,10 @@ void AuctionRunnerJak::submitBid(IdxType bidderIdx, const IdxValPair& itemsBidVa } } -void AuctionRunnerJak::printDebug(void) +void AuctionRunnerJac::printDebug(void) { #ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA sanityCheck(); std::cout << "**********************" << std::endl; std::cout << "Current assignment:" << std::endl; @@ -168,9 +172,10 @@ void AuctionRunnerJak::printDebug(void) //} std::cout << "**********************" << std::endl; #endif +#endif } -void AuctionRunnerJak::flushAssignment(void) +void AuctionRunnerJac::flushAssignment(void) { for(auto& b2g : biddersToItems) { b2g = -1; @@ -181,7 +186,7 @@ void AuctionRunnerJak::flushAssignment(void) //oracle->adjustPrices(); } -void AuctionRunnerJak::runAuction(void) +void AuctionRunnerJac::runAuction(void) { // relative error // choose some initial epsilon @@ -213,15 +218,19 @@ void AuctionRunnerJak::runAuction(void) // decrease epsilon for the next iteration oracle->setEpsilon( oracle->getEpsilon() / epsilonCommonRatio ); if (iterNum > maxIterNum) { +#ifndef FOR_R_TDA std::cerr << "Maximum iteration number exceeded, exiting. Current result is:"; std::cerr << wassersteinDistance << std::endl; std::exit(1); +#else + throw std::runtime_error("Max. iter. exceeded"); +#endif } } while ( notDone ); //printMatching(); } -void AuctionRunnerJak::runAuctionPhase(void) +void AuctionRunnerJac::runAuctionPhase(void) { //std::cout << "Entered runAuctionPhase" << std::endl; //int numUnassignedBidders { 0 }; @@ -266,7 +275,7 @@ void AuctionRunnerJak::runAuctionPhase(void) } // assertion: the matching must be perfect -double AuctionRunnerJak::getDistanceToQthPowerInternal(void) +double AuctionRunnerJac::getDistanceToQthPowerInternal(void) { sanityCheck(); double result = 0.0; @@ -276,19 +285,29 @@ double AuctionRunnerJak::getDistanceToQthPowerInternal(void) auto pB = items[biddersToItems[bIdx]]; result += pow(distLp(pA, pB, internal_p), wassersteinPower); } + wassersteinCost = result; wassersteinDistance = pow(result, 1.0 / wassersteinPower); return result; } -double AuctionRunnerJak::getWassersteinDistance(void) +double AuctionRunnerJac::getWassersteinDistance(void) { runAuction(); return wassersteinDistance; } -void AuctionRunnerJak::sanityCheck(void) +double AuctionRunnerJac::getWassersteinCost(void) +{ + runAuction(); + return wassersteinCost; +} + + + +void AuctionRunnerJac::sanityCheck(void) { #ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA if (biddersToItems.size() != numBidders) { std::cerr << "Wrong size of biddersToItems, must be " << numBidders << ", is " << biddersToItems.size() << std::endl; throw "Wrong size of biddersToItems"; @@ -344,11 +363,13 @@ void AuctionRunnerJak::sanityCheck(void) } } #endif +#endif } -void AuctionRunnerJak::printMatching(void) +void AuctionRunnerJac::printMatching(void) { //#ifdef DEBUG_AUCTION +#ifndef FOR_R_TDA sanityCheck(); for(size_t bIdx = 0; bIdx < biddersToItems.size(); ++bIdx) { if (biddersToItems[bIdx] >= 0) { @@ -359,6 +380,7 @@ void AuctionRunnerJak::printMatching(void) assert(false); } } +#endif //#endif } diff --git a/geom_matching/wasserstein/src/basic_defs.cpp b/geom_matching/wasserstein/src/basic_defs.cpp index d228123..a46e6aa 100644 --- a/geom_matching/wasserstein/src/basic_defs.cpp +++ b/geom_matching/wasserstein/src/basic_defs.cpp @@ -30,6 +30,10 @@ derivative works thereof, in binary and source code form. #include #include "basic_defs_ws.h" +#ifndef FOR_R_TDA +#include +#endif + namespace geom_ws { // Point @@ -43,11 +47,14 @@ bool Point::operator!=(const Point& other) const return !(*this == other); } + +#ifndef FOR_R_TDA std::ostream& operator<<(std::ostream& output, const Point p) { output << "(" << p.x << ", " << p.y << ")"; return output; } +#endif double sqrDist(const Point& a, const Point& b) { @@ -97,6 +104,7 @@ double distLp(const DiagramPoint& a, const DiagramPoint& b, const double p) } +#ifndef FOR_R_TDA std::ostream& operator<<(std::ostream& output, const DiagramPoint p) { if ( p.type == DiagramPoint::DIAG ) { @@ -106,7 +114,7 @@ std::ostream& operator<<(std::ostream& output, const DiagramPoint p) } return output; } - +#endif DiagramPoint::DiagramPoint(double xx, double yy, Type ttype) : x(xx), diff --git a/geom_matching/wasserstein/src/wasserstein.cpp b/geom_matching/wasserstein/src/wasserstein.cpp index cac811b..fc1b662 100644 --- a/geom_matching/wasserstein/src/wasserstein.cpp +++ b/geom_matching/wasserstein/src/wasserstein.cpp @@ -36,7 +36,7 @@ derivative works thereof, in binary and source code form. #ifdef GAUSS_SEIDEL_AUCTION #include "auction_runner_gs.h" #else -#include "auction_runner_jak.h" +#include "auction_runner_jac.h" #endif namespace geom_ws { @@ -50,29 +50,78 @@ double wassersteinDistVec(const std::vector& A, const double _epsFactor) { if (q < 1) { +#ifndef FOR_R_TDA std::cerr << "Wasserstein distance not defined for q = " << q << ", must be >= 1" << std::endl; - throw "Bad q in Wasserstein"; +#endif + throw std::runtime_error("Bad q in Wasserstein " + std::to_string(q)); } if (delta < 0.0) { +#ifndef FOR_R_TDA std::cerr << "Relative error " << delta << ", must be > 0" << std::endl; - throw "Bad delta in Wasserstein"; +#endif + throw std::runtime_error("Bad delta in Wasserstein " + std::to_string(delta)); } if (_initialEpsilon < 0.0) { +#ifndef FOR_R_TDA std::cerr << "Initial epsilon = " << _initialEpsilon << ", must be non-negative" << std::endl; - throw "Bad delta in Wasserstein"; +#endif + throw std::runtime_error("Bad initial epsilon in Wasserstein" + std::to_string(_initialEpsilon)); } if (_epsFactor < 0.0) { +#ifndef FOR_R_TDA std::cerr << "Epsilon factor = " << _epsFactor << ", must be non-negative" << std::endl; - throw "Bad delta in Wasserstein"; +#endif + throw std::runtime_error("Bad epsilon factor in Wasserstein " + std::to_string(_epsFactor)); } + #ifdef GAUSS_SEIDEL_AUCTION AuctionRunnerGS auction(A, B, q, delta, _internal_p, _initialEpsilon, _epsFactor); #else - AuctionRunnerJak auction(A, B, q, delta, _internal_p); + AuctionRunnerJac auction(A, B, q, delta, _internal_p); #endif return auction.getWassersteinDistance(); } +double wassersteinCostVec(const std::vector& A, + const std::vector& B, + const double q, + const double delta, + const double _internal_p, + const double _initialEpsilon, + const double _epsFactor) +{ + if (q < 1) { +#ifndef FOR_R_TDA + std::cerr << "Wasserstein distance not defined for q = " << q << ", must be >= 1" << std::endl; +#endif + throw std::runtime_error("Bad q in Wasserstein " + std::to_string(q)); + } + if (delta < 0.0) { +#ifndef FOR_R_TDA + std::cerr << "Relative error " << delta << ", must be > 0" << std::endl; +#endif + throw std::runtime_error("Bad delta in Wasserstein " + std::to_string(delta)); + } + if (_initialEpsilon < 0.0) { +#ifndef FOR_R_TDA + std::cerr << "Initial epsilon = " << _initialEpsilon << ", must be non-negative" << std::endl; +#endif + throw std::runtime_error("Bad initial epsilon in Wasserstein" + std::to_string(_initialEpsilon)); + } + if (_epsFactor < 0.0) { +#ifndef FOR_R_TDA + std::cerr << "Epsilon factor = " << _epsFactor << ", must be non-negative" << std::endl; +#endif + throw std::runtime_error("Bad epsilon factor in Wasserstein " + std::to_string(_epsFactor)); + } +#ifdef GAUSS_SEIDEL_AUCTION + AuctionRunnerGS auction(A, B, q, delta, _internal_p, _initialEpsilon, _epsFactor); +#else + AuctionRunnerJac auction(A, B, q, delta, _internal_p); +#endif + return auction.getWassersteinCost(); +} + bool readDiagramPointSet(const std::string& fname, std::vector>& result) { return readDiagramPointSet(fname.c_str(), result); @@ -84,7 +133,9 @@ bool readDiagramPointSet(const char* fname, std::vector> x >> y)) { +#ifndef FOR_R_TDA std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl; +#endif return false; } result.push_back(std::make_pair(x,y)); -- cgit v1.2.3