From 4fe495be8dd0233fe961271cd8896e93e7fa2e0c Mon Sep 17 00:00:00 2001 From: Oliver Ruebel Date: Tue, 1 Mar 2022 16:24:37 -0700 Subject: [PATCH] Implement branch decomposition for hierarchical contour tree Co-authored-by: Gunther H. Weber --- .../BranchCompilerApp.cxx | 131 ++++ .../contour_tree_distributed/CMakeLists.txt | 6 + .../ContourTreeApp.cxx | 116 ++- .../hact_test_branch_decomposition.sh | 44 ++ .../testrun_branch_decomposition.sh | 99 +++ vtkm/filter/ContourTreeUniformDistributed.h | 16 +- vtkm/filter/ContourTreeUniformDistributed.hxx | 210 +++++- ...stingContourTreeUniformDistributedFilter.h | 5 +- .../BranchDecompositionBlock.h | 107 +++ .../contourtree_distributed/CMakeLists.txt | 4 + .../CombineHyperSweepBlockFunctor.h | 8 +- ...uteDistributedBranchDecompositionFunctor.h | 238 +++++++ .../ComputeDistributedContourTreeFunctor.h | 4 +- .../DistributedContourTreeBlockData.h | 2 +- .../HierarchicalContourTree.h | 8 + .../HierarchicalVolumetricBranchDecomposer.h | 658 ++++++++++++++++++ .../contourtree_distributed/HyperSweepBlock.h | 4 +- .../hierarchical_contour_tree/CMakeLists.txt | 1 + .../FindSuperArcBetweenNodes.h | 142 ++++ .../CMakeLists.txt | 21 + .../CollapseBranchesPointerDoublingWorklet.h | 94 +++ .../CollapseBranchesWorklet.h | 244 +++++++ .../FindBestSupernodeWorklet.h | 97 +++ ...lBestUpDownByVolumeBestUpDownEdgeWorklet.h | 154 ++++ ...estUpDownByVolumeInitSuperarcListWorklet.h | 104 +++ .../LocalBestUpDownByVolumeWorklet.h | 229 ++++++ ...tricComparatorIndirectGlobalIdComparator.h | 247 +++++++ 27 files changed, 2918 insertions(+), 75 deletions(-) create mode 100644 examples/contour_tree_distributed/BranchCompilerApp.cxx create mode 100755 examples/contour_tree_distributed/hact_test_branch_decomposition.sh create mode 100755 examples/contour_tree_distributed/testrun_branch_decomposition.sh create mode 100644 vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h create mode 100644 vtkm/worklet/contourtree_distributed/ComputeDistributedBranchDecompositionFunctor.h create mode 100644 vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcBetweenNodes.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CMakeLists.txt create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesPointerDoublingWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/FindBestSupernodeWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeInitSuperarcListWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeWorklet.h create mode 100644 vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/SuperArcVolumetricComparatorIndirectGlobalIdComparator.h diff --git a/examples/contour_tree_distributed/BranchCompilerApp.cxx b/examples/contour_tree_distributed/BranchCompilerApp.cxx new file mode 100644 index 000000000..a9cd1afa3 --- /dev/null +++ b/examples/contour_tree_distributed/BranchCompilerApp.cxx @@ -0,0 +1,131 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// +// COMMENTS: +// +// Input is assumed to be a sequence of lines of the form: +// I Global ID of branch root +// II Value of supernode +// III Global ID of supernode +// +// All lines are assumed to have been sorted already. Because of how the +// Unix sort utility operates (textual sort), the most we can assume is that all +// supernodes corresponding to a given branch root are sorted together. +// +// We therefore do simple stream processing, identifying new branches by +// the changes in root ID. +// +//======================================================================================= + +#include +#include + +int main() +{ // main() + // variables tracking the best & worst so far for this extent + vtkm::Id currentBranch = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + // this is slightly tricky, since we don't know the range of the data type + // yet, but we can initialize to 0 for both floats and integers, then test on + // current branch + vtkm::Float64 highValue = 0; + vtkm::Float64 lowValue = 0; + vtkm::Id highEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + vtkm::Id lowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + + // values to read in + vtkm::Id nextBranch; + vtkm::Float64 nextValue; + vtkm::Id nextSupernode; + + while (true) + { // until stream goes bad + // read the next triple + std::cin >> nextBranch >> nextValue >> nextSupernode; + + // test in the middle before processing + if (std::cin.eof()) + break; + + // test to see if the branch is different from the current one + if (nextBranch != currentBranch) + { // new branch + // special test for initial one + if (!vtkm::worklet::contourtree_augmented::NoSuchElement(currentBranch)) + printf("%12llu %12llu\n", highEnd, lowEnd); + + // set the high & low value ends to this one + highValue = nextValue; + lowValue = nextValue; + highEnd = nextSupernode; + lowEnd = nextSupernode; + + // and reset the branch ID + currentBranch = nextBranch; + } // new branch + else + { // existing branch + // test value with simulation of simplicity + if ((nextValue > highValue) || ((nextValue == highValue) && (nextSupernode > highEnd))) + { // new high end + highEnd = nextSupernode; + highValue = nextValue; + } // new high end + // test value with simulation of simplicity + else if ((nextValue < lowValue) || ((nextValue == lowValue) && (nextSupernode < lowEnd))) + { // new low end + lowEnd = nextSupernode; + lowValue = nextValue; + } // new low end + } // existing branch + } // until stream goes bad + + printf("%12llu %12llu\n", highEnd, lowEnd); +} // main() diff --git a/examples/contour_tree_distributed/CMakeLists.txt b/examples/contour_tree_distributed/CMakeLists.txt index 75f3012a4..6c146ac3a 100644 --- a/examples/contour_tree_distributed/CMakeLists.txt +++ b/examples/contour_tree_distributed/CMakeLists.txt @@ -81,10 +81,16 @@ if (VTKm_ENABLE_MPI) target_link_libraries(TreeCompiler vtkm_filter) vtkm_add_target_information(TreeCompiler DROP_UNUSED_SYMBOLS) + add_executable(BranchCompiler BranchCompilerApp.cxx) + target_link_libraries(BranchCompiler vtkm_filter) + vtkm_add_target_information(BranchCompiler DROP_UNUSED_SYMBOLS) + configure_file(split_data_2d.py split_data_2d.py COPYONLY) configure_file(split_data_3d.py split_data_3d.py COPYONLY) configure_file(hact_test.sh hact_test.sh COPYONLY) configure_file(hact_test_volume.sh hact_test_volume.sh COPYONLY) + configure_file(hact_test_branch_decomposition.sh hact_test_branch_decomposition.sh COPYONLY) + configure_file(testrun_branch_decomposition.sh testrun_branch_decomposition.sh COPYONLY) configure_file(testrun.sh testrun.sh COPYONLY) configure_file(testrun_volume.sh testrun_volume.sh COPYONLY) endif() diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index 0278356bf..0ea3ae8e4 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -76,6 +76,7 @@ #include #include #include +#include #include // clang-format off @@ -203,6 +204,19 @@ int main(int argc, char* argv[]) augmentHierarchicalTree = true; } + bool computeHierarchicalVolumetricBranchDecomposition = false; + if (parser.hasOption("--computeVolumeBranchDecomposition")) + { + computeHierarchicalVolumetricBranchDecomposition = true; + if (!augmentHierarchicalTree) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Warn, + "Warning: --computeVolumeBranchDecomposition only " + "allowed augmentation. Enabling --augmentHierarchicalTree option."); + augmentHierarchicalTree = true; + } + } + bool useBoundaryExtremaOnly = true; if (parser.hasOption("--useFullBoundary")) { @@ -261,7 +275,7 @@ int main(int argc, char* argv[]) { if (rank == 0) { - std::cout << "ContourTreeAugmented " << std::endl; + std::cout << "ContourTreeDistributed " << std::endl; std::cout << std::endl; std::cout << " Name of the input data file." << std::endl; std::cout << "The file is expected to be ASCII with either: " << std::endl; @@ -284,6 +298,9 @@ int main(int argc, char* argv[]) << std::endl; std::cout << " and when using only boundary extrema." << std::endl; std::cout << "--augmentHierarchicalTree Augment the hierarchical tree." << std::endl; + std::cout << "--computeVolumeBranchDecomposition Compute the volume branch decomposition. " + << std::endl; + std::cout << " Requries --augmentHierarchicalTree to be set." << std::endl; std::cout << "--preSplitFiles Input data is already pre-split into blocks." << std::endl; std::cout << "--saveDot Save DOT files of the distributed contour tree " << std::endl << " computation (Default=False). " << std::endl; @@ -311,6 +328,9 @@ int main(int argc, char* argv[]) << " mc=" << useMarchingCubes << std::endl << " useFullBoundary=" << !useBoundaryExtremaOnly << std::endl << " saveDot=" << saveDotFiles << std::endl + << " augmentHierarchicalTree=" << augmentHierarchicalTree << std::endl + << " computeVolumetricBranchDecomposition=" + << computeHierarchicalVolumetricBranchDecomposition << std::endl << " saveOutputData=" << saveOutputData << std::endl << " forwardSummary=" << forwardSummary << std::endl << " nblocks=" << numBlocks << std::endl); @@ -842,17 +862,19 @@ int main(int argc, char* argv[]) prevTime = currTime; // Convert the mesh of values into contour tree, pairs of vertex ids - vtkm::filter::ContourTreeUniformDistributed filter(blocksPerDim, - globalSize, - localBlockIndices, - localBlockOrigins, - localBlockSizes, - useBoundaryExtremaOnly, - useMarchingCubes, - augmentHierarchicalTree, - saveDotFiles, - timingsLogLevel, - treeLogLevel); + vtkm::filter::ContourTreeUniformDistributed filter( + blocksPerDim, + globalSize, + localBlockIndices, + localBlockOrigins, + localBlockSizes, + useBoundaryExtremaOnly, + useMarchingCubes, + augmentHierarchicalTree, + computeHierarchicalVolumetricBranchDecomposition, + saveDotFiles, + timingsLogLevel, + treeLogLevel); filter.SetActiveField("values"); // Execute the contour tree analysis @@ -872,35 +894,53 @@ int main(int argc, char* argv[]) { if (augmentHierarchicalTree) { - for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) + if (computeHierarchicalVolumetricBranchDecomposition) { - auto ds = result.GetPartition(ds_no); - vtkm::worklet::contourtree_augmented::IdArrayType supernodes; - ds.GetField("Supernodes").GetData().AsArrayHandle(supernodes); - vtkm::worklet::contourtree_augmented::IdArrayType superarcs; - ds.GetField("Superarcs").GetData().AsArrayHandle(superarcs); - vtkm::worklet::contourtree_augmented::IdArrayType regularNodeGlobalIds; - ds.GetField("RegularNodeGlobalIds").GetData().AsArrayHandle(regularNodeGlobalIds); - vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; - vtkm::worklet::contourtree_augmented::IdArrayType intrinsicVolume; - ds.GetField("IntrinsicVolume").GetData().AsArrayHandle(intrinsicVolume); - vtkm::worklet::contourtree_augmented::IdArrayType dependentVolume; - ds.GetField("DependentVolume").GetData().AsArrayHandle(dependentVolume); + for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) + { + auto ds = result.GetPartition(ds_no); + std::string branchDecompositionFileName = std::string("BranchDecomposition_Rank_") + + std::to_string(static_cast(rank)) + std::string("_Block_") + + std::to_string(static_cast(ds_no)) + std::string(".txt"); - std::string dumpVolumesString = - vtkm::worklet::contourtree_distributed::HierarchicalContourTree::DumpVolumes( - supernodes, - superarcs, - regularNodeGlobalIds, - totalVolume, - intrinsicVolume, - dependentVolume); + std::ofstream treeStream(branchDecompositionFileName.c_str()); + treeStream + << vtkm::worklet::contourtree_distributed::HierarchicalVolumetricBranchDecomposer< + ValueType>::PrintBranches(ds); + } + } + else + { + for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) + { + auto ds = result.GetPartition(ds_no); + vtkm::worklet::contourtree_augmented::IdArrayType supernodes; + ds.GetField("Supernodes").GetData().AsArrayHandle(supernodes); + vtkm::worklet::contourtree_augmented::IdArrayType superarcs; + ds.GetField("Superarcs").GetData().AsArrayHandle(superarcs); + vtkm::worklet::contourtree_augmented::IdArrayType regularNodeGlobalIds; + ds.GetField("RegularNodeGlobalIds").GetData().AsArrayHandle(regularNodeGlobalIds); + vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; + vtkm::worklet::contourtree_augmented::IdArrayType intrinsicVolume; + ds.GetField("IntrinsicVolume").GetData().AsArrayHandle(intrinsicVolume); + vtkm::worklet::contourtree_augmented::IdArrayType dependentVolume; + ds.GetField("DependentVolume").GetData().AsArrayHandle(dependentVolume); - std::string volumesFileName = std::string("TreeWithVolumes_Rank_") + - std::to_string(static_cast(rank)) + std::string("_Block_") + - std::to_string(static_cast(ds_no)) + std::string(".txt"); - std::ofstream treeStream(volumesFileName.c_str()); - treeStream << dumpVolumesString; + std::string dumpVolumesString = + vtkm::worklet::contourtree_distributed::HierarchicalContourTree::DumpVolumes( + supernodes, + superarcs, + regularNodeGlobalIds, + totalVolume, + intrinsicVolume, + dependentVolume); + + std::string volumesFileName = std::string("TreeWithVolumes_Rank_") + + std::to_string(static_cast(rank)) + std::string("_Block_") + + std::to_string(static_cast(ds_no)) + std::string(".txt"); + std::ofstream treeStream(volumesFileName.c_str()); + treeStream << dumpVolumesString; + } } } else diff --git a/examples/contour_tree_distributed/hact_test_branch_decomposition.sh b/examples/contour_tree_distributed/hact_test_branch_decomposition.sh new file mode 100755 index 000000000..34737f97c --- /dev/null +++ b/examples/contour_tree_distributed/hact_test_branch_decomposition.sh @@ -0,0 +1,44 @@ +#!/bin/sh +GTCT_DIR=${GTCT_DIR:-${HOME}/devel/parallel-peak-pruning/ContourTree/SweepAndMergeSerial/out} +RED="" +GREEN="" +NC="" +if [ -t 1 ]; then +# If stdout is a terminal, color Pass and FAIL green and red, respectively + RED=$(tput setaf 1) + GREEN=$(tput setaf 2) + NC=$(tput sgr0) +fi + +echo "Removing previously generated files" +rm *.log *.dat + +echo "Copying target file "$1 "into current directory" +filename=${1##*/} +fileroot=${filename%.txt} + +cp $1 ${filename} + +echo "Splitting data into "$2" x "$2" parts" +./split_data_2d.py ${filename} $2 +rm ${filename} + +echo "Running HACT" +n_parts=$(($2*$2)) +# mpirun -np 4 --oversubscribe ./ContourTree_Distributed --vtkm-device Any --preSplitFiles --saveOutputData --augmentHierarchicalTree --computeVolumeBranchDecomposition --numBlocks=${n_parts} ${fileroot}_part_%d_of_${n_parts}.txt +mpirun -np 2 --oversubscribe ./ContourTree_Distributed --vtkm-device Any --preSplitFiles --saveOutputData --augmentHierarchicalTree --computeVolumeBranchDecomposition --numBlocks=${n_parts} ${fileroot}_part_%d_of_${n_parts}.txt +rm ${fileroot}_part_*_of_${n_parts}.txt + +echo "Compiling Outputs" +sort -u BranchDecomposition_Rank_*.txt > outsort${fileroot}_$2x$2.txt +cat outsort${fileroot}_$2x$2.txt | ./BranchCompiler | sort > bcompile${fileroot}_$2x$2.txt +rm BranchDecomposition_Rank_*.txt outsort${fileroot}_$2x$2.txt + +echo "Diffing" +echo diff bcompile${fileroot}_$2x$2.txt ${GTCT_DIR}/branch_decomposition_volume_${fileroot}.txt +diff bcompile${fileroot}_$2x$2.txt ${GTCT_DIR}/branch_decomposition_volume_${fileroot}.txt + +if test $? -eq 0; then echo "${GREEN}Pass${NC}"; rm bcompile${fileroot}_$2x$2.txt; else echo "${RED}FAIL${NC}"; fi; + +# echo "Generating Dot files" +# ./makedot.sh diff --git a/examples/contour_tree_distributed/testrun_branch_decomposition.sh b/examples/contour_tree_distributed/testrun_branch_decomposition.sh new file mode 100755 index 000000000..5de23883a --- /dev/null +++ b/examples/contour_tree_distributed/testrun_branch_decomposition.sh @@ -0,0 +1,99 @@ +#!/bin/sh + +mkdir -p out +DATA_DIR=${DATA_DIR:-${HOME}/devel/parallel-peak-pruning/Data/2D} + +if [ ! -d $DATA_DIR ]; then + echo "Error: Directory $DATA_DIR does not exist!" + exit 1; +fi; + +echo +echo "Starting Timing Runs" +echo +echo "8x9 Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/8x9test.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/8x9test.txt 4 +# ./hact_test_branch_decomposition.sh $DATA_DIR/8x9test.txt 8 +echo +echo "Vancouver Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vanc.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vanc.txt 4 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vanc.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vanc.txt 16 +echo +echo "Vancouver SWSW Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSW.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSW.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSW.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSW.txt 16 +echo +echo "Vancouver SWNW Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNW.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNW.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNW.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNW.txt 16 +echo +echo "Vancouver SWSE Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSE.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSE.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSE.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWSE.txt 16 +echo +echo "Vancouver SWNE Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNE.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNE.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNE.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSWNE.txt 16 +echo +echo "Vancouver NE Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNE.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNE.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNE.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNE.txt 16 +echo +echo "Vancouver NW Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNW.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNW.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNW.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverNW.txt 16 +echo +echo "Vancouver SE Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSE.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSE.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSE.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSE.txt 16 +echo +echo "Vancouver SW Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSW.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSW.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSW.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/vancouverSW.txt 16 +echo +echo "Icefields Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 16 +# ./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 32 +# ./hact_test_branch_decomposition.sh $DATA_DIR/icefield.txt 64 +echo +echo "GTOPO30 Full Tiny Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 16 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 32 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo_full_tiny.txt 64 +echo +echo "GTOPO30 UK Tile Test Set" +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 2 +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 4 +./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 8 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 16 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 32 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 64 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 128 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 256 +# ./hact_test_branch_decomposition.sh $DATA_DIR/gtopo30w020n40.txt 512 +echo "Done" diff --git a/vtkm/filter/ContourTreeUniformDistributed.h b/vtkm/filter/ContourTreeUniformDistributed.h index 595db046e..caf4befa4 100644 --- a/vtkm/filter/ContourTreeUniformDistributed.h +++ b/vtkm/filter/ContourTreeUniformDistributed.h @@ -119,6 +119,7 @@ public: bool useBoundaryExtremaOnly = true, bool useMarchingCubes = false, bool augmentHierarchicalTree = false, + bool computeHierarchicalVolumetricBranchDecomposition = false, bool saveDotFiles = false, vtkm::cont::LogLevel timingsLogLevel = vtkm::cont::LogLevel::Perf, vtkm::cont::LogLevel treeLogLevel = vtkm::cont::LogLevel::Info); @@ -165,12 +166,20 @@ public: template VTKM_CONT void ComputeVolumeMetric( - vtkmdiy::Master& inputContoourTreeMaster, + vtkmdiy::Master& inputContourTreeMaster, vtkmdiy::DynamicAssigner& assigner, vtkmdiy::RegularSwapPartners& partners, const FieldType&, // dummy parameter to get the type std::stringstream& timingsStream, - std::vector hierarchicalTreeOutputDataSet); + std::vector& hierarchicalTreeOutputDataSet); + + template + VTKM_CONT void ComputeBranchDecomposition(vtkmdiy::Master& inputContourTreeMaster, + vtkmdiy::DynamicAssigner& assigner, + vtkmdiy::RegularSwapPartners& partners, + const FieldType&, // dummy parameter to get the type + std::stringstream& timingsStream, + std::vector& outputDataSet); /// /// Internal helper function that implements the actual functionality of PostExecute @@ -196,6 +205,9 @@ private: /// Augment hierarchical tree bool AugmentHierarchicalTree; + /// Compute the hierarchical volumetric branch decomposition + bool ComputeHierarchicalVolumetricBranchDecomposition; + /// Save dot files for all tree computations bool SaveDotFiles; diff --git a/vtkm/filter/ContourTreeUniformDistributed.hxx b/vtkm/filter/ContourTreeUniformDistributed.hxx index 666342150..afe7ae2a5 100644 --- a/vtkm/filter/ContourTreeUniformDistributed.hxx +++ b/vtkm/filter/ContourTreeUniformDistributed.hxx @@ -66,12 +66,15 @@ // distributed contour tree includes #include #include +#include #include +#include #include #include #include #include #include +#include #include #include #include @@ -208,6 +211,7 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed( bool useBoundaryExtremaOnly, bool useMarchingCubes, bool augmentHierarchicalTree, + bool computeHierarchicalVolumetricBranchDecomposition, bool saveDotFiles, vtkm::cont::LogLevel timingsLogLevel, vtkm::cont::LogLevel treeLogLevel) @@ -215,6 +219,8 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed( , UseBoundaryExtremaOnly(useBoundaryExtremaOnly) , UseMarchingCubes(useMarchingCubes) , AugmentHierarchicalTree(augmentHierarchicalTree) + , ComputeHierarchicalVolumetricBranchDecomposition( + computeHierarchicalVolumetricBranchDecomposition) , SaveDotFiles(saveDotFiles) , TimingsLogLevel(timingsLogLevel) , TreeLogLevel(treeLogLevel) @@ -229,6 +235,13 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed( , LocalInteriorForests(static_cast(localBlockSizes.GetNumberOfValues())) { this->SetOutputFieldName("resultData"); + if (this->ComputeHierarchicalVolumetricBranchDecomposition && !this->AugmentHierarchicalTree) + { + this->AugmentHierarchicalTree = true; + VTKM_LOG_S(vtkm::cont::LogLevel::Error, + "Volumetric branch decmomposition requires augmentation. " + << "Enabling the AugmentHierarchicalTree option."); + } } @@ -626,12 +639,12 @@ inline VTKM_CONT void ContourTreeUniformDistributed::PostExecute( template inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( - vtkmdiy::Master& inputContoourTreeMaster, + vtkmdiy::Master& inputContourTreeMaster, vtkmdiy::DynamicAssigner& assigner, vtkmdiy::RegularSwapPartners& partners, const FieldType&, // dummy parameter to get the type std::stringstream& timingsStream, - std::vector hierarchicalTreeOutputDataSet) + std::vector& hierarchicalTreeOutputDataSet) { // TODO/FIXME: CONSIDER MOVING CONTENTS OF THIS METHOD TO SEPARATE FILTER vtkm::cont::Timer timer; @@ -644,7 +657,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( 1, // Use 1 thread, VTK-M will do the treading -1, // All blocks in memory 0, // No create function - HyperSweepBlock::destroy); + HyperSweepBlock::Destroy); // Log the time to create the DIY master for the hyper sweep timingsStream << " " << std::setw(38) << std::left << "Create DIY Master (Hypersweep)" @@ -654,7 +667,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( // Copy data from hierarchical tree computation to initialize volume computation using DistributedContourTreeBlockData = vtkm::worklet::contourtree_distributed::DistributedContourTreeBlockData; - inputContoourTreeMaster.foreach ( + inputContourTreeMaster.foreach ( [&](DistributedContourTreeBlockData* currInBlock, const vtkmdiy::Master::ProxyWithLink&) { vtkm::Id blockNo = currInBlock->LocalBlockNo; // The block size and origin may be modified during the FanIn so we need to use the @@ -667,7 +680,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( // NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all // communication is via RegularDecomposer, which sets up its own links. No need // to keep the pointer, as DIY will "own" it and delete it when no longer needed. - // NOTE: Since we passed a "destroy" function to DIY master, it will own the local data + // NOTE: Since we passed a "Destroy" function to DIY master, it will own the local data // blocks and delete them when done. hierarchical_hyper_sweep_master.add( currInBlock->GlobalBlockId, @@ -680,13 +693,18 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( new vtkmdiy::Link()); }); - vtkmdiy::fix_links(hierarchical_hyper_sweep_master, assigner); - // Log time to copy the data to the HyperSweepBlock data objects timingsStream << " " << std::setw(38) << std::left << "Initialize Hypersweep Data" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); + vtkmdiy::fix_links(hierarchical_hyper_sweep_master, assigner); + + // Record time to fix the links + timingsStream << " " << std::setw(38) << std::left << "Fix DIY Links (Hypersweep)" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + hierarchical_hyper_sweep_master.foreach ([&](HyperSweepBlock* b, const vtkmdiy::Master::ProxyWithLink&) { std::stringstream localHypersweepTimingsStream; @@ -695,8 +713,8 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( // Create HyperSweeper #ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block " << b->GlobalBlockId); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, + VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); + VTKM_LOG_S(this->TreeLogLevel, b->HierarchicalContourTree.DebugPrint( "Before initializing HierarchicalHyperSweeper", __FILE__, __LINE__)); #endif @@ -731,8 +749,8 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( } #ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block " << b->GlobalBlockId); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, + VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); + VTKM_LOG_S(this->TreeLogLevel, b->HierarchicalContourTree.DebugPrint( "After initializing intrinsic vertex count", __FILE__, __LINE__)); std::ostringstream volumeStream; @@ -742,8 +760,8 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( "Intrinsic Volume", b->IntrinsicVolume, -1, volumeStream); vtkm::worklet::contourtree_augmented::PrintIndices( "Dependent Volume", b->DependentVolume, -1, volumeStream); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, volumeStream.str()); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "FLUSH" << std::endl << std::flush); + VTKM_LOG_S(this->TreeLogLevel, volumeStream.str()); + VTKM_LOG_S(this->TreeLogLevel, "FLUSH" << std::endl << std::flush); #endif // Initialize dependentVolume by copy from intrinsicVolume vtkm::cont::Algorithm::Copy(b->IntrinsicVolume, b->DependentVolume); @@ -758,8 +776,8 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( hyperSweeper.LocalHyperSweep(); #ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block " << b->GlobalBlockId); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, + VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); + VTKM_LOG_S(this->TreeLogLevel, b->HierarchicalContourTree.DebugPrint("After local hypersweep", __FILE__, __LINE__)); #endif // Log the local hypersweep time @@ -805,9 +823,9 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( "DependentVolume", vtkm::cont::Field::Association::WholeMesh, b->DependentVolume); hierarchicalTreeOutputDataSet[b->LocalBlockNo].AddField(dependentVolumeField); #ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block " << b->GlobalBlockId); + VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); VTKM_LOG_S( - vtkm::cont::LogLevel::Info, + this->TreeLogLevel, b->HierarchicalContourTree.DebugPrint("Called from DumpVolumes", __FILE__, __LINE__)); std::ostringstream volumeStream; vtkm::worklet::contourtree_augmented::PrintHeader(b->IntrinsicVolume.GetNumberOfValues(), @@ -816,15 +834,151 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( "Intrinsic Volume", b->IntrinsicVolume, -1, volumeStream); vtkm::worklet::contourtree_augmented::PrintIndices( "Dependent Volume", b->DependentVolume, -1, volumeStream); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, volumeStream.str()); + VTKM_LOG_S(this->TreeLogLevel, volumeStream.str()); #endif // Log the time for adding hypersweep data to the output dataset timingsStream << " " << std::setw(38) << std::left << "Create Output Data (Hypersweep)" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; - timer.Start(); }); } +//----------------------------------------------------------------------------- +template +inline VTKM_CONT void ContourTreeUniformDistributed::ComputeBranchDecomposition( + vtkmdiy::Master& inputContourTreeMaster, + vtkmdiy::DynamicAssigner& assigner, + vtkmdiy::RegularSwapPartners& partners, + const FieldType&, // dummy parameter to get the type + std::stringstream& timingsStream, + std::vector& outputDataSet) +{ + vtkm::cont::Timer timer; + timer.Start(); + + using BranchDecompositionBlock = + vtkm::worklet::contourtree_distributed::BranchDecompositionBlock; + + auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator(); + vtkmdiy::Master branch_decomposition_master(comm, + 1, // Use 1 thread, VTK-M will do the treading + -1, // All blocks in memory + 0, // No create function + BranchDecompositionBlock::Destroy); + + // Log the time to create the DIY master for the hyper sweep + timingsStream << " " << std::setw(38) << std::left + << "Create DIY Master (Branch Decomposition):" << timer.GetElapsedTime() + << " seconds" << std::endl; + timer.Start(); + + // Copy data from hierarchical hypersweep to initialize branch decomposition computation + using DistributedContourTreeBlockData = + vtkm::worklet::contourtree_distributed::DistributedContourTreeBlockData; + inputContourTreeMaster.foreach ( + [&](DistributedContourTreeBlockData* currInBlock, const vtkmdiy::Master::ProxyWithLink&) { + BranchDecompositionBlock* newBlock = new BranchDecompositionBlock( + currInBlock->LocalBlockNo, currInBlock->GlobalBlockId, currInBlock->AugmentedTree); + // NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all + // communication is via RegularDecomposer, which sets up its own links. No need + // to keep the pointer, as DIY will "own" it and delete it when no longer needed. + // NOTE: Since we passed a "Destroy" function to DIY master, it will own the local data + // blocks and delete them when done. + branch_decomposition_master.add(currInBlock->GlobalBlockId, newBlock, new vtkmdiy::Link()); + }); + + // Log time to copy the data to the HyperSweepBlock data objects + timingsStream << " " << std::setw(38) << std::left << "Initialize Branch Decomposition Data" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + + vtkmdiy::fix_links(branch_decomposition_master, assigner); + + // Record time to fix the links + timingsStream << " " << std::setw(38) << std::left << "Fix DIY Links (Branch Decomposition)" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + + branch_decomposition_master.foreach ( + [&](BranchDecompositionBlock* b, const vtkmdiy::Master::ProxyWithLink&) { + // Get intrinsic and dependent volume from data set + vtkm::cont::DataSet ds = outputDataSet[b->LocalBlockNo]; + vtkm::cont::ArrayHandle intrinsicVolume = + ds.GetField("IntrinsicVolume").GetData().AsArrayHandle>(); + vtkm::cont::ArrayHandle dependentVolume = + ds.GetField("DependentVolume").GetData().AsArrayHandle>(); + + // Get global size and compute total volume from it + const auto& globalSize = this->MultiBlockSpatialDecomposition.GlobalSize; + vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; + + // Compute local best up and down paths by volume + b->HierarchicalVolumetricBranchDecomposer.LocalBestUpDownByVolume( + &b->HierarchicalContourTree, intrinsicVolume, dependentVolume, totalVolume); + +#ifdef DEBUG_PRINT + VTKM_LOG_S(this->TreeLogLevel, "Before reduction"); + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader( + b->HierarchicalVolumetricBranchDecomposer.BestUpSupernode.GetNumberOfValues(), rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestUpSupernode", b->HierarchicalVolumetricBranchDecomposer.BestUpSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestDownSupernode", b->HierarchicalVolumetricBranchDecomposer.BestDownSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestUpVolume", b->HierarchicalVolumetricBranchDecomposer.BestUpVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestDownVolume", b->HierarchicalVolumetricBranchDecomposer.BestDownVolume, -1, rs); + VTKM_LOG_S(this->TreeLogLevel, rs.str()); + } +#endif + }); + + timingsStream << " " << std::setw(38) << std::left << "LocalBestUpDownByVolume" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + + // Reduce + // partners for merge over regular block grid + vtkmdiy::reduce( + branch_decomposition_master, + assigner, + partners, + vtkm::worklet::contourtree_distributed::ComputeDistributedBranchDecompositionFunctor< + FieldType>{}); + + timingsStream << " " << std::setw(38) << std::left + << "Exchanging best up/down supernode and volume" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + + branch_decomposition_master.foreach ( + [&](BranchDecompositionBlock* b, const vtkmdiy::Master::ProxyWithLink&) { + b->HierarchicalVolumetricBranchDecomposer.CollapseBranches(&b->HierarchicalContourTree, + b->BranchRoots); + }); + + timingsStream << " " << std::setw(38) << std::left << "CollapseBranches" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + + branch_decomposition_master.foreach ([&](BranchDecompositionBlock* b, + const vtkmdiy::Master::ProxyWithLink&) { +#ifdef DEBUG_PRINT + VTKM_LOG_S(this->TreeLogLevel, + b->HierarchicalVolumetricBranchDecomposer.PrintBranches(&b->HierarchicalContourTree, + b->BranchRoots)); +#endif + vtkm::cont::Field branchRootField( + "BranchRoots", vtkm::cont::Field::Association::WholeMesh, b->BranchRoots); + outputDataSet[b->LocalBlockNo].AddField(branchRootField); + }); + + timingsStream << " " << std::setw(38) << std::left + << "Creating Branch Decomposition Output Data" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); +} //----------------------------------------------------------------------------- template @@ -854,10 +1008,11 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( 1, // Use 1 thread, VTK-M will do the treading -1, // All blocks in memory 0, // No create function (since all blocks in memory) - DistributedContourTreeBlockData::destroy); + DistributedContourTreeBlockData::Destroy); // ... and record time for creating the DIY master - timingsStream << " " << std::setw(38) << std::left << "Create DIY Master" + timingsStream << " " << std::setw(38) << std::left + << "Create DIY Master (Distributed Contour Tree)" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); @@ -946,7 +1101,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( // NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all // communication is via RegularDecomposer, which sets up its own links. No need // to keep the pointer, as DIY will "own" it and delete it when no longer needed. - // NOTE: Since we passed a "destroy" function to DIY master, it will own the local data + // NOTE: Since we passed a "Destroy" function to DIY master, it will own the local data // blocks and delete them when done. master.add(vtkmdiyLocalBlockGids[bi], newBlock, new vtkmdiy::Link); } // for @@ -1020,7 +1175,8 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( vtkmdiy::fix_links(master, assigner); // Record time to fix the links - timingsStream << " " << std::setw(38) << std::left << "Fix DIY Links" + timingsStream << " " << std::setw(38) << std::left + << "Fix DIY Links (Distributed Contour Tree)" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); @@ -1192,7 +1348,7 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( timer.Start(); } - // ******** 4. Create output data set ******** + // ******** 5. Create output data set ******** std::vector hierarchicalTreeOutputDataSet(master.size()); master.foreach ( [&](DistributedContourTreeBlockData* blockData, const vtkmdiy::Master::ProxyWithLink&) { @@ -1277,6 +1433,12 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( master, assigner, partners, FieldType{}, timingsStream, hierarchicalTreeOutputDataSet); } + if (this->ComputeHierarchicalVolumetricBranchDecomposition) + { + this->ComputeBranchDecomposition( + master, assigner, partners, FieldType{}, timingsStream, hierarchicalTreeOutputDataSet); + } + VTKM_LOG_S(this->TimingsLogLevel, std::endl << " ------------ DoPostExecute Timings ------------" << std::endl diff --git a/vtkm/filter/testing/TestingContourTreeUniformDistributedFilter.h b/vtkm/filter/testing/TestingContourTreeUniformDistributedFilter.h index c41def4d0..47fb7add2 100644 --- a/vtkm/filter/testing/TestingContourTreeUniformDistributedFilter.h +++ b/vtkm/filter/testing/TestingContourTreeUniformDistributedFilter.h @@ -292,8 +292,9 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed( // work when only using boundary extrema !useMarchingCubes, useMarchingCubes, - false, - false, + false, // no augmentationa + false, // no branch decmompostion + false, // no save dot vtkm::cont::LogLevel::UserVerboseLast, vtkm::cont::LogLevel::UserVerboseLast); filter.SetActiveField(fieldName); diff --git a/vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h b/vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h new file mode 100644 index 000000000..17044c5b7 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h @@ -0,0 +1,107 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// +// This code is an extension of the algorithm presented in the paper: +// Parallel Peak Pruning for Scalable SMP Contour Tree Computation. +// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens. +// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization +// (LDAV), October 2016, Baltimore, Maryland. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_branchdecompositionblock_h +#define vtk_m_worklet_contourtree_distributed_branchdecompositionblock_h + +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ + +template +struct BranchDecompositionBlock +{ + BranchDecompositionBlock( + vtkm::Id localBlockNo, + int globalBlockId, + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree& + hierarchicalContourTree) + : LocalBlockNo(localBlockNo) + , GlobalBlockId(globalBlockId) + , HierarchicalContourTree(hierarchicalContourTree) + { + // Import/alias for readability + using vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + const vtkm::Id nSupernodes = HierarchicalContourTree.Supernodes.GetNumberOfValues(); + //VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Allocating " << nSupernodes << " valies in BranchRoots"); + BranchRoots.AllocateAndFill(nSupernodes, NO_SUCH_ELEMENT); + } + + // Block metadata TODO/FIXME: Check whether really needed + vtkm::Id LocalBlockNo; + int GlobalBlockId; + + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree& + HierarchicalContourTree; + vtkm::worklet::contourtree_distributed::HierarchicalVolumetricBranchDecomposer< + ContourTreeDataFieldType> + HierarchicalVolumetricBranchDecomposer; + vtkm::cont::ArrayHandle BranchRoots; + + // Destroy function allowing DIY to own blocks and clean them up after use + static void Destroy(void* b) + { + delete static_cast*>(b); + } +}; + +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm +#endif diff --git a/vtkm/worklet/contourtree_distributed/CMakeLists.txt b/vtkm/worklet/contourtree_distributed/CMakeLists.txt index 19ed032ad..f0294879a 100644 --- a/vtkm/worklet/contourtree_distributed/CMakeLists.txt +++ b/vtkm/worklet/contourtree_distributed/CMakeLists.txt @@ -11,12 +11,15 @@ set(headers BoundaryTree.h BoundaryTreeMaker.h + BranchDecompositionBlock.h CombineHyperSweepBlockFunctor.h + ComputeDistributedBranchDecompositionFunctor.h ComputeDistributedContourTreeFunctor.h ContourTreeBlockData.h DistributedContourTreeBlockData.h HierarchicalAugmenter.h HierarchicalAugmenterFunctor.h + HierarchicalVolumetricBranchDecomposer.h HierarchicalContourTree.h HierarchicalHyperSweeper.h HyperSweepBlock.h @@ -36,3 +39,4 @@ add_subdirectory(tree_grafter) add_subdirectory(hierarchical_contour_tree) add_subdirectory(hierarchical_hyper_sweeper) add_subdirectory(hierarchical_augmenter) +add_subdirectory(hierarchical_volumetric_branch_decomposer) diff --git a/vtkm/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h b/vtkm/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h index ee5be0010..aa6945964 100644 --- a/vtkm/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h +++ b/vtkm/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h @@ -123,7 +123,7 @@ struct CobmineHyperSweepBlockFunctor // TODO/FIXME: Is there a way to do an in-place transform without a temporary array? vtkm::cont::Algorithm::Transform( intrinsicVolumeView, incomingIntrinsicVolume, tempSum, vtkm::Sum()); - vtkm::cont::Algorithm::Copy(tempSum, intrinsicVolumeView); + vtkm::cont::ArrayCopy(tempSum, intrinsicVolumeView); auto dependentVolumeView = make_ArrayHandleView(b->DependentVolume, 0, numSupernodesToProcess); @@ -132,7 +132,7 @@ struct CobmineHyperSweepBlockFunctor // TODO/FIXME: Is there a way to do an in-place transform without a temporary array? vtkm::cont::Algorithm::Transform( dependentVolumeView, incomingDependentVolume, tempSum, vtkm::Sum()); - vtkm::cont::Algorithm::Copy(tempSum, dependentVolumeView); + vtkm::cont::ArrayCopy(tempSum, dependentVolumeView); } } @@ -156,9 +156,9 @@ struct CobmineHyperSweepBlockFunctor // without copy. enqueue does not accept ArrayHandleView as input as it // is not trivially copyable vtkm::cont::ArrayHandle sendIntrinsicVolume; - vtkm::cont::Algorithm::Copy(intrinsicVolumeView, sendIntrinsicVolume); + vtkm::cont::ArrayCopy(intrinsicVolumeView, sendIntrinsicVolume); vtkm::cont::ArrayHandle sendDependentVolume; - vtkm::cont::Algorithm::Copy(dependentVolumeView, sendDependentVolume); + vtkm::cont::ArrayCopy(dependentVolumeView, sendDependentVolume); // Send necessary data portions rp.enqueue(target, sendIntrinsicVolume); diff --git a/vtkm/worklet/contourtree_distributed/ComputeDistributedBranchDecompositionFunctor.h b/vtkm/worklet/contourtree_distributed/ComputeDistributedBranchDecompositionFunctor.h new file mode 100644 index 000000000..167db6a41 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/ComputeDistributedBranchDecompositionFunctor.h @@ -0,0 +1,238 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// +// This code is an extension of the algorithm presented in the paper: +// Parallel Peak Pruning for Scalable SMP Contour Tree Computation. +// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens. +// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization +// (LDAV), October 2016, Baltimore, Maryland. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed__h +#define vtk_m_worklet_contourtree_distributed__h + +#include +#include +#include + +// clang-format off +VTKM_THIRDPARTY_PRE_INCLUDE +#include +VTKM_THIRDPARTY_POST_INCLUDE +// clang-format on + +#ifdef DEBUG_PRINT +#define DEBUG_PRINT_COMBINED_BLOCK_IDS +#endif + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ + +template +struct ComputeDistributedBranchDecompositionFunctor +{ + void operator()( + vtkm::worklet::contourtree_distributed::BranchDecompositionBlock* b, + const vtkmdiy::ReduceProxy& rp, // communication proxy + const vtkmdiy::RegularSwapPartners& // partners of the current block (unused) + ) const + { + // Get our rank and DIY id + //const vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank(); + const auto selfid = rp.gid(); + + // Aliases to reduce verbosity + auto& branchDecomposer = b->HierarchicalVolumetricBranchDecomposer; + + std::vector incoming; + rp.incoming(incoming); + for (const int ingid : incoming) + { + // NOTE/IMPORTANT: In each round we should have only one swap partner (despite for-loop here). + // If that assumption does not hold, it will break things. + // NOTE/IMPORTANT: This assumption only holds if the number of blocks is a power of two. + // Otherwise, we may need to process more than one incoming block + if (ingid != selfid) + { +#ifdef DEBUG_PRINT_COMBINED_BLOCK_IDS + int incomingGlobalBlockId; + rp.dequeue(ingid, incomingGlobalBlockId); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + "Combining local block " << b->GlobalBlockId << " with incomoing block " + << incomingGlobalBlockId); +#endif + + // Receive data from swap partner + vtkm::cont::ArrayHandle incomingBestUpVolume; + rp.dequeue(ingid, incomingBestUpVolume); + vtkm::cont::ArrayHandle incomingBestUpSupernode; + rp.dequeue(ingid, incomingBestUpSupernode); + + vtkm::cont::ArrayHandle incomingBestDownVolume; + rp.dequeue(ingid, incomingBestDownVolume); + vtkm::cont::ArrayHandle incomingBestDownSupernode; + rp.dequeue(ingid, incomingBestDownSupernode); + + vtkm::Id prefixLength = vtkm::cont::ArrayGetValue( + 0, b->HierarchicalContourTree.FirstSupernodePerIteration[rp.round() - 1]); + +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Prefix length is " << prefixLength); + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader( + incomingBestUpSupernode.GetNumberOfValues(), rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingBestUpSupernode", incomingBestUpSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingBestDownSupernode", incomingBestDownSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingBestUpVolume", incomingBestUpVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingBestDownVolume", incomingBestDownVolume, -1, rs); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); + } +#endif + // NOTE: We are processing input data from the previous round, hence, get + // first supernide per ieteration from previous round + + // Create 'views' to restrict worklet to relevant portion of arrays + auto bestUpVolumeView = + make_ArrayHandleView(branchDecomposer.BestUpVolume, 0, prefixLength); + auto bestUpSupernodeView = + make_ArrayHandleView(branchDecomposer.BestUpSupernode, 0, prefixLength); + auto bestDownVolumeView = + make_ArrayHandleView(branchDecomposer.BestDownVolume, 0, prefixLength); + auto bestDownSupernodeView = + make_ArrayHandleView(branchDecomposer.BestDownSupernode, 0, prefixLength); + + // Check if swap partner knows a better up /down and update + vtkm::cont::Invoker invoke; + invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + FindBestSupernodeWorklet{}, + incomingBestUpVolume, + incomingBestUpSupernode, + bestUpVolumeView, + bestUpSupernodeView); + + invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + FindBestSupernodeWorklet{}, + incomingBestDownVolume, + incomingBestDownSupernode, + bestDownVolumeView, + bestDownSupernodeView); + +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "After round " << rp.round() - 1); + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader( + branchDecomposer.BestUpSupernode.GetNumberOfValues(), rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestUpSupernode", branchDecomposer.BestUpSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestDownSupernode", branchDecomposer.BestDownSupernode, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestUpVolume", branchDecomposer.BestUpVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BestDownVolume", branchDecomposer.BestDownVolume, -1, rs); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); + } +#endif + } + } + + for (int cc = 0; cc < rp.out_link().size(); ++cc) + { + auto target = rp.out_link().target(cc); + if (target.gid != selfid) + { +#ifdef DEBUG_PRINT_COMBINED_BLOCK_IDS + rp.enqueue(target, b->GlobalBlockId); +#endif + // Determine which portion of up/down volume/supernode to send + vtkm::Id prefixLength = vtkm::cont::ArrayGetValue( + 0, b->HierarchicalContourTree.FirstSupernodePerIteration[rp.round()]); + + // Create 'views' to restrict sending to relevant portion of arrays + auto bestUpVolumeView = + make_ArrayHandleView(branchDecomposer.BestUpVolume, 0, prefixLength); + auto bestUpSupernodeView = + make_ArrayHandleView(branchDecomposer.BestUpSupernode, 0, prefixLength); + auto bestDownVolumeView = + make_ArrayHandleView(branchDecomposer.BestDownVolume, 0, prefixLength); + auto bestDownSupernodeView = + make_ArrayHandleView(branchDecomposer.BestDownSupernode, 0, prefixLength); + + // TODO/FIXME: Check if it is possible to send a portion of the arrays + // without copy. enqueue does not accept ArrayHandleView as input as it + // is not trivially copyable + vtkm::cont::ArrayHandle sendBestUpVolume; + vtkm::cont::Algorithm::Copy(bestUpVolumeView, sendBestUpVolume); + vtkm::cont::ArrayHandle sendBestUpSupernode; + vtkm::cont::Algorithm::Copy(bestUpSupernodeView, sendBestUpSupernode); + vtkm::cont::ArrayHandle sendBestDownVolume; + vtkm::cont::Algorithm::Copy(bestDownVolumeView, sendBestDownVolume); + vtkm::cont::ArrayHandle sendBestDownSupernode; + vtkm::cont::Algorithm::Copy(bestDownSupernodeView, sendBestDownSupernode); + + rp.enqueue(target, sendBestUpVolume); + rp.enqueue(target, sendBestUpSupernode); + rp.enqueue(target, sendBestDownVolume); + rp.enqueue(target, sendBestDownSupernode); + } + } + } +}; + +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h b/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h index 565129706..9541e15b1 100644 --- a/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h +++ b/vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h @@ -177,11 +177,11 @@ public: // save the corresponding .gv file for the contour tree mesh std::string contourTreeMeshFileName = std::string("Rank_") + std::to_string(static_cast(rank)) + std::string("_Block_") + - std::to_string(static_cast(block->BlockIndex)) + "_Round_" + + std::to_string(static_cast(block->LocalBlockNo)) + "_Round_" + std::to_string(rp.round()) + "_Partner_" + std::to_string(ingid) + std::string("_Step_0_Combined_Mesh.gv"); std::string contourTreeMeshLabel = std::string("Block ") + - std::to_string(static_cast(block->BlockIndex)) + " Round " + + std::to_string(static_cast(block->LocalBlockNo)) + " Round " + std::to_string(rp.round()) + " Partner " + std::to_string(ingid) + std::string(" Step 0 Combined Mesh"); std::string contourTreeMeshString = diff --git a/vtkm/worklet/contourtree_distributed/DistributedContourTreeBlockData.h b/vtkm/worklet/contourtree_distributed/DistributedContourTreeBlockData.h index 387770b70..88bdfff6f 100644 --- a/vtkm/worklet/contourtree_distributed/DistributedContourTreeBlockData.h +++ b/vtkm/worklet/contourtree_distributed/DistributedContourTreeBlockData.h @@ -94,7 +94,7 @@ struct DistributedContourTreeBlockData vtkm::worklet::contourtree_distributed::HierarchicalContourTree AugmentedTree; // Destroy function allowing DIY to own blocks and clean them up after use - static void destroy(void* b) + static void Destroy(void* b) { delete static_cast*>(b); } diff --git a/vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h b/vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h index 8f57b9b7b..a39c6b4c7 100644 --- a/vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h +++ b/vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h @@ -75,6 +75,7 @@ #include #include #include +#include #include #include #include @@ -188,6 +189,13 @@ public: this->DataValues); } + /// routine to create a FindSuperArcBetweenNodes object that we can use as an input for worklets to call the function + VTKM_CONT + FindSuperArcBetweenNodes GetFindSuperArcBetweenNodes() const + { + return FindSuperArcBetweenNodes(this->Superarcs); + } + /// routine to initialize the hierarchical tree with the top level tree VTKM_CONT void Initialize(vtkm::Id numRounds, diff --git a/vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h b/vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h new file mode 100644 index 000000000..e9d7c7b87 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h @@ -0,0 +1,658 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//======================================================================================= +// +// Parallel Peak Pruning v. 2.0 +// +// Started June 15, 2017 +// +// Copyright Hamish Carr, University of Leeds +// +// HierarchicalVolumetricBranchDecomposer.h +// +//======================================================================================= +// +// COMMENTS: +// +// This class computes the branch decomposition by volume for a given hierarchical +// contour tree. +// +// It takes as input arrays of dependent and intrinsic volumes for each superarc +// (it needs both, in order to compute the dependent volume at each end of each superarc) +// +// Recall from the non-hierarchical version that in order to compute the branch decomposition, +// we need to choose the "best up" and "best down" superarc for each supernode - i.e. the +// superarc with the largest dependent volume. Since we only wish to compare superarcs that +// meet at a given supernode, we tiebreak by always taking the superarc whose "other" end +// has a higher ID. +// +// Once the best up & best down have been found for each supernode, branches are identified +// with (essentially) a graph connectivity computation. +// +// Conceptually, each superarc is a vertex in a new (temporary) graph. For each supernode, the +// "best up" superarc is connected to the "best down" superarc. This defines a graph in which +// each branch is a connected component. A single path-doubling pass then collects the branches +// +// In the non-hierarchical version, this was done with supernode IDs (I can't remember why), +// with the upper end of each branch being treated as the root node. +// +// To construct the hierarchical branch decomposition, we assume that the hierarchical contour +// tree has already been augmented with all attachment points. If not, the code may produce +// undefined results. +// +// In the first step, we will run a local routine for each rank to determine the best up / down +// as far as the rank knows. We will then do a fan-in swap to determine the best up / down for +// shared vertices. At the end of this step, all ranks will share the knowledge of the best +// up / down superarc, stored as: +// i. the superarc ID, which may be reused on other ranks +// ii. the global ID of the outer end of that superarc, which is unique across all ranks +// iii. the volume dependent on that superarc +// +// In the second stage, each rank will do a local computation of the branches. However, most ranks +// will not have the full set of supernodes / superarcs for each branch, even (or especially) +// for the master branch. It is therefore a bad idea to collapse to the upper end of the branch +// as we did in the non-hierarchical version. +// +// Instead, we will define the root of each component to be the most senior superarc ID. This will +// be canonical, because of the way we construct the hierarchical tree, with low superarc IDs +// occurring at higher levels of the tree, so all shared superarcs are a prefix set. Therefore, +// the most senior superarc ID will always indicate the highest level of the tree through which the +// branch passes, and is safe. Moreover, it is not necessary for each rank to determine the full +// branch, merely the part of the branch that passes through the superarcs it tracks. It may even +// happen that no single rank stores the entire branch, as for example if the global minimum +// and maximum are interior to different ranks. +// +// Note that most senior means testing iteration, round, then ID +// +// RESIZE SEMANTICS: Oliver Ruebel has asked for the semantics of all resize() calls to be annotated +// in order to ease porting to vtkm. These will be flagged with a RESIZE SEMANTICS: comment, and will +// generally fall into several patterns: +// 1. FIXED: Resize() is used to initialize the array size for an array that will never change size +// 2. COMPRESS: Resize() is used after a compression operation (eg remove_if()) so that the +// array size() call does not include the elements removed. This is a standard +// C++ pattern, but could be avoided by storing an explicit element count (curiously, +// the std::vector class does exactly this with logical vs. physical array sizes). +// 3. MULTI-COMPRESS: Resize() may also be used (as in the early stages of the PPP algorithm) to give +// a collapsing array size of working elements. Again, this could on principle by +// avoided with an array count, but is likely to be intricate. +// +//======================================================================================= + + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_h + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEBUG_PRINT +#define DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER +#endif + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ + +/// Facture class for augmenting the hierarchical contour tree to enable computations of measures, e.g., volumne +template +class HierarchicalVolumetricBranchDecomposer +{ // class HierarchicalVolumetricBranchDecomposer +public: + /// we will want arrays for swapping with our partners, holding the best up/down superarc & the corresponding volume + /// the best up/down will be in local supernode IDs initially, but during the swap will need to be global node IDs + vtkm::worklet::contourtree_augmented::IdArrayType BestUpSupernode; + vtkm::worklet::contourtree_augmented::IdArrayType BestDownSupernode; + vtkm::worklet::contourtree_augmented::IdArrayType BestUpVolume; + vtkm::worklet::contourtree_augmented::IdArrayType BestDownVolume; + + /// working arrays - kept at class level to simplify debug print + vtkm::worklet::contourtree_augmented::IdArrayType UpVolume; + vtkm::worklet::contourtree_augmented::IdArrayType DownVolume; + + /// Because we want to use array allocation in main, we need a default constructor + /// Unfortunately, that means we have problems with references, and the workaround is to have + /// the references passed to the individual functions + HierarchicalVolumetricBranchDecomposer(){}; + + /// routines to compute branch decomposition by volume + /// WARNING: we now have two types of hierarchical tree sharing a data structure: + /// I. hierarchical tree without augmentation + /// II. hierarchical tree with augmentation + /// We only expect to call this for II, but it's wiser to make sure that it computes for I as well. + /// Also, this code is substantially identical to ContourTreeMaker::ComputeVolumeBranchDecomposition() + /// except for: + /// A. it has to deal with the round/iteration paradigm of hierarchical trees, and + /// B. Stages III-IV in particular are modified + /// C. Several stages involve fan-ins + /// The principal reason for the modifications in B. is that the old code collapses branches to their maximum + /// which is often a leaf. In the hierarchical version, the leaf will often not be represented on all ranks, so + /// we modify it to collapse towards the "most senior". This will be easiest if we collapse by superarc IDs instead of supernode IDs + /// For C., we have to break the code into separate routines so that the fan-in MPI can be outside this unit. + /// + /// WARNING! WARNING! WARNING! + /// In the non-hierarchical version, the last (virtual root) superarc goes from the highest ID supernode to NO_SUCH_ELEMENT + /// If it was included in the sorts, this could cause problems + /// The (simple) way out of this was to set nSuperarcs = nSupernodes - 1 when copying our temporary list of superarcs + /// that way we don't use it at all. + /// In the hierarchical version, this no longer works, because attachment points may also have virtual superarcs + /// So we either need to compress them out (an extra log step) or ignore them in the later loop. + /// Of the two, compressing them out is safer + /// + /// routine that determines the best upwards/downwards edges at each vertex + /// Unlike the local version, the best might only be stored on another rank + /// so we will compute the locally best up or down, then swap until all ranks choose the same best + + void LocalBestUpDownByVolume( + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree* + hierarchicalTree, + const vtkm::cont::ArrayHandle& intrinsicValues, + const vtkm::cont::ArrayHandle& dependentValues, + vtkm::Id totalVolume); + + /// routine to compute the local set of superarcs that root at a given one + void CollapseBranches(const vtkm::worklet::contourtree_distributed::HierarchicalContourTree< + FieldType>* hierarchicalTree, + vtkm::worklet::contourtree_augmented::IdArrayType& branchRoot); + + // routine to print branches + std::string PrintBranches(const vtkm::worklet::contourtree_distributed::HierarchicalContourTree< + FieldType>* hierarchicalTree, + const vtkm::worklet::contourtree_augmented::IdArrayType& branchRoot); + static std::string PrintBranches(const vtkm::cont::DataSet& ds); + + /// debug routine + std::string DebugPrint(std::string message, const char* fileName, long lineNum); + +private: + /// Used internally to Invoke worklets + vtkm::cont::Invoker Invoke; + +}; // class HierarchicalVolumetricBranchDecomposer + + +template +void HierarchicalVolumetricBranchDecomposer::LocalBestUpDownByVolume( + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree* + hierarchicalTree, + const vtkm::cont::ArrayHandle& intrinsicValues, + const vtkm::cont::ArrayHandle& dependentValues, + vtkm::Id totalVolume) +{ + // LocalBestUpDownByVolume + // STAGE I: Allocate memory for our arrays + vtkm::Id nSupernodes = hierarchicalTree->Supernodes.GetNumberOfValues(); + // WARNING: This differs from the non-hierarchical version by using the full size *WITH* virtual superarcs + vtkm::Id nSuperarcs = hierarchicalTree->Superarcs.GetNumberOfValues(); + + // set up a list of superarcs as Edges for reference in our comparator + vtkm::worklet::contourtree_augmented::EdgePairArray superarcList; + superarcList.Allocate(nSuperarcs); + VTKM_ASSERT(hierarchicalTree->Superarcs.GetNumberOfValues() == superarcList.GetNumberOfValues()); + this->Invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + LocalBestUpDownByVolumeInitSuperarcListWorklet{}, // the worklet + hierarchicalTree->Superarcs, // input + superarcList // output + ); + + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + { + std::stringstream resultStream; + vtkm::worklet::contourtree_augmented::PrintHeader(superarcList.GetNumberOfValues(), + resultStream); + vtkm::worklet::contourtree_augmented::PrintEdgePairArray( + "Superarc List", superarcList, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); + } +#endif + + // create a list of the non-virtual superarcs + vtkm::worklet::contourtree_augmented::IdArrayType actualSuperarcs; + // and fill it up with index values [0, 1, 2 ... nSuperarcs-1] while simultaneously stream compacting the + // values by keeping only those indices where the hierarchicalTree->Superarcs is not NoSuchElement. + vtkm::cont::Algorithm::CopyIf(vtkm::cont::ArrayHandleIndex(nSuperarcs), //input + hierarchicalTree->Superarcs, // stencil + actualSuperarcs, // output target array + vtkm::worklet::contourtree_augmented::NotNoSuchElementPredicate{}); + // NOTE: The behavior here is slightly different from the original implementation, as the original code + // here does not resize actualSuperarcs but keeps it at the full length of nSuperacs and instead + // relies on the nActualSuperarcs parameter. However, the extra values are never used, so compacting + // the array here should be fine. + vtkm::Id nActualSuperarcs = actualSuperarcs.GetNumberOfValues(); + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + { + std::stringstream resultStream; + vtkm::worklet::contourtree_augmented::PrintHeader(nActualSuperarcs, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Actual Superarcs", actualSuperarcs, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); + } +#endif + + // and set up arrays for the best upwards, downwards superarcs at each supernode + // initialize everything to NO_SUCH_ELEMENT for safety (we will test against this, so it's necessary) + // Set up temporary constant arrays for each group of arrays and initalize the arrays + // Initalize the arrays + using vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + this->UpVolume.AllocateAndFill(nSuperarcs, 0); + this->DownVolume.AllocateAndFill(nSuperarcs, 0); + this->BestUpSupernode.AllocateAndFill(nSupernodes, NO_SUCH_ELEMENT); + this->BestDownSupernode.AllocateAndFill(nSupernodes, NO_SUCH_ELEMENT); + this->BestUpVolume.AllocateAndFill(nSupernodes, 0); + this->BestDownVolume.AllocateAndFill(nSupernodes, 0); + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + VTKM_LOG_S(vtkm::cont::LogLevel::Info, DebugPrint("Arrays Allocated", __FILE__, __LINE__)); +#endif + + // STAGE II: Pick the best (largest volume) edge upwards and downwards + // II A. Compute the up / down volumes for indirect sorting + // this is the same in spirit as ContourTreeMaker::ComputeVolumeBranchDecomposition() STAGE II A. + // given that we have already suppressed the non-virtual superarcs + // however, in this case, we need to use the actualSuperarcs array instead of the main array + { + vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + LocalBestUpDownByVolumeBestUpDownEdgeWorklet bestUpDownEdgeWorklet(totalVolume); + // permut input and output arrays here so we can use FieldIn and FieldOut to + // avoid the use of WholeArray access in the worklet + auto permutedHierarchicalTreeSuperarcs = vtkm::cont::make_ArrayHandlePermutation( + actualSuperarcs, hierarchicalTree->Superarcs); // input + auto permutedDependetValues = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, dependentValues); // input + auto permutedIntrinsicValues = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, intrinsicValues); // input + auto permutedUpVolume = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, this->UpVolume); // output + auto permitedDownVolume = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, this->DownVolume); // outout + + this->Invoke(bestUpDownEdgeWorklet, // the worklet + permutedHierarchicalTreeSuperarcs, // input + permutedDependetValues, // input + permutedIntrinsicValues, // input + permutedUpVolume, // output + permitedDownVolume // outout + ); + } + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + VTKM_LOG_S(vtkm::cont::LogLevel::Info, DebugPrint("Volume Arrays Set Up", __FILE__, __LINE__)); + { + std::stringstream resultStream; + vtkm::worklet::contourtree_augmented::PrintHeader(superarcList.GetNumberOfValues(), + resultStream); + vtkm::worklet::contourtree_augmented::PrintEdgePairArray( + "Superarc List", superarcList, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); + } +#endif + // II B. Pick the best downwards volume by sorting on upper vertex then processing by segments (segmented by vertex) + // II B 1. Sort the superarcs by upper vertex + // NB: We reuse the actual superarcs list here - this works because we have indexed the volumes on the underlying superarc ID + // NB 2: Notice that we only sort the "actual" ones - this is to avoid unnecessary resize() calls in vtkm later on + { + vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + SuperArcVolumetricComparatorIndirectGlobalIdComparator + SuperArcVolumetricComparatorIndirectGlobalIdComparator( + this->UpVolume, superarcList, hierarchicalTree->RegularNodeGlobalIds, false); + vtkm::cont::Algorithm::Sort(actualSuperarcs, + SuperArcVolumetricComparatorIndirectGlobalIdComparator); + } + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + { + std::stringstream resultStream; + resultStream + << "Actual Superarc List After Sorting By High End (Full Array, including ignored elements)" + << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nActualSuperarcs, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Actual Superarcs", actualSuperarcs, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); + } +#endif + // II B 2. Per vertex, best superarc writes to the best downward array + { + auto permutedUpVolume = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, this->UpVolume); + this->Invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + LocalBestUpDownByVolumeWorklet{ nActualSuperarcs }, + actualSuperarcs, // input + superarcList, // input + permutedUpVolume, // input + hierarchicalTree->RegularNodeGlobalIds, // input + hierarchicalTree->Supernodes, // input + this->BestDownSupernode, // output + this->BestDownVolume // output + ); + } +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + DebugPrint("BestDownSupernode Written", __FILE__, __LINE__)); +#endif + + // II B 3. Repeat for lower vertex + { + vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + SuperArcVolumetricComparatorIndirectGlobalIdComparator + SuperArcVolumetricComparatorIndirectGlobalIdComparator( + this->DownVolume, superarcList, hierarchicalTree->RegularNodeGlobalIds, true); + vtkm::cont::Algorithm::Sort(actualSuperarcs, + SuperArcVolumetricComparatorIndirectGlobalIdComparator); + } + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + { + std::stringstream resultStream; + resultStream + << "Actual Superarc List After Sorting By Low End (Full Array, including ignored elements)" + << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nActualSuperarcs, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Actual Superarcs", actualSuperarcs, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); + } +#endif + + // II B 2. Per vertex, best superarc writes to the best upward array + { + auto permutedDownVolume = + vtkm::cont::make_ArrayHandlePermutation(actualSuperarcs, this->DownVolume); + this->Invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + LocalBestUpDownByVolumeWorklet{ nActualSuperarcs }, + actualSuperarcs, // input + superarcList, // input + permutedDownVolume, // input + hierarchicalTree->RegularNodeGlobalIds, // input + hierarchicalTree->Supernodes, // input + this->BestUpSupernode, // output + this->BestUpVolume // output + ); + } + +#ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + DebugPrint("Local Best Up/Down Computed", __FILE__, __LINE__)); +#endif +} // LocalBestUpDownByVolume + +template +void HierarchicalVolumetricBranchDecomposer::CollapseBranches( + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree* + hierarchicalTree, + vtkm::worklet::contourtree_augmented::IdArrayType& branchRoot) +{ // CollapseBranches + // initialise the superarcs to be their own branch roots + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + "Initializing branch root with " << branchRoot.GetNumberOfValues() << " values."); + vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(branchRoot.GetNumberOfValues()), + branchRoot); + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Convert to superarc ID."); + // For each supernode, convert the best up into a superarc ID + { + vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + CollapseBranchesWorklet collapseBranchesWorklet; + auto findRegularByGlobal = hierarchicalTree->GetFindRegularByGlobal(); + auto findSuperArcBetweenNodes = hierarchicalTree->GetFindSuperArcBetweenNodes(); + this->Invoke(collapseBranchesWorklet, // the worklet + this->BestUpSupernode, // input + this->BestDownSupernode, // input + findRegularByGlobal, // input ExecutionObject + findSuperArcBetweenNodes, // input ExecutionObject + hierarchicalTree->Regular2Supernode, // input + hierarchicalTree->WhichRound, // input + branchRoot); + } + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Pointer doubling."); + // OK. We've now initialized it, and can use pointer-doubling + // Compute the number of log steps required + vtkm::Id nLogSteps = 1; + for (vtkm::Id shifter = branchRoot.GetNumberOfValues(); shifter != 0; shifter >>= 1) + { + nLogSteps++; + } + + // loop that many times, pointer-doubling + for (vtkm::Id iteration = 0; iteration < nLogSteps; iteration++) + { // per iteration + // loop through the vertices, updating + vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer:: + CollapseBranchesPointerDoublingWorklet collapseBranchesPointerDoublingWorklet; + this->Invoke(collapseBranchesPointerDoublingWorklet, branchRoot); + } // per iteration + +} // CollapseBranches + + +// routine to print branches +template +std::string HierarchicalVolumetricBranchDecomposer::PrintBranches( + const vtkm::worklet::contourtree_distributed::HierarchicalContourTree* + hierarchicalTree, + const vtkm::worklet::contourtree_augmented::IdArrayType& branchRoot) +{ // PrintBranches + // we want to dump out the branches as viewed by this rank. + // most of the processing will be external, so we keep this simple. + // For each end of the superarc, we print out value & global ID prefixed by global ID of the branch root + // The external processing will then sort them to construct segments (as usual) in the array + // then a post-process can find the first and last in each segment & print out the branch + // In order for the sort to work lexicographically, we need to print out in the following order: + // I Branch Root Global ID + // II Supernode Value + // III Supernode Global ID + std::stringstream resultStream; + // loop through the individual superarcs + auto hierarchicalTreeSuperarcsPortal = hierarchicalTree->Superarcs.ReadPortal(); + auto hierarchicalTreeSupernodesPortal = hierarchicalTree->Supernodes.ReadPortal(); + auto hierarchicalTreeRegularNodeGlobalIdsPortal = + hierarchicalTree->RegularNodeGlobalIds.ReadPortal(); + auto hierarchicalTreeDataValuesPortal = hierarchicalTree->DataValues.ReadPortal(); + auto branchRootPortal = branchRoot.ReadPortal(); + + for (vtkm::Id superarc = 0; superarc < hierarchicalTree->Superarcs.GetNumberOfValues(); + superarc++) + { // per superarc + // explicit test for root superarc / attachment points + if (vtkm::worklet::contourtree_augmented::NoSuchElement( + hierarchicalTreeSuperarcsPortal.Get(superarc))) + { + continue; + } + + // now retrieve the branch root's global ID + vtkm::Id branchRootSuperId = branchRootPortal.Get(superarc); + vtkm::Id branchRootRegularId = hierarchicalTreeSupernodesPortal.Get(branchRootSuperId); + vtkm::Id branchRootGlobalId = + hierarchicalTreeRegularNodeGlobalIdsPortal.Get(branchRootRegularId); + + // now retrieve the global ID & value for each end & output them + vtkm::Id superFromRegularId = hierarchicalTreeSupernodesPortal.Get(superarc); + vtkm::Id superFromGlobalId = hierarchicalTreeRegularNodeGlobalIdsPortal.Get(superFromRegularId); + FieldType superFromValue = hierarchicalTreeDataValuesPortal.Get(superFromRegularId); + resultStream << branchRootGlobalId << "\t" << superFromValue << "\t" << superFromGlobalId + << std::endl; + + // now retrieve the global ID & value for each end & output them + vtkm::Id superToRegularId = vtkm::worklet::contourtree_augmented::MaskedIndex( + hierarchicalTreeSuperarcsPortal.Get(superarc)); + vtkm::Id superToGlobalId = hierarchicalTreeRegularNodeGlobalIdsPortal.Get(superToRegularId); + FieldType superToValue = hierarchicalTreeDataValuesPortal.Get(superToRegularId); + resultStream << branchRootGlobalId << "\t" << superToValue << "\t" << superToGlobalId + << std::endl; + } // per superarc + return resultStream.str(); +} // PrintBranches + +// routine to print branches +template +std::string HierarchicalVolumetricBranchDecomposer::PrintBranches( + const vtkm::cont::DataSet& ds) +{ // PrintBranches + // we want to dump out the branches as viewed by this rank. + // most of the processing will be external, so we keep this simple. + // For each end of the superarc, we print out value & global ID prefixed by global ID of the branch root + // The external processing will then sort them to construct segments (as usual) in the array + // then a post-process can find the first and last in each segment & print out the branch + // In order for the sort to work lexicographically, we need to print out in the following order: + // I Branch Root Global ID + // II Supernode Value + // III Supernode Global ID + std::stringstream resultStream; + // loop through the individual superarcs + + + auto hierarchicalTreeSuperarcsAH = + ds.GetField("Superarcs").GetData().AsArrayHandle>(); + auto hierarchicalTreeSuperarcsPortal = hierarchicalTreeSuperarcsAH.ReadPortal(); + vtkm::Id nSuperarcs = hierarchicalTreeSuperarcsAH.GetNumberOfValues(); + auto hierarchicalTreeSupernodesPortal = ds.GetField("Supernodes") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto hierarchicalTreeRegularNodeGlobalIdsPortal = + ds.GetField("RegularNodeGlobalIds") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto hierarchicalTreeDataValuesPortal = ds.GetField("DataValues") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto branchRootPortal = ds.GetField("BranchRoots") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + + for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++) + { // per superarc + // explicit test for root superarc / attachment points + if (vtkm::worklet::contourtree_augmented::NoSuchElement( + hierarchicalTreeSuperarcsPortal.Get(superarc))) + { + continue; + } + + // now retrieve the branch root's global ID + vtkm::Id branchRootSuperId = branchRootPortal.Get(superarc); + vtkm::Id branchRootRegularId = hierarchicalTreeSupernodesPortal.Get(branchRootSuperId); + vtkm::Id branchRootGlobalId = + hierarchicalTreeRegularNodeGlobalIdsPortal.Get(branchRootRegularId); + + // now retrieve the global ID & value for each end & output them + vtkm::Id superFromRegularId = hierarchicalTreeSupernodesPortal.Get(superarc); + vtkm::Id superFromGlobalId = hierarchicalTreeRegularNodeGlobalIdsPortal.Get(superFromRegularId); + FieldType superFromValue = hierarchicalTreeDataValuesPortal.Get(superFromRegularId); + resultStream << branchRootGlobalId << "\t" << superFromValue << "\t" << superFromGlobalId + << std::endl; + + // now retrieve the global ID & value for each end & output them + vtkm::Id superToRegularId = vtkm::worklet::contourtree_augmented::MaskedIndex( + hierarchicalTreeSuperarcsPortal.Get(superarc)); + vtkm::Id superToGlobalId = hierarchicalTreeRegularNodeGlobalIdsPortal.Get(superToRegularId); + FieldType superToValue = hierarchicalTreeDataValuesPortal.Get(superToRegularId); + resultStream << branchRootGlobalId << "\t" << superToValue << "\t" << superToGlobalId + << std::endl; + } // per superarc + return resultStream.str(); +} // PrintBranches + +// debug routine +template +std::string HierarchicalVolumetricBranchDecomposer::DebugPrint(std::string message, + const char* fileName, + long lineNum) +{ // DebugPrint() + std::stringstream resultStream; + resultStream << "----------------------------------------" << std::endl; + resultStream << std::setw(30) << std::left << fileName << ":" << std::right << std::setw(4) + << lineNum << std::endl; + resultStream << std::left << message << std::endl; + resultStream << "Hypersweep Value Array Contains: " << std::endl; + resultStream << "----------------------------------------" << std::endl; + resultStream << std::endl; + + vtkm::worklet::contourtree_augmented::PrintHeader(this->UpVolume.GetNumberOfValues(), + resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Up Volume by SA", this->UpVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Down Volume by SA", this->DownVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Best Down Snode by SN", this->BestDownSupernode, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Best Down Volume by SN", this->BestDownVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Best Up Snode by SN", this->BestUpSupernode, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Best Up Volume by SN", this->BestUpVolume, -1, resultStream); + std::cout << std::endl; + return resultStream.str(); +} // DebugPrint() + +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + + +#endif diff --git a/vtkm/worklet/contourtree_distributed/HyperSweepBlock.h b/vtkm/worklet/contourtree_distributed/HyperSweepBlock.h index 5d4f1cd7b..8d5948dd6 100644 --- a/vtkm/worklet/contourtree_distributed/HyperSweepBlock.h +++ b/vtkm/worklet/contourtree_distributed/HyperSweepBlock.h @@ -68,7 +68,7 @@ template struct HyperSweepBlock { HyperSweepBlock( - const vtkm::Id& localBlockNo, + const vtkm::Id localBlockNo, const int globalBlockId, const vtkm::Id3& origin, const vtkm::Id3& size, @@ -100,7 +100,7 @@ struct HyperSweepBlock vtkm::cont::ArrayHandle DependentVolume; // Destroy function allowing DIY to own blocks and clean them up after use - static void destroy(void* b) + static void Destroy(void* b) { delete static_cast*>(b); } diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/CMakeLists.txt b/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/CMakeLists.txt index f7fce9bf8..417f732a9 100644 --- a/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/CMakeLists.txt +++ b/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/CMakeLists.txt @@ -11,6 +11,7 @@ set(headers FindRegularByGlobal.h FindSuperArcForUnknownNode.h + FindSuperArcBetweenNodes.h PermuteComparator.h InitalizeSuperchildrenWorklet.h ) diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcBetweenNodes.h b/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcBetweenNodes.h new file mode 100644 index 000000000..b20358601 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcBetweenNodes.h @@ -0,0 +1,142 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// +// This code is an extension of the algorithm presented in the paper: +// Parallel Peak Pruning for Scalable SMP Contour Tree Computation. +// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens. +// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization +// (LDAV), October 2016, Baltimore, Maryland. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_find_superarc_between_nodes_h +#define vtk_m_worklet_contourtree_distributed_find_superarc_between_nodes_h + +#include +#include + + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ + + +/// Device implementation of FindSuperArcBetweenNodes for the HierarchicalContourTree +/// Used in the hierarchical branch decomposition +class FindSuperArcBetweenNodesDeviceData +{ +public: + using IndicesPortalType = + typename vtkm::worklet::contourtree_augmented::IdArrayType::ReadPortalType; + + VTKM_CONT + FindSuperArcBetweenNodesDeviceData( + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token, + const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs) + { + // Prepare the arrays for input and store the array portals + // so that they can be used inside a workelt + this->SuperarcsPortal = superarcs.PrepareForInput(device, token); + } + + // routine to find the superarc from one node to another + // it will always be the same ID as one of them if it exists + // if not, it will be NO_SUCH_ELEMENT + VTKM_EXEC + vtkm::Id FindSuperArcBetweenNodes(vtkm::Id firstSupernode, vtkm::Id secondSupernode) const + { // FindSuperArcBetweenNodes() + // if the second is the target of the first's superarc + if (vtkm::worklet::contourtree_augmented::MaskedIndex( + this->SuperarcsPortal.Get(firstSupernode)) == secondSupernode) + { + return firstSupernode; + } + // flip and test the other way + if (vtkm::worklet::contourtree_augmented::MaskedIndex( + this->SuperarcsPortal.Get(secondSupernode)) == firstSupernode) + { + return secondSupernode; + } + // otherwise it fails + return vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + } // FindSuperArcBetweenNodes() + +private: + // Array portals needed by FindSuperArcBetweenNodes + IndicesPortalType SuperarcsPortal; +}; + + +/// ExecutionObject to generate a device object to use FindSuperArcBetweenNodes for the HierarchicalContourTree +class FindSuperArcBetweenNodes : public vtkm::cont::ExecutionObjectBase +{ +public: + /// constructor + VTKM_CONT + FindSuperArcBetweenNodes(const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs) + : Superarcs(superarcs) + { + } + + VTKM_CONT FindSuperArcBetweenNodesDeviceData + PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const + { + return FindSuperArcBetweenNodesDeviceData(device, token, this->Superarcs); + } + +private: + // Array portals needed by FindSuperArcBetweenNodes + vtkm::worklet::contourtree_augmented::IdArrayType Superarcs; +}; + + +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CMakeLists.txt b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CMakeLists.txt new file mode 100644 index 000000000..a7228d0ec --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CMakeLists.txt @@ -0,0 +1,21 @@ +##============================================================================ +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +##============================================================================ + +set(headers + CollapseBranchesPointerDoublingWorklet.h + CollapseBranchesWorklet.h + FindBestSupernodeWorklet.h + LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h + LocalBestUpDownByVolumeInitSuperarcListWorklet.h + LocalBestUpDownByVolumeWorklet.h + SuperArcVolumetricComparatorIndirectGlobalIdComparator.h +) + +vtkm_declare_headers(${headers}) diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesPointerDoublingWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesPointerDoublingWorklet.h new file mode 100644 index 000000000..05489faf3 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesPointerDoublingWorklet.h @@ -0,0 +1,94 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_collapse_branches_pointer_doubling_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_collapse_branches_pointer_doubling_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +class CollapseBranchesPointerDoublingWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + using ControlSignature = void(WholeArrayInOut branchRoot); + using ExecutionSignature = void(InputIndex, _1); + using InputDomain = _1; + + /// Default Constructor + VTKM_EXEC_CONT + CollapseBranchesPointerDoublingWorklet() {} + + /// operator() of the workelt + template + VTKM_EXEC void operator()(const vtkm::Id superarc, + const InOutFieldPortalType& branchRootPortal) const + { + vtkm::Id branchRootVal1 = branchRootPortal.Get(superarc); + vtkm::Id branchRootVal2 = branchRootPortal.Get(branchRootPortal.Get(superarc)); + if (branchRootVal1 != branchRootVal2) + { + branchRootPortal.Set(superarc, branchRootVal2); + } + } // operator()() + + +}; // CollapseBranchesPointerDoublingWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h new file mode 100644 index 000000000..7be4e3d7b --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h @@ -0,0 +1,244 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_collapse_branches_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_collapse_branches_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +class CollapseBranchesWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + using ControlSignature = void( + FieldIn bestUpSupernode, + FieldIn bestDownSupernode, + // Execution objects from the hierarchical tree to use the FindRegularByGlobal function + ExecObject findRegularByGlobal, + // Execution objects from the hierarchical tree to use the FindSuperArcBetweenNodes, function + ExecObject findSuperArcBetweenNodes, + WholeArrayIn hierarchicalTreeRegular2supernode, + WholeArrayIn hierarchicalTreeWhichRound, + WholeArrayInOut branchRoot); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7); + using InputDomain = _1; + + /// Default Constructor + VTKM_EXEC_CONT + CollapseBranchesWorklet() {} + + /// operator() of the workelt + template + VTKM_EXEC void operator()( + const vtkm::Id& supernode, // iteration index + const vtkm::Id& bestUpSupernodeId, // bestUpSupernode[supernode] + const vtkm::Id& bestDownSupernodeId, // bestDownSupernode[supernode] + const ExecObjectType1& findRegularByGlobal, // Execution object to call FindRegularByGlobal + const ExecObjectType2& + findSuperArcBetweenNodes, // Execution object to call FindSuperArcBetweenNodes + const InFieldPortalType& hierarchicalTreeRegular2supernodePortal, + const InFieldPortalType& hierarchicalTreeWhichRoundPortal, + const InOutFieldPortalType& branchRootPortal) const + { + // per supernode + // For each supernode, convert the best up into a superarc ID + // Note that the superarc may not belong to this rank, and that the superarc might be oriented either direction + // So we search for the best up global ID in the hierarchical tree + // If it does not exist, then this superarc does not belong to the rank, and can be ignored + // If it does exist and is a downwards superarc, we now have the correct ID + // If it does exist and is an upwards superarc, then the current supernode must have an ascending arc to it, and we're done + // Also do the same for the best down, then for each supernode, point the higher numbered at the lower + + // if there is no best up, we're at an upper leaf and will not connect up two superarcs anyway, so we can skip the supernode + if (vtkm::worklet::contourtree_augmented::NoSuchElement(bestUpSupernodeId)) + { + return; + } + + // Search for the regular ID of the best up supernode + vtkm::Id bestUpLocalRegularId = findRegularByGlobal.FindRegularByGlobal(bestUpSupernodeId); + + // test to see whether it exists in this rank's hierarchical tree. + if (vtkm::worklet::contourtree_augmented::NoSuchElement(bestUpLocalRegularId)) + { + return; + } + + // do the same for the best down + // Search for the regular ID of the best down supernode + vtkm::Id bestDownLocalRegularId = findRegularByGlobal.FindRegularByGlobal(bestDownSupernodeId); + + // test to see whether it exists in this rank's hierarchical tree. + if (vtkm::worklet::contourtree_augmented::NoSuchElement(bestDownLocalRegularId)) + { + return; + } + + // Convert regular to super ID + vtkm::Id bestUpLocalSupernodeId = + hierarchicalTreeRegular2supernodePortal.Get(bestUpLocalRegularId); + vtkm::Id bestDownLocalSupernodeId = + hierarchicalTreeRegular2supernodePortal.Get(bestDownLocalRegularId); + + // local variable for the superarc IDs + vtkm::Id bestUpSuperarc = + findSuperArcBetweenNodes.FindSuperArcBetweenNodes(bestUpLocalSupernodeId, supernode); + vtkm::Id bestDownSuperarc = + findSuperArcBetweenNodes.FindSuperArcBetweenNodes(bestDownLocalSupernodeId, supernode); + + // right: we now know the local IDs of both. Take the more junior and point it at the more senior - i.e. always orient inbound + // at the root supernode, the virtual root superarc will not be used, so we compare round/iteration/ID of the two superarcs anyway + // WARNING: it might appear that there is potential for loops in the algorithm &/or write collisions + // but there isn't because our superarcs are *ALWAYS* oriented inwards, as long as the test is correct ;-> + + // so to find seniority, &c., we retrieve round number. + // we don't need the iteration number, because a higher iteration (more senior) always has a higher ID for the same round + vtkm::Id bestUpRound = hierarchicalTreeWhichRoundPortal.Get(bestUpSuperarc); + vtkm::Id bestDownRound = hierarchicalTreeWhichRoundPortal.Get(bestDownSuperarc); + + // more senior rounds are higher numbered + if ((bestUpRound > bestDownRound) || + // within each round, more senior iterations are higher numbered + ((bestUpRound == bestDownRound) && (bestUpSuperarc > bestDownSuperarc))) + { // up is more senior + branchRootPortal.Set(bestDownSuperarc, bestUpSuperarc); + } // up is more senior + else // all other cases, go the opposite way - NB: assumes we will never have the same superarc twice + { // down is more senior + branchRootPortal.Set(bestUpSuperarc, bestDownSuperarc); + } // down is more senior + + /* + // This worklet implements the following loop. + for (vtkm::Id supernode = 0; supernode < hierarchicalTree.supernodes.size(); supernode++) + { // per supernode + // For each supernode, convert the best up into a superarc ID + // Note that the superarc may not belong to this rank, and that the superarc might be oriented either direction + // So we search for the best up global ID in the hierarchical tree + // If it does not exist, then this superarc does not belong to the rank, and can be ignored + // If it does exist and is a downwards superarc, we now have the correct ID + // If it does exist and is an upwards superarc, then the current supernode must have an ascending arc to it, and we're done + // Also do the same for the best down, then for each supernode, point the higher numbered at the lower + + vtkm::Id bestUpSupernodeID = bestUpSupernode[supernode]; + + // if there is no best up, we're at an upper leaf and will not connect up two superarcs anyway, so we can skip the supernode + if (noSuchElement(bestUpSupernodeID)) + continue; + + // Search for the regular ID of the best up supernode + vtkm::Id bestUpLocalRegularID = hierarchicalTree.FindRegularByGlobal(bestUpSupernodeID); + + // test to see whether it exists in this rank's hierarchical tree. + if (noSuchElement(bestUpLocalRegularID)) + continue; + + // do the same for the best down + vtkm::Id bestDownSupernodeID = bestDownSupernode[supernode]; + + // Search for the regular ID of the best down supernode + vtkm::Id bestDownLocalRegularID = hierarchicalTree.FindRegularByGlobal(bestDownSupernodeID); + + // test to see whether it exists in this rank's hierarchical tree. + if (noSuchElement(bestDownLocalRegularID)) + continue; + + // Convert regular to super ID + vtkm::Id bestUpLocalSupernodeID = hierarchicalTree.regular2supernode[bestUpLocalRegularID]; + + vtkm::Id bestDownLocalSupernodeID = hierarchicalTree.regular2supernode[bestDownLocalRegularID]; + + // local variable for the superarc IDs + vtkm::Id bestUpSuperarc = hierarchicalTree.FindSuperArcBetweenNodes(bestUpLocalSupernodeID, supernode); + + vtkm::Id bestDownSuperarc = hierarchicalTree.FindSuperArcBetweenNodes(bestDownLocalSupernodeID, supernode); + + // right: we now know the local IDs of both. Take the more junior and point it at the more senior - i.e. always orient inbound + // at the root supernode, the virtual root superarc will not be used, so we compare round/iteration/ID of the two superarcs anyway + // WARNING: it might appear that there is potential for loops in the algorithm &/or write collisions + // but there isn't because our superarcs are *ALWAYS* oriented inwards, as long as the test is correct ;-> + + // so to find seniority, &c., we retrieve round number. + // we don't need the iteration number, because a higher iteration (more senior) always has a higher ID for the same round + vtkm::Id bestUpRound = hierarchicalTree.whichRound[bestUpSuperarc]; + vtkm::Id bestDownRound = hierarchicalTree.whichRound[bestDownSuperarc]; + + // more senior rounds are higher numbered + if ( (bestUpRound > bestDownRound) || + // within each round, more senior iterations are higher numbered + ( (bestUpRound == bestDownRound) && (bestUpSuperarc > bestDownSuperarc) )) + { // up is more senior + branchRoot[bestDownSuperarc] = bestUpSuperarc; + } // up is more senior + else // all other cases, go the opposite way - NB: assumes we will never have the same superarc twice + { // down is more senior + branchRoot[bestUpSuperarc] = bestDownSuperarc; + } // down is more senior + } // per supernode + */ + } // operator()() + + +}; // CollapseBranchesWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/FindBestSupernodeWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/FindBestSupernodeWorklet.h new file mode 100644 index 000000000..77d51df04 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/FindBestSupernodeWorklet.h @@ -0,0 +1,97 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_combine_branch_decomposers_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_combine_branch_decomposers_worklet_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +// Note: We use a template bool paramter for the tiebreak role choice instead of a variable +// so that the if statement can be optimized at compile time and we have fewer tests inside +// the worklet. +template +class FindBestSupernodeWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn incomingBestVolume, + FieldIn incomingBestSupernode, + FieldInOut bestVolume, + FieldInOut bestSupernode); + using ExecutionSignature = void(_1, _2, _3, _4); + + VTKM_EXEC void operator()(vtkm::Id incomingBestVolume, + vtkm::Id incomingBestSupernode, + vtkm::Id& bestVolume, + vtkm::Id& bestSupernode) const + { + // this is the same test as the SuperArcVolumetricComparator, hard-coded since we're not + // dealing with an array here + if ((incomingBestVolume > bestVolume) || + // don't forget the tiebreak rule + ((incomingBestVolume == bestVolume) && + (tieBreakGreaterThan ? (incomingBestSupernode > bestSupernode) + : (incomingBestSupernode < bestSupernode)))) + { // new one is better + bestVolume = incomingBestVolume; + bestSupernode = incomingBestSupernode; + } // new one is better + } // operator()() +}; // FindBestSupernodeWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h new file mode 100644 index 000000000..00483d9fe --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h @@ -0,0 +1,154 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_localbestupdownbyvolume_best_up_down_edge_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_localbestupdownbyvolume_best_up_down_edge_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +/// Worklet used in HierarchicalAugmenter::CopyBaseRegularStructure for +/// finding the superparent for each node needed +class LocalBestUpDownByVolumeBestUpDownEdgeWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + /// NOTE: We require the input arrays (aside form the input domain) to be permutted by the + /// regularNodesNeeded input domain so that we can use FieldIn instead of WholeArrayIn + /// NOTE: We require ArrayHandleView for the output arrays of the range [numExistingRegular:end] so + /// that we can use FieldOut instead of requiring WholeArrayInOut + using ControlSignature = void( + FieldIn + permutedHierarchicalTreeSuperarcs, // hierarchicalTree.Superarcs permuted by actualSuperarcs + FieldIn permutedDependetValues, // dependentValues permuted by actualSuperarcs + FieldIn permutedIntrinsicValues, // intrinsicValues permuted by actualSuperarcs + FieldOut permutedUpVolume, // upVolume permuted by actualSuperarcs + FieldOut permitedDownVolume // downVolume permited by actualSuperarcs + ); + using ExecutionSignature = void(_1, _2, _3, _4, _5); + using InputDomain = _1; + + /// Default Constructor + VTKM_EXEC_CONT + LocalBestUpDownByVolumeBestUpDownEdgeWorklet(const vtkm::Id totalVolume) + : TotalVolume(totalVolume) + { + } + + /// operator() of the workelt + template + VTKM_EXEC void operator()( + const vtkm::Id& + hierarchicalTreeSuperarc, // hierarchicalTree.superarcs[actualSuperarcs[InputIndex]] + const FieldType& dependentValue, // dependentValues[actualSuperarcs[InputIndex]] + const FieldType& intrinsicValue, // intrinsicValues[actualSuperarcs[InputIndex]] + vtkm::Id& upVolume, // upVolume[actualSuperarcs[InputIndex]] + vtkm::Id& downVolume // downVolume[actualSuperarcs[InputIndex]] + ) const + { + // per actual superarc + // retrieve the superarc Id (down via the superarcId FieldIn parameter) + if (vtkm::worklet::contourtree_augmented::IsAscending(hierarchicalTreeSuperarc)) + { // ascending superarc + upVolume = dependentValue; + // at the inner end, dependent volume is the total in the subtree. Then there are vertices along the edge + // itself (intrinsic volume), including the supernode at the outer end + // So, to get the "dependent" volume in the other direction, we start + // with totalVolume - dependent, then subtract (intrinsic - 1) + downVolume = (this->TotalVolume - dependentValue) + (intrinsicValue - 1); + } // ascending superarc + else + { // descending superarc + downVolume = dependentValue; + // at the inner end, dependent volume is the total in the subtree. Then there are vertices along the edge + // itself (intrinsic volume), including the supernode at the outer end + // So, to get the "dependent" volume in the other direction, we start + // with totalVolume - dependent, then subtract (intrinsic - 1) + upVolume = (this->TotalVolume - dependentValue) + (intrinsicValue - 1); + } // descending superarc + + /* // This worklet implements the follwing loop + for (vtkm::Id actualSuperarc = 0; actualSuperarc < nActualSuperarcs; actualSuperarc++) + { // per actual superarc + // retrieve the superarc ID + vtkm::Id superarcID = actualSuperarcs[actualSuperarc]; + if (isAscending(hierarchicalTree.superarcs[superarcID])) + { // ascending superarc + upVolume[superarcID] = dependentValues[superarcID]; + // at the inner end, dependent volume is the total in the subtree. Then there are vertices along the edge itself (intrinsic volume), including the supernode at the outer end + // So, to get the "dependent" volume in the other direction, we start with totalVolume - dependent, then subtract (intrinsic - 1) + downVolume[superarcID] = (totalVolume - dependentValues[superarcID]) + (intrinsicValues[superarcID] - 1); + } // ascending superarc + else + { // descending superarc + downVolume[superarcID] = dependentValues[superarcID]; + // at the inner end, dependent volume is the total in the subtree. Then there are vertices along the edge itself (intrinsic volume), including the supernode at the outer end + // So, to get the "dependent" volume in the other direction, we start with totalVolume - dependent, then subtract (intrinsic - 1) + upVolume[superarcID] = (totalVolume - dependentValues[superarcID]) + (intrinsicValues[superarcID] - 1); + } // descending superarc + } // per superarc + */ + } // operator()() + +private: + vtkm::Id TotalVolume; + +}; // LocalBestUpDownByVolumeBestUpDownEdgeWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeInitSuperarcListWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeInitSuperarcListWorklet.h new file mode 100644 index 000000000..7b42f3145 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeInitSuperarcListWorklet.h @@ -0,0 +1,104 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_localbestupdownbyvolume_init_superarc_list_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_localbestupdownbyvolume_init_superarc_list_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +class LocalBestUpDownByVolumeInitSuperarcListWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn hierarchicalTreeSuperarcs, FieldOut superarcList); + using ExecutionSignature = _2(InputIndex, _1); + using InputDomain = _1; + + VTKM_EXEC vtkm::worklet::contourtree_augmented::EdgePair operator()( + vtkm::Id superarc, // InputIndex in [0, hierarchicalTree.Superarcs.GetNumberOfValues() - 1] + vtkm::Id hierarchicalTreeSuperarc // hierarchicalTree.Superarcs[superarc] + ) const + { + using vtkm::worklet::contourtree_augmented::EdgePair; + using vtkm::worklet::contourtree_augmented::IsAscending; + using vtkm::worklet::contourtree_augmented::MaskedIndex; + + if (IsAscending(hierarchicalTreeSuperarc)) + { + return EdgePair(superarc, MaskedIndex(hierarchicalTreeSuperarc)); + } + else + { + return EdgePair(MaskedIndex(hierarchicalTreeSuperarc), superarc); + } + + /* // This worklet implements the follwing loop + for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++) + { // per superarc + if (isAscending(hierarchicalTree.superarcs[superarc])) + superarcList[superarc] = Edge(superarc, maskedIndex(hierarchicalTree.superarcs[superarc])); + else + superarcList[superarc] = Edge(maskedIndex(hierarchicalTree.superarcs[superarc]), superarc); + } // per superarc + */ + } // operator()() + +}; // LocalBestUpDownByVolumeInitSuperarcListWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeWorklet.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeWorklet.h new file mode 100644 index 000000000..54ff73f68 --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeWorklet.h @@ -0,0 +1,229 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_local_best_updown_by_volume_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_local_best_updown_by_volume_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + +/// Worklet used in HierarchicalAugmenter::CopyBaseRegularStructure for +/// finding the superparent for each node needed + +/// Template parameter is bool indicating whether we are processing up- or +/// down best volumes and corresponding whether we need to use the low or +/// high end of the edge. Note: We make this a template paramter so that +/// the corresponding if statement can already be optimozed away during +/// compile time. +template +class LocalBestUpDownByVolumeWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + using ControlSignature = void( + WholeArrayIn actualSuperacrs, + WholeArrayIn superarcList, // superarc list + FieldIn + permutedUpDownVolume, // upVolumne if IsDown==True, or downVolume if IsDown==False. These are swapped as IsDown refers to the output arrays + WholeArrayIn hierarchicalTreeRegularNodeGlobalIds, + WholeArrayIn hierarchicalTreeSupernodes, + WholeArrayInOut + bestUpDownSupernode, // bestUpSupernode if IsDown==False, or bestDownSupernode if IsDown==True + WholeArrayInOut + bestUpDownVolume // bestUpVolume if IsDown==False, or bestDownVolume if IsDown==True + ); + // TODO: Check if we need WholeArrayInOut here for the output arrays or if WholeArrayOut is sufficient + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7); + using InputDomain = _1; + + /// Default Constructor + VTKM_EXEC_CONT + LocalBestUpDownByVolumeWorklet(const vtkm::Id numActualSuperarcs) + : NumberActualSuperarcs(numActualSuperarcs) + { + } + + /// operator() of the workelt + template + VTKM_EXEC void operator()(const vtkm::Id& actualSuperarcIndex, + const InFieldPortalType1& actualSuperarcsPortal, + const InFieldPortalType2& superarcListPortal, + const vtkm::Id& upDownVolumeValue, // upDownVolume[superarcID] + const InFieldPortalType3& hierarchicalTreeRegularNodeGlobalIdsPortal, + const InFieldPortalType4& hierarchicalTreeSupernodesPortal, + const OutFieldPortalType1& bestUpDownSupernodePortal, + const OutFieldPortalType2& bestUpDownVolumePortal) const + { + // per actual superarc + vtkm::Id superarcId = actualSuperarcsPortal.Get(actualSuperarcIndex); + const vtkm::worklet::contourtree_augmented::EdgePair& edge = superarcListPortal.Get(superarcId); + + if (IsDown) + { + // if it's the last one + if (actualSuperarcIndex == NumberActualSuperarcs - 1) + { // last in array + bestUpDownSupernodePortal.Set(edge.second, + hierarchicalTreeRegularNodeGlobalIdsPortal.Get( + hierarchicalTreeSupernodesPortal.Get(edge.first))); + bestUpDownVolumePortal.Set(edge.second, upDownVolumeValue); + } // last in array + else + { // not the last one + const vtkm::worklet::contourtree_augmented::EdgePair& nextEdge = + superarcListPortal.Get(actualSuperarcsPortal.Get(actualSuperarcIndex + 1)); + // if the next edge belongs to another, we're the highest + if (nextEdge.second != edge.second) + { // last in group + bestUpDownSupernodePortal.Set(edge.second, + hierarchicalTreeRegularNodeGlobalIdsPortal.Get( + hierarchicalTreeSupernodesPortal.Get(edge.first))); + bestUpDownVolumePortal.Set(edge.second, upDownVolumeValue); + } // last in group + } // not the last one + } // if(this->IsDown) + else // Processing the Up arrays. This is essentiall the same, but we need to use the lower end of the edge instead + { + // if it's the last one + if (actualSuperarcIndex == NumberActualSuperarcs - 1) + { // last in array + bestUpDownSupernodePortal.Set(edge.first, + hierarchicalTreeRegularNodeGlobalIdsPortal.Get( + hierarchicalTreeSupernodesPortal.Get(edge.second))); + bestUpDownVolumePortal.Set(edge.first, upDownVolumeValue); + } // last in array + else + { // not the last one + const vtkm::worklet::contourtree_augmented::EdgePair& nextEdge = + superarcListPortal.Get(actualSuperarcsPortal.Get(actualSuperarcIndex + 1)); + // if the next edge belongs to another, we're the highest + if (nextEdge.first != edge.first) + { // best in group + bestUpDownSupernodePortal.Set(edge.first, + hierarchicalTreeRegularNodeGlobalIdsPortal.Get( + hierarchicalTreeSupernodesPortal.Get(edge.second))); + bestUpDownVolumePortal.Set(edge.first, upDownVolumeValue); + } // best in group + } // not the last one + } // else (if(this->isDown)) + /* + // This worklet implements the following loop. Depending on whether we are working with the Up- or DownVolumes + // This function implements one of the following logic + // II B 2. Per vertex, best superarc writes to the best downward array + for (vtkm::Id actualSuperarc = 0; actualSuperarc < nActualSuperarcs; actualSuperarc++) + { // per actual superarc + vtkm::Id superarcID = actualSuperarcs[actualSuperarc]; + Edge &edge = superarcList[superarcID]; + // if it's the last one + if (actualSuperarc == nActualSuperarcs-1) + { // last in array + bestDownSupernode[edge.high] = hierarchicalTree.regularNodeGlobalIDs[hierarchicalTree.supernodes[edge.low]]; + bestDownVolume[edge.high] = upVolume[superarcID]; + } // last in array + else + { // not the last one + Edge &nextEdge = superarcList[actualSuperarcs[actualSuperarc+1]]; + // if the next edge belongs to another, we're the highest + if (nextEdge.high != edge.high) + { // last in group + bestDownSupernode[edge.high] = hierarchicalTree.regularNodeGlobalIDs[hierarchicalTree.supernodes[edge.low]]; + bestDownVolume[edge.high] = upVolume[superarcID]; + } // last in group + } // not the last one + } // per actual superarc + */ + /* // Or in the Up case we have. THe main difference is the we either use the Up or Down Volume and + // and we following here the low end of the edge and top end of the edge in the other case + + for (vtkm::Id actualSuperarc = 0; actualSuperarc < nActualSuperarcs; actualSuperarc++) + { // per actual superarc + vtkm::Id superarcID = actualSuperarcs[actualSuperarc]; + Edge &edge = superarcList[superarcID]; + // if it's the last one + if (actualSuperarc == nActualSuperarcs-1) + { // last in array + bestUpSupernode[edge.low] = hierarchicalTree.regularNodeGlobalIDs[hierarchicalTree.supernodes[edge.high]]; + bestUpVolume[edge.low] = downVolume[superarcID]; + } // last in array + else + { // not the last one + Edge &nextEdge = superarcList[actualSuperarcs[actualSuperarc+1]]; + // if the next edge belongs to another, we're the highest + if (nextEdge.low != edge.low) + { // best in group + bestUpSupernode[edge.low] = hierarchicalTree.regularNodeGlobalIDs[hierarchicalTree.supernodes[edge.high]]; + bestUpVolume[edge.low] = downVolume[superarcID]; + } // best in group + } // not the last one + + } // per actual superarc + */ + } // operator()() + +private: + vtkm::Id NumberActualSuperarcs; + +}; // LocalBestUpDownByVolumeWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/SuperArcVolumetricComparatorIndirectGlobalIdComparator.h b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/SuperArcVolumetricComparatorIndirectGlobalIdComparator.h new file mode 100644 index 000000000..8d92d653b --- /dev/null +++ b/vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/SuperArcVolumetricComparatorIndirectGlobalIdComparator.h @@ -0,0 +1,247 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// +// SuperArcVolumetricComparatorIndirectGlobalID.h - a comparator for sorting superarcs by volume +// Has to take a flag for high end vs. low end sorting +// Also, this version takes supernode IDs rather than global IDs, so has an extra indirection +// +//======================================================================================= +// +// COMMENTS: +// +// A comparator that sorts superarc pairs by: +// 1. ID of low end vertex +// 2. volumetric measure at low end +// 3. global index of upper end, OR +// +// the same for the higher end. +// +// Notice that 2. only applies if two edges share a lower end and have the same volume. +// We then look at the index at the upper end to see which is "furthest" from the low end +// +//======================================================================================= + + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_superarc_volumetric_comparator_indirect_globalid_comparator_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_volumetric_branch_decomposer_superarc_volumetric_comparator_indirect_globalid_comparator_h + +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_volumetric_branch_decomposer +{ + + +/// Implementation of the comparator for the SuperArcVolumetricComparatorIndirectGlobalId ExecutionObject +class SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl +{ +public: + using IdArrayPortalType = + typename vtkm::worklet::contourtree_augmented::IdArrayType::ReadPortalType; + using EdgePairArrayPortalType = + typename vtkm::worklet::contourtree_augmented::EdgePairArray::ReadPortalType; + + // constructor + VTKM_CONT + SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl( + IdArrayPortalType weightPortal, + EdgePairArrayPortalType superarcListPortal, + IdArrayPortalType globalIdPortal, + bool pairsAtLowEnd) + : WeightPortal(weightPortal) + , SuperarcListPortal(superarcListPortal) + , GlobalIdPortal(globalIdPortal) + , PairsAtLowEnd(pairsAtLowEnd) + { // constructor + } // constructor + + // () operator - gets called to do comparison + VTKM_EXEC + bool operator()(const vtkm::Id& left, const vtkm::Id& right) const + { // operator() + // get local references to the edge details + vtkm::worklet::contourtree_augmented::EdgePair edgeLeft = this->SuperarcListPortal.Get(left); + vtkm::worklet::contourtree_augmented::EdgePair edgeRight = this->SuperarcListPortal.Get(right); + + if (this->PairsAtLowEnd) + { // pairs at low end + // test by low end ID + if (edgeLeft.first < edgeRight.first) + { + return true; + } + if (edgeLeft.first > edgeRight.first) + { + return false; + } + + // test by volumetric measure + vtkm::Id weightLeft = this->WeightPortal.Get(left); + vtkm::Id weightRight = this->WeightPortal.Get(right); + if (weightLeft < weightRight) + return true; + if (weightLeft > weightRight) + return false; + + // test by the global ID - we were past a supernode ID, so there's an + // extra level of indirection + vtkm::Id globalIdLeftEdgeSecond = this->GlobalIdPortal.Get(edgeLeft.second); + vtkm::Id globalIdRightEdgeSecond = this->GlobalIdPortal.Get(edgeRight.second); + if (globalIdLeftEdgeSecond < globalIdRightEdgeSecond) + { + return true; + } + if (globalIdLeftEdgeSecond > globalIdRightEdgeSecond) + { + return false; + } + + // fallback + return false; + } // pairs at low end + else + { // pairs at high end + // test by high end ID + if (edgeLeft.second < edgeRight.second) + { + return true; + } + if (edgeLeft.second > edgeRight.second) + { + return false; + } + + // test by volumetric measure + vtkm::Id weightLeft = this->WeightPortal.Get(left); + vtkm::Id weightRight = this->WeightPortal.Get(right); + if (weightLeft < weightRight) + { + return true; + } + if (weightLeft > weightRight) + { + return false; + } + + // test by the global ID - we were past a supernode ID, so there's an + // extra level of indirection + // Note the reversal from above - we want the greatest difference, not + // the greatest value + vtkm::Id globalIdLeftEdgeFirst = this->GlobalIdPortal.Get(edgeLeft.first); + vtkm::Id globalIdRightEdgeFirst = this->GlobalIdPortal.Get(edgeRight.first); + if (globalIdLeftEdgeFirst > globalIdRightEdgeFirst) + { + return true; + } + if (globalIdLeftEdgeFirst < globalIdRightEdgeFirst) + { + return false; + } + + // fallback + return false; + } // pairs at high end + } // operator() + +private: + IdArrayPortalType WeightPortal; + EdgePairArrayPortalType SuperarcListPortal; + IdArrayPortalType GlobalIdPortal; + bool PairsAtLowEnd; +}; // SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl + + +/// Execution object for Compartor used in HierarchicalVolumetricBranchDecomposer::LocalBestUpDownByVolume. +/// The comparator is used to sort superarc pairs by: +/// 1. ID of low end vertex +/// 2. volumetric measure at low end +/// 3. global index of upper end, OR +/// the same for the higher end. Notice that 2. only applies if two edges share a lower end and have the same volume. +/// We then look at the index at the upper end to see which is "furthest" from the low end +class SuperArcVolumetricComparatorIndirectGlobalIdComparator + : public vtkm::cont::ExecutionObjectBase +{ +public: + // constructor - takes vectors as parameters + VTKM_CONT + SuperArcVolumetricComparatorIndirectGlobalIdComparator( + const vtkm::worklet::contourtree_augmented::IdArrayType& weight, + const vtkm::worklet::contourtree_augmented::EdgePairArray& superarcList, + const vtkm::worklet::contourtree_augmented::IdArrayType& globalId, + bool pairsAtLowEnd) + : Weight(weight) + , SuperarcList(superarcList) + , GlobalId(globalId) + , PairsAtLowEnd(pairsAtLowEnd) + { // constructor + } // constructor + + /// Create a SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl object for use in the sort or worklet + VTKM_CONT SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl + PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const + { + return SuperArcVolumetricComparatorIndirectGlobalIdComparatorImpl( + this->Weight.PrepareForInput(device, token), + this->SuperarcList.PrepareForInput(device, token), + this->GlobalId.PrepareForInput(device, token), + this->PairsAtLowEnd); + } + +private: + vtkm::worklet::contourtree_augmented::IdArrayType Weight; + vtkm::worklet::contourtree_augmented::EdgePairArray SuperarcList; + vtkm::worklet::contourtree_augmented::IdArrayType GlobalId; + bool PairsAtLowEnd; +}; // SuperArcVolumetricComparatorIndirectGlobalIdComparator + +} // namespace hierarchical_volumetric_branch_decomposer +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif