Implement branch decomposition for hierarchical contour tree

Co-authored-by: Gunther H. Weber <GHWeber@lbl.gov>
This commit is contained in:
Oliver Ruebel 2022-03-01 16:24:37 -07:00 committed by Gunther H. Weber
parent af5073885e
commit 4fe495be8d
27 changed files with 2918 additions and 75 deletions

@ -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 <iostream>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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()

@ -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()

@ -76,6 +76,7 @@
#include <vtkm/worklet/contourtree_augmented/ProcessContourTree.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h>
#include <vtkm/worklet/contourtree_distributed/TreeCompiler.h>
// 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 <options> <fileName>" << std::endl;
std::cout << "ContourTreeDistributed <options> <fileName>" << std::endl;
std::cout << std::endl;
std::cout << "<fileName> 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<int>(rank)) + std::string("_Block_") +
std::to_string(static_cast<int>(ds_no)) + std::string(".txt");
std::string dumpVolumesString =
vtkm::worklet::contourtree_distributed::HierarchicalContourTree<ValueType>::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<int>(rank)) + std::string("_Block_") +
std::to_string(static_cast<int>(ds_no)) + std::string(".txt");
std::ofstream treeStream(volumesFileName.c_str());
treeStream << dumpVolumesString;
std::string dumpVolumesString =
vtkm::worklet::contourtree_distributed::HierarchicalContourTree<ValueType>::DumpVolumes(
supernodes,
superarcs,
regularNodeGlobalIds,
totalVolume,
intrinsicVolume,
dependentVolume);
std::string volumesFileName = std::string("TreeWithVolumes_Rank_") +
std::to_string(static_cast<int>(rank)) + std::string("_Block_") +
std::to_string(static_cast<int>(ds_no)) + std::string(".txt");
std::ofstream treeStream(volumesFileName.c_str());
treeStream << dumpVolumesString;
}
}
}
else

@ -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

@ -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"

@ -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 <typename FieldType>
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<vtkm::cont::DataSet> hierarchicalTreeOutputDataSet);
std::vector<vtkm::cont::DataSet>& hierarchicalTreeOutputDataSet);
template <typename FieldType>
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<vtkm::cont::DataSet>& 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;

@ -66,12 +66,15 @@
// distributed contour tree includes
#include <vtkm/worklet/contourtree_distributed/BoundaryTree.h>
#include <vtkm/worklet/contourtree_distributed/BoundaryTreeMaker.h>
#include <vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h>
#include <vtkm/worklet/contourtree_distributed/CombineHyperSweepBlockFunctor.h>
#include <vtkm/worklet/contourtree_distributed/ComputeDistributedBranchDecompositionFunctor.h>
#include <vtkm/worklet/contourtree_distributed/ComputeDistributedContourTreeFunctor.h>
#include <vtkm/worklet/contourtree_distributed/DistributedContourTreeBlockData.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalAugmenter.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalAugmenterFunctor.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalHyperSweeper.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h>
#include <vtkm/worklet/contourtree_distributed/HyperSweepBlock.h>
#include <vtkm/worklet/contourtree_distributed/InteriorForest.h>
#include <vtkm/worklet/contourtree_distributed/PrintGraph.h>
@ -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<std::size_t>(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 <typename FieldType>
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<vtkm::cont::DataSet> hierarchicalTreeOutputDataSet)
std::vector<vtkm::cont::DataSet>& 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<FieldType>;
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 <typename FieldType>
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<vtkm::cont::DataSet>& outputDataSet)
{
vtkm::cont::Timer timer;
timer.Start();
using BranchDecompositionBlock =
vtkm::worklet::contourtree_distributed::BranchDecompositionBlock<FieldType>;
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<FieldType>;
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<vtkm::Id> intrinsicVolume =
ds.GetField("IntrinsicVolume").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>();
vtkm::cont::ArrayHandle<vtkm::Id> dependentVolume =
ds.GetField("DependentVolume").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>();
// 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 <typename FieldType, typename StorageType, typename DerivedPolicy>
@ -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<vtkm::cont::DataSet> 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

@ -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);

