diff options
Diffstat (limited to 'geom_bottleneck/bottleneck/src/bottleneck.cpp')
-rw-r--r-- | geom_bottleneck/bottleneck/src/bottleneck.cpp | 51 |
1 files changed, 46 insertions, 5 deletions
diff --git a/geom_bottleneck/bottleneck/src/bottleneck.cpp b/geom_bottleneck/bottleneck/src/bottleneck.cpp index 34f7ed4..92e9dbf 100644 --- a/geom_bottleneck/bottleneck/src/bottleneck.cpp +++ b/geom_bottleneck/bottleneck/src/bottleneck.cpp @@ -100,7 +100,7 @@ double bottleneckDistApprox(DiagramPointSet& A, DiagramPointSet& B, const double } -double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B, std::vector<double>& pairwiseDist) +double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B, std::vector<double>& pairwiseDist, const int decPrecision) { //for(size_t k = 0; k < pairwiseDist.size(); ++k) { //std::cout << "pairwiseDist[" << k << "] = " << std::setprecision(15) << pairwiseDist[k] << std::endl; @@ -111,13 +111,19 @@ double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B bool useRangeSearch = true; double distEpsilon = std::numeric_limits<double>::max(); + double diffThreshold = 0.1; + for(int k = 0; k < decPrecision; ++k) { + diffThreshold /= 10.0; + } for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) { auto diff = pairwiseDist[k+1]- pairwiseDist[k]; - if ( diff > 1.0e-14 and diff < distEpsilon ) { + //std::cout << "diff = " << diff << ", pairwiseDist[k] = " << pairwiseDist[k] << std::endl; + if ( diff > diffThreshold and diff < distEpsilon ) { distEpsilon = diff; } } distEpsilon /= 3.0; + //std::cout << "decPrecision = " << decPrecision << ", distEpsilon = " << distEpsilon << std::endl; BoundMatchOracle oracle(A, B, distEpsilon, useRangeSearch); // binary search @@ -145,6 +151,11 @@ double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) { + return bottleneckDistExact(A, B, 14); +} + +double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B, const int decPrecision) +{ constexpr double epsilon = 0.001; auto interval = bottleneckDistApproxInterval(A, B, epsilon); const double delta = 0.5 * (interval.second - interval.first); @@ -291,7 +302,7 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B) //std::cout << std::setprecision(15) << ddd << std::endl; //} - return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist); + return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist, decPrecision); } double bottleneckDistSlow(DiagramPointSet& A, DiagramPointSet& B) @@ -495,13 +506,27 @@ double bottleneckDistSlow(DiagramPointSet& A, DiagramPointSet& B) */ } +// wrappers bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result) { - return readDiagramPointSet(fname.c_str(), result); + int decPrecision; + return readDiagramPointSet(fname.c_str(), result, decPrecision); } bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result) { + int decPrecision; + return readDiagramPointSet(fname, result, decPrecision); +} + +bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result, int& decPrecision) +{ + return readDiagramPointSet(fname.c_str(), result, decPrecision); +} + +// reading function +bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result, int& decPrecision) +{ size_t lineNumber { 0 }; result.clear(); std::ifstream f(fname); @@ -530,6 +555,22 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double continue; } line = std::string(whiteSpaceFront,whiteSpaceBack); + bool fracPart = false; + int currDecPrecision = 0; + for(auto c : line) { + if (c == '.') { + fracPart = true; + } else if (fracPart) { + if (isdigit(c)) { + currDecPrecision++; + } else { + fracPart = false; + if (currDecPrecision > decPrecision) + decPrecision = currDecPrecision; + currDecPrecision = 0; + } + } + } double x, y; std::istringstream iss(line); if (not(iss >> x >> y)) { @@ -542,7 +583,7 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double result.push_back(std::make_pair(x,y)); } else { #ifndef FOR_R_TDA - std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point: \"" << line << "\"" << std::endl; + //std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point ignored: \"" << line << "\"" << std::endl; #endif } } |