257 std::vector<std::tuple<ftm::idNode, ftm::idNode, double>> *outputMatching
265 int const rootID1 = tree1->
getRoot();
266 int const rootID2 = tree2->
getRoot();
272 std::stack<int> stack;
275 while(!stack.empty()) {
276 int const nIdx = stack.top();
278 preorder1[count] = nIdx;
280 depth1 = std::max((
int)predecessors1[nIdx].size(), depth1);
281 std::vector<ftm::idNode> children;
283 for(
int const cIdx : children) {
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);
294 while(!stack.empty()) {
295 int const nIdx = stack.top();
297 preorder2[count] = nIdx;
299 depth2 = std::max((
int)predecessors2[nIdx].size(), depth2);
300 std::vector<ftm::idNode> children;
302 for(
int const cIdx : children) {
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);
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;
321 std::vector<dataType> memT((nn1 + 1) * (depth1 + 1) * (nn2 + 1)
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;
329 for(
size_t l = 1; l <= predecessors1[preorder1[i]].size(); l++) {
330 int parent1 = predecessors1[preorder1[i]]
331 [predecessors1[preorder1[i]].size() - l];
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);
351 dataType c = std::numeric_limits<dataType>::max();
352 for(
auto child1_mb : children1) {
354 = memT[child1_mb + (l + 1) * dim2 + nn2 * dim3 + 0 * dim4];
355 for(
auto child1 : children1) {
356 if(child1 == child1_mb) {
359 c_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
363 memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4] = c;
367 for(
size_t j = 0; j < nn2; j++) {
368 int curr2 = preorder2[j];
369 std::vector<ftm::idNode> 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];
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);
393 dataType c = std::numeric_limits<dataType>::max();
394 for(
auto child2_mb : children2) {
396 = memT[nn1 + 0 * dim2 + child2_mb * dim3 + (l + 1) * dim4];
397 for(
auto child2 : children2) {
398 if(child2 == child2_mb) {
401 c_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
405 memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4] = c;
410 for(
size_t i = 0; i < nn1; i++) {
411 int curr1 = preorder1[i];
412 std::vector<ftm::idNode> children1;
414 for(
size_t j = 0; j < nn2; j++) {
415 int curr2 = preorder2[j];
416 std::vector<ftm::idNode> children2;
418 for(
size_t l1 = 1; l1 <= predecessors1[preorder1[i]].size(); l1++) {
420 = predecessors1[preorder1[i]]
421 [predecessors1[preorder1[i]].size() - l1];
422 for(
size_t l2 = 1; l2 <= predecessors2[preorder2[j]].size(); l2++) {
424 = predecessors2[preorder2[j]]
425 [predecessors2[preorder2[j]].size() - l2];
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);
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
454 for(
auto child2 : children2) {
455 if(child2 == child2_mb) {
458 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
462 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
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
472 for(
auto child1 : children1) {
473 if(child1 == child1_mb) {
476 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
480 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
486 dataType d = std::numeric_limits<dataType>::max();
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>(
498 memT[child11 + (l1 + 1) * dim2 + child21 * dim3
500 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]);
501 d = std::min<dataType>(
503 memT[child12 + (l1 + 1) * dim2 + child22 * dim3
505 + memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]);
506 d = std::min<dataType>(
508 memT[child11 + (l1 + 1) * dim2 + child22 * dim3
510 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]);
511 d = std::min<dataType>(
513 memT[child12 + (l1 + 1) * dim2 + child21 * dim3
515 + memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]);
517 for(
auto child1_mb : children1) {
518 std::vector<ftm::idNode> topo1_;
521 std::remove(topo1_.begin(), topo1_.end(), child1_mb),
523 for(
auto child2_mb : children2) {
524 std::vector<ftm::idNode> topo2_;
527 std::remove(topo2_.begin(), topo2_.end(), child2_mb),
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];
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);
549 switch(assignmentSolverID_) {
552 assignmentSolver = &solverExhaustive;
556 assignmentSolver = &solverMunkres;
561 assignmentSolver = &solverAuction;
563 assignmentSolver->
setInput(costMatrix);
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);
578 for(
auto child1_mb : children1) {
579 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3
581 for(
auto child1 : children1) {
582 if(child1 == child1_mb) {
585 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
593 for(
auto child2_mb : children2) {
594 dataType d_ = memT[curr1 + l1 * dim2 + child2_mb * dim3
596 for(
auto child2 : children2) {
597 if(child2 == child2_mb) {
600 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
604 memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] = d;
611 std::vector<ftm::idNode> children1;
613 std::vector<ftm::idNode> children2;
617 = memT[children1[0] + 1 * dim2 + children2[0] * dim3 + 1 * dim4];
619 if(computeMapping_ && outputMatching) {
621 outputMatching->clear();
624 std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>>
629 predecessors1, predecessors2, depth1, depth2, memT,
632 for(
auto m : mapping) {
633 if(writeOptimalBranchDecomposition_ && m.first.first >= 0
634 && m.first.second >= 0) {
637 linkedNodes1[m.first.first] = m.first.second;
638 linkedNodes1[m.first.second] = m.first.first;
640 if(writeOptimalBranchDecomposition_ && m.second.first >= 0
641 && m.second.second >= 0) {
644 linkedNodes2[m.second.first] = m.second.second;
645 linkedNodes2[m.second.second] = m.second.first;
647 if(m.first.first == -1)
649 if(m.first.second == -1)
651 if(m.second.first == -1)
653 if(m.second.second == -1)
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,
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,
667 m.second.second, tree1, tree2)
668 : editCost_Shifting<dataType>(m.first.first, m.first.second,
669 m.second.first, m.second.second,
671 matchedCost[m.first.second] = matchedCost[m.first.first];
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]));
680 return squared_ ? std::sqrt(res) : res;
691 std::vector<std::vector<int>> &predecessors1,
692 std::vector<std::vector<int>> &predecessors2,
695 std::vector<dataType> &memT,
696 std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>>
702 int dim2 = (nn1 + 1) * dim1;
703 int dim3 = (depth1 + 1) * dim2;
704 int dim4 = (nn2 + 1) * dim3;
710 std::vector<ftm::idNode> children1;
712 int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
717 mapping.emplace_back(std::make_pair(
718 std::make_pair(curr1, parent1), std::make_pair(-1, -1)));
724 for(
auto child1_mb : children1) {
726 = memT[child1_mb + (l1 + 1) * dim2 + nn2 * dim3 + 0 * dim4];
727 for(
auto child1 : children1) {
728 if(child1 == child1_mb) {
731 c_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
733 if(c_ == memT[curr1 + l1 * dim2 + nn2 * dim3 + 0 * dim4]) {
735 predecessors1, predecessors2, depth1, depth2,
737 for(
auto child1 : children1) {
738 if(child1 == child1_mb) {
742 predecessors1, predecessors2, depth1,
743 depth2, memT, mapping);
748 this->
printErr(
"Mapping traceback not correct.");
756 std::vector<ftm::idNode> children2;
758 int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
763 mapping.emplace_back(std::make_pair(
764 std::make_pair(-1, -1), std::make_pair(curr2, parent2)));
770 for(
auto child2_mb : children2) {
772 = memT[nn1 + 0 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
773 for(
auto child2 : children2) {
774 if(child2 == child2_mb) {
777 c_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
779 if(c_ == memT[nn1 + 0 * dim2 + curr2 * dim3 + l2 * dim4]) {
781 predecessors1, predecessors2, depth1, depth2,
783 for(
auto child2 : children2) {
784 if(child2 == child2_mb) {
788 predecessors1, predecessors2, depth1,
789 depth2, memT, mapping);
794 this->
printErr(
"Mapping traceback not correct.");
798 std::vector<ftm::idNode> children1;
800 std::vector<ftm::idNode> children2;
802 int parent1 = predecessors1[curr1][predecessors1[curr1].size() - l1];
803 int parent2 = predecessors2[curr2][predecessors2[curr2].size() - l2];
813 mapping.emplace_back(std::make_pair(
814 std::make_pair(curr1, parent1), std::make_pair(curr2, parent2)));
820 else if(children1.size() == 0) {
821 for(
auto child2_mb : children2) {
823 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
824 for(
auto child2 : children2) {
825 if(child2 == child2_mb) {
828 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
830 if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
832 predecessors1, predecessors2, depth1, depth2,
834 for(
auto child2 : children2) {
835 if(child2 == child2_mb) {
839 predecessors1, predecessors2, depth1, depth2,
849 else if(children2.size() == 0) {
850 for(
auto child1_mb : children1) {
852 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
853 for(
auto child1 : children1) {
854 if(child1 == child1_mb) {
857 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
859 if(d_ == memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]) {
861 predecessors1, predecessors2, depth1, depth2,
863 for(
auto child1 : children1) {
864 if(child1 == child1_mb) {
868 predecessors1, predecessors2, depth1, depth2,
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
891 + memT[child12 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
894 (l2 + 1), predecessors1, predecessors2, depth1,
895 depth2, memT, mapping);
897 predecessors1, predecessors2, depth1, depth2,
902 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
903 == memT[child12 + (l1 + 1) * dim2 + child22 * dim3
905 + memT[child11 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
908 (l2 + 1), predecessors1, predecessors2, depth1,
909 depth2, memT, mapping);
911 predecessors1, predecessors2, depth1, depth2,
916 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
917 == memT[child11 + (l1 + 1) * dim2 + child22 * dim3
919 + memT[child12 + 1 * dim2 + child21 * dim3 + 1 * dim4]) {
922 (l2 + 1), predecessors1, predecessors2, depth1,
923 depth2, memT, mapping);
925 predecessors1, predecessors2, depth1, depth2,
930 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4]
931 == memT[child12 + (l1 + 1) * dim2 + child21 * dim3
933 + memT[child11 + 1 * dim2 + child22 * dim3 + 1 * dim4]) {
936 (l2 + 1), predecessors1, predecessors2, depth1,
937 depth2, memT, mapping);
939 predecessors1, predecessors2, depth1, depth2,
945 for(
auto child1_mb : children1) {
946 std::vector<ftm::idNode> topo1_;
948 topo1_.erase(std::remove(topo1_.begin(), topo1_.end(), child1_mb),
950 for(
auto child2_mb : children2) {
951 std::vector<ftm::idNode> topo2_;
953 topo2_.erase(std::remove(topo2_.begin(), topo2_.end(), child2_mb),
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];
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);
975 switch(assignmentSolverID_) {
978 assignmentSolver = &solverExhaustive;
982 assignmentSolver = &solverMunkres;
987 assignmentSolver = &solverAuction;
989 assignmentSolver->
setInput(costMatrix);
991 assignmentSolver->
run(matching);
992 dataType d_ = memT[child1_mb + (l1 + 1) * dim2 + child2_mb * dim3
994 for(
auto m : matching)
995 d_ += std::get<2>(m);
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)]
1005 int n2 = std::get<1>(m) <
static_cast<int>(topo2_.size())
1006 ? topo2_[std::get<1>(m)]
1008 if(n1 >= 0 && n2 >= 0)
1010 predecessors1, predecessors2, depth1,
1011 depth2, memT, mapping);
1014 predecessors1, predecessors2, depth1,
1015 depth2, memT, mapping);
1018 predecessors1, predecessors2, depth1,
1019 depth2, memT, mapping);
1030 for(
auto child1_mb : children1) {
1032 = memT[child1_mb + (l1 + 1) * dim2 + curr2 * dim3 + l2 * dim4];
1033 for(
auto child1 : children1) {
1034 if(child1 == child1_mb) {
1037 d_ += memT[child1 + 1 * dim2 + nn2 * dim3 + 0 * dim4];
1039 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
1041 predecessors1, predecessors2, depth1, depth2,
1043 for(
auto child1 : children1) {
1044 if(child1 == child1_mb) {
1048 predecessors1, predecessors2, depth1, depth2,
1058 for(
auto child2_mb : children2) {
1060 = memT[curr1 + l1 * dim2 + child2_mb * dim3 + (l2 + 1) * dim4];
1061 for(
auto child2 : children2) {
1062 if(child2 == child2_mb) {
1065 d_ += memT[nn1 + 0 * dim2 + child2 * dim3 + 1 * dim4];
1067 if(memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] == d_) {
1069 predecessors1, predecessors2, depth1, depth2,
1071 for(
auto child2 : children2) {
1072 if(child2 == child2_mb) {
1076 predecessors1, predecessors2, depth1, depth2,
1082 this->
printErr(
"Mapping traceback not correct");