diff --git a/core/base/constrainedGradientDescent/CMakeLists.txt b/core/base/constrainedGradientDescent/CMakeLists.txt new file mode 100644 index 0000000000..1717feff7f --- /dev/null +++ b/core/base/constrainedGradientDescent/CMakeLists.txt @@ -0,0 +1,14 @@ +ttk_add_base_library(constrainedGradientDescent + SOURCES + ConstrainedGradientDescent.cpp + HEADERS + ConstrainedGradientDescent.h + DEPENDS + common + persistenceDiagram +) + +if(TTK_ENABLE_EIGEN) + target_compile_definitions(constrainedGradientDescent PRIVATE TTK_ENABLE_EIGEN) + target_link_libraries(constrainedGradientDescent PRIVATE Eigen3::Eigen) +endif() \ No newline at end of file diff --git a/core/base/constrainedGradientDescent/ConstrainedGradientDescent.cpp b/core/base/constrainedGradientDescent/ConstrainedGradientDescent.cpp new file mode 100644 index 0000000000..74ea1d97b5 --- /dev/null +++ b/core/base/constrainedGradientDescent/ConstrainedGradientDescent.cpp @@ -0,0 +1,366 @@ +#include +#include +#include + +#ifdef TTK_ENABLE_EIGEN +#include +#include +#endif // TTK_ENABLE_EIGEN + +using namespace ttk; + +void ConstrainedGradientDescent::executeWeightsProjected( + std::vector &hessianList, + std::vector &weights, + const std::vector &grad, + bool maxEigenValue) { + gradientDescentWeights(hessianList, weights, grad, maxEigenValue); + projectionOnSimplex(weights); +} + +void ConstrainedGradientDescent::executeAtoms( + std::vector &DictDiagrams, + const std::vector> &matchings, + const ttk::DiagramType &Barycenter, + const std::vector &gradsLists, + const std::vector &checkerAtomsExt, + std::vector> &projForDiag, + ttk::DiagramType &featuresToAdd, + std::vector> &projLocations, + std::vector> &vectorForProjContrib, + std::vector>> &pairToAddGradList, + ttk::DiagramType &infoToAdd) { + gradientDescentAtoms(DictDiagrams, matchings, Barycenter, gradsLists, + checkerAtomsExt, projForDiag, featuresToAdd, + projLocations, vectorForProjContrib, pairToAddGradList, + infoToAdd); +} + +// simple projection on simplex, aka where a vector has positive elements and +// sum to 1. +void ConstrainedGradientDescent::projectionOnSimplex( + std::vector &weights) { + + int n = weights.size(); + std::vector copy_temp = weights; + std::sort(copy_temp.rbegin(), copy_temp.rend()); + double K = 1.; + double somme_u = copy_temp[0]; + double theta = (somme_u - 1.) / K; + while(K < n && (somme_u + copy_temp[K] - 1.) / (K + 1.) < copy_temp[K]) { + somme_u += copy_temp[K]; + K += 1.; + theta = (somme_u - 1.) / K; + } + for(int i = 0; i < n; ++i) { + weights[i] = std::max(weights[i] - theta, 0.); + } + + double sum = 0.; + for(int i = 0; i < n - 1; ++i) { + weights[i] = trunc(weights[i] * 1e8) / 1e8; + sum += weights[i]; + } + weights[n - 1] = 1. - sum; +} + +void ConstrainedGradientDescent::gradientDescentWeights( + std::vector &hessianList, + std::vector &weights, + const std::vector &grad, + bool maxEigenValue) { + + int n = weights.size(); + double stepWeight; + double L = 0.; +#ifndef TTK_ENABLE_EIGEN + maxEigenValue = false; +#endif // TTK_ENABLE_EIGEN + if(maxEigenValue) { +#ifdef TTK_ENABLE_EIGEN + for(size_t i = 0; i < hessianList.size(); ++i) { + auto &hessian = hessianList[i]; + int m = hessian.size(); + Eigen::MatrixXd H(m, m); + for(size_t j = 0; j < hessian.size(); ++j) { + for(size_t k = 0; k < hessian.size(); ++k) { + H(j, k) = hessian[j][k]; + } + } + Eigen::EigenSolver es; + es.compute(H, false); + Eigen::VectorXcd eigvals = es.eigenvalues(); + L += eigvals.lpNorm(); + } +#endif // TTK_ENABLE_EIGEN + } else { + for(size_t i = 0; i < hessianList.size(); ++i) { + auto &hessian = hessianList[i]; + for(size_t k = 0; k < hessian.size(); ++k) { + double diag = hessian[k][k]; + L += 1. * diag; + } + } + } + + if(L > 0) { + stepWeight = 1. / L; + } else { + stepWeight = 0; + } + + for(int i = 0; i < n; ++i) { + weights[i] = weights[i] - stepWeight * grad[i]; + } +} + +// TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO +// METTRE TIMER POUR VOIR QUOI PARALELLISER +// TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO + +void ConstrainedGradientDescent::gradientDescentAtoms( + std::vector &DictDiagrams, + const std::vector> &matchings, + const ttk::DiagramType &Barycenter, + const std::vector &gradsLists, + const std::vector &checkerAtomsExt, + std::vector> &projForDiag, + ttk::DiagramType &featuresToAdd, + std::vector> &projLocations, + std::vector> &vectorForProjContrib, + std::vector>> &pairToAddGradList, + ttk::DiagramType &infoToAdd) { + + // Here vector of diagramTuple because it is not a persistence diagram per + // say. + // we get the right pairs to update for each barycenter pair. + + std::vector miniBirth(matchings.size()); + for(size_t i = 0; i < matchings.size(); ++i) { + auto &t = DictDiagrams[i][0]; + miniBirth[i] = t.birth.sfValue; + } + + std::vector>> gradBuffersList( + Barycenter.size()); + std::vector> projectionsBuffer(Barycenter.size()); + + for(size_t i = 0; i < Barycenter.size(); ++i) { + projectionsBuffer[i].resize(matchings.size()); + } + + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + gradBuffersList[i].resize(matchings.size()); + } + + std::vector> checker(Barycenter.size()); + for(size_t i = 0; i < checker.size(); ++i) { + checker[i].resize(matchings.size()); + } + std::vector tracker(Barycenter.size(), 0); + std::vector> trackerDiagonal(Barycenter.size()); + std::vector> trackerMatch(Barycenter.size()); + for(size_t i = 0; i < trackerDiagonal.size(); ++i) { + trackerDiagonal[i].resize(matchings.size()); + } + for(size_t i = 0; i < trackerMatch.size(); ++i) { + trackerMatch[i].resize(matchings.size()); + } + + for(size_t i = 0; i < matchings.size(); ++i) { + for(size_t j = 0; j < matchings[i].size(); ++j) { + const auto &t = matchings[i][j]; + // Id in atom + const SimplexId Id1 = std::get<0>(t); + // Id in barycenter + const SimplexId Id2 = std::get<1>(t); + if(Id2 < 0 || static_cast(gradBuffersList.size()) <= Id2 + || static_cast(DictDiagrams[i].size()) <= Id1) { + continue; + } else { + if(Id1 < 0) { + const auto &t3 = Barycenter[Id2]; + auto &point = gradBuffersList[Id2][i]; + const double birthBarycenter = t3.birth.sfValue; + const double deathBarycenter = t3.death.sfValue; + const double birthDeathAtom + = birthBarycenter + (deathBarycenter - birthBarycenter) / 2.; + point[0] = birthDeathAtom; + point[1] = birthDeathAtom; + checker[Id2][i] = i; + if(checker[Id2].size() > matchings.size()) { + std::raise(SIGINT); + } + tracker[Id2] = 1; + trackerMatch[Id2][i] = Id1; + if(static_cast(DictDiagrams[i].size()) <= Id1) { + std::cout << "ID1: " << Id1 << std::endl; + } + trackerDiagonal[Id2][i] = 1; + projectionsBuffer[Id2][i] = birthDeathAtom; + + } else { + const auto &t2 = DictDiagrams[i][Id1]; + auto &point = gradBuffersList[Id2][i]; + const double birthAtom = t2.birth.sfValue; + const double deathAtom = t2.death.sfValue; + point[0] = birthAtom; + point[1] = deathAtom; + checker[Id2][i] = i; + if(checker[Id2].size() > matchings.size()) { + std::raise(SIGINT); + } + tracker[Id2] = 1; + trackerMatch[Id2][i] = Id1; + trackerDiagonal[Id2][i] = 0; + projectionsBuffer[Id2][i] = 0.; + } + } + } + } + + gradBuffersList.insert( + gradBuffersList.end(), pairToAddGradList.begin(), pairToAddGradList.end()); + for(size_t j = 0; j < pairToAddGradList.size(); ++j) { + std::vector temp1(matchings.size(), 1); + std::vector temp2(matchings.size(), -1); + std::vector temp3(matchings.size()); + for(size_t l = 0; l < matchings.size(); ++l) { + temp3[l] = static_cast(l); + } + trackerDiagonal.push_back(temp1); + trackerMatch.push_back(temp1); + checker.push_back(temp3); + tracker.push_back(1); + } + + // #ifdef TTK_ENABLE_OPENMP + // #pragma omp parallel for num_threads(threadNumber_) + // #endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + if(tracker[i] == 0 || checkerAtomsExt[i] == 0) { + + continue; + } else { + + std::vector pos(gradBuffersList[i].size(), 0.); + for(size_t j = 0; j < checker[i].size(); ++j) { + auto &t = gradBuffersList[i][checker[i][j]]; + double birth = t[0]; + double death = t[1]; + pos[j] = death - birth; + } + std::vector pos2(pos.size(), false); + std::vector temp2; + for(size_t p = 0; p < checker[i].size(); ++p) { + auto &t = gradBuffersList[i][checker[i][p]]; + double birth = t[0]; + pos2[p] = birth == 0.; + if(birth > 0.) { + temp2.push_back(birth); + } + } + for(size_t p = 0; p < checker[i].size(); ++p) { + if(pos2[p]) { + + auto &t = gradBuffersList[i][checker[i][p]]; + + t[1] = t[1] - (this->stepAtom) * gradsLists[i][checker[i][p]][1]; + } else if(pos[p] < 1e-7) { + continue; + } else { + + auto &t0 = gradBuffersList[i][checker[i][p]]; + + t0[0] = t0[0] - (this->stepAtom) * gradsLists[i][checker[i][p]][0]; + t0[1] = t0[1] - (this->stepAtom) * gradsLists[i][checker[i][p]][1]; + + if(t0[0] > t0[1]) { + t0[1] = t0[0]; + } + } + } + } + } + + for(size_t i = 0; i < checker.size(); ++i) { + if(tracker[i] == 0 || checkerAtomsExt[i] == 0) { + continue; + } else { + if(i < Barycenter.size()) { + for(size_t j = 0; j < checker[i].size(); ++j) { + auto &trackerTemp = trackerMatch[i][j]; + if(trackerDiagonal[i][j] == 1 || trackerTemp == -1 + || trackerTemp > 10000) { + std::vector projAndIndex(DictDiagrams.size() + 1); + for(size_t m = 0; m < DictDiagrams.size(); ++m) { + projAndIndex[m] = trackerMatch[i][m]; + } + + } else { + auto &index = checker[i][j]; + auto &t2 = gradBuffersList[i][index]; + + auto &t1 = DictDiagrams[index][trackerTemp]; + + if(t2[0] < miniBirth[index]) { + t1.birth.sfValue = miniBirth[index]; + } else { + t1.birth.sfValue = t2[0]; + } + if(t2[1] < t2[0]) { + t1.birth.sfValue = t2[0]; + t1.death.sfValue = t2[0]; + } else { + t1.death.sfValue = t2[1]; + } + } + } + } else { + const auto &infos = infoToAdd[static_cast(i) + - static_cast(Barycenter.size())]; + + for(size_t j = 0; j < checker[i].size(); ++j) { + auto &trackerTemp = trackerMatch[i][j]; + if(trackerDiagonal[i][j] == 1 || trackerTemp == -1 + || trackerTemp > 10000) { + auto &index = checker[i][j]; + std::vector projAndIndex(DictDiagrams.size() + 1); + int atomIndex = static_cast(index); + for(size_t m = 0; m < DictDiagrams.size(); ++m) { + projAndIndex[m] = trackerMatch[i][m]; + } + projAndIndex[DictDiagrams.size()] = atomIndex; + projForDiag.push_back(projAndIndex); + featuresToAdd.push_back(infos); + projLocations.push_back(gradBuffersList[i][index]); + vectorForProjContrib.push_back(gradsLists[i][index]); + } else { + auto &index = checker[i][j]; + auto &t2 = gradBuffersList[i][index]; + auto &t1 = DictDiagrams[index][trackerTemp]; + if(t2[0] < miniBirth[index]) { + t1.birth.sfValue = miniBirth[index]; + } else { + t1.birth.sfValue = t2[0]; + } + if(t2[1] < t2[0]) { + t1.birth.sfValue = t2[0]; + t1.death.sfValue = t2[0]; + } else { + t1.death.sfValue = t2[1]; + } + } + } + } + } + } +} + +void ConstrainedGradientDescent::setStep(double factEquiv) { + this->stepAtom = 1. / (2. * 2. * factEquiv); +} + +void ConstrainedGradientDescent::reduceStep() { + this->stepAtom = this->stepAtom / 2.; +} \ No newline at end of file diff --git a/core/base/constrainedGradientDescent/ConstrainedGradientDescent.h b/core/base/constrainedGradientDescent/ConstrainedGradientDescent.h new file mode 100644 index 0000000000..cc94ffab3b --- /dev/null +++ b/core/base/constrainedGradientDescent/ConstrainedGradientDescent.h @@ -0,0 +1,70 @@ +#pragma once + +#include +#include + +#include +#include +#include + +namespace ttk { + using Matrix = std::vector>; + + class ConstrainedGradientDescent : public Debug { + + public: + ConstrainedGradientDescent() { + this->setDebugMsgPrefix("ConstrainedGradientDescent"); + }; + + void executeWeightsProjected(std::vector &hessianList, + std::vector &weights, + const std::vector &grad, + bool maxEigenValue); + + void executeAtoms( + std::vector &DictDiagrams, + const std::vector> &matchings, + const ttk::DiagramType &Barycenter, + const std::vector &gradsLists, + const std::vector &checkerAtomsExt, + std::vector> &projForDiag, + ttk::DiagramType &featuresToAdd, + std::vector> &projLocations, + std::vector> &vectorForProjContrib, + std::vector>> &pairToAddGradList, + ttk::DiagramType &infoToAdd); + + void setStep(double factEquiv); + void reduceStep(); + // void executeAtoms(std::vector &DictDiagrams); + + // inline void setNbAtoms(const int nbAtoms) { + // NbAtoms = nbAtoms; + //} + + protected: + void projectionOnSimplex(std::vector &weights); + + void gradientDescentWeights(std::vector &hessianList, + std::vector &weights, + const std::vector &grad, + bool maxEigenValue); + + void gradientDescentAtoms( + std::vector &DictDiagrams, + const std::vector> &matchings, + const ttk::DiagramType &Barycenter, + const std::vector &gradsLists, + const std::vector &checkerAtomsExt, + std::vector> &projForDiag, + ttk::DiagramType &featuresToAdd, + std::vector> &projLocations, + std::vector> &vectorForProjContrib, + std::vector>> &pairToAddGradList, + ttk::DiagramType &infoToAdd); + + double stepAtom; + }; + +} // namespace ttk diff --git a/core/base/initDictPersistenceDiagram/CMakeLists.txt b/core/base/initDictPersistenceDiagram/CMakeLists.txt new file mode 100644 index 0000000000..035f4f0831 --- /dev/null +++ b/core/base/initDictPersistenceDiagram/CMakeLists.txt @@ -0,0 +1,11 @@ +ttk_add_base_library(initDictPersistenceDiagram + SOURCES + InitDictPersistenceDiagram.cpp + InitDictRandomly.cpp + HEADERS + InitDictPersistenceDiagram.h + InitDictRandomly.h + DEPENDS + common + persistenceDiagramDistanceMatrix +) diff --git a/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.cpp b/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.cpp new file mode 100644 index 0000000000..56a75ac0e2 --- /dev/null +++ b/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.cpp @@ -0,0 +1,69 @@ +#include + +using namespace ttk; + +void InitDictPersistenceDiagram::execute( + std::vector &DictDiagrams, + const std::vector &datas, + const int nbAtoms, + bool do_min_, + bool do_sad_, + bool do_max_) { + + const size_t nDiags = datas.size(); + + PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nDiags, 0}; + MatrixCalculator.setDos(do_min_, do_sad_, do_max_); + MatrixCalculator.setThreadNumber(this->threadNumber_); + Matrix distMatrix = MatrixCalculator.execute(datas, nInputs); + std::vector allDistsSummed(nDiags); + for(size_t i = 0; i < nDiags; ++i) { + const auto &line = distMatrix[i]; + for(size_t j = 0; j < nDiags; ++j) { + allDistsSummed[j] += line[j]; + } + } + std::vector indices; + int Id1 = std::max_element(allDistsSummed.begin(), allDistsSummed.end()) + - allDistsSummed.begin(); + indices.push_back(Id1); + + for(int i = 1; i < nbAtoms; ++i) { + indices.push_back(getNextIndex(distMatrix, indices)); + } + + std::vector tempDistsSummed(nbAtoms); + for(int i = 0; i < nbAtoms; ++i) { + tempDistsSummed[i] = allDistsSummed[indices[i]]; + } + + int FirstId = std::min_element(tempDistsSummed.begin(), tempDistsSummed.end()) + - tempDistsSummed.begin(); + DictDiagrams.push_back(datas[indices[FirstId]]); + for(int i = 0; i < nbAtoms; ++i) { + if(i == FirstId) { + continue; + } else { + DictDiagrams.push_back(datas[indices[i]]); + } + } +} + +int InitDictPersistenceDiagram::getNextIndex( + const Matrix &distMatrix, const std::vector &indices) const { + std::vector allSumCumul(distMatrix.size(), 0.); + for(size_t k = 0; k < indices.size(); ++k) { + const auto &line = distMatrix[indices[k]]; + for(size_t i = 0; i < distMatrix.size(); ++i) { + if(std::find(indices.begin(), indices.end(), i) != indices.end()) { + allSumCumul[i] += 0.; + } else { + allSumCumul[i] += line[i]; + } + } + } + int newId = std::max_element(allSumCumul.begin(), allSumCumul.end()) + - allSumCumul.begin(); + return newId; +} diff --git a/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.h b/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.h new file mode 100644 index 0000000000..b8a13ef7b6 --- /dev/null +++ b/core/base/initDictPersistenceDiagram/InitDictPersistenceDiagram.h @@ -0,0 +1,41 @@ +#pragma once + +// #include +// #include +// #include +#include +#include +#include + +#include +#include + +namespace ttk { + using Matrix = std::vector>; + + class InitDictPersistenceDiagram : public Debug { + + public: + InitDictPersistenceDiagram() { + this->setDebugMsgPrefix("InitFarBorderDict"); + }; + + void execute(std::vector &DictDiagrams, + const std::vector &datas, + const int nbAtoms, + bool do_min_, + bool do_sad_, + bool do_max_); + + protected: + int getNextIndex(const Matrix &distMatrix, + const std::vector &indices) const; + + int Wasserstein{2}; + double Alpha{1.0}; + double DeltaLim{0.01}; + double Lambda{0}; + size_t MaxNumberOfPairs{20}; + double MinPersistence{0.1}; + }; +} // namespace ttk diff --git a/core/base/initDictPersistenceDiagram/InitDictRandomly.cpp b/core/base/initDictPersistenceDiagram/InitDictRandomly.cpp new file mode 100644 index 0000000000..b2d3b2cc77 --- /dev/null +++ b/core/base/initDictPersistenceDiagram/InitDictRandomly.cpp @@ -0,0 +1,24 @@ +#include +#include + +#include +#include + +using namespace ttk; + +void InitRandomDict::execute(std::vector &DictDiagrams, + const std::vector &datas, + const int nbAtom, + const int seed) { + int nDiags = datas.size(); + DictDiagrams.resize(nbAtom); + std::vector indices(nDiags); + std::iota(indices.begin(), indices.end(), 0); + std::mt19937 random_engine{}; + random_engine.seed(seed); + ttk::shuffle(indices, random_engine); + for(int i = 0; i < nbAtom; ++i) { + const auto &atom = datas[indices[i]]; + DictDiagrams[i] = atom; + } +} diff --git a/core/base/initDictPersistenceDiagram/InitDictRandomly.h b/core/base/initDictPersistenceDiagram/InitDictRandomly.h new file mode 100644 index 0000000000..59784cfd94 --- /dev/null +++ b/core/base/initDictPersistenceDiagram/InitDictRandomly.h @@ -0,0 +1,25 @@ +#pragma once + +#include +#include + +#include +#include +#include + +namespace ttk { + using Matrix = std::vector>; + + class InitRandomDict : public Debug { + + public: + InitRandomDict() { + this->setDebugMsgPrefix("InitRandomDict"); + }; + + void execute(std::vector &DictDiagrams, + const std::vector &datas, + const int nbAtoms, + const int seed); + }; +} // namespace ttk diff --git a/core/base/persistenceDiagramClustering/PDBarycenter.cpp b/core/base/persistenceDiagramClustering/PDBarycenter.cpp index bb4ff80985..198fc7c865 100644 --- a/core/base/persistenceDiagramClustering/PDBarycenter.cpp +++ b/core/base/persistenceDiagramClustering/PDBarycenter.cpp @@ -172,6 +172,9 @@ double ttk::PDBarycenter::updateBarycenter( points_deleted_ = 0; double max_shift = 0; + std::vector weight_diag_matchings( + n_goods); // Total weight for diagonal matchings for each point of the + // barycenter std::vector count_diag_matchings( n_goods); // Number of diagonal matchings for each point of the barycenter std::vector x(n_goods); @@ -181,6 +184,7 @@ double ttk::PDBarycenter::updateBarycenter( std::vector crit_coords_z(n_goods); for(size_t i = 0; i < n_goods; i++) { count_diag_matchings[i] = 0; + weight_diag_matchings[i] = 0.0; x[i] = 0; y[i] = 0; crit_coords_x[i] = 0; @@ -192,22 +196,31 @@ double ttk::PDBarycenter::updateBarycenter( min_prices[j] = std::numeric_limits::max(); } - std::vector - points_to_append; // Will collect bidders linked to diagonal + // Will collect bidders linked to diagonal second value + // indicate the corresponding weight. + std::vector> points_to_append; + // 2. Preprocess the matchings for(size_t j = 0; j < matchings.size(); j++) { + double weight = useCustomWeights_ ? (*customWeights_)[j] : 0.0; for(size_t i = 0; i < matchings[j].size(); i++) { int const bidder_id = std::get<0>(matchings[j][i]); int const good_id = std::get<1>(matchings[j][i]); if(good_id < 0 && bidder_id >= 0) { // Future new barycenter point - points_to_append.push_back(¤t_bidder_diagrams_[j].at(bidder_id)); + points_to_append.emplace_back( + ¤t_bidder_diagrams_[j].at(bidder_id), weight); } else if(good_id >= 0 && bidder_id >= 0) { // Update coordinates (to be divided by the number of diagrams later on) - x[good_id] += current_bidder_diagrams_[j].at(bidder_id).x_; - y[good_id] += current_bidder_diagrams_[j].at(bidder_id).y_; + if(useCustomWeights_) { + x[good_id] += weight * current_bidder_diagrams_[j].at(bidder_id).x_; + y[good_id] += weight * current_bidder_diagrams_[j].at(bidder_id).y_; + } else { + x[good_id] += current_bidder_diagrams_[j].at(bidder_id).x_; + y[good_id] += current_bidder_diagrams_[j].at(bidder_id).y_; + } if(geometrical_factor_ < 1) { const auto &critical_coordinates = current_bidder_diagrams_[j] .at(bidder_id) @@ -220,6 +233,10 @@ double ttk::PDBarycenter::updateBarycenter( // Counting the number of times this barycenter point is linked to the // diagonal count_diag_matchings[good_id] = count_diag_matchings[good_id] + 1; + if(useCustomWeights_) { + weight_diag_matchings[good_id] + = weight_diag_matchings[good_id] + weight; + } } } } @@ -229,20 +246,42 @@ double ttk::PDBarycenter::updateBarycenter( if(count_diag_matchings[i] < n_diagrams) { // Barycenter point i is matched at least to one off-diagonal bidder // 3.1 Compute the arithmetic mean of off-diagonal bidders linked to it - double const x_bar - = x[i] / (double)(n_diagrams - count_diag_matchings[i]); - double const y_bar - = y[i] / (double)(n_diagrams - count_diag_matchings[i]); + double x_bar = 0; + double y_bar = 0; + + if(useCustomWeights_) { + // weighted arithmetic mean of OFF-DIAGONAL points + // = sum_i weight[i]*x[i] / sum_j weight[j] + // if sum_j weight[j] = 0 (matched to off-diagonal points with weight + // zero), then skip the computation + if(weight_diag_matchings[i] < 1) { + x_bar = x[i] / (1 - weight_diag_matchings[i]); + y_bar = y[i] / (1 - weight_diag_matchings[i]); + } + } else { + x_bar = x[i] / (double)(n_diagrams - count_diag_matchings[i]); + y_bar = y[i] / (double)(n_diagrams - count_diag_matchings[i]); + } + // 3.2 Compute the new coordinates of the point (the more linked to the // diagonal it was, the closer to the diagonal it'll be) - double const new_x - = ((double)(n_diagrams - count_diag_matchings[i]) * x_bar - + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.) - / (double)n_diagrams; - double const new_y - = ((double)(n_diagrams - count_diag_matchings[i]) * y_bar - + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.) - / (double)n_diagrams; + double new_x = 0; + double new_y = 0; + + if(useCustomWeights_) { + new_x = (1 - weight_diag_matchings[i]) * x_bar + + weight_diag_matchings[i] * (x_bar + y_bar) / 2.; + new_y = (1 - weight_diag_matchings[i]) * y_bar + + weight_diag_matchings[i] * (x_bar + y_bar) / 2.; + } else { + new_x = ((double)(n_diagrams - count_diag_matchings[i]) * x_bar + + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.) + / (double)n_diagrams; + new_y = ((double)(n_diagrams - count_diag_matchings[i]) * y_bar + + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.) + / (double)n_diagrams; + } + // TODO Weight by persistence double const new_crit_coord_x = crit_coords_x[i] / (double)(n_diagrams - count_diag_matchings[i]); @@ -301,11 +340,17 @@ double ttk::PDBarycenter::updateBarycenter( // 5. Append the new points to the barycenter for(size_t k = 0; k < points_to_append.size(); k++) { points_added_ += 1; - Bidder *b = points_to_append[k]; - double const gx + Bidder *b = points_to_append[k].first; + double gx = (b->x_ + (n_diagrams - 1) * (b->x_ + b->y_) / 2.) / (n_diagrams); - double const gy + double gy = (b->y_ + (n_diagrams - 1) * (b->x_ + b->y_) / 2.) / (n_diagrams); + if(useCustomWeights_) { + // need to retrieve the correct weight here + double weight = points_to_append[k].second; + gx = weight * b->x_ + (1 - weight) * (b->x_ + b->y_) / 2.; + gy = weight * b->y_ + (1 - weight) * (b->x_ + b->y_) / 2.; + } const auto &critical_coordinates = b->GetCriticalCoordinates(); for(size_t j = 0; j < n_diagrams; j++) { Good g = Good(gx, gy, false, barycenter_goods_[j].size()); diff --git a/core/base/persistenceDiagramClustering/PDBarycenter.h b/core/base/persistenceDiagramClustering/PDBarycenter.h index d728f74e93..0ec82839df 100644 --- a/core/base/persistenceDiagramClustering/PDBarycenter.h +++ b/core/base/persistenceDiagramClustering/PDBarycenter.h @@ -187,6 +187,13 @@ namespace ttk { return cost_; } + inline void setUseCustomWeights(bool data) { + useCustomWeights_ = data; + } + inline void setCustomWeights(std::vector *pdata) { + customWeights_ = pdata; + } + protected: // std::vector precision_objective_; std::vector precision_; @@ -226,6 +233,9 @@ namespace ttk { std::vector> current_bidder_ids_; std::vector barycenter_goods_; + std::vector *customWeights_{}; + bool useCustomWeights_{false}; + bool reinit_prices_{true}; bool epsilon_decreases_{true}; bool early_stoppage_{true}; diff --git a/core/base/persistenceDiagramClustering/PDClustering.cpp b/core/base/persistenceDiagramClustering/PDClustering.cpp index ea96108fa8..84a017ab8f 100644 --- a/core/base/persistenceDiagramClustering/PDClustering.cpp +++ b/core/base/persistenceDiagramClustering/PDClustering.cpp @@ -2025,7 +2025,8 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( bool add_points_to_barycenter, bool first_enrichment) { - std::vector new_min_persistence = min_persistence; + std::vector new_min_persistence + = first_enrichment ? previous_min_persistence : min_persistence; if(!do_min_) { new_min_persistence[0] = previous_min_persistence[0]; @@ -2084,31 +2085,37 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( if(do_min_) { for(int i = 0; i < numberOfInputs_; i++) { std::vector persistences; + int number_of_added_pairs = current_bidder_diagrams_min_[i].size(); + double last_added_persistence = number_of_added_pairs + ? current_bidder_diagrams_min_[i] + .at(number_of_added_pairs - 1) + .getPersistence() + : previous_min_persistence[0]; for(size_t j = 0; j < bidder_diagrams_min_[i].size(); j++) { - Bidder const b = bidder_diagrams_min_[i].at(j); - double const persistence = b.getPersistence(); + const auto &b = bidder_diagrams_min_[i].at(j); + const auto persistence = b.getPersistence(); if(persistence >= min_persistence[0] - && persistence <= previous_min_persistence[0]) { + && persistence < last_added_persistence) { candidates_to_be_added_min[i].emplace_back(j); idx_min[i].emplace_back(idx_min[i].size()); persistences.emplace_back(persistence); } } - sort( + std::sort( idx_min[i].begin(), idx_min[i].end(), [&persistences](int &a, int &b) { return ((persistences[a] > persistences[b]) || ((persistences[a] == persistences[b]) && (a > b))); }); - int const size = candidates_to_be_added_min[i].size(); - if(size >= max_points_to_add_min) { - double const last_persistence_added_min - = persistences[idx_min[i][max_points_to_add_min - 1]]; + int size = candidates_to_be_added_min[i].size(); + if(size > 0) { + double last_persistence_added_min + = persistences[idx_min[i][std::min(size, max_points_to_add_min) - 1]]; if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add) // added per diagram if(i == 0) { new_min_persistence[0] = last_persistence_added_min; } else { - if(last_persistence_added_min < new_min_persistence[0]) + if(last_persistence_added_min > new_min_persistence[0]) new_min_persistence[0] = last_persistence_added_min; } } else { // a maxima max_point_to_add added per diagram @@ -2123,31 +2130,37 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( if(do_sad_) { for(int i = 0; i < numberOfInputs_; i++) { std::vector persistences; + int number_of_added_pairs = current_bidder_diagrams_saddle_[i].size(); + double last_added_persistence = number_of_added_pairs + ? current_bidder_diagrams_saddle_[i] + .at(number_of_added_pairs - 1) + .getPersistence() + : previous_min_persistence[1]; for(size_t j = 0; j < bidder_diagrams_saddle_[i].size(); j++) { - Bidder const b = bidder_diagrams_saddle_[i].at(j); - double const persistence = b.getPersistence(); + const auto &b = bidder_diagrams_saddle_[i].at(j); + const auto persistence = b.getPersistence(); if(persistence >= min_persistence[1] - && persistence <= previous_min_persistence[1]) { + && persistence < last_added_persistence) { candidates_to_be_added_sad[i].emplace_back(j); idx_sad[i].emplace_back(idx_sad[i].size()); persistences.emplace_back(persistence); } } - sort( + std::sort( idx_sad[i].begin(), idx_sad[i].end(), [&persistences](int &a, int &b) { return ((persistences[a] > persistences[b]) || ((persistences[a] == persistences[b]) && (a > b))); }); - int const size = candidates_to_be_added_sad[i].size(); - if(size >= max_points_to_add_sad) { - double const last_persistence_added_sad - = persistences[idx_sad[i][max_points_to_add_sad - 1]]; + int size = candidates_to_be_added_sad[i].size(); + if(size > 0) { + double last_persistence_added_sad + = persistences[idx_sad[i][std::min(size, max_points_to_add_sad) - 1]]; if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add) // added per diagram if(i == 0) { new_min_persistence[1] = last_persistence_added_sad; } else { - if(last_persistence_added_sad < new_min_persistence[1]) + if(last_persistence_added_sad > new_min_persistence[1]) new_min_persistence[1] = last_persistence_added_sad; } } else { // a maxima max_point_to_add added per diagram @@ -2161,31 +2174,38 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( if(do_max_) { for(int i = 0; i < numberOfInputs_; i++) { std::vector persistences; + int number_of_added_pairs = current_bidder_diagrams_max_[i].size(); + double last_added_persistence = number_of_added_pairs + ? current_bidder_diagrams_max_[i] + .at(number_of_added_pairs - 1) + .getPersistence() + : previous_min_persistence[2]; for(size_t j = 0; j < bidder_diagrams_max_[i].size(); j++) { - Bidder const b = bidder_diagrams_max_[i].at(j); - double const persistence = b.getPersistence(); + const auto &b = bidder_diagrams_max_[i].at(j); + const auto persistence = b.getPersistence(); + if(persistence >= min_persistence[2] - && persistence <= previous_min_persistence[2]) { + && persistence < last_added_persistence) { candidates_to_be_added_max[i].emplace_back(j); idx_max[i].emplace_back(idx_max[i].size()); persistences.emplace_back(persistence); } } - sort( + std::sort( idx_max[i].begin(), idx_max[i].end(), [&persistences](int &a, int &b) { return ((persistences[a] > persistences[b]) || ((persistences[a] == persistences[b]) && (a > b))); }); - int const size = candidates_to_be_added_max[i].size(); - if(size >= max_points_to_add_max) { - double const last_persistence_added_max - = persistences[idx_max[i][max_points_to_add_max - 1]]; + int size = candidates_to_be_added_max[i].size(); + if(size > 0) { + double last_persistence_added_max + = persistences[idx_max[i][std::min(size, max_points_to_add_max) - 1]]; if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add) // added per diagram if(i == 0) { new_min_persistence[2] = last_persistence_added_max; } else { - if(last_persistence_added_max < new_min_persistence[2]) + if(last_persistence_added_max > new_min_persistence[2]) new_min_persistence[2] = last_persistence_added_max; } } else { // a maxima max_point_to_add added per diagram @@ -2205,8 +2225,8 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( for(int j = 0; j < std::min(max_points_to_add_min, size); j++) { Bidder b = bidder_diagrams_min_[i].at( candidates_to_be_added_min[i][idx_min[i][j]]); - double const persistence = b.getPersistence(); - if(persistence >= new_min_persistence[0]) { + double persistence = b.getPersistence(); + if(persistence >= new_min_persistence[0] or first_enrichment) { b.id_ = current_bidder_diagrams_min_[i].size(); b.setPositionInAuction(current_bidder_diagrams_min_[i].size()); b.setDiagonalPrice(initial_diagonal_prices[0][i]); @@ -2265,8 +2285,8 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( for(int j = 0; j < std::min(max_points_to_add_sad, size); j++) { Bidder b = bidder_diagrams_saddle_[i].at( candidates_to_be_added_sad[i][idx_sad[i][j]]); - double const persistence = b.getPersistence(); - if(persistence >= new_min_persistence[1]) { + double persistence = b.getPersistence(); + if(persistence >= new_min_persistence[1] or first_enrichment) { b.id_ = current_bidder_diagrams_saddle_[i].size(); b.setPositionInAuction(current_bidder_diagrams_saddle_[i].size()); b.setDiagonalPrice(initial_diagonal_prices[1][i]); @@ -2321,8 +2341,8 @@ std::vector ttk::PDClustering::enrichCurrentBidderDiagrams( for(int j = 0; j < std::min(max_points_to_add_max, size); j++) { Bidder b = bidder_diagrams_max_[i].at( candidates_to_be_added_max[i][idx_max[i][j]]); - double const persistence = b.getPersistence(); - if(persistence >= new_min_persistence[2]) { + double persistence = b.getPersistence(); + if(persistence >= new_min_persistence[2] or first_enrichment) { b.id_ = current_bidder_diagrams_max_[i].size(); b.setPositionInAuction(current_bidder_diagrams_max_[i].size()); b.setDiagonalPrice(initial_diagonal_prices[2][i]); @@ -2396,6 +2416,16 @@ void ttk::PDClustering::initializeBarycenterComputers( barycenter_computer_min_[c].setDeterministic(true); barycenter_computer_min_[c].setGeometricalFactor(geometrical_factor_); barycenter_computer_min_[c].setDebugLevel(debugLevel_); + if(useCustomWeights_) { + if(customWeights_ != nullptr + and customWeights_->size() == diagrams_c.size()) { + barycenter_computer_min_[c].setUseCustomWeights(useCustomWeights_); + barycenter_computer_min_[c].setCustomWeights(customWeights_); + } else { + printErr("Weights incorrectly initialized, returning to balanced " + "barycenter"); + } + } barycenter_computer_min_[c].setNumberOfInputs(diagrams_c.size()); barycenter_computer_min_[c].setCurrentBidders(diagrams_c); barycenter_computer_min_[c].setNonMatchingWeight(nonMatchingWeight_); @@ -2418,7 +2448,16 @@ void ttk::PDClustering::initializeBarycenterComputers( barycenter_computer_sad_[c].setDebugLevel(debugLevel_); barycenter_computer_sad_[c].setNumberOfInputs(diagrams_c.size()); barycenter_computer_sad_[c].setCurrentBidders(diagrams_c); - barycenter_computer_sad_[c].setNonMatchingWeight(nonMatchingWeight_); + if(useCustomWeights_) { + if(customWeights_ != nullptr + and customWeights_->size() == diagrams_c.size()) { + barycenter_computer_sad_[c].setUseCustomWeights(useCustomWeights_); + barycenter_computer_sad_[c].setCustomWeights(customWeights_); + } else { + printErr("Weights incorrectly initialized, returning to balanced " + "barycenter"); + } + } std::vector barycenter_goods(clustering_[c].size()); for(size_t i_diagram = 0; i_diagram < clustering_[c].size(); @@ -2447,7 +2486,16 @@ void ttk::PDClustering::initializeBarycenterComputers( barycenter_computer_max_[c].setDebugLevel(debugLevel_); barycenter_computer_max_[c].setNumberOfInputs(diagrams_c.size()); barycenter_computer_max_[c].setCurrentBidders(diagrams_c); - barycenter_computer_max_[c].setNonMatchingWeight(nonMatchingWeight_); + if(useCustomWeights_) { + if(customWeights_ != nullptr + and customWeights_->size() == diagrams_c.size()) { + barycenter_computer_max_[c].setUseCustomWeights(useCustomWeights_); + barycenter_computer_max_[c].setCustomWeights(customWeights_); + } else { + printErr("Weights incorrectly initialized, returning to balanced " + "barycenter"); + } + } std::vector barycenter_goods(clustering_[c].size()); for(size_t i_diagram = 0; i_diagram < clustering_[c].size(); diff --git a/core/base/persistenceDiagramClustering/PDClustering.h b/core/base/persistenceDiagramClustering/PDClustering.h index 40be499601..2a3d224d61 100644 --- a/core/base/persistenceDiagramClustering/PDClustering.h +++ b/core/base/persistenceDiagramClustering/PDClustering.h @@ -242,6 +242,13 @@ namespace ttk { this->printMsg(msg.str()); } + inline void setUseCustomWeights(bool data) { + useCustomWeights_ = data; + } + inline void setCustomWeights(std::vector *pdata) { + customWeights_ = pdata; + } + protected: std::vector barycenter_computer_min_{}; std::vector barycenter_computer_sad_{}; @@ -322,5 +329,8 @@ namespace ttk { std::vector distanceToCentroid_{}; int n_iterations_; + + std::vector *customWeights_{}; + bool useCustomWeights_{false}; }; } // namespace ttk diff --git a/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.cpp b/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.cpp index 825eba5ae1..0fbb471694 100644 --- a/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.cpp +++ b/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.cpp @@ -1,5 +1,57 @@ #include +void ttk::computeWeightedBarycenter( + std::vector &intermediateDiagrams, + std::vector &weights, + DiagramType &barycenter, + std::vector> &matchings, + const ttk::Debug &dbg, + const bool ProgBarycenter) { + + std::vector final_centroids; + std::vector>> all_matchings; + + bool useCustomWeights = true; + + if(weights.size() == 0) { + dbg.printMsg("Uniform weights"); + useCustomWeights = false; + } else if(weights.size() != intermediateDiagrams.size()) { + dbg.printErr("Number of weights != Number of inputs"); + dbg.printErr("Defaulting to uniform barycenter"); + useCustomWeights = false; + } else { + std::stringstream msg; + double sum = 0.0; + // msg << "Weights: "; + for(auto w : weights) { + // msg << " " << w; + sum += w; + } + // msg << " | sum: " << sum; + // dbg.printMsg(msg.str()); + if(sum != 1) { + dbg.printWrn("Sum of weights different from 1"); + // dbg.printErr("Defaulting to uniform barycenter"); + // useCustomWeights = false; + } + } + + PersistenceDiagramClustering barycenterComputer{}; + barycenterComputer.setForceUseOfAlgorithm(true); + barycenterComputer.setUseCustomWeights(useCustomWeights); + barycenterComputer.setUseInterruptible(true); + barycenterComputer.setUseProgressive(ProgBarycenter); + barycenterComputer.setCustomWeights(&weights); + barycenterComputer.setDebugLevel(0); + barycenterComputer.setTimeLimit(10); + barycenterComputer.setThreadNumber(1); + barycenterComputer.execute( + intermediateDiagrams, final_centroids, all_matchings); + barycenter = std::move(final_centroids[0]); + matchings = std::move(all_matchings[0]); +} + std::vector ttk::PersistenceDiagramClustering::execute( std::vector &intermediateDiagrams, std::vector &final_centroids, @@ -121,6 +173,10 @@ std::vector ttk::PersistenceDiagramClustering::execute( KMeans.setDiagrams(&data_min, &data_sad, &data_max); KMeans.setDos(do_min, do_sad, do_max); KMeans.setNonMatchingWeight(NonMatchingWeight); + + KMeans.setUseCustomWeights(UseCustomWeights); + KMeans.setCustomWeights(CustomWeights); + inv_clustering = KMeans.execute(final_centroids, all_matchings_per_type_and_cluster); std::vector> centroids_sizes = KMeans.get_centroids_sizes(); diff --git a/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.h b/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.h index 7788a50301..b33aa229bc 100644 --- a/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.h +++ b/core/base/persistenceDiagramClustering/PersistenceDiagramClustering.h @@ -53,6 +53,25 @@ namespace ttk { return this->distances; } + inline void setTimeLimit(const double tl) { + TimeLimit = tl; + } + inline void setForceUseOfAlgorithm(bool data) { + ForceUseOfAlgorithm = data; + } + inline void setUseInterruptible(bool data) { + UseInterruptible = data; + } + inline void setUseProgressive(bool data) { + UseProgressive = data; + } + inline void setUseCustomWeights(bool data) { + UseCustomWeights = data; + } + inline void setCustomWeights(std::vector *pdata) { + CustomWeights = pdata; + } + protected: // Critical pairs used for clustering // 0:min-saddles ; 1:saddles-saddles ; 2:sad-max ; else : all @@ -82,6 +101,17 @@ namespace ttk { int points_added_; int points_deleted_; + + std::vector *CustomWeights{}; + bool UseCustomWeights{false}; }; + void + computeWeightedBarycenter(std::vector &intermediateDiagrams, + std::vector &weights, + DiagramType &barycenter, + std::vector> &matchings, + const ttk::Debug &dbg, + const bool ProgBarycenter); + } // namespace ttk diff --git a/core/base/persistenceDiagramDictionary/CMakeLists.txt b/core/base/persistenceDiagramDictionary/CMakeLists.txt new file mode 100644 index 0000000000..2603be829e --- /dev/null +++ b/core/base/persistenceDiagramDictionary/CMakeLists.txt @@ -0,0 +1,17 @@ +ttk_add_base_library(persistenceDiagramDictionary + SOURCES + PersistenceDiagramDictionary.cpp + HEADERS + PersistenceDiagramDictionary.h + DEPENDS + common + persistenceDiagramAuction + constrainedGradientDescent + persistenceDiagramClustering + initDictPersistenceDiagram +) + +if(TTK_ENABLE_EIGEN) + target_compile_definitions(persistenceDiagramDictionary PRIVATE TTK_ENABLE_EIGEN) + target_link_libraries(persistenceDiagramDictionary PRIVATE Eigen3::Eigen) +endif() \ No newline at end of file diff --git a/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.cpp b/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.cpp new file mode 100644 index 0000000000..904a91e7ed --- /dev/null +++ b/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.cpp @@ -0,0 +1,1724 @@ +#include "PersistenceDiagramUtils.h" +#include +#include +#ifdef TTK_ENABLE_EIGEN +#include +#include +#endif + +#include + +using namespace ttk; + +void PersistenceDiagramDictionary::execute( + std::vector &intermediateDiagrams, + const std::vector &intermediateAtoms, + std::vector &dictDiagrams, + std::vector> &vectorWeights, + const int seed, + const int numAtom, + std::vector &lossTab, + std::vector> &allLosses, + double percent) { + + if(!ProgApproach_) { + // Regular approach + for(size_t i = 0; i < intermediateDiagrams.size(); ++i) { + if(sortedForTest_) { + auto &diag = intermediateDiagrams[i]; + std::sort(diag.begin(), diag.end(), + [](const PersistencePair &t1, const PersistencePair &t2) { + return (t1.death.sfValue - t1.birth.sfValue) + > (t2.death.sfValue - t2.birth.sfValue); + }); + } + } + + bool doCompression = false; + if(CompressionMode_) { + doCompression = true; + } + + std::vector trueBidderDiagramMin( + intermediateDiagrams.size()); + std::vector trueBidderDiagramSad( + intermediateDiagrams.size()); + std::vector trueBidderDiagramMax( + intermediateDiagrams.size()); + + std::vector> histoVectorWeights( + intermediateDiagrams.size()); + std::vector histoDictDiagrams(numAtom); + this->maxLag2_ = 5; + + Timer tm_init{}; + bool preWeightOpt = false; + this->printMsg("Regular approach:"); + initDictionary(dictDiagrams, intermediateDiagrams, intermediateAtoms, + numAtom, this->do_min_, this->do_sad_, this->do_max_, seed, + percent); + this->printMsg("Initialization computed ", 1, tm_init.getElapsedTime(), + threadNumber_, debug::LineMode::NEW); + + controlAtomsSize(intermediateDiagrams, dictDiagrams); + + method(intermediateDiagrams, dictDiagrams, vectorWeights, numAtom, lossTab, + allLosses, histoVectorWeights, histoDictDiagrams, preWeightOpt, + percent, doCompression); + } else { + // Multi scale approach + bool doCompression = false; + for(size_t i = 0; i < intermediateDiagrams.size(); ++i) { + auto &diag = intermediateDiagrams[i]; + std::sort(diag.begin(), diag.end(), + [](const PersistencePair &t1, const PersistencePair &t2) { + return (t1.death.sfValue - t1.birth.sfValue) + > (t2.death.sfValue - t2.birth.sfValue); + }); + } + + int start = 20; + double stop = percent; + std::vector percentages; + for(int value = start; value > stop; value -= 5) { + percentages.push_back(static_cast(value) / 100.); + } + percentages.push_back(static_cast(percent) / 100.); + std::vector> histoVectorWeights( + intermediateDiagrams.size()); + std::vector histoDictDiagrams(numAtom); + std::vector dataTemp(intermediateDiagrams.size()); + bool preWeightOpt = true; + Timer tm_method{}; + for(size_t j = 0; j < 1; ++j) { + double percentage = percentages[j]; + this->maxLag2_ = 0; + this->printMsg("First step multi-scale approach:"); + for(size_t i = 0; i < intermediateDiagrams.size(); ++i) { + auto &diag = intermediateDiagrams[i]; + auto &t = diag[0]; + double maxPers = getMaxPers(diag); + dataTemp[i].push_back(t); + for(size_t p = 1; p < diag.size(); ++p) { + auto &t2 = diag[p]; + if(percentage * maxPers <= (t2.death.sfValue - t2.birth.sfValue)) { + dataTemp[i].push_back(t2); + } else { + continue; + } + } + } + + Timer tm_init{}; + initDictionary(dictDiagrams, dataTemp, intermediateAtoms, numAtom, + this->do_min_, this->do_sad_, this->do_max_, seed, + percent); + this->printMsg("Initialization computed ", 1, tm_init.getElapsedTime(), + threadNumber_, debug::LineMode::NEW); + + method(dataTemp, dictDiagrams, vectorWeights, numAtom, lossTab, allLosses, + histoVectorWeights, histoDictDiagrams, preWeightOpt, percent, + doCompression); + } + + int counter = 0; + for(size_t j = 1; j < percentages.size(); ++j) { + if(j == percentages.size() - static_cast(1) && CompressionMode_) { + doCompression = true; + } + double percentage = percentages[j]; + double previousPerc = percentages[j - 1]; + if(j < percentages.size() - 1) { + this->maxLag2_ = 0; + } else { + this->maxLag2_ = 5; + } + if(j < percentages.size() - 1) { + this->printMsg("New step multi-scale approach:"); + } else { + this->printMsg("Final step multi-scale approach:"); + } + + for(size_t i = 0; i < intermediateDiagrams.size(); ++i) { + auto &diag = intermediateDiagrams[i]; + double maxPers = getMaxPers(diag); + for(size_t p = 0; p < diag.size(); ++p) { + auto &t2 = diag[p]; + if(percentage * maxPers <= (t2.death.sfValue - t2.birth.sfValue) + && (t2.death.sfValue - t2.birth.sfValue) + < previousPerc * maxPers) { + dataTemp[i].push_back(t2); + counter += 1; + } + } + } + if(counter == 0) { + continue; + } + method(dataTemp, dictDiagrams, vectorWeights, numAtom, lossTab, allLosses, + histoVectorWeights, histoDictDiagrams, preWeightOpt, percent, + doCompression); + } + + this->printMsg( + "Total time", 1.0, tm_method.getElapsedTime(), threadNumber_); + } +} + +void PersistenceDiagramDictionary::method( + const std::vector &intermediateDiagrams, + std::vector &dictDiagrams, + std::vector> &vectorWeights, + const int numAtom, + std::vector &lossTab, + std::vector> &allLosses, + std::vector> &histoVectorWeights, + std::vector &histoDictDiagrams, + bool preWeightOpt, + double percent, + bool doCompression) { + + Timer tm{}; + + bool doOptimizeAtoms = false; + bool doOptimizeWeights = false; + + if(OptimizeWeights_) { + printMsg("Weight Optimization activated"); + doOptimizeWeights = true; + } else { + printWrn("Weight Optimization desactivated"); + } + if(OptimizeAtoms_) { + doOptimizeAtoms = true; + printMsg("Atom Optimization activated"); + } else { + printWrn("Atom Optimization desactivated"); + } + + const auto nDiags = intermediateDiagrams.size(); + + // tracking the original indices + std::vector inputDiagramsMin(nDiags); + std::vector inputDiagramsSad(nDiags); + std::vector inputDiagramsMax(nDiags); + + std::vector bidderDiagramsMin{}; + std::vector bidderDiagramsSad{}; + std::vector bidderDiagramsMax{}; + + std::vector> originIndexDatasMin(nDiags); + std::vector> originIndexDatasSad(nDiags); + std::vector> originIndexDatasMax(nDiags); + + gettingBidderDiagrams( + intermediateDiagrams, inputDiagramsMin, inputDiagramsSad, inputDiagramsMax, + bidderDiagramsMin, bidderDiagramsSad, bidderDiagramsMax, + originIndexDatasMin, originIndexDatasSad, originIndexDatasMax, true); + + std::vector barycentersList(nDiags); + std::vector>> allMatchingsAtoms( + nDiags); + std::vector> matchingsDatasMin(nDiags); + std::vector> matchingsDatasSad(nDiags); + std::vector> matchingsDatasMax(nDiags); + ConstrainedGradientDescent gradActor; + double loss; + int lag = 0; + int lag2 = 0; + int lag3 = 0; + int lagLimit = 10; + int minEpoch = 20; + int maxEpoch = MaxEpoch_; + bool cond = true; + int epoch = 0; + int nbEpochPrevious = lossTab.size(); + std::vector initSizes(dictDiagrams.size()); + for(size_t j = 0; j < dictDiagrams.size(); ++j) { + initSizes[j] = dictDiagrams[j].size(); + } + + double factEquiv = static_cast(numAtom); + double step = 1. / (2. * 2. * factEquiv); + gradActor.setStep(factEquiv); + + // BUFFERS + std::vector> bufferHistoAllEpochLife(dictDiagrams.size()); + std::vector> bufferHistoAllBoolLife(dictDiagrams.size()); + std::vector> bufferCheckUnderDiag(dictDiagrams.size()); + std::vector> bufferCheckDiag(dictDiagrams.size()); + std::vector> bufferCheckAboveGlobal(dictDiagrams.size()); + + std::vector> histoAllEpochLife(dictDiagrams.size()); + std::vector> histoAllBoolLife(dictDiagrams.size()); + std::vector> checkUnderDiag(dictDiagrams.size()); + std::vector> checkDiag(dictDiagrams.size()); + std::vector> checkAboveGlobal(dictDiagrams.size()); + + std::vector allLossesAtEpoch(nDiags, 0.); + std::vector trueAllLossesAtEpoch(nDiags, 0.); + + while(epoch < maxEpoch && cond) { + loss = 0.; + Timer tm_it{}; +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + //////////////////////////////WEIGHTS/////////////////////////////////// + for(size_t i = 0; i < nDiags; ++i) { + auto &barycenter = barycentersList[i]; + std::vector &weight = vectorWeights[i]; + std::vector> &matchings + = allMatchingsAtoms[i]; + computeWeightedBarycenter( + dictDiagrams, weight, barycenter, matchings, *this, ProgBarycenter_); + } + this->printMsg( + "Computed 1st Barycenters for epoch " + std::to_string(epoch), + epoch / static_cast(maxEpoch), tm_it.getElapsedTime(), + threadNumber_, debug::LineMode::NEW, debug::Priority::DETAIL); + + // "====================BARYCENTER FINISHED======================"); + + // tracking the original indices + std::vector barycentersListMin(nDiags); + std::vector barycentersListSad(nDiags); + std::vector barycentersListMax(nDiags); + + std::vector bidderBarycentersListMin{}; + std::vector bidderBarycentersListSad{}; + std::vector bidderBarycentersListMax{}; + + std::vector> originIndexBarysMin(nDiags); + std::vector> originIndexBarysSad(nDiags); + std::vector> originIndexBarysMax(nDiags); + + computeAllDistances( + barycentersList, nDiags, barycentersListMin, barycentersListSad, + barycentersListMax, bidderBarycentersListMin, bidderBarycentersListSad, + bidderBarycentersListMax, originIndexBarysMin, originIndexBarysSad, + originIndexBarysMax, bidderDiagramsMin, bidderDiagramsMax, + bidderDiagramsSad, matchingsDatasMin, matchingsDatasMax, + matchingsDatasSad, allLossesAtEpoch, true); + + for(size_t p = 0; p < nDiags; ++p) { + loss += allLossesAtEpoch[p]; + } + + for(size_t p = 0; p < nDiags; ++p) { + if(!ProgApproach_) { + allLosses[p].push_back(allLossesAtEpoch[p]); + } else { + allLosses[p].push_back(trueAllLossesAtEpoch[p]); + } + } + + lossTab.push_back(loss); + + printMsg( + " Epoch " + std::to_string(epoch) + ", loss = " + std::to_string(loss), 1, + threadNumber_, ttk::debug::LineMode::REPLACE); + + if(preWeightOpt && OptimizeAtoms_) { + if(epoch < 10) { + // Pre optimization of barycentric weights + doOptimizeAtoms = false; + } else { + doOptimizeAtoms = true; + } + } + + if(loss < 1e-7) { + cond = false; + doOptimizeWeights = false; + doOptimizeAtoms = false; + } + + double mini + = *std::min_element(lossTab.begin() + nbEpochPrevious, lossTab.end() - 1); + if(loss <= mini) { + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + const auto atom = dictDiagrams[p]; + histoDictDiagrams[p] = atom; + + const auto &histoEpochAtom = histoAllEpochLife[p]; + const auto &histoBoolAtom = histoAllBoolLife[p]; + const auto &boolUnderDiag = checkUnderDiag[p]; + const auto &boolDiag = checkDiag[p]; + const auto &boolAboveGlobal = checkAboveGlobal[p]; + + bufferHistoAllEpochLife[p] = histoEpochAtom; + bufferHistoAllBoolLife[p] = histoBoolAtom; + bufferCheckUnderDiag[p] = boolUnderDiag; + bufferCheckDiag[p] = boolDiag; + bufferCheckAboveGlobal[p] = boolAboveGlobal; + } + for(size_t p = 0; p < nDiags; ++p) { + const auto weights = vectorWeights[p]; + histoVectorWeights[p] = weights; + } + lag = 0; + lag3 = 0; + } else { + if(epoch > minEpoch) { + lag += 1; + } else { + lag = 0; + } + } + + if((epoch > minEpoch) + && (lossTab[epoch + nbEpochPrevious] + / lossTab[epoch + nbEpochPrevious - 1] + > 0.99)) { + if(lossTab[epoch + nbEpochPrevious] + < lossTab[epoch + nbEpochPrevious - 1]) { + if(lag2 == this->maxLag2_) { + this->printMsg("Loss not decreasing enough"); + if(StopCondition_) { + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + const auto atom = histoDictDiagrams[p]; + dictDiagrams[p] = atom; + } + for(size_t p = 0; p < nDiags; ++p) { + const auto weights = histoVectorWeights[p]; + vectorWeights[p] = weights; + } + doOptimizeWeights = false; + doOptimizeAtoms = false; + cond = false; + } + } else { + lag2 += 1; + } + } + } else { + lag2 = 0; + } + + if(epoch > minEpoch && lag > lagLimit) { + + if(StopCondition_) { + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + const auto atom = histoDictDiagrams[p]; + dictDiagrams[p] = atom; + } + for(size_t p = 0; p < nDiags; ++p) { + const auto weights = histoVectorWeights[p]; + vectorWeights[p] = weights; + } + this->printMsg("Minimum not passed"); + doOptimizeWeights = false; + doOptimizeAtoms = false; + + cond = false; + } + } + + if(cond && (epoch > 0) && (lossTab[epoch + nbEpochPrevious] > 2. * mini)) { + lag3 += 1; + lag2 = 0; + + if(lag3 > 2) { + this->printMsg("Loss increasing too much, reducing step and recompute " + "Barycenters and matchings"); + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + const auto atom = histoDictDiagrams[p]; + dictDiagrams[p] = atom; + + const auto &bufferHistoEpochAtom = bufferHistoAllEpochLife[p]; + const auto &bufferHistoBoolAtom = bufferHistoAllBoolLife[p]; + const auto &bufferBoolUnderDiag = bufferCheckUnderDiag[p]; + const auto &bufferBoolDiag = bufferCheckDiag[p]; + const auto &bufferBoolAboveGlobal = bufferCheckAboveGlobal[p]; + + histoAllEpochLife[p] = bufferHistoEpochAtom; + histoAllBoolLife[p] = bufferHistoBoolAtom; + checkUnderDiag[p] = bufferBoolUnderDiag; + checkDiag[p] = bufferBoolDiag; + checkAboveGlobal[p] = bufferBoolAboveGlobal; + } + + for(size_t p = 0; p < nDiags; ++p) { + const auto weights = histoVectorWeights[p]; + vectorWeights[p] = weights; + } + gradActor.reduceStep(); + step = step / 2.; + lag3 = 0; + + barycentersList.clear(); + barycentersList.resize(nDiags); + allMatchingsAtoms.clear(); + allMatchingsAtoms.resize(nDiags); + + barycentersListMin.clear(); + barycentersListSad.clear(); + barycentersListMax.clear(); + + barycentersListMin.resize(nDiags); + barycentersListSad.resize(nDiags); + barycentersListMax.resize(nDiags); + + bidderBarycentersListMin.clear(); + bidderBarycentersListSad.clear(); + bidderBarycentersListMax.clear(); + + originIndexBarysMin.clear(); + originIndexBarysSad.clear(); + originIndexBarysMax.clear(); + + originIndexBarysMin.resize(nDiags); + originIndexBarysSad.resize(nDiags); + originIndexBarysMax.resize(nDiags); + + matchingsDatasMin.clear(); + matchingsDatasSad.clear(); + matchingsDatasMax.clear(); + + matchingsDatasMin.resize(nDiags); + matchingsDatasSad.resize(nDiags); + matchingsDatasMax.resize(nDiags); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + //////////////////////////////WEIGHTS/////////////////////////////////// + for(size_t i = 0; i < nDiags; ++i) { + auto &barycenter = barycentersList[i]; + std::vector &weight = vectorWeights[i]; + std::vector> &matchings + = allMatchingsAtoms[i]; + computeWeightedBarycenter(dictDiagrams, weight, barycenter, matchings, + *this, ProgBarycenter_); + } + + computeAllDistances( + barycentersList, nDiags, barycentersListMin, barycentersListSad, + barycentersListMax, bidderBarycentersListMin, + bidderBarycentersListSad, bidderBarycentersListMax, + originIndexBarysMin, originIndexBarysSad, originIndexBarysMax, + bidderDiagramsMin, bidderDiagramsMax, bidderDiagramsSad, + matchingsDatasMin, matchingsDatasMax, matchingsDatasSad, + allLossesAtEpoch, false); + } + } + + std::vector> allHessianLists(nDiags); + std::vector> gradWeightsList(nDiags); + Timer tm_opt1{}; + + // WEIGHT OPTIMIZATION + if(doOptimizeWeights) { +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nDiags; ++i) { + auto &gradWeights = gradWeightsList[i]; + const auto &matchingsAtoms = allMatchingsAtoms[i]; + const auto &Barycenter = barycentersList[i]; + const auto &Data = intermediateDiagrams[i]; + std::vector &hessianList = allHessianLists[i]; + const std::vector &matchingsMin + = matchingsDatasMin[i]; + const std::vector &matchingsMax + = matchingsDatasMax[i]; + const std::vector &matchingsSad + = matchingsDatasSad[i]; + const std::vector &indexBaryMin = originIndexBarysMin[i]; + const std::vector &indexBarySad = originIndexBarysSad[i]; + const std::vector &indexBaryMax = originIndexBarysMax[i]; + const std::vector &indexDataMin = originIndexDatasMin[i]; + const std::vector &indexDataSad = originIndexDatasSad[i]; + const std::vector &indexDataMax = originIndexDatasMax[i]; + std::vector &weights = vectorWeights[i]; + computeGradientWeights(gradWeights, hessianList, dictDiagrams, + matchingsAtoms, Barycenter, Data, matchingsMin, + matchingsMax, matchingsSad, indexBaryMin, + indexBaryMax, indexBarySad, indexDataMin, + indexDataMax, indexDataSad, doOptimizeAtoms); + gradActor.executeWeightsProjected( + hessianList, weights, gradWeights, MaxEigenValue_); + } + + this->printMsg("Computed 1st opt for epoch " + std::to_string(epoch), + epoch / static_cast(maxEpoch), + tm_opt1.getElapsedTime(), threadNumber_, + debug::LineMode::NEW, debug::Priority::DETAIL); + } + + for(size_t p = 0; p < nDiags; ++p) { + allLossesAtEpoch[p] = 0.; + trueAllLossesAtEpoch[p] = 0.; + } + barycentersList.clear(); + barycentersList.resize(nDiags); + allMatchingsAtoms.clear(); + allMatchingsAtoms.resize(nDiags); + ////////////////////////////////ATOM//////////////////////////////////////// + if(doOptimizeAtoms) { + Timer tm_it2{}; +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nDiags; ++i) { + auto &barycenter = barycentersList[i]; + std::vector &weight = vectorWeights[i]; + std::vector> &matchings + = allMatchingsAtoms[i]; + computeWeightedBarycenter( + dictDiagrams, weight, barycenter, matchings, *this, ProgBarycenter_); + } + this->printMsg( + "Computed 2nd Barycenters for epoch " + std::to_string(epoch), + epoch / static_cast(maxEpoch), tm_it2.getElapsedTime(), + threadNumber_, debug::LineMode::NEW, debug::Priority::DETAIL); + + barycentersListMin.clear(); + barycentersListSad.clear(); + barycentersListMax.clear(); + + barycentersListMin.resize(nDiags); + barycentersListSad.resize(nDiags); + barycentersListMax.resize(nDiags); + + bidderBarycentersListMin.clear(); + bidderBarycentersListSad.clear(); + bidderBarycentersListMax.clear(); + + originIndexBarysMin.clear(); + originIndexBarysSad.clear(); + originIndexBarysMax.clear(); + + originIndexBarysMin.resize(nDiags); + originIndexBarysSad.resize(nDiags); + originIndexBarysMax.resize(nDiags); + + matchingsDatasMin.clear(); + matchingsDatasSad.clear(); + matchingsDatasMax.clear(); + + matchingsDatasMin.resize(nDiags); + matchingsDatasSad.resize(nDiags); + matchingsDatasMax.resize(nDiags); + + computeAllDistances( + barycentersList, nDiags, barycentersListMin, barycentersListSad, + barycentersListMax, bidderBarycentersListMin, bidderBarycentersListSad, + bidderBarycentersListMax, originIndexBarysMin, originIndexBarysSad, + originIndexBarysMax, bidderDiagramsMin, bidderDiagramsMax, + bidderDiagramsSad, matchingsDatasMin, matchingsDatasMax, + matchingsDatasSad, allLossesAtEpoch, false); + + std::vector>>> + allPairToAddToGradList(nDiags); + std::vector allInfoToAdd(nDiags); + std::vector> gradsAtomsList(nDiags); + std::vector> checkerAtomsList(nDiags); + + bool doDimReduct = false; + if(DimReductMode_ && numAtom <= 3) { + doDimReduct = true; + } else { + doDimReduct = false; + } + + Timer tm_opt2{}; +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nDiags; ++i) { + auto &pairToAddGradList = allPairToAddToGradList[i]; + auto &infoToAdd = allInfoToAdd[i]; + auto &gradsAtoms = gradsAtomsList[i]; + auto &checkerAtoms = checkerAtomsList[i]; + const auto &Barycenter = barycentersList[i]; + const auto &Data = intermediateDiagrams[i]; + const std::vector &matchingsMin + = matchingsDatasMin[i]; + const std::vector &matchingsMax + = matchingsDatasMax[i]; + const std::vector &matchingsSad + = matchingsDatasSad[i]; + const std::vector &indexBaryMin = originIndexBarysMin[i]; + const std::vector &indexBarySad = originIndexBarysSad[i]; + const std::vector &indexBaryMax = originIndexBarysMax[i]; + const std::vector &indexDataMin = originIndexDatasMin[i]; + const std::vector &indexDataSad = originIndexDatasSad[i]; + const std::vector &indexDataMax = originIndexDatasMax[i]; + const std::vector &weights = vectorWeights[i]; + computeGradientAtoms( + gradsAtoms, weights, Barycenter, Data, matchingsMin, matchingsMax, + matchingsSad, indexBaryMin, indexBaryMax, indexBarySad, indexDataMin, + indexDataMax, indexDataSad, checkerAtoms, pairToAddGradList, + infoToAdd, doDimReduct); + } + + std::vector maxiDeath(numAtom); + std::vector minBirth(numAtom); + for(int j = 0; j < numAtom; ++j) { + auto &atom = dictDiagrams[j]; + auto &temp = atom[0]; + maxiDeath[j] = temp.death.sfValue; + minBirth[j] = temp.birth.sfValue; + } + + std::vector>> allProjectionsList(nDiags); + std::vector allFeaturesToAdd(nDiags); + std::vector>> allProjLocations(nDiags); + std::vector>> + allVectorForProjContributions(nDiags); + std::vector allTrueFeaturesToAdd(numAtom); + std::vector>> allTrueProj(numAtom); + std::vector>> allTrueProjLoc(numAtom); + + for(size_t i = 0; i < nDiags; ++i) { + auto &pairToAddGradList = allPairToAddToGradList[i]; + auto &infoToAdd = allInfoToAdd[i]; + auto &projForDiag = allProjectionsList[i]; + auto &vectorForProjContrib = allVectorForProjContributions[i]; + auto &featuresToAdd = allFeaturesToAdd[i]; + auto &projLocations = allProjLocations[i]; + auto &gradsAtoms = gradsAtomsList[i]; + const auto &matchingsAtoms = allMatchingsAtoms[i]; + const auto &Barycenter = barycentersList[i]; + const auto &checkerAtoms = checkerAtomsList[i]; + gradActor.executeAtoms( + dictDiagrams, matchingsAtoms, Barycenter, gradsAtoms, checkerAtoms, + projForDiag, featuresToAdd, projLocations, vectorForProjContrib, + pairToAddGradList, infoToAdd); + } + + if(CreationFeatures_) { + for(size_t i = 0; i < nDiags; ++i) { + auto &projForDiag = allProjectionsList[i]; + auto &featuresToAdd = allFeaturesToAdd[i]; + auto &projLocations = allProjLocations[i]; + auto &vectorForProjContrib = allVectorForProjContributions[i]; + for(size_t j = 0; j < projForDiag.size(); ++j) { + auto &t = featuresToAdd[j]; + std::array &pair = projLocations[j]; + std::vector &vectorContrib = vectorForProjContrib[j]; + std::vector &projAndIndex = projForDiag[j]; + std::vector proj(numAtom); + for(int m = 0; m < numAtom; ++m) { + proj[m] = projAndIndex[m]; + } + int atomIndex = static_cast(projAndIndex[numAtom]); + bool lenNull = allTrueProj[atomIndex].size() == 0; + if(lenNull) { + const CriticalType c1 = t.birth.type; + const CriticalType c2 = t.death.type; + const SimplexId idTemp = t.dim; + pair[0] = pair[0] - step * vectorContrib[0]; + pair[1] = pair[1] - step * vectorContrib[1]; + if(pair[0] > pair[1]) { + pair[1] = pair[0]; + } + if(pair[0] < minBirth[atomIndex]) { + continue; + } + if(pair[1] - pair[0] < 1e-7) { + continue; + } + + PersistencePair newPair{CriticalVertex{0, c1, pair[0], {}}, + CriticalVertex{0, c2, pair[1], {}}, + idTemp, true}; + allTrueProj[atomIndex].push_back(proj); + allTrueFeaturesToAdd[atomIndex].push_back(newPair); + allTrueProjLoc[atomIndex].push_back(pair); + } else { + bool ralph = true; + size_t index = 0; + if(Fusion_) { + for(size_t n = 0; n < allTrueProj[atomIndex].size(); ++n) { + auto &projStocked = allTrueProj[atomIndex][n]; + auto &projLocStocked = allTrueProj[atomIndex][n]; + double distance + = sqrt(pow((pair[0] - projLocStocked[0]), 2) + + pow((pair[1] - projLocStocked[1]), 2)); + if(proj == projStocked && distance < 1e-3) { + ralph = false; + index = n; + break; + } + } + } + if(ralph) { + + const CriticalType c1 = t.birth.type; + const CriticalType c2 = t.death.type; + const SimplexId idTemp = t.dim; + pair[0] = pair[0] - step * vectorContrib[0]; + pair[1] = pair[1] - step * vectorContrib[1]; + if(pair[0] > pair[1]) { + pair[1] = pair[0]; + } + if(pair[0] < minBirth[atomIndex]) { + continue; + } + if(pair[1] - pair[0] < 1e-7) { + continue; + } + + PersistencePair newPair{CriticalVertex{0, c1, pair[0], {}}, + CriticalVertex{0, c2, pair[1], {}}, + idTemp, true}; + allTrueProj[atomIndex].push_back(proj); + allTrueFeaturesToAdd[atomIndex].push_back(newPair); + allTrueProjLoc[atomIndex].push_back(pair); + + } else { + auto &tReal = allTrueFeaturesToAdd[atomIndex][index]; + tReal.birth.sfValue + = tReal.birth.sfValue - step * vectorContrib[0]; + tReal.death.sfValue + = tReal.death.sfValue - step * vectorContrib[1]; + if(tReal.birth.sfValue > tReal.death.sfValue) { + tReal.death.sfValue = tReal.birth.sfValue; + } + } + } + } + } + + if(!doCompression) { + for(int i = 0; i < numAtom; ++i) { + auto &atom = dictDiagrams[i]; + auto &histoEpochAtom = histoAllEpochLife[i]; + auto &histoBoolAtom = histoAllBoolLife[i]; + auto &boolUnderDiag = checkUnderDiag[i]; + auto &boolDiag = checkDiag[i]; + auto initSize = initSizes[i]; + auto &boolAboveGlobal = checkAboveGlobal[i]; + if(histoEpochAtom.size() > 0) { + for(size_t j = 0; j < histoEpochAtom.size(); ++j) { + auto &t = atom[initSize + j]; + histoEpochAtom[j] += 1; + histoBoolAtom[j] = t.death.sfValue - t.birth.sfValue + < 0.1 * (percent / 100.) * maxiDeath[i]; + boolDiag[j] = t.death.sfValue - t.birth.sfValue < 1e-6; + boolUnderDiag[j] = t.death.sfValue < t.birth.sfValue; + boolAboveGlobal[j] = t.birth.sfValue > maxiDeath[i]; + } + } + } + + for(int i = 0; i < numAtom; ++i) { + auto &atom = dictDiagrams[i]; + auto &histoEpochAtom = histoAllEpochLife[i]; + auto &histoBoolAtom = histoAllBoolLife[i]; + auto &boolUnderDiag = checkUnderDiag[i]; + auto &boolDiag = checkDiag[i]; + auto &trueFeaturesToAdd = allTrueFeaturesToAdd[i]; + auto &boolAboveGlobal = checkAboveGlobal[i]; + for(size_t j = 0; j < trueFeaturesToAdd.size(); ++j) { + auto &t = trueFeaturesToAdd[j]; + atom.push_back(t); + histoEpochAtom.push_back(0); + histoBoolAtom.push_back(t.death.sfValue - t.birth.sfValue + < 0.1 * (percent / 100.) * maxiDeath[i]); + boolDiag.push_back(t.death.sfValue - t.birth.sfValue < 1e-6); + boolUnderDiag.push_back(t.death.sfValue < t.birth.sfValue); + boolAboveGlobal.push_back(t.birth.sfValue > maxiDeath[i]); + } + } + + std::vector> allIndicesToDelete(numAtom); + for(int i = 0; i < numAtom; ++i) { + auto &indicesAtomToDelete = allIndicesToDelete[i]; + auto &histoEpochAtom = histoAllEpochLife[i]; + auto &histoBoolAtom = histoAllBoolLife[i]; + auto &boolUnderDiag = checkUnderDiag[i]; + auto &boolDiag = checkDiag[i]; + auto &boolAboveGlobal = checkAboveGlobal[i]; + for(size_t j = 0; j < histoEpochAtom.size(); ++j) { + if(boolUnderDiag[j] || boolDiag[j] || boolAboveGlobal[j] + || (histoEpochAtom[j] > 5 && histoBoolAtom[j])) { + indicesAtomToDelete.push_back(j); + } + } + } + + for(int i = 0; i < numAtom; ++i) { + auto &atom = dictDiagrams[i]; + auto &histoEpochAtom = histoAllEpochLife[i]; + auto &histoBoolAtom = histoAllBoolLife[i]; + auto &indicesAtomToDelete = allIndicesToDelete[i]; + auto &boolUnderDiag = checkUnderDiag[i]; + auto &boolDiag = checkDiag[i]; + auto &boolAboveGlobal = checkAboveGlobal[i]; + auto initSize = initSizes[i]; + // size_t initAtomSize = atom.size() -histoEpochAtom.size(); + if(static_cast(indicesAtomToDelete.size()) > 0) { + for(int j = static_cast(indicesAtomToDelete.size()) - 1; + j >= 0; j--) { + atom.erase(atom.begin() + initSize + indicesAtomToDelete[j]); + histoEpochAtom.erase(histoEpochAtom.begin() + + indicesAtomToDelete[j]); + histoBoolAtom.erase(histoBoolAtom.begin() + + indicesAtomToDelete[j]); + boolUnderDiag.erase(boolUnderDiag.begin() + + indicesAtomToDelete[j]); + boolDiag.erase(boolDiag.begin() + indicesAtomToDelete[j]); + boolAboveGlobal.erase(boolAboveGlobal.begin() + + indicesAtomToDelete[j]); + } + } else { + continue; + } + } + } else { + for(int i = 0; i < numAtom; ++i) { + auto &atom = dictDiagrams[i]; + auto &trueFeaturesToAdd = allTrueFeaturesToAdd[i]; + for(size_t j = 0; j < trueFeaturesToAdd.size(); ++j) { + auto &t = trueFeaturesToAdd[j]; + atom.push_back(t); + } + } + controlAtomsSize(intermediateDiagrams, dictDiagrams); + } + } else { + if(doCompression) { + if(epoch > -1) { + controlAtomsSize(intermediateDiagrams, dictDiagrams); + } + } + } + // Deleting unallowed pairs: + for(int i = 0; i < numAtom; ++i) { + auto &atom = dictDiagrams[i]; + auto &globalPair = atom[0]; + for(size_t j = 0; j < atom.size(); ++j) { + auto &t = atom[j]; + + if(t.death.sfValue > globalPair.death.sfValue) { + t.death.sfValue = globalPair.death.sfValue; + } + + if(t.birth.sfValue > t.death.sfValue) { + t.death.sfValue = t.birth.sfValue; + } + + if(t.birth.sfValue < globalPair.birth.sfValue) { + t.birth.sfValue = globalPair.birth.sfValue; + } + } + } + this->printMsg("Computed 2nd opt for epoch " + std::to_string(epoch), + epoch / static_cast(maxEpoch), + tm_opt2.getElapsedTime(), threadNumber_, + debug::LineMode::NEW, debug::Priority::DETAIL); + + // ATOM OPTIMIZATION + } + epoch += 1; + barycentersList.clear(); + barycentersList.resize(nDiags); + allMatchingsAtoms.clear(); + allMatchingsAtoms.resize(nDiags); + } + printMsg( + " Epoch " + std::to_string(epoch) + ", loss = " + std::to_string(loss), 1, + threadNumber_); + + printMsg("Loss returned " + + std::to_string(*std::min_element( + lossTab.begin() + nbEpochPrevious, lossTab.end())) + + " at Epoch " + + std::to_string( + std::min_element(lossTab.begin() + nbEpochPrevious, lossTab.end()) + - lossTab.begin())); + + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + auto atom = histoDictDiagrams[p]; + dictDiagrams[p] = atom; + } + for(size_t p = 0; p < nDiags; ++p) { + auto weights = histoVectorWeights[p]; + vectorWeights[p] = weights; + } + + this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_); +} + +double PersistenceDiagramDictionary::distVect( + const std::vector &vec1, const std::vector &vec2) const { + + double dist = 0.; + for(size_t i = 0; i < vec1.size(); ++i) { + dist = dist + (vec1[i] - vec2[i]) * (vec1[i] - vec2[i]); + } + return std::sqrt(dist); +} + +double PersistenceDiagramDictionary::getMostPersistent( + const std::vector &bidder_diags) const { + + double max_persistence = 0; + + for(unsigned int i = 0; i < bidder_diags.size(); ++i) { + for(size_t j = 0; j < bidder_diags[i].size(); ++j) { + const double persistence = bidder_diags[i].at(j).getPersistence(); + if(persistence > max_persistence) { + max_persistence = persistence; + } + } + } + + return max_persistence; +} + +double PersistenceDiagramDictionary::computeDistance( + const BidderDiagram &D1, + const BidderDiagram &D2, + std::vector &matching) const { + + GoodDiagram D2_bis{}; + for(size_t i = 0; i < D2.size(); i++) { + const Bidder &b = D2.at(i); + Good g(b.x_, b.y_, b.isDiagonal(), D2_bis.size()); + g.SetCriticalCoordinates(b.coords_); + g.setPrice(0); + D2_bis.emplace_back(g); + } + + PersistenceDiagramAuction auction( + this->Wasserstein, this->Alpha, this->Lambda, this->DeltaLim, true); + auction.BuildAuctionDiagrams(D1, D2_bis); + double loss = auction.run(matching); + return loss; +} + +void PersistenceDiagramDictionary::computeGradientWeights( + std::vector &gradWeights, + std::vector &hessianList, + const std::vector &dictDiagrams, + const std::vector> &matchingsAtoms, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsMin, + const std::vector &matchingsMax, + const std::vector &matchingsSad, + const std::vector &indexBaryMin, + const std::vector &indexBaryMax, + const std::vector &indexBarySad, + const std::vector &indexDataMin, + const std::vector &indexDataMax, + const std::vector &indexDataSad, + const bool doOptimizeAtoms) const { + + // initialization + std::vector>> gradBuffersList( + Barycenter.size()); + std::vector>> pairToAddGradList; + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + gradBuffersList[i].resize(matchingsAtoms.size()); + } + std::vector> directions(Barycenter.size()); + std::vector> dataAssigned(Barycenter.size()); + gradWeights.resize(dictDiagrams.size()); + for(size_t i = 0; i < dictDiagrams.size(); ++i) { + gradWeights[i] = 0.; + } + + std::vector> checker(Barycenter.size()); + for(size_t j = 0; j < Barycenter.size(); ++j) { + checker[j].resize(matchingsAtoms.size()); + } + std::vector tracker(Barycenter.size(), 0); + std::vector tracker2(Barycenter.size(), 0); + + // computing gradients + for(size_t i = 0; i < matchingsAtoms.size(); ++i) { + for(size_t j = 0; j < matchingsAtoms[i].size(); ++j) { + const ttk::MatchingType &t = matchingsAtoms[i][j]; + // Id in atom + const SimplexId Id1 = std::get<0>(t); + // Id in barycenter + const SimplexId Id2 = std::get<1>(t); + // if(Id2 < 0) { + if(Id2 < 0 || static_cast(gradBuffersList.size()) <= Id2 + || static_cast(dictDiagrams[i].size()) <= Id1) { + continue; + } else if(Id1 < 0) { + const PersistencePair &t3 = Barycenter[Id2]; + auto &point = gradBuffersList[Id2][i]; + const double birthBarycenter = t3.birth.sfValue; + const double deathBarycenter = t3.death.sfValue; + const double birth_death_atom + = birthBarycenter + (deathBarycenter - birthBarycenter) / 2.; + point[0] = birth_death_atom; + point[1] = birth_death_atom; + checker[Id2][i] = i; + tracker[Id2] = 1; + } else { + const PersistencePair &t2 = dictDiagrams[i][Id1]; + auto &point = gradBuffersList[Id2][i]; + const double birth_atom = t2.birth.sfValue; + const double death_atom = t2.death.sfValue; + point[0] = birth_atom; + point[1] = death_atom; + checker[Id2][i] = i; + tracker[Id2] = 1; + } + } + } + + // Compute directions for min diagram + computeDirectionsGradWeight(matchingsAtoms, Barycenter, newData, matchingsMin, + indexBaryMin, indexDataMin, pairToAddGradList, + directions, dataAssigned, tracker2, + doOptimizeAtoms); + + // Compute directions for max diagram + computeDirectionsGradWeight(matchingsAtoms, Barycenter, newData, matchingsMax, + indexBaryMax, indexDataMax, pairToAddGradList, + directions, dataAssigned, tracker2, + doOptimizeAtoms); + + // Compute directions for sad diagram + computeDirectionsGradWeight(matchingsAtoms, Barycenter, newData, matchingsSad, + indexBarySad, indexDataSad, pairToAddGradList, + directions, dataAssigned, tracker2, + doOptimizeAtoms); + + std::vector temp(pairToAddGradList.size(), 1); + gradBuffersList.insert( + gradBuffersList.end(), pairToAddGradList.begin(), pairToAddGradList.end()); + tracker2.insert(tracker2.end(), temp.begin(), temp.end()); + tracker.insert(tracker.end(), temp.begin(), temp.end()); + std::vector temp2(matchingsAtoms.size()); + for(size_t j = 0; j < matchingsAtoms.size(); ++j) { + temp2.push_back(static_cast(j)); + } + for(size_t j = 0; j < pairToAddGradList.size(); ++j) { + checker.push_back(temp2); + } + + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + const auto &data_point = dataAssigned[i]; + for(size_t j = 0; j < checker[i].size(); ++j) { + auto &point = gradBuffersList[i][checker[i][j]]; + point[0] -= data_point[0]; + point[1] -= data_point[1]; + } + } + + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + if(tracker[i] == 0 || tracker2[i] == 0) { + continue; + } else { + for(size_t j = 0; j < checker[i].size(); ++j) { + const auto &point = gradBuffersList[i][checker[i][j]]; + const auto &direction = directions[i]; + gradWeights[checker[i][j]] + += -2 * (point[0] * direction[0] + point[1] * direction[1]); + } + } + } + hessianList.resize(gradBuffersList.size()); + for(size_t i = 0; i < gradBuffersList.size(); ++i) { + Matrix &hessian = hessianList[i]; + hessian.resize(checker[i].size()); + for(size_t j = 0; j < checker[i].size(); ++j) { + auto &line = hessian[j]; + line.resize(checker[i].size()); + const auto &point = gradBuffersList[i][checker[i][j]]; + for(size_t q = 0; q < checker[i].size(); ++q) { + const auto &point_temp = gradBuffersList[i][checker[i][q]]; + line[q] = point[0] * point_temp[0] + point[1] * point_temp[1]; + } + } + } +} + +void PersistenceDiagramDictionary::computeGradientAtoms( + std::vector &gradsAtoms, + const std::vector &weights, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsMin, + const std::vector &matchingsMax, + const std::vector &matchingsSad, + const std::vector &indexBaryMin, + const std::vector &indexBaryMax, + const std::vector &indexBarySad, + const std::vector &indexDataMin, + const std::vector &indexDataMax, + const std::vector &indexDataSad, + std::vector &checker, + std::vector>> &pairToAddGradList, + std::vector &infoToAdd, + bool doDimReduct) const { + + gradsAtoms.resize(Barycenter.size()); + for(size_t i = 0; i < Barycenter.size(); ++i) { + gradsAtoms[i].resize(weights.size()); + } + + checker.resize(Barycenter.size()); + for(size_t i = 0; i < Barycenter.size(); ++i) { + checker[i] = 0; + } + std::vector> directions(Barycenter.size()); + + // Compute directions for min diagram + computeDirectionsGradAtoms(gradsAtoms, Barycenter, weights, newData, + matchingsMin, indexBaryMin, indexDataMin, + pairToAddGradList, directions, checker, infoToAdd, + doDimReduct); + + // Compute directions for max diagram + computeDirectionsGradAtoms(gradsAtoms, Barycenter, weights, newData, + matchingsMax, indexBaryMax, indexDataMax, + pairToAddGradList, directions, checker, infoToAdd, + doDimReduct); + + // Compute directions for sad diagram + computeDirectionsGradAtoms(gradsAtoms, Barycenter, weights, newData, + matchingsSad, indexBarySad, indexDataSad, + pairToAddGradList, directions, checker, infoToAdd, + doDimReduct); + + for(size_t i = 0; i < Barycenter.size(); ++i) { + if(checker[i] == 0) { + continue; + } else { + for(size_t j = 0; j < weights.size(); ++j) { + std::vector temp(2); + const std::vector &direction = directions[i]; + temp[0] = -2 * weights[j] * direction[0]; + temp[1] = -2 * weights[j] * direction[1]; + gradsAtoms[i][j] = temp; + } + } + } +} + +void PersistenceDiagramDictionary::setBidderDiagrams( + const size_t nInputs, + std::vector &inputDiagrams, + std::vector &bidder_diags) const { + + bidder_diags.resize(nInputs); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nInputs; i++) { + auto &diag = inputDiagrams[i]; + auto &bidders = bidder_diags[i]; + + for(size_t j = 0; j < diag.size(); j++) { + // Add bidder to bidders + Bidder b(diag[j], j, this->Lambda); + b.setPositionInAuction(bidders.size()); + bidders.emplace_back(b); + if(b.isDiagonal() || b.x_ == b.y_) { + this->printMsg("Diagonal point in diagram " + std::to_string(i) + "!", + ttk::debug::Priority::DETAIL); + } + } + } +} + +int PersistenceDiagramDictionary::initDictionary( + std::vector &dictDiagrams, + const std::vector &datas, + const std::vector &inputAtoms, + const int nbAtom, + bool do_min, + bool do_sad, + bool do_max, + int seed, + double percent) { + switch(this->BackEnd) { + case BACKEND::INPUT_ATOMS: { + if(static_cast(inputAtoms.size()) != nbAtom) { + printErr("number of atoms don't match, going with " + + std::to_string(inputAtoms.size()) + " atoms"); + } + for(size_t i = 0; i < inputAtoms.size(); i++) { + const auto &t = inputAtoms[i]; + double maxPers = getMaxPers(t); + + dictDiagrams.push_back(t); + dictDiagrams[i].erase( + std::remove_if(dictDiagrams[i].begin(), dictDiagrams[i].end(), + [maxPers](ttk::PersistencePair &p) { + return (p.death.sfValue - p.birth.sfValue) + < 0. * maxPers; + }), + dictDiagrams[i].end()); + } + break; + } + + case BACKEND::BORDER_INIT: { + ttk::InitDictPersistenceDiagram initializer; + initializer.setThreadNumber(this->threadNumber_); + initializer.execute(dictDiagrams, datas, nbAtom, do_min, do_sad, do_max); + break; + } + + case BACKEND::RANDOM_INIT: + ttk::InitRandomDict{}.execute(dictDiagrams, datas, nbAtom, seed); // TODO + break; + + case BACKEND::FIRST_DIAGS: { + for(int i = 0; i < nbAtom; ++i) { + const auto &t = datas[i]; + dictDiagrams.push_back(t); + } + break; + } + case BACKEND::GREEDY_INIT: { + dictDiagrams.resize(datas.size()); + for(size_t i = 0; i < datas.size(); ++i) { + dictDiagrams[i] = datas[i]; + } + while(static_cast(dictDiagrams.size()) != nbAtom) { + std::vector allEnergy(dictDiagrams.size()); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t j = 0; j < dictDiagrams.size(); ++j) { + std::vector dictTemp; + std::vector dataAlone; + std::vector lossTabTemp; + std::vector> allLossesTemp(1); + + for(size_t p = 0; p < dictDiagrams.size(); ++p) { + if(p != j) { + dictTemp.push_back(dictDiagrams[p]); + } else { + dataAlone.push_back(dictDiagrams[p]); + } + } + std::vector> weightsTemp(1); + std::vector weights( + dictTemp.size(), 1. / (static_cast(dictTemp.size()) * 1.)); + weightsTemp[0] = weights; + std::vector> histoVectorWeights(1); + std::vector histoDictDiagrams(dictTemp.size()); + bool doCompression = false; + this->method(dataAlone, dictTemp, weightsTemp, + static_cast(dictTemp.size()), lossTabTemp, + allLossesTemp, histoVectorWeights, histoDictDiagrams, + false, percent, doCompression); + double min_loss + = *std::min_element(lossTabTemp.begin(), lossTabTemp.end()); + allEnergy[j] = min_loss; + } + int index = std::min_element(allEnergy.begin(), allEnergy.end()) + - allEnergy.begin(); + dictDiagrams.erase(dictDiagrams.begin() + index); + } + break; + } + default: + break; + } + return 0; +} + +void PersistenceDiagramDictionary::gettingBidderDiagrams( + const std::vector &intermediateDiagrams, + std::vector &inputDiagramsMin, + std::vector &inputDiagramsSad, + std::vector &inputDiagramsMax, + std::vector &bidderDiagramsMin, + std::vector &bidderDiagramsSad, + std::vector &bidderDiagramsMax, + std::vector> &originIndexMin, + std::vector> &originIndexSad, + std::vector> &originIndexMax, + bool insertOriginIndexMode) const { + + size_t nDiags = intermediateDiagrams.size(); + + // Create diagrams for min, saddle and max persistence pairs +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nDiags; i++) { + const auto &CTDiagram = intermediateDiagrams[i]; + + for(size_t j = 0; j < CTDiagram.size(); ++j) { + const PersistencePair &t = CTDiagram[j]; + const ttk::CriticalType nt1 = t.birth.type; + const ttk::CriticalType nt2 = t.death.type; + const double pers = t.persistence(); + if(pers > 0) { + if(nt1 == CriticalType::Local_minimum + && nt2 == CriticalType::Local_maximum) { + inputDiagramsMin[i].emplace_back(t); + if(insertOriginIndexMode) { + originIndexMin[i].push_back(j); + } + } else { + if(nt1 == CriticalType::Local_maximum + || nt2 == CriticalType::Local_maximum) { + inputDiagramsMax[i].emplace_back(t); + if(insertOriginIndexMode) { + originIndexMax[i].push_back(j); + } + } + if(nt1 == CriticalType::Local_minimum + || nt2 == CriticalType::Local_minimum) { + inputDiagramsMin[i].emplace_back(t); + if(insertOriginIndexMode) { + originIndexMin[i].push_back(j); + } + } + if((nt1 == CriticalType::Saddle1 && nt2 == CriticalType::Saddle2) + || (nt1 == CriticalType::Saddle2 + && nt2 == CriticalType::Saddle1)) { + inputDiagramsSad[i].emplace_back(t); + if(insertOriginIndexMode) { + originIndexSad[i].push_back(j); + } + } + } + } + } + } + + if(this->do_min_) { + setBidderDiagrams(nDiags, inputDiagramsMin, bidderDiagramsMin); + } + if(this->do_sad_) { + setBidderDiagrams(nDiags, inputDiagramsSad, bidderDiagramsSad); + } + if(this->do_max_) { + setBidderDiagrams(nDiags, inputDiagramsMax, bidderDiagramsMax); + } +} + +double PersistenceDiagramDictionary::getMaxPers(const ttk::DiagramType &data) { + double maxPers = 0.; + for(size_t j = 0; j < data.size(); ++j) { + auto &t = data[j]; + double pers = t.death.sfValue - t.birth.sfValue; + if(pers > maxPers) { + maxPers = pers; + } + } + + return maxPers; +} + +void PersistenceDiagramDictionary::controlAtomsSize( + const std::vector &intermediateDiagrams, + std::vector &dictDiagrams) const { + + size_t m = dictDiagrams.size(); + int globalSize = 0; + + for(size_t j = 0; j < intermediateDiagrams.size(); ++j) { + auto &data = intermediateDiagrams[j]; + globalSize += static_cast(data.size()); + } + + int dictSize = 0; + + for(size_t j = 0; j < m; ++j) { + auto &atom = dictDiagrams[j]; + dictSize += static_cast(atom.size()); + } + + if(static_cast(dictSize) + > (1. / this->CompressionFactor) * static_cast(globalSize)) { + double factor + = (1. / this->CompressionFactor) + * (static_cast(globalSize) / static_cast(dictSize)); + std::vector> tempDictPersistencePairs(m); + + for(size_t j = 0; j < m; ++j) { + auto &temp = tempDictPersistencePairs[j]; + auto &atom = dictDiagrams[j]; + temp.resize(atom.size()); + for(size_t p = 0; p < atom.size(); ++p) { + auto &t = atom[p]; + temp[p] = t.persistence(); + } + } + + for(size_t j = 0; j < m; ++j) { + auto &temp = tempDictPersistencePairs[j]; + std::sort(temp.begin(), temp.end(), std::greater()); + } + + std::vector persThresholds(m); + for(size_t j = 0; j < m; ++j) { + auto &temp = tempDictPersistencePairs[j]; + int index + = static_cast(floor(factor * static_cast(temp.size()))); + persThresholds[j] = temp[index]; + } + + for(size_t j = 0; j < m; ++j) { + double persThreshold = persThresholds[j]; + auto &atom = dictDiagrams[j]; + atom.erase(std::remove_if(atom.begin(), atom.end(), + [persThreshold](ttk::PersistencePair &t) { + return (t.death.sfValue - t.birth.sfValue) + < persThreshold; + }), + atom.end()); + } + } +} + +void PersistenceDiagramDictionary::computeDirectionsGradWeight( + const std::vector> &matchingsAtoms, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsCritType, + const std::vector &indexBaryCritType, + const std::vector &indexDataCritType, + std::vector>> &pairToAddGradList, + std::vector> &directions, + std::vector> &dataAssigned, + std::vector &tracker2, + const bool doOptimizeAtoms) const { + + size_t m = matchingsCritType.size(); + // int k = 0; + for(size_t i = 0; i < m; ++i) { + const ttk::MatchingType &t = matchingsCritType[i]; + // Id in newData + const SimplexId Id1 = std::get<0>(t); + // Id in barycenter + const SimplexId Id2 = std::get<1>(t); + if(Id2 < 0) { + + if(Id1 < 0) { + continue; + } else { + if(doOptimizeAtoms && CreationFeatures_ && ProgApproach_) { + const PersistencePair &t2 = newData[indexDataCritType[Id1]]; + const double birthData = t2.birth.sfValue; + const double deathData = t2.death.sfValue; + const double birthDeathBarycenter + = birthData + (deathData - birthData) / 2.; + std::array direction; + direction[0] = birthData - birthDeathBarycenter; + direction[1] = deathData - birthDeathBarycenter; + + std::vector> newPairs(matchingsAtoms.size()); + for(size_t j = 0; j < matchingsAtoms.size(); ++j) { + std::array pair{ + birthDeathBarycenter, birthDeathBarycenter}; + newPairs[j] = pair; + } + // pair to add later from the diagonal + pairToAddGradList.push_back(newPairs); + dataAssigned.push_back({birthData, deathData}); + directions.push_back(direction); + } else { + continue; + } + } + } else { + const PersistencePair &t3 = Barycenter[indexBaryCritType[Id2]]; + const double birthBarycenter = t3.birth.sfValue; + const double deathBarycenter = t3.death.sfValue; + auto &direction = directions[indexBaryCritType[Id2]]; + if(Id1 < 0) { + // If matching on the diagonal + const double birthDeathData + = birthBarycenter + (deathBarycenter - birthBarycenter) / 2.; + direction[0] = birthDeathData - birthBarycenter; + direction[1] = birthDeathData - deathBarycenter; + dataAssigned[indexBaryCritType[Id2]] = {birthDeathData, birthDeathData}; + + } else { + const PersistencePair &t2 = newData[indexDataCritType[Id1]]; + const double birthData = t2.birth.sfValue; + const double deathData = t2.death.sfValue; + direction[0] = birthData - birthBarycenter; + direction[1] = deathData - deathBarycenter; + dataAssigned[indexBaryCritType[Id2]] = {birthData, deathData}; + } + tracker2[indexBaryCritType[Id2]] = 1; + } + } +} + +void PersistenceDiagramDictionary::computeDirectionsGradAtoms( + std::vector &gradsAtoms, + const ttk::DiagramType &Barycenter, + const std::vector &weights, + const ttk::DiagramType &newData, + const std::vector &matchingsCritType, + const std::vector &indexBaryCritType, + const std::vector &indexDataCritType, + std::vector>> &pairToAddGradList, + std::vector> &directions, + std::vector &checker, + std::vector &infoToAdd, + const bool doDimReduct) const { + + for(size_t i = 0; i < matchingsCritType.size(); ++i) { + const ttk::MatchingType &t = matchingsCritType[i]; + // Id in newData + const SimplexId Id1 = std::get<0>(t); + // Id in barycenter + const SimplexId Id2 = std::get<1>(t); + int maxWeights + = std::max_element(weights.begin(), weights.end()) - weights.begin(); + if(Id2 < 0) { + if(Id1 < 0) { + continue; + } else { + if(CreationFeatures_) { + const PersistencePair &t2 = newData[indexDataCritType[Id1]]; + const double birthData = t2.birth.sfValue; + const double deathData = t2.death.sfValue; + const double birthDeathBarycenter + = birthData + (deathData - birthData) / 2.; + std::vector direction(2); + direction[0] = birthData - birthDeathBarycenter; + direction[1] = deathData - birthDeathBarycenter; + std::vector> temp3(weights.size()); + std::vector temp2(2); + if(CompressionMode_ && !doDimReduct) { + for(size_t j = 0; j < weights.size(); ++j) { + if(j == static_cast(maxWeights)) { + temp2[0] = -2 * weights[j] * direction[0]; + temp2[1] = -2 * weights[j] * direction[1]; + } else { + temp2[0] = 0.; + temp2[1] = 0.; + } + temp3[j] = temp2; + } + } else { + for(size_t j = 0; j < weights.size(); ++j) { + temp2[0] = -2 * weights[j] * direction[0]; + temp2[1] = -2 * weights[j] * direction[1]; + temp3[j] = temp2; + } + } + gradsAtoms.push_back(temp3); + checker.push_back(1); + std::vector> newPairs(weights.size()); + for(size_t j = 0; j < weights.size(); ++j) { + std::array pair{ + birthDeathBarycenter, birthDeathBarycenter}; + newPairs[j] = pair; + } + // pair to add later from the diagonal + pairToAddGradList.push_back(newPairs); + infoToAdd.push_back(t2); + + } else { + continue; + } + } + } else { + const PersistencePair &t3 = Barycenter[indexBaryCritType[Id2]]; + const double birthBarycenter = t3.birth.sfValue; + const double deathBarycenter = t3.death.sfValue; + std::vector direction(2); + if(Id1 < 0) { + const double birthDeathData + = birthBarycenter + (deathBarycenter - birthBarycenter) / 2.; + direction[0] = birthDeathData - birthBarycenter; + direction[1] = birthDeathData - deathBarycenter; + + } else { + const PersistencePair &t2 = newData[indexDataCritType[Id1]]; + const double birthData = t2.birth.sfValue; + const double deathData = t2.death.sfValue; + direction[0] = birthData - birthBarycenter; + direction[1] = deathData - deathBarycenter; + } + directions[indexBaryCritType[Id2]] = std::move(direction); + checker[indexBaryCritType[Id2]] = 1; + } + } +} + +void PersistenceDiagramDictionary::computeAllDistances( + std::vector &barycentersList, + const size_t nDiags, + std::vector &barycentersListMin, + std::vector &barycentersListSad, + std::vector &barycentersListMax, + std::vector &bidderBarycentersListMin, + std::vector &bidderBarycentersListSad, + std::vector &bidderBarycentersListMax, + std::vector> &originIndexBarysMin, + std::vector> &originIndexBarysSad, + std::vector> &originIndexBarysMax, + std::vector &bidderDiagramsMin, + std::vector &bidderDiagramsMax, + std::vector &bidderDiagramsSad, + std::vector> &matchingsDatasMin, + std::vector> &matchingsDatasMax, + std::vector> &matchingsDatasSad, + std::vector &allLossesAtEpoch, + bool firstDistComputation) const { + + gettingBidderDiagrams(barycentersList, barycentersListMin, barycentersListSad, + barycentersListMax, bidderBarycentersListMin, + bidderBarycentersListSad, bidderBarycentersListMax, + originIndexBarysMin, originIndexBarysSad, + originIndexBarysMax, true); + + // Compute distance and matchings +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(size_t i = 0; i < nDiags; ++i) { + std::vector matchingMin; + std::vector matchingSad; + std::vector matchingMax; + + if(this->do_min_) { + auto &barycentermin = bidderBarycentersListMin[i]; + auto &datamin = bidderDiagramsMin[i]; + + if(firstDistComputation) { + allLossesAtEpoch[i] + += computeDistance(datamin, barycentermin, matchingMin); + + } else { + computeDistance(datamin, barycentermin, matchingMin); + } + } + if(this->do_max_) { + auto &barycentermax = bidderBarycentersListMax[i]; + auto &datamax = bidderDiagramsMax[i]; + + if(firstDistComputation) { + allLossesAtEpoch[i] + += computeDistance(datamax, barycentermax, matchingMax); + + } else { + computeDistance(datamax, barycentermax, matchingMax); + } + } + if(this->do_sad_) { + auto &barycentersListad = bidderBarycentersListSad[i]; + auto &datasad = bidderDiagramsSad[i]; + + if(firstDistComputation) { + allLossesAtEpoch[i] + += computeDistance(datasad, barycentersListad, matchingSad); + + } else { + computeDistance(datasad, barycentersListad, matchingSad); + } + } + matchingsDatasMin[i] = std::move(matchingMin); + matchingsDatasSad[i] = std::move(matchingSad); + matchingsDatasMax[i] = std::move(matchingMax); + } +} \ No newline at end of file diff --git a/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.h b/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.h new file mode 100644 index 0000000000..6c449ddb9f --- /dev/null +++ b/core/base/persistenceDiagramDictionary/PersistenceDiagramDictionary.h @@ -0,0 +1,245 @@ +/// \ingroup base +/// \class ttk::PersistenceDiagramDictionary +/// \author Keanu Sisouk +/// \author Pierre Guillou +/// \date Mai 2023 +/// +/// \b Related \b publication \n +/// "Wasserstein Dictionaries of Persistence Diagrams" \n +/// Keanu Sisouk, Julie Delon and Julien Tierny \n +/// IEEE Transactions on Visualization and Computer Graphics, 2023. +/// +/// \sa PersistenceDiagramDictionary + +#pragma once + +#include "PersistenceDiagramUtils.h" +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace ttk { + class PersistenceDiagramDictionary : virtual public Debug { + + public: + PersistenceDiagramDictionary() { + this->setDebugMsgPrefix("PersistenceDiagramDictionary"); + } + + void execute(std::vector &intermediateDiagrams, + const std::vector &intermediateAtoms, + std::vector &dictDiagrams, + std::vector> &vectorWeights, + const int seed, + const int numAtom, + std::vector &lossTab, + std::vector> &allLosses, + double percent); + + void method(const std::vector &intermediateDiagrams, + std::vector &dictDiagrams, + std::vector> &vectorWeights, + const int numAtom, + std::vector &lossTab, + std::vector> &allLosses, + std::vector> &histoVectorWeights, + std::vector &histoDictDiagrams, + bool preWeightOpt, + double percent, + bool doCompression); + + enum class BACKEND { + BORDER_INIT = 0, + RANDOM_INIT = 1, + FIRST_DIAGS = 2, + INPUT_ATOMS = 3, + GREEDY_INIT = 4 + }; + + inline void setUseDimReduct(bool data) { + DimReductMode_ = data; + } + + inline void setUseProgApproach(bool data) { + ProgApproach_ = data; + } + + inline void setDos(const bool min, const bool sad, const bool max) { + do_min_ = min; + do_sad_ = sad; + do_max_ = max; + } + + inline void setMinPersistence_(const double data) { + MinPersistence_ = data; + } + + protected: + BACKEND BackEnd{BACKEND::BORDER_INIT}; + double distVect(const std::vector &vec1, + const std::vector &vec2) const; + + double + getMostPersistent(const std::vector &bidder_diags) const; + double computeDistance(const BidderDiagram &D1, + const BidderDiagram &D2, + std::vector &matching) const; + + void computeGradientWeights( + std::vector &gradWeights, + std::vector &hessianList, + const std::vector &dictDiagrams, + const std::vector> &matchingsAtoms, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsMin, + const std::vector &matchingsMax, + const std::vector &matchingsSad, + const std::vector &indexBaryMin, + const std::vector &indexBaryMax, + const std::vector &indexBarySad, + const std::vector &indexDataMin, + const std::vector &indexDataMax, + const std::vector &indexDataSad, + const bool doOptimizeAtoms) const; + + void computeGradientAtoms( + std::vector &gradsAtoms, + const std::vector &weights, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsMin, + const std::vector &matchingsMax, + const std::vector &matchingsSad, + const std::vector &indexBaryMin, + const std::vector &indexBaryMax, + const std::vector &indexBarySad, + const std::vector &indexDataMin, + const std::vector &indexDataMax, + const std::vector &indexDataSad, + std::vector &checker, + std::vector>> &pairToAddGradList, + ttk::DiagramType &infoToAdd, + bool doDimReduct) const; + + void computeDirectionsGradWeight( + const std::vector> &matchingsAtoms, + const ttk::DiagramType &Barycenter, + const ttk::DiagramType &newData, + const std::vector &matchingsCritType, + const std::vector &indexBaryCritType, + const std::vector &indexDataCritType, + std::vector>> &pairToAddGradList, + std::vector> &directions, + std::vector> &data_assigned, + std::vector &tracker2, + const bool doOptimizeAtoms) const; + + void computeDirectionsGradAtoms( + std::vector &gradsAtoms, + const ttk::DiagramType &Barycenter, + const std::vector &weights, + const ttk::DiagramType &newData, + const std::vector &matchingsCritType, + const std::vector &indexBaryCritType, + const std::vector &indexDataCritType, + std::vector>> &pairToAddGradList, + std::vector> &directions, + std::vector &checker, + std::vector &infoToAdd, + const bool doOptimizeAtoms) const; + + void setBidderDiagrams(const size_t nInputs, + std::vector &inputDiagrams, + std::vector &bidder_diags) const; + + int initDictionary(std::vector &dictDiagrams, + const std::vector &datas, + const std::vector &inputAtoms, + const int nbAtom, + bool do_min_, + bool do_sad_, + bool do_max_, + int seed, + double percent); + + void gettingBidderDiagrams( + const std::vector &intermediateDiagrams, + std::vector &inputDiagramsMin, + std::vector &inputDiagramsSad, + std::vector &inputDiagramsMax, + std::vector &bidderDiagramsMin, + std::vector &bidderDiagramsSad, + std::vector &bidderDiagramsMax, + std::vector> &originIndexMin, + std::vector> &originIndexSad, + std::vector> &originIndexMax, + bool insertOriginIndexMode) const; + + void computeAllDistances( + std::vector &barycentersList, + const size_t nDiag, + std::vector &barycentersListMin, + std::vector &barycentersListSad, + std::vector &barycentersListMax, + std::vector &bidderBarycentersListMin, + std::vector &bidderBarycentersListSad, + std::vector &bidderBarycentersListMax, + std::vector> &originIndexBarysMin, + std::vector> &originIndexBarysSad, + std::vector> &originIndexBarysMax, + std::vector &bidderDiagramsMin, + std::vector &bidderDiagramsMax, + std::vector &bidderDiagramsSad, + std::vector> &matchingsDatasMin, + std::vector> &matchingsDatasMax, + std::vector> &matchingsDatasSad, + std::vector &allLossesAtEpoch, + bool firstDistComputation) const; + + void controlAtomsSize( + const std::vector &intermediateDiagrams, + std::vector &dictDiagrams) const; + + double getMaxPers(const ttk::DiagramType &data); + + int Wasserstein{2}; + double Alpha{1.0}; + double DeltaLim{0.01}; + // lambda : 0<=lambda<=1 + // parametrizes the point used for the physical (critical) coordinates of + // the persistence paired lambda = 1 : extremum (min if pair min-sad, max if + // pair sad-max) lambda = 0 : saddle (bad stability) lambda = 1/2 : middle + // of the 2 critical points of the pair + double Lambda{1.0}; + size_t MaxNumberOfPairs{20}; + double MinPersistence_{0.1}; + + double CompressionFactor{1.5}; + bool do_min_{true}, do_sad_{true}, do_max_{true}; + + int maxLag2_; + + int MaxEpoch_; + bool MaxEigenValue_{true}; + bool OptimizeWeights_{true}; + bool OptimizeAtoms_{true}; + + bool CreationFeatures_{true}; + bool Fusion_{false}; + bool ProgBarycenter_{false}; + + bool sortedForTest_{false}; + bool ProgApproach_{false}; + bool StopCondition_{true}; + + bool CompressionMode_{false}; + bool DimReductMode_{false}; + }; +} // namespace ttk diff --git a/core/base/persistenceDiagramDictionaryDecoding/CMakeLists.txt b/core/base/persistenceDiagramDictionaryDecoding/CMakeLists.txt new file mode 100644 index 0000000000..b166bd7d77 --- /dev/null +++ b/core/base/persistenceDiagramDictionaryDecoding/CMakeLists.txt @@ -0,0 +1,13 @@ +ttk_add_base_library(persistenceDiagramDictionaryDecoding + SOURCES + PersistenceDiagramDictionaryDecoding.cpp + HEADERS + PersistenceDiagramDictionaryDecoding.h + DEPENDS + common + persistenceDiagram + persistenceDiagramClustering + persistenceDiagramDistanceMatrix + persistenceDiagramDictionary + dimensionReduction + ) diff --git a/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.cpp b/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.cpp new file mode 100644 index 0000000000..8e538184c5 --- /dev/null +++ b/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.cpp @@ -0,0 +1,279 @@ +#include "DimensionReduction.h" +#include "PersistenceDiagramUtils.h" +#include +#include + +void ttk::PersistenceDiagramDictionaryDecoding::execute( + std::vector &dictDiagrams, + std::vector> &vectorWeights, + std::vector &Barycenters) const { + + Timer tm{}; + + std::vector>> AllMatchingsAtoms( + Barycenters.size()); + + for(size_t i = 0; i < Barycenters.size(); ++i) { + auto &barycenter = Barycenters[i]; + auto &weight = vectorWeights[i]; + auto &matchings = AllMatchingsAtoms[i]; + computeWeightedBarycenter( + dictDiagrams, weight, barycenter, matchings, *this, ProgBarycenter); + } + + this->printMsg( + "Computed barycenters", 1.0, tm.getElapsedTime(), this->threadNumber_); +} + +void ttk::PersistenceDiagramDictionaryDecoding::computeAtomsCoordinates( + std::vector &atoms, + const std::vector> &vectorWeights, + std::vector> &coords, + std::vector> &trueCoords, + std::vector &xVector, + std::vector &yVector, + std::vector &zVector, + const double spacing, + const size_t nAtoms) const { + + if(nAtoms == 2) { + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nAtoms, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(2); + const auto distMatrix = MatrixCalculator.execute(atoms, nInputs); + coords[0][0] = 0.; + trueCoords[0][0] = 0.; + coords[0][1] = 0.; + trueCoords[0][0] = 0.; + coords[1][0] = spacing * distMatrix[0][1]; + trueCoords[1][0] = distMatrix[0][1]; + trueCoords[1][1] = 0.; + + } else if(nAtoms == 3) { + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nAtoms, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(3); + std::vector> distMatrix + = MatrixCalculator.execute(atoms, nInputs); + coords[0][0] = 0.; + trueCoords[0][0] = 0.; + coords[0][1] = 0.; + trueCoords[0][1] = 0.; + coords[1][0] = spacing * distMatrix[0][1]; + trueCoords[1][0] = distMatrix[0][1]; + coords[1][1] = 0.; + trueCoords[1][1] = 0.; + double distOpposed = distMatrix[2][1]; + double firstDist = distMatrix[0][1]; + double distAdja = distMatrix[0][2]; + double alpha = std::acos( + (distOpposed * distOpposed - firstDist * firstDist - distAdja * distAdja) + / (-2. * firstDist * distAdja)); + coords[2][0] = spacing * distAdja * std::cos(alpha); + trueCoords[2][0] = distAdja * std::cos(alpha); + coords[2][1] = spacing * distAdja * std::sin(alpha); + trueCoords[2][1] = distAdja * std::sin(alpha); + + } else if(nAtoms == 4) { + switch(this->ProjMet) { + case BACKEND::DICTIONARY: { + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nAtoms, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(3); + std::vector> distMatrix + = MatrixCalculator.execute(atoms, nInputs); + coords[0][0] = 0.; + trueCoords[0][0] = 0.; + coords[0][1] = 0.; + trueCoords[0][1] = 0.; + coords[1][0] = spacing * distMatrix[0][1]; + trueCoords[1][0] = distMatrix[0][1]; + coords[1][1] = 0.; + trueCoords[1][1] = 0.; + double distOpposed = distMatrix[2][1]; + double firstDist = distMatrix[0][1]; + double distAdja = distMatrix[0][2]; + double alpha = std::acos((distOpposed * distOpposed + - firstDist * firstDist - distAdja * distAdja) + / (-2. * firstDist * distAdja)); + coords[2][0] = spacing * distAdja * std::cos(alpha); + trueCoords[2][0] = distAdja * std::cos(alpha); + coords[2][1] = spacing * distAdja * std::sin(alpha); + trueCoords[2][1] = distAdja * std::sin(alpha); + double firstHeight = distMatrix[3][0]; + double secondHeight = distMatrix[3][1]; + double thirdHeight = distMatrix[3][2]; + trueCoords[3][0] = (pow(firstHeight, 2) - pow(secondHeight, 2) + + pow(trueCoords[1][0], 2)) + / (2 * trueCoords[1][0]); + trueCoords[3][1] + = (pow(firstHeight, 2) - pow(thirdHeight, 2) + + pow(trueCoords[2][0], 2) + pow(trueCoords[2][1], 2) + - 2 * trueCoords[3][0] * trueCoords[2][0]) + / (2 * trueCoords[2][1]); + trueCoords[3][2] + = std::sqrt(pow(firstHeight, 2) - pow(trueCoords[3][0], 2) + - pow(trueCoords[3][1], 2)); + + break; + } + + case BACKEND::MDS: { + ttk::DimensionReduction DimProjector; + DimProjector.setIsInputDistanceMatrix(true); + DimProjector.setInputNumberOfComponents(3); + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nAtoms, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(4); + std::vector> distMatrix + = MatrixCalculator.execute(atoms, nInputs); + int nRow = distMatrix.size(); + std::vector matrixForProjector; + for(int i = 0; i < nRow; ++i) { + for(int j = 0; j < nRow; ++j) { + matrixForProjector.push_back(distMatrix[j][i]); + } + } + std::vector> coordsAtom; + DimProjector.execute(coordsAtom, matrixForProjector, nRow, nRow); + + for(size_t i = 0; i < 3; ++i) { + for(size_t j = 0; j < nAtoms; ++j) { + if(i == 0) { + trueCoords[j][0] = coordsAtom[0][j]; + } else if(i == 1) { + trueCoords[j][1] = coordsAtom[1][j]; + } else { + trueCoords[j][2] = coordsAtom[2][j]; + } + } + } + + break; + } + } + } else { + switch(this->ProjMet) { + case BACKEND::DICTIONARY: { + + std::vector dictDiagrams; + std::vector intermediateAtoms; + std::vector lossTab; + std::vector> allLosses(nAtoms); + const int seed = 0; + const int m = 3; + std::vector> tempWeights(nAtoms); + for(size_t i = 0; i < tempWeights.size(); ++i) { + std::vector weights(m, 1. / (m * 1.)); + tempWeights[i] = std::move(weights); + } + ttk::PersistenceDiagramDictionary DictionaryEncoder; + DictionaryEncoder.setUseDimReduct(false); + DictionaryEncoder.setUseProgApproach(true); + DictionaryEncoder.execute(atoms, atoms, dictDiagrams, tempWeights, seed, + m, lossTab, allLosses, 0.); + std::vector> tempCoords(3); + std::vector> tempTrueCoords(3); + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{3, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(3); + std::vector> distMatrix + = MatrixCalculator.execute(dictDiagrams, nInputs); + tempCoords[0][0] = 0.; + tempTrueCoords[0][0] = 0.; + tempCoords[0][1] = 0.; + tempTrueCoords[0][1] = 0.; + tempCoords[1][0] = spacing * distMatrix[0][1]; + tempTrueCoords[1][0] = distMatrix[0][1]; + tempCoords[1][1] = 0.; + tempTrueCoords[0][1] = 0.; + double distOpposed = distMatrix[2][1]; + double firstDist = distMatrix[0][1]; + double distAdja = distMatrix[0][2]; + double alpha = std::acos((distOpposed * distOpposed + - firstDist * firstDist - distAdja * distAdja) + / (-2. * firstDist * distAdja)); + tempCoords[2][0] = spacing * distAdja * std::cos(alpha); + tempTrueCoords[2][0] = distAdja * std::cos(alpha); + tempCoords[2][1] = spacing * distAdja * std::sin(alpha); + tempTrueCoords[2][1] = distAdja * std::sin(alpha); + for(int i = 0; i < 2; ++i) { + for(size_t j = 0; j < nAtoms; ++j) { + double temp = 0.; + for(int iAtom = 0; iAtom < 3; ++iAtom) { + if(i == 0) { + temp += tempWeights[j][iAtom] * tempTrueCoords[iAtom][0]; + trueCoords[j][0] = temp; + } else { + temp += tempWeights[j][iAtom] * tempTrueCoords[iAtom][1]; + trueCoords[j][1] = temp; + } + } + } + } + + break; + } + + case BACKEND::MDS: { + ttk::DimensionReduction DimProjector; + DimProjector.setIsInputDistanceMatrix(true); + ttk::PersistenceDiagramDistanceMatrix MatrixCalculator; + std::array nInputs{nAtoms, 0}; + MatrixCalculator.setDos(true, true, true); + MatrixCalculator.setThreadNumber(3); + std::vector> distMatrix + = MatrixCalculator.execute(atoms, nInputs); + int nRow = distMatrix.size(); + std::vector matrixForProjector; + for(int i = 0; i < nRow; ++i) { + for(int j = 0; j < nRow; ++j) { + matrixForProjector.push_back(distMatrix[j][i]); + } + } + std::vector> coordsAtom; + DimProjector.execute(coordsAtom, matrixForProjector, nRow, nRow); + + for(size_t i = 0; i < 2; ++i) { + for(size_t j = 0; j < nAtoms; ++j) { + if(i == 0) { + trueCoords[j][0] = coordsAtom[0][j]; + } else { + trueCoords[j][1] = coordsAtom[1][j]; + } + } + } + + break; + } + } + } + size_t nDiags = vectorWeights.size(); + + for(int i = 0; i < 3; ++i) { + for(size_t j = 0; j < nDiags; ++j) { + double temp = 0.; + for(size_t iAtom = 0; iAtom < nAtoms; ++iAtom) { + if(i == 0) { + temp += vectorWeights[j][iAtom] * trueCoords[iAtom][0]; + } else if(i == 1) { + temp += vectorWeights[j][iAtom] * trueCoords[iAtom][1]; + } else { + temp += vectorWeights[j][iAtom] * trueCoords[iAtom][2]; + } + } + if(i == 0) { + xVector[j] = temp; + } else if(i == 1) { + yVector[j] = temp; + } else { + zVector[j] = temp; + } + } + } +} \ No newline at end of file diff --git a/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.h b/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.h new file mode 100644 index 0000000000..2f1b6df03c --- /dev/null +++ b/core/base/persistenceDiagramDictionaryDecoding/PersistenceDiagramDictionaryDecoding.h @@ -0,0 +1,67 @@ +/// \ingroup base +/// \class ttk::PersistenceDiagramDictionaryDecoding +/// \author Keanu Sisouk +/// \author Pierre Guillou +/// \date Mai 2023 +/// +/// \brief TTK processing package for the computation of a Dictionary +/// of Persistence Diagrams and barycentric weights to approximate +/// an ensemble of Persistence Diagrams. +/// +/// \b Related \b publication \n +/// "Wasserstein Dictionaries of Persistence Diagrams" \n +/// Keanu Sisouk, Julie Delon and Julien Tierny \n +/// IEEE Transactions on Visualization and Computer Graphics, 2023. +/// +/// \sa PersistenceDiagramDictionary +#pragma once + +// ttk common includes +#include +#include +#include +#include +#include + +#include + +namespace ttk { + using Matrice = std::vector>; + using VectorMatchingTuple = std::vector; + + /** + * The PersistenceDiagramDictionaryDecoding class provides methods to compute + * for each vertex of a triangulation the average scalar value of itself and + * its direct neighbors. + */ + class PersistenceDiagramDictionaryDecoding : virtual public Debug { + + public: + enum class BACKEND { MDS = 0, DICTIONARY = 1 }; + + PersistenceDiagramDictionaryDecoding() { + this->setDebugMsgPrefix("PersistenceDiagramDictionaryDecoding"); + } + + void execute(std::vector &dictDiagrams, + std::vector> &vectorWeights, + std::vector &Barycenters) const; + + protected: + BACKEND ProjMet{BACKEND::MDS}; + bool ProgBarycenter{false}; + + void computeAtomsCoordinates( + std::vector &atoms, + const std::vector> &vectorWeights, + std::vector> &coords, + std::vector> &trueCoords, + std::vector &xVector, + std::vector &yVector, + std::vector &zVector, + const double spacing, + const size_t nAtoms) const; + + }; // PersistenceDiagramDictionaryDecoding class + +} // namespace ttk diff --git a/core/vtk/ttkPersistenceDiagramDictionary/CMakeLists.txt b/core/vtk/ttkPersistenceDiagramDictionary/CMakeLists.txt new file mode 100644 index 0000000000..73e96de098 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionary/CMakeLists.txt @@ -0,0 +1 @@ +ttk_add_vtk_module() diff --git a/core/vtk/ttkPersistenceDiagramDictionary/ttk.module b/core/vtk/ttkPersistenceDiagramDictionary/ttk.module new file mode 100644 index 0000000000..18692ace94 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionary/ttk.module @@ -0,0 +1,10 @@ +NAME + ttkPersistenceDiagramDictionary +SOURCES + ttkPersistenceDiagramDictionary.cpp +HEADERS + ttkPersistenceDiagramDictionary.h +DEPENDS + persistenceDiagramDictionary + ttkAlgorithm + ttkPersistenceDiagram diff --git a/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.cpp b/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.cpp new file mode 100644 index 0000000000..22414453c2 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.cpp @@ -0,0 +1,221 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(ttkPersistenceDiagramDictionary); + +ttkPersistenceDiagramDictionary::ttkPersistenceDiagramDictionary() { + SetNumberOfInputPorts(2); + SetNumberOfOutputPorts(2); +} + +int ttkPersistenceDiagramDictionary::FillInputPortInformation( + int port, vtkInformation *info) { + if(port == 0) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkMultiBlockDataSet"); + return 1; + } else if(port == 1) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkMultiBlockDataSet"); + info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1); + return 1; + } + return 0; +} + +int ttkPersistenceDiagramDictionary::FillOutputPortInformation( + int port, vtkInformation *info) { + if(port == 0) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkMultiBlockDataSet"); + return 1; + } else if(port == 1) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable"); + return 1; + } else { + return 0; + } +} + +// to adapt if your wrapper does not inherit from vtkDataSetAlgorithm +int ttkPersistenceDiagramDictionary::RequestData( + vtkInformation * /*request*/, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + ttk::Memory m; + + // Get input data + auto blocks = vtkMultiBlockDataSet::GetData(inputVector[0], 0); + auto atomBlocks = vtkMultiBlockDataSet::GetData(inputVector[1], 0); + + // Flat storage for diagrams extracted from blocks + std::vector inputDiagrams; + + std::vector inputAtoms; + + int numInputAtoms = 0; + if(atomBlocks != nullptr) { + numInputAtoms = atomBlocks->GetNumberOfBlocks(); + inputAtoms.resize(numInputAtoms); + for(int i = 0; i < numInputAtoms; ++i) { + inputAtoms[i] + = vtkUnstructuredGrid::SafeDownCast(atomBlocks->GetBlock(i)); + } + } + if(BackEnd == BACKEND::INPUT_ATOMS) { + AtomNumber_ = numInputAtoms; + } + if(blocks != nullptr) { + int numInputs = blocks->GetNumberOfBlocks(); + inputDiagrams.resize(numInputs); + for(int i = 0; i < numInputs; ++i) { + inputDiagrams[i] = vtkUnstructuredGrid::SafeDownCast(blocks->GetBlock(i)); + } + } + + const int numAtom = this->GetAtomNumber_(); + printMsg("Number of atoms: " + ttk::debug::output::YELLOW + + ttk::debug::output::UNDERLINED + std::to_string(numAtom) + + ttk::debug::output::ENDCOLOR + ttk::debug::output::ENDCOLOR); + + // total number of diagrams + const int nDiags = inputDiagrams.size(); + + // Sanity check + for(const auto vtu : inputDiagrams) { + if(vtu == nullptr) { + this->printErr("Input diagrams are not all vtkUnstructuredGrid"); + return 0; + } + } + + // Set output + auto outputDgm = vtkMultiBlockDataSet::GetData(outputVector, 0); + auto outputWeights = vtkTable::GetData(outputVector, 1); + + outputDgm->SetNumberOfBlocks(numAtom); + + if(BackEnd == BACKEND::INPUT_ATOMS) { + for(int i = 0; i < numAtom; ++i) { + vtkNew vtu; + vtu->DeepCopy(inputAtoms[i]); + outputDgm->SetBlock(i, vtu); + } + } else { + for(int i = 0; i < numAtom; ++i) { + vtkNew vtu; + vtu->DeepCopy(inputDiagrams[i]); + outputDgm->SetBlock(i, vtu); + } + } + + std::vector intermediateDiagrams(nDiags); + std::vector intermediateAtoms(numInputAtoms); + + double max_dimension_total = 0.0; + for(int i = 0; i < nDiags; ++i) { + + const auto ret + = VTUToDiagram(intermediateDiagrams[i], inputDiagrams[i], *this); + + double maxPers = this->getMaxPers(intermediateDiagrams[i]); + double percentage = this->Percent_; + + intermediateDiagrams[i].erase( + std::remove_if(intermediateDiagrams[i].begin(), + intermediateDiagrams[i].end(), + [maxPers, percentage](ttk::PersistencePair &t) { + return (t.death.sfValue - t.birth.sfValue) + < (percentage / 100.) * maxPers; + }), + intermediateDiagrams[i].end()); + + if(ret != 0) { + this->printErr("Could not read Persistence Diagram"); + return 0; + } + if(max_dimension_total < maxPers) { + max_dimension_total = maxPers; + } + } + for(int i = 0; i < numInputAtoms; ++i) { + const auto ret = VTUToDiagram(intermediateAtoms[i], inputAtoms[i], *this); + if(ret != 0) { + this->printErr("Could not read Persistence Diagram"); + return 0; + } + } + + std::vector dictDiagrams; + const int seed = this->GetSeed_(); + + std::vector> vectorWeights(nDiags); + for(size_t i = 0; i < vectorWeights.size(); ++i) { + std::vector weights(numAtom, 1. / (numAtom * 1.)); + vectorWeights[i] = std::move(weights); + } + + std::vector lossTab; + std::vector> allLosses(nDiags); + this->execute(intermediateDiagrams, intermediateAtoms, dictDiagrams, + vectorWeights, seed, numAtom, lossTab, allLosses, + this->Percent_); + // zero-padd column name to keep Row Data columns ordered + outputWeights->SetNumberOfRows(nDiags); + + const auto zeroPad + = [](std::string &colName, const size_t numberCols, const size_t colIdx) { + std::string max{std::to_string(numberCols - 1)}; + std::string cur{std::to_string(colIdx)}; + std::string zer(max.size() - cur.size(), '0'); + colName.append(zer).append(cur); + }; + for(int i = 0; i < numAtom; ++i) { + std::string name{"Atom"}; + zeroPad(name, numAtom, i); + // name + vtkNew col{}; + col->SetNumberOfValues(nDiags); + col->SetName(name.c_str()); + for(int j = 0; j < nDiags; ++j) { + col->SetValue(j, vectorWeights[j][i]); + } + col->Modified(); + outputWeights->AddColumn(col); + } + + vtkNew fd{}; + fd->CopyStructure(inputDiagrams[0]->GetFieldData()); + fd->SetNumberOfTuples(nDiags); + for(int i = 0; i < nDiags; ++i) { + fd->SetTuple(i, 0, inputDiagrams[i]->GetFieldData()); + } + + for(int i = 0; i < fd->GetNumberOfArrays(); ++i) { + outputWeights->AddColumn(fd->GetAbstractArray(i)); + } + + vtkNew dummy{}; + + for(int i = 0; i < numAtom; ++i) { + vtkNew vtu; + ttk::DiagramType &diagram = dictDiagrams[i]; + DiagramToVTU(vtu, diagram, dummy, *this, 3, false); + outputDgm->SetBlock(i, vtu); + } + return 1; +} diff --git a/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.h b/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.h new file mode 100644 index 0000000000..7eacf0c271 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionary/ttkPersistenceDiagramDictionary.h @@ -0,0 +1,177 @@ +/// \ingroup vtk +/// \class ttkPersistenceDiagramDictionary +/// \author Keanu Sisouk +/// \author Pierre Guillou +/// \date Mai 2023 +/// +/// \brief TTK processing package for the computation of a Dictionary +/// of Persistence Diagrams and barycentric weights to approximate +/// an ensemble of Persistence Diagrams. +/// +/// \b Related \b publication \n +/// "Wasserstein Dictionaries of Persistence Diagrams" \n +/// Keanu Sisouk, Julie Delon and Julien Tierny \n +/// IEEE Transactions on Visualization and Computer Graphics, 2023. +/// +/// \sa PersistenceDiagramDictionary + +#pragma once + +// VTK includes +#include +#include +#include +#include + +// VTK Module +#include + +// ttk code includes + +#include +#include +#include + +class TTKPERSISTENCEDIAGRAMDICTIONARY_EXPORT ttkPersistenceDiagramDictionary + : public ttkAlgorithm, + protected ttk::PersistenceDiagramDictionary { + +private: + int AtomNumber_{3}; + int Seed_{0}; + double Percent_{0}; + +public: + // enum class BACKEND{BORDER_INIT = 0 , RANDOM_INIT = 1 , FIRST_DIAGS = 2}; + + static ttkPersistenceDiagramDictionary *New(); + + vtkTypeMacro(ttkPersistenceDiagramDictionary, ttkAlgorithm); + + void SetWassersteinMetric(const std::string &data) { + Wasserstein = (data == "inf") ? -1 : stoi(data); + Modified(); + } + std::string GetWassersteinMetric() { + return Wasserstein == -1 ? "inf" : std::to_string(Wasserstein); + } + + void SetAntiAlpha(double data) { + data = 1 - data; + if(data > 0 && data <= 1) { + Alpha = data; + } else if(data > 1) { + Alpha = 1; + } else { + Alpha = 0.001; + } + Modified(); + } + + vtkGetMacro(Alpha, double); + + vtkSetMacro(Percent_, double); + vtkGetMacro(Percent_, double); + + vtkSetMacro(OptimizeWeights_, int); + vtkGetMacro(OptimizeWeights_, int); + + vtkSetMacro(OptimizeAtoms_, int); + vtkGetMacro(OptimizeAtoms_, int); + + vtkSetMacro(MaxEigenValue_, int); + vtkGetMacro(MaxEigenValue_, int); + + vtkSetMacro(Fusion_, int); + vtkGetMacro(Fusion_, int); + + vtkSetMacro(ProgBarycenter_, int); + vtkGetMacro(ProgBarycenter_, int); + + vtkSetMacro(MaxEpoch_, int); + vtkGetMacro(MaxEpoch_, int); + + vtkSetMacro(ProgApproach_, int); + vtkGetMacro(ProgApproach_, int); + + vtkSetMacro(StopCondition_, int); + vtkGetMacro(StopCondition_, int); + + vtkSetMacro(CompressionMode_, int); + vtkGetMacro(CompressionMode_, int); + + vtkSetMacro(DimReductMode_, int); + vtkGetMacro(DimReductMode_, int); + + vtkSetMacro(sortedForTest_, int); + vtkGetMacro(sortedForTest_, int); + + vtkSetMacro(CreationFeatures_, int); + vtkGetMacro(CreationFeatures_, int); + + vtkSetMacro(AtomNumber_, int); + vtkGetMacro(AtomNumber_, int); + + vtkSetMacro(Seed_, int); + vtkGetMacro(Seed_, int); + + vtkSetMacro(DeltaLim, double); + vtkGetMacro(DeltaLim, double); + + vtkSetMacro(Lambda, double); + vtkGetMacro(Lambda, double); + + ttkSetEnumMacro(BackEnd, BACKEND); + vtkGetEnumMacro(BackEnd, BACKEND); + + vtkSetMacro(CompressionFactor, double); + vtkGetMacro(CompressionFactor, double); + + void SetPairType(const int data) { + switch(data) { + case(0): + this->setDos(true, false, false); + break; + case(1): + this->setDos(false, true, false); + break; + case(2): + this->setDos(false, false, true); + break; + default: + this->setDos(true, true, true); + break; + } + Modified(); + } + int GetPairType() { + if(do_min_ && do_sad_ && do_max_) { + return -1; + } else if(do_min_) { + return 0; + } else if(do_sad_) { + return 1; + } else if(do_max_) { + return 2; + } + return -1; + } + + vtkSetMacro(MaxNumberOfPairs, unsigned int); + vtkGetMacro(MaxNumberOfPairs, unsigned int); + + vtkSetMacro(MinPersistence_, double); + vtkGetMacro(MinPersistence_, double); + +protected: + ttkPersistenceDiagramDictionary(); + ~ttkPersistenceDiagramDictionary() override = default; + + // BACKEND BackEnd{BACKEND::BORDER_INIT}; + int FillInputPortInformation(int port, vtkInformation *info) override; + int FillOutputPortInformation(int port, vtkInformation *info) override; + + int RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; +}; diff --git a/core/vtk/ttkPersistenceDiagramDictionary/vtk.module b/core/vtk/ttkPersistenceDiagramDictionary/vtk.module new file mode 100644 index 0000000000..c19f2074af --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionary/vtk.module @@ -0,0 +1,4 @@ +NAME + ttkPersistenceDiagramDictionary +DEPENDS + ttkAlgorithm diff --git a/core/vtk/ttkPersistenceDiagramDictionaryDecoding/CMakeLists.txt b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/CMakeLists.txt new file mode 100644 index 0000000000..73e96de098 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/CMakeLists.txt @@ -0,0 +1 @@ +ttk_add_vtk_module() diff --git a/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttk.module b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttk.module new file mode 100644 index 0000000000..5db94f0fe0 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttk.module @@ -0,0 +1,10 @@ +NAME + ttkPersistenceDiagramDictionaryDecoding +SOURCES + ttkPersistenceDiagramDictionaryDecoding.cpp +HEADERS + ttkPersistenceDiagramDictionaryDecoding.h +DEPENDS + persistenceDiagramDictionaryDecoding + ttkAlgorithm + ttkPersistenceDiagram diff --git a/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.cpp b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.cpp new file mode 100644 index 0000000000..a3999c356a --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.cpp @@ -0,0 +1,329 @@ +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +vtkStandardNewMacro(ttkPersistenceDiagramDictionaryDecoding); + +ttkPersistenceDiagramDictionaryDecoding:: + ttkPersistenceDiagramDictionaryDecoding() { + this->SetNumberOfInputPorts(2); + this->SetNumberOfOutputPorts(2); +} + +int ttkPersistenceDiagramDictionaryDecoding::FillInputPortInformation( + int port, vtkInformation *info) { + if(port == 0) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkMultiBlockDataSet"); + return 1; + } else if(port == 1) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkTable"); + return 1; + } else { + return 0; + } +} + +int ttkPersistenceDiagramDictionaryDecoding::FillOutputPortInformation( + int port, vtkInformation *info) { + if(port == 0) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkMultiBlockDataSet"); + return 1; + } else if(port == 1) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable"); + return 1; + } else { + return 0; + } +} + +int ttkPersistenceDiagramDictionaryDecoding::RequestData( + vtkInformation * /*request*/, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + + const auto blocks = vtkMultiBlockDataSet::GetData(inputVector[0]); + const auto weightsVTK = vtkTable::GetData(inputVector[1]); + + std::vector inputDiagrams; + + if(blocks != nullptr) { + int numInputs = blocks->GetNumberOfBlocks(); + inputDiagrams.resize(numInputs); + for(int i = 0; i < numInputs; ++i) { + inputDiagrams[i] = vtkUnstructuredGrid::SafeDownCast(blocks->GetBlock(i)); + } + } + + const size_t nDiags = inputDiagrams.size(); + + std::vector dictDiagrams(nDiags); + for(size_t i = 0; i < nDiags; ++i) { + vtkNew vtu; + vtu->DeepCopy(inputDiagrams[i]); + auto &atom = dictDiagrams[i]; + const auto ret = VTUToDiagram(atom, vtu, *this); + if(ret != 0) { + this->printWrn("Could not read Persistence Diagram"); + } + } + + // Sanity check + for(const auto vtu : inputDiagrams) { + if(vtu == nullptr) { + this->printErr("Input diagrams are not all vtkUnstructuredGrid"); + return 0; + } + } + + const auto zeroPad + = [](std::string &colName, const size_t numberCols, const size_t colIdx) { + std::string max{std::to_string(numberCols - 1)}; + std::string cur{std::to_string(colIdx)}; + std::string zer(max.size() - cur.size(), '0'); + colName.append(zer).append(cur); + }; + + vtkNew temp; + temp->DeepCopy(weightsVTK); + std::vector inputWeights; + int numWeights = temp->GetNumberOfRows(); + for(int i = 0; i < temp->GetNumberOfColumns(); ++i) { + std::cout << temp->GetColumnName(i) << "\n"; + } + + if(weightsVTK != nullptr) { + inputWeights.resize(nDiags); + for(size_t i = 0; i < nDiags; ++i) { + std::string name{"Atom"}; + zeroPad(name, nDiags, i); + inputWeights[i] + = vtkDataArray::SafeDownCast(temp->GetColumnByName(name.c_str())); + } + } + + std::vector> vectorWeights(numWeights); + for(int i = 0; i < numWeights; ++i) { + std::vector &t1 = vectorWeights[i]; + for(size_t j = 0; j < nDiags; ++j) { + double weight = inputWeights[j]->GetTuple1(i); + t1.push_back(weight); + } + } + std::vector Barycenters(numWeights); + + if(!ComputePoints) { + this->execute(dictDiagrams, vectorWeights, Barycenters); + } + + auto outputDgm = vtkMultiBlockDataSet::GetData(outputVector, 0); + auto outputCoordinates = vtkTable::GetData(outputVector, 1); + outputDgm->SetNumberOfBlocks(numWeights); + outputCoordinates->SetNumberOfRows(numWeights); + + outputDiagrams(outputDgm, outputCoordinates, Barycenters, dictDiagrams, + weightsVTK, vectorWeights, Spacing, 1); + + // Get input object from input vector + // Note: has to be a vtkDataSet as required by FillInputPortInformation + + // make a SHALLOW copy of the input + // outputDataSet->ShallowCopy(inputDataSet); + + // add to the output point data the computed output array + // outputDataSet->GetPointData()->AddArray(outputArray); + + // return success + return 1; +} + +void ttkPersistenceDiagramDictionaryDecoding::outputDiagrams( + vtkMultiBlockDataSet *output, + vtkTable *outputCoordinates, + const std::vector &diags, + std::vector &atoms, + vtkTable *weightsVTK, + const std::vector> &weights, + const double spacing, + const double maxPersistence) const { + + const auto nDiags = diags.size(); + const auto nAtoms = atoms.size(); + + ttk::SimplexId n_existing_blocks = ShowAtoms ? nAtoms : 0; + + output->SetNumberOfBlocks(nDiags + n_existing_blocks); + std::vector> coords(nAtoms); + std::vector> trueCoords(nAtoms); + std::vector xVector(nDiags); + std::vector yVector(nDiags); + std::vector zVector(nDiags, 0.); + vtkNew dummy{}; + + computeAtomsCoordinates(atoms, weights, coords, trueCoords, xVector, yVector, + zVector, spacing, nAtoms); + + if(nAtoms == 2) { + + if(ShowAtoms) { + for(size_t i = 0; i < nAtoms; ++i) { + double X = coords[i][0]; + double Y = coords[i][1]; + vtkNew vtu{}; + DiagramToVTU(vtu, atoms[i], dummy, *this, 3, false); + TranslateDiagram(vtu, std::array{X, Y, 0.0}); + output->SetBlock(i, vtu); + } + } + + } else if(nAtoms == 3) { + + if(ShowAtoms) { + for(size_t i = 0; i < nAtoms; ++i) { + double X = coords[i][0]; + double Y = coords[i][1]; + vtkNew vtu{}; + DiagramToVTU(vtu, atoms[i], dummy, *this, 3, false); + TranslateDiagram(vtu, std::array{X, Y, 0.0}); + output->SetBlock(i, vtu); + } + } + + } else { + + if(ShowAtoms) { + for(size_t i = 0; i < nAtoms; ++i) { + const auto angle + = 2.0 * M_PI * static_cast(i) / static_cast(nAtoms); + double X = spacing * maxPersistence * std::cos(angle); + double Y = spacing * maxPersistence * std::sin(angle); + vtkNew vtu{}; + DiagramToVTU(vtu, atoms[i], dummy, *this, 3, false); + TranslateDiagram(vtu, std::array{X, Y, 0.0}); + output->SetBlock(i, vtu); + } + } + } + + int numDiags = diags.size(); + + for(int i = 0; i < numDiags; ++i) { + vtkNew vtu{}; + if(!ComputePoints) { + DiagramToVTU(vtu, diags[i], dummy, *this, 3, false); + } + + double X = 0; + double Y = 0; + for(size_t iAtom = 0; iAtom < nAtoms; ++iAtom) { + X += weights[i][iAtom] * coords[iAtom][0]; + Y += weights[i][iAtom] * coords[iAtom][1]; + } + + TranslateDiagram(vtu, std::array{X, Y, 0.0}); + output->SetBlock(i + n_existing_blocks, vtu); + } + + for(size_t i = 0; i < 3; ++i) { + vtkNew col{}; + col->SetNumberOfValues(nDiags); + std::string name; + + if(i == 0) { + name = "X"; + } else if(i == 1) { + name = "Y"; + } else { + name = "Z"; + } + col->SetName(name.c_str()); + for(size_t j = 0; j < nDiags; ++j) { + if(i == 0) { + col->SetValue(j, xVector[j]); + } else if(i == 1) { + col->SetValue(j, yVector[j]); + } else { + col->SetValue(j, zVector[j]); + } + } + col->Modified(); + outputCoordinates->AddColumn(col); + } + + const auto zeroPad + = [](std::string &colName, const size_t numberCols, const size_t colIdx) { + std::string max{std::to_string(numberCols - 1)}; + std::string cur{std::to_string(colIdx)}; + std::string zer(max.size() - cur.size(), '0'); + colName.append(zer).append(cur); + }; + + vtkNew temp; + temp->DeepCopy(weightsVTK); + for(int i = 0; i < temp->GetNumberOfColumns(); ++i) { + int test = 0; + const auto array = temp->GetColumn(i); + + for(size_t j = 0; j < nDiags; ++j) { + std::string name{"Atom"}; + zeroPad(name, nDiags, j); + if(strcmp(name.c_str(), temp->GetColumnName(i)) == 0) { + test += 1; + } + } + if(test > 0) { + continue; + } + outputCoordinates->AddColumn(array); + } + + for(size_t i = 0; i < nAtoms; ++i) { + vtkNew row{}; + row->SetNumberOfValues(outputCoordinates->GetNumberOfColumns()); + for(int j = 0; j < outputCoordinates->GetNumberOfColumns(); ++j) { + if(strcmp(outputCoordinates->GetColumnName(j), "X") == 0) { + row->SetValue(j, trueCoords[i][0]); + } else if(strcmp(outputCoordinates->GetColumnName(j), "Y") == 0) { + row->SetValue(j, trueCoords[i][1]); + } else if(strcmp(outputCoordinates->GetColumnName(j), "Z") == 0) { + row->SetValue(j, trueCoords[i][2]); + } else if(strcmp(outputCoordinates->GetColumnName(j), "ClusterID") == 0) { + row->SetValue(j, -1); + } else { + continue; + } + } + row->Modified(); + outputCoordinates->InsertNextRow(row); + } +} + +double ttkPersistenceDiagramDictionaryDecoding::getMaxPersistence( + const ttk::DiagramType &diagram) const { + + double maxPersistence{0}; + for(size_t i = 0; i < diagram.size(); ++i) { + const auto &t = diagram[i]; + const double &pers = t.persistence(); + maxPersistence = std::max(pers, maxPersistence); + } + return maxPersistence; +} diff --git a/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.h b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.h new file mode 100644 index 0000000000..aedfea5d57 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/ttkPersistenceDiagramDictionaryDecoding.h @@ -0,0 +1,83 @@ +/// \ingroup vtk +/// \class ttkPersistenceDiagramDictionaryDecoding +/// \author Keanu Sisouk +/// \author Pierre Guillou +/// \date Mai 2023 +/// +/// \brief TTK processing package for the computation of a Dictionary +/// of Persistence Diagrams and barycentric weights to approximate +/// an ensemble of Persistence Diagrams. +/// +/// \b Related \b publication \n +/// "Wasserstein Dictionaries of Persistence Diagrams" \n +/// Keanu Sisouk, Julie Delon and Julien Tierny \n +/// IEEE Transactions on Visualization and Computer Graphics, 2023. +/// +/// \sa PersistenceDiagramDictionaryDecoding + +#pragma once + +#include "ttkMacros.h" +#include +#include +#include +#include +#include +// VTK Module +#include + +// VTK Includes +#include + +// TTK Base Includes +#include + +class TTKPERSISTENCEDIAGRAMDICTIONARYDECODING_EXPORT + ttkPersistenceDiagramDictionaryDecoding + : public ttkAlgorithm, + protected ttk::PersistenceDiagramDictionaryDecoding { + +private: +public: + static ttkPersistenceDiagramDictionaryDecoding *New(); + vtkTypeMacro(ttkPersistenceDiagramDictionaryDecoding, ttkAlgorithm); + + vtkGetMacro(Spacing, double); + vtkSetMacro(Spacing, double); + + vtkGetMacro(ShowAtoms, int); + vtkSetMacro(ShowAtoms, int); + + vtkGetMacro(ProgBarycenter, int); + vtkSetMacro(ProgBarycenter, int); + + ttkSetEnumMacro(ProjMet, BACKEND); + vtkGetEnumMacro(ProjMet, BACKEND); + +protected: + ttkPersistenceDiagramDictionaryDecoding(); + ~ttkPersistenceDiagramDictionaryDecoding() override = default; + + int FillInputPortInformation(int port, vtkInformation *info) override; + + int FillOutputPortInformation(int port, vtkInformation *info) override; + + void outputDiagrams(vtkMultiBlockDataSet *output, + vtkTable *output_coordinates, + const std::vector &diags, + std::vector &atoms, + vtkTable *weights_vtk, + const std::vector> &weights, + const double spacing, + const double maxPersistence) const; + + double getMaxPersistence(const ttk::DiagramType &diagram) const; + + int RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + + double Spacing{}; + int ShowAtoms{1}; + bool ComputePoints{false}; +}; diff --git a/core/vtk/ttkPersistenceDiagramDictionaryDecoding/vtk.module b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/vtk.module new file mode 100644 index 0000000000..ef10b07793 --- /dev/null +++ b/core/vtk/ttkPersistenceDiagramDictionaryDecoding/vtk.module @@ -0,0 +1,4 @@ +NAME + ttkPersistenceDiagramDictionaryDecoding +DEPENDS + ttkAlgorithm diff --git a/paraview/xmls/PersistenceDiagramDictionary.xml b/paraview/xmls/PersistenceDiagramDictionary.xml new file mode 100644 index 0000000000..a084a8ffc6 --- /dev/null +++ b/paraview/xmls/PersistenceDiagramDictionary.xml @@ -0,0 +1,303 @@ + + + + + This filter encodes new persistence diagrams by a dictionary and a vector of number. + + Related publication: "Progressive Wasserstein Barycenters of Persistence Diagrams" + Jules Vidal, Joseph Budin, Julien Tierny + IEEE Transactions on Visualization and Computer Graphics. + Proceedings of IEEE VIS 2019. + + See also PersistenceDiagram, BottleneckDistance + + + + + + + + + + + + + + + + Data-set to process. + + + + + + + + + + + + + + + + + Input Atoms + + + + + + + + + + + + + Backend for the computation of the dictionary of the persistence diagrams. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Performs the weight optimization + + + + + + + Performs the atoms optimization + + + + + + + + + Set strategy step to Max Eigen Value + + + + + + + Performs fusion + + + + + + + + Performs progressivity for computing barycenters + + + + + + + + + Activate progressive approach + + + + + + + + Activate the interrupt condition + + + + + + + Sort each diagrams by persistence for test + + + + + + + Performs the creation of features + + + + + + + Activate compression mode + + + + + + + Activate dimension reduction mode + + + + + + + + + + + + + + + + + + + + + + + + + ${DEBUG_WIDGETS} + + + + + + + + diff --git a/paraview/xmls/PersistenceDiagramDictionaryDecoding.xml b/paraview/xmls/PersistenceDiagramDictionaryDecoding.xml new file mode 100644 index 0000000000..ce6de7d2ee --- /dev/null +++ b/paraview/xmls/PersistenceDiagramDictionaryDecoding.xml @@ -0,0 +1,142 @@ + + + + + This filter encodes new persistence diagrams by a dictionary and a vector of number. + + Related publication: "Progressive Wasserstein Barycenters of Persistence Diagrams" + Jules Vidal, Joseph Budin, Julien Tierny + IEEE Transactions on Visualization and Computer Graphics. + Proceedings of IEEE VIS 2019. + + See also PersistenceDiagram, BottleneckDistance + + + + + + + + + + + + + + + Data-set to process. + + + + + + + + + + + + + + + + Data-set to process. + + + + + + + + + + + + + + + + + + + Backend for the computation of the 2D projection when the size of the dictionary is greater than 4. + + + + + + > + + Ajust the spacing for the display of diagrams. + + + + + > + + + Includes the atoms of the dictionnary in the output representation + + + + + + > + + + Compute the barycenters using progressivity + + + + + + + + + + + ${DEBUG_WIDGETS} + + + + + + + +