@ -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 <vtkm/worklet/contourtree/Types.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalVolumetricBranchDecomposer.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_distributed
{
template <typename ContourTreeDataFieldType>
struct BranchDecompositionBlock
{
BranchDecompositionBlock(
vtkm::Id localBlockNo,
int globalBlockId,
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<ContourTreeDataFieldType>&
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<ContourTreeDataFieldType>&
HierarchicalContourTree;
vtkm::worklet::contourtree_distributed::HierarchicalVolumetricBranchDecomposer<
ContourTreeDataFieldType>
HierarchicalVolumetricBranchDecomposer;
vtkm::cont::ArrayHandle<vtkm::Id> BranchRoots;
// Destroy function allowing DIY to own blocks and clean them up after use
static void Destroy(void* b)
{
delete static_cast<BranchDecompositionBlock<ContourTreeDataFieldType>*>(b);
}
};
} // namespace contourtree_distributed
} // namespace worklet
} // namespace vtkm
#endif

@ -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)

@ -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<vtkm::Id> sendIntrinsicVolume;
vtkm::cont::Algorithm::Copy(intrinsicVolumeView, sendIntrinsicVolume);
vtkm::cont::ArrayCopy(intrinsicVolumeView, sendIntrinsicVolume);
vtkm::cont::ArrayHandle<vtkm::Id> sendDependentVolume;
vtkm::cont::Algorithm::Copy(dependentVolumeView, sendDependentVolume);
vtkm::cont::ArrayCopy(dependentVolumeView, sendDependentVolume);
// Send necessary data portions
rp.enqueue(target, sendIntrinsicVolume);

