diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index b3ff96568..69b27b90d 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -220,6 +220,12 @@ int main(int argc, char* argv[]) } } + vtkm::Id presimplifyThreshold = 0; // Do not presimplify the hierachical contour tree + if (parser.hasOption("--presimplifyThreshold")) + { + presimplifyThreshold = std::stoi(parser.getOption("--presimplifyThreshold")); + } + bool useBoundaryExtremaOnly = true; if (parser.hasOption("--useFullBoundary")) { @@ -368,6 +374,10 @@ int main(int argc, char* argv[]) std::cout << "--eps= Floating point offset awary from the critical point. (default=0.00001)" << std::endl; + std::cout << "--presimplifyThreshold Integer volume threshold for presimplifying the tree" + << std::endl; + std::cout << " Default value is 0, indicating no presimplification" + << 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; @@ -411,6 +421,7 @@ int main(int argc, char* argv[]) << " augmentHierarchicalTree=" << augmentHierarchicalTree << std::endl << " computeVolumetricBranchDecomposition=" << computeHierarchicalVolumetricBranchDecomposition << std::endl + << " presimplifyThreshold=" << presimplifyThreshold << std::endl << " saveOutputData=" << saveOutputData << std::endl << " forwardSummary=" << forwardSummary << std::endl << " nblocks=" << numBlocks << std::endl @@ -648,6 +659,10 @@ int main(int argc, char* argv[]) filter.SetUseBoundaryExtremaOnly(useBoundaryExtremaOnly); filter.SetUseMarchingCubes(useMarchingCubes); filter.SetAugmentHierarchicalTree(augmentHierarchicalTree); + if (presimplifyThreshold > 0) + { + filter.SetPresimplifyThreshold(presimplifyThreshold); + } filter.SetSaveDotFiles(saveDotFiles); filter.SetActiveField("values"); @@ -923,6 +938,7 @@ int main(int argc, char* argv[]) close(save_err); } + std::cout << "DONE!!!" << std::endl; // Finalize and finish the execution MPI_Finalize(); return EXIT_SUCCESS; diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx index 9992fd64a..5832d9e98 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx +++ b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx @@ -169,6 +169,7 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed(vtkm::cont::LogLeve : UseBoundaryExtremaOnly(true) , UseMarchingCubes(false) , AugmentHierarchicalTree(false) + , PresimplifyThreshold(0) , SaveDotFiles(false) , TimingsLogLevel(timingsLogLevel) , TreeLogLevel(treeLogLevel) @@ -585,7 +586,10 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( vtkmdiy::RegularSwapPartners& partners, const FieldType&, // dummy parameter to get the type std::stringstream& timingsStream, - std::vector& hierarchicalTreeOutputDataSet) + const vtkm::cont::PartitionedDataSet& input, + bool useAugmentedTree, + std::vector>& intrinsicVolumes, + std::vector>& dependentVolumes) { // TODO/FIXME: CONSIDER MOVING CONTENTS OF THIS METHOD TO SEPARATE FILTER vtkm::cont::Timer timer; @@ -610,30 +614,35 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( inputContourTreeMaster.foreach ( [&](DistributedContourTreeBlockData* currInBlock, const vtkmdiy::Master::ProxyWithLink&) { vtkm::Id blockNo = currInBlock->LocalBlockNo; - const vtkm::cont::DataSet& currDS = hierarchicalTreeOutputDataSet[blockNo]; + //const vtkm::cont::DataSet& currDS = hierarchicalTreeOutputDataSet[blockNo]; + auto currOriginalBlock = input.GetPartition(static_cast(blockNo)); // The block size and origin may be modified during the FanIn so we need to use the // size and origin from the original decomposition instead of looking it up in the currInBlock vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart; - currDS.GetCellSet().CastAndCallForTypes( + currOriginalBlock.GetCellSet().CastAndCallForTypes( vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), pointDimensions, globalPointDimensions, globalPointIndexStart); - // 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. - hierarchical_hyper_sweep_master.add( - currInBlock->GlobalBlockId, - new HyperSweepBlock(blockNo, - currInBlock->GlobalBlockId, - globalPointIndexStart, - pointDimensions, - globalPointDimensions, - *currInBlock->HierarchicalAugmenter.AugmentedTree), - new vtkmdiy::Link()); + + // If we are pre-simplifying the tree then we need to use the base tree and if we compute the + // final volume, then we need to use the augmented tree + auto hierarchicalTreeToProcess = useAugmentedTree + ? currInBlock->HierarchicalAugmenter.AugmentedTree + : currInBlock->HierarchicalAugmenter.BaseTree; + hierarchical_hyper_sweep_master.add(currInBlock->GlobalBlockId, + new HyperSweepBlock(blockNo, + currInBlock->GlobalBlockId, + globalPointIndexStart, + pointDimensions, + globalPointDimensions, + *hierarchicalTreeToProcess), + new vtkmdiy::Link()); }); // Log time to copy the data to the HyperSweepBlock data objects @@ -675,7 +684,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler idRelabeler{ b->Origin, b->Size, b->GlobalSize }; - if (b->GlobalSize[2] <= 1) { vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation2DFreudenthal mesh( @@ -690,7 +698,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( hyperSweeper.InitializeIntrinsicVertexCount( b->HierarchicalContourTree, mesh, idRelabeler, b->IntrinsicVolume); } - #ifdef DEBUG_PRINT VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); VTKM_LOG_S(this->TreeLogLevel, @@ -719,7 +726,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( hyperSweeper.LocalHyperSweep(); #ifdef DEBUG_PRINT - VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); + VTKM_LOG_S(this->Tree LogLevel, "Block " << b->GlobalBlockId); VTKM_LOG_S(this->TreeLogLevel, b->HierarchicalContourTree.DebugPrint("After local hypersweep", __FILE__, __LINE__)); #endif @@ -755,16 +762,14 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); - // Print & add to output data set - //std::vector hierarchicalTreeAndVolumeOutputDataSet(localDataBlocks.size()); + // Add the intrinsic and dependent volumes to the output vectors + intrinsicVolumes.resize(inputContourTreeMaster.size()); + dependentVolumes.resize(inputContourTreeMaster.size()); hierarchical_hyper_sweep_master.foreach ( [&](HyperSweepBlock* b, const vtkmdiy::Master::ProxyWithLink&) { - vtkm::cont::Field intrinsicVolumeField( - "IntrinsicVolume", vtkm::cont::Field::Association::WholeDataSet, b->IntrinsicVolume); - hierarchicalTreeOutputDataSet[b->LocalBlockNo].AddField(intrinsicVolumeField); - vtkm::cont::Field dependentVolumeField( - "DependentVolume", vtkm::cont::Field::Association::WholeDataSet, b->DependentVolume); - hierarchicalTreeOutputDataSet[b->LocalBlockNo].AddField(dependentVolumeField); + intrinsicVolumes[b->LocalBlockNo] = b->IntrinsicVolume; + dependentVolumes[b->LocalBlockNo] = b->DependentVolume; + #ifdef DEBUG_PRINT VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId); VTKM_LOG_S( @@ -779,9 +784,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric( "Dependent Volume", b->DependentVolume, -1, volumeStream); 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; }); } @@ -1135,14 +1137,60 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); + + // Compute the volume for pre-simplification if we want to pre-simplify + // The dependent volumes from the unaugemented hierarchical tree are used for the pre-simplification + // as part of HierarchicalAugmenter.Initialize. + std::vector> unaugmentedDependentVolumes; + if (this->PresimplifyThreshold > 0) + { + // we don't need the unaugemented intrinsic Volumes for the pre-simplification, so we + // use a local variable that is being deleted automatically after the context + std::vector> unaugmentedIntrinsicVolumes; + // Compute the volume for the base hierarchical tree before augmentation in order to allow for pre-simplification. + this->ComputeVolumeMetric( + master, + assigner, + partners, + FieldType{}, + timingsStream, + input, + false, // use the unaugmented hierarchical tree (i.e., the base tree) for the volume computation + unaugmentedIntrinsicVolumes, + unaugmentedDependentVolumes); + timingsStream << " " << std::setw(38) << std::left << "Compute Volume for Presimplication" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + timer.Start(); + } + // ******** 3. Augment the hierarchical tree if requested ******** if (this->AugmentHierarchicalTree) { - master.foreach ( - [](DistributedContourTreeBlockData* blockData, const vtkmdiy::Master::ProxyWithLink&) { - blockData->HierarchicalAugmenter.Initialize( - blockData->GlobalBlockId, &blockData->HierarchicalTree, &blockData->AugmentedTree); - }); + vtkm::Id localPresimplifyThreshold = this->PresimplifyThreshold; + master.foreach ([globalPointDimensions, localPresimplifyThreshold, unaugmentedDependentVolumes]( + DistributedContourTreeBlockData* blockData, + const vtkmdiy::Master::ProxyWithLink&) { + // if we don't presimplify then use a NULL pointer for the dependent volume used for pre-simplification + vtkm::worklet::contourtree_augmented::IdArrayType* volumeArrayForPresimplifiction = NULL; + // if we presimplify then get a pointer for the dependent volume for the current block + if (localPresimplifyThreshold > 0) + { + volumeArrayForPresimplifiction = + const_cast( + &unaugmentedDependentVolumes[blockData->LocalBlockNo]); + } + // Initialize the hierarchical augmenter + blockData->HierarchicalAugmenter.Initialize( + blockData->GlobalBlockId, + &blockData->HierarchicalTree, + &blockData->AugmentedTree, + blockData->BlockOrigin, // Origin of the data block + blockData->BlockSize, // Extends of the data block + globalPointDimensions, // global point dimensions + volumeArrayForPresimplifiction, // DependentVolume if we computed it or NULL if no presimplification is used + localPresimplifyThreshold // presimplify if threshold is > 0 + ); + }); timingsStream << " " << std::setw(38) << std::left << "Initalize Hierarchical Trees" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; @@ -1259,8 +1307,34 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( if (this->AugmentHierarchicalTree) { - this->ComputeVolumeMetric( - master, assigner, partners, FieldType{}, timingsStream, hierarchicalTreeOutputDataSet); + std::vector> augmentedIntrinsicVolumes; + std::vector> augmentedDependentVolumes; + this->ComputeVolumeMetric(master, + assigner, + partners, + FieldType{}, + timingsStream, + input, + true, // use the augmented tree + augmentedIntrinsicVolumes, + augmentedDependentVolumes); + timer.Start(); + + master.foreach ( + [&](DistributedContourTreeBlockData* blockData, const vtkmdiy::Master::ProxyWithLink&) { + // Add the intrinsic and dependent volumes to the output data set + vtkm::cont::Field intrinsicVolumeField("IntrinsicVolume", + vtkm::cont::Field::Association::WholeDataSet, + augmentedIntrinsicVolumes[blockData->LocalBlockNo]); + hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(intrinsicVolumeField); + vtkm::cont::Field dependentVolumeField("DependentVolume", + vtkm::cont::Field::Association::WholeDataSet, + augmentedDependentVolumes[blockData->LocalBlockNo]); + hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(dependentVolumeField); + // Log the time for adding hypersweep data to the output dataset + timingsStream << " " << std::setw(38) << std::left << "Add Volume Output Data" + << ": " << timer.GetElapsedTime() << " seconds" << std::endl; + }); } VTKM_LOG_S(this->TimingsLogLevel, diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h index 23610971a..525fa4814 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h +++ b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h @@ -118,6 +118,11 @@ public: this->AugmentHierarchicalTree = augmentHierarchicalTree; } + VTKM_CONT void SetPresimplifyThreshold(vtkm::Id presimplifyThreshold) + { + this->PresimplifyThreshold = presimplifyThreshold; + } + VTKM_CONT void SetBlockIndices(vtkm::Id3 blocksPerDim, const vtkm::cont::ArrayHandle& localBlockIndices) { @@ -127,6 +132,8 @@ public: VTKM_CONT bool GetAugmentHierarchicalTree() { return this->AugmentHierarchicalTree; } + VTKM_CONT vtkm::Id GetPresimplifyThreshold() { return this->PresimplifyThreshold; } + VTKM_CONT void SetSaveDotFiles(bool saveDotFiles) { this->SaveDotFiles = saveDotFiles; } VTKM_CONT bool GetSaveDotFiles() { return this->SaveDotFiles; } @@ -166,7 +173,10 @@ private: vtkmdiy::RegularSwapPartners& partners, const FieldType&, // dummy parameter to get the type std::stringstream& timingsStream, - std::vector& hierarchicalTreeOutputDataSet); + const vtkm::cont::PartitionedDataSet& input, + bool useAugmentedTree, + std::vector>& intrinsicVolumes, + std::vector>& dependentVolumes); /// /// Internal helper function that implements the actual functionality of PostExecute @@ -188,6 +198,9 @@ private: /// Augment hierarchical tree bool AugmentHierarchicalTree; + /// Threshold to use for volume pre-simplification + vtkm::Id PresimplifyThreshold; + /// Save dot files for all tree computations bool SaveDotFiles; diff --git a/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h b/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h index 35a8b8609..034be5d6e 100644 --- a/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h +++ b/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h @@ -506,12 +506,19 @@ inline void HierarchicalVolumetricBranchDecomposer::CollapseBranches( vtkm::worklet::contourtree_distributed::FindSuperArcBetweenNodes findSuperArcBetweenNodes{ hierarchicalTreeSuperarcs }; + // Get the number of rounds + auto numRoundsArray = hierarchicalTreeDataSet.GetField("NumRounds") + .GetData() + .AsArrayHandle>(); + vtkm::Id numRounds = vtkm::cont::ArrayGetValue(0, numRoundsArray); using vtkm::worklet::scalar_topology::hierarchical_volumetric_branch_decomposer:: CollapseBranchesWorklet; - this->Invoke(CollapseBranchesWorklet{}, // the worklet + CollapseBranchesWorklet collapseBranchesWorklet(numRounds); + this->Invoke(collapseBranchesWorklet, // the worklet this->BestUpSupernode, // input this->BestDownSupernode, // input + hierarchicalTreeSuperarcs, // input findRegularByGlobal, // input ExecutionObject findSuperArcBetweenNodes, // input ExecutionObject hierarchicalTreeRegular2Supernode, // input diff --git a/vtkm/filter/scalar_topology/worklet/branch_decomposition/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h b/vtkm/filter/scalar_topology/worklet/branch_decomposition/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h index 498f7eb5f..d16e8b108 100644 --- a/vtkm/filter/scalar_topology/worklet/branch_decomposition/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/branch_decomposition/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h @@ -65,6 +65,7 @@ public: using ControlSignature = void( FieldIn bestUpSupernode, FieldIn bestDownSupernode, + FieldIn superarcs, // Execution objects from the hierarchical tree to use the FindRegularByGlobal function ExecObject findRegularByGlobal, // Execution objects from the hierarchical tree to use the FindSuperArcBetweenNodes, function @@ -72,12 +73,15 @@ public: WholeArrayIn hierarchicalTreeRegular2supernode, WholeArrayIn hierarchicalTreeWhichRound, WholeArrayInOut branchRoot); - using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8); using InputDomain = _1; /// Default Constructor VTKM_EXEC_CONT - CollapseBranchesWorklet() {} + CollapseBranchesWorklet(vtkm::Id numRounds) + : NumRounds(numRounds) + { + } /// operator() of the workelt template NumRounds) && + (vtkm::worklet::contourtree_augmented::NoSuchElement(superarcsId))) + { + return; + } + // 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)) { @@ -233,6 +250,8 @@ public: */ } // operator()() +private: + vtkm::Id NumRounds; }; // CollapseBranchesWorklet diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h index e1d02f68b..dbac3b6c2 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h @@ -86,12 +86,20 @@ constexpr vtkm::Id IS_REGULAR = static_cast(2); constexpr vtkm::Id IS_SADDLE = static_cast(3); constexpr vtkm::Id IS_ATTACHMENT = static_cast(4); +// TERMINAL_ELEMENT is primarily used for optimisation of memory access during pointer doubling operations +// We now need to distinguish between a supernode and superarc when sorting by superarc(node) IDs +// This only (at present) comes up when processing attachment points, which have null superarcs, so it +// is reasonable to reuse TERMINAL_ELEMENT for this purpose. However, we give it a separate macro name with +// the same value to aid comprehension +constexpr vtkm::Id TRANSFER_TO_SUPERARC = TERMINAL_ELEMENT; + // clang-format on using IdArrayType = vtkm::cont::ArrayHandle; using EdgePair = vtkm::Pair; // here EdgePair.first=low and EdgePair.second=high using EdgePairArray = vtkm::cont::ArrayHandle; // Array of edge pairs + // inline functions for retrieving flags or index VTKM_EXEC_CONT inline bool NoSuchElement(vtkm::Id flaggedIndex) @@ -143,6 +151,12 @@ inline bool NoFlagsSet(vtkm::Id flaggedIndex) return (flaggedIndex & ~INDEX_MASK) == 0; } // NoFlagsSet +// Helper function: to check that the TRANSFER_TO_SUPERARC flag is set +VTKM_EXEC_CONT +inline bool TransferToSuperarc(vtkm::Id flaggedIndex) +{ // transferToSuperarc() + return ((flaggedIndex & TRANSFER_TO_SUPERARC) != 0); +} // transferToSuperarc() // Debug helper function: Assert that an index array has no element with any flags set template @@ -225,6 +239,23 @@ inline std::string FlagString(vtkm::Id flaggedIndex) return fString; } // FlagString() + +// == comparison operator for edges +inline bool edgeEqual(const EdgePair& LHS, const EdgePair& RHS) +{ // operator == + + if (LHS.first != RHS.first) + { + return false; + } + if (LHS.second != RHS.second) + { + return false; + } + return true; +} // operator == + + class EdgeDataHeight { public: diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalAugmenter.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalAugmenter.h index a9ad55f7f..13db70bfb 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalAugmenter.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalAugmenter.h @@ -101,6 +101,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -127,6 +130,11 @@ template class HierarchicalAugmenter { // class HierarchicalAugmenter public: + /// base mesh variable needs to determine whether a vertex is inside or outside of the block + vtkm::Id3 MeshBlockOrigin; + vtkm::Id3 MeshBlockSize; + vtkm::Id3 MeshGlobalSize; + /// the tree that it hypersweeps over vtkm::worklet::contourtree_distributed::HierarchicalContourTree* BaseTree; /// the tree that it is building @@ -173,7 +181,7 @@ public: /// these are essentially temporary local variables, but are placed here to make the DebugPrint() /// more comprehensive. They will be allocated where used - /// this one makes a list of attachment Ids and is used in sevral different places, so resize it when done + /// this one makes a list of attachment Ids and is used in several different places, so resize it when done vtkm::worklet::contourtree_augmented::IdArrayType AttachmentIds; /// tracks segments of attachment points by round vtkm::worklet::contourtree_augmented::IdArrayType FirstAttachmentPointInRound; @@ -198,7 +206,12 @@ public: void Initialize( vtkm::Id blockId, vtkm::worklet::contourtree_distributed::HierarchicalContourTree* inBaseTree, - vtkm::worklet::contourtree_distributed::HierarchicalContourTree* inAugmentedTree); + vtkm::worklet::contourtree_distributed::HierarchicalContourTree* inAugmentedTree, + vtkm::Id3 meshBlockOrigin, + vtkm::Id3 meshBockSize, + vtkm::Id3 meshGlobalSize, + vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray = NULL, + vtkm::Id presimplifyThreshold = 0); /// routine to prepare the set of attachment points to transfer void PrepareOutAttachmentPoints(vtkm::Id round); @@ -220,7 +233,7 @@ public: /// transfer level of superstructure with insertions void CopySuperstructure(); /// reset the super Ids in the hyperstructure to the new values - void UpdateHyperstructure(); + void UpdateHyperstructure(vtkm::Id roundNumber); /// copy the remaining base level regular nodes void CopyBaseRegularStructure(); @@ -249,22 +262,36 @@ template void HierarchicalAugmenter::Initialize( vtkm::Id blockId, vtkm::worklet::contourtree_distributed::HierarchicalContourTree* baseTree, - vtkm::worklet::contourtree_distributed::HierarchicalContourTree* augmentedTree) + vtkm::worklet::contourtree_distributed::HierarchicalContourTree* augmentedTree, + vtkm::Id3 meshBlockOrigin, + vtkm::Id3 meshBockSize, + vtkm::Id3 meshGlobalSize, + vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray, + vtkm::Id presimplifyThreshold) { // Initialize() // copy the parameters for use this->BlockId = blockId; this->BaseTree = baseTree; this->AugmentedTree = augmentedTree; + this->MeshBlockOrigin = meshBlockOrigin; + this->MeshBlockSize = meshBockSize; + this->MeshGlobalSize = meshGlobalSize; // now construct a list of all attachment points on the block + // except those under the presimplify threshold. The presimplification is + // handled in the IsAttachementPointPredicate + // to do this, we construct an index array with all supernode ID's that satisfy: // 1. superparent == NO_SUCH_ELEMENT (i.e. root of interior tree) // 2. round < nRounds (except the top level, where 1. indicates the tree root) - // initalize AttachementIds + // initialize AttachementIds { vtkm::worklet::contourtree_distributed::hierarchical_augmenter::IsAttachementPointPredicate - isAttachementPointPredicate( - this->BaseTree->Superarcs, this->BaseTree->WhichRound, this->BaseTree->NumRounds); + isAttachementPointPredicate(this->BaseTree->Superarcs, + this->BaseTree->WhichRound, + this->BaseTree->NumRounds, + volumeArray, + presimplifyThreshold); auto tempSupernodeIndex = vtkm::cont::ArrayHandleIndex(this->BaseTree->Supernodes.GetNumberOfValues()); vtkm::cont::Algorithm::CopyIf( @@ -476,7 +503,11 @@ void HierarchicalAugmenter::ReleaseSwapArrays() } // ReleaseSwapArrays() -// routine to reconstruct a hierarchical tree using the augmenting supernodes +// routine to reconstruct a hierarchical tree using the augmenting supernodes. +// Allowing pre-simplification needs the superstructure and hyperstructure to be done one layer +// at a time. This means lifting the code from CopySuperstructure up to here CopyHyperstructure() +// can actually be left the way it is - copying the entire hyperstructure, as the augmentation +// process doesn't change it. It might be sensible to change it anyway, just to make the code better organised template void HierarchicalAugmenter::BuildAugmentedTree() { // BuildAugmentedTree() @@ -484,11 +515,22 @@ void HierarchicalAugmenter::BuildAugmentedTree() this->PrepareAugmentedTree(); // 2. Copy the hyperstructure, using the old super IDs for now this->CopyHyperstructure(); - // 3. Copy the superstructure, inserting additional points as we do - this->CopySuperstructure(); - // 4. Update the hyperstructure to use the new super IDs - this->UpdateHyperstructure(); - // 5. Copy the remaining regular structure at the bottom level, setting up the regular sort order in the process + // 3. Copy superstructure one round at a time, updating the hyperstructure as well + // (needed to permit search for superarcs) + // Loop from the top down: + for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--) + { // per round + // start by retrieving list of old supernodes from the tree (except for attachment points) + this->RetrieveOldSupernodes(roundNumber); + // since we know the number of attachment points, we can now allocate space for the level + // and set up arrays for sorting the supernodes + this->ResizeArrays(roundNumber); + // now we create the superarcs for the round in the new tree + this->CreateSuperarcs(roundNumber); + // finally, we update the hyperstructure for the round in the new tree + this->UpdateHyperstructure(roundNumber); + } // per round + // 4. Copy the remaining regular structure at the bottom level, setting up the regular sort order in the process this->CopyBaseRegularStructure(); } // BuildAugmentedTree() @@ -623,53 +665,91 @@ void HierarchicalAugmenter::CopyHyperstructure() this->AugmentedTree->FirstHypernodePerIteration[roundNumber]); } // per round + // WARNING 28/05/2023: Since this resize is for the full hyperstructure, it should be safe to put here. + // Unless of course, anything relies on the sizes: but they were 0, so it is unlikely. + // A search for hyperarcs.size() & hypernodes.size() in this unit confirmed that nothing uses them. + // Nevertheless, set them all to NO_SUCH_ELEMENT out of paranoia + // 5. Reset hypernodes, hyperarcs and superchildren using supernode IDs + // The hyperstructure is unchanged, but uses old supernode IDs + vtkm::worklet::contourtree_augmented::ResizeVector( + this->AugmentedTree->Hypernodes, // resize array + this->BaseTree->Hypernodes.GetNumberOfValues(), // new size + vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT // set all elements to this value + ); + vtkm::worklet::contourtree_augmented::ResizeVector( + this->AugmentedTree->Hyperarcs, // resize array + this->BaseTree->Hyperarcs.GetNumberOfValues(), // new size + vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT // set all elements to this value + ); + vtkm::worklet::contourtree_augmented::ResizeVector( + this->AugmentedTree->Superchildren, // resize array + this->BaseTree->Superchildren.GetNumberOfValues(), // new size + 0 // set all elements to this value + ); + #ifdef DEBUG_PRINT VTKM_LOG_S(vtkm::cont::LogLevel::Info, DebugPrint("Hyperstructure Copied", __FILE__, __LINE__)); #endif } // CopyHyperstructure() +// WARNING 28/05/2023: deleted and moved up to AugmentTree() routine // transfer level of superstructure with insertions -template -void HierarchicalAugmenter::CopySuperstructure() -{ // CopySuperstructure() - - // Loop from the top down: - for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--) - { // per round - // start by retrieving list of old supernodes from the tree (except for attachment points) - this->RetrieveOldSupernodes(roundNumber); - // since we know the number of attachment points, we can now allocate space for the level - // and set up arrays for sorting the supernodes - this->ResizeArrays(roundNumber); - // now we create the superarcs for the round in the new tree - this->CreateSuperarcs(roundNumber); - } // per round -#ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, - this->DebugPrint("Superstructure Copied", __FILE__, __LINE__)); -#endif -} // CopySuperstructure() +//template +//void HierarchicalAugmenter::CopySuperstructure() +//{ // CopySuperstructure() +// +// // Loop from the top down: +// for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--) +// { // per round +// // start by retrieving list of old supernodes from the tree (except for attachment points) +// this->RetrieveOldSupernodes(roundNumber); +// // since we know the number of attachment points, we can now allocate space for the level +// // and set up arrays for sorting the supernodes +// this->ResizeArrays(roundNumber); +// // now we create the superarcs for the round in the new tree +// this->CreateSuperarcs(roundNumber); +// } // per round +//#ifdef DEBUG_PRINT +// VTKM_LOG_S(vtkm::cont::LogLevel::Info, +// this->DebugPrint("Superstructure Copied", __FILE__, __LINE__)); +//#endif +//} // CopySuperstructure() // reset the super IDs in the hyperstructure to the new values template -void HierarchicalAugmenter::UpdateHyperstructure() +void HierarchicalAugmenter::UpdateHyperstructure(vtkm::Id roundNumber) { // UpdateHyperstructure() - - // 5. Reset hypernodes, hyperarcs and superchildren using supernode IDs - // The hyperstructure is unchanged, but uses old supernode IDs + // now that the superstructure is known, we can find the new supernode IDs for all + // of the old hypernodes at this level and update.Wwe want to update the entire round + // at once, so we would like to use the firstHypernodePerIteration array. + vtkm::Id startIndex = + vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstHypernodePerIteration[roundNumber]); + vtkm::Id stopIndex = vtkm::cont::ArrayGetValue( + vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations), + this->AugmentedTree->FirstHypernodePerIteration[roundNumber]); + vtkm::Id selectSize = stopIndex - startIndex; { vtkm::worklet::contourtree_distributed::hierarchical_augmenter:: UpdateHyperstructureSetHyperarcsAndNodesWorklet updateHyperstructureSetHyperarcsAndNodesWorklet; - this->Invoke( - updateHyperstructureSetHyperarcsAndNodesWorklet, - this->BaseTree->Hypernodes, // input - this->BaseTree->Hyperarcs, // input - this->NewSupernodeIds, // input - this->AugmentedTree->Hypernodes, // output (the array is automatically resized here) - this->AugmentedTree->Hyperarcs // output (the array is automatically resized here) + // Create ArrayHandleViews of the subrange of the input and output arrays we need to process + auto baseTreeHypernodesView = + vtkm::cont::make_ArrayHandleView(this->BaseTree->Hypernodes, startIndex, selectSize); + auto baseTreeHyperarcsView = + vtkm::cont::make_ArrayHandleView(this->BaseTree->Hyperarcs, startIndex, selectSize); + auto augmentedTreeHypernodesView = + vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hypernodes, startIndex, selectSize); + auto augmentedTreeHyperarcsView = + vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hyperarcs, startIndex, selectSize); + // Invoke the worklet with the subset view of our arrays + this->Invoke(updateHyperstructureSetHyperarcsAndNodesWorklet, + baseTreeHypernodesView, // input + baseTreeHyperarcsView, // input + this->NewSupernodeIds, // input + augmentedTreeHypernodesView, // output (the array is automatically resized here) + augmentedTreeHyperarcsView // output (the array is automatically resized here) ); } @@ -679,10 +759,20 @@ void HierarchicalAugmenter::UpdateHyperstructure() vtkm::worklet::contourtree_distributed::hierarchical_augmenter:: UpdateHyperstructureSetSuperchildrenWorklet updateHyperstructureSetSuperchildrenWorklet( this->AugmentedTree->Supernodes.GetNumberOfValues()); - this->Invoke( - updateHyperstructureSetSuperchildrenWorklet, - this->AugmentedTree->Hypernodes, // input - this->AugmentedTree->Superchildren // output (the array is automatically resized here) + // As above, we need to create views of the relevant subranges of our arrays + auto augmentedTreeSuperarcsView = + vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs, startIndex, selectSize); + vtkm::Id extraSelectSize = + ((startIndex + selectSize) < this->AugmentedTree->Hyperparents.GetNumberOfValues()) + ? selectSize + 1 + : selectSize; + auto augmentedTreeHyperparentsView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Hyperparents, startIndex, extraSelectSize); + this->Invoke(updateHyperstructureSetSuperchildrenWorklet, + this->AugmentedTree->Hypernodes, // input + augmentedTreeSuperarcsView, // input + augmentedTreeHyperparentsView, // inpu + this->AugmentedTree->Superchildren // output ); } } // UpdateHyperstructure() @@ -728,12 +818,13 @@ void HierarchicalAugmenter::CopyBaseRegularStructure() vtkm::worklet::contourtree_augmented::IdArrayType tempRegularNodesNeeded; // create the worklet vtkm::worklet::contourtree_distributed::hierarchical_augmenter:: - FindSuperparentForNecessaryNodesWorklet findSuperparentForNecessaryNodesWorklet; + FindSuperparentForNecessaryNodesWorklet findSuperparentForNecessaryNodesWorklet( + this->MeshBlockOrigin, this->MeshBlockSize, this->MeshGlobalSize); // Get a FindRegularByGlobal and FindSuperArcForUnknownNode execution object for our worklet auto findRegularByGlobal = this->AugmentedTree->GetFindRegularByGlobal(); auto findSuperArcForUnknownNode = this->AugmentedTree->GetFindSuperArcForUnknownNode(); - // excute the worklet + // execute the worklet this->Invoke(findSuperparentForNecessaryNodesWorklet, // the worklet to call // inputs this->BaseTree->RegularNodeGlobalIds, // input domain @@ -831,6 +922,12 @@ void HierarchicalAugmenter::CopyBaseRegularStructure() ); } + + // Reset the number of regular nodes in round 0 + vtkm::Id regularNodesInRound0 = + numTotalRegular - this->AugmentedTree->NumRegularNodesInRound.ReadPortal().Get(1); + this->AugmentedTree->NumRegularNodesInRound.WritePortal().Set(0, regularNodesInRound0); + // Finally, we resort the regular node sort order { vtkm::worklet::contourtree_distributed::PermuteComparator // hierarchical_contour_tree:: @@ -850,23 +947,22 @@ void HierarchicalAugmenter::RetrieveOldSupernodes(vtkm::Id roundNumbe // TODO PERFORMANCE STATISTICS: // the # of supernodes at each level minus the # of kept supernodes gives us the # of attachment points we lose at this level // in addition to this, the firstAttachmentPointInRound array gives us the # of attachment points we gain at this level + // + // COMMENT: In adding presimplification, this will have to change. + // Previously, it made the hard assumption that all attachment points were transferred & used that to suppress + // them. Now it can do that no longer, which gives us a choice. We can pass in the threshold & volume array + // and test here, but that's not ideal as it does the same test in multiple places. Alternatively, we can + // look up whether the supernode is already present in the structure, which has an associated search cost. + // BUT, we have an array called newSupernodeIDs for the purpose already, so that's how we'll do it. + vtkm::Id supernodeIndexBase = vtkm::cont::ArrayGetValue(0, this->BaseTree->FirstSupernodePerIteration[roundNumber]); vtkm::cont::ArrayHandleCounting supernodeIdVals( supernodeIndexBase, // start 1, // step this->BaseTree->NumSupernodesInRound.ReadPortal().Get(roundNumber)); - // the test for whether to keep it is: - // a1. at the top level, keep everything - if (!(roundNumber < this->BaseTree->NumRounds)) - { - vtkm::cont::Algorithm::Copy(supernodeIdVals, this->KeptSupernodes); - } - // a2. at lower levels, keep them if the superarc is NO_SUCH_ELEMENT - else { // Reset this-KeptSupernodes to the right size and initalize with NO_SUCH_ELEMENT. - // TODO: Check if a simple free and allocate without initalizing the array is sufficient vtkm::cont::Algorithm::Copy( // Create const array to copy vtkm::cont::ArrayHandleConstant( @@ -882,7 +978,7 @@ void HierarchicalAugmenter::RetrieveOldSupernodes(vtkm::Id roundNumbe supernodeIdVals, // Stencil with baseTree->superarcs[supernodeID] vtkm::cont::make_ArrayHandleView( - this->BaseTree->Superarcs, supernodeIndexBase, this->KeptSupernodes.GetNumberOfValues()), + this->NewSupernodeIds, supernodeIndexBase, this->KeptSupernodes.GetNumberOfValues()), // And the CopyIf compresses the array to eliminate unnecssary elements // save to this->KeptSupernodes this->KeptSupernodes, @@ -908,6 +1004,10 @@ template void HierarchicalAugmenter::ResizeArrays(vtkm::Id roundNumber) { // ResizeArrays() // at this point, we know how many supernodes are kept from the same level of the old tree + // TODO/WARNING 23/07/2023: + // Actually, this is no longer true. If the same vertex is an attachment point for two adjacent blocks (i.e. it is on the boundary), it is entirely possible + // for one block to add it at a higher level than the other. To preclude this, we will need to edit RetrieveOldSupernodes() + // we can also find out how many supernodes are being inserted, which gives us the correct amount to expand by, saving a double resize() call // note that some of these arrays could probably be resized later, but it's cleaner this way // also note that if it becomes a problem, we could resize all of the arrays to baseTree->supernodes.size() + # of attachmentPoints as an over-estimate @@ -1013,6 +1113,9 @@ void HierarchicalAugmenter::ResizeArrays(vtkm::Id roundNumber) #endif // b. Transfer attachment points for level into new supernode array + // NOTE: this means the set of attachment points that we have determined by swapping + // need to be inserted onto a superarc at this level. All of them should be from + // lower levels originally, but are being moved up to this level for insertion // to copy them in, we use the existing array of attachment point IDs by round { vtkm::Id firstAttachmentPointInRoundCurrent = @@ -1061,7 +1164,10 @@ void HierarchicalAugmenter::ResizeArrays(vtkm::Id roundNumber) __LINE__)); #endif - // Now we copy in the kept supernodes + // Now we copy in the kept supernodes: this used to mean only the non-attachment points + // now it includes the attachment points at this level that the simplification removed + // so they need to be put back where they were + // However, that means that all of them do exist in the base tree, so we can copy from there { auto oldRegularIdArr = vtkm::cont::make_ArrayHandlePermutation(this->KeptSupernodes, // index @@ -1135,15 +1241,24 @@ void HierarchicalAugmenter::ResizeArrays(vtkm::Id roundNumber) ResizeArraysBuildNewSupernodeIdsWorklet resizeArraysBuildNewSupernodeIdsWorklet( numSupernodesAlready); auto supernodeIndex = vtkm::cont::ArrayHandleIndex(this->SupernodeSorter.GetNumberOfValues()); - auto supernodeIdSetPermuted = - vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet); + // auto supernodeIdSetPermuted = + // vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet); + auto globalRegularIdSetPermuted = + vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet); + auto findRegularByGlobal = this->BaseTree->GetFindRegularByGlobal(); this->Invoke( resizeArraysBuildNewSupernodeIdsWorklet, supernodeIndex, // input domain. We only need the index because supernodeIdSetPermuted already does the permute - supernodeIdSetPermuted, // input input supernodeIDSet permuted by supernodeSorter to allow for FieldIn + // supernodeIdSetPermuted, // input input supernodeIDSet permuted by supernodeSorter to allow for FieldIn + globalRegularIdSetPermuted, + findRegularByGlobal, + this->BaseTree->Regular2Supernode, this ->NewSupernodeIds // output/input (both are necessary since not all values will be overwritten) ); + + // Add const ExecObjectType1& findRegularByGlobal, + // Add baseTree->regular2supernode } #ifdef DEBUG_PRINT @@ -1161,167 +1276,250 @@ template void HierarchicalAugmenter::CreateSuperarcs(vtkm::Id roundNumber) { // CreateSuperarcs() // retrieve the ID number of the first supernode at this level +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + DebugPrint(std::string("Round ") + std::to_string(roundNumber) + + std::string(" Starting CreateSuperarcs()"), + __FILE__, + __LINE__)); +#endif + + vtkm::Id currNumIterations = + vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations); vtkm::Id numSupernodesAlready = vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstSupernodePerIteration[roundNumber]); - // e. Connect superarcs for the level & set hyperparents & superchildren count, whichRound, whichIteration, super2hypernode - { // START scope for e. to delete temporary variables - // Note: The original PPP algorithm performed all operations listed in this block - // in a single parralel for loop. Many of those operations were smart array - // copies. So to simplfy the worklet and to make more effective use of - // VTKm algorithm, a large number of copy operations have been extracted from - // the loop and are performed here via combinations of fancy array handles and - // vtkm::cont::Algorithm::Copy operations. + // e. Connect superarcs for the level & set hyperparents & superchildren count, whichRound, whichIteration, super2hypernode + // 24/05/2023: Expansion of comment to help debug. + // At this point, we know that all higher rounds are correctly constructed, and that any attachment points that survived simplification + // have already been inserted in a higher round. - // Define the new supernode and regular Id. Both are actually the same here since we are - // augmenting the tree here, but for clarity we define them as separate variables. - // At all levels above 0, we used to keep regular vertices in case - // they are attachment points. After augmentation, we don't need to. - // Instead, at all levels above 0, the regular nodes in each round - // are identical to the supernodes. In order to avoid confusion, - // we will copy the Id into a separate variable - vtkm::cont::ArrayHandleCounting newSupernodeId( - numSupernodesAlready, // start - static_cast(1), // step - this->SupernodeSorter.GetNumberOfValues() // number of values - ); - auto newRegularId = newSupernodeId; - // define the superparentOldSuperId - auto permutedSuperparentSet = - vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SuperparentSet); - auto superparentOldSuperId = vtkm::cont::make_ArrayHandleTransform( - permutedSuperparentSet, vtkm::worklet::contourtree_augmented::MaskedIndexFunctor()); + // The sort should have resulted in the supernodes being segmented along old superarcs. Most supernodes should be in a segment of length + // 1, and should be their own superparent in the sort array. But we can't readily test that, because other supernodes may also have them + // as the superparent. - // set the supernode's regular Id. Set: augmentedTree->supernodes[newSupernodeID] = newRegularID; - auto permutedAugmentedTreeSupernodes = - vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->Supernodes); - vtkm::cont::Algorithm::Copy(newRegularId, permutedAugmentedTreeSupernodes); + // This loop will principally determine the superarc for each supernode. For this, the rules break down to: + // 1. If the supernode is the global root, connect it nowhere + // 2. If the supernode is the last in the set of all supernodes in this round, treat it as the end of a segment + // 3. If the supernode is the last in a segment by superarc, connect it to the target of its superparent in the old tree, using the new supernode ID + // 4. Otherwise, connect to the new supernode ID of the next supernode in the segment - // Run the worklet for more complex operations - { // START block for CreateSuperarcsWorklet - // create the worklet - vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsWorklet - createSuperarcsWorklet( - numSupernodesAlready, - this->BaseTree->NumRounds, - vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations), - roundNumber, - this->AugmentedTree->Supernodes.GetNumberOfValues()); - // create fancy arrays needed to allow use of FieldIn for worklet parameters - auto permutedGlobalRegularIdSet = - vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet); - auto augmentedTreeSuperarcsView = - vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs, - numSupernodesAlready, - this->SupernodeSorter.GetNumberOfValues()); - auto augmentedTreeSuper2HypernodeView = - vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Super2Hypernode, - numSupernodesAlready, - this->SupernodeSorter.GetNumberOfValues()); - // invoke the worklet - this->Invoke(createSuperarcsWorklet, // the worklet - this->SupernodeSorter, // input domain - this->SuperparentSet, // input - this->BaseTree->Superarcs, // input - this->NewSupernodeIds, // input - this->BaseTree->Supernodes, // input - this->BaseTree->RegularNodeGlobalIds, // input - permutedGlobalRegularIdSet, // input - this->BaseTree->Super2Hypernode, // input - this->BaseTree->WhichIteration, // input - augmentedTreeSuperarcsView, // output - this->AugmentedTree->FirstSupernodePerIteration[roundNumber], // input/output - augmentedTreeSuper2HypernodeView // output - ); - } // END block for CreateSuperarcsWorklet + // In each case, we will need to preserve the ascending / descending flag - // setting the hyperparent is straightforward since the hyperstructure is preserved - // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that - // Set: augmentedTree->hyperparents[newSupernodeID] = baseTree->hyperparents[superparentOldSuperID]; - auto permutedAugmentedTreeHyperparents = - vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->Hyperparents); - auto permutedBaseTreeHyperparents = - vtkm::cont::make_ArrayHandlePermutation(superparentOldSuperId, this->BaseTree->Hyperparents); - vtkm::cont::Algorithm::Copy(permutedBaseTreeHyperparents, permutedAugmentedTreeHyperparents); + // We will also have to set the first supernode per iteration - if possible, in a separate loop - // which round and iteration carry over - // Set: augmentedTree->whichRound[newSupernodeID] = baseTree->whichRound[superparentOldSuperID]; - auto permutedAugmentedTreeWhichRound = - vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->WhichRound); - auto permutedBaseTreeWhichRound = - vtkm::cont::make_ArrayHandlePermutation(superparentOldSuperId, this->BaseTree->WhichRound); - vtkm::cont::Algorithm::Copy(permutedBaseTreeWhichRound, permutedAugmentedTreeWhichRound); + // We need this to determine which supernodes are inserted and which are attached (see below) + vtkm::Id numInsertedSupernodes = + (vtkm::cont::ArrayGetValue(roundNumber + 1, this->FirstAttachmentPointInRound) - + vtkm::cont::ArrayGetValue(roundNumber, this->FirstAttachmentPointInRound)); - // Set: augmentedTree->whichIteration[newSupernodeID] = baseTree->whichIteration[superparentOldSuperID]; - auto permutedAugmentedTreeWhichIteration = - vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->WhichIteration); - auto permutedBaseTreeWhichIterationPortal = vtkm::cont::make_ArrayHandlePermutation( - superparentOldSuperId, this->BaseTree->WhichIteration); - vtkm::cont::Algorithm::Copy(permutedBaseTreeWhichIterationPortal, - permutedAugmentedTreeWhichIteration); + { // START Call CreateSuperarcsWorklet (scope to delete temporary variables) + // create the worklet + vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsWorklet< + FieldType> + createSuperarcsWorklet( + numSupernodesAlready, this->BaseTree->NumRounds, numInsertedSupernodes, roundNumber); - // now we deal with the regular-sized arrays. In the following supernodeSetIndex is simply supernodeSorterPortal.Get(supernode); - // copy the global regular Id and data value - // Set: augmentedTree->regularNodeGlobalIDs[newRegularID] = globalRegularIDSet[supernodeSetIndex]; - auto permutedAugmentedTreeRegularNodeGlobalIds = vtkm::cont::make_ArrayHandlePermutation( - newRegularId, this->AugmentedTree->RegularNodeGlobalIds); + // create fancy arrays needed to allow use of FieldIn for worklet parameters + auto permutedSupernodeIdSet = + vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet); + auto permutedBaseTreeSuper2Hypernode = vtkm::cont::make_ArrayHandlePermutation( + permutedSupernodeIdSet, + this->BaseTree + ->Super2Hypernode); // Get this->BaseTree->Super2hypernode[oldSupernodeId]; as FieldIn for the worklet auto permutedGlobalRegularIdSet = vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet); - vtkm::cont::Algorithm::Copy(permutedGlobalRegularIdSet, - permutedAugmentedTreeRegularNodeGlobalIds); - // SetL augmentedTree->dataValues[newRegularID] = dataValueSet[supernodeSetIndex]; - auto permutedAugmentedTreeDataValues = - vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->DataValues); auto permutedDataValueSet = vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->DataValueSet); - vtkm::cont::Algorithm::Copy(permutedDataValueSet, permutedAugmentedTreeDataValues); + auto oldRegularId = + vtkm::cont::make_ArrayHandlePermutation(permutedSupernodeIdSet, this->BaseTree->Supernodes); + auto oldSuperFrom = + vtkm::cont::make_ArrayHandlePermutation(oldRegularId, this->BaseTree->Superparents); - // the sort order will be dealt with later - // since all of these nodes are supernodes, they will be their own superparent, which means that: - // a. the regular2node can be set immediately - // Set: augmentedTree->regular2supernode[newRegularID] = newSupernodeID; - auto permutedAugmentedTreeRegular2Supernode = - vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->Regular2Supernode); - vtkm::cont::Algorithm::Copy(newSupernodeId, permutedAugmentedTreeRegular2Supernode); + // Create a view of the range of this->AugmentedTree->Superarcs that will be updated by the worklet + // so that we can use FieldOut as the type in the worklet + auto augmentedTreeSuperarcsView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Superarcs, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeHyperparentsView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Hyperparents, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeSuper2HypernodeView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Super2Hypernode, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeWhichRoundView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->WhichRound, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeWhichIterationView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->WhichIteration, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeRegularNodeGlobalIdsView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->RegularNodeGlobalIds, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeDataValuesView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->DataValues, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeRegular2SupernodeView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Regular2Supernode, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); + auto augmentedTreeSuperparentsView = vtkm::cont::make_ArrayHandleView( + this->AugmentedTree->Superparents, + numSupernodesAlready, // start view here + this->SupernodeSorter.GetNumberOfValues() // select this many values + ); - // b. as can the superparent - // Set: augmentedTree->superparents[newRegularID] = newSupernodeID; - auto permutedAugmentedTreeSuperparents = - vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->Superparents); - vtkm::cont::Algorithm::Copy(newSupernodeId, permutedAugmentedTreeSuperparents); - } // END scope for e. to delete temporary variables + // Required execution objects to call other functions + auto findSuperArcForUnknownNode = this->AugmentedTree->GetFindSuperArcForUnknownNode(); - // We have one last bit of cleanup to do. If there were attachment points, then the round in which they transfer has been removed - // While it is possible to turn this into a null round, it is better to reduce the iteration count by one and resize the arrays - // To do this, we access the *LAST* element written and check to see whether it is in the final iteration (according to the base tree) + // Execution object used to encapsulate data form the BaseTree to avoid the limit of 20 input parameters per worklet + auto createSuperarcsDataExecObj = + vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsDataExec( + this->BaseTree->Hyperparents, + this->BaseTree->WhichRound, + this->BaseTree->WhichIteration, + this->BaseTree->Superarcs, + this->BaseTree->Hypernodes, + this->SuperparentSet, + this->NewSupernodeIds); + + // invoke the worklet + this->Invoke(createSuperarcsWorklet, // the worklet + // inputs + this->SupernodeSorter, // input domain (WholeArrayIn) + permutedSupernodeIdSet, // input (FieldIn) + permutedBaseTreeSuper2Hypernode, // input (FieldIn) + permutedGlobalRegularIdSet, // input (FieldIn) + permutedDataValueSet, // input (FieldIn) + oldSuperFrom, // input (FieldIn) + findSuperArcForUnknownNode, // input (Execution object) + createSuperarcsDataExecObj, // input (Execution object with BaseTreeData + // Outputs + this->AugmentedTree->Supernodes, // input/output + augmentedTreeSuperarcsView, // output + augmentedTreeHyperparentsView, // output + augmentedTreeSuper2HypernodeView, // output + augmentedTreeWhichRoundView, // output + augmentedTreeWhichIterationView, // output + augmentedTreeRegularNodeGlobalIdsView, // output + augmentedTreeDataValuesView, // output + augmentedTreeRegular2SupernodeView, // output + augmentedTreeSuperparentsView // output + ); + + } // END Call CreateSuperarcsWorklet + +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + DebugPrint(std::string("Round ") + std::to_string(roundNumber) + + std::string(" Details Filled in For Supernodes "), + __FILE__, + __LINE__)); +#endif + + // Now, in order to set the first supernode per iteration, we do an additional loop + // We are guaranteed that all supernodes at this level are implicitly sorted by iteration, so we test for ends of segments + // NOTE that we do this after the previous loop, since we depend on a value that it has set + { + vtkm::cont::ArrayHandleCounting tempIndex(1, 1, currNumIterations - 1); + vtkm::worklet::contourtree_distributed::hierarchical_augmenter:: + CreateSuperarcsSetFirstSupernodePerIterationWorklet + createSuperarcsSetFirstSupernodePerIterationWorklet(numSupernodesAlready); + vtkm::cont::ArrayHandleIndex tempSupernodeIndex(this->SupernodeSorter.GetNumberOfValues()); + this->Invoke(createSuperarcsSetFirstSupernodePerIterationWorklet, + tempSupernodeIndex, + this->AugmentedTree->WhichIteration, // input + this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // input/output + ); + } + + // since there's an extra entry in the firstSupernode array as a sentinel, set it + vtkm::worklet::contourtree_augmented::IdArraySetValue( + currNumIterations, // index + this->AugmentedTree->Supernodes.GetNumberOfValues(), // new value + this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // array + ); + + // The following loop should be safe in parallel since there should never be two zeros in sequence, i.e., the next + // entry after a zero will always be valid, regardless of execution order + // This was added because in rare cases there are no supernodes transferred in an iteration, for example because there + // are no available upper leaves to prune. If this is case, we are guaranteed that there will be available lower leaves + // so the next iteration will have a non-zero number. We had a major bug from this, and it's cropped back up in the + // Hierarchical Augmentation, so I'm expanding the comment just in case. + { + vtkm::cont::ArrayHandleCounting tempIndex(1, 1, currNumIterations - 1); + vtkm::worklet::contourtree_distributed::hierarchical_augmenter:: + CreateSuperarcsUpdateFirstSupernodePerIterationWorklet + createSuperarcsUpdateFirstSupernodePerIterationWorklet; + this->Invoke(createSuperarcsUpdateFirstSupernodePerIterationWorklet, + tempIndex, // input index + this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // input/output + ); + } + + // We have one last bit of cleanup to do. If there were attachment points, + // then the round in which they transfer has been removed + // While it is possible to turn this into a null round, it is better to + // reduce the iteration count by one and resize the arrays + // To do this, we access the *LAST* element written and check to see whether + // it is in the final iteration (according to the base tree) // But there might be *NO* supernodes in the round, so we check first - vtkm::Id iterationArraySize = - vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations); - if (iterationArraySize > 0) + if (currNumIterations > 0) { // at least one iteration vtkm::Id lastSupernodeThisLevel = this->AugmentedTree->Supernodes.GetNumberOfValues() - 1; vtkm::Id lastIterationThisLevel = vtkm::worklet::contourtree_augmented::MaskedIndex( vtkm::cont::ArrayGetValue(lastSupernodeThisLevel, this->AugmentedTree->WhichIteration)); - // if there were no attachment points, it will be in the last iteration: if there were attachment points, it will be in the previous one - if (lastIterationThisLevel < iterationArraySize - 1) + // if there were no attachment points, it will be in the last iteration: if there were + // attachment points, it will be in the previous one + if (lastIterationThisLevel < currNumIterations - 1) { // attachment point round was removed // decrement the iteration count (still with an extra element as sentinel) + vtkm::Id iterationArraySize = currNumIterations; + // decrease iterations by 1. I.e,: augmentedTree->nIterations[roundNo]--; vtkm::worklet::contourtree_augmented::IdArraySetValue( - roundNumber, iterationArraySize - 1, this->AugmentedTree->NumIterations); + roundNumber, // index + currNumIterations - 1, // new value + this->AugmentedTree->NumIterations // array + ); // shrink the supernode array this->AugmentedTree->FirstSupernodePerIteration[roundNumber].Allocate( iterationArraySize, vtkm::CopyFlag::On); // shrink array but keep values vtkm::worklet::contourtree_augmented::IdArraySetValue( - iterationArraySize - 1, - this->AugmentedTree->Supernodes.GetNumberOfValues(), - this->AugmentedTree->FirstSupernodePerIteration[roundNumber]); + iterationArraySize - 1, // index + this->AugmentedTree->Supernodes.GetNumberOfValues(), // new value + this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // array + ); // for the hypernode array, the last iteration is guaranteed not to have hyperarcs by construction // so the last iteration will already have the correct sentinel value, and we just need to shrink the array - this->AugmentedTree->FirstHypernodePerIteration[roundNumber].Allocate( + this->AugmentedTree->FirstSupernodePerIteration[roundNumber].Allocate( iterationArraySize, vtkm::CopyFlag::On); // shrink array but keep values } // attachment point round was removed } // at least one iteration + +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + DebugPrint(std::string("Round ") + std::to_string(roundNumber) + + std::string(" Superarcs Created "), + __FILE__, + __LINE__)); +#endif + // in the interests of debug, we resize the sorting array to zero here, // even though we will re-resize them in the next function this->SupernodeSorter.ReleaseResources(); diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalContourTree.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalContourTree.h index acd5eca2b..a996d4292 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalContourTree.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalContourTree.h @@ -1054,6 +1054,13 @@ void HierarchicalContourTree::AddToVTKMDataSet(vtkm::cont::DataSet& d ds.AddField(firstSupernodePerIterationOffsetsField); // TODO/FIXME: It seems we may only need the counts for the first iteration, so check, which // information we actually need. + // Add the number of rounds as an array of length 1 + vtkm::cont::ArrayHandle tempNumRounds; + tempNumRounds.Allocate(1); + vtkm::worklet::contourtree_augmented::IdArraySetValue(0, this->NumRounds, tempNumRounds); + vtkm::cont::Field numRoundsField( + "NumRounds", vtkm::cont::Field::Association::WholeDataSet, tempNumRounds); + ds.AddField(numRoundsField); } } // namespace contourtree_distributed diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalHyperSweeper.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalHyperSweeper.h index da7126a9b..22a20f6c2 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalHyperSweeper.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalHyperSweeper.h @@ -91,7 +91,9 @@ #include #include #include +#include #include +#include namespace vtkm @@ -530,6 +532,9 @@ void HierarchicalHyperSweeper::ComputeSupe } // scope ComputeSuperarcTransferWeightsWorklet // 5. Now we need to sort the transfer targets into contiguous segments + // TODO / WARNING 11/07/2023 + // We have now got a flag of ATTACHMENT_POINT_TRANSFER whose effect is to separate out transfers to + // the superarc from transfers to the supernode { // create view of superSortPermute[firstSupernode, lastSupernode) for sorting vtkm::cont::ArrayHandleView @@ -578,6 +583,16 @@ void HierarchicalHyperSweeper::TransferWei vtkm::Id firstSupernode, vtkm::Id lastSupernode) { // TransferWeights() + // TODO / WARNING 11/07/2023 + // This code was originally written on the assumption that the hierarchical tree had been augmented by the attachment points. + // As a result, it assumed that no attachment points remained. + // It is now being used for partially augmented versions due to pre-simplification, for which the correct treatment is to + // transfer not as dependent weight, but as intrinsic weight. Note that this ONLY applies to attachment points: if the + // subtree attaches at a proper supernode in the ancestor level, it should still be treated as dependent weight. The logic + // behind this is that an attachment point is regular with respect to the superarc along which it inserts. Adding it as + // dependent weight means that it is treated as *OUTSIDE* the superarc in the reverse sweep (or equivalent computation) + // Treating it as dependent weight means that both ends of the superarc end up with the correct value. + // 7. Now perform a segmented prefix sum vtkm::Id numSupernodesToProcess = lastSupernode - firstSupernode; // Same as std::partial_sum(valuePrefixSum.begin() + firstSupernode, valuePrefixSum.begin() + lastSupernode, valuePrefixSum.begin() + firstSupernode); @@ -602,6 +617,12 @@ void HierarchicalHyperSweeper::TransferWei // 7a. and 7b. { + // TODO / WARNING 11/07/2023 + // Before dealing with attachment points, we just transferred by segments. We now have + // the possibility of transferring some weight at an attachment point, + // and some not. To avoid write conflicts, we treat this as two passes: one for attachment + // points, one for all others. This means duplicating 7a/7b, sadly. + // 7a. Find the RHE of each group and transfer the prefix sum weight // Note that we do not compute the transfer weight separately, we add it in place instead // Instantiate the worklet @@ -634,7 +655,7 @@ void HierarchicalHyperSweeper::TransferWei auto valuePrefixSumPreviousValueView = vtkm::cont::make_ArrayHandleView( this->ValuePrefixSum, firstSupernode, numSupernodesToProcess - 1); - // 7b. Now find the LHE of each group and subtract out the prior weight + // 7b (non-attachment). Now find the LHE of each group and subtract out the prior weight. vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper:: TransferWeightsUpdateLHEWorklet transferWeightsUpdateLHEWorklet; this->Invoke(transferWeightsUpdateLHEWorklet, @@ -643,6 +664,43 @@ void HierarchicalHyperSweeper::TransferWei valuePrefixSumPreviousValueView, this->DependentValues); } + + // 7a (attachment). Find the RHE of each group and transfer the prefix sum weight + // Note that we do not compute the transfer weight separately, we add it in place instead + { + auto supernodeIndex = + vtkm::cont::make_ArrayHandleCounting(firstSupernode, vtkm::Id{ 1 }, numSupernodesToProcess); + auto valuePrefixSumView = vtkm::cont::make_ArrayHandleView( + this->ValuePrefixSum, firstSupernode, numSupernodesToProcess); + + vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper:: + TransferWeightsUpdateRHEWorkletRound2 transferWeightsUpdateRHEWorkletRound2(lastSupernode); + // Invoke the worklet + this->Invoke(transferWeightsUpdateRHEWorkletRound2, // worklet + supernodeIndex, // input counting array [firstSupernode, lastSupernode) + this->SortedTransferTarget, + valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode) + this->IntrinsicValues, + this->DependentValues); + } + // 7b (i). Now find the LHE of each group and subtract out the prior weight. + { + auto sortedTransferTargetView = vtkm::cont::make_ArrayHandleView( + this->SortedTransferTarget, firstSupernode + 1, numSupernodesToProcess - 1); + auto sortedTransferTargetShiftedView = vtkm::cont::make_ArrayHandleView( + this->SortedTransferTarget, firstSupernode, numSupernodesToProcess - 1); + auto valuePrefixSumPreviousValueView = vtkm::cont::make_ArrayHandleView( + this->ValuePrefixSum, firstSupernode, numSupernodesToProcess - 1); + + vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper:: + TransferWeightsUpdateLHEWorkletRound2 transferWeightsUpdateLHEWorkletRound2; + this->Invoke(transferWeightsUpdateLHEWorkletRound2, + sortedTransferTargetView, + sortedTransferTargetShiftedView, + valuePrefixSumPreviousValueView, + this->IntrinsicValues, + this->DependentValues); + } } // TransferWeights() diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CMakeLists.txt b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CMakeLists.txt index 56a3b61c1..79a50d0a8 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CMakeLists.txt @@ -23,6 +23,9 @@ set(headers AttachmentAndSupernodeComparator.h ResizeArraysBuildNewSupernodeIdsWorklet.h CreateSuperarcsWorklet.h + CreateSuperarcsData.h + CreateSuperarcsSetFirstSupernodePerIterationWorklet.h + CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h HierarchicalAugmenterInOutData.h ) diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsData.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsData.h new file mode 100644 index 000000000..3881a2ec5 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsData.h @@ -0,0 +1,169 @@ +//============================================================================ +// 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) +//============================================================================== + +// This header contains an Execution Object used to pass a arrays to the +// CreateSuperarcsWorklet to overcome the limitation of 20 input parameters for a worklet + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_superarcs_data_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_superarcs_data_h + +#include + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_augmenter +{ + + +class CreateSuperarcsData +{ +public: + // Sort indicies types + using IndicesPortalType = vtkm::worklet::contourtree_augmented::IdArrayType::ReadPortalType; + + VTKM_EXEC_CONT + CreateSuperarcsData() {} + + VTKM_CONT + CreateSuperarcsData( + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHyperparents, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichRound, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichIteration, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeSuperarcs, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHypernodes, + const vtkm::worklet::contourtree_augmented::IdArrayType& superparentSet, + const vtkm::worklet::contourtree_augmented::IdArrayType& newSupernodeIds, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) + { + this->BaseTreeHyperparents = baseTreeHyperparents.PrepareForInput(device, token); + this->BaseTreeWhichRound = baseTreeWhichRound.PrepareForInput(device, token); + this->BaseTreeWhichIteration = baseTreeWhichIteration.PrepareForInput(device, token); + this->BaseTreeSuperarcs = baseTreeSuperarcs.PrepareForInput(device, token); + this->BaseTreeHypernodes = baseTreeHypernodes.PrepareForInput(device, token); + this->SuperparentSet = superparentSet.PrepareForInput(device, token); + this->NewSupernodeIds = newSupernodeIds.PrepareForInput(device, token); + } + +public: + IndicesPortalType BaseTreeHyperparents; + IndicesPortalType BaseTreeWhichRound; + IndicesPortalType BaseTreeWhichIteration; + IndicesPortalType BaseTreeSuperarcs; + IndicesPortalType BaseTreeHypernodes; + IndicesPortalType SuperparentSet; + IndicesPortalType NewSupernodeIds; +}; + + +class CreateSuperarcsDataExec : public vtkm::cont::ExecutionObjectBase +{ +public: + VTKM_EXEC_CONT + CreateSuperarcsDataExec( + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHyperparents, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichRound, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichIteration, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeSuperarcs, + const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHypernodes, + const vtkm::worklet::contourtree_augmented::IdArrayType& superparentSet, + const vtkm::worklet::contourtree_augmented::IdArrayType& newSupernodeIds) + : BaseTreeHyperparents(baseTreeHyperparents) + , BaseTreeWhichRound(baseTreeWhichRound) + , BaseTreeWhichIteration(baseTreeWhichIteration) + , BaseTreeSuperarcs(baseTreeSuperarcs) + , BaseTreeHypernodes(baseTreeHypernodes) + , SuperparentSet(superparentSet) + , NewSupernodeIds(newSupernodeIds) + { + } + + VTKM_CONT + CreateSuperarcsData PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const + { + return CreateSuperarcsData(BaseTreeHyperparents, + BaseTreeWhichRound, + BaseTreeWhichIteration, + BaseTreeSuperarcs, + BaseTreeHypernodes, + SuperparentSet, + NewSupernodeIds, + device, + token); + } + +private: + // Whole array data used from the BaseTree in CreateSuperarcsWorklet + const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeHyperparents; + const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeWhichRound; + const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeWhichIteration; + const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeSuperarcs; + const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeHypernodes; + const vtkm::worklet::contourtree_augmented::IdArrayType& SuperparentSet; + const vtkm::worklet::contourtree_augmented::IdArrayType& NewSupernodeIds; +}; + + + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsSetFirstSupernodePerIterationWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsSetFirstSupernodePerIterationWorklet.h new file mode 100644 index 000000000..134555100 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsSetFirstSupernodePerIterationWorklet.h @@ -0,0 +1,148 @@ +//============================================================================ +// 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_augmenter_create_auperarcs_set_first_supernode_per_iteration_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_set_first_supernode_per_iteration_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_augmenter +{ + +/// Worklet used in HierarchicalAugmenter::UpdateHyperstructure to set the hyperarcs and hypernodes +class CreateSuperarcsSetFirstSupernodePerIterationWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + using ControlSignature = void(FieldIn supernodeIndex, + WholeArrayIn augmentedTreeWhichIteration, + WholeArrayInOut augmentedTreeFirstSupernodePerIteration); + using ExecutionSignature = void(_1, _2, _3); + using InputDomain = _1; + + // Default Constructor + VTKM_EXEC_CONT + CreateSuperarcsSetFirstSupernodePerIterationWorklet(vtkm::Id numSupernodesAlready) + : NumSupernodesAlready(numSupernodesAlready) + { + } + + template + VTKM_EXEC void operator()( + const vtkm::Id& supernode, // index in supernodeSorter + // const vtkm::Id& supernodeSetindex, // supernodeSorter[supernode] + const InFieldPortalType& augmentedTreeWhichIterationPortal, + const InOutFieldPortalType& augmentedTreeFirstSupernodePerIterationPortal) const + { // operator()() + // per supernode in the set + // retrieve the index from the sorting index array (Done on input)(NOT USED) + // indexType supernodeSetIndex = supernodeSorter[supernode]; + + // work out the new supernode ID + vtkm::Id newSupernodeId = this->NumSupernodesAlready + supernode; + + // The 0th element sets the first element in the zeroth iteration + if (supernode == 0) + { + augmentedTreeFirstSupernodePerIterationPortal.Set(0, newSupernodeId); + } + // otherwise, mismatch to the left identifies a new iteration + else + { + if (vtkm::worklet::contourtree_augmented::MaskedIndex( + augmentedTreeWhichIterationPortal.Get(newSupernodeId)) != + vtkm::worklet::contourtree_augmented::MaskedIndex( + augmentedTreeWhichIterationPortal.Get(newSupernodeId - 1))) + { // start of segment + augmentedTreeFirstSupernodePerIterationPortal.Set( + vtkm::worklet::contourtree_augmented::MaskedIndex( + augmentedTreeWhichIterationPortal.Get(newSupernodeId)), + newSupernodeId); + } // start of segmen + } + + /* + #pragma omp parallel for + for (indexType supernode = 0; supernode < supernodeSorter.size(); supernode++) + { // per supernode in the set + // retrieve the index from the sorting index array + indexType supernodeSetIndex = supernodeSorter[supernode]; + + // work out the new supernode ID + indexType newSupernodeID = nSupernodesAlready + supernode; + + // The 0th element sets the first element in the zeroth iteration + if (supernode == 0) + augmentedTree->firstSupernodePerIteration[roundNo][0] = newSupernodeID; + // otherwise, mismatch to the left identifies a new iteration + else + { + if (augmentedTree->whichIteration[newSupernodeID] != augmentedTree->whichIteration[newSupernodeID-1]) + augmentedTree->firstSupernodePerIteration[roundNo][maskedIndex(augmentedTree->whichIteration[newSupernodeID])] = newSupernodeID; + } + } // per supernode in the set + */ + + } // operator()() + +private: + vtkm::Id NumSupernodesAlready; + + +}; // CreateSuperarcsSetFirstSupernodePerIterationWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h new file mode 100644 index 000000000..386c74399 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h @@ -0,0 +1,104 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +//============================================================================ +// Copyright (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// (3) Neither the name of the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +// OF THE POSSIBILITY OF SUCH DAMAGE. +// +//============================================================================= +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_update_first_supernode_per_iteration_worklet_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_update_first_supernode_per_iteration_worklet_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_augmenter +{ + +/// Worklet used in HierarchicalAugmenter::UpdateHyperstructure to set the hyperarcs and hypernodes +class CreateSuperarcsUpdateFirstSupernodePerIterationWorklet : public vtkm::worklet::WorkletMapField +{ +public: + /// Control signature for the worklet + using ControlSignature = void(FieldIn indexArray, + WholeArrayInOut augmentedTreeFirstSupernodePerIteration); + using ExecutionSignature = void(_1, _2); + using InputDomain = _1; + + // Default Constructor + VTKM_EXEC_CONT + CreateSuperarcsUpdateFirstSupernodePerIterationWorklet() {} + + template + VTKM_EXEC void operator()( + const vtkm::Id& iteration, + const InOutFieldPortalType& augmentedTreeFirstSupernodePerIterationPortal) const + { // operator()() + if (augmentedTreeFirstSupernodePerIterationPortal.Get(iteration) == 0) + { + augmentedTreeFirstSupernodePerIterationPortal.Set( + iteration, augmentedTreeFirstSupernodePerIterationPortal.Get(iteration + 1)); + } + + /* + for (indexType iteration = 1; iteration < augmentedTree->nIterations[roundNo]; ++iteration) + { + if (augmentedTree->firstSupernodePerIteration[roundNo][iteration] == 0) + { + augmentedTree->firstSupernodePerIteration[roundNo][iteration] = + augmentedTree->firstSupernodePerIteration[roundNo][iteration+1]; + } + } + */ + } // operator()() +}; // CreateSuperarcsUpdateFirstSupernodePerIterationWorklet + +} // namespace hierarchical_augmenter +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsWorklet.h index 32bcc6037..28097d36c 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsWorklet.h @@ -61,6 +61,7 @@ namespace hierarchical_augmenter /// Worklet used to implement the main part of HierarchicalAugmenter::CreateSuperarcs /// Connect superarcs for the level & set hyperparents & superchildren count, whichRound, /// whichIteration, super2hypernode +template class CreateSuperarcsWorklet : public vtkm::worklet::WorkletMapField { public: @@ -69,342 +70,559 @@ public: /// Control signature for the worklet /// @param[in] supernodeSorter input domain. We need access to InputIndex and InputIndex+1, /// therefore this is a WholeArrayIn transfer. - /// @param[in] superparentSet WholeArrayIn because we need access to superparentSet[supernodeSorter[InputIndex]] - /// and superparentSet[supernodeSorter[InputIndex+1]]. - /// @param[in] baseTreeSuperarcs WholeArrayIn because we need access to baseTreeSuperarcsPortal.Get(superparentOldSuperId) - /// While this could be done with fancy array magic, it would require a sequence of multiple - /// fancy arrays and would likely not be cheaper then computing things in the worklet. - /// @param[in] newSupernodeIds WholeArrayIn because we need to access newSupernodeIdsPortal.Get(oldTargetSuperId) - /// where oldTargetSuperId is the unmasked baseTreeSuperarcsPortal.Get(superparentOldSuperId) - /// @param[in] baseTreeSupernodes WholeArrayIn because we need to access baseTreeSupernodesPortal.Get(superparentOldSuperId); - /// @param[in] baseTreeRegularNodeGlobalIds WholeArrayIn because we need to access - /// baseTreeRegularNodeGlobalIdsPortal.Get(superparentOldSuperId); - /// @param[in] globalRegularIdSet FieldInd. Permute globalRegularIdSet with supernodeSorter in order to allow this to be a FieldIn. - /// @param[in] baseTreeSuper2Hypernode WholeArrayIn because we need to access - /// baseTreeSuper2HypernodePortal.Get(superparentOldSuperId) - /// @param[in] baseTreeWhichIteration WholeArrayIn because we need to access baseTreeWhichIterationPortal.Get(superparentOldSuperId) - /// and baseTreeWhichIterationPortal.Get(superparentOldSuperId+1) - /// @param[in] augmentedTreeSuperarcsView output view of this->AugmentedTree->Superarcs with + /// @param[in] supernodeIdSetPermuted Field in of supernodeIdSet permuted by the supernodeSorter array to + /// allow us to use FieldIn + /// @param[in] baseTreeHyperparents whole array input of the BaseTree->Hyperparents + /// @param[in] baseTreeSuper2HypernodePermuted baseTreeSuper2Hypernode permuted by supernodeIdSetPermuted in order + /// to extract baseTree->super2hypernode[oldSupernodeID]; + /// @param[in] baseTreeWhichRound + /// @param[in] baseTreeWhichIteration + /// @param[in] globalRegularIdSetPermuted Field in of globalRegularIdSet permuted by supernodeSorter array to + /// allow use of FieldIn + /// @param[in] dataValueSetPermuted Field in of dataValyeSet permuted by supernodeSorter array to + /// allow use of FieldIn + /// @param[in] baseTreeSuperarcs BaseTree->Superarcs + /// @param[in] oldSuperFrom Permuted baseTree->SuperParents[ baseTree->Supernodes[ supernodeIdSetPermuted ] ] + /// @param[in] baseTreeHypernodes + /// @param[out, in] augmentedTreeSupernodes this->AugmentedTree->Supernodes array + /// @param[out] augmentedTreeSuperarcsView output view of this->AugmentedTree->Superarcs with /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs, /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). /// By using this view allows us to do this one as a FieldOut and it effectively the /// same as accessing the array at the newSuppernodeId location. - /// @param[in] augmentedTreeFirstSupernodePerIteration WholeArrayInOut because we need to update multiple locations. - /// In is used to preseve original values. Set to augmentedTree->firstSupernodePerIteration[roundNumber]. - /// @param[in] augmentedTreeSuper2hypernode FieldOut. Output view of this->AugmentedTree->Super2Hypernode + /// @param[out] augmentedTreeHyperparentsView output view of this->AugmentedTree->Hyperparents with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hyperparents, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newSuppernodeId location. + /// @param[out] augmentedTreeSuper2HypernodeView output view of this->AugmentedTree->Super2Hypernode with /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Super2Hypernode, /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). /// By using this view allows us to do this one as a FieldOut and it effectively the /// same as accessing the array at the newSuppernodeId location. + /// @param[out] augmentedTreeWhichRoundView output view of this->AugmentedTree->WhichRound with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->WhichRound, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newSuppernodeId location. + /// @param[out] augmentedTreeWhichIterationView output view of this->AugmentedTree->WhichIteration with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->WhichIteration, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newSuppernodeId location. + /// @param[out] augmentedTreeRegularNodeGlobalIdsView output view of this->AugmentedTree->RegularNodeGlobalIds with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->RegularNodeGlobalIds, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newRegularId location. + /// @param[out] augmentedTreeDataValuesView output view of this->AugmentedTree->DataValues with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->DataValues, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newRegularId location. + /// @param[out] augmentedTreeRegular2SupernodeView output view of this->AugmentedTree->Regular2Supernode with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Regular2Supernode, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newRegularId location. + /// @param[out] augmentedTreeSuperparentsView output view of this->AugmentedTree->Superparents with + /// vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superparents, + /// numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()). + /// By using this view allows us to do this one as a FieldOut and it effectively the + /// same as accessing the array at the newRegularId location. + using ControlSignature = void( + // Inputs WholeArrayIn supernodeSorter, - WholeArrayIn superparentSet, // input - WholeArrayIn baseTreeSuperarcs, // input - WholeArrayIn newSupernodeIds, // input - WholeArrayIn baseTreeSupernodes, // input - WholeArrayIn baseTreeRegularNodeGlobalIds, // input - FieldIn globalRegularIdSet, // input - WholeArrayIn baseTreeSuper2Hypernode, // input - WholeArrayIn baseTreeWhichIteration, // input - FieldOut augmentedTreeSuperarcsView, // output - WholeArrayInOut augmentedTreeFirstSupernodePerIteration, // input/output - FieldOut augmentedTreeSuper2hypernode // ouput - ); - using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12); + FieldIn supernodeIdSetPermuted, + FieldIn baseTreeSuper2HypernodePermuted, + FieldIn globalRegularIdSetPermuted, + FieldIn dataValueSetPermuted, + FieldIn + oldSuperFrom, // baseTree->SuperParents[ baseTree->Supernodes[ supernodeIdSetPermuted ] ] + ExecObject findSuperArcForUnknownNode, + ExecObject createSuperarcsData, + // Outputs + WholeArrayInOut augmentedTreeSupernodes, + FieldOut augmentedTreeSuperarcsView, + FieldOut augmentedTreeHyperparensView, + FieldOut augmentedTreeSuper2Hypernode, + FieldOut augmentedTreeWhichRoundView, + FieldOut augmentedTreeWhichIterationView, + FieldOut augmentedTreeRegularNodeGlobalIdsView, + FieldOut augmentedTreeDataValuesView, + FieldOut augmentedTreeRegular2SupernodeView, + FieldOut augmentedTreeSuperparentsViews); + using ExecutionSignature = void(InputIndex, + _1, + _2, + _3, + _4, + _5, + _6, + _7, + _8, + _9, + _10, + _11, + _12, + _13, + _14, + _15, + _16, + _17, + _18); + // using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, _16, _17, _18, _19, _20); using InputDomain = _1; /// Default Constructor /// @param[in] numSupernodesAlready Set to vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstSupernodePerIteration[roundNumber]); /// @param[in] baseTreeNumRounds Set to this->BaseTree->NumRounds - /// @param[in] augmentedTreeNumIterations Set to vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations); - /// @param[in] roundNumber Set the current round - /// @param[in] numAugmentedTreeSupernodes Set to augmentedTreeSupernodes this->AugmentedTree->Supernodes.GetNumberOfValues(); + /// @param[in] numInsertedSupernodes Set to numInsertedSupernodes VTKM_EXEC_CONT CreateSuperarcsWorklet(const vtkm::Id& numSupernodesAlready, const vtkm::Id& baseTreeNumRounds, - const vtkm::Id& augmentedTreeNumIterations, - const vtkm::Id& roundNumber, - const vtkm::Id& numAugmentedTreeSupernodes) + const vtkm::Id& numInsertedSupernodes, + const vtkm::Id& roundNo) : NumSupernodesAlready(numSupernodesAlready) , BaseTreeNumRounds(baseTreeNumRounds) - , AugmentedTreeNumIterations(augmentedTreeNumIterations) - , RoundNumber(roundNumber) - , NumAugmentedTreeSupernodes(numAugmentedTreeSupernodes) + , NumInsertedSupernodes(numInsertedSupernodes) + , RoundNo(roundNo) { } /// operator() of the workelt - template + template + //typename InOutDataFieldPortalType> + //typename InOutFieldPortalType, + //typename ExecObjectTypeData, + //typename ExecObjType> VTKM_EXEC void operator()( + // Inputs const vtkm::Id& supernode, // InputIndex of supernodeSorter const InFieldPortalType& supernodeSorterPortal, - const InFieldPortalType& superparentSetPortal, - const InFieldPortalType& baseTreeSuperarcsPortal, - const InFieldPortalType& newSupernodeIdsPortal, - const InFieldPortalType& baseTreeSupernodesPortal, - const InFieldPortalType& baseTreeRegularNodeGlobalIdsPortal, - const vtkm::Id& globalRegularIdSetValue, - const InFieldPortalType& baseTreeSuper2HypernodePortal, - const InFieldPortalType& baseTreeWhichIterationPortal, - vtkm::Id& augmentedTreeSuperarcsValue, // same as augmentedTree->superarcs[newSupernodeId] - const InOutFieldPortalType& - augmentedTreeFirstSupernodePerIterationPortal, // augmentedTree->firstSupernodePerIteration[roundNumber] - vtkm::Id& augmentedTreeSuper2hypernodeValue) const + const vtkm::Id& oldSupernodeId, // supernodeIdSet[supernodeSorterPortal.Get(supernode)] + const vtkm::Id& + baseTreeSuper2HypernodeOldSupernodeIdValue, // baseTree->super2hypernode[oldSupernodeID]; + const vtkm::Id& + globalRegularIdSetValue, // globalRegularIdSet[supernodeSorterPortal.Get(supernode)]]; + const FieldType& dataValueSetValue, // dataValueSet[supernodeSorterPortal.Get(supernode)]]; + const vtkm::Id& + oldSuperFromValue, // baseTree->SuperParents[ baseTree->Supernodes[ supernodeIdSetPermuted ] ] + const ExecObjType& + findSuperArcForUnknownNode, // Execution object to call FindSuperArcForUnknownNode + const ExecObjectTypeData& + createSuperarcsData, // Execution object of collect BaseTree data array + // Outputs + const InOutFieldPortalType& augmentedTreeSupernodesPortal, + vtkm::Id& augmentedTreeSuperarcsValue, // set value for AugmentedTree->Auperarcs[newSupernodeId] + vtkm::Id& + augmentedTreeHyperparentsValue, // set value for AugmentedTree->Hyperparents[newSupernodeId] + vtkm::Id& + augmentedTreeSuper2HypernodeValue, // set value for AugmentedTree->Super2Hypernode[newSupernodeId] + vtkm::Id& augmentedTreeWhichRoundValue, // AugmentedTree->WhichRound[newSupernodeId] + vtkm::Id& augmentedTreeWhichIterationValue, // AugmentedTree->WhichIteration[newSupernodeId] + vtkm::Id& + augmentedTreeRegularNodeGlobalIdsValue, // AugmentedTree->RegularNodeGlobalIds[newRegularID] + FieldType& augmentedTreeDataValuesValue, // AugmentedTree->DataValues[newRegularID] + vtkm::Id& augmentedTreeRegular2SupernodeValue, // AugmentedTree->Regular2Supernode[newRegularID] + vtkm::Id& augmentedTreeSuperparentsValue // AugmentedTree->Superparents[newRegularID] + ) const { // per supernode in the set // retrieve the index from the sorting index array vtkm::Id supernodeSetIndex = supernodeSorterPortal.Get(supernode); - // work out the new supernode Id. We have this defined on the outside as a fancy array handle, - // however, using the fancy handle here would not really make a performance differnce and - // computing it here is more readable + // work out the new supernode ID vtkm::Id newSupernodeId = this->NumSupernodesAlready + supernode; - // NOTE: The newRegularId is no longer needed here since all parts - // that used it in the worklet have been moved outside - // vtkm::Id newRegularId = newSupernodeId; + // and the old supernode ID + // vtkm::Id oldSupernodeId = supernodeIDSet[supernodeSetIndex]; Extracted on call and provides as input - // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs - // // setting the supernode's regular Id is now trivial - // augmentedTreeSupernodesPortal.Set(newSupernodeId, newRegularId); + // At all levels above 0, we used to keep regular vertices in case they are attachment points. + // After augmentation, we don't need to. + // Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes + // In order to avoid confusion, we will copy the ID into a separate variable + vtkm::Id newRegularId = newSupernodeId; - // retrieve the ascending flag from the superparent - vtkm::Id superparentSetVal = superparentSetPortal.Get(supernodeSetIndex); - // get the ascending flag from the parent - bool superarcAscends = vtkm::worklet::contourtree_augmented::IsAscending(superparentSetVal); - // strip the ascending flag from the superparent. - vtkm::Id superparentOldSuperId = - vtkm::worklet::contourtree_augmented::MaskedIndex(superparentSetVal); + // setting the supernode's regular ID is now trivial + augmentedTreeSupernodesPortal.Set(newSupernodeId, newRegularId); - // setting the superarc is done the usual way. Our sort routine has ended up - // with the supernodes arranged in either ascending or descending order - // inwards along the parent superarc (as expressed by the superparent Id). - // Each superarc except the last in the segment points to the next one: - // the last one points to the target of the original superarc. - // first test to see if we're the last in the array - if (supernode == supernodeSorterPortal.GetNumberOfValues() - 1) - { // last in the array - // special case for root of entire tree at end of top level - if (RoundNumber == this->BaseTreeNumRounds) + // retrieve the old superID of the superparent. This is slightly tricky, as we have four classes of supernodes: + // 1. the root of the entire tree + // 2. attachment points not being inserted. In this case, the supernode ID is stored in the superparentSet + // array, not the superparent for insertion purposes + // 3. attachment points being inserted. In this case, the superparent is stored in the superparentSet array + // 4. "ordinary" supernodes, where the superparent is the same as the supernode ID anyway + // + // Note that an attachment point gets inserted into a parent superarc. But the attachment point itself has + // a NULL superarc, because it's only a virtual insertion. + // This means that such an attachment superarc cannot be the superparent of any other attachment point + // It is therefore reasonable to deal with 1. & 2 separately. 3. & 4. then combine together + + // first we test for the root of the tree + if ((this->RoundNo == BaseTreeNumRounds) && + (supernode == supernodeSorterPortal.GetNumberOfValues() - 1)) + { // root of the tree + // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the root is in every tree + // set the super arrays + augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + // hyperstructure carries over, so we use the same hyperparent as before + augmentedTreeHyperparentsValue = createSuperarcsData.BaseTreeHyperparents.Get(oldSupernodeId); + // and set the hypernode ID + augmentedTreeSuper2HypernodeValue = baseTreeSuper2HypernodeOldSupernodeIdValue; + // and the round and iteration + augmentedTreeWhichRoundValue = createSuperarcsData.BaseTreeWhichRound.Get(oldSupernodeId); + augmentedTreeWhichIterationValue = + createSuperarcsData.BaseTreeWhichIteration.Get(oldSupernodeId); + // and set the relevant regular arrays + augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue; + augmentedTreeDataValuesValue = dataValueSetValue; + // for the root, these always point to itself + augmentedTreeRegular2SupernodeValue = newSupernodeId; + augmentedTreeSuperparentsValue = newSupernodeId; + } // root of the tree + // now deal with unsimplified attachment points, which we can identify because they were in the "kept" batch, not the "inserted" batch, + // and this is given away by the index into the set of supernodes to be added + // and the fact that the superarc is NO_SUCH_ELEMENT + else if ((supernodeSetIndex >= this->NumInsertedSupernodes) && + (vtkm::worklet::contourtree_augmented::NoSuchElement( + createSuperarcsData.BaseTreeSuperarcs.Get(oldSupernodeId)))) + { // preserved attachment point + // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the supernode came from this block originally + // set the superarc to NO_SUCH_ELEMENT, as before + augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + // hyperstructure carries over, so we use the same hyperparent as before + augmentedTreeHyperparentsValue = createSuperarcsData.BaseTreeHyperparents.Get(oldSupernodeId); + // attachment points are never hypernodes anyway, so set it directly + augmentedTreeSuper2HypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + // and the round and iteration + augmentedTreeWhichRoundValue = createSuperarcsData.BaseTreeWhichRound.Get(oldSupernodeId); + augmentedTreeWhichIterationValue = + createSuperarcsData.BaseTreeWhichIteration.Get(oldSupernodeId); + // and set the relevant regular arrays + augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue; + augmentedTreeDataValuesValue = dataValueSetValue; + // for a preserved attachment point, this always points to itself + augmentedTreeRegular2SupernodeValue = newSupernodeId; + // the superparent is the tricky one, as the old one may have been broken up by insertions at a higher level + + // Here, what we need to do is a search in the augmented tree to find which superarc to attach to. This is necessary + // because the old superarc it attached to may have been broken up. + // We are guaranteed that there is one, and that it only uses the higher levels of the augmented tree, + // so the fact that we are partially constructed doesn't get in the way. To do this, we need supernodes + // known to be in the higher level that are above and below the supernode. + // Since the point was an attachment point in the base tree, that means that there is a higher round superarc + // it inserts into. Moreover, the algorithm ALWAYS inserts a supernode at or above its original round, so + // we can guarantee that both ends of the parent are in the higher levels. Which means we only need to work + // out which end is higher. + + // oldSuperFrom is provided as input and extracted as FieldIn on call. oldRegularId is not needed here + // indexType oldRegularID = baseTree->supernodes[oldSupernodeID]; + // indexType oldSuperFrom = baseTree->superparents[oldRegularID]; + // indexType oldSuperTo = baseTree->superarcs[oldSuperFrom]; + vtkm::Id oldSuperToValue = createSuperarcsData.BaseTreeSuperarcs.Get(oldSuperFromValue); + + // retrieve the ascending flag + bool ascendingSuperarc = vtkm::worklet::contourtree_augmented::IsAscending(oldSuperToValue); + // and mask out the flags + vtkm::Id oldSuperToMaskedIndex = + vtkm::worklet::contourtree_augmented::MaskedIndex(oldSuperToValue); + + // since we haven't set up the regular search array yet, we can't use that + // instead, we know that the two supernodes must be in the new tree, so we retrieve their new super IDs + // and convert them to regular + + // retrieve their new super IDs + vtkm::Id newSuperFrom = createSuperarcsData.NewSupernodeIds.Get(oldSuperFromValue); + vtkm::Id newSuperTo = createSuperarcsData.NewSupernodeIds.Get(oldSuperToMaskedIndex); + + // convert to regular IDs (which is what the FindSuperArcForUnknownNode() routine assumes) + vtkm::Id newRegularFrom = augmentedTreeSupernodesPortal.Get(newSuperFrom); + vtkm::Id newRegularTo = augmentedTreeSupernodesPortal.Get(newSuperTo); + + // the new superparent after the search + vtkm::Id newSuperparentId = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + + // depending on the ascending flag + if (ascendingSuperarc) { - augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + newSuperparentId = findSuperArcForUnknownNode.FindSuperArcForUnknownNode( + globalRegularIdSetValue, dataValueSetValue, newRegularTo, newRegularFrom); } else - { // not the tree root - // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag) - vtkm::Id oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex( - baseTreeSuperarcsPortal.Get(superparentOldSuperId)); - // convert to a new supernode Id - vtkm::Id newTargetSuperId = newSupernodeIdsPortal.Get(oldTargetSuperId); - // add the ascending flag back in and store in the array - augmentedTreeSuperarcsValue = newTargetSuperId | - (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00); - } // not the tree root - // since there's an extra entry in the firstSupernode array as a sentinel, set it - augmentedTreeFirstSupernodePerIterationPortal.Set(this->AugmentedTreeNumIterations, - NumAugmentedTreeSupernodes); - } // last in the array - else if (superparentOldSuperId != - vtkm::worklet::contourtree_augmented::MaskedIndex( - superparentSetPortal.Get(supernodeSorterPortal.Get(supernode + 1)))) - { // last in the segment - // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag) - vtkm::Id oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex( - baseTreeSuperarcsPortal.Get(superparentOldSuperId)); - // convert to a new supernode Id - vtkm::Id newTargetSuperId = newSupernodeIdsPortal.Get(oldTargetSuperId); - // add the ascending flag back in and store in the array - augmentedTreeSuperarcsValue = newTargetSuperId | - (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00); + { + newSuperparentId = findSuperArcForUnknownNode.FindSuperArcForUnknownNode( + globalRegularIdSetValue, dataValueSetValue, newRegularFrom, newRegularTo); + } - // since we're the last in the segment, we check to see if we are at the end of an iteration - vtkm::Id iterationNumber = vtkm::worklet::contourtree_augmented::MaskedIndex( - baseTreeWhichIterationPortal.Get(superparentOldSuperId)); - vtkm::Id iterationNumberOfNext = vtkm::worklet::contourtree_augmented::MaskedIndex( - baseTreeWhichIterationPortal.Get(superparentOldSuperId + 1)); + // attachment points use the superparent to store the superarc they insert onto + augmentedTreeSuperparentsValue = newSuperparentId; - if (iterationNumber != iterationNumberOfNext) - { // boundary of iterations - // If so, we set the "firstSupernodePerIteration" for the next - augmentedTreeFirstSupernodePerIterationPortal.Set(iterationNumberOfNext, - newSupernodeId + 1); - } // boundary of iterations - } // last in the segment + } // preserved attachment point else - { // not last in the segment - // the target is always the next one, so just store it with the ascending flag - augmentedTreeSuperarcsValue = (newSupernodeId + 1) | - (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00); - } // not last in the segment + { // raised attachment point or "ordinary" supernodes + // Since all of the superparents must be in the base tree, we can now retrieve the target + vtkm::Id superparentOldSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex( + createSuperarcsData.SuperparentSet.Get(supernodeSetIndex)); + vtkm::Id oldTargetSuperId = createSuperarcsData.BaseTreeSuperarcs.Get(superparentOldSuperId); + // and break it into a target and flags + bool ascendingSuperarc = vtkm::worklet::contourtree_augmented::IsAscending(oldTargetSuperId); + // NOTE: if the target was NO_SUCH_ELEMENT, this will hold 0 + oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex(oldTargetSuperId); - // set the first supernode in the first iteration to the beginning of the round - augmentedTreeFirstSupernodePerIterationPortal.Set(0, this->NumSupernodesAlready); + // and another boolean for whether we are the last element in a segment + bool isLastInSegment = false; - - // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs - // // setting the hyperparent is straightforward since the hyperstructure is preserved - // // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that - // augmentedTreeHyperparentsPortal.Set(newSupernodeId, baseTreeHyperparentsPortal.Get(superparentOldSuperId)); - - // NOTE: This part could potentially be made a separate worklet but it does not seem necessary - // similarly, the super2hypernode should carry over, but it's harder to test because of the attachment points which - // do not have valid old supernode Ids. Instead, we check their superparent's regular global Id against them: if it - // matches, then it must be the start of the superarc, in which case it does have an old Id, and we can then use the - // existing hypernode Id - vtkm::Id superparentOldRegularId = baseTreeSupernodesPortal.Get(superparentOldSuperId); - vtkm::Id superparentGlobalId = baseTreeRegularNodeGlobalIdsPortal.Get(superparentOldRegularId); - // Here: globalRegularIdSetValue is the same as globalRegularIdSetPortal.Get(supernodeSetIndex) - if (superparentGlobalId == globalRegularIdSetValue) - { - // augmentedTreeSuper2hypernodePortal.Set(newSupernodeId, baseTreeSuper2HypernodePortal.Get(superparentOldSuperId)); - augmentedTreeSuper2hypernodeValue = baseTreeSuper2HypernodePortal.Get(superparentOldSuperId); - } - else - { - // augmentedTreeSuper2hypernodePortal.Set(newSupernodeId, vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT); - augmentedTreeSuper2hypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; - } - - // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs - // // which round and iteration carry over - // augmentedTreeWhichRoundPortal.Set(newSupernodeId, baseTreeWhichRoundPortal.Get(superparentOldSuperId)); - // augmentedTreeWhichIterationPortal.Set(newSupernodeId, baseTreeWhichIterationPortal.Get(superparentOldSuperId)); - - // now we deal with the regular-sized arrays - - // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs - // // copy the global regular Id and data value - // augmentedTreeRegularNodeGlobalIdsPortal.Set(newRegularId, globalRegularIdSetPortal.Get(supernodeSetIndex)); - // augmentedTreeDataValuesPortal.Set(newRegularId, dataValueSetPortal.Get(supernodeSetIndex)); - - // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs - // // the sort order will be dealt with later - // // since all of these nodes are supernodes, they will be their own superparent, which means that: - // // a. the regular2node can be set immediately - // augmentedTreeRegular2SupernodePortal.Set(newRegularId, newSupernodeId); - // // b. as can the superparent - // augmentedTreeSuperparentsPortal.Set(newRegularId, newSupernodeId); - - // In serial this worklet implements the following operation - /* - for (vtkm::Id supernode = 0; supernode < supernodeSorter.size(); supernode++) - { // per supernode in the set - // retrieve the index from the sorting index array - vtkm::Id supernodeSetIndex = supernodeSorter[supernode]; - - // work out the new supernode ID - vtkm::Id newSupernodeID = numSupernodesAlready + supernode; - - // At all levels above 0, we used to keep regular vertices in case they are attachment points. After augmentation, we don't need to. - // Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes - // In order to avoid confusion, we will copy the ID into a separate variable - vtkm::Id newRegularID = newSupernodeID; - - // setting the supernode's regular ID is now trivial - augmentedTree->supernodes [newSupernodeID] = newRegularID; - - // retrieve the ascending flag from the superparent - bool superarcAscends = isAscending(superparentSet[supernodeSetIndex]); - - // strip the ascending flag from the superparent - vtkm::Id superparentOldSuperID = maskedIndex(superparentSet[supernodeSetIndex]); + // end of the entire array counts as last in segment + if (supernode == supernodeSorterPortal.GetNumberOfValues() - 1) + { + isLastInSegment = true; + } + // otherwise, check for a mismatch in the sorting superparent which indicates the end of a segment + else if (vtkm::worklet::contourtree_augmented::MaskedIndex( + createSuperarcsData.SuperparentSet.Get(supernodeSetIndex)) != + vtkm::worklet::contourtree_augmented::MaskedIndex( + createSuperarcsData.SuperparentSet.Get(supernodeSorterPortal.Get(supernode + 1)))) + { + isLastInSegment = true; + } // setting the superarc is done the usual way. Our sort routine has ended up with the supernodes arranged in either ascending or descending order // inwards along the parent superarc (as expressed by the superparent ID). Each superarc except the last in the segment points to the next one: // the last one points to the target of the original superarc. - // first test to see if we're the last in the array - if (supernode == supernodeSorter.size() - 1) - { // last in the array - // special case for root of entire tree at end of top level - if (roundNumber == baseTree->nRounds) - { - augmentedTree->superarcs[newSupernodeID] = NO_SUCH_ELEMENT; - } - else - { // not the tree root - // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag) - vtkm::Id oldTargetSuperID = maskedIndex(baseTree->superarcs[superparentOldSuperID]); - // convert to a new supernode ID - vtkm::Id newTargetSuperID = newSupernodeIDs[oldTargetSuperID]; - // add the ascending flag back in and store in the array - augmentedTree->superarcs[newSupernodeID] = newTargetSuperID | (superarcAscends ? IS_ASCENDING : 0x00); - } // not the tree root - // since there's an extra entry in the firstSupernode array as a sentinel, set it - augmentedTree->firstSupernodePerIteration[roundNumber][augmentedTree->nIterations[roundNumber]] = augmentedTree->supernodes.size(); - } // last in the array - else if (superparentOldSuperID != maskedIndex(superparentSet[supernodeSorter[supernode+1]])) - { // last in the segment - // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag) - vtkm::Id oldTargetSuperID = maskedIndex(baseTree->superarcs[superparentOldSuperID]); - // convert to a new supernode ID - vtkm::Id newTargetSuperID = newSupernodeIDs[oldTargetSuperID]; - // add the ascending flag back in and store in the array - augmentedTree->superarcs[newSupernodeID] = newTargetSuperID | (superarcAscends ? IS_ASCENDING : 0x00); - - // since we're the last in the segment, we check to see if we are at the end of an iteration - vtkm::Id iterationNumber = maskedIndex(baseTree->whichIteration[superparentOldSuperID]); - vtkm::Id iterationNumberOfNext = maskedIndex(baseTree->whichIteration[superparentOldSuperID + 1]); - - if (iterationNumber != iterationNumberOfNext) - { // boundary of iterations - // If so, we set the "firstSupernodePerIteration" for the next - augmentedTree->firstSupernodePerIteration[roundNumber][iterationNumberOfNext] = newSupernodeID + 1; - } // boundary of iterations - } // last in the segment + if (isLastInSegment) + { // last in segment + // we take the old target of the superarc (in old supernode IDs) and convert it to a new supernode ID + augmentedTreeSuperarcsValue = createSuperarcsData.NewSupernodeIds.Get(oldTargetSuperId) | + (ascendingSuperarc ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00); + } // last in segment else - { // not last in the segment + { // not last in segment // the target is always the next one, so just store it with the ascending flag - augmentedTree->superarcs[newSupernodeID] = (newSupernodeID+1) | (superarcAscends ? IS_ASCENDING : 0x00); - } // not last in the segment + augmentedTreeSuperarcsValue = (newSupernodeId + 1) | + (ascendingSuperarc ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00); + } // not last in segment - // set the first supernode in the first iteration to the beginning of the round - augmentedTree->firstSupernodePerIteration[roundNumber][0] = numSupernodesAlready; + // first we identify the hyperarc on which the superarc sits + // this will be visible in the old base tree, since hyperstructure carries over + vtkm::Id oldHyperparent = createSuperarcsData.BaseTreeHyperparents.Get(superparentOldSuperId); - // setting the hyperparent is straightforward since the hyperstructure is preserved - // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that - augmentedTree->hyperparents [newSupernodeID] = baseTree->hyperparents [superparentOldSuperID]; + // hyperstructure carries over, so we use the same hyperparent as the superparent + augmentedTreeHyperparentsValue = oldHyperparent; - // similarly, the super2hypernode should carry over, but it's harder to test because of the attachment points which - // do not have valid old supernode IDs. Instead, we check their superparent's regular global ID against them: if it - // matches, then it must be the start of the superarc, in which case it does have an old ID, and we can then use the - // existing hypernode ID - vtkm::Id superparentOldRegularID = baseTree->supernodes[superparentOldSuperID]; - vtkm::Id superparentGlobalID = baseTree->regularNodeGlobalIDs[superparentOldRegularID]; - if (superparentGlobalID == globalRegularIDSet[supernodeSetIndex]) + // retrieve the hyperparent's old supernode ID & convert to a new one, then test it + if (createSuperarcsData.NewSupernodeIds.Get( + createSuperarcsData.BaseTreeHypernodes.Get(oldHyperparent)) == newSupernodeId) { - augmentedTree->super2hypernode [newSupernodeID] = baseTree->super2hypernode[superparentOldSuperID]; + augmentedTreeSuper2HypernodeValue = oldHyperparent; } else { - augmentedTree->super2hypernode [newSupernodeID] = NO_SUCH_ELEMENT; + augmentedTreeSuper2HypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; } - // which round and iteration carry over - augmentedTree->whichRound [newSupernodeID] = baseTree->whichRound[superparentOldSuperID]; - augmentedTree->whichIteration [newSupernodeID] = baseTree->whichIteration[superparentOldSuperID]; + // round and iteration are set from the superparent, since we are raising to its level + augmentedTreeWhichRoundValue = + createSuperarcsData.BaseTreeWhichRound.Get(superparentOldSuperId); + augmentedTreeWhichIterationValue = + createSuperarcsData.BaseTreeWhichIteration.Get(superparentOldSuperId); + // and set the relevant regular arrays + augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue; + augmentedTreeDataValuesValue = dataValueSetValue; + // for all supernodes, this points to itself + augmentedTreeRegular2SupernodeValue = newSupernodeId; + // and since we're inserted, so does this + augmentedTreeSuperparentsValue = newSupernodeId; + } // raised attachment point or "ordinary" supernodes - // now we deal with the regular-sized arrays + /* Original PPP2 code + for (indexType supernode = 0; supernode < supernodeSorter.size(); supernode++) + { // per supernode in the set + // retrieve the index from the sorting index array + indexType supernodeSetIndex = supernodeSorter[supernode]; - // copy the global regular ID and data value - augmentedTree->regularNodeGlobalIDs [newRegularID] = globalRegularIDSet[supernodeSetIndex]; - augmentedTree->dataValues [newRegularID] = dataValueSet[supernodeSetIndex]; + // work out the new supernode ID + indexType newSupernodeID = nSupernodesAlready + supernode; - // the sort order will be dealt with later - // since all of these nodes are supernodes, they will be their own superparent, which means that: - // a. the regular2node can be set immediately - augmentedTree->regular2supernode [newRegularID] = newSupernodeID; - // b. as can the superparent - augmentedTree->superparents [newRegularID] = newSupernodeID; - } // per supernode in the set + // and the old supernode ID + // NB: May be NO_SUCH_ELEMENT if not already in the base tree + indexType oldSupernodeID = supernodeIDSet[supernodeSetIndex]; + + // At all levels above 0, we used to keep regular vertices in case they are attachment points. After augmentation, we don't need to. + // Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes + // In order to avoid confusion, we will copy the ID into a separate variable + indexType newRegularID = newSupernodeID; + + // setting the supernode's regular ID is now trivial + augmentedTree->supernodes [newSupernodeID] = newRegularID; + + // retrieve the old superID of the superparent. This is slightly tricky, as we have four classes of supernodes: + // 1. the root of the entire tree + // 2. attachment points not being inserted. In this case, the supernode ID is stored in the superparentSet array, not the superparent for insertion purposes + // 3. attachment points being inserted. In this case, the superparent is stored in the superparentSet array + // 4. "ordinary" supernodes, where the superparent is the same as the supernode ID anyway + // + // Note that an attachment point gets inserted into a parent superarc. But the attachment point itself has a NULL superarc, because it's only a virtual insertion + // This means that such an attachment superarc cannot be the superparent of any other attachment point + // It is therefore reasonable to deal with 1. & 2 separately. 3. & 4. then combine together + + // first we test for the root of the tree + if ((roundNo == baseTree->nRounds) && (supernode == supernodeSorter.size() - 1)) + { // root of the tree + // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the root is in every tree + // set the super arrays + augmentedTree->superarcs [newSupernodeID] = NO_SUCH_ELEMENT; + // hyperstructure carries over, so we use the same hyperparent as before + augmentedTree->hyperparents [newSupernodeID] = baseTree->hyperparents[oldSupernodeID]; + // and set the hypernode ID + augmentedTree->super2hypernode [newSupernodeID] = baseTree->super2hypernode[oldSupernodeID]; + // and the round and iteration + augmentedTree->whichRound [newSupernodeID] = baseTree->whichRound[oldSupernodeID]; + augmentedTree->whichIteration [newSupernodeID] = baseTree->whichIteration[oldSupernodeID]; + // and set the relevant regular arrays + augmentedTree->regularNodeGlobalIDs [newRegularID] = globalRegularIDSet[supernodeSetIndex]; + augmentedTree->dataValues [newRegularID] = dataValueSet[supernodeSetIndex]; + // for the root, these always point to itself + augmentedTree->regular2supernode [newRegularID] = newSupernodeID; + augmentedTree->superparents [newRegularID] = newSupernodeID; + } // root of the tree + // now deal with unsimplified attachment points, which we can identify because they were in the "kept" batch, not the "inserted" batch, + // and this is given away by the index into the set of supernodes to be added + // and the fact that the superarc is NO_SUCH_ELEMENT + else if ((supernodeSetIndex >= nInsertedSupernodes) && (noSuchElement(baseTree->superarcs[oldSupernodeID]))) + { // preserved attachment point + // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the supernode came from this block originally + // set the superarc to NO_SUCH_ELEMENT, as before + augmentedTree->superarcs [newSupernodeID] = NO_SUCH_ELEMENT; + // hyperstructure carries over, so we use the same hyperparent as before + augmentedTree->hyperparents [newSupernodeID] = baseTree->hyperparents[oldSupernodeID]; + // attachment points are never hypernodes anyway, so set it directly + augmentedTree->super2hypernode [newSupernodeID] = NO_SUCH_ELEMENT; + // and the round and iteration + augmentedTree->whichRound [newSupernodeID] = baseTree->whichRound[oldSupernodeID]; + augmentedTree->whichIteration [newSupernodeID] = baseTree->whichIteration[oldSupernodeID]; + // and set the relevant regular arrays + augmentedTree->regularNodeGlobalIDs [newRegularID] = globalRegularIDSet[supernodeSetIndex]; + augmentedTree->dataValues [newRegularID] = dataValueSet[supernodeSetIndex]; + // for a preserved attachment point, this always points to itself + augmentedTree->regular2supernode [newRegularID] = newSupernodeID; + // the superparent is the tricky one, as the old one may have been broken up by insertions at a higher level + + // Here, what we need to do is a search in the augmented tree to find which superarc to attach to. This is necessary + // because the old superarc it attached to may have been broken up. + // We are guaranteed that there is one, and that it only uses the higher levels of the augmented tree, + // so the fact that we are partially constructed doesn't get in the way. To do this, we need supernodes + // known to be in the higher level that are above and below the supernode. + // Since the point was an attachment point in the base tree, that means that there is a higher round superarc + // it inserts into. Moreover, the algorithm ALWAYS inserts a supernode at or above its original round, so + // we can guarantee that both ends of the parent are in the higher levels. Which means we only need to work + // out which end is higher. + indexType oldRegularID = baseTree->supernodes[oldSupernodeID]; + indexType oldSuperFrom = baseTree->superparents[oldRegularID]; + indexType oldSuperTo = baseTree->superarcs[oldSuperFrom]; + // retrieve the ascending flag + bool ascendingSuperarc = isAscending(oldSuperTo); + // and mask out the flags + oldSuperTo = maskedIndex(oldSuperTo); + + // since we haven't set up the regular search array yet, we can't use that + // instead, we know that the two supernodes must be in the new tree, so we retrieve their new super IDs + // and convert them to regular + + // retrieve their new super IDs + indexType newSuperFrom = newSupernodeIDs[oldSuperFrom]; + indexType newSuperTo = newSupernodeIDs[oldSuperTo]; + + // convert to regular IDs (which is what the FindSuperArcForUnknownNode() routine assumes) + indexType newRegularFrom = augmentedTree->supernodes[newSuperFrom]; + indexType newRegularTo = augmentedTree->supernodes[newSuperTo]; + + // the new superparent after the search + indexType newSuperparentID = NO_SUCH_ELEMENT; + + // depending on the ascending flag + if (ascendingSuperarc) + newSuperparentID = augmentedTree->FindSuperArcForUnknownNode(globalRegularIDSet[supernodeSetIndex],dataValueSet[supernodeSetIndex], newRegularTo, newRegularFrom); + else + newSuperparentID = augmentedTree->FindSuperArcForUnknownNode(globalRegularIDSet[supernodeSetIndex],dataValueSet[supernodeSetIndex], newRegularFrom, newRegularTo); + + // attachment points use the superparent to store the superarc they insert onto + augmentedTree->superparents [newRegularID] = newSuperparentID; + + } // preserved attachment point + else + { // raised attachment point or "ordinary" supernodes + // Since all of the superparents must be in the base tree, we can now retrieve the target + indexType superparentOldSuperID = maskedIndex(superparentSet[supernodeSetIndex]); + indexType oldTargetSuperID = baseTree->superarcs[superparentOldSuperID]; + // and break it into a target and flags + bool ascendingSuperarc = isAscending(oldTargetSuperID); + // NOTE: if the target was NO_SUCH_ELEMENT, this will hold 0 + oldTargetSuperID = maskedIndex(oldTargetSuperID); + + // and another boolean for whether we are the last element in a segment + bool isLastInSegment = false; + + // end of the entire array counts as last in segment + if (supernode == supernodeSorter.size() - 1) + isLastInSegment = true; + // otherwise, check for a mismatch in the sorting superparent which indicates the end of a segment + else if (maskedIndex(superparentSet[supernodeSetIndex]) != maskedIndex(superparentSet[supernodeSorter[supernode+1]])) + isLastInSegment = true; + + // setting the superarc is done the usual way. Our sort routine has ended up with the supernodes arranged in either ascending or descending order + // inwards along the parent superarc (as expressed by the superparent ID). Each superarc except the last in the segment points to the next one: + // the last one points to the target of the original superarc. + if (isLastInSegment) + { // last in segment + // we take the old target of the superarc (in old supernode IDs) and convert it to a new supernode ID + augmentedTree->superarcs[newSupernodeID] = newSupernodeIDs[oldTargetSuperID] | (ascendingSuperarc ? IS_ASCENDING : 0x00); + } // last in segment + else + { // not last in segment + // the target is always the next one, so just store it with the ascending flag + augmentedTree->superarcs[newSupernodeID] = (newSupernodeID+1) | (ascendingSuperarc ? IS_ASCENDING : 0x00); + } // not last in segment + + // first we identify the hyperarc on which the superarc sits + // this will be visible in the old base tree, since hyperstructure carries over + indexType oldHyperparent = baseTree->hyperparents[superparentOldSuperID]; + + // hyperstructure carries over, so we use the same hyperparent as the superparent + augmentedTree->hyperparents [newSupernodeID] = oldHyperparent; + // retrieve the hyperparent's old supernode ID & convert to a new one, then test it + if (newSupernodeIDs[baseTree->hypernodes[oldHyperparent]] == newSupernodeID) + augmentedTree->super2hypernode [newSupernodeID] = oldHyperparent; + else + augmentedTree->super2hypernode [newSupernodeID] = NO_SUCH_ELEMENT; + + // round and iteration are set from the superparent, since we are raising to its level + augmentedTree->whichRound [newSupernodeID] = baseTree->whichRound[superparentOldSuperID]; + augmentedTree->whichIteration [newSupernodeID] = baseTree->whichIteration[superparentOldSuperID]; + // and set the relevant regular arrays + augmentedTree->regularNodeGlobalIDs [newRegularID] = globalRegularIDSet[supernodeSetIndex]; + augmentedTree->dataValues [newRegularID] = dataValueSet[supernodeSetIndex]; + // for all supernodes, this points to itself + augmentedTree->regular2supernode [newRegularID] = newSupernodeID; + // and since we're inserted, so does this + augmentedTree->superparents [newRegularID] = newSupernodeID; + } // raised attachment point or "ordinary" supernodes + + } // per supernode in the set */ + + } // operator()() private: const vtkm::Id NumSupernodesAlready; const vtkm::Id BaseTreeNumRounds; - const vtkm::Id AugmentedTreeNumIterations; - const vtkm::Id RoundNumber; - const vtkm::Id NumAugmentedTreeSupernodes; + const vtkm::Id NumInsertedSupernodes; + const vtkm::Id RoundNo; }; // CreateSuperarcsWorklet diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h index f4c745478..27af1e05a 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h @@ -81,9 +81,17 @@ public: using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9); using InputDomain = _1; + /// Default Constructor VTKM_EXEC_CONT - FindSuperparentForNecessaryNodesWorklet() {} + FindSuperparentForNecessaryNodesWorklet(vtkm::Id3 meshBlockOrigin, + vtkm::Id3 meshBlockSize, + vtkm::Id3 meshGlobalSize) + : MeshBlockOrigin(meshBlockOrigin) + , MeshBlockSize(meshBlockSize) + , MeshGlobalSize(meshGlobalSize) + { + } /// operator() of the workelt template IsInMesh(globalRegularId)) + { + // Set to NO_SUCH_ELEMENT by default. By doing this in the worklet we an avoid having to + // initialize the output arrays first and we can use FieldIn instead of FieldInOut + regularSuperparentsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + regularNodesNeededValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + } // if it fails this test, then it's already in tree - if (vtkm::worklet::contourtree_augmented::NoSuchElement(newRegularId)) + else if (vtkm::worklet::contourtree_augmented::NoSuchElement(newRegularId)) { // not yet in tree // since it's not in the tree, we want to find where it belongs // to do so, we need to find an "above" and "below" node for it. Since it exists in the old tree, it belongs to a superarc, and we @@ -200,7 +216,80 @@ public: */ } // operator()() -}; // FindSuperparentForNecessaryNodesWorklet + +private: + // Mesh data + vtkm::Id3 MeshBlockOrigin; + vtkm::Id3 MeshBlockSize; + vtkm::Id3 MeshGlobalSize; + + VTKM_EXEC + bool IsInMesh(vtkm::Id globalId) const + { // IsInMesh() + if (this->MeshGlobalSize[2] > 1) // 3D + { + // convert from global ID to global coords + vtkm::Id globalSliceSize = this->MeshGlobalSize[0] * this->MeshGlobalSize[1]; + vtkm::Id globalSlice = globalId / globalSliceSize; + vtkm::Id globalRow = globalId / this->MeshGlobalSize[0]; + vtkm::Id globalCol = globalId % this->MeshGlobalSize[0]; + + // test validity + if (globalSlice < this->MeshBlockOrigin[2]) + { + return false; + } + if (globalSlice >= this->MeshBlockOrigin[2] + this->MeshBlockSize[2]) + { + return false; + } + if (globalRow < this->MeshBlockOrigin[1]) + { + return false; + } + if (globalRow >= this->MeshBlockOrigin[1] + this->MeshBlockSize[1]) + { + return false; + } + if (globalCol < this->MeshBlockOrigin[0]) + { + return false; + } + if (globalCol >= this->MeshBlockOrigin[0] + this->MeshBlockSize[0]) + { + return false; + } + // it's in the block - return true + return true; + } // end if 3D + else // 2D mesh + { + // convert from global ID to global coords + vtkm::Id globalRow = globalId / this->MeshGlobalSize[0]; + vtkm::Id globalCol = globalId % this->MeshGlobalSize[0]; + + // test validity + if (globalRow < this->MeshBlockOrigin[1]) + { + return false; + } + if (globalRow >= this->MeshBlockOrigin[1] + this->MeshBlockSize[1]) + { + return false; + } + if (globalCol < this->MeshBlockOrigin[0]) + { + return false; + } + if (globalCol >= this->MeshBlockOrigin[0] + this->MeshBlockSize[0]) + { + return false; + } + // it's in the block - return true + return true; + } + } // IsInMesh() +}; // FindSuperparentForNecessaryNodesWorklet } // namespace hierarchical_augmenter } // namespace contourtree_distributed diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/IsAttachementPointPredicate.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/IsAttachementPointPredicate.h index dd79e4ce7..a85ccfd94 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/IsAttachementPointPredicate.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/IsAttachementPointPredicate.h @@ -81,27 +81,50 @@ public: const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs, const vtkm::worklet::contourtree_augmented::IdArrayType& whichRound, const vtkm::Id numRounds, + vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray, + vtkm::Id presimplifyThreshold, vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) : SuperarcsPortal(superarcs.PrepareForInput(device, token)) , WhichRoundPortal(whichRound.PrepareForInput(device, token)) , NumRounds(numRounds) + , PresimplifyThreshold(presimplifyThreshold) { // constructor + this->Presimplify = ((volumeArray != NULL) && (presimplifyThreshold > 0)); + // If we presimplify then store the volumeArray. Otherwise we don't need to volume array and we + // set it to another portal, just to make sure the variable is being initalized with something + this->VolumeArrayPortal = + this->Presimplify ? volumeArray->PrepareForInput(device, token) : this->WhichRoundPortal; } // constructor // () operator - gets called to do comparison VTKM_EXEC bool operator()(const vtkm::Id& supernode) const { // operator() - return ( - vtkm::worklet::contourtree_augmented::NoSuchElement(this->SuperarcsPortal.Get(supernode)) && - (this->WhichRoundPortal.Get(supernode) < this->NumRounds)); + // an attachment point is defined by having no superarc (NO_SUCH_ELEMENT) and not being in + // the final round (where this indicates the global root) + bool predicate = + (vtkm::worklet::contourtree_augmented::NoSuchElement(this->SuperarcsPortal.Get(supernode)) && + (this->WhichRoundPortal.Get(supernode) < this->NumRounds)); + // if we pass this check then we need to also check that the supernode passes the pre-simplification threshold + if (predicate && this->Presimplify) + { + // suppress if it's volume is at or below the threshold + if (this->VolumeArrayPortal.Get(supernode) <= this->PresimplifyThreshold) + { // below threshold + predicate = false; // do not keep attachement point below the simplification threshold + } // below threshold + } + return predicate; } // operator() private: IdPortalType SuperarcsPortal; IdPortalType WhichRoundPortal; const vtkm::Id NumRounds; + bool Presimplify; + IdPortalType VolumeArrayPortal; + vtkm::Id PresimplifyThreshold; }; // IsAttachementPointPredicateImpl @@ -113,24 +136,35 @@ public: VTKM_CONT IsAttachementPointPredicate(const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs, const vtkm::worklet::contourtree_augmented::IdArrayType& whichRound, - const vtkm::Id numRounds) + const vtkm::Id numRounds, + vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray = NULL, + vtkm::Id presimplifyThreshold = 0) : Superarcs(superarcs) , WhichRound(whichRound) , NumRounds(numRounds) + , VolumeArray(volumeArray) + , PresimplifyThreshold(presimplifyThreshold) { } VTKM_CONT IsAttachementPointPredicateImpl PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const { - return IsAttachementPointPredicateImpl( - this->Superarcs, this->WhichRound, this->NumRounds, device, token); + return IsAttachementPointPredicateImpl(this->Superarcs, + this->WhichRound, + this->NumRounds, + this->VolumeArray, + this->PresimplifyThreshold, + device, + token); } private: vtkm::worklet::contourtree_augmented::IdArrayType Superarcs; vtkm::worklet::contourtree_augmented::IdArrayType WhichRound; const vtkm::Id NumRounds; + vtkm::worklet::contourtree_augmented::IdArrayType* VolumeArray; + vtkm::Id PresimplifyThreshold; }; // IsAttachementPointPredicate } // namespace hierarchical_augmenter diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/ResizeArraysBuildNewSupernodeIdsWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/ResizeArraysBuildNewSupernodeIdsWorklet.h index 0903f85f4..927854844 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/ResizeArraysBuildNewSupernodeIdsWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/ResizeArraysBuildNewSupernodeIdsWorklet.h @@ -65,12 +65,16 @@ public: /// Control signature for the worklet using ControlSignature = void( FieldIn supernodeIndex, // input domain ArrayHandleIndex(SupernodeSorter.GetNumberOfValues()) + // FieldIn supernodeIdSetPermuted, // input supernodeIDSet permuted by supernodeSorter to allow for FieldIn FieldIn - supernodeIdSetPermuted, // input supernodeIDSet permuted by supernodeSorter to allow for FieldIn + globalRegularIdSet, // input globalRegularIdSet permuted by supernodeSorter to allow for FieldIn + ExecObject findRegularByGlobal, + WholeArrayIn baseTreeRegular2Supernode, WholeArrayInOut newSupernodeIds // output/input (both are necessary since not all valyes will be overwritten) ); - using ExecutionSignature = void(_1, _2, _3); + + using ExecutionSignature = void(_1, _2, _3, _4, _5); using InputDomain = _1; // Default Constructor @@ -80,10 +84,13 @@ public: { } - template + template VTKM_EXEC void operator()( - const vtkm::Id& supernode, // InputIndex of supernodeSorter - const vtkm::Id& oldSupernodeId, // same as supernodeIDSet[supernodeSetIndex]; + const vtkm::Id& supernode, // InputIndex of supernodeSorter + // const vtkm::Id& oldSupernodeId, // same as supernodeIDSet[supernodeSetIndex]; + const vtkm::Id& globalRegularIdSetValue, // same as globalRegularIDSet[supernodeSetIndex] + const ExecObjectType& findRegularByGlobal, + const InFieldPortalType& baseTreeRegular2SupernodePortal, const InOutFieldPortalType& newSupernodeIdsPortal) const { // per supernode @@ -96,6 +103,17 @@ public: // that if it came from another block it will be set to NO_SUCH_ELEMENT // vtkm::Id oldSupernodeId set on input since we use ArrayHandlePermutation to // shuffle supernodeIDSet by supernodeSorter; + // TODO/WARNING: Logic error in that comment for presimplified trees, but not for the original version. See RetrieveOldSupernodes() for why. + + // TODO/WARNING: We substitute a search in the old hierarchical tree for the supernode. If it is present, then we fill in it's entry in the + + // newSupernodeIDs array. If not, we're happy. + vtkm::Id oldRegularId = findRegularByGlobal.FindRegularByGlobal(globalRegularIdSetValue); + vtkm::Id oldSupernodeId = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT; + if (!vtkm::worklet::contourtree_augmented::NoSuchElement(oldRegularId)) + { + oldSupernodeId = baseTreeRegular2SupernodePortal.Get(oldRegularId); + } // and write to the lookup array if (!vtkm::worklet::contourtree_augmented::NoSuchElement(oldSupernodeId)) diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/UpdateHyperstructureSetSuperchildrenWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/UpdateHyperstructureSetSuperchildrenWorklet.h index c81154f9e..f0b67fb7e 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/UpdateHyperstructureSetSuperchildrenWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/UpdateHyperstructureSetSuperchildrenWorklet.h @@ -66,10 +66,12 @@ class UpdateHyperstructureSetSuperchildrenWorklet : public vtkm::worklet::Workle public: /// Control signature for the worklet using ControlSignature = void( - WholeArrayIn augmentedTreeHypernodes, // input (we need both this and the next value) - FieldOut augmentedTreeSuperchildren // output + WholeArrayIn augmentedTreeHypernodes, // input (we need both this and the next value) + FieldIn augmentedTreeSuperarcs, // input + WholeArrayIn augmentedTreeHyperparents, // input + WholeArrayInOut augmentedTreeSuperchildren // output ); - using ExecutionSignature = void(InputIndex, _1, _2); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4); using InputDomain = _1; // Default Constructor @@ -80,46 +82,58 @@ public: } - template - VTKM_EXEC void operator()( - const vtkm::Id& hypernode, - const InFieldPortalType& augmentedTreeHypernodesPortal, - vtkm::Id& - augmentedTreeSuperchildrenValue // same as augmentedTree->superchildren[InputIndex] = ... - ) const + template + VTKM_EXEC void operator()(const vtkm::Id& supernode, + const InFieldPortalType1& augmentedTreeHypernodesPortal, + const vtkm::Id& augmentedTreeSuperarcsValue, + const InFieldPortalType2& augmentedTreeHyperparentsPortal, + const OutFieldPortalType& augmentedTreeSuperchildrenPortal) const { - // per hypernode - // retrieve the new superId - vtkm::Id superId = augmentedTreeHypernodesPortal.Get(hypernode); - // and the next one over - vtkm::Id nextSuperId; - if (hypernode == augmentedTreeHypernodesPortal.GetNumberOfValues() - 1) + // per supernode + // attachment points have NULL superarcs and are skipped + if (vtkm::worklet::contourtree_augmented::NoSuchElement(augmentedTreeSuperarcsValue)) { - nextSuperId = this->AugmentedTreeNumSupernodes; + return; } - else - { - nextSuperId = augmentedTreeHypernodesPortal.Get(hypernode + 1); - } - // the difference is the number of superchildren - augmentedTreeSuperchildrenValue = nextSuperId - superId; + // we are now guaranteed to have a valid hyperparent + vtkm::Id hyperparent = augmentedTreeHyperparentsPortal.Get(supernode); + vtkm::Id hyperparentSuperId = augmentedTreeHypernodesPortal.Get(hyperparent); + // we could be at the end of the array, so test explicitly + if (supernode == this->AugmentedTreeNumSupernodes - 1) + { + // this means that we are the end of the segment and can subtract the hyperparent's super ID to get the number of superchildren + augmentedTreeSuperchildrenPortal.Set(hyperparent, + this->AugmentedTreeNumSupernodes - hyperparentSuperId); + } + // otherwise, if our hyperparent is different from our neighbor's, we are the end of the segment + else if (hyperparent != augmentedTreeHyperparentsPortal.Get(supernode + 1)) + { + // again, subtract to get the number + augmentedTreeSuperchildrenPortal.Set( + hyperparent, this->AugmentedTreeNumSupernodes + 1 - hyperparentSuperId); + } + // per supernode // In serial this worklet implements the following operation /* - for (vtkm::Id hypernode = 0; hypernode < augmentedTree->hypernodes.size(); hypernode++) - { // per hypernode - // retrieve the new super ID - vtkm::Id superID = augmentedTree->hypernodes[hypernode]; - // and the next one over - vtkm::Id nextSuperID; - if (hypernode == augmentedTree->hypernodes.size() - 1) - nextSuperID = augmentedTree->supernodes.size(); - else - nextSuperID = augmentedTree->hypernodes[hypernode+1]; - // the difference is the number of superchildren - augmentedTree->superchildren[hypernode] = nextSuperID - superID; - } // per hypernode + for (indexType supernode = augmentedTree->firstSupernodePerIteration[roundNo][0]; supernode < augmentedTree->firstSupernodePerIteration[roundNo][augmentedTree->nIterations[roundNo]]; supernode++) + { // per supernode + // attachment points have NULL superarcs and are skipped + if (noSuchElement(augmentedTree->superarcs[supernode])) + continue; + // we are now guaranteed to have a valid hyperparent + indexType hyperparent = augmentedTree->hyperparents[supernode]; + indexType hyperparentSuperID = augmentedTree->hypernodes[hyperparent]; + // we could be at the end of the array, so test explicitly + if (supernode == augmentedTree->supernodes.size() - 1) + // this means that we are the end of the segment and can subtract the hyperparent's super ID to get the number of superchildren + augmentedTree->superchildren[hyperparent] = augmentedTree->supernodes.size() - hyperparentSuperID; + // otherwise, if our hyperparent is different from our neighbor's, we are the end of the segment + else if (hyperparent != augmentedTree->hyperparents[supernode + 1]) + // again, subtract to get the number + augmentedTree->superchildren[hyperparent] = supernode + 1 - hyperparentSuperID; + } // per supernode */ } // operator()() diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/CMakeLists.txt b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/CMakeLists.txt index eb302fd40..18b089ddb 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/CMakeLists.txt @@ -17,6 +17,8 @@ set(headers TransferTargetComperator.h TransferWeightsUpdateRHEWorklet.h TransferWeightsUpdateLHEWorklet.h + TransferWeightsUpdateRHEWorkletRound2.h + TransferWeightsUpdateLHEWorkletRound2.h ) vtkm_declare_headers(${headers}) diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h index b885cb116..93882ac2b 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h @@ -116,7 +116,8 @@ public: else { // attachment point // set the transfer target - transferTarget = hierarchicalTreeSuperparentsPortal.Get(supernodeRegularId); + transferTarget = hierarchicalTreeSuperparentsPortal.Get(supernodeRegularId) | + vtkm::worklet::contourtree_augmented::TRANSFER_TO_SUPERARC; } // attachment point } // null superarc else diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h index 0d5c1ca98..68cc08b77 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h @@ -74,10 +74,11 @@ public: FieldIn sortedTransferTargetShiftedView, FieldIn valuePrefixSumShiftedView, WholeArrayInOut dependentValuesPortal); - using ExecutionSgnature = void(_1, _2, _3, _4); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4); template - VTKM_EXEC void operator()(const vtkm::Id& sortedTransferTargetValue, + VTKM_EXEC void operator()(const vtkm::Id& supernode, + const vtkm::Id& sortedTransferTargetValue, const vtkm::Id& sortedTransferTargetPreviousValue, const vtkm::Id& valuePrefixSumPreviousValue, InOutPortalType& dependentValuesPortal) const @@ -88,30 +89,33 @@ public: { return; } - if (sortedTransferTargetValue != sortedTransferTargetPreviousValue) + + // we need to separate out the flag for attachment points + bool superarcTransfer = + vtkm::worklet::contourtree_augmented::TransferToSuperarc(sortedTransferTargetValue); + vtkm::Id superarcOrNodeId = + vtkm::worklet::contourtree_augmented::MaskedIndex(sortedTransferTargetValue); + + // ignore the transfers for attachment points + if (superarcTransfer) { - auto originalValue = dependentValuesPortal.Get(sortedTransferTargetValue); + return; + } + + // the LHE at 0 is special - it subtracts zero. In practice, since NO_SUCH_ELEMENT will sort low, this will never + // occur, but let's keep the logic strict + auto originalValue = dependentValuesPortal.Get(superarcOrNodeId); + if (supernode == + 0) // i.e., supernode == firstSupernode but we view the arrays when we call the worklet + { // LHE 0 + dependentValuesPortal.Set(superarcOrNodeId, originalValue - 0); + } + else if (sortedTransferTargetValue != sortedTransferTargetPreviousValue) + { // LHE not 0 + dependentValuesPortal.Set(sortedTransferTargetValue, originalValue - valuePrefixSumPreviousValue); } - - - // In serial this worklet implements the following operation - /* - for (vtkm::Id supernode = firstSupernode + 1; supernode < lastSupernode; supernode++) - { // per supernode - // ignore any that point at NO_SUCH_ELEMENT - if (noSuchElement(sortedTransferTarget[supernode])) - continue; - - // the LHE at 0 is special - it subtracts zero. In practice, since NO_SUCH_ELEMENT will sort low, this will never - // occur, but let's keep the logic strict - if (sortedTransferTarget[supernode] != sortedTransferTarget[supernode-1]) - { // LHE not 0 - dependentValues[sortedTransferTarget[supernode]] -= valuePrefixSum[supernode-1]; - } // LHE not 0 - } // per supernode - */ } // operator()() }; // TransferWeightsUpdateLHEWorklet diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorkletRound2.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorkletRound2.h new file mode 100644 index 000000000..485572eef --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorkletRound2.h @@ -0,0 +1,132 @@ +//============================================================================ +// 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_hierarchical_hyper_sweeper_transfer_weights_update_lhe_worklet_round2_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_lhe_worklet_round2_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_hyper_sweeper +{ + +/// Worklet used in HierarchicalHyperSweeper.TransferWeights(...) to implement +/// step 7b. Now find the LHE of each group and subtract out the prior weight +class TransferWeightsUpdateLHEWorkletRound2 : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn sortedTransferTargetPortal, + FieldIn sortedTransferTargetShiftedView, + FieldIn valuePrefixSumShiftedView, + WholeArrayInOut intrinsicValues, + WholeArrayInOut dependentValues); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5); + + template + VTKM_EXEC void operator()(const vtkm::Id& supernode, + const vtkm::Id& sortedTransferTargetValue, + const vtkm::Id& sortedTransferTargetPreviousValue, + const vtkm::Id& valuePrefixSumPreviousValue, + InOutPortalType& intrinsicValuesPortal, + InOutPortalType& dependentValuesPortal) const + { + // per supernode + // ignore any that point at NO_SUCH_ELEMENT + if (vtkm::worklet::contourtree_augmented::NoSuchElement(sortedTransferTargetValue)) + { + return; + } + + // we need to separate out the flag for attachment points + bool superarcTransfer = + vtkm::worklet::contourtree_augmented::TransferToSuperarc(sortedTransferTargetValue); + vtkm::Id superarcOrNodeId = + vtkm::worklet::contourtree_augmented::MaskedIndex(sortedTransferTargetValue); + + // ignore the transfers for attachment points + if (superarcTransfer) + { + return; + } + + // the LHE at 0 is special - it subtracts zero. In practice, since NO_SUCH_ELEMENT will sort low, this will never + // occur, but let's keep the logic strict + auto originalIntrinsicValue = dependentValuesPortal.Get(superarcOrNodeId); + auto originalDependentValue = dependentValuesPortal.Get(superarcOrNodeId); + if (supernode == + 0) // i.e., supernode == firstSupernode but we view the arrays when we call the worklet + { // LHE 0 + intrinsicValuesPortal.Set(superarcOrNodeId, originalIntrinsicValue - 0); + dependentValuesPortal.Set(superarcOrNodeId, originalDependentValue - 0); + } + else if (sortedTransferTargetValue != sortedTransferTargetPreviousValue) + { // LHE not 0 + intrinsicValuesPortal.Set(sortedTransferTargetValue, + originalIntrinsicValue - valuePrefixSumPreviousValue); + dependentValuesPortal.Set(sortedTransferTargetValue, + originalDependentValue - valuePrefixSumPreviousValue); + } + } // operator()() +}; // TransferWeightsUpdateLHEWorklet + +} // namespace hierarchical_hyper_sweeper +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h index f1265b88a..108a019e0 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h @@ -100,8 +100,18 @@ public: if ((supernode == this->LastSupernode - 1) || (transferTarget != sortedTransferTargetPortal.Get(supernode + 1))) { // RHE of segment - auto originalValue = dependentValuesPortal.Get(transferTarget); - dependentValuesPortal.Set(transferTarget, originalValue + valuePrefixSum); + // we need to separate out the flag for attachment points + bool superarcTransfer = + vtkm::worklet::contourtree_augmented::TransferToSuperarc(transferTarget); + vtkm::Id superarcOrNodeId = + vtkm::worklet::contourtree_augmented::MaskedIndex(transferTarget); + // we ignore attachment points + if (superarcTransfer) + { + return; + } + auto originalValue = dependentValuesPortal.Get(superarcOrNodeId); + dependentValuesPortal.Set(superarcOrNodeId, originalValue + valuePrefixSum); } // RHE of segment } @@ -113,10 +123,19 @@ public: if (noSuchElement(sortedTransferTarget[supernode])) continue; - // the RHE of each segment transfers its weight (including all irrelevant prefixes) if ((supernode == lastSupernode - 1) || (sortedTransferTarget[supernode] != sortedTransferTarget[supernode+1])) { // RHE of segment - dependentValues[sortedTransferTarget[supernode]] += valuePrefixSum[supernode]; + // TODO / WARNING 11/07/2023 + // we need to separate out the flag for attachment points + bool superarcTransfer = transferToSuperarc(sortedTransferTarget[supernode]); + indexType superarcOrNodeID = maskedIndex(sortedTransferTarget[supernode]); + + // we ignore attachment points + if (superarcTransfer) + continue; + + // transfer as dependent weight + dependentValues[superarcOrNodeID] += valuePrefixSum[supernode]; } // RHE of segment } // per supernode */ diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorkletRound2.h b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorkletRound2.h new file mode 100644 index 000000000..d4e064e40 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorkletRound2.h @@ -0,0 +1,161 @@ +//============================================================================ +// 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_hierarchical_hyper_sweeper_transfer_weights_update_rhe_worklet_round2_h +#define vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_rhe_worklet_round2_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace contourtree_distributed +{ +namespace hierarchical_hyper_sweeper +{ + +/// Worklet used in HierarchicalHyperSweeper.TransferWeights(...) to implement +/// step 7a. Find the RHE of each group and transfer the prefix sum weight. +/// Note that we do not compute the transfer weight separately, we add it in place instead +class TransferWeightsUpdateRHEWorkletRound2 : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = + void(FieldIn supernodeIndex, // input counting array [firstSupernode, lastSupernode) + WholeArrayIn sortedTransferTarget, + FieldIn valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode) + WholeArrayInOut intrinsicValuesPortal, + WholeArrayInOut dependentValuesPortal); + using ExecutionSignature = void(_1, _2, _3, _4, _5); + + // Default Constructor + VTKM_EXEC_CONT + TransferWeightsUpdateRHEWorkletRound2(const vtkm::Id& lastSupernode) + : LastSupernode(lastSupernode) + { + } + + template + VTKM_EXEC void operator()(const vtkm::Id& supernode, + const InPortalType& sortedTransferTargetPortal, + const vtkm::Id& valuePrefixSum, // same as valuePrefixSum[supernode] + OutPortalType& intrinsicValuesPortal, + OutPortalType& dependentValuesPortal) const + { + // per supernode + // ignore any that point at NO_SUCH_ELEMENT + vtkm::Id transferTarget = sortedTransferTargetPortal.Get(supernode); + if (!vtkm::worklet::contourtree_augmented::NoSuchElement(transferTarget)) + { + // the RHE of each segment transfers its weight (including all irrelevant prefixes) + if ((supernode == this->LastSupernode - 1) || + (transferTarget != sortedTransferTargetPortal.Get(supernode + 1))) + { // RHE of segment + // we need to separate out the flag for attachment points + bool superarcTransfer = + vtkm::worklet::contourtree_augmented::TransferToSuperarc(transferTarget); + vtkm::Id superarcOrNodeId = + vtkm::worklet::contourtree_augmented::MaskedIndex(transferTarget); + // we ignore attachment points + if (superarcTransfer) + { + return; + } + // we modify both intrinsic and dependent values + // we modify both intrinsic and dependent values + auto originalIntrinsicValue = intrinsicValuesPortal.Get(superarcOrNodeId); + intrinsicValuesPortal.Set(superarcOrNodeId, originalIntrinsicValue + valuePrefixSum); + auto originalDependentValue = dependentValuesPortal.Get(superarcOrNodeId); + dependentValuesPortal.Set(superarcOrNodeId, originalDependentValue + valuePrefixSum); + } // RHE of segment + } + + // In serial this worklet implements the following operation + /* + for (indexType supernode = firstSupernode; supernode < lastSupernode; supernode++) + { // per supernode + // ignore any that point at NO_SUCH_ELEMENT + if (noSuchElement(sortedTransferTarget[supernode])) + continue; + + // the RHE of each segment transfers its weight (including all irrelevant prefixes) + if ((supernode == lastSupernode - 1) || (sortedTransferTarget[supernode] != sortedTransferTarget[supernode+1])) + { // RHE of segment + // we need to separate out the flag for attachment points + bool superarcTransfer = transferToSuperarc(sortedTransferTarget[supernode]); + indexType superarcOrNodeID = maskedIndex(sortedTransferTarget[supernode]); + // ignore the transfers for non-attachment points + if (!superarcTransfer) + continue; + + // we modify both intrinsic and dependent values + intrinsicValues[superarcOrNodeID] += valuePrefixSum[supernode]; + dependentValues[superarcOrNodeID] += valuePrefixSum[supernode]; + } // RHE of segment + + } // per supernode + */ + } // operator()() + +private: + const vtkm::Id LastSupernode; + +}; // TransferWeightsUpdateRHEWorklet + +} // namespace hierarchical_hyper_sweeper +} // namespace contourtree_distributed +} // namespace worklet +} // namespace vtkm + +#endif