TTK
Loading...
Searching...
No Matches
BranchMappingDistance.h
Go to the documentation of this file.
1
13
14#pragma once
15
16#include <set>
17#include <vector>
18
19#include <algorithm>
20#include <cfloat>
21#include <chrono>
22#include <cmath>
23#include <iostream>
24#include <limits>
25#include <set>
26#include <stack>
27#include <tuple>
28#include <vector>
29
30// ttk common includes
31#include "MergeTreeBase.h"
32#include <AssignmentAuction.h>
34#include <AssignmentMunkres.h>
35#include <Debug.h>
36#include <FTMTree_MT.h>
37
38namespace ttk {
39
40 class BranchMappingDistance : virtual public Debug, public MergeTreeBase {
41
42 private:
43 int baseMetric_ = 0;
44 int assignmentSolverID_ = 0;
45 bool squared_ = false;
46 bool computeMapping_ = false;
47 bool writeOptimalBranchDecomposition_ = false;
48
49 bool preprocess_ = true;
50 bool saveTree_ = false;
51
52 template <class dataType>
53 inline dataType editCost_Wasserstein1(int n1,
54 int p1,
55 int n2,
56 int p2,
57 ftm::FTMTree_MT *tree1,
58 ftm::FTMTree_MT *tree2) {
59 dataType d;
60 if(n1 < 0) {
61 dataType b1 = tree2->getValue<dataType>(n2);
62 dataType d1 = tree2->getValue<dataType>(p2);
63 dataType b2 = (b1 + d1) * 0.5;
64 dataType d2 = (b1 + d1) * 0.5;
65 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
66 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
67 d = db + dd;
68 } else if(n2 < 0) {
69 dataType b1 = tree1->getValue<dataType>(n1);
70 dataType d1 = tree1->getValue<dataType>(p1);
71 dataType b2 = (b1 + d1) * 0.5;
72 dataType d2 = (b1 + d1) * 0.5;
73 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
74 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
75 d = db + dd;
76 } else {
77 dataType b1 = tree1->getValue<dataType>(n1);
78 dataType d1 = tree1->getValue<dataType>(p1);
79 dataType b2 = tree2->getValue<dataType>(n2);
80 dataType d2 = tree2->getValue<dataType>(p2);
81 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
82 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
83 d = db + dd;
84 }
85 return squared_ ? d * d : d;
86 }
87
88 template <class dataType>
89 inline dataType editCost_Wasserstein2(int n1,
90 int p1,
91 int n2,
92 int p2,
93 ftm::FTMTree_MT *tree1,
94 ftm::FTMTree_MT *tree2) {
95 dataType d;
96 if(n1 < 0) {
97 dataType b1 = tree2->getValue<dataType>(n2);
98 dataType d1 = tree2->getValue<dataType>(p2);
99 dataType b2 = (b1 + d1) * 0.5;
100 dataType d2 = (b1 + d1) * 0.5;
101 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
102 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
103 d = std::sqrt(db * db + dd * dd);
104 } else if(n2 < 0) {
105 dataType b1 = tree1->getValue<dataType>(n1);
106 dataType d1 = tree1->getValue<dataType>(p1);
107 dataType b2 = (b1 + d1) * 0.5;
108 dataType d2 = (b1 + d1) * 0.5;
109 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
110 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
111 d = std::sqrt(db * db + dd * dd);
112 } else {
113 dataType b1 = tree1->getValue<dataType>(n1);
114 dataType d1 = tree1->getValue<dataType>(p1);
115 dataType b2 = tree2->getValue<dataType>(n2);
116 dataType d2 = tree2->getValue<dataType>(p2);
117 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
118 dataType dd = d1 > d2 ? d1 - d2 : d2 - d1;
119 d = std::sqrt(db * db + dd * dd);
120 }
121 return squared_ ? d * d : d;
122 }
123
124 template <class dataType>
125 inline dataType editCost_Persistence(int n1,
126 int p1,
127 int n2,
128 int p2,
129 ftm::FTMTree_MT *tree1,
130 ftm::FTMTree_MT *tree2) {
131 dataType d;
132 if(n1 < 0) {
133 dataType b1 = tree2->getValue<dataType>(n2);
134 dataType d1 = tree2->getValue<dataType>(p2);
135 d = d1 > b1 ? d1 - b1 : b1 - d1;
136 } else if(n2 < 0) {
137 dataType b1 = tree1->getValue<dataType>(n1);
138 dataType d1 = tree1->getValue<dataType>(p1);
139 d = d1 > b1 ? d1 - b1 : b1 - d1;
140 } else {
141 dataType b1 = tree1->getValue<dataType>(n1);
142 dataType d1 = tree1->getValue<dataType>(p1);
143 dataType b2 = tree2->getValue<dataType>(n2);
144 dataType d2 = tree2->getValue<dataType>(p2);
145 dataType dist1 = d1 > b1 ? d1 - b1 : b1 - d1;
146 dataType dist2 = d2 > b2 ? d2 - b2 : b2 - d2;
147 d = dist1 > dist2 ? dist1 - dist2 : dist2 - dist1;
148 }
149 return squared_ ? d * d : d;
150 }
151
152 template <class dataType>
153 inline dataType editCost_Shifting(int n1,
154 int p1,
155 int n2,
156 int p2,
157 ftm::FTMTree_MT *tree1,
158 ftm::FTMTree_MT *tree2) {
159 dataType d;
160 if(n1 < 0) {
161 dataType b1 = tree2->getValue<dataType>(n2);
162 dataType d1 = tree2->getValue<dataType>(p2);
163 d = d1 > b1 ? d1 - b1 : b1 - d1;
164 } else if(n2 < 0) {
165 dataType b1 = tree1->getValue<dataType>(n1);
166 dataType d1 = tree1->getValue<dataType>(p1);
167 d = d1 > b1 ? d1 - b1 : b1 - d1;
168 } else {
169 dataType b1 = tree1->getValue<dataType>(n1);
170 dataType d1 = tree1->getValue<dataType>(p1);
171 dataType b2 = tree2->getValue<dataType>(n2);
172 dataType d2 = tree2->getValue<dataType>(p2);
173 dataType pers1 = d1 > b1 ? d1 - b1 : b1 - d1;
174 dataType pers2 = d2 > b2 ? d2 - b2 : b2 - d2;
175 dataType db = b1 > b2 ? b1 - b2 : b2 - b1;
176 dataType dp = pers1 > pers2 ? pers1 - pers2 : pers2 - pers1;
177 d = db + dp;
178 }
179 return squared_ ? d * d : d;
180 }
181
182 public:
184 this->setDebugMsgPrefix(
185 "MergeTreeDistance"); // inherited from Debug: prefix will be printed at
186 // the beginning of every msg
187 }
188 ~BranchMappingDistance() override = default;
189
190 void setBaseMetric(int m) {
191 baseMetric_ = m;
192 }
193
194 void setAssignmentSolver(int assignmentSolver) {
195 assignmentSolverID_ = assignmentSolver;
196 }
197
198 void setSquared(bool s) {
199 squared_ = s;
200 }
201
202 void setComputeMapping(bool m) {
203 computeMapping_ = m;
204 }
205
206 void setWriteBD(bool w) {
207 writeOptimalBranchDecomposition_ = w;
208 }
209
210 void setPreprocess(bool p) {
211 preprocess_ = p;
212 }
213
214 void setSaveTree(bool save) {
215 saveTree_ = save;
216 }
217
218 template <class dataType>
219 dataType execute(
222 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>> *outputMatching
223 = nullptr) {
224
225 ftm::MergeTree<dataType> mTree1Copy;
226 ftm::MergeTree<dataType> mTree2Copy;
227 if(saveTree_) {
228 mTree1Copy = ftm::copyMergeTree<dataType>(mTree1);
229 mTree2Copy = ftm::copyMergeTree<dataType>(mTree2);
230 }
231 ftm::MergeTree<dataType> &mTree1Int = (saveTree_ ? mTree1Copy : mTree1);
232 ftm::MergeTree<dataType> &mTree2Int = (saveTree_ ? mTree2Copy : mTree2);
233 ftm::FTMTree_MT *tree1 = &(mTree1Int.tree);
234 ftm::FTMTree_MT *tree2 = &(mTree2Int.tree);
235
236 // optional preprocessing
237 if(preprocess_) {
238 treesNodeCorr_.resize(2);
240 mTree1Int, epsilonTree1_, epsilon2Tree1_, epsilon3Tree1_, false,
241 useMinMaxPair_, cleanTree_, treesNodeCorr_[0], true, true);
243 mTree2Int, epsilonTree2_, epsilon2Tree2_, epsilon3Tree2_, false,
244 useMinMaxPair_, cleanTree_, treesNodeCorr_[1], true, true);
245 }
246
247 tree1 = &(mTree1Int.tree);
248 tree2 = &(mTree2Int.tree);
249
250 return computeDistance<dataType>(tree1, tree2, outputMatching);
251 }
252
253 template <class dataType>
255 ftm::FTMTree_MT *tree1,
256 ftm::FTMTree_MT *tree2,
257 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>> *outputMatching
258 = nullptr) {
259
260 // compute preorder of both trees (necessary for bottom-up dynamic
261 // programming)
262
263 std::vector<std::vector<int>> predecessors1(tree1->getNumberOfNodes());
264 std::vector<std::vector<int>> predecessors2(tree2->getNumberOfNodes());
265 int const rootID1 = tree1->getRoot();
266 int const rootID2 = tree2->getRoot();
267 std::vector<int> preorder1(tree1->getNumberOfNodes());
268 std::vector<int> preorder2(tree2->getNumberOfNodes());
269
270 int depth1 = 0;
271 int depth2 = 0;
272 std::stack<int> stack;
273 stack.push(rootID1);
274 int count = tree1->getNumberOfNodes() - 1;
275 while(!stack.empty()) {
276 int const nIdx = stack.top();
277 stack.pop();
278 preorder1[count] = nIdx;
279 count--;
280 depth1 = std::max((int)predecessors1[nIdx].size(), depth1);
281 std::vector<ftm::idNode> children;
282 tree1->getChildren(nIdx, children);
283 for(int const cIdx : children) {
284 stack.push(cIdx);
285 predecessors1[cIdx].reserve(predecessors1[nIdx].size() + 1);
286 predecessors1[cIdx].insert(predecessors1[cIdx].end(),
287 predecessors1[nIdx].begin(),
288 predecessors1[nIdx].end());
289 predecessors1[cIdx].push_back(nIdx);
290 }
291 }
292 stack.push(rootID2);
293 count = tree2->getNumberOfNodes() - 1;
294 while(!stack.empty()) {
295 int const nIdx = stack.top();
296 stack.pop();
297 preorder2[count] = nIdx;
298 count--;
299 depth2 = std::max((int)predecessors2[nIdx].size(), depth2);
300 std::vector<ftm::idNode> children;
301 tree2->getChildren(nIdx, children);
302 for(int const cIdx : children) {
303 stack.push(cIdx);
304 predecessors2[cIdx].reserve(predecessors2[nIdx].size() + 1);
305 predecessors2[cIdx].insert(predecessors2[cIdx].end(),
306 predecessors2[nIdx].begin(),
307 predecessors2[nIdx].end());
308 predecessors2[cIdx].push_back(nIdx);
309 }
310 }
311
312 // initialize memoization tables
313
314 size_t nn1 = tree1->getNumberOfNodes();
315 size_t nn2 = tree2->getNumberOfNodes();
316 size_t const dim1 = 1;
317 size_t const dim2 = (nn1 + 1) * dim1;
318 size_t const dim3 = (depth1 + 1) * dim2;
319 size_t const dim4 = (nn2 + 1) * dim3;
320
321 std::vector<dataType> memT((nn1 + 1) * (depth1 + 1) * (nn2 + 1)
322 * (depth2 + 1));
323
324 memT[nn1 + 0 * dim2 + nn2 * dim3 + 0 * dim4] = 0;
325 for(size_t i = 0; i < nn1; i++) {
326 int curr1 = preorder1[i];
327 std::vector<ftm::idNode> children1;
328 tree1->getChildren(curr1, children1);
329 for(size_t l = 1; l <= predecessors1[preorder1[i]].size(); l++) {
330 int parent1 = predecessors1[preorder1[i]]
331 [predecessors1[preorder1[i]].size() - l];
332
333 //-----------------------------------------------------------------------
334 // If first subtree has only one branch, return deletion cost of this
335 // branch
336 if(tree1->getNumberOfChildren(curr1) == 0) {
337 memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4]
338 = this->baseMetric_ == 0 ? editCost_Wasserstein1<dataType>(
339 curr1, parent1, -1, -1, tree1, tree2)
340 : this->baseMetric_ == 1 ? editCost_Wasserstein2<dataType>(
341 curr1, parent1, -1, -1, tree1, tree2)
342 : this->baseMetric_ == 2
343 ? editCost_Persistence<dataType>(
344 curr1, parent1, -1, -1, tree1, tree2)
345 : editCost_Shifting<dataType>(
346 curr1, parent1, -1, -1, tree1, tree2);
347 }
348 //-----------------------------------------------------------------------
349 // If first subtree has more than one branch, try all decompositions
350 else {
351 dataType c = std::numeric_limits<dataType>::max();
352 for(auto child1_mb : children1) {
353 dataType c_
354 = memT[child1_mb + (l + 1) * dim2 + nn2 * dim3 + 0 * dim4];
355 for(auto child1 : children1) {
356 if(child1 == child1_mb) {
357 continue;
358 }
359 c_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
360 }
361 c = std::min(c, c_);
362 }
363 memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4] = c;
364 }
365 }
366 }
367 for(size_t j = 0; j < nn2; j++) {
368 int curr2 = preorder2[j];
369 std::vector<ftm::idNode> children2;
370 tree2->getChildren(curr2, children2);
371 for(size_t l = 1; l <= predecessors2[preorder2[j]].size(); l++) {
372 int parent2 = predecessors2[preorder2[j]]
373 [predecessors2[preorder2[j]].size() - l];
374
375 //-----------------------------------------------------------------------
376 // If first subtree has only one branch, return deletion cost of this
377 // branch
378 if(tree2->getNumberOfChildren(curr2) == 0) {
379 memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4]
380 = this->baseMetric_ == 0 ? editCost_Wasserstein1<dataType>(
381 -1, -1, curr2, parent2, tree1, tree2)
382 : this->baseMetric_ == 1 ? editCost_Wasserstein2<dataType>(
383 -1, -1, curr2, parent2, tree1, tree2)
384 : this->baseMetric_ == 2
385 ? editCost_Persistence<dataType>(
386 -1, -1, curr2, parent2, tree1, tree2)
387 : editCost_Shifting<dataType>(
388 -1, -1, curr2, parent2, tree1, tree2);
389 }
390 //-----------------------------------------------------------------------
391 // If first subtree has more than one branch, try all decompositions
392 else {
393 dataType c = std::numeric_limits<dataType>::max();
394 for(auto child2_mb : children2) {
395 dataType c_
396 = memT[nn1 + 0 * dim2 + child2_mb * dim3 + (l + 1) * dim4];
397 for(auto child2 : children2) {
398 if(child2 == child2_mb) {
399 continue;
400 }
401 c_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
402 }
403 c = std::min(c, c_);
404 }
405 memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4] = c;
406 }
407 }
408 }
409
410 for(size_t i = 0; i < nn1; i++) {
411 int curr1 = preorder1[i];
412 std::vector<ftm::idNode> children1;
413 tree1->getChildren(curr1, children1);
414 for(size_t j = 0; j < nn2; j++) {
415 int curr2 = preorder2[j];
416 std::vector<ftm::idNode> children2;
417 tree2->getChildren(curr2, children2);
418 for(size_t l1 = 1; l1 <= predecessors1[preorder1[i]].size(); l1++) {
419 int parent1
420 = predecessors1[preorder1[i]]
421 [predecessors1[preorder1[i]].size() - l1];
422 for(size_t l2 = 1; l2 <= predecessors2[preorder2[j]].size(); l2++) {
423 int parent2
424 = predecessors2[preorder2[j]]
425 [predecessors2[preorder2[j]].size() - l2];
426
427 //===============================================================================
428 // If both trees not empty, find optimal edit operation
429
430 //---------------------------------------------------------------------------
431 // If both trees only have one branch, return edit cost between
432 // the two branches
433 if(tree1->getNumberOfChildren(curr1) == 0
434 and tree2->getNumberOfChildren(curr2) == 0) {
435 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
436 = this->baseMetric_ == 0 ? editCost_Wasserstein1<dataType>(
437 curr1, parent1, curr2, parent2, tree1, tree2)
438 : this->baseMetric_ == 1 ? editCost_Wasserstein2<dataType>(
439 curr1, parent1, curr2, parent2, tree1, tree2)
440 : this->baseMetric_ == 2
441 ? editCost_Persistence<dataType>(
442 curr1, parent1, curr2, parent2, tree1, tree2)
443 : editCost_Shifting<dataType>(
444 curr1, parent1, curr2, parent2, tree1, tree2);
445 }
446 //---------------------------------------------------------------------------
447 // If first tree only has one branch, try all decompositions of
448 // second tree
449 else if(children1.size() == 0) {
450 dataType d = std::numeric_limits<dataType>::max();
451 for(auto child2_mb : children2) {
452 dataType d_ = memT[curr1 + l1 * dim2 + child2_mb * dim3
453 + (l2 + 1) * dim4];
454 for(auto child2 : children2) {
455 if(child2 == child2_mb) {
456 continue;
457 }
458 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
459 }
460 d = std::min(d, d_);
461 }
462 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
463 }
464 //---------------------------------------------------------------------------
465 // If second tree only has one branch, try all decompositions of
466 // first tree
467 else if(children2.size() == 0) {
468 dataType d = std::numeric_limits<dataType>::max();
469 for(auto child1_mb : children1) {
470 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3
471 + l2 * dim4];
472 for(auto child1 : children1) {
473 if(child1 == child1_mb) {
474 continue;
475 }
476 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
477 }
478 d = std::min(d, d_);
479 }
480 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
481 }
482 //---------------------------------------------------------------------------
483 // If both trees have more than one branch, try all decompositions
484 // of both trees
485 else {
486 dataType d = std::numeric_limits<dataType>::max();
487 //-----------------------------------------------------------------------
488 // Try all possible main branches of first tree (child1_mb) and
489 // all possible main branches of second tree (child2_mb) Then
490 // try all possible matchings of subtrees
491 if(children1.size() == 2 && children2.size() == 2) {
492 int const child11 = children1[0];
493 int const child12 = children1[1];
494 int const child21 = children2[0];
495 int const child22 = children2[1];
496 d = std::min<dataType>(
497 d,
498 memT[child11 + (l1 + 1) * dim2 + child21 * dim3
499 + (l2 + 1) * dim4]
500 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]);
501 d = std::min<dataType>(
502 d,
503 memT[child12 + (l1 + 1) * dim2 + child22 * dim3
504 + (l2 + 1) * dim4]
505 + memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]);
506 d = std::min<dataType>(
507 d,
508 memT[child11 + (l1 + 1) * dim2 + child22 * dim3
509 + (l2 + 1) * dim4]
510 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]);
511 d = std::min<dataType>(
512 d,
513 memT[child12 + (l1 + 1) * dim2 + child21 * dim3
514 + (l2 + 1) * dim4]
515 + memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]);
516 } else {
517 for(auto child1_mb : children1) {
518 std::vector<ftm::idNode> topo1_;
519 tree1->getChildren(curr1, topo1_);
520 topo1_.erase(
521 std::remove(topo1_.begin(), topo1_.end(), child1_mb),
522 topo1_.end());
523 for(auto child2_mb : children2) {
524 std::vector<ftm::idNode> topo2_;
525 tree2->getChildren(curr2, topo2_);
526 topo2_.erase(
527 std::remove(topo2_.begin(), topo2_.end(), child2_mb),
528 topo2_.end());
529
530 auto f = [&](unsigned r, unsigned c) {
531 int const c1 = r < topo1_.size() ? topo1_[r] : -1;
532 int const c2 = c < topo2_.size() ? topo2_[c] : -1;
533 return memT[c1 + 1 * dim2 + c2 * dim3 + 1 * dim4];
534 };
535 int size = std::max(topo1_.size(), topo2_.size()) + 1;
536 auto costMatrix = std::vector<std::vector<dataType>>(
537 size, std::vector<dataType>(size, 0));
538 std::vector<MatchingType> matching;
539 for(int r = 0; r < size; r++) {
540 for(int c = 0; c < size; c++) {
541 costMatrix[r][c] = f(r, c);
542 }
543 }
544
545 AssignmentSolver<dataType> *assignmentSolver;
546 AssignmentExhaustive<dataType> solverExhaustive;
547 AssignmentMunkres<dataType> solverMunkres;
548 AssignmentAuction<dataType> solverAuction;
549 switch(assignmentSolverID_) {
550 case 1:
551 solverExhaustive = AssignmentExhaustive<dataType>();
552 assignmentSolver = &solverExhaustive;
553 break;
554 case 2:
555 solverMunkres = AssignmentMunkres<dataType>();
556 assignmentSolver = &solverMunkres;
557 break;
558 case 0:
559 default:
560 solverAuction = AssignmentAuction<dataType>();
561 assignmentSolver = &solverAuction;
562 }
563 assignmentSolver->setInput(costMatrix);
564 assignmentSolver->setBalanced(true);
565 assignmentSolver->run(matching);
566 dataType d_ = memT[child1_mb + (l1 + 1) * dim2
567 + child2_mb * dim3 + (l2 + 1) * dim4];
568 for(auto m : matching)
569 d_ += std::get<2>(m);
570 d = std::min(d, d_);
571 }
572 }
573 }
574 //-----------------------------------------------------------------------
575 // Try to continue main branch on one child of first tree and
576 // delete all other subtrees Then match continued branch to
577 // current branch in second tree
578 for(auto child1_mb : children1) {
579 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3
580 + l2 * dim4];
581 for(auto child1 : children1) {
582 if(child1 == child1_mb) {
583 continue;
584 }
585 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
586 }
587 d = std::min(d, d_);
588 }
589 //-----------------------------------------------------------------------
590 // Try to continue main branch on one child of second tree and
591 // delete all other subtrees Then match continued branch to
592 // current branch in first tree
593 for(auto child2_mb : children2) {
594 dataType d_ = memT[curr1 + l1 * dim2 + child2_mb * dim3
595 + (l2 + 1) * dim4];
596 for(auto child2 : children2) {
597 if(child2 == child2_mb) {
598 continue;
599 }
600 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
601 }
602 d = std::min(d, d_);
603 }
604 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
605 }
606 }
607 }
608 }
609 }
610
611 std::vector<ftm::idNode> children1;
612 tree1->getChildren(rootID1, children1);
613 std::vector<ftm::idNode> children2;
614 tree2->getChildren(rootID2, children2);
615
616 dataType res
617 = memT[children1[0] + 1 * dim2 + children2[0] * dim3 + 1 * dim4];
618
619 if(computeMapping_ && outputMatching) {
620
621 outputMatching->clear();
622 std::vector<int> matchedNodes(tree1->getNumberOfNodes(), -1);
623 std::vector<dataType> matchedCost(tree1->getNumberOfNodes(), -1);
624 std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>>
625 mapping;
626 std::vector<int> linkedNodes1(tree1->getNumberOfNodes(), -1);
627 std::vector<int> linkedNodes2(tree2->getNumberOfNodes(), -1);
628 traceMapping_branch(tree1, tree2, children1[0], 1, children2[0], 1,
629 predecessors1, predecessors2, depth1, depth2, memT,
630 mapping);
631 // dataType cost_mapping = 0;
632 for(auto m : mapping) {
633 if(writeOptimalBranchDecomposition_ && m.first.first >= 0
634 && m.first.second >= 0) {
635 tree1->getNode(m.first.first)->setOrigin(m.first.second);
636 tree1->getNode(m.first.second)->setOrigin(m.first.first);
637 linkedNodes1[m.first.first] = m.first.second;
638 linkedNodes1[m.first.second] = m.first.first;
639 }
640 if(writeOptimalBranchDecomposition_ && m.second.first >= 0
641 && m.second.second >= 0) {
642 tree2->getNode(m.second.first)->setOrigin(m.second.second);
643 tree2->getNode(m.second.second)->setOrigin(m.second.first);
644 linkedNodes2[m.second.first] = m.second.second;
645 linkedNodes2[m.second.second] = m.second.first;
646 }
647 if(m.first.first == -1)
648 continue;
649 if(m.first.second == -1)
650 continue;
651 if(m.second.first == -1)
652 continue;
653 if(m.second.second == -1)
654 continue;
655 matchedNodes[m.first.first] = m.second.first;
656 matchedNodes[m.first.second] = m.second.second;
657 matchedCost[m.first.first]
658 = this->baseMetric_ == 0 ? editCost_Wasserstein1<dataType>(
659 m.first.first, m.first.second, m.second.first, m.second.second,
660 tree1, tree2)
661 : this->baseMetric_ == 1 ? editCost_Wasserstein2<dataType>(
662 m.first.first, m.first.second, m.second.first,
663 m.second.second, tree1, tree2)
664 : this->baseMetric_ == 2
665 ? editCost_Persistence<dataType>(m.first.first, m.first.second,
666 m.second.first,
667 m.second.second, tree1, tree2)
668 : editCost_Shifting<dataType>(m.first.first, m.first.second,
669 m.second.first, m.second.second,
670 tree1, tree2);
671 matchedCost[m.first.second] = matchedCost[m.first.first];
672 }
673 for(ftm::idNode i = 0; i < matchedNodes.size(); i++) {
674 if(matchedNodes[i] >= 0)
675 outputMatching->emplace_back(
676 std::make_tuple(i, matchedNodes[i], matchedCost[i]));
677 }
678 }
679
680 return squared_ ? std::sqrt(res) : res;
681 }
682
683 template <class dataType>
685 ftm::FTMTree_MT *tree1,
686 ftm::FTMTree_MT *tree2,
687 int curr1,
688 int l1,
689 int curr2,
690 int l2,
691 std::vector<std::vector<int>> &predecessors1,
692 std::vector<std::vector<int>> &predecessors2,
693 int depth1,
694 int depth2,
695 std::vector<dataType> &memT,
696 std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>>
697 &mapping) {
698
699 int nn1 = tree1->getNumberOfNodes();
700 int nn2 = tree2->getNumberOfNodes();
701 int dim1 = 1;
702 int dim2 = (nn1 + 1) * dim1;
703 int dim3 = (depth1 + 1) * dim2;
704 int dim4 = (nn2 + 1) * dim3;
705
706 //===============================================================================
707 // If second tree empty, track optimal branch decomposition of first tree
708
709 if(curr2 == nn2) {
710 std::vector<ftm::idNode> children1;
711 tree1->getChildren(curr1, children1);
712 int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
713 //-----------------------------------------------------------------------
714 // If first subtree has only one branch, return deletion cost of this
715 // branch
716 if(tree1->getNumberOfChildren(curr1) == 0) {
717 mapping.emplace_back(std::make_pair(
718 std::make_pair(curr1, parent1), std::make_pair(-1, -1)));
719 return;
720 }
721 //-----------------------------------------------------------------------
722 // If first subtree has more than one branch, try all decompositions
723 else {
724 for(auto child1_mb : children1) {
725 dataType c_
726 = memT[child1_mb + (l1 + 1) * dim2 + nn2 * dim3 + 0 * dim4];
727 for(auto child1 : children1) {
728 if(child1 == child1_mb) {
729 continue;
730 }
731 c_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
732 }
733 if(c_ == memT[curr1 + l1 * dim2 + nn2 * dim3 + 0 * dim4]) {
734 traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), nn2, 0,
735 predecessors1, predecessors2, depth1, depth2,
736 memT, mapping);
737 for(auto child1 : children1) {
738 if(child1 == child1_mb) {
739 continue;
740 }
741 traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
742 predecessors1, predecessors2, depth1,
743 depth2, memT, mapping);
744 }
745 return;
746 }
747 }
748 this->printErr("Mapping traceback not correct.");
749 }
750 }
751
752 //===============================================================================
753 // If first tree empty, track optimal branch decomposition of second tree
754
755 if(curr1 == nn1) {
756 std::vector<ftm::idNode> children2;
757 tree2->getChildren(curr2, children2);
758 int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
759 //-----------------------------------------------------------------------
760 // If first subtree has only one branch, return deletion cost of this
761 // branch
762 if(tree2->getNumberOfChildren(curr2) == 0) {
763 mapping.emplace_back(std::make_pair(
764 std::make_pair(-1, -1), std::make_pair(curr2, parent2)));
765 return;
766 }
767 //-----------------------------------------------------------------------
768 // If first subtree has more than one branch, try all decompositions
769 else {
770 for(auto child2_mb : children2) {
771 dataType c_
772 = memT[nn1 + 0 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
773 for(auto child2 : children2) {
774 if(child2 == child2_mb) {
775 continue;
776 }
777 c_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
778 }
779 if(c_ == memT[nn1 + 0 * dim2 + curr2 * dim3 + l2 * dim4]) {
780 traceMapping_branch(tree1, tree2, nn1, 0, child2_mb, (l2 + 1),
781 predecessors1, predecessors2, depth1, depth2,
782 memT, mapping);
783 for(auto child2 : children2) {
784 if(child2 == child2_mb) {
785 continue;
786 }
787 traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
788 predecessors1, predecessors2, depth1,
789 depth2, memT, mapping);
790 }
791 return;
792 }
793 }
794 this->printErr("Mapping traceback not correct.");
795 }
796 }
797
798 std::vector<ftm::idNode> children1;
799 tree1->getChildren(curr1, children1);
800 std::vector<ftm::idNode> children2;
801 tree2->getChildren(curr2, children2);
802 int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
803 int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
804
805 //===============================================================================
806 // If both trees not empty, find optimal edit operation
807
808 //---------------------------------------------------------------------------
809 // If both trees only have one branch, return edit cost between
810 // the two branches
811 if(tree1->getNumberOfChildren(curr1) == 0
812 and tree2->getNumberOfChildren(curr2) == 0) {
813 mapping.emplace_back(std::make_pair(
814 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2)));
815 return;
816 }
817 //---------------------------------------------------------------------------
818 // If first tree only has one branch, try all decompositions of
819 // second tree
820 else if(children1.size() == 0) {
821 for(auto child2_mb : children2) {
822 dataType d_
823 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
824 for(auto child2 : children2) {
825 if(child2 == child2_mb) {
826 continue;
827 }
828 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
829 }
830 if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
831 traceMapping_branch(tree1, tree2, curr1, l1, child2_mb, (l2 + 1),
832 predecessors1, predecessors2, depth1, depth2,
833 memT, mapping);
834 for(auto child2 : children2) {
835 if(child2 == child2_mb) {
836 continue;
837 }
838 traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
839 predecessors1, predecessors2, depth1, depth2,
840 memT, mapping);
841 }
842 return;
843 }
844 }
845 }
846 //---------------------------------------------------------------------------
847 // If second tree only has one branch, try all decompositions of
848 // first tree
849 else if(children2.size() == 0) {
850 for(auto child1_mb : children1) {
851 dataType d_
852 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
853 for(auto child1 : children1) {
854 if(child1 == child1_mb) {
855 continue;
856 }
857 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
858 }
859 if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
860 traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), curr2, l2,
861 predecessors1, predecessors2, depth1, depth2,
862 memT, mapping);
863 for(auto child1 : children1) {
864 if(child1 == child1_mb) {
865 continue;
866 }
867 traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
868 predecessors1, predecessors2, depth1, depth2,
869 memT, mapping);
870 }
871 return;
872 }
873 }
874 }
875 //---------------------------------------------------------------------------
876 // If both trees have more than one branch, try all decompositions
877 // of both trees
878 else {
879 //-----------------------------------------------------------------------
880 // Try all possible main branches of first tree (child1_mb) and
881 // all possible main branches of second tree (child2_mb) Then
882 // try all possible matchings of subtrees
883 if(children1.size() == 2 && children2.size() == 2) {
884 int child11 = children1[0];
885 int child12 = children1[1];
886 int child21 = children2[0];
887 int child22 = children2[1];
888 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
889 == memT[child11 + (l1 + 1) * dim2 + child21 * dim3
890 + (l2 + 1) * dim4]
891 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
892
893 traceMapping_branch(tree1, tree2, child11, (l1 + 1), child21,
894 (l2 + 1), predecessors1, predecessors2, depth1,
895 depth2, memT, mapping);
896 traceMapping_branch(tree1, tree2, child12, 1, child22, 1,
897 predecessors1, predecessors2, depth1, depth2,
898 memT, mapping);
899
900 return;
901 }
902 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
903 == memT[child12 + (l1 + 1) * dim2 + child22 * dim3
904 + (l2 + 1) * dim4]
905 + memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
906
907 traceMapping_branch(tree1, tree2, child12, (l1 + 1), child22,
908 (l2 + 1), predecessors1, predecessors2, depth1,
909 depth2, memT, mapping);
910 traceMapping_branch(tree1, tree2, child11, 1, child21, 1,
911 predecessors1, predecessors2, depth1, depth2,
912 memT, mapping);
913
914 return;
915 }
916 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
917 == memT[child11 + (l1 + 1) * dim2 + child22 * dim3
918 + (l2 + 1) * dim4]
919 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
920
921 traceMapping_branch(tree1, tree2, child11, (l1 + 1), child22,
922 (l2 + 1), predecessors1, predecessors2, depth1,
923 depth2, memT, mapping);
924 traceMapping_branch(tree1, tree2, child12, 1, child21, 1,
925 predecessors1, predecessors2, depth1, depth2,
926 memT, mapping);
927
928 return;
929 }
930 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
931 == memT[child12 + (l1 + 1) * dim2 + child21 * dim3
932 + (l2 + 1) * dim4]
933 + memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
934
935 traceMapping_branch(tree1, tree2, child12, (l1 + 1), child21,
936 (l2 + 1), predecessors1, predecessors2, depth1,
937 depth2, memT, mapping);
938 traceMapping_branch(tree1, tree2, child11, 1, child22, 1,
939 predecessors1, predecessors2, depth1, depth2,
940 memT, mapping);
941
942 return;
943 }
944 } else {
945 for(auto child1_mb : children1) {
946 std::vector<ftm::idNode> topo1_;
947 tree1->getChildren(curr1, topo1_);
948 topo1_.erase(std::remove(topo1_.begin(), topo1_.end(), child1_mb),
949 topo1_.end());
950 for(auto child2_mb : children2) {
951 std::vector<ftm::idNode> topo2_;
952 tree2->getChildren(curr2, topo2_);
953 topo2_.erase(std::remove(topo2_.begin(), topo2_.end(), child2_mb),
954 topo2_.end());
955
956 auto f = [&](unsigned r, unsigned c) {
957 int c1 = r < topo1_.size() ? topo1_[r] : -1;
958 int c2 = c < topo2_.size() ? topo2_[c] : -1;
959 return memT[c1 + 1 * dim2 + c2 * dim3 + 1 * dim4];
960 };
961 int size = std::max(topo1_.size(), topo2_.size()) + 1;
962 auto costMatrix = std::vector<std::vector<dataType>>(
963 size, std::vector<dataType>(size, 0));
964 std::vector<MatchingType> matching;
965 for(int r = 0; r < size; r++) {
966 for(int c = 0; c < size; c++) {
967 costMatrix[r][c] = f(r, c);
968 }
969 }
970
971 AssignmentSolver<dataType> *assignmentSolver;
972 AssignmentExhaustive<dataType> solverExhaustive;
973 AssignmentMunkres<dataType> solverMunkres;
974 AssignmentAuction<dataType> solverAuction;
975 switch(assignmentSolverID_) {
976 case 1:
977 solverExhaustive = AssignmentExhaustive<dataType>();
978 assignmentSolver = &solverExhaustive;
979 break;
980 case 2:
981 solverMunkres = AssignmentMunkres<dataType>();
982 assignmentSolver = &solverMunkres;
983 break;
984 case 0:
985 default:
986 solverAuction = AssignmentAuction<dataType>();
987 assignmentSolver = &solverAuction;
988 }
989 assignmentSolver->setInput(costMatrix);
990 assignmentSolver->setBalanced(true);
991 assignmentSolver->run(matching);
992 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + child2_mb * dim3
993 + (l2 + 1) * dim4];
994 for(auto m : matching)
995 d_ += std::get<2>(m);
996
997 if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
999 tree1, tree2, child1_mb, (l1 + 1), child2_mb, (l2 + 1),
1000 predecessors1, predecessors2, depth1, depth2, memT, mapping);
1001 for(auto m : matching) {
1002 int n1 = std::get<0>(m) < static_cast<int>(topo1_.size())
1003 ? topo1_[std::get<0>(m)]
1004 : -1;
1005 int n2 = std::get<1>(m) < static_cast<int>(topo2_.size())
1006 ? topo2_[std::get<1>(m)]
1007 : -1;
1008 if(n1 >= 0 && n2 >= 0)
1009 traceMapping_branch(tree1, tree2, n1, 1, n2, 1,
1010 predecessors1, predecessors2, depth1,
1011 depth2, memT, mapping);
1012 else if(n1 >= 0)
1013 traceMapping_branch(tree1, tree2, n1, 1, nn2, 0,
1014 predecessors1, predecessors2, depth1,
1015 depth2, memT, mapping);
1016 else if(n2 >= 0)
1017 traceMapping_branch(tree1, tree2, nn1, 0, n2, 1,
1018 predecessors1, predecessors2, depth1,
1019 depth2, memT, mapping);
1020 }
1021 return;
1022 }
1023 }
1024 }
1025 }
1026 //-----------------------------------------------------------------------
1027 // Try to continue main branch on one child of first tree and
1028 // delete all other subtrees Then match continued branch to
1029 // current branch in second tree
1030 for(auto child1_mb : children1) {
1031 dataType d_
1032 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
1033 for(auto child1 : children1) {
1034 if(child1 == child1_mb) {
1035 continue;
1036 }
1037 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
1038 }
1039 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
1040 traceMapping_branch(tree1, tree2, child1_mb, (l1 + 1), curr2, l2,
1041 predecessors1, predecessors2, depth1, depth2,
1042 memT, mapping);
1043 for(auto child1 : children1) {
1044 if(child1 == child1_mb) {
1045 continue;
1046 }
1047 traceMapping_branch(tree1, tree2, child1, 1, nn2, 0,
1048 predecessors1, predecessors2, depth1, depth2,
1049 memT, mapping);
1050 }
1051 return;
1052 }
1053 }
1054 //-----------------------------------------------------------------------
1055 // Try to continue main branch on one child of second tree and
1056 // delete all other subtrees Then match continued branch to
1057 // current branch in first tree
1058 for(auto child2_mb : children2) {
1059 dataType d_
1060 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
1061 for(auto child2 : children2) {
1062 if(child2 == child2_mb) {
1063 continue;
1064 }
1065 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
1066 }
1067 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
1068 traceMapping_branch(tree1, tree2, curr1, l1, child2_mb, (l2 + 1),
1069 predecessors1, predecessors2, depth1, depth2,
1070 memT, mapping);
1071 for(auto child2 : children2) {
1072 if(child2 == child2_mb) {
1073 continue;
1074 }
1075 traceMapping_branch(tree1, tree2, nn1, 0, child2, 1,
1076 predecessors1, predecessors2, depth1, depth2,
1077 memT, mapping);
1078 }
1079 return;
1080 }
1081 }
1082 this->printErr("Mapping traceback not correct");
1083 }
1084 }
1085 };
1086
1087} // namespace ttk
virtual int run(std::vector< MatchingType > &matchings)=0
virtual int setInput(std::vector< std::vector< dataType > > &C_)
virtual void setBalanced(bool balanced)
void traceMapping_branch(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2, int curr1, int l1, int curr2, int l2, std::vector< std::vector< int > > &predecessors1, std::vector< std::vector< int > > &predecessors2, int depth1, int depth2, std::vector< dataType > &memT, std::vector< std::pair< std::pair< int, int >, std::pair< int, int > > > &mapping)
dataType computeDistance(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > *outputMatching=nullptr)
~BranchMappingDistance() override=default
dataType execute(ftm::MergeTree< dataType > &mTree1, ftm::MergeTree< dataType > &mTree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > *outputMatching=nullptr)
void setAssignmentSolver(int assignmentSolver)
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
void preprocessingPipeline(ftm::MergeTree< dataType > &mTree, double epsilonTree, double epsilon2Tree, double epsilon3Tree, bool branchDecompositionT, bool useMinMaxPairT, bool cleanTreeT, double persistenceThreshold, std::vector< int > &nodeCorr, bool deleteInconsistentNodes=true, bool removeMergedSaddles=false)
std::vector< std::vector< int > > treesNodeCorr_
Node * getNode(idNode nodeId) const
Definition FTMTree_MT.h:393
const scalarType & getValue(SimplexId nodeId) const
Definition FTMTree_MT.h:339
void getChildren(idNode nodeId, std::vector< idNode > &res) const
idNode getNumberOfNodes() const
Definition FTMTree_MT.h:389
idNode getRoot() const
int getNumberOfChildren(idNode nodeId) const
void setOrigin(SimplexId linked)
Definition FTMNode.h:72
MergeTree< dataType > copyMergeTree(const ftm::FTMTree_MT *tree, bool doSplitMultiPersPairs=false)
unsigned int idNode
Node index in vect_nodes_.
TTK base package defining the standard types.
T end(std::pair< T, T > &p)
Definition ripser.cpp:503
T begin(std::pair< T, T > &p)
Definition ripser.cpp:499
ftm::FTMTree_MT tree
Definition FTMTree_MT.h:906