@ -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 <vtkm/Types.h>
#include <vtkm/worklet/contourtree_distributed/BranchDecompositionBlock.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/FindBestSupernodeWorklet.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/diy.h>
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 <typename ContourTreeDataFieldType>
struct ComputeDistributedBranchDecompositionFunctor
{
void operator()(
vtkm::worklet::contourtree_distributed::BranchDecompositionBlock<ContourTreeDataFieldType>* 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<int> 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<vtkm::Id> incomingBestUpVolume;
rp.dequeue(ingid, incomingBestUpVolume);
vtkm::cont::ArrayHandle<vtkm::Id> incomingBestUpSupernode;
rp.dequeue(ingid, incomingBestUpSupernode);
vtkm::cont::ArrayHandle<vtkm::Id> incomingBestDownVolume;
rp.dequeue(ingid, incomingBestDownVolume);
vtkm::cont::ArrayHandle<vtkm::Id> 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<true>{},
incomingBestUpVolume,
incomingBestUpSupernode,
bestUpVolumeView,
bestUpSupernodeView);
invoke(vtkm::worklet::contourtree_distributed::hierarchical_volumetric_branch_decomposer::
FindBestSupernodeWorklet<false>{},
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<vtkm::Id> sendBestUpVolume;
vtkm::cont::Algorithm::Copy(bestUpVolumeView, sendBestUpVolume);
vtkm::cont::ArrayHandle<vtkm::Id> sendBestUpSupernode;
vtkm::cont::Algorithm::Copy(bestUpSupernodeView, sendBestUpSupernode);
vtkm::cont::ArrayHandle<vtkm::Id> sendBestDownVolume;
vtkm::cont::Algorithm::Copy(bestDownVolumeView, sendBestDownVolume);
vtkm::cont::ArrayHandle<vtkm::Id> 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

@ -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<int>(rank)) + std::string("_Block_") +
std::to_string(static_cast<int>(block->BlockIndex)) + "_Round_" +
std::to_string(static_cast<int>(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<int>(block->BlockIndex)) + " Round " +
std::to_string(static_cast<int>(block->LocalBlockNo)) + " Round " +
std::to_string(rp.round()) + " Partner " + std::to_string(ingid) +
std::string(" Step 0 Combined Mesh");
std::string contourTreeMeshString =

@ -94,7 +94,7 @@ struct DistributedContourTreeBlockData
vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType> 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<DistributedContourTreeBlockData<FieldType>*>(b);
}

@ -75,6 +75,7 @@
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindRegularByGlobal.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcBetweenNodes.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcForUnknownNode.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/InitalizeSuperchildrenWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/PermuteComparator.h>
@ -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,

@ -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 <iomanip>
#include <string>
#include <vtkm/worklet/contourtree_augmented/NotNoSuchElementPredicate.h>
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_distributed/PrintGraph.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesPointerDoublingWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeBestUpDownEdgeWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeInitSuperarcListWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/LocalBestUpDownByVolumeWorklet.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_volumetric_branch_decomposer/SuperArcVolumetricComparatorIndirectGlobalIdComparator.h>
#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 <typename FieldType>
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<FieldType>(){};
/// 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<FieldType>*
hierarchicalTree,
const vtkm::cont::ArrayHandle<vtkm::Id>& intrinsicValues,
const vtkm::cont::ArrayHandle<vtkm::Id>& 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 <typename FieldType>
void HierarchicalVolumetricBranchDecomposer<FieldType>::LocalBestUpDownByVolume(
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>*
hierarchicalTree,
const vtkm::cont::ArrayHandle<vtkm::Id>& intrinsicValues,
const vtkm::cont::ArrayHandle<vtkm::Id>& 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<true>{ 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<false>{ 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 <typename FieldType>
void HierarchicalVolumetricBranchDecomposer<FieldType>::CollapseBranches(
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>*
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 <typename FieldType>
std::string HierarchicalVolumetricBranchDecomposer<FieldType>::PrintBranches(
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>*
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 <typename FieldType>
std::string HierarchicalVolumetricBranchDecomposer<FieldType>::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<vtkm::cont::ArrayHandle<vtkm::Id>>();
auto hierarchicalTreeSuperarcsPortal = hierarchicalTreeSuperarcsAH.ReadPortal();
vtkm::Id nSuperarcs = hierarchicalTreeSuperarcsAH.GetNumberOfValues();
auto hierarchicalTreeSupernodesPortal = ds.GetField("Supernodes")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
.ReadPortal();
auto hierarchicalTreeRegularNodeGlobalIdsPortal =
ds.GetField("RegularNodeGlobalIds")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
.ReadPortal();
auto hierarchicalTreeDataValuesPortal = ds.GetField("DataValues")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<FieldType>>()
.ReadPortal();
auto branchRootPortal = ds.GetField("BranchRoots")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>()
.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 <typename FieldType>
std::string HierarchicalVolumetricBranchDecomposer<FieldType>::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

@ -68,7 +68,7 @@ template <typename ContourTreeDataFieldType>
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<vtkm::Id> 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<HyperSweepBlock<ContourTreeDataFieldType>*>(b);
}

@ -11,6 +11,7 @@
set(headers
FindRegularByGlobal.h
FindSuperArcForUnknownNode.h
FindSuperArcBetweenNodes.h
PermuteComparator.h
InitalizeSuperchildrenWorklet.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 <vtkm/Types.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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

@ -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})

@ -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 <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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 <typename InOutFieldPortalType>
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

@ -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 <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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 <typename ExecObjectType1,
typename ExecObjectType2,
typename InFieldPortalType,
typename InOutFieldPortalType>
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

@ -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 <vtkm/worklet/WorkletMapField.h>
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 <bool tieBreakGreaterThan>
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

@ -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 <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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 <typename FieldType>
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

@ -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 <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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

@ -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 <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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 <bool IsDown>
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 <typename InFieldPortalType1,
typename InFieldPortalType2,
typename InFieldPortalType3,
typename InFieldPortalType4,
typename OutFieldPortalType1,
typename OutFieldPortalType2>
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

@ -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 <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
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<FieldType>::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