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.

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 <cstdlib>			// standard lib includes
 #include <cmath>			// math includes
-#include <iostream>			// I/O streams
 #include <cstring>			// C-style strings
 #include <vector>			 
 #include <assert.h>			 
@@ -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 <iostream>			// 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 <ANN/ANNx.h>					// 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 <utility>                      // for std::pair
 #include <ANN/ANNx.h>					// 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 <vector>
+#include <stdexcept>
 #include <math.h>
 #include <cstddef>
 #include <unordered_map>
 #include <unordered_set>
-#include <iostream>
 #include <string>
 #include <assert.h>
 
 #include "def_debug_bt.h"
 
+#ifndef FOR_R_TDA
+#include <iostream>
+#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<DiagramPoint, DiagramPointHash>::iterator end() { return points.end(); }
     std::unordered_set<DiagramPoint, DiagramPointHash>::const_iterator cbegin() const { return points.cbegin(); }
     std::unordered_set<DiagramPoint, DiagramPointHash>::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<class PairIterator> DiagramPointSet(PairIterator first, PairIterator last);
     template<class PairIterator> 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 <iostream>
+//#include <iostream>
 #include <fstream>
 #include <vector>
 #include <algorithm>
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 <ciso646>                      // make VS more conformal
 #endif
 
+#include <stdexcept>
 #include <cstdlib>						// C standard lib defs
 #include <ANN/ANNx.h>					// all ANN includes
 #include <ANN/ANNperf.h>				// 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 <ANN/ANNperf.h>				// 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 <ANN/ANNperf.h>				// 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 <algorithm>
 #include <cfloat>
+#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<double> 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<double> 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<std::pair<double, double
     result.clear();
     std::ifstream f(fname);
     if (!f.good()) {
+#ifndef FOR_R_TDA
         std::cerr << "Cannot open file " << fname << std::endl;
+#endif
         return false;
     }
     std::string line;
@@ -541,7 +533,9 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double
         double x, y;
         std::istringstream iss(line);
         if (not(iss >> 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 <http://www.gnu.org/licenses/>.
 
 */
 
-#include <iostream>
 #include <assert.h>
+#include "def_debug_bt.h"
 #include "bound_match.h"
+#ifndef FOR_R_TDA
+#include <iostream>
+#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<AuctionOracle> oracle;
 #ifdef KEEP_UNASSIGNED_ORDERED
diff --git a/geom_matching/wasserstein/include/auction_runner_jac.h b/geom_matching/wasserstein/include/auction_runner_jac.h
index 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<DiagramPoint>& A, const std::vector<DiagramPoint>& B, const double q,  const double _delta, const double _internal_p);
+    AuctionRunnerJac(const std::vector<DiagramPoint>& A, const std::vector<DiagramPoint>& 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<IdxValPair> bidTable;
     // to get the 2 best items
     std::unique_ptr<AuctionOracle> 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 <cstddef>
 #include <unordered_map>
 #include <unordered_set>
-#include <iostream>
 #include <string>
 #include <assert.h>
 
@@ -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 <boost/serialization/access.hpp>
 #include <boost/serialization/base_object.hpp>
 
-#include <iostream>
+//#include <iostream>
 #include <fstream>
 #include <string>
 #include <sstream>
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 <stack>
 
 #include "../parallel/tbb.h"
+#include "def_debug_ws.h"
 
 template<class T>
 dnn::KDTree<T>::
@@ -127,13 +128,8 @@ search(PointHandle q, ResultsFunctor& rf) const
     // TODO: use tbb::scalable_allocator for the queue
     std::queue<KDTreeNode>  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<class T>
@@ -290,6 +285,7 @@ void
 dnn::KDTree<T>::
 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 <iostream>
+//#include <iostream>
 #include <vector>
 
 #include <boost/range.hpp>
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<DiagramPoint>& 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<DiagramPoint>& A, 
+                          const std::vector<DiagramPoint>& B, 
+                          const double q, 
+                          const double delta,
+                          const double _internal_p = std::numeric_limits<double>::infinity(),
+                          const double _initialEpsilon = 0.0,
+                          const double _epsFactor = 0.0);
+
 
 // compare as multisets
 template<class PairContainer>
@@ -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<class PairContainer>
+double wassersteinCost(PairContainer& A, PairContainer& B, const double q, const double delta, const double _internal_p = std::numeric_limits<double>::infinity(), const double _initialEpsilon = 0.0, const double _epsFactor = 0.0)
+{
+    if (areEqual(A, B)) {
+        return 0.0;
+    }
+
+    std::vector<DiagramPoint> 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<DiagramPoint>& 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<double> weightVec;
@@ -134,8 +132,10 @@ 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,
     // do not update queues.
@@ -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<double> weightVec;
@@ -319,8 +321,10 @@ 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() ) {
         // lazy: record the moment we updated the price of the items,
@@ -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<DiagramPoint>& _bidde
     }
     DnnTraits traits;
     traits.internal_p = internal_p;
-    //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl;
     kdtree = new dnn::KDTree<DnnTraits>(traits, dnnPointHandles, wassersteinPower);
 
     size_t dnnItemIdxAll { 0 };
@@ -658,7 +660,6 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector<DiagramPoint>& _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<DiagramPoint>& _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<DnnTraits>(traits, dnnPointHandlesAll, wassersteinPower);
     
     size_t handleIdx {0};
@@ -677,13 +677,9 @@ AuctionOracleKDTree::AuctionOracleKDTree(const std::vector<DiagramPoint>& _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<IdxValPair> 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<D
         dnnPointHandles.push_back(&dnnPoints[i]);
     }
     DnnTraits traits;
-    //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl;
     kdtree = new dnn::KDTree<DnnTraits>(traits, dnnPointHandles, wassersteinPower);
     
     size_t handleIdx {0};
@@ -997,12 +980,9 @@ AuctionOracleKDTreeRestricted::AuctionOracleKDTreeRestricted(const std::vector<D
         }
     }
     //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 << "AuctionOracleKDTreeRestricted: weightAdjConst = " << weightAdjConst << std::endl;
