diff --git a/CMake/merge_tree_planar_layout.xml b/CMake/merge_tree_planar_layout.xml
index bf4b57cd3f..333b3035d1 100644
--- a/CMake/merge_tree_planar_layout.xml
+++ b/CMake/merge_tree_planar_layout.xml
@@ -17,6 +17,25 @@ panel_visibility="advanced">
+
+
+
+
+
+
+
+
+
+
+
diff --git a/CMake/merge_tree_preprocess.xml b/CMake/merge_tree_preprocess.xml
index 4b7c7bfc53..56e3cabc9a 100644
--- a/CMake/merge_tree_preprocess.xml
+++ b/CMake/merge_tree_preprocess.xml
@@ -20,6 +20,14 @@ panel_visibility="advanced">
mode="visibility"
property="Backend"
value="2" />
+
+
@@ -81,6 +89,14 @@ default_values="5">
mode="visibility"
property="Backend"
value="2" />
+
+
@@ -177,6 +193,14 @@ default_values="0">
mode="visibility"
property="Backend"
value="2" />
+
+
diff --git a/CMakePresets.json b/CMakePresets.json
index 732d7dd1aa..d720fb236f 100644
--- a/CMakePresets.json
+++ b/CMakePresets.json
@@ -1,41 +1,61 @@
{
- "version": 2,
- "configurePresets": [
- {
- "name": "TTK-Default",
- "hidden": true,
- "binaryDir": "build",
- "generator": "Ninja",
- "cacheVariables": {
- "TTK_ENABLE_DOUBLE_TEMPLATING": "ON"
- }
- },
- {
- "name": "TTK-Release",
- "inherits": "TTK-Default",
- "cacheVariables": {
- "CMAKE_BUILD_TYPE": "Release",
- "TTK_ENABLE_CPU_OPTIMIZATION": "OFF",
- "TTK_ENABLE_MPI": "OFF",
- "TTK_ENABLE_KAMIKAZE": "ON"
- }
- },
- {
- "name": "TTK-PerformanceBenchmark",
- "inherits": "TTK-Default",
- "cacheVariables": {
- "CMAKE_BUILD_TYPE": "Release",
- "TTK_ENABLE_CPU_OPTIMIZATION": "ON",
- "TTK_ENABLE_KAMIKAZE": "ON"
- }
- },
- {
- "name": "TTK-Debug",
- "inherits": "TTK-Default",
- "cacheVariables": {
- "CMAKE_BUILD_TYPE": "Debug",
- "TTK_ENABLE_KAMIKAZE": "OFF"
- }
- }
- ]
-}
+ "version": 6,
+ "configurePresets": [
+ {
+ "name": "TTK-Default",
+ "hidden": true,
+ "binaryDir": "build",
+ "generator": "Unix Makefiles",
+ "cacheVariables": {
+ "TTK_ENABLE_DOUBLE_TEMPLATING": "OFF"
+ }
+ },
+ {
+ "name": "TTK-Release",
+ "inherits": "TTK-Default",
+ "binaryDir": "build/release",
+ "cacheVariables": {
+ "CMAKE_BUILD_TYPE": "Release",
+ "TTK_ENABLE_CPU_OPTIMIZATION": "OFF",
+ "TTK_ENABLE_MPI": "OFF",
+ "TTK_ENABLE_KAMIKAZE": "ON",
+ "TTK_ENABLE_WEBSOCKETPP": "OFF",
+ "WEBSOCKETPP_DIR": "/home/wetzels/ttk/ttk-jonas/websocketpp/install/lib/cmake/websocketpp"
+ }
+ },
+ {
+ "name": "TTK-PerformanceBenchmark",
+ "inherits": "TTK-Default",
+ "cacheVariables": {
+ "CMAKE_BUILD_TYPE": "Release",
+ "TTK_ENABLE_CPU_OPTIMIZATION": "ON",
+ "TTK_ENABLE_KAMIKAZE": "ON"
+ }
+ },
+ {
+ "name": "TTK-Debug",
+ "inherits": "TTK-Default",
+ "binaryDir": "build/debug",
+ "cacheVariables": {
+ "CMAKE_BUILD_TYPE": "Debug",
+ "TTK_ENABLE_KAMIKAZE": "OFF"
+ }
+ }
+ ],
+ "buildPresets": [
+ {
+ "name": "TTK-Release",
+ "description": "",
+ "displayName": "",
+ "configurePreset": "TTK-Release",
+ "jobs": 4
+ },
+ {
+ "name": "TTK-Debug",
+ "description": "",
+ "displayName": "",
+ "configurePreset": "TTK-Debug",
+ "jobs": 4
+ }
+ ]
+}
\ No newline at end of file
diff --git a/core/base/ftmTree/FTMTreeUtils_Template.h b/core/base/ftmTree/FTMTreeUtils_Template.h
index 119befe4d0..587cb492f0 100644
--- a/core/base/ftmTree/FTMTreeUtils_Template.h
+++ b/core/base/ftmTree/FTMTreeUtils_Template.h
@@ -39,21 +39,63 @@ namespace ttk {
double threshold,
std::vector &excludeLower,
std::vector &excludeHigher) {
- dataType rootPers = this->getNodePersistence(this->getRoot());
- if(threshold > 1)
- threshold /= 100.0;
+ idNode treeRoot = this->getRoot();
+ dataType rootValue = this->getValue(treeRoot);
+ dataType lowestNodeValue
+ = this->getValue(this->getLowestNode(treeRoot));
+ dataType rootPers
+ = (rootValue > lowestNodeValue ? rootValue - lowestNodeValue
+ : lowestNodeValue - rootValue);
+ threshold /= 100.0;
threshold = rootPers * threshold;
- auto pers = this->getNodePersistence(nodeId);
-
- // Excluded pairs
- bool isExcluded = false;
- if(excludeLower.size() == excludeHigher.size())
- for(unsigned i = 0; i < excludeLower.size(); ++i) {
- isExcluded |= (pers > rootPers * excludeLower[i] / 100.0
- and pers < rootPers * excludeHigher[i] / 100.0);
- }
- return pers > threshold and not isExcluded;
+ auto isImportantPairOneNode = [&](idNode node) {
+ auto pers = this->getNodePersistence(node);
+
+ // Excluded pairs
+ bool isExcluded = false;
+ if(excludeLower.size() == excludeHigher.size())
+ for(unsigned i = 0; i < excludeLower.size(); ++i) {
+ isExcluded |= (pers > rootPers * excludeLower[i] / 100.0
+ and pers < rootPers * excludeHigher[i] / 100.0);
+ }
+ return (pers > threshold and not isExcluded);
+ };
+
+ if(isImportantPairOneNode(nodeId))
+ return true;
+
+ // Test if it is a parent of an important pair (for not persistence based
+ // branch decomposition)
+ idNode saddleNode
+ = (this->isLeaf(nodeId) ? this->getNode(nodeId)->getOrigin() : nodeId);
+ idNode leafNode
+ = (this->isLeaf(nodeId) ? nodeId : this->getNode(nodeId)->getOrigin());
+ std::queue queue;
+ std::vector nodeDone(this->getNumberOfNodes(), false);
+ queue.emplace(leafNode);
+ while(!queue.empty()) {
+ idNode node = queue.front();
+ queue.pop();
+
+ if(isImportantPairOneNode(node))
+ return true;
+
+ nodeDone[node] = true;
+
+ idNode parent = this->getParentSafe(node);
+ if(parent != saddleNode)
+ if(not nodeDone[parent])
+ queue.emplace(parent);
+
+ std::vector children;
+ this->getChildren(node, children);
+ for(auto child : children)
+ if(not nodeDone[child])
+ queue.emplace(child);
+ }
+
+ return false;
}
template
diff --git a/core/base/mergeTreeClustering/BranchMappingDistance.h b/core/base/mergeTreeClustering/BranchMappingDistance.h
index ceb9d1d932..de3e47c835 100644
--- a/core/base/mergeTreeClustering/BranchMappingDistance.h
+++ b/core/base/mergeTreeClustering/BranchMappingDistance.h
@@ -28,6 +28,7 @@
#include
// ttk common includes
+#include "MergeTreeBase.h"
#include
#include
#include
@@ -36,12 +37,17 @@
namespace ttk {
- class BranchMappingDistance : virtual public Debug {
+ class BranchMappingDistance : virtual public Debug, public MergeTreeBase {
private:
int baseMetric_ = 0;
int assignmentSolverID_ = 0;
bool squared_ = false;
+ bool computeMapping_ = false;
+ bool writeOptimalBranchDecomposition_ = false;
+
+ bool preprocess_ = true;
+ bool saveTree_ = false;
template
inline dataType editCost_Wasserstein1(int n1,
@@ -193,11 +199,66 @@ namespace ttk {
squared_ = s;
}
+ void setComputeMapping(bool m) {
+ computeMapping_ = m;
+ }
+
+ void setWriteBD(bool w) {
+ writeOptimalBranchDecomposition_ = w;
+ }
+
+ void setPreprocess(bool p) {
+ preprocess_ = p;
+ }
+
+ void setSaveTree(bool save) {
+ saveTree_ = save;
+ }
+
template
- dataType editDistance_branch(ftm::FTMTree_MT *tree1,
- ftm::FTMTree_MT *tree2) {
+ dataType execute(
+ ftm::MergeTree &mTree1,
+ ftm::MergeTree &mTree2,
+ std::vector> *outputMatching
+ = nullptr) {
- // initialize memoization tables
+ ftm::MergeTree mTree1Copy;
+ ftm::MergeTree mTree2Copy;
+ if(saveTree_) {
+ mTree1Copy = ftm::copyMergeTree(mTree1);
+ mTree2Copy = ftm::copyMergeTree(mTree2);
+ }
+ ftm::MergeTree &mTree1Int = (saveTree_ ? mTree1Copy : mTree1);
+ ftm::MergeTree &mTree2Int = (saveTree_ ? mTree2Copy : mTree2);
+ ftm::FTMTree_MT *tree1 = &(mTree1Int.tree);
+ ftm::FTMTree_MT *tree2 = &(mTree2Int.tree);
+
+ // optional preprocessing
+ if(preprocess_) {
+ treesNodeCorr_.resize(2);
+ preprocessingPipeline(
+ mTree1Int, epsilonTree1_, epsilon2Tree1_, epsilon3Tree1_, false,
+ useMinMaxPair_, cleanTree_, treesNodeCorr_[0], true, true);
+ preprocessingPipeline(
+ mTree2Int, epsilonTree2_, epsilon2Tree2_, epsilon3Tree2_, false,
+ useMinMaxPair_, cleanTree_, treesNodeCorr_[1], true, true);
+ }
+
+ tree1 = &(mTree1Int.tree);
+ tree2 = &(mTree2Int.tree);
+
+ return computeDistance(tree1, tree2, outputMatching);
+ }
+
+ template
+ dataType computeDistance(
+ ftm::FTMTree_MT *tree1,
+ ftm::FTMTree_MT *tree2,
+ std::vector> *outputMatching
+ = nullptr) {
+
+ // compute preorder of both trees (necessary for bottom-up dynamic
+ // programming)
std::vector> predecessors1(tree1->getNumberOfNodes());
std::vector> predecessors2(tree2->getNumberOfNodes());
@@ -248,6 +309,8 @@ namespace ttk {
}
}
+ // initialize memoization tables
+
size_t nn1 = tree1->getNumberOfNodes();
size_t nn2 = tree2->getNumberOfNodes();
size_t const dim1 = 1;
@@ -452,12 +515,14 @@ namespace ttk {
+ memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]);
} else {
for(auto child1_mb : children1) {
- auto topo1_ = children1;
+ std::vector topo1_;
+ tree1->getChildren(curr1, topo1_);
topo1_.erase(
std::remove(topo1_.begin(), topo1_.end(), child1_mb),
topo1_.end());
for(auto child2_mb : children2) {
- auto topo2_ = children2;
+ std::vector topo2_;
+ tree2->getChildren(curr2, topo2_);
topo2_.erase(
std::remove(topo2_.begin(), topo2_.end(), child2_mb),
topo2_.end());
@@ -551,8 +616,472 @@ namespace ttk {
dataType res
= memT[children1[0] + 1 * dim2 + children2[0] * dim3 + 1 * dim4];
+ if(computeMapping_ && outputMatching) {
+
+ outputMatching->clear();
+ std::vector matchedNodes(tree1->getNumberOfNodes(), -1);
+ std::vector matchedCost(tree1->getNumberOfNodes(), -1);
+ std::vector, std::pair>>
+ mapping;
+ std::vector linkedNodes1(tree1->getNumberOfNodes(), -1);
+ std::vector linkedNodes2(tree2->getNumberOfNodes(), -1);
+ traceMapping_branch(tree1, tree2, children1[0], 1, children2[0], 1,
+ predecessors1, predecessors2, depth1, depth2, memT,
+ mapping);
+ // dataType cost_mapping = 0;
+ for(auto m : mapping) {
+ if(writeOptimalBranchDecomposition_ && m.first.first >= 0
+ && m.first.second >= 0) {
+ tree1->getNode(m.first.first)->setOrigin(m.first.second);
+ tree1->getNode(m.first.second)->setOrigin(m.first.first);
+ linkedNodes1[m.first.first] = m.first.second;
+ linkedNodes1[m.first.second] = m.first.first;
+ }
+ if(writeOptimalBranchDecomposition_ && m.second.first >= 0
+ && m.second.second >= 0) {
+ tree2->getNode(m.second.first)->setOrigin(m.second.second);
+ tree2->getNode(m.second.second)->setOrigin(m.second.first);
+ linkedNodes2[m.second.first] = m.second.second;
+ linkedNodes2[m.second.second] = m.second.first;
+ }
+ if(m.first.first == -1)
+ continue;
+ if(m.first.second == -1)
+ continue;
+ if(m.second.first == -1)
+ continue;
+ if(m.second.second == -1)
+ continue;
+ matchedNodes[m.first.first] = m.second.first;
+ matchedNodes[m.first.second] = m.second.second;
+ matchedCost[m.first.first]
+ = this->baseMetric_ == 0 ? editCost_Wasserstein1(
+ m.first.first, m.first.second, m.second.first, m.second.second,
+ tree1, tree2)
+ : this->baseMetric_ == 1 ? editCost_Wasserstein2(
+ m.first.first, m.first.second, m.second.first,
+ m.second.second, tree1, tree2)
+ : this->baseMetric_ == 2
+ ? editCost_Persistence(m.first.first, m.first.second,
+ m.second.first,
+ m.second.second, tree1, tree2)
+ : editCost_Shifting(m.first.first, m.first.second,
+ m.second.first, m.second.second,
+ tree1, tree2);
+ matchedCost[m.first.second] = matchedCost[m.first.first];
+ }
+ for(ftm::idNode i = 0; i < matchedNodes.size(); i++) {
+ if(matchedNodes[i] >= 0)
+ outputMatching->emplace_back(
+ std::make_tuple(i, matchedNodes[i], matchedCost[i]));
+ }
+ }
+
return squared_ ? std::sqrt(res) : res;
}
+
+ template
+ void traceMapping_branch(
+ ftm::FTMTree_MT *tree1,
+ ftm::FTMTree_MT *tree2,
+ int curr1,
+ int l1,
+ int curr2,
+ int l2,
+ std::vector> &predecessors1,
+ std::vector> &predecessors2,
+ int depth1,
+ int depth2,
+ std::vector &memT,
+ std::vector, std::pair>>
+ &mapping) {
+
+ int nn1 = tree1->getNumberOfNodes();
+ int nn2 = tree2->getNumberOfNodes();
+ int dim1 = 1;
+ int dim2 = (nn1 + 1) * dim1;
+ int dim3 = (depth1 + 1) * dim2;
+ int dim4 = (nn2 + 1) * dim3;
+
+ //===============================================================================
+ // If second tree empty, track optimal branch decomposition of first tree
+
+ if(curr2 == nn2) {
+ std::vector children1;
+ tree1->getChildren(curr1, children1);
+ int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
+ //-----------------------------------------------------------------------
+ // If first subtree has only one branch, return deletion cost of this
+ // branch
+ if(tree1->getNumberOfChildren(curr1) == 0) {
+ mapping.emplace_back(std::make_pair(
+ std::make_pair(curr1, parent1), std::make_pair(-1, -1)));
+ return;
+ }
+ //-----------------------------------------------------------------------
+ // If first subtree has more than one branch, try all decompositions
+ else {
+ for(auto child1_mb : children1) {
+ dataType c_
+ = memT[child1_mb + (l1 + 1) * dim2 + nn2 * dim3 + 0 * dim4];
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ c_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
+ }
+ if(c_ == memT[curr1 + l1 * dim2 + nn2 * dim3 + 0 * dim4]) {
+ traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), nn2, 0,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
+ predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ }
+ return;
+ }
+ }
+ this->printErr("Mapping traceback not correct.");
+ }
+ }
+
+ //===============================================================================
+ // If first tree empty, track optimal branch decomposition of second tree
+
+ if(curr1 == nn1) {
+ std::vector children2;
+ tree2->getChildren(curr2, children2);
+ int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
+ //-----------------------------------------------------------------------
+ // If first subtree has only one branch, return deletion cost of this
+ // branch
+ if(tree2->getNumberOfChildren(curr2) == 0) {
+ mapping.emplace_back(std::make_pair(
+ std::make_pair(-1, -1), std::make_pair(curr2, parent2)));
+ return;
+ }
+ //-----------------------------------------------------------------------
+ // If first subtree has more than one branch, try all decompositions
+ else {
+ for(auto child2_mb : children2) {
+ dataType c_
+ = memT[nn1 + 0 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ c_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
+ }
+ if(c_ == memT[nn1 + 0 * dim2 + curr2 * dim3 + l2 * dim4]) {
+ traceMapping_branch(tree1, tree2, nn1, 0, child2_mb, (l2 + 1),
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
+ predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ }
+ return;
+ }
+ }
+ this->printErr("Mapping traceback not correct.");
+ }
+ }
+
+ std::vector children1;
+ tree1->getChildren(curr1, children1);
+ std::vector children2;
+ tree2->getChildren(curr2, children2);
+ int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
+ int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
+
+ //===============================================================================
+ // If both trees not empty, find optimal edit operation
+
+ //---------------------------------------------------------------------------
+ // If both trees only have one branch, return edit cost between
+ // the two branches
+ if(tree1->getNumberOfChildren(curr1) == 0
+ and tree2->getNumberOfChildren(curr2) == 0) {
+ mapping.emplace_back(std::make_pair(
+ std::make_pair(curr1, parent1), std::make_pair(curr2, parent2)));
+ return;
+ }
+ //---------------------------------------------------------------------------
+ // If first tree only has one branch, try all decompositions of
+ // second tree
+ else if(children1.size() == 0) {
+ for(auto child2_mb : children2) {
+ dataType d_
+ = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
+ }
+ if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
+ traceMapping_branch(tree1, tree2, curr1, l1, child2_mb, (l2 + 1),
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ }
+ return;
+ }
+ }
+ }
+ //---------------------------------------------------------------------------
+ // If second tree only has one branch, try all decompositions of
+ // first tree
+ else if(children2.size() == 0) {
+ for(auto child1_mb : children1) {
+ dataType d_
+ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
+ }
+ if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
+ traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), curr2, l2,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ }
+ return;
+ }
+ }
+ }
+ //---------------------------------------------------------------------------
+ // If both trees have more than one branch, try all decompositions
+ // of both trees
+ else {
+ //-----------------------------------------------------------------------
+ // Try all possible main branches of first tree (child1_mb) and
+ // all possible main branches of second tree (child2_mb) Then
+ // try all possible matchings of subtrees
+ if(children1.size() == 2 && children2.size() == 2) {
+ int child11 = children1[0];
+ int child12 = children1[1];
+ int child21 = children2[0];
+ int child22 = children2[1];
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
+ == memT[child11 + (l1 + 1) * dim2 + child21 * dim3
+ + (l2 + 1) * dim4]
+ + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
+
+ traceMapping_branch(tree1, tree2, child11, (l1 + 1), child21,
+ (l2 + 1), predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ traceMapping_branch(tree1, tree2, child12, 1, child22, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+
+ return;
+ }
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
+ == memT[child12 + (l1 + 1) * dim2 + child22 * dim3
+ + (l2 + 1) * dim4]
+ + memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
+
+ traceMapping_branch(tree1, tree2, child12, (l1 + 1), child22,
+ (l2 + 1), predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ traceMapping_branch(tree1, tree2, child11, 1, child21, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+
+ return;
+ }
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
+ == memT[child11 + (l1 + 1) * dim2 + child22 * dim3
+ + (l2 + 1) * dim4]
+ + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
+
+ traceMapping_branch(tree1, tree2, child11, (l1 + 1), child22,
+ (l2 + 1), predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ traceMapping_branch(tree1, tree2, child12, 1, child21, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+
+ return;
+ }
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
+ == memT[child12 + (l1 + 1) * dim2 + child21 * dim3
+ + (l2 + 1) * dim4]
+ + memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
+
+ traceMapping_branch(tree1, tree2, child12, (l1 + 1), child21,
+ (l2 + 1), predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ traceMapping_branch(tree1, tree2, child11, 1, child22, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+
+ return;
+ }
+ } else {
+ for(auto child1_mb : children1) {
+ std::vector topo1_;
+ tree1->getChildren(curr1, topo1_);
+ topo1_.erase(std::remove(topo1_.begin(), topo1_.end(), child1_mb),
+ topo1_.end());
+ for(auto child2_mb : children2) {
+ std::vector topo2_;
+ tree2->getChildren(curr2, topo2_);
+ topo2_.erase(std::remove(topo2_.begin(), topo2_.end(), child2_mb),
+ topo2_.end());
+
+ auto f = [&](unsigned r, unsigned c) {
+ int c1 = r < topo1_.size() ? topo1_[r] : -1;
+ int c2 = c < topo2_.size() ? topo2_[c] : -1;
+ return memT[c1 + 1 * dim2 + c2 * dim3 + 1 * dim4];
+ };
+ int size = std::max(topo1_.size(), topo2_.size()) + 1;
+ auto costMatrix = std::vector>(
+ size, std::vector(size, 0));
+ std::vector matching;
+ for(int r = 0; r < size; r++) {
+ for(int c = 0; c < size; c++) {
+ costMatrix[r][c] = f(r, c);
+ }
+ }
+
+ AssignmentSolver *assignmentSolver;
+ AssignmentExhaustive solverExhaustive;
+ AssignmentMunkres solverMunkres;
+ AssignmentAuction solverAuction;
+ switch(assignmentSolverID_) {
+ case 1:
+ solverExhaustive = AssignmentExhaustive();
+ assignmentSolver = &solverExhaustive;
+ break;
+ case 2:
+ solverMunkres = AssignmentMunkres();
+ assignmentSolver = &solverMunkres;
+ break;
+ case 0:
+ default:
+ solverAuction = AssignmentAuction();
+ assignmentSolver = &solverAuction;
+ }
+ assignmentSolver->setInput(costMatrix);
+ assignmentSolver->setBalanced(true);
+ assignmentSolver->run(matching);
+ dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + child2_mb * dim3
+ + (l2 + 1) * dim4];
+ for(auto m : matching)
+ d_ += std::get<2>(m);
+
+ if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
+ traceMapping_branch(
+ tree1, tree2, child1_mb, (l1 + 1), child2_mb, (l2 + 1),
+ predecessors1, predecessors2, depth1, depth2, memT, mapping);
+ for(auto m : matching) {
+ int n1 = std::get<0>(m) < static_cast(topo1_.size())
+ ? topo1_[std::get<0>(m)]
+ : -1;
+ int n2 = std::get<1>(m) < static_cast(topo2_.size())
+ ? topo2_[std::get<1>(m)]
+ : -1;
+ if(n1 >= 0 && n2 >= 0)
+ traceMapping_branch(tree1, tree2, n1, 1, n2, 1,
+ predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ else if(n1 >= 0)
+ traceMapping_branch(tree1, tree2, n1, 1, nn2, 0,
+ predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ else if(n2 >= 0)
+ traceMapping_branch(tree1, tree2, nn1, 0, n2, 1,
+ predecessors1, predecessors2, depth1,
+ depth2, memT, mapping);
+ }
+ return;
+ }
+ }
+ }
+ }
+ //-----------------------------------------------------------------------
+ // Try to continue main branch on one child of first tree and
+ // delete all other subtrees Then match continued branch to
+ // current branch in second tree
+ for(auto child1_mb : children1) {
+ dataType d_
+ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
+ }
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
+ traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), curr2, l2,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child1 : children1) {
+ if(child1 == child1_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ }
+ return;
+ }
+ }
+ //-----------------------------------------------------------------------
+ // Try to continue main branch on one child of second tree and
+ // delete all other subtrees Then match continued branch to
+ // current branch in first tree
+ for(auto child2_mb : children2) {
+ dataType d_
+ = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
+ }
+ if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
+ traceMapping_branch(tree1, tree2, curr1, l1, child2_mb, (l2 + 1),
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ for(auto child2 : children2) {
+ if(child2 == child2_mb) {
+ continue;
+ }
+ traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
+ predecessors1, predecessors2, depth1, depth2,
+ memT, mapping);
+ }
+ return;
+ }
+ }
+ this->printErr("Mapping traceback not correct");
+ }
+ }
};
} // namespace ttk
\ No newline at end of file
diff --git a/core/base/mergeTreeClustering/MergeTreeBarycenter.h b/core/base/mergeTreeClustering/MergeTreeBarycenter.h
index f74f1f6a28..b226aaff92 100644
--- a/core/base/mergeTreeClustering/MergeTreeBarycenter.h
+++ b/core/base/mergeTreeClustering/MergeTreeBarycenter.h
@@ -1,6 +1,7 @@
/// \ingroup base
/// \class MergeTreeBarycenter
/// \author Mathieu Pont (mathieu.pont@lip6.fr)
+/// \author Florian Wetzels (wetzels@cs.uni-kl.de)
/// \date 2021.
///
/// This module defines the %MergeTreeBarycenter class that computes
@@ -12,6 +13,12 @@
/// Proc. of IEEE VIS 2021.\n
/// IEEE Transactions on Visualization and Computer Graphics, 2021
+/// \b Related \b publication \n
+/// "Merge Tree Geodesics and Barycenters with Path Mappings" \n
+/// F. Wetzels, M. Pont, J. Tierny and C. Garth.\n
+/// Proc. of IEEE VIS 2023.\n
+/// IEEE Transactions on Visualization and Computer Graphics, 2024
+
#pragma once
#include
@@ -22,6 +29,10 @@
#include "MergeTreeBase.h"
#include "MergeTreeDistance.h"
+#include "PathMappingDistance.h"
+
+#include
+#include
namespace ttk {
@@ -49,6 +60,14 @@ namespace ttk {
bool preprocess_ = true;
bool postprocess_ = true;
+ int pathMetric_ = 0;
+ int baseModule_ = 0;
+ bool useMedianBarycenter_ = false;
+ bool useFixedInit_ = false;
+ // bool useEarlyOut_ = true;
+ int fixedInitNumber_ = 0;
+ int iterationLimit_ = 100;
+
// Output
std::vector finalDistances_;
@@ -119,6 +138,34 @@ namespace ttk {
return finalDistances_;
}
+ void setBaseModule(int m) {
+ baseModule_ = m;
+ }
+
+ void setPathMetric(int m) {
+ pathMetric_ = m;
+ }
+
+ void setUseMedianBarycenter(bool useMedian) {
+ useMedianBarycenter_ = useMedian;
+ }
+
+ void setUseFixedInit(bool useFixedInit) {
+ useFixedInit_ = useFixedInit;
+ }
+
+ // void setUseEarlyOut(bool useEarlyOut) {
+ // useEarlyOut_ = useEarlyOut;
+ // }
+
+ void setFixedInitNumber(int fixedInitNumber) {
+ fixedInitNumber_ = fixedInitNumber;
+ }
+
+ void setIterationLimit(int l) {
+ iterationLimit_ = l;
+ }
+
/**
* Implementation of the algorithm.
*/
@@ -140,9 +187,13 @@ namespace ttk {
for(unsigned int i = 0; i < trees.size(); ++i)
for(unsigned int j = i + 1; j < trees.size(); ++j) {
std::vector> matching;
+ std::vector,
+ std::pair>>
+ matching_path;
dataType distance;
- computeOneDistance(trees[i], trees2[j], matching, distance,
- useDoubleInput, isFirstInput);
+ computeOneDistance(trees[i], trees2[j], matching,
+ matching_path, distance, useDoubleInput,
+ isFirstInput);
distanceMatrix[i][j] = distance;
distanceMatrix[j][i] = distance;
}
@@ -273,9 +324,19 @@ namespace ttk {
void initBarycenterTree(std::vector &trees,
ftm::MergeTree &baryTree,
bool distMinimizer = true) {
- int const bestIndex
- = getBestInitTreeIndex(trees, distMinimizer);
- baryTree = ftm::copyMergeTree(trees[bestIndex], true);
+ int bestIndex;
+ if(useFixedInit_) {
+ if(fixedInitNumber_ >= 0 && fixedInitNumber_ < (int)trees.size())
+ bestIndex = fixedInitNumber_;
+ else
+ bestIndex = 0;
+ } else
+ bestIndex = getBestInitTreeIndex(trees, distMinimizer);
+ // bestIndex = 10;
+ // baryTree = ftm::copyMergeTree(trees[bestIndex], true);
+ baryTree
+ = ftm::copyMergeTree(trees[bestIndex], baseModule_ != 2);
+ // ftm::FTMTree_MT* bt = &(baryTree.tree);
limitSizeBarycenter(baryTree, trees);
}
@@ -687,6 +748,294 @@ namespace ttk {
ftm::cleanMergeTree(baryMergeTree);
}
+ template
+ void updateBarycenterTree_path(
+ std::vector &trees,
+ ftm::MergeTree &baryMergeTree,
+ std::vector &alphas,
+ std::vector,
+ std::pair>>>
+ &matchings) {
+ ftm::FTMTree_MT *baryTree = &(baryMergeTree.tree);
+ double alphaSum = 0;
+ for(unsigned int i = 0; i < trees.size(); ++i)
+ alphaSum += alphas[i];
+ bool joinTrees = trees[0]->isJoinTree();
+ int oldSize = baryTree->getNumberOfNodes();
+
+ // compute matched and unmatched nodes for all trees and barycenter
+ std::vector baryNodesMatched(baryTree->getNumberOfNodes(), false);
+ std::vector> treeNodesMatched(trees.size());
+ for(unsigned int i = 0; i < trees.size(); i++) {
+ if(alphas[i] == 0)
+ continue;
+ treeNodesMatched[i].resize(trees[i]->getNumberOfNodes(), false);
+ for(auto match : matchings[i]) {
+ baryNodesMatched[match.first.first] = true;
+ baryNodesMatched[match.first.second] = true;
+ treeNodesMatched[i][match.second.first] = true;
+ treeNodesMatched[i][match.second.second] = true;
+ }
+ }
+ // compute size of new barycenter tree
+ int newSize = oldSize;
+ for(unsigned int i = 0; i < treeNodesMatched.size(); i++) {
+ if(alphas[i] == 0 || useMedianBarycenter_)
+ continue;
+ for(unsigned int j = 0; j < treeNodesMatched[i].size(); j++) {
+ if(!treeNodesMatched[i][j])
+ newSize++;
+ }
+ }
+
+ // Create new barycenter tree
+ ftm::MergeTree baryMergeTreeNew
+ = ftm::createEmptyMergeTree(newSize);
+ // newScalars.resize(newSize);
+ // ftm::setTreeScalars(baryMergeTreeNew, newScalars);
+ ftm::FTMTree_MT *baryTreeNew = &(baryMergeTreeNew.tree);
+
+ // Copy the old tree structure
+ baryTreeNew->copyMergeTreeStructure(baryTree);
+
+ // delete not-matched nodes in barycenter
+ for(ftm::idNode i = 0; i < baryTree->getNumberOfNodes(); i++) {
+ if(not baryNodesMatched[i]) {
+ baryTreeNew->getNode(i)->setOrigin(-1);
+ baryTreeNew->deleteNode(i);
+ }
+ }
+
+ // relabel paths
+ std::vector> parentEdgeLengths(
+ baryTree->getNumberOfNodes());
+ for(unsigned int i = 0; i < trees.size(); i++) {
+ if(alphas[i] == 0)
+ continue;
+ auto tree = trees[i];
+ for(auto match : matchings[i]) {
+ dataType bv1 = baryTree->getValue(match.first.first);
+ dataType bv2 = baryTree->getValue(match.first.second);
+ dataType tv1 = tree->getValue(match.second.first);
+ dataType tv2 = tree->getValue(match.second.second);
+ dataType pathRangeB = bv1 > bv2 ? bv1 - bv2 : bv2 - bv1;
+ dataType pathRangeT = tv1 > tv2 ? tv1 - tv2 : tv2 - tv1;
+ ftm::idNode currB = baryTreeNew->getParentSafe(match.first.first);
+ ftm::idNode lastB = match.first.first;
+ while(lastB != match.first.second) {
+ dataType currValueB = baryTree->getValue(currB);
+ dataType lastValueB = baryTree->getValue(lastB);
+ dataType relativeValueB = lastValueB > currValueB
+ ? lastValueB - currValueB
+ : currValueB - lastValueB;
+ relativeValueB = relativeValueB / pathRangeB;
+ if(useMedianBarycenter_)
+ parentEdgeLengths[lastB].emplace_back(relativeValueB
+ * pathRangeT);
+ else
+ parentEdgeLengths[lastB].emplace_back(relativeValueB * pathRangeT
+ * alphas[i]);
+ // continue iteration
+ lastB = currB;
+ currB = baryTreeNew->getParentSafe(currB);
+ }
+ }
+ }
+ std::queue q;
+ q.push(baryTreeNew->getRoot());
+ // std::vector newScalars(baryTree->getNumberOfNodes(),0);
+ std::vector newScalars(newSize, 0);
+ newScalars[baryTreeNew->getRoot()]
+ = baryTree->getValue(baryTree->getRoot());
+ while(!q.empty()) {
+ auto curr = q.front();
+ q.pop();
+ std::vector children;
+ baryTreeNew->getChildren(curr, children);
+ for(auto child : children) {
+ q.emplace(child);
+ if(useMedianBarycenter_) {
+ auto m = parentEdgeLengths[child].begin()
+ + parentEdgeLengths[child].size() / 2;
+ std::nth_element(parentEdgeLengths[child].begin(), m,
+ parentEdgeLengths[child].end());
+ auto medianEdgeLength
+ = parentEdgeLengths[child][parentEdgeLengths[child].size() / 2];
+ newScalars[child]
+ = newScalars[curr]
+ + (joinTrees ? -medianEdgeLength : medianEdgeLength);
+ } else {
+ dataType avgEdgeLength = 0;
+ for(auto l : parentEdgeLengths[child]) {
+ avgEdgeLength += l;
+ }
+ // avgEdgeLength =
+ // avgEdgeLength/static_cast(trees.size());
+ avgEdgeLength = avgEdgeLength / alphaSum;
+ newScalars[child]
+ = newScalars[curr] + (joinTrees ? -avgEdgeLength : avgEdgeLength);
+ }
+ }
+ }
+ setTreeScalars(baryMergeTreeNew, newScalars);
+
+ // insert new nodes
+ int currSize = oldSize;
+ for(unsigned int i = 0; i < trees.size(); i++) {
+ if(alphas[i] == 0 || useMedianBarycenter_)
+ continue;
+ auto tree = trees[i];
+ std::vector newIndices(tree->getNumberOfNodes(), -1);
+ for(auto match : matchings[i]) {
+ dataType bv1 = baryTreeNew->getValue(match.first.first);
+ dataType bv2 = baryTreeNew->getValue(match.first.second);
+ dataType tv1 = tree->getValue(match.second.first);
+ dataType tv2 = tree->getValue(match.second.second);
+ dataType pathRangeB = bv1 > bv2 ? bv1 - bv2 : bv2 - bv1;
+ dataType pathRangeT = tv1 > tv2 ? tv1 - tv2 : tv2 - tv1;
+ ftm::idNode currB = baryTreeNew->getParentSafe(match.first.first);
+ ftm::idNode currT = tree->getParentSafe(match.second.first);
+ ftm::idNode lastB = match.first.first;
+ ftm::idNode lastT = match.second.first;
+ ftm::idNode lastNode = lastB;
+ while(currB != match.first.second || currT != match.second.second) {
+ dataType currValueB = baryTreeNew->getValue(currB);
+ dataType currValueT = tree->getValue(currT);
+ dataType relativeValueB
+ = bv1 > bv2 ? bv1 - currValueB : currValueB - bv1;
+ dataType relativeValueT
+ = tv1 > tv2 ? tv1 - currValueT : currValueT - tv1;
+ relativeValueB = relativeValueB / pathRangeB;
+ relativeValueT = relativeValueT / pathRangeT;
+ // if next node in barycenter, ignore
+ if(relativeValueB < relativeValueT) {
+ // continue iteration
+ lastB = currB;
+ currB = baryTreeNew->getParentSafe(currB);
+ lastNode = lastB;
+ }
+ // if next node in tree, add nodes
+ else if(relativeValueB > relativeValueT) {
+ q = std::queue();
+ std::vector currChildren;
+ tree->getChildren(currT, currChildren);
+ newIndices[currT] = currSize; // newScalars.size();
+ currSize++;
+ ftm::idNode nI = newIndices[currT];
+ // newScalars.emplace_back(tree->getValue(currT));
+ // newScalars.emplace_back(bv1 + (joinTrees ? relativeValueT *
+ // pathRangeB : - relativeValueT * pathRangeB));
+ newScalars[nI] = bv1
+ + (joinTrees ? relativeValueT * pathRangeB
+ : -relativeValueT * pathRangeB);
+ baryTreeNew->makeNode(nI);
+ baryTreeNew->setParent(nI, currB);
+ baryTreeNew->deleteParent(lastNode);
+ baryTreeNew->setParent(lastNode, nI);
+ baryTreeNew->getNode(nI)->setOrigin(-1);
+ std::vector nodesWithoutLink;
+ // baryTreeNew->getNode(nI)->setOrigin(newIndices[tree->getNode(currT)->getOrigin()]);
+ lastNode = newIndices[currT];
+ for(auto child : currChildren) {
+ if(child == lastT)
+ continue;
+ q.emplace(child);
+ newIndices[child] = currSize; // newScalars.size();
+ currSize++;
+ nI = newIndices[child];
+ // newScalars.emplace_back(tree->getValue(child));
+ dataType edgeLength
+ = (joinTrees ? tree->getValue(currT)
+ - tree->getValue(child)
+ : tree->getValue(child)
+ - tree->getValue(currT));
+ // newScalars.emplace_back(newScalars[newIndices[currT]] +
+ // (joinTrees ? - edgeLength * (alphas[i]/alphaSum) : edgeLength
+ // * (alphas[i]/alphaSum)));
+ newScalars[nI]
+ = newScalars[newIndices[currT]]
+ + (joinTrees ? -edgeLength * (alphas[i] / alphaSum)
+ : edgeLength * (alphas[i] / alphaSum));
+ baryTreeNew->makeNode(nI);
+ baryTreeNew->setParent(nI, newIndices[currT]);
+ baryTreeNew->getNode(nI)->setOrigin(-1);
+ if(tree->getNumberOfChildren(child) == 0
+ && newIndices[tree->getNode(child)->getOrigin()] >= 0) {
+ ftm::idNode ln
+ = newIndices[tree->getNode(child)->getOrigin()];
+ baryTreeNew->getNode(nI)->setOrigin(ln);
+ baryTreeNew->getNode(ln)->setOrigin(nI);
+ } else {
+ nodesWithoutLink.push_back(nI);
+ }
+ }
+ while(!q.empty()) {
+ auto currNode = q.front();
+ q.pop();
+ currChildren.clear();
+ tree->getChildren(currNode, currChildren);
+ for(auto child : currChildren) {
+ q.emplace(child);
+ newIndices[child] = currSize; // newScalars.size();
+ currSize++;
+ nI = newIndices[child];
+ // newScalars.emplace_back(tree->getValue(child));
+ dataType edgeLength
+ = (joinTrees ? tree->getValue(currNode)
+ - tree->getValue(child)
+ : tree->getValue(child)
+ - tree->getValue(currNode));
+ // newScalars.emplace_back(newScalars[newIndices[currNode]] +
+ // (joinTrees ? - edgeLength * (alphas[i]/alphaSum) :
+ // edgeLength * (alphas[i]/alphaSum)));
+ newScalars[nI]
+ = newScalars[newIndices[currNode]]
+ + (joinTrees ? -edgeLength * (alphas[i] / alphaSum)
+ : edgeLength * (alphas[i] / alphaSum));
+ baryTreeNew->makeNode(nI);
+ baryTreeNew->getNode(nI)->setOrigin(-1);
+ baryTreeNew->setParent(nI, newIndices[currNode]);
+ if(tree->getNumberOfChildren(child) == 0
+ && newIndices[tree->getNode(child)->getOrigin()] >= 0) {
+ ftm::idNode ln
+ = newIndices[tree->getNode(child)->getOrigin()];
+ baryTreeNew->getNode(nI)->setOrigin(ln);
+ baryTreeNew->getNode(ln)->setOrigin(nI);
+ } else {
+ nodesWithoutLink.push_back(nI);
+ }
+ }
+ }
+ // std::cout <<
+ // baryTreeNew->getNode(newIndices[currT])->getOrigin() << " " <<
+ // nodesWithoutLink.size() << std::endl;
+ if(baryTreeNew->getNode(newIndices[currT])->getOrigin() < 0) {
+ baryTreeNew->getNode(newIndices[currT])
+ ->setOrigin(nodesWithoutLink[0]);
+ }
+ for(ftm::idNode n : nodesWithoutLink) {
+ baryTreeNew->getNode(n)->setOrigin(newIndices[currT]);
+ }
+ // continue iteration
+ lastT = currT;
+ currT = tree->getParentSafe(currT);
+ } else {
+ // this should not happen
+ printErr("Impossible Matching behaviour.");
+ lastB = currB;
+ lastT = currT;
+ currB = baryTreeNew->getParentSafe(currB);
+ currT = tree->getParentSafe(currT);
+ }
+ }
+ }
+ setTreeScalars(baryMergeTreeNew, newScalars);
+ }
+
+ ftm::cleanMergeTree(baryMergeTreeNew, true);
+ baryMergeTree = baryMergeTreeNew;
+ }
+
template
void updateBarycenterTree(
std::vector &trees,
@@ -703,37 +1052,53 @@ namespace ttk {
// ------------------------------------------------------------------------
// Assignment
// ------------------------------------------------------------------------
+
template
void computeOneDistance(
ftm::FTMTree_MT *tree,
ftm::FTMTree_MT *baryTree,
std::vector> &matching,
+ std::vector,
+ std::pair>>
+ &matching_path,
dataType &distance,
bool useDoubleInput = false,
bool isFirstInput = true) {
// Timer t_distance;
- MergeTreeDistance mergeTreeDistance;
- mergeTreeDistance.setDebugLevel(std::min(debugLevel_, 2));
- mergeTreeDistance.setPreprocess(false);
- mergeTreeDistance.setPostprocess(false);
- mergeTreeDistance.setBranchDecomposition(true);
- mergeTreeDistance.setNormalizedWasserstein(normalizedWasserstein_);
- mergeTreeDistance.setKeepSubtree(keepSubtree_);
- mergeTreeDistance.setAssignmentSolver(assignmentSolverID_);
- mergeTreeDistance.setIsCalled(true);
- mergeTreeDistance.setThreadNumber(this->threadNumber_);
- mergeTreeDistance.setDistanceSquaredRoot(true); // squared root
- mergeTreeDistance.setNodePerTask(nodePerTask_);
- if(useDoubleInput) {
- double const weight = mixDistancesMinMaxPairWeight(isFirstInput);
- mergeTreeDistance.setMinMaxPairWeight(weight);
+ if(baseModule_ == 2) {
+ PathMappingDistance pathDistance;
+ pathDistance.setDebugLevel(std::min(debugLevel_, 2));
+ pathDistance.setPreprocess(false);
+ pathDistance.setAssignmentSolver(assignmentSolverID_);
+ pathDistance.setThreadNumber(this->threadNumber_);
+ pathDistance.setDistanceSquaredRoot(false); // squared root
+ pathDistance.setComputeMapping(true);
+ distance = pathDistance.computeDistance(
+ baryTree, tree, &matching, &matching_path);
+ } else {
+ MergeTreeDistance mergeTreeDistance;
+ mergeTreeDistance.setDebugLevel(std::min(debugLevel_, 2));
+ mergeTreeDistance.setPreprocess(false);
+ mergeTreeDistance.setPostprocess(false);
+ mergeTreeDistance.setBranchDecomposition(true);
+ mergeTreeDistance.setNormalizedWasserstein(normalizedWasserstein_);
+ mergeTreeDistance.setKeepSubtree(keepSubtree_);
+ mergeTreeDistance.setAssignmentSolver(assignmentSolverID_);
+ mergeTreeDistance.setIsCalled(true);
+ mergeTreeDistance.setThreadNumber(this->threadNumber_);
+ mergeTreeDistance.setDistanceSquaredRoot(true); // squared root
+ mergeTreeDistance.setNodePerTask(nodePerTask_);
+ if(useDoubleInput) {
+ double const weight = mixDistancesMinMaxPairWeight(isFirstInput);
+ mergeTreeDistance.setMinMaxPairWeight(weight);
+ }
+ /*if(progressiveBarycenter_){
+ mergeTreeDistance.setAuctionNoRounds(1);
+ mergeTreeDistance.setAuctionEpsilonDiviser(NoIteration-1);
+ }*/
+ distance = mergeTreeDistance.computeDistance(
+ baryTree, tree, matching);
}
- /*if(progressiveBarycenter_){
- mergeTreeDistance.setAuctionNoRounds(1);
- mergeTreeDistance.setAuctionEpsilonDiviser(NoIteration-1);
- }*/
- distance
- = mergeTreeDistance.computeDistance(baryTree, tree, matching);
std::stringstream ss, ss2;
ss << "distance tree : " << distance;
printMsg(ss.str(), debug::Priority::VERBOSE);
@@ -749,11 +1114,15 @@ namespace ttk {
ftm::FTMTree_MT *tree,
ftm::MergeTree &baryMergeTree,
std::vector> &matching,
+ std::vector,
+ std::pair>>
+ &matching_path,
dataType &distance,
bool useDoubleInput = false,
bool isFirstInput = true) {
computeOneDistance(tree, &(baryMergeTree.tree), matching,
- distance, useDoubleInput, isFirstInput);
+ matching_path, distance, useDoubleInput,
+ isFirstInput);
}
template
@@ -761,29 +1130,82 @@ namespace ttk {
ftm::MergeTree &baryMergeTree,
ftm::MergeTree &baryMergeTree2,
std::vector> &matching,
+ std::vector,
+ std::pair>>
+ &matching_path,
dataType &distance,
bool useDoubleInput = false,
bool isFirstInput = true) {
computeOneDistance(&(baryMergeTree.tree), baryMergeTree2,
- matching, distance, useDoubleInput,
+ matching, matching_path, distance,
+ useDoubleInput, isFirstInput);
+ }
+
+ template
+ void computeOneDistance(
+ ftm::FTMTree_MT *tree,
+ ftm::FTMTree_MT *baryTree,
+ std::vector> &matching,
+ dataType &distance,
+ bool useDoubleInput = false,
+ bool isFirstInput = true) {
+ std::vector,
+ std::pair>>
+ matching_path;
+ computeOneDistance(tree, baryTree, matching, matching_path,
+ distance, useDoubleInput, isFirstInput);
+ }
+
+ template
+ void computeOneDistance(
+ ftm::FTMTree_MT *tree,
+ ftm::MergeTree &baryMergeTree,
+ std::vector> &matching,
+ dataType &distance,
+ bool useDoubleInput = false,
+ bool isFirstInput = true) {
+ std::vector,
+ std::pair>>
+ matching_path;
+ computeOneDistance(tree, &(baryMergeTree.tree), matching,
+ matching_path, distance, useDoubleInput,
isFirstInput);
}
+ template
+ void computeOneDistance(
+ ftm::MergeTree &baryMergeTree,
+ ftm::MergeTree &baryMergeTree2,
+ std::vector> &matching,
+ dataType &distance,
+ bool useDoubleInput = false,
+ bool isFirstInput = true) {
+ std::vector,
+ std::pair>>
+ matching_path;
+ computeOneDistance(&(baryMergeTree.tree), baryMergeTree2,
+ matching, matching_path, distance,
+ useDoubleInput, isFirstInput);
+ }
+
template
void assignment(
std::vector &trees,
ftm::MergeTree &baryMergeTree,
std::vector>>
&matchings,
+ std::vector,
+ std::pair>>>
+ &matchings_path,
std::vector &distances,
bool useDoubleInput = false,
bool isFirstInput = true) {
if(not isCalled_)
- assignmentPara(trees, baryMergeTree, matchings, distances,
- useDoubleInput, isFirstInput);
+ assignmentPara(trees, baryMergeTree, matchings, matchings_path,
+ distances, useDoubleInput, isFirstInput);
else
- assignmentTask(trees, baryMergeTree, matchings, distances,
- useDoubleInput, isFirstInput);
+ assignmentTask(trees, baryMergeTree, matchings, matchings_path,
+ distances, useDoubleInput, isFirstInput);
}
template
@@ -792,6 +1214,9 @@ namespace ttk {
ftm::MergeTree &baryMergeTree,
std::vector>>
&matchings,
+ std::vector,
+ std::pair>>>
+ &matchings_path,
std::vector &distances,
bool useDoubleInput = false,
bool isFirstInput = true) {
@@ -801,8 +1226,8 @@ namespace ttk {
{
#pragma omp single nowait
#endif
- assignmentTask(trees, baryMergeTree, matchings, distances,
- useDoubleInput, isFirstInput);
+ assignmentTask(trees, baryMergeTree, matchings, matchings_path,
+ distances, useDoubleInput, isFirstInput);
#ifdef TTK_ENABLE_OPENMP4
} // pragma omp parallel
#endif
@@ -814,17 +1239,20 @@ namespace ttk {
ftm::MergeTree &baryMergeTree,
std::vector>>
&matchings,
+ std::vector,
+ std::pair>>>
+ &matchings_path,
std::vector &distances,
bool useDoubleInput = false,
bool isFirstInput = true) {
for(unsigned int i = 0; i < trees.size(); ++i)
#ifdef TTK_ENABLE_OPENMP4
#pragma omp task firstprivate(i) UNTIED() \
- shared(baryMergeTree, matchings, distances)
+ shared(baryMergeTree, matchings, matchings_path, distances)
#endif
computeOneDistance(trees[i], baryMergeTree, matchings[i],
- distances[i], useDoubleInput,
- isFirstInput);
+ matchings_path[i], distances[i],
+ useDoubleInput, isFirstInput);
#ifdef TTK_ENABLE_OPENMP4
#pragma omp taskwait
#endif
@@ -911,6 +1339,9 @@ namespace ttk {
std::vector &alphas,
std::vector>>
&finalMatchings,
+ std::vector,
+ std::pair>>>
+ &finalMatchings_path,
bool finalAsgnDoubleInput = false,
bool finalAsgnFirstInput = true) {
Timer t_bary;
@@ -939,7 +1370,11 @@ namespace ttk {
dataType minFrechet = std::numeric_limits::max();
int cptBlocked = 0;
int NoIteration = 0;
- while(not converged) {
+ std::stringstream energySequence;
+ int minBarySize = std::numeric_limits::max();
+ int maxBarySize = 0;
+ while(not converged
+ && (iterationLimit_ < 0 || NoIteration < iterationLimit_)) {
++NoIteration;
printMsg(debug::Separator::L2);
@@ -950,9 +1385,13 @@ namespace ttk {
// --- Assignment
std::vector>>
matchings(trees.size());
+ std::vector,
+ std::pair>>>
+ matchings_path(trees.size());
std::vector distances(trees.size(), -1);
Timer t_assignment;
- assignment(trees, baryMergeTree, matchings, distances);
+ assignment(
+ trees, baryMergeTree, matchings, matchings_path, distances);
Timer t_addDeletedNodes;
if(progressiveBarycenter_)
addScaledDeletedNodesCost(
@@ -965,7 +1404,13 @@ namespace ttk {
// --- Update
Timer t_update;
- updateBarycenterTree(trees, baryMergeTree, alphas, matchings);
+ if(baseModule_ == 2) {
+ updateBarycenterTree_path(
+ trees, baryMergeTree, alphas, matchings_path);
+ } else {
+ updateBarycenterTree(
+ trees, baryMergeTree, alphas, matchings);
+ }
auto t_update_time = t_update.getElapsedTime();
baryTree = &(baryMergeTree.tree);
printMsg("Update", 1, t_update_time, this->threadNumber_,
@@ -973,28 +1418,40 @@ namespace ttk {
// --- Check convergence
dataType currentFrechetEnergy = 0;
- for(unsigned int i = 0; i < trees.size(); ++i)
+ dataType currentFrechetEnergy2 = 0;
+ for(unsigned int i = 0; i < trees.size(); ++i) {
+ currentFrechetEnergy2 += alphas[i] * distances[i];
currentFrechetEnergy += alphas[i] * distances[i] * distances[i];
+ }
auto frechetDiff
= std::abs((double)(frechetEnergy - currentFrechetEnergy));
converged = (frechetDiff <= tol_);
converged = converged and (not progressiveBarycenter_ or treesUnscaled);
frechetEnergy = currentFrechetEnergy;
tol_ = frechetEnergy / 125.0;
+ energySequence << currentFrechetEnergy << std::endl;
- std::stringstream ss4;
+ std::stringstream ss4, ss5;
auto barycenterTime = t_bary.getElapsedTime() - addDeletedNodesTime_;
printMsg("Total", 1, barycenterTime, this->threadNumber_,
debug::LineMode::NEW, debug::Priority::INFO);
- printBaryStats(baryTree);
+ printBaryStats(baryTree, debug::Priority::INFO);
ss4 << "Frechet energy : " << frechetEnergy;
+ ss5 << "Frechet energy non-squared: " << currentFrechetEnergy2;
printMsg(ss4.str());
+ printMsg(ss5.str());
+
+ if((int)baryTree->getNumberOfNodes() > maxBarySize)
+ maxBarySize = baryTree->getNumberOfNodes();
+ if((int)baryTree->getNumberOfNodes() < minBarySize)
+ minBarySize = baryTree->getNumberOfNodes();
minFrechet = std::min(minFrechet, frechetEnergy);
if(not converged and (not progressiveBarycenter_ or treesUnscaled)) {
cptBlocked = (minFrechet < frechetEnergy) ? cptBlocked + 1 : 0;
converged = (cptBlocked >= 10);
}
+ // if(!useEarlyOut_) converged = false;
// --- Persistence scaling
if(progressiveBarycenter_) {
@@ -1004,28 +1461,54 @@ namespace ttk {
}
}
+ // std::ofstream energyFile;
+ // energyFile.open("/home/wetzels/ttk/energy.txt");
+ // energyFile << energySequence.str();
+ // energyFile.close();
+
// Final processing
printMsg(debug::Separator::L2);
printMsg("Final assignment");
std::vector distances(trees.size(), -1);
- assignment(trees, baryMergeTree, finalMatchings, distances,
- finalAsgnDoubleInput, finalAsgnFirstInput);
+ if(baseModule_ == 2) {
+ assignment(
+ trees, baryMergeTree, finalMatchings, finalMatchings_path, distances);
+ } else {
+ assignment(trees, baryMergeTree, finalMatchings,
+ finalMatchings_path, distances,
+ finalAsgnDoubleInput, finalAsgnFirstInput);
+ }
for(auto dist : distances)
finalDistances_.push_back(dist);
dataType currentFrechetEnergy = 0;
- for(unsigned int i = 0; i < trees.size(); ++i)
+ dataType currentFrechetEnergy2 = 0;
+ for(unsigned int i = 0; i < trees.size(); ++i) {
+ currentFrechetEnergy2 += alphas[i] * distances[i];
currentFrechetEnergy += alphas[i] * distances[i] * distances[i];
+ }
+ auto barycenterTime = t_bary.getElapsedTime() - addDeletedNodesTime_;
std::stringstream ss, ss2;
ss << "Frechet energy : " << currentFrechetEnergy;
+ ss2 << "Frechet energy non-squared: " << currentFrechetEnergy2;
printMsg(ss.str());
- auto barycenterTime = t_bary.getElapsedTime() - addDeletedNodesTime_;
+ printMsg(ss2.str());
printMsg("Total", 1, barycenterTime, this->threadNumber_,
- debug::LineMode::NEW, debug::Priority::INFO);
+ debug::LineMode::NEW, debug::Priority::PERFORMANCE);
// std::cout << "Bary Distance Time = " << allDistanceTime_ << std::endl;
- if(trees.size() == 2 and not isCalled_)
+ std::stringstream ssIt;
+ ssIt << "Number of iterations: " << NoIteration;
+ printMsg(ssIt.str(), debug::Priority::PERFORMANCE);
+ std::stringstream ssMin;
+ ssMin << "Min barycenter bize: " << minBarySize;
+ printMsg(ssMin.str(), debug::Priority::PERFORMANCE);
+ std::stringstream ssMax;
+ ssMax << "Max barycenter bize: " << maxBarySize;
+ printMsg(ssMax.str(), debug::Priority::PERFORMANCE);
+
+ if(trees.size() == 2 and not isCalled_ && baseModule_ != 2)
verifyBarycenterTwoTrees(
trees, baryMergeTree, finalMatchings, distances);
@@ -1043,17 +1526,21 @@ namespace ttk {
std::vector &alphas,
std::vector>>
&finalMatchings,
+ std::vector,
+ std::pair>>>
+ &finalMatchings_path,
ftm::MergeTree &baryMergeTree,
bool finalAsgnDoubleInput = false,
bool finalAsgnFirstInput = true) {
// --- Preprocessing
if(preprocess_) {
treesNodeCorr_.resize(trees.size());
- for(unsigned int i = 0; i < trees.size(); ++i)
- preprocessingPipeline(trees[i], epsilonTree2_,
- epsilon2Tree2_, epsilon3Tree2_,
- branchDecomposition_, useMinMaxPair_,
- cleanTree_, treesNodeCorr_[i]);
+ for(unsigned int i = 0; i < trees.size(); ++i) {
+ preprocessingPipeline(
+ trees[i], epsilonTree2_, epsilon2Tree2_, epsilon3Tree2_,
+ branchDecomposition_, useMinMaxPair_, cleanTree_, treesNodeCorr_[i],
+ true, baseModule_ == 2);
+ }
printTreesStats(trees);
}
@@ -1064,7 +1551,16 @@ namespace ttk {
// --- Execute
computeBarycenter(treesT, baryMergeTree, alphas, finalMatchings,
- finalAsgnDoubleInput, finalAsgnFirstInput);
+ finalMatchings_path, finalAsgnDoubleInput,
+ finalAsgnFirstInput);
+
+ if(baseModule_ == 2) {
+ ftm::FTMTree_MT *baryTree = &(baryMergeTree.tree);
+ for(ftm::idNode node = 0; node < baryTree->getNumberOfNodes(); node++) {
+ baryTree->getNode(node)->setOrigin(-1);
+ }
+ preprocessTree(baryTree, false);
+ }
// --- Postprocessing
if(postprocess_) {
@@ -1082,6 +1578,49 @@ namespace ttk {
}
}
+ template
+ void execute(
+ std::vector> &trees,
+ std::vector &alphas,
+ std::vector>>
+ &finalMatchings,
+ ftm::MergeTree &baryMergeTree,
+ bool finalAsgnDoubleInput = false,
+ bool finalAsgnFirstInput = true) {
+
+ std::vector,
+ std::pair>>>
+ finalMatchings_path;
+ execute(trees, alphas, finalMatchings, finalMatchings_path,
+ baryMergeTree, finalAsgnDoubleInput,
+ finalAsgnFirstInput);
+ }
+
+ template
+ void execute(
+ std::vector> &trees,
+ std::vector>>
+ &finalMatchings,
+ std::vector,
+ std::pair>>>
+ &finalMatchings_path,
+ ftm::MergeTree &baryMergeTree,
+ bool finalAsgnDoubleInput = false,
+ bool finalAsgnFirstInput = true) {
+ std::vector alphas;
+ if(trees.size() != 2) {
+ for(unsigned int i = 0; i < trees.size(); ++i)
+ alphas.push_back(1.0 / trees.size());
+ } else {
+ alphas.push_back(alpha_);
+ alphas.push_back(1 - alpha_);
+ }
+
+ execute(trees, alphas, finalMatchings, finalMatchings_path,
+ baryMergeTree, finalAsgnDoubleInput,
+ finalAsgnFirstInput);
+ }
+
template
void execute(
std::vector> &trees,
@@ -1099,8 +1638,13 @@ namespace ttk {
alphas.push_back(1 - alpha_);
}
- execute(trees, alphas, finalMatchings, baryMergeTree,
- finalAsgnDoubleInput, finalAsgnFirstInput);
+ std::vector,
+ std::pair>>>
+ finalMatchings_path;
+
+ execute(trees, alphas, finalMatchings, finalMatchings_path,
+ baryMergeTree, finalAsgnDoubleInput,
+ finalAsgnFirstInput);
}
// ------------------------------------------------------------------------
@@ -1204,8 +1748,11 @@ namespace ttk {
&finalMatchings,
std::vector distances) {
std::vector> matching;
+ std::vector,
+ std::pair>>
+ matching_path;
dataType distance;
- computeOneDistance(trees[0], trees[1], matching, distance);
+ computeOneDistance(trees[0], trees[1], matching, matching_path, distance);
if(distance != (distances[0] + distances[1])) {
std::stringstream ss, ss2, ss3, ss4;
ss << "distance T1 T2 : " << distance;
diff --git a/core/base/mergeTreeClustering/MergeTreeBase.h b/core/base/mergeTreeClustering/MergeTreeBase.h
index 64a1b0b49f..d13463307a 100644
--- a/core/base/mergeTreeClustering/MergeTreeBase.h
+++ b/core/base/mergeTreeClustering/MergeTreeBase.h
@@ -636,7 +636,8 @@ namespace ttk {
bool cleanTreeT,
double persistenceThreshold,
std::vector &nodeCorr,
- bool deleteInconsistentNodes = true) {
+ bool deleteInconsistentNodes = true,
+ bool removeMergedSaddles = false) {
Timer t_proc;
ftm::FTMTree_MT *tree = &(mTree.tree);
@@ -650,8 +651,18 @@ namespace ttk {
std::vector> treeNodeMerged(
tree->getNumberOfNodes());
if(not isPersistenceDiagram_ or convertToDiagram_) {
- if(epsilonTree != 0)
+ if(epsilonTree != 0) {
mergeSaddle(tree, epsilonTree, treeNodeMerged);
+ if(removeMergedSaddles) {
+ for(unsigned int j = 0; j < treeNodeMerged.size(); j++) {
+ for(auto k : treeNodeMerged[j]) {
+ auto nodeToDelete = tree->getNode(k)->getOrigin();
+ tree->getNode(k)->setOrigin(j);
+ tree->getNode(nodeToDelete)->setOrigin(-1);
+ }
+ }
+ }
+ }
}
// - Compute branch decomposition
@@ -701,11 +712,12 @@ namespace ttk {
bool useMinMaxPairT,
bool cleanTreeT,
std::vector