/*
Copyrigth 2015, D. Morozov, M. Kerber, A. Nigmetov
This file is part of GeomBottleneck.
GeomBottleneck is free software: you can redistribute it and/or modify
it under the terms of the Lesser GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
GeomBottleneck is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Lesser GNU General Public License for more details.
You should have received a copy of the Lesser GNU General Public License
along with GeomBottleneck. If not, see .
*/
#include
#include
#include
#include
#include "bottleneck.h"
//#include "test_dist_calc.h"
namespace geom_bt {
// return the interval (distMin, distMax) such that:
// a) actual bottleneck distance between A and B is contained in the interval
// b) if the interval is not (0,0), then (distMax - distMin) / distMin < epsilon
std::pair bottleneckDistApproxInterval(DiagramPointSet& A, DiagramPointSet& B, const double epsilon)
{
// empty diagrams are not considered as error
if (A.empty() and B.empty())
return std::make_pair(0.0, 0.0);
// link diagrams A and B by adding projections
addProjections(A, B);
// TODO: think about that!
// we need one threshold for checking if the distance is 0,
// another one for the oracle!
constexpr double epsThreshold { 1.0e-10 };
std::pair result { 0.0, 0.0 };
bool useRangeSearch { true };
// construct an oracle
BoundMatchOracle oracle(A, B, epsThreshold, useRangeSearch);
// check for distance = 0
if (oracle.isMatchLess(2*epsThreshold)) {
return result;
}
// get a 3-approximation of maximal distance between A and B
// as a starting value for probe distance
double distProbe { getFurthestDistance3Approx(A, B) };
// aliases for result components
double& distMin {result.first};
double& distMax {result.second};
if ( oracle.isMatchLess(distProbe) ) {
// distProbe is an upper bound,
// find lower bound with binary search
do {
distMax = distProbe;
distProbe /= 2.0;
} while (oracle.isMatchLess(distProbe));
distMin = distProbe;
} else {
// distProbe is a lower bound,
// find upper bound with exponential search
do {
distMin = distProbe;
distProbe *= 2.0;
} while (!oracle.isMatchLess(distProbe));
distMax = distProbe;
}
// bounds are found, perform binary search
//std::cout << "Bounds found, distMin = " << distMin << ", distMax = " << distMax << ", ratio = " << ( distMax - distMin ) / distMin << std::endl ;
distProbe = ( distMin + distMax ) / 2.0;
while ( ( distMax - distMin ) / distMin >= epsilon ) {
if (oracle.isMatchLess(distProbe)) {
distMax = distProbe;
} else {
distMin = distProbe;
}
distProbe = ( distMin + distMax ) / 2.0;
}
return result;
}
// get approximate distance,
// see bottleneckDistApproxInterval
double bottleneckDistApprox(DiagramPointSet& A, DiagramPointSet& B, const double epsilon)
{
auto interval = bottleneckDistApproxInterval(A, B, epsilon);
return interval.second;
}
double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B, std::vector& pairwiseDist)
{
//for(size_t k = 0; k < pairwiseDist.size(); ++k) {
//std::cout << "pairwiseDist[" << k << "] = " << std::setprecision(15) << pairwiseDist[k] << std::endl;
//}
// trivial case: we have only one candidate
if (pairwiseDist.size() == 1)
return pairwiseDist[0];
bool useRangeSearch = true;
double distEpsilon = std::numeric_limits::max();
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 ) {
distEpsilon = diff;
}
}
distEpsilon /= 3.0;
BoundMatchOracle oracle(A, B, distEpsilon, useRangeSearch);
// binary search
size_t iterNum {0};
size_t idxMin {0}, idxMax {pairwiseDist.size() - 1};
size_t idxMid;
while(idxMax > idxMin) {
idxMid = static_cast(floor(idxMin + idxMax) / 2.0);
//std::cout << "while begin: min = " << idxMin << ", idxMax = " << idxMax << ", idxMid = " << idxMid << ", testing d = " << std::setprecision(15) << pairwiseDist[idxMid] << std::endl;
iterNum++;
// not A[imid] < dist <=> A[imid] >= dist <=> A[imid[ >= dist + eps
if (oracle.isMatchLess(pairwiseDist[idxMid] + distEpsilon / 2.0)) {
//std::cout << "isMatchLess = true" << std::endl;
idxMax = idxMid;
} else {
//std::cout << "isMatchLess = false " << std::endl;
idxMin = idxMid + 1;
}
//std::cout << "while end: idxMin = " << idxMin << ", idxMax = " << idxMax << ", idxMid = " << idxMid << std::endl;
}
idxMid = static_cast(floor(idxMin + idxMax) / 2.0);
return pairwiseDist[idxMid];
}
double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B)
{
constexpr double epsilon = 0.001;
auto interval = bottleneckDistApproxInterval(A, B, epsilon);
const double delta = 0.5 * (interval.second - interval.first);
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;
if ( delta == 0 ) {
return interval.first;
}
// copy points from A to a vector
// todo: get rid of this?
std::vector pointsA;
pointsA.reserve(A.size());
for(const auto& ptA : A) {
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::sort(killDist.begin(), killDist.end());
//for(auto d : killDist) {
//std::cout << d << std::endl;
//}
//std::cout << "*************" << std::endl;
// in this vector we store the distances between the points
// that are candidates to realize
std::vector pairwiseDist;
{
// vector to store centers of vertical stripes
// two for each point in A and the id of the corresponding point
std::vector> xCentersVec;
xCentersVec.reserve(2 * pointsA.size());
for(auto ptA : pointsA) {
xCentersVec.push_back(std::make_pair(ptA.getRealX() - approxDist, ptA));
xCentersVec.push_back(std::make_pair(ptA.getRealX() + approxDist, ptA));
}
// lambda to compare pairs w.r.t coordinate
auto compLambda = [](std::pair a, std::pair b)
{ return a.first < b.first; };
std::sort(xCentersVec.begin(), xCentersVec.end(), compLambda);
//std::cout << "xCentersVec.size = " << xCentersVec.size() << std::endl;
//for(auto p = xCentersVec.begin(); p!= xCentersVec.end(); ++p) {
//if (p->second.id == 200) {
//std::cout << "index of 200: " << p - xCentersVec.begin() << std::endl;
//}
//}
//std::vector
// todo: sort points in B, reduce search range in lower and upper bounds
for(auto ptB : B) {
// iterator to the first stripe such that ptB lies to the left
// from its right boundary (x_B <= x_j + \delta iff x_j >= x_B - \delta
auto itStart = std::lower_bound(xCentersVec.begin(),
xCentersVec.end(),
std::make_pair(ptB.getRealX() - delta, ptB),
compLambda);
//if (ptB.id == 236) {
//std::cout << itStart - xCentersVec.begin() << std::endl;
//}
for(auto iterA = itStart; iterA < xCentersVec.end(); ++iterA) {
//if (ptB.id == 236) {
//std::cout << "consider " << iterA->second << std::endl;
//}
if ( ptB.getRealX() < iterA->first - delta) {
// from that moment x_B >= x_j - delta
// is violated: x_B no longer lies to right from the left
// boundary of current stripe
//if (ptB.id == 236) {
//std::cout << "break" << std::endl;
//}
break;
}
// we're here => ptB lies in vertical stripe,
// check if distance fits into the interval we've found
double pwDist = distLInf(iterA->second, ptB);
//if (ptB.id == 236) {
//std::cout << pwDist << std::endl;
//}
//std::cout << 1000*minDist << " <= " << 1000*pwDist << " <= " << 1000*maxDist << std::endl;
if (pwDist >= minDist and pwDist <= maxDist) {
pairwiseDist.push_back(pwDist);
}
}
}
}
{
// for y
// vector to store centers of vertical stripes
// two for each point in A and the id of the corresponding point
std::vector> yCentersVec;
yCentersVec.reserve(2 * pointsA.size());
for(auto ptA : pointsA) {
yCentersVec.push_back(std::make_pair(ptA.getRealY() - approxDist, ptA));
yCentersVec.push_back(std::make_pair(ptA.getRealY() + approxDist, ptA));
}
// lambda to compare pairs w.r.t coordinate
auto compLambda = [](std::pair a, std::pair b)
{ return a.first < b.first; };
std::sort(yCentersVec.begin(), yCentersVec.end(), compLambda);
// std::cout << "Sorted vector of y-centers:" << std::endl;
//for(auto coordPtPair : yCentersVec) {
//std::cout << coordPtPair.first << ", id = " << coordPtPair.second.id << std::endl;
//}
/*std::cout << "End of sorted vector of y-centers:" << std::endl;*/
//std::vector
// todo: sort points in B, reduce search range in lower and upper bounds
for(auto ptB : B) {
auto itStart = std::lower_bound(yCentersVec.begin(),
yCentersVec.end(),
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) {
break;
}
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::sort(pairwiseDist.begin(), pairwiseDist.end());
//for(auto ddd : pairwiseDist) {
//std::cerr << std::setprecision(15) << ddd << std::endl;
//}
return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist);
}
double bottleneckDistSlow(DiagramPointSet& A, DiagramPointSet& B)
{
// use range search when building the layer graph
bool useRangeSearch { true };
// find maximum of min. distances for each point,
// use this value as lower bound for bottleneck distance
bool useHeurMinIdx { true };
// find matching in a greedy manner to
// get an upper bound for a bottleneck distance
bool useHeurGreedyMatching { false };
// use successive multiplication of idxMin with 2 to get idxMax
bool goUpToFindIdxMax { false };
//
goUpToFindIdxMax = goUpToFindIdxMax and !useHeurGreedyMatching;
if (!useHeurGreedyMatching) {
long int N = 3 * (A.size() / 2 ) * (B.size() / 2);
std::vector pairwiseDist;
pairwiseDist.reserve(N);
double maxMinDist {0.0};
for(auto& p_A : A) {
double minDist { std::numeric_limits::max() };
for(auto& p_B : B) {
if (p_A.type != DiagramPoint::DIAG or p_B.type != DiagramPoint::DIAG) {
double d = distLInf(p_A, p_B);
pairwiseDist.push_back(d);
if (useHeurMinIdx and p_A.type != DiagramPoint::DIAG) {
if (d < minDist)
minDist = d;
}
}
}
if (useHeurMinIdx and DiagramPoint::DIAG != p_A.type and minDist > maxMinDist) {
maxMinDist = minDist;
}
}
std::sort(pairwiseDist.begin(), pairwiseDist.end());
double distEpsilon = std::numeric_limits::max();
for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) {
auto diff = pairwiseDist[k+1]- pairwiseDist[k];
if ( diff > 1.0e-10 and diff < distEpsilon ) {
distEpsilon = diff;
}
}
distEpsilon /= 3.0;
BoundMatchOracle oracle(A, B, distEpsilon, useRangeSearch);
// binary search
size_t iterNum {0};
size_t idxMin {0}, idxMax {pairwiseDist.size() - 1};
if (useHeurMinIdx) {
auto maxMinIter = std::equal_range(pairwiseDist.begin(), pairwiseDist.end(), maxMinDist);
assert(maxMinIter.first != pairwiseDist.end());
idxMin = maxMinIter.first - pairwiseDist.begin();
//std::cout << "maxMinDist = " << maxMinDist << ", idxMin = " << idxMin << ", d = " << pairwiseDist[idxMin] << std::endl;
}
if (goUpToFindIdxMax) {
if ( pairwiseDist.size() == 1) {
return pairwiseDist[0];
}
idxMax = std::max(idxMin, 1);
while (!oracle.isMatchLess(pairwiseDist[idxMax])) {
//std::cout << "entered while" << std::endl;
idxMin = idxMax;
if (2*idxMax > pairwiseDist.size() -1) {
idxMax = pairwiseDist.size() - 1;
break;
} else {
idxMax *= 2;
}
}
//std::cout << "size = " << pairwiseDist.size() << ", idxMax = " << idxMax << ", pw[max] = " << pairwiseDist[idxMax] << std::endl;
}
size_t idxMid { (idxMin + idxMax) / 2 };
while(idxMax > idxMin) {
iterNum++;
if (oracle.isMatchLess(pairwiseDist[idxMid])) {
idxMax = idxMid;
} else {
if (idxMax - idxMin == 1)
idxMin++;
else
idxMin = idxMid;
}
idxMid = (idxMin + idxMax) / 2;
}
return pairwiseDist[idxMid];
} else {
// with greeedy matching
long int N = A.size() * B.size();
std::vector pairwiseDist;
pairwiseDist.reserve(N);
double maxMinDist {0.0};
size_t idxA{0}, idxB{0};
for(auto p_A : A) {
double minDist { std::numeric_limits::max() };
idxB = 0;
for(auto p_B : B) {
double d = distLInf(p_A, p_B);
pairwiseDist.push_back( std::make_pair(d, std::make_pair(idxA, idxB) ) );
if (useHeurMinIdx and p_A.type != DiagramPoint::DIAG) {
if (d < minDist)
minDist = d;
}
idxB++;
}
if (useHeurMinIdx and DiagramPoint::DIAG != p_A.type and minDist > maxMinDist) {
maxMinDist = minDist;
}
idxA++;
}
auto compLambda = [](DistVerticesPair a, DistVerticesPair b)
{ return a.first < b.first;};
std::sort(pairwiseDist.begin(),
pairwiseDist.end(),
compLambda);
double distEpsilon = std::numeric_limits::max();
for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) {
auto diff = pairwiseDist[k+1].first - pairwiseDist[k].first;
if ( diff > 1.0e-10 and diff < distEpsilon ) {
distEpsilon = diff;
}
}
distEpsilon /= 3.0;
BoundMatchOracle oracle(A, B, distEpsilon, useRangeSearch);
// construct greedy matching
size_t numVert { A.size() };
size_t numMatched { 0 };
std::unordered_set aTobMatched, bToaMatched;
aTobMatched.reserve(numVert);
bToaMatched.reserve(numVert);
size_t distVecIdx {0};
while( numMatched < numVert) {
auto vertPair = pairwiseDist[distVecIdx++].second;
//std::cout << "distVecIdx = " << distVecIdx << ", matched: " << numMatched << " out of " << numVert << std::endl;
//std::cout << "vertex A idx = " << vertPair.first << ", B idx: " << vertPair.second << " out of " << numVert << std::endl;
if ( aTobMatched.count(vertPair.first) == 0 and
bToaMatched.count(vertPair.second) == 0 ) {
aTobMatched.insert(vertPair.first);
bToaMatched.insert(vertPair.second);
numMatched++;
}
}
size_t idxMax = distVecIdx-1;
//std::cout << "idxMax = " << idxMax << ", size = " << pairwiseDist.size() << std::endl;
// binary search
size_t iterNum {0};
size_t idxMin {0};
if (useHeurMinIdx) {
auto maxMinIter = std::equal_range(pairwiseDist.begin(),
pairwiseDist.end(),
std::make_pair(maxMinDist, std::make_pair(0,0)),
compLambda);
assert(maxMinIter.first != pairwiseDist.end());
idxMin = maxMinIter.first - pairwiseDist.begin();
//std::cout << "maxMinDist = " << maxMinDist << ", idxMin = " << idxMin << ", d = " << pairwiseDist[idxMin].first << std::endl;
}
size_t idxMid { (idxMin + idxMax) / 2 };
while(idxMax > idxMin) {
iterNum++;
if (oracle.isMatchLess(pairwiseDist[idxMid].first)) {
idxMax = idxMid;
} else {
if (idxMax - idxMin == 1)
idxMin++;
else
idxMin = idxMid;
}
idxMid = (idxMin + idxMax) / 2;
}
return pairwiseDist[idxMid].first;
}
// stats
/*
// count number of edges
// pairwiseDist is sorted, add edges of the same length
int edgeNumber {idxMid};
while(pairwiseDist[edgeNumber + 1] == pairwiseDist[edgeNumber])
edgeNumber++;
// add edges between diagonal points
edgeNumber += N / 3;
// output stats
std::cout << idxMid << "\t" << N;
std::cout << "\t" << iterNum;
std::cout << "\t" << A.size() + B.size();
std::cout << "\t" << edgeNumber << "\t";
std::cout << (double)(edgeNumber) / (double)(A.size() + B.size()) << std::endl;
*/
}
bool readDiagramPointSet(const std::string& fname, std::vector>& result)
{
return readDiagramPointSet(fname.c_str(), result);
}
bool readDiagramPointSet(const char* fname, std::vector>& result)
{
size_t lineNumber { 0 };
result.clear();
std::ifstream f(fname);
if (!f.good()) {
std::cerr << "Cannot open file " << fname << std::endl;
return false;
}
std::string line;
while(std::getline(f, line)) {
lineNumber++;
// process comments: remove everything after hash
auto hashPos = line.find_first_of("#", 0);
if( std::string::npos != hashPos) {
line = std::string(line.begin(), line.begin() + hashPos);
}
if (line.empty()) {
continue;
}
// trim whitespaces
auto whiteSpaceFront = std::find_if_not(line.begin(),line.end(),isspace);
auto whiteSpaceBack = std::find_if_not(line.rbegin(),line.rend(),isspace).base();
if (whiteSpaceBack <= whiteSpaceFront) {
// line consists of spaces only - move to the next line
continue;
}
line = std::string(whiteSpaceFront,whiteSpaceBack);
double x, y;
std::istringstream iss(line);
if (not(iss >> x >> y)) {
std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl;
return false;
}
result.push_back(std::make_pair(x,y));
}
f.close();
return true;
}
}