Added serial (via DIY) unit test for distributed contour tree filter.

This commit is contained in:
Gunther H. Weber 2020-09-24 12:44:44 -07:00
parent ca36a6a420
commit 7a77b1ad89
6 changed files with 759 additions and 168 deletions

@ -59,7 +59,7 @@ find_package(VTKm REQUIRED QUIET)
# MPI
####################################
if (VTKm_ENABLE_MPI)
add_executable(ContourTree_Distributed ContourTreeApp.cxx TreeCompiler.cxx)
add_executable(ContourTree_Distributed ContourTreeApp.cxx)
target_link_libraries(ContourTree_Distributed vtkm_filter vtkm_io MPI::MPI_CXX)
vtkm_add_target_information(ContourTree_Distributed
MODIFY_CUDA_FLAGS

@ -75,6 +75,7 @@
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/ProcessContourTree.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_distributed/TreeCompiler.h>
#ifdef ENABLE_SET_NUM_THREADS
#include "tbb/task_scheduler_init.h"
@ -98,7 +99,6 @@ VTKM_THIRDPARTY_POST_INCLUDE
#include <utility>
#include <vector>
#include "TreeCompiler.h"
#define PRESPLIT_FILE
@ -868,7 +868,7 @@ int main(int argc, char* argv[])
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
TreeCompiler treeCompiler;
vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
char fname[256];
std::snprintf(fname,

@ -1,138 +0,0 @@
#ifndef _TREECOMPILER_H_
#define _TREECOMPILER_H_
#include <iostream>
#include <vtkm/Types.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
// FIXME/HACK: Define here for compatibility with PPP TreeCompiler
typedef double dataType;
typedef unsigned long long indexType;
// small class for storing the contour arcs
class Edge
{ // Edge
public:
indexType low, high;
// constructor - defaults to -1
Edge(vtkm::Id Low = -1, vtkm::Id High = -1)
: low(Low)
, high(High)
{
}
}; // Edge
// comparison operator <
inline bool operator<(const Edge LHS, const Edge RHS)
{ // operator <
if (LHS.low < RHS.low)
return true;
if (LHS.low > RHS.low)
return false;
if (LHS.high < RHS.high)
return true;
if (LHS.high > RHS.high)
return false;
return false;
} // operator <
// a helper class which stores a single supernode inserted onto a superarc
class SupernodeOnSuperarc
{ // class SupernodeOnSuperarc
public:
// the global ID of the supernode
indexType globalID;
// the data value stored at the supernode
dataType dataValue;
// the low and high ends of the superarc it is on (may be itself)
indexType lowEnd, highEnd;
// constructor
SupernodeOnSuperarc(indexType GlobalID = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
dataType DataValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType LowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType HighEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)
: globalID(GlobalID)
, dataValue(DataValue)
, lowEnd(LowEnd)
, highEnd(HighEnd)
{ // constructor
} // constructor
}; // class SupernodeOnSuperarc
// overloaded comparison operator
// primary sort is by superarc (low, high),
// then secondary sort on datavalue
// tertiary on globalID to implement simulated simplicity
inline bool operator<(const SupernodeOnSuperarc& left, const SupernodeOnSuperarc& right)
{ // < operator
// simple lexicographic sort
if (left.lowEnd < right.lowEnd)
return true;
if (left.lowEnd > right.lowEnd)
return false;
if (left.highEnd < right.highEnd)
return true;
if (left.highEnd > right.highEnd)
return false;
if (left.dataValue < right.dataValue)
return true;
if (left.dataValue > right.dataValue)
return false;
if (left.globalID < right.globalID)
return true;
if (left.globalID > right.globalID)
return false;
// fall-through (shouldn't happen, but)
// if they're the same, it's false
return false;
} // < operator
// stream output
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node);
// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node);
// the class that compiles the contour tree
class TreeCompiler
{ // class TreeCompiler
public:
// we want a vector of supernodes on superarcs
std::vector<SupernodeOnSuperarc> supernodes;
// and a vector of Edges (the output)
std::vector<Edge> superarcs;
// default constructor sets it to empty
TreeCompiler()
{ // constructor
// clear out the supernode array
supernodes.resize(0);
// and the superarc array
superarcs.resize(0);
} // constructor
// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
void AddHierarchicalTree(const vtkm::cont::DataSet& addedTree);
// routine to compute the actual superarcs
void ComputeSuperarcs();
// routine to print the superarcs
void PrintSuperarcs();
// routine to write out binary file
void WriteBinary(FILE* outFile);
// routine to read in binary file & append to contents
void ReadBinary(FILE* inFile);
}; // class TreeCompiler
#endif

@ -19,6 +19,7 @@ set(unit_tests
UnitTestContourFilterNormals.cxx
UnitTestContourTreeUniformFilter.cxx
UnitTestContourTreeUniformAugmentedFilter.cxx
UnitTestContourTreeUniformDistributedFilter.cxx
UnitTestCoordinateSystemTransform.cxx
UnitTestCrossProductFilter.cxx
UnitTestDotProductFilter.cxx

@ -0,0 +1,519 @@
//============================================================================
// 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)
//==============================================================================
#include <vtkm/filter/ContourTreeUniformDistributed.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_distributed/TreeCompiler.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/VTKDataSetReader.h>
namespace
{
// numberOf Blocks must be a power of 2
vtkm::Id3 ComputeNumberOfBlocksPerAxis(vtkm::Id3 globalSize, vtkm::Id numberOfBlocks)
{
std::cout << "GlobalSize: " << globalSize << " numberOfBlocks:" << numberOfBlocks << " -> ";
// Inefficient way to compute log2 of numberOfBlocks, i.e., number of total splits
int numSplits = 0;
bool isPowerOfTwo = true;
while (numberOfBlocks > 1)
{
if (numberOfBlocks % 2 != 0)
{
isPowerOfTwo = false;
break;
}
numberOfBlocks /= 2;
++numSplits;
}
if (isPowerOfTwo)
{
vtkm::Id3 splitsPerAxis{ 0, 0, 0 };
while (numSplits > 0)
{
// Find split axis as axis with largest extent
vtkm::Id splitAxis = 0;
for (vtkm::Id d = 1; d < 3; ++d)
if (globalSize[d] > globalSize[splitAxis])
splitAxis = d;
// Split in half along that axis
std::cout << splitAxis << " " << globalSize << std::endl;
VTKM_ASSERT(globalSize[splitAxis] > 1);
++splitsPerAxis[splitAxis];
globalSize[splitAxis] /= 2;
--numSplits;
}
std::cout << "splitsPerAxis: " << splitsPerAxis;
vtkm::Id3 blocksPerAxis;
for (vtkm::Id d = 0; d < 3; ++d)
blocksPerAxis[d] = 1 << splitsPerAxis[d];
std::cout << " blocksPerAxis: " << blocksPerAxis << std::endl;
return blocksPerAxis;
}
else
{
std::cout << "numberOfBlocks is not a power of two. Splitting along longest axis" << std::endl;
vtkm::Id splitAxis = 0;
for (vtkm::Id d = 1; d < 3; ++d)
if (globalSize[d] > globalSize[splitAxis])
splitAxis = d;
vtkm::Id3 blocksPerAxis{ 1, 1, 1 };
blocksPerAxis[splitAxis] = numberOfBlocks;
return blocksPerAxis;
}
}
std::tuple<vtkm::Id3, vtkm::Id3, vtkm::Id3> ComputeBlockExtents(vtkm::Id3 globalSize,
vtkm::Id3 blocksPerAxis,
vtkm::Id blockNo)
{
std::cout << "Block " << blockNo;
vtkm::Id3 blockIndex, blockOrigin, blockSize;
for (vtkm::Id d = 0; d < 3; ++d)
{
blockIndex[d] = blockNo % blocksPerAxis[d];
blockNo /= blocksPerAxis[d];
float dx = float(globalSize[d] - 1) / float(blocksPerAxis[d]);
blockOrigin[d] = vtkm::Id(blockIndex[d] * dx);
vtkm::Id maxIdx =
blockIndex[d] < blocksPerAxis[d] - 1 ? vtkm::Id((blockIndex[d] + 1) * dx) : globalSize[d] - 1;
blockSize[d] = maxIdx - blockOrigin[d] + 1;
}
std::cout << " -> " << blockIndex << " " << blockOrigin << " " << blockSize << std::endl;
return std::make_tuple(blockIndex, blockOrigin, blockSize);
}
vtkm::cont::DataSet CreateSubDataSet(const vtkm::cont::DataSet& ds,
vtkm::Id3 blockOrigin,
vtkm::Id3 blockSize,
const std::string& fieldName)
{
vtkm::Id3 globalSize;
ds.GetCellSet().CastAndCall(vtkm::worklet::contourtree_augmented::GetPointDimensions(),
globalSize);
const vtkm::Id nOutValues = blockSize[0] * blockSize[1] * blockSize[2];
const auto inDataArrayHandle = ds.GetPointField(fieldName).GetData();
vtkm::cont::ArrayHandle<vtkm::Id> copyIdsArray;
copyIdsArray.Allocate(nOutValues);
auto copyIdsPortal = copyIdsArray.WritePortal();
vtkm::Id3 outArrIdx;
for (outArrIdx[2] = 0; outArrIdx[2] < blockSize[2]; ++outArrIdx[2])
for (outArrIdx[1] = 0; outArrIdx[1] < blockSize[1]; ++outArrIdx[1])
for (outArrIdx[0] = 0; outArrIdx[0] < blockSize[0]; ++outArrIdx[0])
{
vtkm::Id3 inArrIdx = outArrIdx + blockOrigin;
vtkm::Id inIdx = (inArrIdx[2] * globalSize[1] + inArrIdx[1]) * globalSize[0] + inArrIdx[0];
vtkm::Id outIdx =
(outArrIdx[2] * blockSize[1] + outArrIdx[1]) * blockSize[0] + outArrIdx[0];
VTKM_ASSERT(inIdx >= 0 && inIdx < inDataArrayHandle.GetNumberOfValues());
VTKM_ASSERT(outIdx >= 0 && outIdx < nOutValues);
copyIdsPortal.Set(outIdx, inIdx);
}
std::cout << copyIdsPortal.GetNumberOfValues() << std::endl;
vtkm::cont::ArrayHandle<vtkm::Float32> inputArrayHandle;
ds.GetPointField(fieldName).GetData().CopyTo(inputArrayHandle);
auto permutedInArray = make_ArrayHandlePermutation(copyIdsArray, inputArrayHandle);
vtkm::cont::ArrayHandle<vtkm::Float32> outputArrayHandle;
vtkm::cont::ArrayCopy(permutedInArray, outputArrayHandle);
outputArrayHandle.SyncControlArray();
VTKM_ASSERT(outputArrayHandle.GetNumberOfValues() == nOutValues);
auto rp = outputArrayHandle.ReadPortal();
for (vtkm::Id i = 0; i < nOutValues; ++i)
std::cout << rp.Get(i) << " ";
std::cout << std::endl;
vtkm::cont::DataSetBuilderUniform dsb;
if (globalSize[2] <= 1) // 2D Data Set
{
vtkm::Id2 dimensions{ blockSize[0], blockSize[1] };
vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
dataSet.AddPointField(fieldName, outputArrayHandle);
return dataSet;
}
else
{
vtkm::cont::DataSet dataSet = dsb.Create(blockSize);
dataSet.AddPointField(fieldName, outputArrayHandle);
return dataSet;
}
}
std::vector<vtkm::worklet::contourtree_distributed::Edge> ReadGroundTruthContourTree(
std::string filename)
{
std::ifstream ct_file(filename);
vtkm::Id val1, val2;
std::vector<vtkm::worklet::contourtree_distributed::Edge> result;
while (ct_file >> val1 >> val2)
{
result.push_back(vtkm::worklet::contourtree_distributed::Edge(val1, val2));
}
std::sort(result.begin(), result.end());
return result;
}
class TestContourTreeUniformDistributedFilter
{
public:
// numberOfBlocks should be a power of 2
vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(const vtkm::cont::DataSet& ds,
std::string fieldName,
bool useMarchingCubes,
int numberOfBlocks,
int rank = 0,
int numberOfRanks = 1) const
{
// Get dimensions of data set
vtkm::Id3 globalSize;
ds.GetCellSet().CastAndCall(vtkm::worklet::contourtree_augmented::GetPointDimensions(),
globalSize);
// Determine split
vtkm::Id3 blocksPerAxis = ComputeNumberOfBlocksPerAxis(globalSize, numberOfBlocks);
VTKM_ASSERT(numberOfBlocks % numberOfRanks == 0);
vtkm::Id blocksPerRank = numberOfBlocks / numberOfRanks;
vtkm::cont::PartitionedDataSet pds;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksPerRank);
localBlockOrigins.Allocate(blocksPerRank);
localBlockSizes.Allocate(blocksPerRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
for (vtkm::Id blockNo = 0; blockNo < blocksPerRank; ++blockNo)
{
vtkm::Id3 blockOrigin, blockSize, blockIndex;
std::tie(blockIndex, blockOrigin, blockSize) =
ComputeBlockExtents(globalSize, blocksPerAxis, rank * blocksPerRank + blockNo);
pds.AppendPartition(CreateSubDataSet(ds, blockOrigin, blockSize, fieldName));
localBlockOriginsPortal.Set(blockNo, blockOrigin);
localBlockSizesPortal.Set(blockNo, blockSize);
localBlockIndicesPortal.Set(blockNo, blockIndex);
}
vtkm::filter::ContourTreeUniformDistributed filter(blocksPerAxis,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes,
useMarchingCubes);
filter.SetActiveField(fieldName);
std::cout << "Executing filter" << std::endl;
// Execute the contour tree analysis
return filter.Execute(pds);
}
void TestContourTreeUniformDistributed8x9(int nBlocks) const
{
std::cout << "Testing ContourTreeUniformDistributed on 2D 8x9 data set divided into " << nBlocks
<< " blocks." << std::endl;
vtkm::cont::DataSet in_ds = vtkm::cont::testing::MakeTestDataSet().Make2DUniformDataSet3();
vtkm::cont::PartitionedDataSet result =
this->RunContourTreeDUniformDistributed(in_ds, "pointvar", false, nBlocks, 0, 1);
vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
}
treeCompiler.ComputeSuperarcs();
// Print the contour tree we computed
std::cout << "Computed Contour Tree" << std::endl;
treeCompiler.PrintSuperarcs();
// Print the expected contour tree
std::cout << "Expected Contour Tree" << std::endl;
std::cout << " 10 20" << std::endl;
std::cout << " 20 34" << std::endl;
std::cout << " 20 38" << std::endl;
std::cout << " 20 61" << std::endl;
std::cout << " 23 34" << std::endl;
std::cout << " 24 34" << std::endl;
std::cout << " 50 61" << std::endl;
std::cout << " 61 71" << std::endl;
using Edge = vtkm::worklet::contourtree_distributed::Edge;
VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 8),
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 10, 20 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 20, 34 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 20, 38 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 20, 61 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 23, 34 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 24, 34 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 50, 61 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 61, 71 },
"Wrong result for ContourTreeUniformDistributed filter");
}
void TestContourTreeUniformDistributed5x6x7(int nBlocks, bool marchingCubes) const
{
std::cout << "Testing ContourTreeUniformDistributed with "
<< (marchingCubes ? "marching cubes" : "Freudenthal");
std::cout << " mesh connectivity on 3D 5x6x7 data set divided into " << nBlocks << " blocks."
<< std::endl;
vtkm::cont::DataSet in_ds = vtkm::cont::testing::MakeTestDataSet().Make3DUniformDataSet4();
vtkm::cont::PartitionedDataSet result =
this->RunContourTreeDUniformDistributed(in_ds, "pointvar", marchingCubes, nBlocks, 0, 1);
vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
}
treeCompiler.ComputeSuperarcs();
// Print the contour tree we computed
std::cout << "Computed Contour Tree" << std::endl;
treeCompiler.PrintSuperarcs();
// Print the expected contour tree
using Edge = vtkm::worklet::contourtree_distributed::Edge;
std::cout << "Expected Contour Tree" << std::endl;
if (!marchingCubes)
{
std::cout << " 0 112" << std::endl;
std::cout << " 71 72" << std::endl;
std::cout << " 72 78" << std::endl;
std::cout << " 72 101" << std::endl;
std::cout << " 101 112" << std::endl;
std::cout << " 101 132" << std::endl;
std::cout << " 107 112" << std::endl;
std::cout << " 131 132" << std::endl;
std::cout << " 132 138" << std::endl;
VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 9),
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 0, 112 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 71, 72 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 72, 78 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 72, 101 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 101, 112 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 101, 132 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 107, 112 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 131, 132 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[8] == Edge{ 132, 138 },
"Wrong result for ContourTreeUniformDistributed filter");
}
else
{
std::cout << " 0 118" << std::endl;
std::cout << " 31 41" << std::endl;
std::cout << " 41 43" << std::endl;
std::cout << " 41 56" << std::endl;
std::cout << " 56 67" << std::endl;
std::cout << " 56 91" << std::endl;
std::cout << " 62 67" << std::endl;
std::cout << " 67 118" << std::endl;
std::cout << " 81 91" << std::endl;
std::cout << " 91 93" << std::endl;
std::cout << " 118 124" << std::endl;
VTKM_TEST_ASSERT(test_equal(treeCompiler.superarcs.size(), 11),
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[0] == Edge{ 0, 118 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[1] == Edge{ 31, 41 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[2] == Edge{ 41, 43 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[3] == Edge{ 41, 56 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[4] == Edge{ 56, 67 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[5] == Edge{ 56, 91 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[6] == Edge{ 62, 67 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[7] == Edge{ 67, 118 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[8] == Edge{ 81, 91 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[9] == Edge{ 91, 93 },
"Wrong result for ContourTreeUniformDistributed filter");
VTKM_TEST_ASSERT(treeCompiler.superarcs[10] == Edge{ 118, 124 },
"Wrong result for ContourTreeUniformDistributed filter");
}
}
void TestContourTreeFile(std::string ds_filename,
std::string fieldName,
std::string gtct_filename,
int nBlocks,
bool marchingCubes = false) const
{
std::cout << "Testing ContourTreeUniformDistributed with "
<< (marchingCubes ? "marching cubes" : "Freudenthal");
std::cout << " mesh connectivity on \"" << ds_filename << "\" divided into " << nBlocks
<< " blocks." << std::endl;
vtkm::io::VTKDataSetReader reader(ds_filename);
vtkm::cont::DataSet ds;
try
{
ds = reader.ReadDataSet();
}
catch (vtkm::io::ErrorIO& e)
{
std::string message("Error reading: ");
message += ds_filename;
message += ", ";
message += e.GetMessage();
VTKM_TEST_FAIL(message.c_str());
}
std::vector<vtkm::worklet::contourtree_distributed::Edge> groundTruthSuperarcs =
ReadGroundTruthContourTree(gtct_filename);
vtkm::cont::PartitionedDataSet result =
this->RunContourTreeDUniformDistributed(ds, fieldName, marchingCubes, nBlocks, 0, 1);
vtkm::worklet::contourtree_distributed::TreeCompiler treeCompiler;
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)
{
treeCompiler.AddHierarchicalTree(result.GetPartition(ds_no));
}
treeCompiler.ComputeSuperarcs();
if (groundTruthSuperarcs.size() < 50)
{
std::cout << "Computed Contour Tree" << std::endl;
treeCompiler.PrintSuperarcs();
// Print the expected contour tree
std::cout << "Expected Contour Tree" << std::endl;
vtkm::worklet::contourtree_distributed::TreeCompiler::PrintSuperarcArray(
groundTruthSuperarcs);
}
else
{
std::cout << "Not printing computed and expected contour tree due to size." << std::endl;
}
VTKM_TEST_ASSERT(treeCompiler.superarcs == groundTruthSuperarcs,
"Test failed for data set " + ds_filename);
}
void operator()() const
{
using vtkm::cont::testing::Testing;
this->TestContourTreeUniformDistributed8x9(2);
// this->TestContourTreeUniformDistributed8x9(3); CRASH???
this->TestContourTreeUniformDistributed8x9(4);
this->TestContourTreeUniformDistributed8x9(8);
this->TestContourTreeUniformDistributed8x9(16);
this->TestContourTreeFile(Testing::DataPath("rectilinear/vanc.vtk"),
"var",
Testing::DataPath("rectilinear/vanc.ct_txt"),
2);
this->TestContourTreeFile(Testing::DataPath("rectilinear/vanc.vtk"),
"var",
Testing::DataPath("rectilinear/vanc.ct_txt"),
4);
this->TestContourTreeFile(Testing::DataPath("rectilinear/vanc.vtk"),
"var",
Testing::DataPath("rectilinear/vanc.ct_txt"),
8);
this->TestContourTreeFile(Testing::DataPath("rectilinear/vanc.vtk"),
"var",
Testing::DataPath("rectilinear/vanc.ct_txt"),
16);
this->TestContourTreeUniformDistributed5x6x7(2, false);
this->TestContourTreeUniformDistributed5x6x7(4, false);
this->TestContourTreeUniformDistributed5x6x7(8, false);
this->TestContourTreeUniformDistributed5x6x7(16, false);
//this->TestContourTreeUniformDistributed5x6x7(32, false); // Hang???
#if 0
this->TestContourTreeUniformDistributed5x6x7(2, true);
this->TestContourTreeUniformDistributed5x6x7(4, true);
this->TestContourTreeUniformDistributed5x6x7(8, true);
this->TestContourTreeUniformDistributed5x6x7(16, true);
// this->TestContourTreeUniformDistributed5x6x7(32, true); // Hang???
#endif
}
};
}
int UnitTestContourTreeUniformDistributedFilter(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestContourTreeUniformDistributedFilter(), argc, argv);
}

@ -1,11 +1,209 @@
#include "TreeCompiler.h"
#include <iomanip>
#include <vtkm/cont/DataSet.h>
//============================================================================
// 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)
//==============================================================================
#define PRINT_WIDTH 12
#ifndef _TREECOMPILER_H_
#define _TREECOMPILER_H_
#include <iomanip>
#include <iostream>
#include <vtkm/Types.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_distributed
{
// Possibly change the following when comapring to PPP prototype
constexpr int PRINT_WIDTH = 12;
using dataType = vtkm::FloatDefault;
using indexType = vtkm::Id;
// small class for storing the contour arcs
class Edge
{ // Edge
public:
indexType low, high;
// constructor - defaults to -1
Edge(vtkm::Id Low = -1, vtkm::Id High = -1)
: low(Low)
, high(High)
{
}
}; // Edge
// comparison operator <
inline bool operator<(const Edge& LHS, const Edge& RHS)
{ // operator <
#if 0
if (LHS.low < RHS.low) return true;
if (LHS.low > RHS.low) return false;
if (LHS.high < RHS.high) return true;
if (LHS.high > RHS.high) return false;
#endif
if (std::min(LHS.low, LHS.high) < std::min(RHS.low, RHS.high))
return true;
else if (std::min(LHS.low, LHS.high) > std::min(RHS.low, RHS.high))
return false;
if (std::max(LHS.low, LHS.high) < std::max(RHS.low, RHS.high))
return true;
else if (std::max(LHS.low, LHS.high) > std::max(RHS.low, RHS.high))
return false;
return false;
} // operator <
// comparison operator ==
inline bool operator==(const Edge& LHS, const Edge& RHS)
{ // operator ==
return (LHS.low == RHS.low && LHS.high == RHS.high) ||
(LHS.low == RHS.high && LHS.high == RHS.low);
} // operator ==
// a helper class which stores a single supernode inserted onto a superarc
class SupernodeOnSuperarc
{ // class SupernodeOnSuperarc
public:
// the global ID of the supernode
indexType globalID;
// the data value stored at the supernode
dataType dataValue;
// the low and high ends of the superarc it is on (may be itself)
indexType lowEnd, highEnd;
// constructor
SupernodeOnSuperarc(indexType GlobalID = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
dataType DataValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType LowEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT,
indexType HighEnd = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)
: globalID(GlobalID)
, dataValue(DataValue)
, lowEnd(LowEnd)
, highEnd(HighEnd)
{ // constructor
} // constructor
}; // class SupernodeOnSuperarc
// overloaded comparison operator
// primary sort is by superarc (low, high),
// then secondary sort on datavalue
// tertiary on globalID to implement simulated simplicity
inline bool operator<(const SupernodeOnSuperarc& left, const SupernodeOnSuperarc& right)
{ // < operator
// simple lexicographic sort
if (left.lowEnd < right.lowEnd)
return true;
if (left.lowEnd > right.lowEnd)
return false;
if (left.highEnd < right.highEnd)
return true;
if (left.highEnd > right.highEnd)
return false;
if (left.dataValue < right.dataValue)
return true;
if (left.dataValue > right.dataValue)
return false;
if (left.globalID < right.globalID)
return true;
if (left.globalID > right.globalID)
return false;
// fall-through (shouldn't happen, but)
// if they're the same, it's false
return false;
} // < operator
// stream output
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node)
std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node);
// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node);
// the class that compiles the contour tree
class TreeCompiler
{ // class TreeCompiler
public:
// we want a vector of supernodes on superarcs
std::vector<SupernodeOnSuperarc> supernodes;
// and a vector of Edges (the output)
std::vector<Edge> superarcs;
// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
void AddHierarchicalTree(const vtkm::cont::DataSet& addedTree);
// routine to compute the actual superarcs
void ComputeSuperarcs();
// routine to print a superarcs array in our format
static void PrintSuperarcArray(const std::vector<Edge>& superarc_array);
// routine to print the superarcs
void PrintSuperarcs() const;
// routine to write out binary file
void WriteBinary(FILE* outFile) const;
// routine to read in binary file & append to contents
void ReadBinary(FILE* inFile);
}; // class TreeCompiler
// stream output
inline std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node)
{ // stream output
outStream << node.lowEnd << " " << node.highEnd << " " << node.dataValue << " " << node.globalID
<< std::endl;
@ -13,7 +211,7 @@ std::ostream& operator<<(std::ostream& outStream, SupernodeOnSuperarc& node)
} // stream output
// stream input
std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node)
inline std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node)
{ // stream input
inStream >> node.lowEnd >> node.highEnd >> node.dataValue >> node.globalID;
return inStream;
@ -21,7 +219,7 @@ std::istream& operator>>(std::istream& inStream, SupernodeOnSuperarc& node)
// routine to add a known hierarchical tree to it
// note that this DOES NOT finalise - we don't want too many sorts
void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
inline void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
{ // TreeCompiler::AddHierarchicalTree()
// Copy relevant tree content to STL arrays
vtkm::cont::VariantArrayHandle dataValues_array = addedTree.GetField("DataValues").GetData();
@ -59,7 +257,8 @@ void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
superparents_handle.SyncControlArray(); //Forces values to get updated if copy happened on GPU
// loop through all of the supernodes in the hierarchical tree
for (indexType supernode = 0; supernode < added_tree_supernodes.size(); supernode++)
for (indexType supernode = 0; supernode < static_cast<indexType>(added_tree_supernodes.size());
supernode++)
{ // per supernode
// retrieve the regular ID for the supernode
indexType regularId = added_tree_supernodes[supernode];
@ -129,7 +328,7 @@ void TreeCompiler::AddHierarchicalTree(const vtkm::cont::DataSet& addedTree)
} // TreeCompiler::AddHierarchicalTree()
// routine to compute the actual superarcs
void TreeCompiler::ComputeSuperarcs()
inline void TreeCompiler::ComputeSuperarcs()
{ // TreeCompiler::ComputeSuperarcs()
// first we sort the vector
std::sort(supernodes.begin(), supernodes.end());
@ -140,7 +339,8 @@ void TreeCompiler::ComputeSuperarcs()
// this is because we know a priori that the last one is the last supernode on a superarc
// and would fail the test inside the loop. By putting it in the loop test, we avoid having
// to have an explicit if statement inside the loop
for (indexType supernode = 0; supernode < supernodes.size() - 1; supernode++)
for (indexType supernode = 0; supernode < static_cast<vtkm::Id>(supernodes.size() - 1);
supernode++)
{ // loop through supernodes
// this is actually painfully simple: if the (lowEnd, highEnd) don't match the next one,
// then we're at the end of the group and do nothing. Otherwise, we link to the next one
@ -161,30 +361,32 @@ void TreeCompiler::ComputeSuperarcs()
} // TreeCompiler::ComputeSuperarcs()
// routine to print the superarcs
void TreeCompiler::PrintSuperarcs()
{ // TreeCompiler::PrintSuperarcs()
std::cout << "============" << std::endl;
std::cout << "Contour Tree" << std::endl;
for (indexType superarc = 0; superarc < superarcs.size(); superarc++)
inline void TreeCompiler::PrintSuperarcArray(const std::vector<Edge>& superarc_array)
{ // TreeCompiler::PrintSuperarcArray()
for (indexType superarc = 0; superarc < static_cast<indexType>(superarc_array.size()); superarc++)
{ // per superarc
if (superarcs[superarc].low < superarcs[superarc].high)
if (superarc_array[superarc].low < superarc_array[superarc].high)
{ // order by ID not value
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << " ";
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << std::endl;
std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].low << " ";
std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].high << std::endl;
} // order by ID not value
else
{ // order by ID not value
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].high << " ";
std::cout << std::setw(PRINT_WIDTH) << superarcs[superarc].low << std::endl;
std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].high << " ";
std::cout << std::setw(PRINT_WIDTH) << superarc_array[superarc].low << std::endl;
} // order by ID not value
} // per superarc
} // TreeCompiler::PrintSuperarcs()
} // TreeCompiler::PrintSuperarcArray()
inline void TreeCompiler::PrintSuperarcs() const
{
PrintSuperarcArray(this->superarcs);
}
// routine to write out binary file
void TreeCompiler::WriteBinary(FILE* outFile)
inline void TreeCompiler::WriteBinary(FILE* outFile) const
{ // WriteBinary()
// do a bulk write of the entire contents
// no error checking, no type checking, no nothing
@ -192,7 +394,7 @@ void TreeCompiler::WriteBinary(FILE* outFile)
} // WriteBinary()
// routine to read in binary file and append
void TreeCompiler::ReadBinary(FILE* inFile)
inline void TreeCompiler::ReadBinary(FILE* inFile)
{ // ReadBinary()
// use fseek to jump to the end
fseek(inFile, 0, SEEK_END);
@ -217,15 +419,16 @@ void TreeCompiler::ReadBinary(FILE* inFile)
} // ReadBinary()
// stream output - just dumps the supernodeonsuperarcs
std::ostream& operator<<(std::ostream& outStream, TreeCompiler& tree)
inline std::ostream& operator<<(std::ostream& outStream, TreeCompiler& tree)
{ // stream output
for (indexType supernode = 0; supernode < tree.supernodes.size(); supernode++)
for (indexType supernode = 0; supernode < static_cast<indexType>(tree.supernodes.size());
supernode++)
outStream << tree.supernodes[supernode];
return outStream;
} // stream output
// stream input - reads in the supernodeonsuperarcs & appends them
std::istream& operator>>(std::istream& inStream, TreeCompiler& tree)
inline std::istream& operator>>(std::istream& inStream, TreeCompiler& tree)
{ // stream input
while (!inStream.eof())
{
@ -237,3 +440,9 @@ std::istream& operator>>(std::istream& inStream, TreeCompiler& tree)
tree.supernodes.resize(tree.supernodes.size() - 1);
return inStream;
} // stream input
} // namespace contourtree_distributed
} // namespace worklet
} // namespace vtkm
#endif