TTK
Loading...
Searching...
No Matches
PathMappingDistance.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 PathMappingDistance : 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
48 bool preprocess_ = true;
49 bool saveTree_ = false;
50
51 template <class dataType>
52 inline dataType editCost_Persistence(int n1,
53 int p1,
54 int n2,
55 int p2,
56 ftm::FTMTree_MT *tree1,
57 ftm::FTMTree_MT *tree2) {
58 dataType d;
59 if(n1 < 0) {
60 dataType b1 = tree2->getValue<dataType>(n2);
61 dataType d1 = tree2->getValue<dataType>(p2);
62 d = (d1 > b1) ? (d1 - b1) : (b1 - d1);
63 } else if(n2 < 0) {
64 dataType b1 = tree1->getValue<dataType>(n1);
65 dataType d1 = tree1->getValue<dataType>(p1);
66 d = (d1 > b1) ? (d1 - b1) : (b1 - d1);
67 } else {
68 dataType b1 = tree1->getValue<dataType>(n1);
69 dataType d1 = tree1->getValue<dataType>(p1);
70 dataType b2 = tree2->getValue<dataType>(n2);
71 dataType d2 = tree2->getValue<dataType>(p2);
72 dataType dist1 = (d1 > b1) ? (d1 - b1) : (b1 - d1);
73 dataType dist2 = (d2 > b2) ? (d2 - b2) : (b2 - d2);
74 d = (dist1 > dist2) ? (dist1 - dist2) : (dist2 - dist1);
75 }
76 return squared_ ? d * d : d;
77 }
78
79 template <class dataType>
80 void traceMapping_path(
81 ftm::FTMTree_MT *tree1,
82 ftm::FTMTree_MT *tree2,
83 int curr1,
84 int l1,
85 int curr2,
86 int l2,
87 std::vector<std::vector<int>> &predecessors1,
88 std::vector<std::vector<int>> &predecessors2,
89 int depth1,
90 int depth2,
91 std::vector<dataType> &memT,
92 std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
93 std::pair<ftm::idNode, ftm::idNode>>> &mapping) {
94
95 //===============================================================================
96 // If both trees not empty, find optimal edit operation
97 std::vector<ftm::idNode> children1;
98 tree1->getChildren(curr1, children1);
99 std::vector<ftm::idNode> children2;
100 tree2->getChildren(curr2, children2);
101 int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
102 int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
103
104 size_t const nn1 = tree1->getNumberOfNodes();
105 size_t const nn2 = tree2->getNumberOfNodes();
106 size_t const dim1 = 1;
107 size_t const dim2 = (nn1 + 1) * dim1;
108 size_t const dim3 = (depth1 + 1) * dim2;
109 size_t const dim4 = (nn2 + 1) * dim3;
110
111 //---------------------------------------------------------------------------
112 // If both trees only have one branch, return edit cost between
113 // the two branches
114 if(tree1->getNumberOfChildren(curr1) == 0
115 and tree2->getNumberOfChildren(curr2) == 0) {
116 mapping.emplace_back(
117 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2));
118 return;
119 }
120 //---------------------------------------------------------------------------
121 // If first tree only has one branch, try all decompositions of
122 // second tree
123 else if(tree1->getNumberOfChildren(curr1) == 0) {
124 for(auto child2_mb : children2) {
125 dataType d_
126 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
127 for(auto child2 : children2) {
128 if(child2 == child2_mb) {
129 continue;
130 }
131 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
132 }
133 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
134 traceMapping_path(tree1, tree2, curr1, l1, child2_mb, l2 + 1,
135 predecessors1, predecessors2, depth1, depth2,
136 memT, mapping);
137 return;
138 }
139 }
140 }
141 //---------------------------------------------------------------------------
142 // If second tree only has one branch, try all decompositions of
143 // first tree
144 else if(tree2->getNumberOfChildren(curr2) == 0) {
145 dataType d = std::numeric_limits<dataType>::max();
146 for(auto child1_mb : children1) {
147 dataType d_
148 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
149 for(auto child1 : children1) {
150 if(child1 == child1_mb) {
151 continue;
152 }
153 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
154 }
155 d = std::min(d, d_);
156 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
157 traceMapping_path(tree1, tree2, child1_mb, l1 + 1, curr2, l2,
158 predecessors1, predecessors2, depth1, depth2,
159 memT, mapping);
160 return;
161 }
162 }
163 }
164 //---------------------------------------------------------------------------
165 // If both trees have more than one branch, try all decompositions
166 // of both trees
167 else {
168 //-----------------------------------------------------------------------
169 // Try all possible main branches of first tree (child1_mb) and
170 // all possible main branches of second tree (child2_mb) Then
171 // try all possible matchings of subtrees
172 if(tree1->getNumberOfChildren(curr1) == 2
173 && tree2->getNumberOfChildren(curr2) == 2) {
174 int child11 = children1[0];
175 int child12 = children1[1];
176 int child21 = children2[0];
177 int child22 = children2[1];
178 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
179 == memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]
180 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]
181 + editCost_Persistence<dataType>(
182 curr1, parent1, curr2, parent2, tree1, tree2)) {
183 mapping.emplace_back(
184 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2));
185 traceMapping_path(tree1, tree2, child11, 1, child21, 1,
186 predecessors1, predecessors2, depth1, depth2,
187 memT, mapping);
188 traceMapping_path(tree1, tree2, child12, 1, child22, 1,
189 predecessors1, predecessors2, depth1, depth2,
190 memT, mapping);
191 return;
192 }
193 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
194 == memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]
195 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]
196 + editCost_Persistence<dataType>(
197 curr1, parent1, curr2, parent2, tree1, tree2)) {
198 mapping.emplace_back(
199 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2));
200 traceMapping_path(tree1, tree2, child11, 1, child22, 1,
201 predecessors1, predecessors2, depth1, depth2,
202 memT, mapping);
203 traceMapping_path(tree1, tree2, child12, 1, child21, 1,
204 predecessors1, predecessors2, depth1, depth2,
205 memT, mapping);
206 return;
207 }
208 } else {
209 auto f = [&](int r, int c) {
210 size_t const c1
211 = r < tree1->getNumberOfChildren(curr1) ? children1[r] : nn1;
212 size_t const c2
213 = c < tree2->getNumberOfChildren(curr2) ? children2[c] : nn2;
214 int const l1_ = c1 == nn1 ? 0 : 1;
215 int const l2_ = c2 == nn2 ? 0 : 1;
216 return memT[c1 + l1_ * dim2 + c2 * dim3 + l2_ * dim4];
217 };
218 int size = std::max(tree1->getNumberOfChildren(curr1),
219 tree2->getNumberOfChildren(curr2))
220 + 1;
221 auto costMatrix = std::vector<std::vector<dataType>>(
222 size, std::vector<dataType>(size, 0));
223 std::vector<MatchingType> matching;
224 for(int r = 0; r < size; r++) {
225 for(int c = 0; c < size; c++) {
226 costMatrix[r][c] = f(r, c);
227 }
228 }
229
230 AssignmentSolver<dataType> *assignmentSolver;
231 AssignmentExhaustive<dataType> solverExhaustive;
232 AssignmentMunkres<dataType> solverMunkres;
233 AssignmentAuction<dataType> solverAuction;
234 switch(assignmentSolverID_) {
235 case 1:
236 solverExhaustive = AssignmentExhaustive<dataType>();
237 assignmentSolver = &solverExhaustive;
238 break;
239 case 2:
240 solverMunkres = AssignmentMunkres<dataType>();
241 assignmentSolver = &solverMunkres;
242 break;
243 case 0:
244 default:
245 solverAuction = AssignmentAuction<dataType>();
246 assignmentSolver = &solverAuction;
247 }
248 assignmentSolver->setInput(costMatrix);
249 assignmentSolver->setBalanced(true);
250 assignmentSolver->run(matching);
251 dataType d_ = editCost_Persistence<dataType>(
252 curr1, parent1, curr2, parent2, tree1, tree2);
253 for(auto m : matching)
254 d_ += std::get<2>(m);
255 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
256 mapping.emplace_back(
257 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2));
258 for(auto m : matching) {
259 int n1 = std::get<0>(m) < tree1->getNumberOfChildren(curr1)
260 ? children1[std::get<0>(m)]
261 : -1;
262 int n2 = std::get<1>(m) < tree2->getNumberOfChildren(curr2)
263 ? children2[std::get<1>(m)]
264 : -1;
265 if(n1 >= 0 && n2 >= 0)
266 traceMapping_path(tree1, tree2, n1, 1, n2, 1, predecessors1,
267 predecessors2, depth1, depth2, memT, mapping);
268 }
269 return;
270 }
271 }
272 //-----------------------------------------------------------------------
273 // Try to continue main branch on one child of first tree and
274 // delete all other subtrees Then match continued branch to
275 // current branch in second tree
276 for(auto child1_mb : children1) {
277 dataType d_
278 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
279 for(auto child1 : children1) {
280 if(child1 == child1_mb) {
281 continue;
282 }
283 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
284 }
285 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
286 traceMapping_path(tree1, tree2, child1_mb, l1 + 1, curr2, l2,
287 predecessors1, predecessors2, depth1, depth2,
288 memT, mapping);
289 return;
290 }
291 }
292 //-----------------------------------------------------------------------
293 // Try to continue main branch on one child of second tree and
294 // delete all other subtrees Then match continued branch to
295 // current branch in first tree
296 for(auto child2_mb : children2) {
297 dataType d_
298 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
299 for(auto child2 : children2) {
300 if(child2 == child2_mb) {
301 continue;
302 }
303 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
304 }
305 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
306 traceMapping_path(tree1, tree2, curr1, l1, child2_mb, l2 + 1,
307 predecessors1, predecessors2, depth1, depth2,
308 memT, mapping);
309 return;
310 }
311 }
312 }
313 }
314
315 public:
317 this->setDebugMsgPrefix(
318 "MergeTreeDistance"); // inherited from Debug: prefix will be printed at
319 // the beginning of every msg
320 }
321 ~PathMappingDistance() override = default;
322
323 void setBaseMetric(int m) {
324 baseMetric_ = m;
325 }
326
327 void setAssignmentSolver(int assignmentSolver) {
328 assignmentSolverID_ = assignmentSolver;
329 }
330
331 void setSquared(bool s) {
332 squared_ = s;
333 }
334
335 void setComputeMapping(bool m) {
336 computeMapping_ = m;
337 }
338
339 void setPreprocess(bool p) {
340 preprocess_ = p;
341 }
342
343 void setSaveTree(bool save) {
344 saveTree_ = save;
345 }
346
347 template <class dataType>
349 ftm::FTMTree_MT *tree1,
350 ftm::FTMTree_MT *tree2,
351 std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
352 std::pair<ftm::idNode, ftm::idNode>>>
353 *outputMatching) {
354
355 // compute preorder of both trees (necessary for bottom-up dynamic
356 // programming)
357
358 std::vector<std::vector<int>> predecessors1(tree1->getNumberOfNodes());
359 std::vector<std::vector<int>> predecessors2(tree2->getNumberOfNodes());
360 int const rootID1 = tree1->getRoot();
361 int const rootID2 = tree2->getRoot();
362 std::vector<int> preorder1(tree1->getNumberOfNodes());
363 std::vector<int> preorder2(tree2->getNumberOfNodes());
364
365 int depth1 = 0;
366 int depth2 = 0;
367 std::stack<int> stack;
368 stack.push(rootID1);
369 int count = tree1->getNumberOfNodes() - 1;
370 while(!stack.empty()) {
371 int const nIdx = stack.top();
372 stack.pop();
373 preorder1[count] = nIdx;
374 count--;
375 depth1 = std::max((int)predecessors1[nIdx].size(), depth1);
376 std::vector<ftm::idNode> children;
377 tree1->getChildren(nIdx, children);
378 for(int const cIdx : children) {
379 stack.push(cIdx);
380 predecessors1[cIdx].reserve(predecessors1[nIdx].size() + 1);
381 predecessors1[cIdx].insert(predecessors1[cIdx].end(),
382 predecessors1[nIdx].begin(),
383 predecessors1[nIdx].end());
384 predecessors1[cIdx].push_back(nIdx);
385 }
386 }
387 stack.push(rootID2);
388 count = tree2->getNumberOfNodes() - 1;
389 while(!stack.empty()) {
390 int const nIdx = stack.top();
391 stack.pop();
392 preorder2[count] = nIdx;
393 count--;
394 depth2 = std::max((int)predecessors2[nIdx].size(), depth2);
395 std::vector<ftm::idNode> children;
396 tree2->getChildren(nIdx, children);
397 for(int const cIdx : children) {
398 stack.push(cIdx);
399 predecessors2[cIdx].reserve(predecessors2[nIdx].size() + 1);
400 predecessors2[cIdx].insert(predecessors2[cIdx].end(),
401 predecessors2[nIdx].begin(),
402 predecessors2[nIdx].end());
403 predecessors2[cIdx].push_back(nIdx);
404 }
405 }
406
407 // initialize memoization tables
408
409 size_t nn1 = tree1->getNumberOfNodes();
410 size_t nn2 = tree2->getNumberOfNodes();
411 size_t const dim1 = 1;
412 size_t const dim2 = (nn1 + 1) * dim1;
413 size_t const dim3 = (depth1 + 1) * dim2;
414 size_t const dim4 = (nn2 + 1) * dim3;
415
416 // std::cout << (nn1 + 1) * (depth1 + 1) * (nn2 + 1) * (depth2 + 1) *
417 // sizeof(dataType) << std::endl;
418 std::vector<dataType> memT((nn1 + 1) * (depth1 + 1) * (nn2 + 1)
419 * (depth2 + 1));
420
421 memT[nn1 + 0 * dim2 + nn2 * dim3 + 0 * dim4] = 0;
422 for(size_t i = 0; i < nn1; i++) {
423 int curr1 = preorder1[i];
424 std::vector<ftm::idNode> children1;
425 tree1->getChildren(curr1, children1);
426 for(size_t l = 1; l <= predecessors1[preorder1[i]].size(); l++) {
427 int parent1 = predecessors1[preorder1[i]]
428 [predecessors1[preorder1[i]].size() - l];
429
430 //-----------------------------------------------------------------------
431 // Delete curr path and full subtree rooted in path
432 memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4]
433 = editCost_Persistence<dataType>(
434 curr1, parent1, -1, -1, tree1, tree2);
435 for(auto child1 : children1) {
436 memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4]
437 += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
438 }
439 }
440 }
441 for(size_t j = 0; j < nn2; j++) {
442 int curr2 = preorder2[j];
443 std::vector<ftm::idNode> children2;
444 tree2->getChildren(curr2, children2);
445 for(size_t l = 1; l <= predecessors2[preorder2[j]].size(); l++) {
446 int parent2 = predecessors2[preorder2[j]]
447 [predecessors2[preorder2[j]].size() - l];
448
449 //-----------------------------------------------------------------------
450 // Delete curr path and full subtree rooted in path
451 memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4]
452 = editCost_Persistence<dataType>(
453 -1, -1, curr2, parent2, tree1, tree2);
454 for(auto child2 : children2) {
455 memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4]
456 += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
457 }
458 }
459 }
460
461 for(size_t i = 0; i < nn1; i++) {
462 int curr1 = preorder1[i];
463 std::vector<ftm::idNode> children1;
464 tree1->getChildren(curr1, children1);
465 for(size_t j = 0; j < nn2; j++) {
466 int curr2 = preorder2[j];
467 std::vector<ftm::idNode> children2;
468 tree2->getChildren(curr2, children2);
469 for(size_t l1 = 1; l1 <= predecessors1[preorder1[i]].size(); l1++) {
470 int parent1
471 = predecessors1[preorder1[i]]
472 [predecessors1[preorder1[i]].size() - l1];
473 for(size_t l2 = 1; l2 <= predecessors2[preorder2[j]].size(); l2++) {
474 int parent2
475 = predecessors2[preorder2[j]]
476 [predecessors2[preorder2[j]].size() - l2];
477
478 //===============================================================================
479 // If both trees not empty, find optimal edit operation
480
481 //---------------------------------------------------------------------------
482 // If both trees only have one branch, return edit cost between
483 // the two branches
484 if(tree1->getNumberOfChildren(curr1) == 0
485 and tree2->getNumberOfChildren(curr2) == 0) {
486 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
487 = editCost_Persistence<dataType>(
488 curr1, parent1, curr2, parent2, tree1, tree2);
489 }
490 //---------------------------------------------------------------------------
491 // If first tree only has one branch, try all decompositions of
492 // second tree
493 else if(tree1->getNumberOfChildren(curr1) == 0) {
494 dataType d = std::numeric_limits<dataType>::max();
495 for(auto child2_mb : children2) {
496 dataType d_ = memT[curr1 + l1 * dim2 + child2_mb * dim3
497 + (l2 + 1) * dim4];
498 for(auto child2 : children2) {
499 if(child2 == child2_mb) {
500 continue;
501 }
502 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
503 }
504 d = std::min(d, d_);
505 }
506 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
507 }
508 //---------------------------------------------------------------------------
509 // If second tree only has one branch, try all decompositions of
510 // first tree
511 else if(tree2->getNumberOfChildren(curr2) == 0) {
512 dataType d = std::numeric_limits<dataType>::max();
513 for(auto child1_mb : children1) {
514 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3
515 + l2 * dim4];
516 for(auto child1 : children1) {
517 if(child1 == child1_mb) {
518 continue;
519 }
520 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
521 }
522 d = std::min(d, d_);
523 }
524 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
525 }
526 //---------------------------------------------------------------------------
527 // If both trees have more than one branch, try all decompositions
528 // of both trees
529 else {
530 dataType d = std::numeric_limits<dataType>::max();
531 //-----------------------------------------------------------------------
532 // Try all possible main branches of first tree (child1_mb) and
533 // all possible main branches of second tree (child2_mb) Then
534 // try all possible matchings of subtrees
535 if(tree1->getNumberOfChildren(curr1) == 2
536 && tree2->getNumberOfChildren(curr2) == 2) {
537 int const child11 = children1[0];
538 int const child12 = children1[1];
539 int const child21 = children2[0];
540 int const child22 = children2[1];
541 d = std::min<dataType>(
542 d, memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]
543 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]
544 + editCost_Persistence<dataType>(
545 curr1, parent1, curr2, parent2, tree1, tree2));
546 d = std::min<dataType>(
547 d, memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]
548 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]
549 + editCost_Persistence<dataType>(
550 curr1, parent1, curr2, parent2, tree1, tree2));
551 } else {
552 auto f = [&](int r, int c) {
553 size_t const c1 = r < tree1->getNumberOfChildren(curr1)
554 ? children1[r]
555 : nn1;
556 size_t const c2 = c < tree2->getNumberOfChildren(curr2)
557 ? children2[c]
558 : nn2;
559 int const l1_ = c1 == nn1 ? 0 : 1;
560 int const l2_ = c2 == nn2 ? 0 : 1;
561 return memT[c1 + l1_ * dim2 + c2 * dim3 + l2_ * dim4];
562 };
563 int size = std::max(tree1->getNumberOfChildren(curr1),
564 tree2->getNumberOfChildren(curr2))
565 + 1;
566 auto costMatrix = std::vector<std::vector<dataType>>(
567 size, std::vector<dataType>(size, 0));
568 std::vector<MatchingType> matching;
569 for(int r = 0; r < size; r++) {
570 for(int c = 0; c < size; c++) {
571 costMatrix[r][c] = f(r, c);
572 }
573 }
574
575 AssignmentSolver<dataType> *assignmentSolver;
576 AssignmentExhaustive<dataType> solverExhaustive;
577 AssignmentMunkres<dataType> solverMunkres;
578 AssignmentAuction<dataType> solverAuction;
579 switch(assignmentSolverID_) {
580 case 1:
581 solverExhaustive = AssignmentExhaustive<dataType>();
582 assignmentSolver = &solverExhaustive;
583 break;
584 case 2:
585 solverMunkres = AssignmentMunkres<dataType>();
586 assignmentSolver = &solverMunkres;
587 break;
588 case 0:
589 default:
590 solverAuction = AssignmentAuction<dataType>();
591 assignmentSolver = &solverAuction;
592 }
593 assignmentSolver->setInput(costMatrix);
594 assignmentSolver->setBalanced(true);
595 assignmentSolver->run(matching);
596 dataType d_ = editCost_Persistence<dataType>(
597 curr1, parent1, curr2, parent2, tree1, tree2);
598 for(auto m : matching)
599 d_ += std::get<2>(m);
600 d = std::min(d, d_);
601 }
602 //-----------------------------------------------------------------------
603 // Try to continue main branch on one child of first tree and
604 // delete all other subtrees Then match continued branch to
605 // current branch in second tree
606 for(auto child1_mb : children1) {
607 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3
608 + l2 * dim4];
609 for(auto child1 : children1) {
610 if(child1 == child1_mb) {
611 continue;
612 }
613 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
614 }
615 d = std::min(d, d_);
616 }
617 //-----------------------------------------------------------------------
618 // Try to continue main branch on one child of second tree and
619 // delete all other subtrees Then match continued branch to
620 // current branch in first tree
621 for(auto child2_mb : children2) {
622 dataType d_ = memT[curr1 + l1 * dim2 + child2_mb * dim3
623 + (l2 + 1) * dim4];
624 for(auto child2 : children2) {
625 if(child2 == child2_mb) {
626 continue;
627 }
628 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
629 }
630 d = std::min(d, d_);
631 }
632 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
633 }
634 }
635 }
636 }
637 }
638
639 std::vector<ftm::idNode> children1;
640 tree1->getChildren(rootID1, children1);
641 std::vector<ftm::idNode> children2;
642 tree2->getChildren(rootID2, children2);
643
644 dataType res
645 = memT[children1[0] + 1 * dim2 + children2[0] * dim3 + 1 * dim4];
646
647 if(computeMapping_ && outputMatching) {
648
649 outputMatching->clear();
650 traceMapping_path(tree1, tree2, children1[0], 1, children2[0], 1,
651 predecessors1, predecessors2, depth1, depth2, memT,
652 *outputMatching);
653 }
654
655 return squared_ ? std::sqrt(res) : res;
656 }
657
658 template <class dataType>
661 std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
662 std::pair<ftm::idNode, ftm::idNode>>>
663 *outputMatching) {
664
665 ftm::MergeTree<dataType> mTree1Copy;
666 ftm::MergeTree<dataType> mTree2Copy;
667 if(saveTree_) {
668 mTree1Copy = ftm::copyMergeTree<dataType>(mTree1);
669 mTree2Copy = ftm::copyMergeTree<dataType>(mTree2);
670 }
671 ftm::MergeTree<dataType> &mTree1Int = (saveTree_ ? mTree1Copy : mTree1);
672 ftm::MergeTree<dataType> &mTree2Int = (saveTree_ ? mTree2Copy : mTree2);
673 ftm::FTMTree_MT *tree1 = &(mTree1Int.tree);
674 ftm::FTMTree_MT *tree2 = &(mTree2Int.tree);
675
676 // optional preprocessing
677 if(preprocess_) {
678 treesNodeCorr_.resize(2);
682 true, true);
686 true, true);
687 }
688
689 tree1 = &(mTree1Int.tree);
690 tree2 = &(mTree2Int.tree);
691
692 return computeDistance<dataType>(tree1, tree2, outputMatching);
693 }
694
695 template <class dataType>
696 dataType
698 ftm::FTMTree_MT *tree2,
699 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>
700 *outputMatching) {
701
702 std::vector<int> matchedNodes(tree1->getNumberOfNodes(), -1);
703 std::vector<double> matchedCost(tree1->getNumberOfNodes(), -1);
704 std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
705 std::pair<ftm::idNode, ftm::idNode>>>
706 mapping;
707 dataType res = computeDistance<dataType>(tree1, tree2, &mapping);
708 if(computeMapping_ && outputMatching) {
709 outputMatching->clear();
710 for(auto m : mapping) {
711 matchedNodes[m.first.first] = m.second.first;
712 matchedNodes[m.first.second] = m.second.second;
713 matchedCost[m.first.first] = editCost_Persistence<dataType>(
714 m.first.first, m.first.second, m.second.first, m.second.second,
715 tree1, tree2);
716 if(m.first.second == tree1->getRoot()) {
717 matchedCost[m.first.second] = matchedCost[m.first.first];
718 }
719 }
720 for(ftm::idNode i = 0; i < matchedNodes.size(); i++) {
721 if(matchedNodes[i] >= 0) {
722 outputMatching->emplace_back(
723 std::make_tuple(i, matchedNodes[i], matchedCost[i]));
724 }
725 }
726 }
727
728 return res;
729 }
730
731 template <class dataType>
733 ftm::FTMTree_MT *tree1,
734 ftm::FTMTree_MT *tree2,
735 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>> *outputMatching,
736 std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
737 std::pair<ftm::idNode, ftm::idNode>>>
738 *outputMatching_path) {
739
740 std::vector<int> matchedNodes(tree1->getNumberOfNodes(), -1);
741 std::vector<double> matchedCost(tree1->getNumberOfNodes(), -1);
742 dataType res
743 = computeDistance<dataType>(tree1, tree2, outputMatching_path);
744 if(computeMapping_ && outputMatching) {
745 outputMatching->clear();
746 for(auto m : *outputMatching_path) {
747 matchedNodes[m.first.first] = m.second.first;
748 matchedNodes[m.first.second] = m.second.second;
749 matchedCost[m.first.first] = editCost_Persistence<dataType>(
750 m.first.first, m.first.second, m.second.first, m.second.second,
751 tree1, tree2);
752 if(m.first.second == tree1->getRoot()) {
753 matchedCost[m.first.second] = matchedCost[m.first.first];
754 }
755 }
756 for(ftm::idNode i = 0; i < matchedNodes.size(); i++) {
757 if(matchedNodes[i] >= 0) {
758 outputMatching->emplace_back(
759 std::make_tuple(i, matchedNodes[i], matchedCost[i]));
760 }
761 }
762 }
763
764 return res;
765 }
766
767 template <class dataType>
770 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>
771 *outputMatching) {
772
773 ftm::MergeTree<dataType> mTree1Copy;
774 ftm::MergeTree<dataType> mTree2Copy;
775 if(saveTree_) {
776 mTree1Copy = ftm::copyMergeTree<dataType>(mTree1);
777 mTree2Copy = ftm::copyMergeTree<dataType>(mTree2);
778 }
779 ftm::MergeTree<dataType> &mTree1Int = (saveTree_ ? mTree1Copy : mTree1);
780 ftm::MergeTree<dataType> &mTree2Int = (saveTree_ ? mTree2Copy : mTree2);
781 ftm::FTMTree_MT *tree1 = &(mTree1Int.tree);
782 ftm::FTMTree_MT *tree2 = &(mTree2Int.tree);
783
784 // optional preprocessing
785 if(preprocess_) {
786 treesNodeCorr_.resize(2);
788 mTree1Int, epsilonTree1_, epsilon2Tree1_, epsilon3Tree1_, false,
789 useMinMaxPair_, cleanTree_, treesNodeCorr_[0], true, true);
791 mTree2Int, epsilonTree2_, epsilon2Tree2_, epsilon3Tree2_, false,
792 useMinMaxPair_, cleanTree_, treesNodeCorr_[1], true, true);
793 }
794
795 tree1 = &(mTree1Int.tree);
796 tree2 = &(mTree2Int.tree);
797
798 return computeDistance<dataType>(tree1, tree2, outputMatching);
799 }
800
801 template <class dataType>
804 tree1, tree2,
805 (std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
806 std::pair<ftm::idNode, ftm::idNode>>> *)nullptr);
807 }
808
809 template <class dataType>
811 ftm::MergeTree<dataType> &mTree2) {
812
813 ftm::MergeTree<dataType> mTree1Copy;
814 ftm::MergeTree<dataType> mTree2Copy;
815 if(saveTree_) {
816 mTree1Copy = ftm::copyMergeTree<dataType>(mTree1);
817 mTree2Copy = ftm::copyMergeTree<dataType>(mTree2);
818 }
819 ftm::MergeTree<dataType> &mTree1Int = (saveTree_ ? mTree1Copy : mTree1);
820 ftm::MergeTree<dataType> &mTree2Int = (saveTree_ ? mTree2Copy : mTree2);
821 ftm::FTMTree_MT *tree1 = &(mTree1Int.tree);
822 ftm::FTMTree_MT *tree2 = &(mTree2Int.tree);
823
824 // optional preprocessing
825 if(preprocess_) {
826 treesNodeCorr_.resize(2);
828 mTree1Int, epsilonTree1_, epsilon2Tree1_, epsilon3Tree1_, false,
829 useMinMaxPair_, cleanTree_, treesNodeCorr_[0], true, true);
831 mTree2Int, epsilonTree2_, epsilon2Tree2_, epsilon3Tree2_, false,
832 useMinMaxPair_, cleanTree_, treesNodeCorr_[1], true, true);
833 }
834
835 tree1 = &(mTree1Int.tree);
836 tree2 = &(mTree2Int.tree);
837
839 tree1, tree2,
840 (std::vector<std::pair<std::pair<ftm::idNode, ftm::idNode>,
841 std::pair<ftm::idNode, ftm::idNode>>> *)nullptr);
842 }
843 };
844} // 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 setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
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_
dataType execute(ftm::MergeTree< dataType > &mTree1, ftm::MergeTree< dataType > &mTree2)
void setAssignmentSolver(int assignmentSolver)
dataType computeDistance(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2)
~PathMappingDistance() override=default
dataType execute(ftm::MergeTree< dataType > &mTree1, ftm::MergeTree< dataType > &mTree2, std::vector< std::pair< std::pair< ftm::idNode, ftm::idNode >, std::pair< ftm::idNode, ftm::idNode > > > *outputMatching)
dataType computeDistance(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > *outputMatching)
dataType computeDistance(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > *outputMatching, std::vector< std::pair< std::pair< ftm::idNode, ftm::idNode >, std::pair< ftm::idNode, ftm::idNode > > > *outputMatching_path)
dataType computeDistance(ftm::FTMTree_MT *tree1, ftm::FTMTree_MT *tree2, std::vector< std::pair< std::pair< ftm::idNode, ftm::idNode >, std::pair< ftm::idNode, ftm::idNode > > > *outputMatching)
dataType execute(ftm::MergeTree< dataType > &mTree1, ftm::MergeTree< dataType > &mTree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > *outputMatching)
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
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