Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compute correct histogram for distributed multiblock cases.

Browse Source 
Histogram filter now correctly computes the histogram for multiblock
input or when operating in a distributed environment. It employs
reduction to reduce local/per block results to generate the final
output.
This commit is contained in:
Utkarsh Ayachit 2018-04-09 15:10:30 -04:00
parent 0af0271f59
commit 89041e4f8e
6 changed files with 397 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
vtkm/filter/Histogram.h
View File

@ -27,9 +27,10 @@ namespace vtkm
{
namespace filter
{
/// \brief Construct the histogram of a given Field
///
/// Construct a histogram with a default of 10 bins
/// Construct a histogram with a default of 10 bins.
///
class Histogram : public vtkm::filter::FilterField<Histogram>
{
@ -44,14 +45,26 @@ public:
VTKM_CONT
vtkm::Id GetNumberOfBins() const { return this->NumberOfBins; }
//Returns the bin delta of the last computed field, be it from DoExecute
//or from MapField
//@{
/// Get/Set the range to use to generate the histogram. If range is set to
/// empty, the field's global range (computed using `vtkm::cont::FieldRangeGlobalCompute`)
/// will be used.
VTKM_CONT
void SetRange(const vtkm::Range& range) { this->Range = range; }
VTKM_CONT
const vtkm::Range& GetRange() const { return this->Range; }
//@}
/// Returns the bin delta of the last computed field.
VTKM_CONT
vtkm::Float64 GetBinDelta() const { return this->BinDelta; }
//Returns the the min and max values for that last computed field
/// Returns the range used for most recent execute. If `SetRange` is used to
/// specify and non-empty range, then this will be same as the range after
/// the `Execute` call.
VTKM_CONT
vtkm::Range GetDataRange() const { return this->DataRange; }
vtkm::Range GetComputedRange() const { return this->ComputedRange; }
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
@ -60,10 +73,25 @@ public:
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
const DeviceAdapter& tag);
//@{
/// when operating on vtkm::cont::MultiBlock, we
/// want to do processing across ranks as well. Just adding pre/post handles
/// for the same does the trick.
template <typename DerivedPolicy>
VTKM_CONT void PreExecute(const vtkm::cont::MultiBlock& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
template <typename DerivedPolicy>
VTKM_CONT void PostExecute(const vtkm::cont::MultiBlock& input,
vtkm::cont::MultiBlock& output,
const vtkm::filter::PolicyBase<DerivedPolicy>&);
//@}
private:
vtkm::Id NumberOfBins;
vtkm::Float64 BinDelta;
vtkm::Range DataRange;
vtkm::Range ComputedRange;
vtkm::Range Range;
};
template <>

231
vtkm/filter/Histogram.hxx
View File

@ -18,20 +18,175 @@
// this software.
//============================================================================
#include <vtkm/filter/internal/CreateResult.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/FieldHistogram.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/AssignerMultiBlock.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/FieldRangeGlobalCompute.h>
#include <vtkm/cont/diy/Serialization.h>
#include <vtkm/filter/internal/CreateResult.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/Configure.h>
#include VTKM_DIY(diy/decomposition.hpp)
#include VTKM_DIY(diy/master.hpp)
#include VTKM_DIY(diy/partners/broadcast.hpp)
#include VTKM_DIY(diy/partners/swap.hpp)
#include VTKM_DIY(diy/reduce.hpp)
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace vtkm
{
namespace filter
{
namespace detail
{
class DistributedHistogram
{
class Reducer
{
public:
void operator()(vtkm::cont::ArrayHandle<vtkm::Id>* result,
const diy::ReduceProxy& srp,
const diy::RegularMergePartners&) const
{
const auto selfid = srp.gid();
// 1. dequeue.
std::vector<int> incoming;
srp.incoming(incoming);
for (const int gid : incoming)
{
if (gid != selfid)
{
vtkm::cont::ArrayHandle<vtkm::Id> in;
srp.dequeue(gid, in);
if (result->GetNumberOfValues() == 0)
{
*result = in;
}
else
{
vtkm::cont::Algorithm::Transform(*result, in, *result, vtkm::Add());
}
}
}
// 2. enqueue
for (int cc = 0; cc < srp.out_link().size(); ++cc)
{
auto target = srp.out_link().target(cc);
if (target.gid != selfid)
{
srp.enqueue(target, *result);
}
}
}
};
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> LocalBlocks;
public:
DistributedHistogram(vtkm::Id numLocalBlocks)
: LocalBlocks(static_cast<size_t>(numLocalBlocks))
{
}
void SetLocalHistogram(vtkm::Id index, const vtkm::cont::ArrayHandle<vtkm::Id>& bins)
{
this->LocalBlocks[static_cast<size_t>(index)] = bins;
}
void SetLocalHistogram(vtkm::Id index, const vtkm::cont::Field& field)
{
this->SetLocalHistogram(index, field.GetData().Cast<vtkm::cont::ArrayHandle<vtkm::Id>>());
}
vtkm::cont::ArrayHandle<vtkm::Id> ReduceAll() const
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
const vtkm::Id numLocalBlocks = static_cast<vtkm::Id>(this->LocalBlocks.size());
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
if (comm.size() == 1 && numLocalBlocks <= 1)
{
// no reduction necessary.
return numLocalBlocks == 0 ? ArrayType() : this->LocalBlocks[0];
}
diy::Master master(
comm,
/*threads*/ 1,
/*limit*/ -1,
[]() -> void* { return new vtkm::cont::ArrayHandle<vtkm::Id>(); },
[](void* ptr) { delete static_cast<vtkm::cont::ArrayHandle<vtkm::Id>*>(ptr); });
vtkm::cont::AssignerMultiBlock assigner(numLocalBlocks);
diy::RegularDecomposer<diy::DiscreteBounds> decomposer(
/*dims*/ 1, diy::interval(0, assigner.nblocks() - 1), assigner.nblocks());
decomposer.decompose(comm.rank(), assigner, master);
assert(static_cast<vtkm::Id>(master.size()) == numLocalBlocks);
for (vtkm::Id cc = 0; cc < numLocalBlocks; ++cc)
{
*master.block<ArrayType>(static_cast<int>(cc)) = this->LocalBlocks[static_cast<size_t>(cc)];
}
diy::RegularMergePartners partners(decomposer, /*k=*/2);
// reduce to block-0.
diy::reduce(master, assigner, partners, Reducer());
ArrayType result;
if (master.local(0))
{
result = *master.block<ArrayType>(master.lid(0));
}
this->Broadcast(result);
return result;
}
private:
void Broadcast(vtkm::cont::ArrayHandle<vtkm::Id>& data) const
{
// broadcast to all ranks (and not blocks).
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
if (comm.size() > 1)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
diy::Master master(
comm,
/*threads*/ 1,
/*limit*/ -1,
[]() -> void* { return new vtkm::cont::ArrayHandle<vtkm::Id>(); },
[](void* ptr) { delete static_cast<vtkm::cont::ArrayHandle<vtkm::Id>*>(ptr); });
diy::ContiguousAssigner assigner(comm.size(), comm.size());
diy::RegularDecomposer<diy::DiscreteBounds> decomposer(
1, diy::interval(0, comm.size() - 1), comm.size());
decomposer.decompose(comm.rank(), assigner, master);
assert(master.size() == 1); // number of local blocks should be 1 per rank.
*master.block<ArrayType>(0) = data;
diy::RegularBroadcastPartners partners(decomposer, /*k=*/2);
diy::reduce(master, assigner, partners, Reducer());
data = *master.block<ArrayType>(0);
}
}
};
} // namespace detail
//-----------------------------------------------------------------------------
inline VTKM_CONT Histogram::Histogram()
: NumberOfBins(10)
, BinDelta(0)
, DataRange()
, ComputedRange()
, Range()
{
this->SetOutputFieldName("histogram");
}
@ -39,9 +194,9 @@ inline VTKM_CONT Histogram::Histogram()
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
inline VTKM_CONT vtkm::cont::DataSet Histogram::DoExecute(
const vtkm::cont::DataSet& inDataSet,
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::FieldMetadata&,
const vtkm::filter::PolicyBase<DerivedPolicy>&,
const DeviceAdapter& device)
{
@ -49,14 +204,70 @@ inline VTKM_CONT vtkm::cont::DataSet Histogram::DoExecute(
T delta;
vtkm::worklet::FieldHistogram worklet;
worklet.Run(field, this->NumberOfBins, this->DataRange, delta, binArray, device);
if (this->ComputedRange.IsNonEmpty())
{
worklet.Run(field,
this->NumberOfBins,
static_cast<T>(this->ComputedRange.Min),
static_cast<T>(this->ComputedRange.Max),
delta,
binArray,
device);
}
else
{
worklet.Run(field, this->NumberOfBins, this->ComputedRange, delta, binArray, device);
}
this->BinDelta = static_cast<vtkm::Float64>(delta);
return internal::CreateResult(inDataSet,
binArray,
this->GetOutputFieldName(),
fieldMetadata.GetAssociation(),
fieldMetadata.GetCellSetName());
vtkm::cont::DataSet output;
vtkm::cont::Field rfield(
this->GetOutputFieldName(), vtkm::cont::Field::ASSOC_WHOLE_MESH, binArray);
output.AddField(rfield);
return output;
}
//-----------------------------------------------------------------------------
template <typename DerivedPolicy>
inline VTKM_CONT void Histogram::PreExecute(const vtkm::cont::MultiBlock& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
if (this->Range.IsNonEmpty())
{
this->ComputedRange = this->Range;
}
else
{
auto handle = vtkm::cont::FieldRangeGlobalCompute(
input, this->GetActiveFieldName(), this->GetActiveFieldAssociation());
if (handle.GetNumberOfValues() != 1)
{
throw vtkm::cont::ErrorFilterExecution("expecting scalar field.");
}
this->ComputedRange = handle.GetPortalConstControl().Get(0);
}
}
//-----------------------------------------------------------------------------
template <typename DerivedPolicy>
inline VTKM_CONT void Histogram::PostExecute(const vtkm::cont::MultiBlock&,
vtkm::cont::MultiBlock& result,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
// iterate and compute histogram for each local block.
detail::DistributedHistogram helper(result.GetNumberOfBlocks());
for (vtkm::Id cc = 0; cc < result.GetNumberOfBlocks(); ++cc)
{
auto& ablock = result.GetBlock(cc);
helper.SetLocalHistogram(cc, ablock.GetField(this->GetOutputFieldName()));
}
vtkm::cont::DataSet output;
vtkm::cont::Field rfield(
this->GetOutputFieldName(), vtkm::cont::Field::ASSOC_WHOLE_MESH, helper.ReduceAll());
output.AddField(rfield);
result = vtkm::cont::MultiBlock(output);
}
}
} // namespace vtkm::filter

1
vtkm/filter/testing/CMakeLists.txt
View File

@ -40,6 +40,7 @@ set(unit_tests
UnitTestMaskFilter.cxx
UnitTestMaskPointsFilter.cxx
UnitTestMultiBlockFilters.cxx
UnitTestMultiBlockHistogramFilter.cxx
UnitTestNDEntropyFilter.cxx
UnitTestNDHistogramFilter.cxx
UnitTestPointAverageFilter.cxx

8
vtkm/filter/testing/UnitTestHistogramFilter.cxx
View File

@ -311,28 +311,28 @@ void TestHistogram()
histogram.SetActiveField("p_poisson");
auto result = histogram.Execute(ds);
delta = histogram.GetBinDelta();
range = histogram.GetDataRange();
range = histogram.GetComputedRange();
VerifyHistogram(result, histogram.GetNumberOfBins(), range, delta);
histogram.SetNumberOfBins(100);
histogram.SetActiveField("p_normal");
result = histogram.Execute(ds);
delta = histogram.GetBinDelta();
range = histogram.GetDataRange();
range = histogram.GetComputedRange();
VerifyHistogram(result, histogram.GetNumberOfBins(), range, delta, false);
histogram.SetNumberOfBins(1);
histogram.SetActiveField("p_chiSquare");
result = histogram.Execute(ds);
delta = histogram.GetBinDelta();
range = histogram.GetDataRange();
range = histogram.GetComputedRange();
VerifyHistogram(result, histogram.GetNumberOfBins(), range, delta);
histogram.SetNumberOfBins(1000000);
histogram.SetActiveField("p_uniform");
result = histogram.Execute(ds);
delta = histogram.GetBinDelta();
range = histogram.GetDataRange();
range = histogram.GetComputedRange();
VerifyHistogram(result, histogram.GetNumberOfBins(), range, delta, false);
} // TestFieldHistogram

9
vtkm/filter/testing/UnitTestMultiBlockFilters.cxx
View File

@ -33,7 +33,6 @@
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/CellAverage.h>
#include <vtkm/filter/Histogram.h>
template <typename T>
@ -85,7 +84,7 @@ vtkm::cont::MultiBlock MultiBlockBuilder(std::size_t BlockNum, std::string Field
}
template <typename D>
void Result_Verify(const vtkm::cont::MultiBlock& Result,
D Filter,
D& Filter,
const vtkm::cont::MultiBlock& Blocks,
std::string FieldName)
{
@ -124,12 +123,6 @@ void TestMultiBlockFilters()
vtkm::cont::MultiBlock result;
vtkm::cont::MultiBlock Blocks;
Blocks = MultiBlockBuilder<vtkm::Float64>(BlockNum, "cellvar");
vtkm::filter::Histogram histogram;
histogram.SetActiveField("cellvar");
result = histogram.Execute(Blocks);
Result_Verify(result, histogram, Blocks, std::string("cellvar"));
Blocks = MultiBlockBuilder<vtkm::Id>(BlockNum, "pointvar");
vtkm::filter::CellAverage cellAverage;
cellAverage.SetOutputFieldName("average");

136
vtkm/filter/testing/UnitTestMultiBlockHistogramFilter.cxx Normal file
View File

@ -0,0 +1,136 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <vtkm/filter/Histogram.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/MultiBlock.h>
#include <vtkm/cont/testing/Testing.h>
#include <algorithm>
#include <numeric>
#include <random>
#include <utility>
namespace
{
static unsigned int uid = 1;
template <typename T>
vtkm::cont::ArrayHandle<T> CreateArrayHandle(T min, T max, vtkm::Id numVals)
{
std::mt19937 gen(uid++);
std::uniform_real_distribution<double> dis(static_cast<double>(min), static_cast<double>(max));
vtkm::cont::ArrayHandle<T> handle;
handle.Allocate(numVals);
std::generate(vtkm::cont::ArrayPortalToIteratorBegin(handle.GetPortalControl()),
vtkm::cont::ArrayPortalToIteratorEnd(handle.GetPortalControl()),
[&]() { return static_cast<T>(dis(gen)); });
return handle;
}
template <typename T, int size>
vtkm::cont::ArrayHandle<vtkm::Vec<T, size>> CreateArrayHandle(const vtkm::Vec<T, size>& min,
const vtkm::Vec<T, size>& max,
vtkm::Id numVals)
{
std::mt19937 gen(uid++);
std::uniform_real_distribution<double> dis[size];
for (int cc = 0; cc < size; ++cc)
{
dis[cc] = std::uniform_real_distribution<double>(static_cast<double>(min[cc]),
static_cast<double>(max[cc]));
}
vtkm::cont::ArrayHandle<T> handle;
handle.Allocate(numVals);
std::generate(vtkm::cont::ArrayPortalToIteratorBegin(handle.GetPortalControl()),
vtkm::cont::ArrayPortalToIteratorEnd(handle.GetPortalControl()),
[&]() {
vtkm::Vec<T, size> val;
for (int cc = 0; cc < size; ++cc)
{
val[cc] = static_cast<T>(dis[cc](gen));
}
return val;
});
return handle;
}
template <typename T>
void AddField(vtkm::cont::DataSet& dataset,
const T& min,
const T& max,
vtkm::Id numVals,
const std::string& name,
vtkm::cont::Field::AssociationEnum assoc = vtkm::cont::Field::ASSOC_POINTS)
{
auto ah = CreateArrayHandle(min, max, numVals);
dataset.AddField(vtkm::cont::Field(name, assoc, ah));
}
}
static void TestMultiBlockHistogram()
{
// init random seed.
std::srand(100);
vtkm::cont::MultiBlock mb;
vtkm::cont::DataSet block0;
AddField<double>(block0, 0.0, 100.0, 1024, "double");
mb.AddBlock(block0);
vtkm::cont::DataSet block1;
AddField<int>(block1, 100, 1000, 1024, "double");
mb.AddBlock(block1);
vtkm::cont::DataSet block2;
AddField<double>(block2, 100.0, 500.0, 1024, "double");
mb.AddBlock(block2);
vtkm::filter::Histogram histogram;
histogram.SetActiveField("double");
auto result = histogram.Execute(mb);
VTKM_TEST_ASSERT(result.GetNumberOfBlocks() == 1, "Expecting 1 block.");
auto bins =
result.GetBlock(0).GetField("histogram").GetData().Cast<vtkm::cont::ArrayHandle<vtkm::Id>>();
VTKM_TEST_ASSERT(bins.GetNumberOfValues() == 10, "Expecting 10 bins.");
auto count = std::accumulate(vtkm::cont::ArrayPortalToIteratorBegin(bins.GetPortalConstControl()),
vtkm::cont::ArrayPortalToIteratorEnd(bins.GetPortalConstControl()),
vtkm::Id(0),
vtkm::Add());
VTKM_TEST_ASSERT(count == 1024 * 3, "Expecting 3072 values");
std::cout << "Values [" << count << "] =";
for (int cc = 0; cc < 10; ++cc)
{
std::cout << " " << bins.GetPortalConstControl().Get(cc);
}
std::cout << std::endl;
};
int UnitTestMultiBlockHistogramFilter(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestMultiBlockHistogram);
}