-    //std::cout << "AuctionOracleKDTreeRestricted constructor done" << std::endl;
 }
 
 DebugOptimalBid AuctionOracleKDTreeRestricted::getOptimalBidDebug(IdxType bidderIdx)
@@ -1010,9 +990,6 @@ DebugOptimalBid AuctionOracleKDTreeRestricted::getOptimalBidDebug(IdxType bidder
     DebugOptimalBid result;
     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.
@@ -1064,7 +1041,6 @@ DebugOptimalBid AuctionOracleKDTreeRestricted::getOptimalBidDebug(IdxType bidder
         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 };
         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 <assert.h>
+#include <stdexcept>
 #include <algorithm>
 #include <functional>
 #include <iterator>
@@ -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,12 +149,14 @@ 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;
             notDone = true;
@@ -157,22 +164,27 @@ 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 );
         }
         // 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<int>(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<int>(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<DiagramPoint>& A, const std::vector<DiagramPoint>& B, const double q, const double _delta, const double _internal_p) :
+AuctionRunnerJac::AuctionRunnerJac(const std::vector<DiagramPoint>& A, const std::vector<DiagramPoint>& 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<DiagramPoint>& A, const std
     oracle = std::unique_ptr<AuctionOracle>(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<IdxType>(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<AuctionOracleKDTree*>(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 <set>
 #include "basic_defs_ws.h"
 
+#ifndef FOR_R_TDA
+#include <iostream>
+#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<DiagramPoint>& 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<DiagramPoint>& A, 
+                          const std::vector<DiagramPoint>& 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<std::pair<double, double>>& result)
 {
     return readDiagramPointSet(fname.c_str(), result);
@@ -84,7 +133,9 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double
     result.clear();
     std::ifstream f(fname);
     if (!f.good()) {
+#ifndef FOR_R_TDA
         std::cerr << "Cannot open file " << fname << std::endl;
+#endif
         return false;
     }
     std::string line;
@@ -109,7 +160,9 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double
         double x, y;
         std::istringstream iss(line);
         if (not(iss >> 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));