Merge branch 'master' into kdtree_docs

This commit is contained in:
Li-Ta Lo 2017-09-18 13:19:55 -06:00
commit 1eaaf342ef
33 changed files with 2561 additions and 45 deletions

@ -28,7 +28,7 @@ add_subdirectory(contour_tree)
add_subdirectory(cosmotools)
add_subdirectory(demo)
add_subdirectory(dynamic_dispatcher)
#add_subdirectory(game_of_life)
add_subdirectory(game_of_life)
add_subdirectory(hello_world)
add_subdirectory(isosurface)
add_subdirectory(multi_backend)

@ -32,7 +32,7 @@
#include <stdexcept>
#include <string>
typedef vtkm::Vec<vtkm::Float32, 3> FloatVec3;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
void TestCosmoCenterFinder(const char* fileName)
{
@ -129,9 +129,6 @@ int main(int argc, char* argv[])
return 1;
}
std::cout << "Device Adapter Name: " << vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetName()
<< std::endl;
TestCosmoCenterFinder(argv[1]);
return 0;

@ -32,7 +32,7 @@
#include <stdexcept>
#include <string>
typedef vtkm::Vec<vtkm::Float32, 3> FloatVec3;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
void TestCosmoHaloFinder(const char* fileName)
{
@ -123,9 +123,6 @@ int main(int argc, char* argv[])
return 1;
}
std::cout << "Device Adapter Name: " << vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetName()
<< std::endl;
TestCosmoHaloFinder(argv[1]);
return 0;

@ -168,7 +168,7 @@ void RunTest(const std::string& fname,
else
{
vtkm::worklet::Streamline streamline;
streamline.Run(rk4, seedArray, fieldArray, numSteps, DeviceAdapter());
streamline.Run(rk4, seedArray, numSteps, DeviceAdapter());
}
auto t1 = std::chrono::high_resolution_clock::now() - t0;

@ -203,7 +203,7 @@ struct DeviceAdapterAlgorithm
VTKM_CONT static void ReduceByKey(const vtkm::cont::ArrayHandle<T, CKeyIn>& keys,
const vtkm::cont::ArrayHandle<U, CValIn>& values,
vtkm::cont::ArrayHandle<T, CKeyOut>& keys_output,
vtkm::cont::ArrayHandle<T, CValOut>& values_output,
vtkm::cont::ArrayHandle<U, CValOut>& values_output,
BinaryFunctor binary_functor);
/// \brief Compute an inclusive prefix sum operation on the input ArrayHandle.

@ -60,6 +60,31 @@ public:
input.PrepareForInput(vtkm::cont::DeviceAdapterTagTBB()), initialValue, binary_functor);
}
template <typename T,
typename U,
class CKeyIn,
class CValIn,
class CKeyOut,
class CValOut,
class BinaryFunctor>
VTKM_CONT static void ReduceByKey(const vtkm::cont::ArrayHandle<T, CKeyIn>& keys,
const vtkm::cont::ArrayHandle<U, CValIn>& values,
vtkm::cont::ArrayHandle<T, CKeyOut>& keys_output,
vtkm::cont::ArrayHandle<U, CValOut>& values_output,
BinaryFunctor binary_functor)
{
vtkm::Id inputSize = keys.GetNumberOfValues();
VTKM_ASSERT(inputSize == values.GetNumberOfValues());
vtkm::Id outputSize =
tbb::ReduceByKeyPortals(keys.PrepareForInput(DeviceAdapterTagTBB()),
values.PrepareForInput(DeviceAdapterTagTBB()),
keys_output.PrepareForOutput(inputSize, DeviceAdapterTagTBB()),
values_output.PrepareForOutput(inputSize, DeviceAdapterTagTBB()),
binary_functor);
keys_output.Shrink(outputSize);
values_output.Shrink(outputSize);
}
template <typename T, class CIn, class COut>
VTKM_CONT static T ScanInclusive(const vtkm::cont::ArrayHandle<T, CIn>& input,
vtkm::cont::ArrayHandle<T, COut>& output)

@ -27,6 +27,9 @@
#include <vtkm/cont/internal/FunctorsGeneral.h>
#include <vtkm/exec/internal/ErrorMessageBuffer.h>
#include <algorithm>
#include <sstream>
VTKM_THIRDPARTY_PRE_INCLUDE
#if defined(VTKM_MSVC)
@ -195,6 +198,359 @@ VTKM_CONT static T ReducePortals(InputPortalType inputPortal,
}
}
// Define this to print out timing information from the reduction and join
// operations in the tbb ReduceByKey algorithm:
//#define _VTKM_DEBUG_TBB_RBK
template <typename KeysInPortalType,
typename ValuesInPortalType,
typename KeysOutPortalType,
typename ValuesOutPortalType,
class BinaryOperationType>
struct ReduceByKeyBody
{
using KeyType = typename KeysInPortalType::ValueType;
using ValueType = typename ValuesInPortalType::ValueType;
struct Range
{
vtkm::Id InputBegin;
vtkm::Id InputEnd;
vtkm::Id OutputBegin;
vtkm::Id OutputEnd;
VTKM_EXEC_CONT
Range()
: InputBegin(-1)
, InputEnd(-1)
, OutputBegin(-1)
, OutputEnd(-1)
{
}
VTKM_EXEC_CONT
Range(vtkm::Id inputBegin, vtkm::Id inputEnd, vtkm::Id outputBegin, vtkm::Id outputEnd)
: InputBegin(inputBegin)
, InputEnd(inputEnd)
, OutputBegin(outputBegin)
, OutputEnd(outputEnd)
{
this->AssertSane();
}
VTKM_EXEC_CONT
void AssertSane() const
{
VTKM_ASSERT("Input begin precedes end" && this->InputBegin <= this->InputEnd);
VTKM_ASSERT("Output begin precedes end" && this->OutputBegin <= this->OutputEnd);
VTKM_ASSERT("Output not past input" && this->OutputBegin <= this->InputBegin &&
this->OutputEnd <= this->InputEnd);
VTKM_ASSERT("Output smaller than input" &&
(this->OutputEnd - this->OutputBegin) <= (this->InputEnd - this->InputBegin));
}
VTKM_EXEC_CONT
bool IsNext(const Range& next) const { return this->InputEnd == next.InputBegin; }
};
KeysInPortalType KeysInPortal;
ValuesInPortalType ValuesInPortal;
KeysOutPortalType KeysOutPortal;
ValuesOutPortalType ValuesOutPortal;
BinaryOperationType BinaryOperation;
Range Ranges;
#ifdef _VTKM_DEBUG_TBB_RBK
double ReduceTime;
double JoinTime;
#endif
VTKM_CONT
ReduceByKeyBody(const KeysInPortalType& keysInPortal,
const ValuesInPortalType& valuesInPortal,
const KeysOutPortalType& keysOutPortal,
const ValuesOutPortalType& valuesOutPortal,
BinaryOperationType binaryOperation)
: KeysInPortal(keysInPortal)
, ValuesInPortal(valuesInPortal)
, KeysOutPortal(keysOutPortal)
, ValuesOutPortal(valuesOutPortal)
, BinaryOperation(binaryOperation)
#ifdef _VTKM_DEBUG_TBB_RBK
, ReduceTime(0)
, JoinTime(0)
#endif
{
}
VTKM_EXEC_CONT
ReduceByKeyBody(const ReduceByKeyBody& body, ::tbb::split)
: KeysInPortal(body.KeysInPortal)
, ValuesInPortal(body.ValuesInPortal)
, KeysOutPortal(body.KeysOutPortal)
, ValuesOutPortal(body.ValuesOutPortal)
, BinaryOperation(body.BinaryOperation)
#ifdef _VTKM_DEBUG_TBB_RBK
, ReduceTime(0)
, JoinTime(0)
#endif
{
}
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC
void operator()(const ::tbb::blocked_range<vtkm::Id>& range)
{
#ifdef _VTKM_DEBUG_TBB_RBK
::tbb::tick_count startTime = ::tbb::tick_count::now();
#endif // _VTKM_DEBUG_TBB_RBK
if (range.empty())
{
return;
}
bool firstRun = this->Ranges.OutputBegin < 0; // First use of this body object
if (firstRun)
{
this->Ranges.InputBegin = range.begin();
}
else
{
// Must be a continuation of the previous input range:
VTKM_ASSERT(this->Ranges.InputEnd == range.begin());
}
this->Ranges.InputEnd = range.end();
this->Ranges.AssertSane();
using KeysInIteratorsType = vtkm::cont::ArrayPortalToIterators<KeysInPortalType>;
using ValuesInIteratorsType = vtkm::cont::ArrayPortalToIterators<ValuesInPortalType>;
using KeysOutIteratorsType = vtkm::cont::ArrayPortalToIterators<KeysOutPortalType>;
using ValuesOutIteratorsType = vtkm::cont::ArrayPortalToIterators<ValuesOutPortalType>;
KeysInIteratorsType keysInIters(this->KeysInPortal);
ValuesInIteratorsType valuesInIters(this->ValuesInPortal);
KeysOutIteratorsType keysOutIters(this->KeysOutPortal);
ValuesOutIteratorsType valuesOutIters(this->ValuesOutPortal);
using KeysInIteratorType = typename KeysInIteratorsType::IteratorType;
using ValuesInIteratorType = typename ValuesInIteratorsType::IteratorType;
using KeysOutIteratorType = typename KeysOutIteratorsType::IteratorType;
using ValuesOutIteratorType = typename ValuesOutIteratorsType::IteratorType;
KeysInIteratorType keysIn = keysInIters.GetBegin();
ValuesInIteratorType valuesIn = valuesInIters.GetBegin();
KeysOutIteratorType keysOut = keysOutIters.GetBegin();
ValuesOutIteratorType valuesOut = valuesOutIters.GetBegin();
vtkm::Id readPos = range.begin();
const vtkm::Id readEnd = range.end();
// Determine output index. If we're reusing the body, pick up where the
// last block left off. If not, use the input range.
vtkm::Id writePos;
if (firstRun)
{
this->Ranges.OutputBegin = range.begin();
this->Ranges.OutputEnd = range.begin();
writePos = range.begin();
}
else
{
writePos = this->Ranges.OutputEnd;
}
this->Ranges.AssertSane();
// We're either writing at the end of a previous block, or at the input
// location. Either way, the write position will never be greater than
// the read position.
VTKM_ASSERT(writePos <= readPos);
// Initialize reduction variables:
BinaryOperationType functor(this->BinaryOperation);
KeyType currentKey = keysIn[readPos];
ValueType currentValue = valuesIn[readPos];
++readPos;
// If the start of the current range continues a previous key block,
// initialize with the previous result and decrement the write index.
VTKM_ASSERT(firstRun || writePos > 0);
if (!firstRun && keysOut[writePos - 1] == currentKey)
{
// Ensure that we'll overwrite the continued key values:
--writePos;
// Update our accumulator with the partial value:
currentValue = functor(valuesOut[writePos], currentValue);
}
// Special case: single value in range
if (readPos >= readEnd)
{
keysOut[writePos] = currentKey;
valuesOut[writePos] = currentValue;
++writePos;
this->Ranges.OutputEnd = writePos;
return;
}
for (;;)
{
while (readPos < readEnd && currentKey == keysIn[readPos])
{
currentValue = functor(currentValue, valuesIn[readPos]);
++readPos;
}
VTKM_ASSERT(writePos <= readPos);
keysOut[writePos] = currentKey;
valuesOut[writePos] = currentValue;
++writePos;
if (readPos < readEnd)
{
currentKey = keysIn[readPos];
currentValue = valuesIn[readPos];
++readPos;
continue;
}
break;
}
this->Ranges.OutputEnd = writePos;
#ifdef _VTKM_DEBUG_TBB_RBK
::tbb::tick_count endTime = ::tbb::tick_count::now();
double time = (endTime - startTime).seconds();
this->ReduceTime += time;
std::ostringstream out;
out << "Reduced " << range.size() << " key/value pairs in " << time << "s. "
<< "InRange: " << this->Ranges.InputBegin << " " << this->Ranges.InputEnd << " "
<< "OutRange: " << this->Ranges.OutputBegin << " " << this->Ranges.OutputEnd << "\n";
std::cerr << out.str();
#endif
}
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC_CONT
void join(const ReduceByKeyBody& rhs)
{
using KeysIteratorsType = vtkm::cont::ArrayPortalToIterators<KeysOutPortalType>;
using ValuesIteratorsType = vtkm::cont::ArrayPortalToIterators<ValuesOutPortalType>;
using KeysIteratorType = typename KeysIteratorsType::IteratorType;
using ValuesIteratorType = typename ValuesIteratorsType::IteratorType;
#ifdef _VTKM_DEBUG_TBB_RBK
::tbb::tick_count startTime = ::tbb::tick_count::now();
#endif
this->Ranges.AssertSane();
rhs.Ranges.AssertSane();
// Ensure that we're joining two consecutive subsets of the input:
VTKM_ASSERT(this->Ranges.IsNext(rhs.Ranges));
KeysIteratorsType keysIters(this->KeysOutPortal);
ValuesIteratorsType valuesIters(this->ValuesOutPortal);
KeysIteratorType keys = keysIters.GetBegin();
ValuesIteratorType values = valuesIters.GetBegin();
const vtkm::Id dstBegin = this->Ranges.OutputEnd;
const vtkm::Id lastDstIdx = this->Ranges.OutputEnd - 1;
vtkm::Id srcBegin = rhs.Ranges.OutputBegin;
const vtkm::Id srcEnd = rhs.Ranges.OutputEnd;
// Merge boundaries if needed:
if (keys[srcBegin] == keys[lastDstIdx])
{
values[lastDstIdx] = this->BinaryOperation(values[lastDstIdx], values[srcBegin]);
++srcBegin; // Don't copy the key/value we just reduced
}
// move data:
if (srcBegin != dstBegin && srcBegin != srcEnd)
{
// Sanity check:
VTKM_ASSERT(srcBegin < srcEnd);
// Not necessary for the copy call to be safe, but if the src range
// overlaps with the dst range there's a problem with the algorithm:
VTKM_ASSERT(dstBegin + (srcEnd - srcBegin) <= srcBegin);
std::copy(keys + srcBegin, keys + srcEnd, keys + dstBegin);
std::copy(values + srcBegin, values + srcEnd, values + dstBegin);
}
this->Ranges.InputEnd = rhs.Ranges.InputEnd;
this->Ranges.OutputEnd += srcEnd - srcBegin;
this->Ranges.AssertSane();
#ifdef _VTKM_DEBUG_TBB_RBK
::tbb::tick_count endTime = ::tbb::tick_count::now();
double time = (endTime - startTime).seconds();
this->JoinTime += rhs.JoinTime + time;
std::ostringstream out;
out << "Joined " << (srcEnd - srcBegin) << " rhs values into body in " << time << "s. "
<< "InRange: " << this->Ranges.InputBegin << " " << this->Ranges.InputEnd << " "
<< "OutRange: " << this->Ranges.OutputBegin << " " << this->Ranges.OutputEnd << "\n";
std::cerr << out.str();
#endif
}
};
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename KeysInPortalType,
typename ValuesInPortalType,
typename KeysOutPortalType,
typename ValuesOutPortalType,
typename BinaryOperationType>
VTKM_CONT vtkm::Id ReduceByKeyPortals(KeysInPortalType keysInPortal,
ValuesInPortalType valuesInPortal,
KeysOutPortalType keysOutPortal,
ValuesOutPortalType valuesOutPortal,
BinaryOperationType binaryOperation)
{
const vtkm::Id inputLength = keysInPortal.GetNumberOfValues();
VTKM_ASSERT(inputLength == valuesInPortal.GetNumberOfValues());
if (inputLength == 0)
{
return 0;
}
using ValueType = typename ValuesInPortalType::ValueType;
using WrappedBinaryOp = internal::WrappedBinaryOperator<ValueType, BinaryOperationType>;
WrappedBinaryOp wrappedBinaryOp(binaryOperation);
ReduceByKeyBody<KeysInPortalType,
ValuesInPortalType,
KeysOutPortalType,
ValuesOutPortalType,
WrappedBinaryOp>
body(keysInPortal, valuesInPortal, keysOutPortal, valuesOutPortal, wrappedBinaryOp);
::tbb::blocked_range<vtkm::Id> range(0, inputLength, TBB_GRAIN_SIZE);
#ifdef _VTKM_DEBUG_TBB_RBK
std::cerr << "\n\nTBB ReduceByKey:\n";
#endif
::tbb::parallel_reduce(range, body);
#ifdef _VTKM_DEBUG_TBB_RBK
std::cerr << "Total reduce time: " << body.ReduceTime << "s\n";
std::cerr << "Total join time: " << body.JoinTime << "s\n";
std::cerr << "\nend\n";
#endif
body.Ranges.AssertSane();
VTKM_ASSERT(body.Ranges.InputBegin == 0 && body.Ranges.InputEnd == inputLength &&
body.Ranges.OutputBegin == 0 && body.Ranges.OutputEnd <= inputLength);
return body.Ranges.OutputEnd;
}
#ifdef _VTKM_DEBUG_TBB_RBK
#undef _VTKM_DEBUG_TBB_RBK
#endif
template <class InputPortalType, class OutputPortalType, class BinaryOperationType>
struct ScanInclusiveBody
{

@ -40,6 +40,8 @@ set(headers
MarchingCubes.h
Mask.h
MaskPoints.h
NDEntropy.h
NDHistogram.h
PointAverage.h
PointElevation.h
PolicyBase.h
@ -74,6 +76,8 @@ set(header_template_sources
MarchingCubes.hxx
Mask.hxx
MaskPoints.hxx
NDEntropy.hxx
NDHistogram.hxx
PointAverage.hxx
PointElevation.hxx
SurfaceNormals.hxx

66
vtkm/filter/NDEntropy.h Normal file

@ -0,0 +1,66 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_filter_NDEntropy_h
#define vtk_m_filter_NDEntropy_h
#include <vtkm/filter/FilterDataSet.h>
namespace vtkm
{
namespace filter
{
/// \brief Clean a mesh to an unstructured grid
///
/// This filter takes a data set and essentially copies it into a new data set.
/// The newly constructed data set will have the same cells as the input and
/// the topology will be stored in a \c CellSetExplicit<>. The filter will also
/// optionally remove all unused points.
///
/// Note that the result of \c CleanGrid is not necessarily smaller than the
/// input. For example, "cleaning" a data set with a \c CellSetStructured
/// topology will actually result in a much larger data set.
///
/// \todo Add a feature to merge points that are coincident or within a
/// tolerance.
///
class NDEntropy : public vtkm::filter::FilterDataSet<NDEntropy>
{
public:
VTKM_CONT
NDEntropy();
VTKM_CONT
void AddFieldAndBin(const std::string& fieldName, vtkm::Id numOfBins);
template <typename Policy, typename Device>
VTKM_CONT vtkm::filter::Result DoExecute(const vtkm::cont::DataSet& inData,
vtkm::filter::PolicyBase<Policy> policy,
Device);
private:
std::vector<vtkm::Id> NumOfBins;
std::vector<std::string> FieldNames;
};
}
} // namespace vtkm::filter
#include <vtkm/filter/NDEntropy.hxx>
#endif //vtk_m_filter_NDEntropy_h

70
vtkm/filter/NDEntropy.hxx Normal file

@ -0,0 +1,70 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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 <vector>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/NDimsEntropy.h>
namespace vtkm
{
namespace filter
{
inline VTKM_CONT NDEntropy::NDEntropy()
{
}
void NDEntropy::AddFieldAndBin(const std::string& fieldName, vtkm::Id numOfBins)
{
this->FieldNames.push_back(fieldName);
this->NumOfBins.push_back(numOfBins);
}
template <typename Policy, typename Device>
inline VTKM_CONT vtkm::filter::Result NDEntropy::DoExecute(
const vtkm::cont::DataSet& inData,
vtkm::filter::PolicyBase<Policy> vtkmNotUsed(policy),
Device device)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
vtkm::worklet::NDimsEntropy ndEntropy;
ndEntropy.SetNumOfDataPoints(inData.GetField(0).GetData().GetNumberOfValues(), device);
// Add field one by one
// (By using AddFieldAndBin(), the length of FieldNames and NumOfBins must be the same)
for (size_t i = 0; i < FieldNames.size(); i++)
{
ndEntropy.AddField(inData.GetField(FieldNames[i]).GetData(), NumOfBins[i], device);
}
// Run worklet to calculate multi-variate entropy
vtkm::Float64 entropy = ndEntropy.Run(device);
vtkm::cont::DataSet outputData;
std::vector<vtkm::Float64> entropyHandle;
entropyHandle.push_back(entropy);
outputData.AddField(vtkm::cont::Field("Entropy", vtkm::cont::Field::ASSOC_POINTS, entropyHandle));
//return outputData;
return vtkm::filter::Result(outputData);
}
}
}

78
vtkm/filter/NDHistogram.h Normal file

@ -0,0 +1,78 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_filter_NDHistogram_h
#define vtk_m_filter_NDHistogram_h
#include <vtkm/filter/FilterDataSet.h>
namespace vtkm
{
namespace filter
{
/// \brief Clean a mesh to an unstructured grid
///
/// This filter takes a data set and essentially copies it into a new data set.
/// The newly constructed data set will have the same cells as the input and
/// the topology will be stored in a \c CellSetExplicit<>. The filter will also
/// optionally remove all unused points.
///
/// Note that the result of \c CleanGrid is not necessarily smaller than the
/// input. For example, "cleaning" a data set with a \c CellSetStructured
/// topology will actually result in a much larger data set.
///
/// \todo Add a feature to merge points that are coincident or within a
/// tolerance.
///
class NDHistogram : public vtkm::filter::FilterDataSet<NDHistogram>
{
public:
VTKM_CONT
NDHistogram();
VTKM_CONT
void AddFieldAndBin(const std::string& fieldName, vtkm::Id numOfBins);
template <typename Policy, typename Device>
VTKM_CONT vtkm::filter::Result DoExecute(const vtkm::cont::DataSet& inData,
vtkm::filter::PolicyBase<Policy> policy,
Device);
// This index is the field position in FieldNames
// (or the input _fieldName string vector of SetFields() Function)
VTKM_CONT
vtkm::Float64 GetBinDelta(size_t fieldIdx);
// This index is the field position in FieldNames
// (or the input _fieldName string vector of SetFields() Function)
VTKM_CONT
vtkm::Range GetDataRange(size_t fieldIdx);
private:
std::vector<vtkm::Id> NumOfBins;
std::vector<std::string> FieldNames;
std::vector<vtkm::Float64> BinDeltas;
std::vector<vtkm::Range> DataRanges; //Min Max of the field
};
}
} // namespace vtkm::filter
#include <vtkm/filter/NDHistogram.hxx>
#endif //vtk_m_filter_NDHistogram_h

@ -0,0 +1,91 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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 <vector>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/NDimsHistogram.h>
namespace vtkm
{
namespace filter
{
inline VTKM_CONT NDHistogram::NDHistogram()
{
}
void NDHistogram::AddFieldAndBin(const std::string& fieldName, vtkm::Id numOfBins)
{
this->FieldNames.push_back(fieldName);
this->NumOfBins.push_back(numOfBins);
}
vtkm::Float64 NDHistogram::GetBinDelta(size_t fieldIdx)
{
return BinDeltas[fieldIdx];
}
vtkm::Range NDHistogram::GetDataRange(size_t fieldIdx)
{
return DataRanges[fieldIdx];
}
template <typename Policy, typename Device>
inline VTKM_CONT vtkm::filter::Result NDHistogram::DoExecute(
const vtkm::cont::DataSet& inData,
vtkm::filter::PolicyBase<Policy> vtkmNotUsed(policy),
Device device)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
vtkm::worklet::NDimsHistogram ndHistogram;
// Set the number of data points
ndHistogram.SetNumOfDataPoints(inData.GetField(0).GetData().GetNumberOfValues(), device);
// Add field one by one
// (By using AddFieldAndBin(), the length of FieldNames and NumOfBins must be the same)
for (size_t i = 0; i < FieldNames.size(); i++)
{
vtkm::Range rangeField;
vtkm::Float64 deltaField;
ndHistogram.AddField(
inData.GetField(FieldNames[i]).GetData(), NumOfBins[i], rangeField, deltaField, device);
DataRanges.push_back(rangeField);
BinDeltas.push_back(deltaField);
}
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> binIds;
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
ndHistogram.Run(binIds, freqs, device);
vtkm::cont::DataSet outputData;
for (size_t i = 0; i < binIds.size(); i++)
{
outputData.AddField(
vtkm::cont::Field(FieldNames[i], vtkm::cont::Field::ASSOC_POINTS, binIds[i]));
}
outputData.AddField(vtkm::cont::Field("Frequency", vtkm::cont::Field::ASSOC_POINTS, freqs));
//return outputData;
return vtkm::filter::Result(outputData);
}
}
}

@ -36,6 +36,8 @@ set(unit_tests
UnitTestMaskFilter.cxx
UnitTestMaskPointsFilter.cxx
UnitTestMultiBlockFilters.cxx
UnitTestNDEntropyFilter.cxx
UnitTestNDHistogramFilter.cxx
UnitTestPointAverageFilter.cxx
UnitTestPointElevationFilter.cxx
UnitTestSurfaceNormalsFilter.cxx

@ -0,0 +1,212 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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/NDEntropy.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
// Make testing dataset with three fields(variables), each one has 1000 values
vtkm::cont::DataSet MakeTestDataSet()
{
vtkm::cont::DataSet dataSet;
const int xVerts = 20;
const int yVerts = 50;
const int nVerts = xVerts * yVerts;
vtkm::Float32 fieldA[nVerts] = {
8, 10, 9, 8, 14, 11, 12, 9, 19, 7, 8, 11, 7, 10, 11, 11, 11, 6, 8, 8, 7, 15, 9, 7,
8, 10, 9, 10, 10, 12, 7, 6, 14, 10, 14, 10, 7, 11, 13, 9, 13, 11, 10, 10, 12, 12, 7, 12,
10, 11, 12, 8, 13, 9, 5, 12, 11, 9, 5, 9, 12, 9, 6, 10, 11, 9, 9, 11, 9, 7, 7, 18,
16, 13, 12, 8, 10, 11, 9, 8, 17, 3, 15, 15, 9, 10, 10, 8, 10, 9, 7, 9, 8, 10, 13, 9,
7, 11, 7, 10, 13, 10, 11, 9, 10, 7, 10, 6, 12, 6, 9, 7, 6, 12, 12, 9, 12, 12, 11, 6,
1, 12, 8, 13, 14, 8, 8, 10, 7, 7, 6, 7, 5, 11, 6, 11, 13, 8, 13, 5, 9, 12, 7, 11,
10, 15, 11, 9, 7, 12, 15, 7, 8, 7, 12, 8, 21, 16, 13, 11, 10, 14, 12, 11, 12, 14, 7, 11,
7, 12, 16, 8, 10, 8, 9, 7, 8, 7, 13, 13, 11, 15, 7, 7, 6, 11, 7, 12, 12, 13, 14, 11,
13, 11, 11, 9, 15, 8, 6, 11, 12, 10, 11, 7, 6, 14, 11, 10, 12, 5, 8, 9, 11, 15, 11, 10,
17, 14, 9, 10, 10, 12, 11, 13, 13, 12, 11, 7, 8, 10, 7, 11, 10, 5, 8, 10, 13, 13, 12, 6,
10, 7, 13, 8, 11, 7, 10, 7, 8, 7, 14, 16, 9, 11, 8, 11, 9, 15, 11, 10, 10, 12, 7, 7,
11, 7, 5, 17, 9, 11, 11, 11, 10, 17, 10, 15, 7, 11, 12, 16, 9, 8, 11, 14, 9, 22, 8, 8,
8, 13, 12, 12, 1, 14, 15, 6, 15, 8, 11, 16, 14, 8, 6, 9, 8, 9, 9, 10, 8, 6, 13, 8,
6, 12, 11, 12, 13, 8, 6, 6, 5, 6, 10, 9, 11, 12, 14, 12, 10, 11, 10, 10, 8, 13, 8, 11,
7, 13, 13, 12, 12, 13, 15, 4, 9, 16, 7, 9, 8, 10, 6, 9, 11, 12, 6, 7, 14, 6, 4, 15,
5, 18, 9, 9, 11, 12, 9, 5, 6, 7, 15, 6, 11, 14, 8, 12, 6, 9, 5, 9, 14, 9, 12, 6,
9, 14, 11, 12, 12, 13, 15, 9, 8, 7, 13, 12, 7, 13, 6, 9, 10, 10, 10, 9, 11, 5, 9, 13,
16, 9, 10, 8, 9, 6, 13, 12, 8, 12, 9, 12, 17, 8, 11, 10, 8, 7, 11, 7, 13, 13, 10, 14,
11, 9, 6, 6, 14, 16, 5, 9, 13, 11, 12, 7, 4, 6, 9, 11, 11, 10, 12, 9, 7, 13, 8, 8,
12, 5, 10, 7, 11, 11, 10, 10, 14, 6, 8, 8, 3, 12, 16, 11, 11, 7, 6, 12, 11, 5, 9, 12,
9, 13, 7, 8, 9, 9, 12, 7, 9, 8, 12, 11, 6, 10, 6, 7, 6, 11, 10, 8, 9, 8, 4, 19,
12, 6, 10, 9, 6, 12, 9, 14, 7, 8, 11, 7, 7, 12, 13, 9, 13, 12, 8, 6, 10, 17, 19, 10,
10, 13, 5, 11, 8, 10, 8, 16, 12, 6, 6, 7, 10, 9, 12, 8, 5, 10, 7, 18, 9, 12, 10, 4,
9, 9, 15, 15, 6, 7, 7, 11, 12, 4, 8, 18, 5, 12, 12, 11, 10, 14, 9, 9, 10, 8, 10, 8,
10, 9, 9, 4, 10, 12, 5, 13, 6, 9, 7, 5, 12, 8, 11, 10, 9, 17, 9, 9, 8, 11, 18, 11,
10, 9, 4, 13, 10, 15, 5, 10, 9, 7, 7, 8, 10, 6, 6, 19, 10, 16, 7, 7, 9, 10, 10, 13,
10, 10, 14, 13, 12, 8, 7, 13, 12, 11, 13, 12, 9, 8, 6, 8, 10, 3, 8, 8, 12, 12, 13, 13,
10, 5, 10, 7, 13, 7, 9, 5, 13, 7, 10, 8, 13, 11, 17, 9, 6, 14, 10, 10, 13, 9, 15, 8,
15, 9, 12, 11, 12, 8, 3, 9, 8, 10, 12, 8, 14, 13, 12, 11, 12, 9, 18, 10, 13, 7, 4, 4,
11, 8, 3, 7, 9, 10, 12, 7, 11, 21, 9, 7, 8, 9, 10, 10, 11, 9, 15, 13, 21, 12, 8, 11,
9, 10, 11, 9, 17, 8, 9, 8, 14, 6, 13, 9, 8, 11, 12, 12, 12, 11, 6, 13, 7, 9, 11, 15,
17, 17, 11, 10, 7, 8, 11, 8, 6, 9, 13, 7, 9, 6, 5, 10, 7, 16, 16, 9, 7, 6, 14, 8,
13, 16, 7, 7, 10, 11, 6, 10, 9, 9, 8, 14, 11, 9, 11, 9, 10, 11, 9, 8, 14, 11, 7, 12,
11, 8, 9, 9, 10, 11, 11, 10, 9, 6, 6, 11, 16, 10, 7, 6, 6, 13, 18, 8, 12, 11, 14, 13,
8, 8, 10, 17, 17, 6, 6, 10, 18, 5, 8, 11, 6, 6, 14, 10, 9, 6, 11, 6, 13, 12, 10, 6,
9, 9, 9, 13, 7, 17, 10, 14, 10, 9, 10, 10, 11, 10, 11, 15, 13, 6, 12, 19, 10, 12, 12, 15,
13, 10, 10, 13, 11, 13, 13, 17, 6, 5, 6, 7, 6, 9, 13, 11, 8, 12, 9, 6, 10, 16, 11, 12,
5, 12, 14, 13, 13, 16, 11, 6, 12, 12, 15, 8, 7, 11, 8, 5, 10, 8, 9, 11, 9, 12, 10, 5,
12, 11, 9, 6, 14, 12, 10, 11, 9, 6, 7, 12, 8, 12, 8, 15, 9, 8, 7, 9, 3, 6, 14, 7,
8, 11, 9, 10, 12, 9, 10, 9, 8, 6, 12, 11, 6, 8, 9, 8, 15, 11, 7, 18, 12, 11, 10, 13,
11, 11, 10, 7, 9, 8, 8, 11, 11, 13, 6, 12, 13, 16, 11, 11, 5, 12, 14, 15, 9, 14, 15, 6,
8, 7, 6, 8, 9, 19, 7, 12, 11, 8, 14, 12, 10, 9, 3, 7
};
vtkm::Float32 fieldB[nVerts] = {
24, 19, 28, 19, 25, 28, 25, 22, 27, 26, 35, 26, 30, 28, 24, 23, 21, 31, 20, 11, 21, 22, 14, 25,
20, 24, 24, 21, 24, 29, 26, 21, 32, 29, 23, 28, 31, 25, 23, 30, 18, 24, 22, 25, 33, 24, 22, 23,
21, 17, 20, 28, 30, 18, 20, 32, 25, 24, 32, 15, 27, 24, 27, 19, 30, 27, 17, 24, 29, 23, 22, 19,
24, 19, 28, 24, 25, 24, 25, 30, 24, 31, 30, 27, 25, 25, 25, 15, 29, 23, 29, 29, 21, 25, 35, 24,
28, 10, 31, 23, 22, 22, 22, 33, 29, 27, 18, 27, 27, 24, 20, 20, 21, 29, 23, 31, 23, 23, 22, 23,
30, 27, 28, 31, 16, 29, 25, 19, 33, 28, 25, 24, 15, 27, 37, 29, 15, 19, 14, 19, 24, 23, 30, 29,
35, 22, 19, 26, 26, 14, 24, 30, 32, 23, 30, 29, 26, 27, 25, 23, 17, 26, 32, 29, 20, 17, 21, 23,
22, 20, 36, 12, 26, 23, 15, 29, 24, 22, 26, 33, 24, 23, 20, 26, 22, 17, 26, 26, 34, 22, 26, 17,
23, 18, 29, 27, 21, 29, 28, 29, 24, 25, 28, 19, 18, 21, 23, 23, 27, 25, 24, 25, 24, 25, 21, 25,
21, 27, 23, 20, 29, 15, 28, 30, 24, 27, 17, 23, 16, 21, 25, 17, 27, 28, 21, 13, 19, 27, 16, 30,
31, 25, 30, 17, 17, 25, 26, 22, 21, 17, 24, 17, 25, 22, 27, 14, 27, 24, 27, 25, 26, 31, 21, 23,
30, 30, 22, 19, 23, 22, 23, 25, 24, 25, 24, 28, 26, 30, 18, 25, 30, 37, 27, 34, 28, 34, 25, 10,
25, 22, 35, 30, 24, 32, 24, 34, 19, 29, 26, 16, 27, 17, 26, 23, 27, 25, 26, 21, 31, 21, 28, 15,
32, 24, 23, 23, 18, 15, 22, 25, 16, 25, 31, 26, 25, 28, 24, 26, 23, 25, 33, 20, 27, 28, 24, 29,
32, 20, 24, 20, 19, 32, 24, 6, 24, 21, 26, 18, 15, 30, 19, 26, 22, 30, 35, 23, 22, 30, 20, 22,
18, 30, 28, 25, 16, 25, 27, 30, 18, 24, 30, 28, 20, 19, 20, 28, 21, 24, 15, 33, 20, 18, 20, 36,
30, 26, 25, 18, 28, 27, 31, 31, 15, 26, 16, 22, 27, 14, 17, 27, 27, 22, 32, 30, 22, 34, 22, 25,
20, 22, 26, 29, 28, 33, 18, 23, 20, 20, 27, 24, 28, 21, 25, 27, 25, 19, 19, 25, 19, 32, 29, 27,
23, 21, 28, 33, 23, 23, 28, 26, 31, 19, 21, 29, 21, 27, 23, 32, 24, 26, 21, 28, 28, 24, 17, 31,
27, 21, 19, 32, 28, 23, 30, 23, 29, 15, 26, 26, 15, 20, 25, 26, 27, 31, 21, 23, 23, 33, 28, 19,
23, 22, 22, 25, 27, 17, 23, 17, 25, 28, 26, 30, 32, 31, 19, 25, 25, 19, 23, 29, 27, 23, 34, 22,
13, 21, 32, 10, 20, 33, 21, 17, 29, 31, 14, 24, 23, 19, 19, 22, 17, 26, 37, 26, 22, 26, 38, 29,
29, 27, 30, 20, 31, 14, 32, 32, 24, 23, 23, 18, 21, 31, 24, 20, 28, 15, 21, 25, 25, 20, 30, 25,
22, 21, 21, 25, 24, 25, 18, 23, 28, 30, 20, 27, 27, 19, 10, 32, 24, 20, 29, 26, 25, 20, 25, 29,
28, 24, 32, 26, 22, 19, 23, 27, 27, 29, 20, 25, 21, 30, 28, 31, 24, 19, 23, 19, 19, 18, 30, 18,
16, 24, 20, 20, 30, 25, 29, 25, 31, 21, 28, 31, 24, 26, 27, 21, 24, 23, 26, 18, 32, 26, 28, 26,
24, 26, 29, 30, 22, 20, 24, 28, 25, 29, 20, 21, 22, 15, 30, 27, 33, 26, 22, 32, 30, 31, 20, 19,
24, 26, 27, 31, 17, 17, 33, 27, 16, 27, 27, 22, 27, 19, 24, 21, 17, 24, 28, 23, 26, 24, 19, 26,
20, 24, 22, 19, 22, 21, 21, 28, 29, 39, 19, 16, 25, 29, 31, 22, 22, 29, 26, 22, 22, 22, 26, 23,
23, 23, 30, 25, 25, 25, 27, 29, 18, 33, 21, 12, 22, 29, 12, 20, 35, 22, 34, 28, 18, 29, 21, 20,
24, 33, 24, 26, 23, 34, 31, 25, 31, 22, 35, 21, 20, 29, 27, 22, 30, 22, 27, 23, 22, 32, 16, 19,
27, 22, 24, 27, 21, 33, 25, 25, 19, 28, 20, 27, 21, 25, 28, 20, 27, 22, 21, 20, 26, 30, 33, 23,
20, 24, 17, 23, 28, 35, 14, 23, 22, 28, 28, 26, 25, 18, 20, 28, 28, 22, 13, 24, 22, 20, 30, 26,
26, 18, 22, 20, 23, 24, 20, 27, 34, 28, 18, 24, 34, 33, 25, 33, 37, 21, 20, 31, 19, 23, 29, 22,
21, 24, 19, 27, 19, 32, 25, 23, 33, 26, 33, 27, 29, 30, 19, 22, 30, 19, 18, 24, 25, 17, 31, 19,
31, 26, 22, 23, 28, 28, 25, 24, 19, 19, 27, 28, 23, 21, 29, 26, 31, 22, 22, 25, 16, 29, 21, 22,
23, 25, 22, 21, 22, 19, 27, 26, 28, 30, 22, 21, 24, 22, 23, 26, 28, 22, 18, 25, 23, 27, 31, 19,
15, 29, 20, 19, 27, 25, 21, 29, 22, 24, 25, 17, 36, 29, 22, 22, 24, 28, 27, 22, 26, 31, 29, 31,
18, 25, 23, 16, 37, 27, 21, 31, 25, 24, 20, 23, 28, 33, 24, 21, 26, 20, 18, 31, 20, 24, 23, 19,
27, 17, 23, 23, 20, 26, 28, 23, 26, 31, 25, 31, 19, 32, 26, 18, 19, 29, 20, 21, 15, 25, 27, 29,
22, 22, 22, 26, 23, 22, 23, 29, 28, 20, 21, 22, 20, 22, 27, 25, 23, 32, 23, 20, 31, 20, 27, 26,
34, 20, 22, 36, 21, 29, 25, 20, 21, 22, 29, 29, 25, 22, 24, 22
};
vtkm::Float32 fieldC[nVerts] = {
3, 1, 4, 6, 5, 4, 8, 7, 2, 9, 2, 0, 0, 4, 3, 2, 5, 2, 3, 6, 3, 8, 3, 4,
3, 3, 2, 7, 2, 10, 9, 6, 1, 1, 4, 7, 3, 3, 1, 4, 4, 3, 9, 4, 4, 7, 3, 2,
4, 7, 3, 3, 2, 10, 1, 6, 2, 2, 3, 8, 3, 3, 6, 9, 4, 1, 4, 3, 16, 7, 0, 1,
8, 7, 13, 3, 5, 0, 3, 8, 10, 3, 5, 5, 1, 5, 2, 1, 3, 2, 5, 3, 4, 3, 3, 3,
3, 1, 13, 2, 3, 1, 2, 7, 3, 4, 1, 2, 5, 4, 4, 4, 2, 6, 3, 2, 7, 8, 1, 3,
4, 1, 2, 0, 1, 6, 1, 8, 8, 1, 1, 4, 2, 1, 4, 3, 5, 4, 6, 4, 2, 3, 8, 8,
3, 3, 3, 4, 5, 8, 8, 16, 7, 12, 4, 3, 14, 8, 3, 12, 5, 0, 5, 3, 5, 2, 9, 2,
9, 4, 1, 0, 0, 4, 4, 6, 3, 4, 11, 2, 4, 7, 4, 2, 1, 9, 4, 3, 2, 5, 1, 5,
3, 8, 2, 8, 1, 8, 0, 4, 1, 3, 2, 1, 2, 3, 2, 1, 8, 5, 4, 1, 9, 9, 1, 3,
5, 0, 1, 6, 10, 8, 3, 12, 3, 4, 4, 7, 1, 3, 6, 4, 4, 6, 1, 4, 7, 5, 6, 11,
6, 5, 2, 7, 2, 5, 3, 5, 6, 3, 6, 2, 1, 10, 8, 3, 7, 0, 2, 6, 9, 3, 11, 3,
2, 5, 1, 4, 6, 10, 9, 1, 4, 3, 7, 12, 3, 10, 0, 2, 11, 2, 1, 0, 4, 1, 2, 16,
5, 17, 7, 8, 2, 10, 10, 3, 1, 3, 2, 2, 4, 8, 4, 3, 2, 4, 4, 6, 8, 6, 2, 3,
2, 4, 2, 4, 7, 10, 5, 3, 5, 2, 4, 6, 9, 3, 1, 1, 1, 1, 4, 2, 2, 7, 4, 9,
2, 3, 5, 6, 2, 5, 1, 6, 5, 7, 8, 3, 7, 2, 2, 8, 6, 2, 10, 2, 1, 4, 5, 1,
1, 1, 5, 6, 1, 1, 4, 5, 4, 2, 4, 3, 2, 7, 19, 4, 7, 2, 7, 5, 2, 5, 3, 8,
4, 6, 7, 2, 0, 0, 2, 12, 6, 2, 2, 3, 5, 9, 4, 9, 2, 2, 7, 8, 3, 3, 10, 6,
3, 2, 1, 6, 2, 4, 6, 3, 5, 8, 2, 3, 6, 14, 0, 3, 6, 5, 2, 7, 0, 3, 8, 5,
3, 2, 2, 5, 1, 3, 12, 11, 16, 2, 1, 3, 7, 3, 1, 6, 4, 3, 12, 5, 1, 3, 1, 4,
9, 1, 3, 3, 4, 4, 6, 7, 7, 5, 2, 4, 2, 3, 2, 2, 6, 4, 2, 2, 3, 5, 1, 4,
9, 1, 0, 7, 6, 4, 3, 3, 7, 3, 3, 6, 2, 7, 9, 3, 1, 16, 5, 4, 3, 6, 3, 2,
5, 2, 2, 4, 3, 1, 3, 3, 6, 3, 5, 9, 1, 10, 1, 7, 2, 2, 6, 7, 3, 5, 3, 7,
2, 2, 2, 2, 6, 4, 3, 2, 5, 5, 3, 15, 4, 2, 7, 7, 4, 3, 3, 5, 1, 2, 9, 0,
5, 7, 12, 2, 4, 8, 5, 7, 8, 3, 2, 2, 18, 1, 7, 2, 2, 1, 3, 3, 3, 7, 1, 9,
8, 4, 3, 7, 6, 4, 5, 2, 0, 5, 1, 5, 10, 4, 2, 8, 2, 2, 0, 5, 6, 4, 5, 0,
1, 5, 11, 3, 3, 4, 4, 2, 3, 5, 1, 6, 5, 7, 2, 2, 5, 7, 4, 8, 4, 1, 1, 7,
2, 3, 9, 6, 13, 1, 5, 4, 6, 2, 4, 11, 2, 5, 5, 1, 4, 1, 4, 7, 1, 5, 8, 3,
1, 10, 9, 13, 1, 7, 2, 9, 4, 3, 3, 10, 12, 2, 0, 4, 6, 5, 5, 1, 4, 7, 2, 12,
7, 6, 5, 0, 6, 4, 4, 12, 1, 3, 10, 1, 9, 2, 2, 2, 1, 5, 5, 6, 9, 6, 4, 1,
11, 6, 9, 3, 2, 7, 1, 7, 4, 3, 0, 3, 1, 12, 17, 2, 1, 6, 4, 4, 2, 1, 5, 5,
3, 2, 2, 4, 6, 5, 4, 6, 11, 3, 12, 6, 3, 6, 3, 0, 6, 3, 7, 4, 8, 5, 14, 5,
1, 9, 4, 6, 5, 3, 9, 3, 1, 1, 0, 3, 7, 3, 5, 1, 6, 2, 2, 6, 2, 12, 1, 0,
6, 3, 3, 5, 4, 7, 2, 2, 15, 7, 3, 10, 4, 2, 6, 3, 4, 8, 3, 1, 5, 5, 5, 4,
3, 7, 3, 4, 5, 5, 2, 4, 2, 5, 1, 12, 5, 6, 3, 2, 8, 5, 2, 3, 11, 11, 6, 5,
0, 3, 3, 9, 4, 2, 11, 1, 5, 3, 5, 6, 3, 6, 4, 2, 4, 10, 11, 3, 3, 4, 1, 1,
1, 3, 5, 5, 1, 1, 4, 1, 5, 1, 6, 8, 6, 4, 6, 7, 6, 3, 5, 3, 6, 6, 6, 4,
0, 6, 3, 1, 2, 4, 2, 6, 1, 1, 1, 2, 2, 4, 7, 2, 6, 2, 5, 7, 6, 4, 6, 3,
1, 4, 5, 1, 4, 6, 2, 3, 0, 6, 11, 2, 9, 2, 6, 4, 5, 6, 2, 19, 2, 10, 4, 2,
3, 3, 11, 7, 3, 3, 1, 5, 3, 6, 4, 3, 0, 6, 6, 6, 4, 2, 5, 2, 2, 2, 6, 10,
4, 9, 3, 7, 7, 0, 6, 8, 5, 2, 3, 2, 3, 3, 3, 1, 6, 1, 8, 2, 5, 3, 6, 11,
5, 7, 2, 6, 7, 3, 4, 1, 0, 5, 8, 3, 2, 9, 3, 1, 2, 3, 3, 9, 5, 6, 5, 1,
4, 5, 6, 7, 6, 1, 5, 1, 6, 6, 2, 6, 7, 2, 4, 6
};
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vtkm::Id3(xVerts, yVerts, 1));
dataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", coordinates));
// Set point scalars
dataSet.AddField(vtkm::cont::Field("fieldA", vtkm::cont::Field::ASSOC_POINTS, fieldA, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldB", vtkm::cont::Field::ASSOC_POINTS, fieldB, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldC", vtkm::cont::Field::ASSOC_POINTS, fieldC, nVerts));
return dataSet;
}
void RunTest()
{
vtkm::cont::DataSet ds = MakeTestDataSet();
vtkm::filter::NDEntropy ndEntropyFilter;
ndEntropyFilter.AddFieldAndBin("fieldA", 10);
ndEntropyFilter.AddFieldAndBin("fieldB", 10);
ndEntropyFilter.AddFieldAndBin("fieldC", 10);
//Entropy is just a single number
//It is stored in "Entropy" (type vtkm::Float64) field in the return dataset
//The field has only one element, the entropy
vtkm::filter::Result resultData = ndEntropyFilter.Execute(ds);
vtkm::cont::DataSet& outputData = resultData.GetDataSet();
vtkm::cont::ArrayHandle<vtkm::Float64> entropyHandle;
outputData.GetField("Entropy").GetData().CopyTo(entropyHandle);
vtkm::Float64 e = entropyHandle.GetPortalControl().Get(0);
VTKM_TEST_ASSERT(fabs(e - 7.457857) < 0.001,
"N-Dimentional entropy filter calculation is incorrect");
}
} // anonymous namespace
int UnitTestNDEntropyFilter(int, char* [])
{
return vtkm::cont::testing::Testing::Run(RunTest);
}

@ -0,0 +1,129 @@
//============================================================================
// 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 2016 Sandia Corporation.
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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/NDHistogram.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
vtkm::cont::DataSet MakeTestDataSet()
{
vtkm::cont::DataSet dataSet;
const int nVerts = 100;
vtkm::Float32 fieldA[nVerts] = { 8, 10, 9, 8, 14, 11, 12, 9, 19, 7, 8, 11, 7, 10, 11,
11, 11, 6, 8, 8, 7, 15, 9, 7, 8, 10, 9, 10, 10, 12,
7, 6, 14, 10, 14, 10, 7, 11, 13, 9, 13, 11, 10, 10, 12,
12, 7, 12, 10, 11, 12, 8, 13, 9, 5, 12, 11, 9, 5, 9,
12, 9, 6, 10, 11, 9, 9, 11, 9, 7, 7, 18, 16, 13, 12,
8, 10, 11, 9, 8, 17, 3, 15, 15, 9, 10, 10, 8, 10, 9,
7, 9, 8, 10, 13, 9, 7, 11, 7, 10 };
vtkm::Float32 fieldB[nVerts] = { 24, 19, 28, 19, 25, 28, 25, 22, 27, 26, 35, 26, 30, 28, 24,
23, 21, 31, 20, 11, 21, 22, 14, 25, 20, 24, 24, 21, 24, 29,
26, 21, 32, 29, 23, 28, 31, 25, 23, 30, 18, 24, 22, 25, 33,
24, 22, 23, 21, 17, 20, 28, 30, 18, 20, 32, 25, 24, 32, 15,
27, 24, 27, 19, 30, 27, 17, 24, 29, 23, 22, 19, 24, 19, 28,
24, 25, 24, 25, 30, 24, 31, 30, 27, 25, 25, 25, 15, 29, 23,
29, 29, 21, 25, 35, 24, 28, 10, 31, 23 };
vtkm::Float32 fieldC[nVerts] = {
3, 1, 4, 6, 5, 4, 8, 7, 2, 9, 2, 0, 0, 4, 3, 2, 5, 2, 3, 6, 3, 8, 3, 4, 3,
3, 2, 7, 2, 10, 9, 6, 1, 1, 4, 7, 3, 3, 1, 4, 4, 3, 9, 4, 4, 7, 3, 2, 4, 7,
3, 3, 2, 10, 1, 6, 2, 2, 3, 8, 3, 3, 6, 9, 4, 1, 4, 3, 16, 7, 0, 1, 8, 7, 13,
3, 5, 0, 3, 8, 10, 3, 5, 5, 1, 5, 2, 1, 3, 2, 5, 3, 4, 3, 3, 3, 3, 1, 13, 2
};
// Set point scalars
dataSet.AddField(vtkm::cont::Field("fieldA", vtkm::cont::Field::ASSOC_POINTS, fieldA, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldB", vtkm::cont::Field::ASSOC_POINTS, fieldB, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldC", vtkm::cont::Field::ASSOC_POINTS, fieldC, nVerts));
return dataSet;
}
void RunTest()
{
vtkm::cont::DataSet ds = MakeTestDataSet();
vtkm::filter::NDHistogram ndHistFilter;
ndHistFilter.AddFieldAndBin("fieldA", 4);
ndHistFilter.AddFieldAndBin("fieldB", 4);
ndHistFilter.AddFieldAndBin("fieldC", 4);
//The return data set contains fieldNames.size() + 1 fields
//The first "fieldNames.size()"" fields are the binId arrays for inputs field
//And their order and field names are the same as the order and name in fieldNames
//The name of last fields in the dataset is "Frequency"
//This field contains the all freqncys of the N-Dims histogram
//The result histogram is stored in sparse representation
//(Do not store and return zero frequency bins)
//All fields in return dataset must have the same length
//So, e.g. (FieldA[i], FieldB[i], FieldC[i], Frequency[i] ) is a bin in the histogram
//First three numbers are binID for FieldA, FieldB, FieldC
//Frequency[i] is frequency for this bin
vtkm::filter::Result resultData = ndHistFilter.Execute(ds);
vtkm::cont::DataSet& outputData = resultData.GetDataSet();
// Ground truth ND histogram
vtkm::Id gtNonSparseBins = 33;
vtkm::Id gtIdx0[33] = { 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3 };
vtkm::Id gtIdx1[33] = { 1, 1, 2, 3, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3,
0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 1, 1, 2, 2, 2, 3 };
vtkm::Id gtIdx2[33] = { 0, 1, 1, 0, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 3,
0, 0, 1, 0, 1, 2, 3, 0, 1, 2, 0, 2, 0, 1, 2, 1 };
vtkm::Id gtFreq[33] = { 1, 1, 1, 3, 2, 1, 1, 6, 6, 3, 17, 8, 2, 6, 2, 1, 2,
1, 1, 4, 11, 4, 1, 1, 3, 3, 1, 1, 1, 1, 1, 2, 1 };
// Check result
vtkm::Id nonSparseBins = outputData.GetField(0).GetData().GetNumberOfValues();
VTKM_TEST_ASSERT(nonSparseBins == gtNonSparseBins, "Incorrect ND-histogram Filter results");
vtkm::cont::ArrayHandle<vtkm::Id> binId0;
outputData.GetField("fieldA").GetData().CopyTo(binId0);
vtkm::cont::ArrayHandle<vtkm::Id> binId1;
outputData.GetField("fieldB").GetData().CopyTo(binId1);
vtkm::cont::ArrayHandle<vtkm::Id> binId2;
outputData.GetField("fieldC").GetData().CopyTo(binId2);
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
outputData.GetField("Frequency").GetData().CopyTo(freqs);
for (int i = 0; i < nonSparseBins; i++)
{
vtkm::Id idx0 = binId0.GetPortalControl().Get(i);
vtkm::Id idx1 = binId1.GetPortalControl().Get(i);
vtkm::Id idx2 = binId2.GetPortalControl().Get(i);
vtkm::Id f = freqs.GetPortalControl().Get(i);
VTKM_TEST_ASSERT(idx0 == gtIdx0[i] && idx1 == gtIdx1[i] && idx2 == gtIdx2[i] && f == gtFreq[i],
"Incorrect ND-histogram Filter results");
}
}
} // anonymous namespace
int UnitTestNDHistogramFilter(int, char* [])
{
return vtkm::cont::testing::Testing::Run(RunTest);
}

@ -1562,9 +1562,9 @@ namespace detail
}
template<class... Ts>
#ifdef BRIGAND_COMP_MSVC_2013
using is_set = typename detail::is_set_impl<range<int, 0, sizeof...(Ts)>, detail::msvc_quali_ref<Ts>...>::type;
using is_set = typename detail::is_set_impl<range<int, 0, static_cast<int>(sizeof...(Ts))>, detail::msvc_quali_ref<Ts>...>::type;
#else
using is_set = typename detail::is_set_impl<range<int, 0, sizeof...(Ts)>, Ts...>::type;
using is_set = typename detail::is_set_impl<range<int, 0, static_cast<int>(sizeof...(Ts))>, Ts...>::type;
#endif
}
#include <type_traits>

@ -46,6 +46,9 @@ set(headers
MarchingCubesDataTables.h
Mask.h
MaskPoints.h
NDimsEntropy.h
NDimsHistogram.h
NDimsHistMarginalization.h
ParticleAdvection.h
PointAverage.h
PointElevation.h

103
vtkm/worklet/NDimsEntropy.h Normal file

@ -0,0 +1,103 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_NDimsEntropy_h
#define vtk_m_worklet_NDimsEntropy_h
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/NDimsHistogram.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/histogram/ComputeNDEntropy.h>
#include <vtkm/cont/Field.h>
namespace vtkm
{
namespace worklet
{
class NDimsEntropy
{
public:
template <typename DeviceAdapter>
void SetNumOfDataPoints(vtkm::Id _numDataPoints, DeviceAdapter device)
{
NumDataPoints = _numDataPoints;
NdHistogram.SetNumOfDataPoints(_numDataPoints, device);
}
// Add a field and the bin for this field
// Return: rangeOfRange is min max value of this array
// binDelta is delta of a bin
template <typename HandleType, typename DeviceAdapter>
void AddField(const HandleType& fieldArray, vtkm::Id numberOfBins, DeviceAdapter device)
{
vtkm::Range range;
vtkm::Float64 delta;
NdHistogram.AddField(fieldArray, numberOfBins, range, delta, device);
}
// Execute the entropy computation filter given
// fields and number of bins of each fields
// Returns:
// Entropy (log2) of the field of the data
template <typename DeviceAdapter>
vtkm::Float64 Run(DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> binIds;
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
NdHistogram.Run(binIds, freqs, device);
///// calculate sum of frequency of the histogram /////
vtkm::Id initFreqSumValue = 0;
vtkm::Id freqSum = DeviceAlgorithms::Reduce(freqs, initFreqSumValue, vtkm::Sum());
///// calculate information content of each bin using self-define worklet /////
vtkm::cont::ArrayHandle<vtkm::Float64> informationContent;
vtkm::worklet::histogram::SetBinInformationContent binWorklet(
static_cast<vtkm::Float64>(freqSum));
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::SetBinInformationContent,
DeviceAdapter>
setBinInformationContentDispatcher(binWorklet);
setBinInformationContentDispatcher.Invoke(freqs, informationContent);
///// calculate entropy by summing up information conetent of all bins /////
vtkm::Float64 initEntropyValue = 0;
vtkm::Float64 entropy =
DeviceAlgorithms::Reduce(informationContent, initEntropyValue, vtkm::Sum());
return entropy;
}
private:
vtkm::worklet::NDimsHistogram NdHistogram;
vtkm::Id NumDataPoints;
};
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_NDimsEntropy_h

@ -0,0 +1,219 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_NDimsHistMarginalization_h
#define vtk_m_worklet_NDimsHistMarginalization_h
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/worklet/histogram/MarginalizeNDHistogram.h>
#include <vtkm/worklet/histogram/MarginalizeNDHistogram.h>
#include <vtkm/worklet/histogram/MarginalizeNDHistogram.h>
#include <vtkm/cont/Field.h>
namespace vtkm
{
namespace worklet
{
class NDimsHistMarginalization
{
public:
// Execute the histogram (conditional) marginalization,
// given the multi-variable histogram(binId, freqIn)
// , marginalVariable and marginal condition
// Input arguments:
// binId, freqsIn: input ND-histogram in the fashion of sparse representation
// (definition of binId and frqIn please refer to NDimsHistogram.h),
// (binId.size() is the number of variables)
// numberOfBins: number of bins of each variable (length of numberOfBins must be the same as binId.size() )
// marginalVariables: length is the same as number of variables.
// 1 indicates marginal variable, otherwise 0.
// conditionFunc: The Condition function for non-marginal variable.
// This func takes two arguments (vtkm::Id var, vtkm::Id binId) and return bool
// var is index of variable and binId is bin index in the varaiable var
// return true indicates considering this bin into final marginal histogram
// more details can refer to example in UnitTestNDimsHistMarginalization.cxx
// marginalBinId, marginalFreqs: return marginalized histogram in the fashion of sparse representation
// the definition is the same as (binId and freqsIn)
template <typename BinaryCompare, typename DeviceAdapter>
void Run(const std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& binId,
vtkm::cont::ArrayHandle<vtkm::Id>& freqsIn,
vtkm::cont::ArrayHandle<vtkm::Id>& numberOfBins,
vtkm::cont::ArrayHandle<bool>& marginalVariables,
BinaryCompare conditionFunc,
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& marginalBinId,
vtkm::cont::ArrayHandle<vtkm::Id>& marginalFreqs,
DeviceAdapter vtkmNotUsed(device))
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
//total variables
vtkm::Id numOfVariable = static_cast<vtkm::Id>(binId.size());
const vtkm::Id numberOfValues = freqsIn.GetPortalConstControl().GetNumberOfValues();
vtkm::cont::ArrayHandleConstant<vtkm::Id> constant0Array(0, numberOfValues);
vtkm::cont::ArrayHandle<vtkm::Id> bin1DIndex;
DeviceAlgorithms::Copy(constant0Array, bin1DIndex);
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
DeviceAlgorithms::Copy(freqsIn, freqs);
vtkm::Id numMarginalVariables = 0; //count num of marginal variables
for (vtkm::Id i = 0; i < numOfVariable; i++)
{
if (marginalVariables.GetPortalConstControl().Get(i) == true)
{
// Worklet to calculate 1D index for marginal variables
numMarginalVariables++;
const vtkm::Id nFieldBins = numberOfBins.GetPortalControl().Get(i);
vtkm::worklet::histogram::To1DIndex binWorklet(nFieldBins);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::To1DIndex, DeviceAdapter>
to1DIndexDispatcher(binWorklet);
size_t vecIndex = static_cast<size_t>(i);
to1DIndexDispatcher.Invoke(binId[vecIndex], bin1DIndex, bin1DIndex);
}
else
{ //non-marginal variable
// Worklet to set the frequency of entities which does not meet the condition
// to 0 on non-marginal variables
vtkm::worklet::histogram::ConditionalFreq<BinaryCompare> conditionalFreqWorklet{
conditionFunc
};
conditionalFreqWorklet.setVar(i);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::ConditionalFreq<BinaryCompare>,
DeviceAdapter>
cfDispatcher(conditionalFreqWorklet);
size_t vecIndex = static_cast<size_t>(i);
cfDispatcher.Invoke(binId[vecIndex], freqs, freqs);
}
}
// Sort the freq array for counting by key(1DIndex)
DeviceAlgorithms::SortByKey(bin1DIndex, freqs);
// Add frequency within same 1d index bin (this get a nonSparse representation)
vtkm::cont::ArrayHandle<vtkm::Id> nonSparseMarginalFreqs;
DeviceAlgorithms::ReduceByKey(
bin1DIndex, freqs, bin1DIndex, nonSparseMarginalFreqs, vtkm::Add());
// Convert to sparse representation(remove all zero freqncy entities)
vtkm::cont::ArrayHandle<vtkm::Id> sparseMarginal1DBinId;
DeviceAlgorithms::CopyIf(bin1DIndex, nonSparseMarginalFreqs, sparseMarginal1DBinId);
DeviceAlgorithms::CopyIf(nonSparseMarginalFreqs, nonSparseMarginalFreqs, marginalFreqs);
//convert back to multi variate binId
marginalBinId.resize(static_cast<size_t>(numMarginalVariables));
vtkm::Id marginalVarIdx = numMarginalVariables - 1;
for (vtkm::Id i = numOfVariable - 1; i >= 0; i--)
{
if (marginalVariables.GetPortalConstControl().Get(i) == true)
{
const vtkm::Id nFieldBins = numberOfBins.GetPortalControl().Get(i);
vtkm::worklet::histogram::ConvertHistBinToND binWorklet(nFieldBins);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::ConvertHistBinToND,
DeviceAdapter>
ConvertHistBinToNDDispatcher(binWorklet);
size_t vecIndex = static_cast<size_t>(marginalVarIdx);
ConvertHistBinToNDDispatcher.Invoke(
sparseMarginal1DBinId, sparseMarginal1DBinId, marginalBinId[vecIndex]);
marginalVarIdx--;
}
}
} //Run()
// Execute the histogram marginalization WITHOUT CONDITION,
// Please refer to the other Run() functions for the definition of input arguments.
template <typename DeviceAdapter>
void Run(const std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& binId,
vtkm::cont::ArrayHandle<vtkm::Id>& freqsIn,
vtkm::cont::ArrayHandle<vtkm::Id>& numberOfBins,
vtkm::cont::ArrayHandle<bool>& marginalVariables,
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& marginalBinId,
vtkm::cont::ArrayHandle<vtkm::Id>& marginalFreqs,
DeviceAdapter vtkmNotUsed(device))
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
//total variables
vtkm::Id numOfVariable = static_cast<vtkm::Id>(binId.size());
const vtkm::Id numberOfValues = freqsIn.GetPortalConstControl().GetNumberOfValues();
vtkm::cont::ArrayHandleConstant<vtkm::Id> constant0Array(0, numberOfValues);
vtkm::cont::ArrayHandle<vtkm::Id> bin1DIndex;
DeviceAlgorithms::Copy(constant0Array, bin1DIndex);
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
DeviceAlgorithms::Copy(freqsIn, freqs);
vtkm::Id numMarginalVariables = 0; //count num of marginal varaibles
for (vtkm::Id i = 0; i < numOfVariable; i++)
{
if (marginalVariables.GetPortalConstControl().Get(i) == true)
{
// Worklet to calculate 1D index for marginal variables
numMarginalVariables++;
const vtkm::Id nFieldBins = numberOfBins.GetPortalControl().Get(i);
vtkm::worklet::histogram::To1DIndex binWorklet(nFieldBins);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::To1DIndex, DeviceAdapter>
to1DIndexDispatcher(binWorklet);
size_t vecIndex = static_cast<size_t>(i);
to1DIndexDispatcher.Invoke(binId[vecIndex], bin1DIndex, bin1DIndex);
}
}
// Sort the freq array for counting by key (1DIndex)
DeviceAlgorithms::SortByKey(bin1DIndex, freqs);
// Add frequency within same 1d index bin
DeviceAlgorithms::ReduceByKey(bin1DIndex, freqs, bin1DIndex, marginalFreqs, vtkm::Add());
//convert back to multi variate binId
marginalBinId.resize(static_cast<size_t>(numMarginalVariables));
vtkm::Id marginalVarIdx = numMarginalVariables - 1;
for (vtkm::Id i = numOfVariable - 1; i >= 0; i--)
{
if (marginalVariables.GetPortalConstControl().Get(i) == true)
{
const vtkm::Id nFieldBins = numberOfBins.GetPortalControl().Get(i);
vtkm::worklet::histogram::ConvertHistBinToND binWorklet(nFieldBins);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::ConvertHistBinToND,
DeviceAdapter>
ConvertHistBinToNDDispatcher(binWorklet);
size_t vecIndex = static_cast<size_t>(marginalVarIdx);
ConvertHistBinToNDDispatcher.Invoke(bin1DIndex, bin1DIndex, marginalBinId[vecIndex]);
marginalVarIdx--;
}
}
} //Run()
};
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_NDimsHistMarginalization_h

@ -0,0 +1,127 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_NDimsHistogram_h
#define vtk_m_worklet_NDimsHistogram_h
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ErrorBadValue.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/worklet/histogram/ComputeNDHistogram.h>
#include <vtkm/cont/Field.h>
namespace vtkm
{
namespace worklet
{
class NDimsHistogram
{
public:
template <typename DeviceAdapter>
void SetNumOfDataPoints(vtkm::Id _numDataPoints, DeviceAdapter vtkmNotUsed(device))
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
NumDataPoints = _numDataPoints;
// Initialize bin1DIndex array
vtkm::cont::ArrayHandleConstant<vtkm::Id> constant0Array(0, NumDataPoints);
DeviceAlgorithms::Copy(constant0Array, Bin1DIndex);
}
// Add a field and the bin number for this field
// Return: rangeOfRange is min max value of this array
// binDelta is delta of a bin
template <typename HandleType, typename DeviceAdapter>
void AddField(const HandleType& fieldArray,
vtkm::Id numberOfBins,
vtkm::Range& rangeOfValues,
vtkm::Float64& binDelta,
DeviceAdapter vtkmNotUsed(device))
{
NumberOfBins.push_back(numberOfBins);
if (fieldArray.GetNumberOfValues() != NumDataPoints)
{
throw vtkm::cont::ErrorBadValue("Array lengths does not match");
}
else
{
CastAndCall(fieldArray.ResetTypeList(vtkm::TypeListTagScalarAll()),
vtkm::worklet::histogram::ComputeBins<DeviceAdapter>(
Bin1DIndex, numberOfBins, rangeOfValues, binDelta));
}
}
// Execute N-Dim histogram worklet to get N-Dims histogram from input fields
// Input arguments:
// binId: returned bin id of NDims-histogram, binId has n arrays, if length of fieldName is n
// freqs: returned frequncy(count) array
// Note: the ND-histogram is returned in the fashion of sparse representation.
// (no zero freqncy in freqs array)
// the length of all arrays in binId and freqs array must be the same
// if the length of fieldNames is n (compute a n-dimensional hisotgram)
// freqs[i] is the freqncy of the bin with bin Ids{ binId[0][i], binId[1][i], ... binId[n-1][i] }
template <typename DeviceAdapter>
void Run(std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& binId,
vtkm::cont::ArrayHandle<vtkm::Id>& freqs,
DeviceAdapter vtkmNotUsed(device))
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
binId.resize(NumberOfBins.size());
// Sort the resulting bin(1D) array for counting
DeviceAlgorithms::Sort(Bin1DIndex);
// Count frequency of each bin
vtkm::cont::ArrayHandleConstant<vtkm::Id> constArray(1, NumDataPoints);
DeviceAlgorithms::ReduceByKey(Bin1DIndex, constArray, Bin1DIndex, freqs, vtkm::Add());
//convert back to multi variate binId
for (vtkm::Id i = static_cast<vtkm::Id>(NumberOfBins.size()) - 1; i >= 0; i--)
{
const vtkm::Id nFieldBins = NumberOfBins[static_cast<size_t>(i)];
vtkm::worklet::histogram::ConvertHistBinToND binWorklet(nFieldBins);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::ConvertHistBinToND, DeviceAdapter>
ConvertHistBinToNDDispatcher(binWorklet);
size_t vectorId = static_cast<size_t>(i);
ConvertHistBinToNDDispatcher.Invoke(Bin1DIndex, Bin1DIndex, binId[vectorId]);
}
}
private:
std::vector<vtkm::Id> NumberOfBins;
vtkm::cont::ArrayHandle<vtkm::Id> Bin1DIndex;
vtkm::Id NumDataPoints;
};
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_NDimsHistogram_h

@ -80,8 +80,12 @@ public:
vtkm::cont::ArrayHandleConstant<vtkm::Id> init(0, numSeeds);
stepsTaken.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(init, stepsTaken);
worklet.Run(it, pts, nSteps, status, stepsTaken);
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(statusOK, status);
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
return res;
@ -94,7 +98,7 @@ public:
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsTaken,
const vtkm::cont::ArrayHandle<vtkm::Id>& stepsAlreadyTaken,
const vtkm::Id& nSteps,
const DeviceAdapter&)
{
@ -102,9 +106,18 @@ public:
FieldType,
DeviceAdapter>
worklet;
vtkm::cont::ArrayHandle<vtkm::Id> status;
worklet.Run(it, pts, nSteps, status, stepsTaken);
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken, status;
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
//Allocate status and steps arrays.
stepsTaken.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(stepsAlreadyTaken, stepsTaken);
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(statusOK, status);
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
return res;

@ -0,0 +1,27 @@
##============================================================================
## 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 Sandia Corporation.
## Copyright 2014 UT-Battelle, LLC.
## Copyright 2014 Los Alamos National Security.
##
## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
## 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.
##============================================================================
set(headers
ComputeNDEntropy.h
ComputeNDHistogram.h
MarginalizeNDHistogram.h
)
vtkm_declare_headers(${headers})

@ -0,0 +1,61 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_ComputeNDEntropy_h
#define vtk_m_worklet_ComputeNDEntropy_h
#include <vtkm/worklet/DispatcherMapField.h>
namespace vtkm
{
namespace worklet
{
namespace histogram
{
// For each bin, calculate its information content (log2)
class SetBinInformationContent : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> freq, FieldOut<> informationContent);
typedef void ExecutionSignature(_1, _2);
vtkm::Float64 FreqSum;
VTKM_CONT
SetBinInformationContent(vtkm::Float64 _freqSum)
: FreqSum(_freqSum)
{
}
template <typename FreqType>
VTKM_EXEC void operator()(const FreqType& freq, vtkm::Float64& informationContent) const
{
vtkm::Float64 p = ((vtkm::Float64)freq) / FreqSum;
if (p > 0)
informationContent = -1 * p * vtkm::Log2(p);
else
informationContent = 0;
}
};
}
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_ComputeNDEntropy_h

@ -0,0 +1,157 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_ComputeNDHistogram_h
#define vtk_m_worklet_ComputeNDHistogram_h
#include <vtkm/worklet/DispatcherMapField.h>
namespace vtkm
{
namespace worklet
{
namespace histogram
{
// GCC creates false positive warnings for signed/unsigned char* operations.
// This occurs because the values are implicitly casted up to int's for the
// operation, and than casted back down to char's when return.
// This causes a false positive warning, even when the values is within
// the value types range
#if defined(VTKM_GCC)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
#endif // gcc
template <typename T>
T compute_delta(T fieldMinValue, T fieldMaxValue, vtkm::Id num)
{
typedef vtkm::VecTraits<T> VecType;
const T fieldRange = fieldMaxValue - fieldMinValue;
return fieldRange / static_cast<typename VecType::ComponentType>(num);
}
#if defined(VTKM_GCC)
#pragma GCC diagnostic pop
#endif // gcc
// For each value set the bin it should be in
template <typename FieldType>
class SetHistogramBin : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> value, FieldIn<> binIndexIn, FieldOut<> binIndexOut);
typedef void ExecutionSignature(_1, _2, _3);
typedef _1 InputDomain;
vtkm::Id numberOfBins;
vtkm::Float64 minValue;
vtkm::Float64 delta;
VTKM_CONT
SetHistogramBin(vtkm::Id numberOfBins0, vtkm::Float64 minValue0, vtkm::Float64 delta0)
: numberOfBins(numberOfBins0)
, minValue(minValue0)
, delta(delta0)
{
}
VTKM_EXEC
void operator()(const FieldType& value, const vtkm::Id& binIndexIn, vtkm::Id& binIndexOut) const
{
vtkm::Id localBinIdx = static_cast<vtkm::Id>((value - minValue) / delta);
if (localBinIdx < 0)
localBinIdx = 0;
else if (localBinIdx >= numberOfBins)
localBinIdx = numberOfBins - 1;
binIndexOut = binIndexIn * numberOfBins + localBinIdx;
}
};
template <typename DeviceAdapter>
class ComputeBins
{
public:
VTKM_CONT
ComputeBins(vtkm::cont::ArrayHandle<vtkm::Id>& _bin1DIdx,
vtkm::Id& _numOfBins,
vtkm::Range& _minMax,
vtkm::Float64& _binDelta)
: Bin1DIdx(_bin1DIdx)
, NumOfBins(_numOfBins)
, MinMax(_minMax)
, BinDelta(_binDelta)
{
}
template <typename T, typename Storage>
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<T, Storage>& field) const
{
typedef vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
const vtkm::Vec<T, 2> initValue(field.GetPortalConstControl().Get(0));
vtkm::Vec<T, 2> minMax = Algorithm::Reduce(field, initValue, vtkm::MinAndMax<T>());
MinMax.Min = static_cast<vtkm::Float64>(minMax[0]);
MinMax.Max = static_cast<vtkm::Float64>(minMax[1]);
BinDelta = compute_delta(MinMax.Min, MinMax.Max, NumOfBins);
SetHistogramBin<T> binWorklet(NumOfBins, MinMax.Min, BinDelta);
vtkm::worklet::DispatcherMapField<vtkm::worklet::histogram::SetHistogramBin<T>, DeviceAdapter>
setHistogramBinDispatcher(binWorklet);
setHistogramBinDispatcher.Invoke(field, Bin1DIdx, Bin1DIdx);
}
private:
vtkm::cont::ArrayHandle<vtkm::Id>& Bin1DIdx;
vtkm::Id& NumOfBins;
vtkm::Range& MinMax;
vtkm::Float64& BinDelta;
};
// Convert N-dims bin index into 1D index
class ConvertHistBinToND : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> bin1DIndexIn,
FieldOut<> bin1DIndexOut,
FieldOut<> oneVariableIndexOut);
typedef void ExecutionSignature(_1, _2, _3);
typedef _1 InputDomain;
vtkm::Id numberOfBins;
VTKM_CONT
ConvertHistBinToND(vtkm::Id numberOfBins0)
: numberOfBins(numberOfBins0)
{
}
VTKM_EXEC
void operator()(const vtkm::Id& bin1DIndexIn,
vtkm::Id& bin1DIndexOut,
vtkm::Id& oneVariableIndexOut) const
{
oneVariableIndexOut = bin1DIndexIn % numberOfBins;
bin1DIndexOut = (bin1DIndexIn - oneVariableIndexOut) / numberOfBins;
}
};
}
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_ComputeNDHistogram_h

@ -0,0 +1,88 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//============================================================================
#ifndef vtk_m_worklet_MarginalizeNDHistogram_h
#define vtk_m_worklet_MarginalizeNDHistogram_h
#include <vtkm/worklet/DispatcherMapField.h>
namespace vtkm
{
namespace worklet
{
namespace histogram
{
// Set freq of the entity, which does not meet the condition, to 0
template <class BinaryCompare>
class ConditionalFreq : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
ConditionalFreq(BinaryCompare _bop)
: bop(_bop)
{
}
VTKM_CONT
void setVar(vtkm::Id _var) { var = _var; }
BinaryCompare bop;
vtkm::Id var;
typedef void ControlSignature(FieldIn<>, FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_EXEC
void operator()(const vtkm::Id& binIdIn, const vtkm::Id& freqIn, vtkm::Id& freqOut) const
{
if (bop(var, binIdIn))
freqOut = freqIn;
else
freqOut = 0;
}
};
// Convert multiple indices to 1D index
class To1DIndex : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> bin, FieldIn<> binIndexIn, FieldOut<> binIndexOut);
typedef void ExecutionSignature(_1, _2, _3);
typedef _1 InputDomain;
vtkm::Id numberOfBins;
VTKM_CONT
To1DIndex(vtkm::Id numberOfBins0)
: numberOfBins(numberOfBins0)
{
}
VTKM_EXEC
void operator()(const vtkm::Id& bin, const vtkm::Id& binIndexIn, vtkm::Id& binIndexOut) const
{
binIndexOut = binIndexIn * numberOfBins + bin;
}
};
}
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_MarginalizeNDHistogram_h

@ -57,7 +57,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -101,7 +101,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -210,7 +210,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -378,7 +378,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);

@ -79,46 +79,59 @@ public:
}
VTKM_EXEC
FieldType GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
FieldType stepLength,
vtkm::Vec<FieldType, 3>& velocity) const
ParticleStatus GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
FieldType& stepLength,
vtkm::Vec<FieldType, 3>& velocity) const
{
this->CheckStep(inpos, stepLength, velocity);
ParticleStatus status = this->CheckStep(inpos, stepLength, velocity);
if (status != ParticleStatus::STATUS_OK)
{
stepLength += this->Tolerance;
return status;
}
FieldType magnitude = vtkm::Magnitude(velocity);
vtkm::Vec<FieldType, 3> dir = velocity / magnitude;
vtkm::Vec<FieldType, 3> dirBounds;
this->Evaluator.GetBoundary(dir, dirBounds);
/*Add a fraction just push the particle beyond the bounds*/
FieldType hx = (std::abs(dirBounds[0] - inpos[0]) + this->Tolerance) / std::abs(velocity[0]);
FieldType hy = (std::abs(dirBounds[1] - inpos[1]) + this->Tolerance) / std::abs(velocity[1]);
FieldType hz = (std::abs(dirBounds[2] - inpos[2]) + this->Tolerance) / std::abs(velocity[2]);
return vtkm::Min(hx, vtkm::Min(hy, hz));
FieldType hx = (vtkm::Abs(dirBounds[0] - inpos[0]) + this->Tolerance) / vtkm::Abs(velocity[0]);
FieldType hy = (vtkm::Abs(dirBounds[1] - inpos[1]) + this->Tolerance) / vtkm::Abs(velocity[1]);
FieldType hz = (vtkm::Abs(dirBounds[2] - inpos[2]) + this->Tolerance) / vtkm::Abs(velocity[2]);
stepLength = vtkm::Min(hx, vtkm::Min(hy, hz));
return status;
}
VTKM_EXEC
ParticleStatus PushOutOfDomain(vtkm::Vec<FieldType, 3> inpos,
ParticleStatus PushOutOfDomain(const vtkm::Vec<FieldType, 3>& inpos,
vtkm::Id numSteps,
vtkm::Vec<FieldType, 3>& outpos) const
{
ParticleStatus status;
outpos = inpos;
numSteps = (numSteps == 0) ? 1 : numSteps;
FieldType totalTime = static_cast<FieldType>(numSteps) * this->StepLength;
FieldType timeFraction = totalTime * this->Tolerance;
FieldType stepLength = this->StepLength / 2;
vtkm::Vec<FieldType, 3> velocity;
vtkm::Vec<FieldType, 3> velocity, currentVelocity;
this->CheckStep(inpos, 0.0f, currentVelocity);
if (this->ShortStepsSupported)
{
do
{
ParticleStatus status = this->CheckStep(inpos, stepLength, velocity);
status = this->CheckStep(inpos, stepLength, velocity);
if (status == ParticleStatus::STATUS_OK)
{
inpos = inpos + stepLength * velocity;
outpos = outpos + stepLength * velocity;
currentVelocity = velocity;
}
stepLength = stepLength / 2;
} while (stepLength > timeFraction);
}
stepLength = GetEscapeStepLength(inpos, stepLength, velocity);
outpos = inpos + stepLength * velocity;
status = GetEscapeStepLength(inpos, stepLength, velocity);
if (status != ParticleStatus::STATUS_OK)
outpos = outpos + stepLength * currentVelocity;
else
outpos = outpos + stepLength * velocity;
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
}

@ -114,20 +114,16 @@ private:
void run(vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& statusArray,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& stepsTaken)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
ParticleWorkletDispatchType;
typedef vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag> ParticleType;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleType = vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
//Allocate status and steps arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> ok(ParticleStatus::STATUS_OK, numSeeds);
statusArray.Allocate(numSeeds);
DeviceAlgorithm::Copy(ok, statusArray);
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleType particles(seedArray, stepsTaken, statusArray, maxSteps);
//Invoke particle advection worklet
ParticleAdvectWorkletType particleWorklet(integrator);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
particleWorkletDispatch.Invoke(idxArray, particles);

@ -130,9 +130,17 @@ public:
SetBit(idx, TERMINATED);
}
VTKM_EXEC
void SetExitedSpatialBoundary(const vtkm::Id& idx) { SetBit(idx, EXITED_SPATIAL_BOUNDARY); }
void SetExitedSpatialBoundary(const vtkm::Id& idx)
{
ClearBit(idx, STATUS_OK);
SetBit(idx, EXITED_SPATIAL_BOUNDARY);
}
VTKM_EXEC
void SetExitedTemporalBoundary(const vtkm::Id& idx) { SetBit(idx, EXITED_TEMPORAL_BOUNDARY); }
void SetExitedTemporalBoundary(const vtkm::Id& idx)
{
ClearBit(idx, STATUS_OK);
SetBit(idx, EXITED_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
void SetError(const vtkm::Id& idx)
{

@ -38,6 +38,9 @@ set(unit_tests
UnitTestMarchingCubes.cxx
UnitTestMask.cxx
UnitTestMaskPoints.cxx
UnitTestNDimsEntropy.cxx
UnitTestNDimsHistogram.cxx
UnitTestNDimsHistMarginalization.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx

@ -0,0 +1,213 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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/worklet/NDimsEntropy.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
// Make testing dataset with three fields(variables), each one has 1000 values
vtkm::cont::DataSet MakeTestDataSet()
{
vtkm::cont::DataSet dataSet;
const int xVerts = 20;
const int yVerts = 50;
const int nVerts = xVerts * yVerts;
vtkm::Float32 fieldA[nVerts] = {
8, 10, 9, 8, 14, 11, 12, 9, 19, 7, 8, 11, 7, 10, 11, 11, 11, 6, 8, 8, 7, 15, 9, 7,
8, 10, 9, 10, 10, 12, 7, 6, 14, 10, 14, 10, 7, 11, 13, 9, 13, 11, 10, 10, 12, 12, 7, 12,
10, 11, 12, 8, 13, 9, 5, 12, 11, 9, 5, 9, 12, 9, 6, 10, 11, 9, 9, 11, 9, 7, 7, 18,
16, 13, 12, 8, 10, 11, 9, 8, 17, 3, 15, 15, 9, 10, 10, 8, 10, 9, 7, 9, 8, 10, 13, 9,
7, 11, 7, 10, 13, 10, 11, 9, 10, 7, 10, 6, 12, 6, 9, 7, 6, 12, 12, 9, 12, 12, 11, 6,
1, 12, 8, 13, 14, 8, 8, 10, 7, 7, 6, 7, 5, 11, 6, 11, 13, 8, 13, 5, 9, 12, 7, 11,
10, 15, 11, 9, 7, 12, 15, 7, 8, 7, 12, 8, 21, 16, 13, 11, 10, 14, 12, 11, 12, 14, 7, 11,
7, 12, 16, 8, 10, 8, 9, 7, 8, 7, 13, 13, 11, 15, 7, 7, 6, 11, 7, 12, 12, 13, 14, 11,
13, 11, 11, 9, 15, 8, 6, 11, 12, 10, 11, 7, 6, 14, 11, 10, 12, 5, 8, 9, 11, 15, 11, 10,
17, 14, 9, 10, 10, 12, 11, 13, 13, 12, 11, 7, 8, 10, 7, 11, 10, 5, 8, 10, 13, 13, 12, 6,
10, 7, 13, 8, 11, 7, 10, 7, 8, 7, 14, 16, 9, 11, 8, 11, 9, 15, 11, 10, 10, 12, 7, 7,
11, 7, 5, 17, 9, 11, 11, 11, 10, 17, 10, 15, 7, 11, 12, 16, 9, 8, 11, 14, 9, 22, 8, 8,
8, 13, 12, 12, 1, 14, 15, 6, 15, 8, 11, 16, 14, 8, 6, 9, 8, 9, 9, 10, 8, 6, 13, 8,
6, 12, 11, 12, 13, 8, 6, 6, 5, 6, 10, 9, 11, 12, 14, 12, 10, 11, 10, 10, 8, 13, 8, 11,
7, 13, 13, 12, 12, 13, 15, 4, 9, 16, 7, 9, 8, 10, 6, 9, 11, 12, 6, 7, 14, 6, 4, 15,
5, 18, 9, 9, 11, 12, 9, 5, 6, 7, 15, 6, 11, 14, 8, 12, 6, 9, 5, 9, 14, 9, 12, 6,
9, 14, 11, 12, 12, 13, 15, 9, 8, 7, 13, 12, 7, 13, 6, 9, 10, 10, 10, 9, 11, 5, 9, 13,
16, 9, 10, 8, 9, 6, 13, 12, 8, 12, 9, 12, 17, 8, 11, 10, 8, 7, 11, 7, 13, 13, 10, 14,
11, 9, 6, 6, 14, 16, 5, 9, 13, 11, 12, 7, 4, 6, 9, 11, 11, 10, 12, 9, 7, 13, 8, 8,
12, 5, 10, 7, 11, 11, 10, 10, 14, 6, 8, 8, 3, 12, 16, 11, 11, 7, 6, 12, 11, 5, 9, 12,
9, 13, 7, 8, 9, 9, 12, 7, 9, 8, 12, 11, 6, 10, 6, 7, 6, 11, 10, 8, 9, 8, 4, 19,
12, 6, 10, 9, 6, 12, 9, 14, 7, 8, 11, 7, 7, 12, 13, 9, 13, 12, 8, 6, 10, 17, 19, 10,
10, 13, 5, 11, 8, 10, 8, 16, 12, 6, 6, 7, 10, 9, 12, 8, 5, 10, 7, 18, 9, 12, 10, 4,
9, 9, 15, 15, 6, 7, 7, 11, 12, 4, 8, 18, 5, 12, 12, 11, 10, 14, 9, 9, 10, 8, 10, 8,
10, 9, 9, 4, 10, 12, 5, 13, 6, 9, 7, 5, 12, 8, 11, 10, 9, 17, 9, 9, 8, 11, 18, 11,
10, 9, 4, 13, 10, 15, 5, 10, 9, 7, 7, 8, 10, 6, 6, 19, 10, 16, 7, 7, 9, 10, 10, 13,
10, 10, 14, 13, 12, 8, 7, 13, 12, 11, 13, 12, 9, 8, 6, 8, 10, 3, 8, 8, 12, 12, 13, 13,
10, 5, 10, 7, 13, 7, 9, 5, 13, 7, 10, 8, 13, 11, 17, 9, 6, 14, 10, 10, 13, 9, 15, 8,
15, 9, 12, 11, 12, 8, 3, 9, 8, 10, 12, 8, 14, 13, 12, 11, 12, 9, 18, 10, 13, 7, 4, 4,
11, 8, 3, 7, 9, 10, 12, 7, 11, 21, 9, 7, 8, 9, 10, 10, 11, 9, 15, 13, 21, 12, 8, 11,
9, 10, 11, 9, 17, 8, 9, 8, 14, 6, 13, 9, 8, 11, 12, 12, 12, 11, 6, 13, 7, 9, 11, 15,
17, 17, 11, 10, 7, 8, 11, 8, 6, 9, 13, 7, 9, 6, 5, 10, 7, 16, 16, 9, 7, 6, 14, 8,
13, 16, 7, 7, 10, 11, 6, 10, 9, 9, 8, 14, 11, 9, 11, 9, 10, 11, 9, 8, 14, 11, 7, 12,
11, 8, 9, 9, 10, 11, 11, 10, 9, 6, 6, 11, 16, 10, 7, 6, 6, 13, 18, 8, 12, 11, 14, 13,
8, 8, 10, 17, 17, 6, 6, 10, 18, 5, 8, 11, 6, 6, 14, 10, 9, 6, 11, 6, 13, 12, 10, 6,
9, 9, 9, 13, 7, 17, 10, 14, 10, 9, 10, 10, 11, 10, 11, 15, 13, 6, 12, 19, 10, 12, 12, 15,
13, 10, 10, 13, 11, 13, 13, 17, 6, 5, 6, 7, 6, 9, 13, 11, 8, 12, 9, 6, 10, 16, 11, 12,
5, 12, 14, 13, 13, 16, 11, 6, 12, 12, 15, 8, 7, 11, 8, 5, 10, 8, 9, 11, 9, 12, 10, 5,
12, 11, 9, 6, 14, 12, 10, 11, 9, 6, 7, 12, 8, 12, 8, 15, 9, 8, 7, 9, 3, 6, 14, 7,
8, 11, 9, 10, 12, 9, 10, 9, 8, 6, 12, 11, 6, 8, 9, 8, 15, 11, 7, 18, 12, 11, 10, 13,
11, 11, 10, 7, 9, 8, 8, 11, 11, 13, 6, 12, 13, 16, 11, 11, 5, 12, 14, 15, 9, 14, 15, 6,
8, 7, 6, 8, 9, 19, 7, 12, 11, 8, 14, 12, 10, 9, 3, 7
};
vtkm::Float32 fieldB[nVerts] = {
24, 19, 28, 19, 25, 28, 25, 22, 27, 26, 35, 26, 30, 28, 24, 23, 21, 31, 20, 11, 21, 22, 14, 25,
20, 24, 24, 21, 24, 29, 26, 21, 32, 29, 23, 28, 31, 25, 23, 30, 18, 24, 22, 25, 33, 24, 22, 23,
21, 17, 20, 28, 30, 18, 20, 32, 25, 24, 32, 15, 27, 24, 27, 19, 30, 27, 17, 24, 29, 23, 22, 19,
24, 19, 28, 24, 25, 24, 25, 30, 24, 31, 30, 27, 25, 25, 25, 15, 29, 23, 29, 29, 21, 25, 35, 24,
28, 10, 31, 23, 22, 22, 22, 33, 29, 27, 18, 27, 27, 24, 20, 20, 21, 29, 23, 31, 23, 23, 22, 23,
30, 27, 28, 31, 16, 29, 25, 19, 33, 28, 25, 24, 15, 27, 37, 29, 15, 19, 14, 19, 24, 23, 30, 29,
35, 22, 19, 26, 26, 14, 24, 30, 32, 23, 30, 29, 26, 27, 25, 23, 17, 26, 32, 29, 20, 17, 21, 23,
22, 20, 36, 12, 26, 23, 15, 29, 24, 22, 26, 33, 24, 23, 20, 26, 22, 17, 26, 26, 34, 22, 26, 17,
23, 18, 29, 27, 21, 29, 28, 29, 24, 25, 28, 19, 18, 21, 23, 23, 27, 25, 24, 25, 24, 25, 21, 25,
21, 27, 23, 20, 29, 15, 28, 30, 24, 27, 17, 23, 16, 21, 25, 17, 27, 28, 21, 13, 19, 27, 16, 30,
31, 25, 30, 17, 17, 25, 26, 22, 21, 17, 24, 17, 25, 22, 27, 14, 27, 24, 27, 25, 26, 31, 21, 23,
30, 30, 22, 19, 23, 22, 23, 25, 24, 25, 24, 28, 26, 30, 18, 25, 30, 37, 27, 34, 28, 34, 25, 10,
25, 22, 35, 30, 24, 32, 24, 34, 19, 29, 26, 16, 27, 17, 26, 23, 27, 25, 26, 21, 31, 21, 28, 15,
32, 24, 23, 23, 18, 15, 22, 25, 16, 25, 31, 26, 25, 28, 24, 26, 23, 25, 33, 20, 27, 28, 24, 29,
32, 20, 24, 20, 19, 32, 24, 6, 24, 21, 26, 18, 15, 30, 19, 26, 22, 30, 35, 23, 22, 30, 20, 22,
18, 30, 28, 25, 16, 25, 27, 30, 18, 24, 30, 28, 20, 19, 20, 28, 21, 24, 15, 33, 20, 18, 20, 36,
30, 26, 25, 18, 28, 27, 31, 31, 15, 26, 16, 22, 27, 14, 17, 27, 27, 22, 32, 30, 22, 34, 22, 25,
20, 22, 26, 29, 28, 33, 18, 23, 20, 20, 27, 24, 28, 21, 25, 27, 25, 19, 19, 25, 19, 32, 29, 27,
23, 21, 28, 33, 23, 23, 28, 26, 31, 19, 21, 29, 21, 27, 23, 32, 24, 26, 21, 28, 28, 24, 17, 31,
27, 21, 19, 32, 28, 23, 30, 23, 29, 15, 26, 26, 15, 20, 25, 26, 27, 31, 21, 23, 23, 33, 28, 19,
23, 22, 22, 25, 27, 17, 23, 17, 25, 28, 26, 30, 32, 31, 19, 25, 25, 19, 23, 29, 27, 23, 34, 22,
13, 21, 32, 10, 20, 33, 21, 17, 29, 31, 14, 24, 23, 19, 19, 22, 17, 26, 37, 26, 22, 26, 38, 29,
29, 27, 30, 20, 31, 14, 32, 32, 24, 23, 23, 18, 21, 31, 24, 20, 28, 15, 21, 25, 25, 20, 30, 25,
22, 21, 21, 25, 24, 25, 18, 23, 28, 30, 20, 27, 27, 19, 10, 32, 24, 20, 29, 26, 25, 20, 25, 29,
28, 24, 32, 26, 22, 19, 23, 27, 27, 29, 20, 25, 21, 30, 28, 31, 24, 19, 23, 19, 19, 18, 30, 18,
16, 24, 20, 20, 30, 25, 29, 25, 31, 21, 28, 31, 24, 26, 27, 21, 24, 23, 26, 18, 32, 26, 28, 26,
24, 26, 29, 30, 22, 20, 24, 28, 25, 29, 20, 21, 22, 15, 30, 27, 33, 26, 22, 32, 30, 31, 20, 19,
24, 26, 27, 31, 17, 17, 33, 27, 16, 27, 27, 22, 27, 19, 24, 21, 17, 24, 28, 23, 26, 24, 19, 26,
20, 24, 22, 19, 22, 21, 21, 28, 29, 39, 19, 16, 25, 29, 31, 22, 22, 29, 26, 22, 22, 22, 26, 23,
23, 23, 30, 25, 25, 25, 27, 29, 18, 33, 21, 12, 22, 29, 12, 20, 35, 22, 34, 28, 18, 29, 21, 20,
24, 33, 24, 26, 23, 34, 31, 25, 31, 22, 35, 21, 20, 29, 27, 22, 30, 22, 27, 23, 22, 32, 16, 19,
27, 22, 24, 27, 21, 33, 25, 25, 19, 28, 20, 27, 21, 25, 28, 20, 27, 22, 21, 20, 26, 30, 33, 23,
20, 24, 17, 23, 28, 35, 14, 23, 22, 28, 28, 26, 25, 18, 20, 28, 28, 22, 13, 24, 22, 20, 30, 26,
26, 18, 22, 20, 23, 24, 20, 27, 34, 28, 18, 24, 34, 33, 25, 33, 37, 21, 20, 31, 19, 23, 29, 22,
21, 24, 19, 27, 19, 32, 25, 23, 33, 26, 33, 27, 29, 30, 19, 22, 30, 19, 18, 24, 25, 17, 31, 19,
31, 26, 22, 23, 28, 28, 25, 24, 19, 19, 27, 28, 23, 21, 29, 26, 31, 22, 22, 25, 16, 29, 21, 22,
23, 25, 22, 21, 22, 19, 27, 26, 28, 30, 22, 21, 24, 22, 23, 26, 28, 22, 18, 25, 23, 27, 31, 19,
15, 29, 20, 19, 27, 25, 21, 29, 22, 24, 25, 17, 36, 29, 22, 22, 24, 28, 27, 22, 26, 31, 29, 31,
18, 25, 23, 16, 37, 27, 21, 31, 25, 24, 20, 23, 28, 33, 24, 21, 26, 20, 18, 31, 20, 24, 23, 19,
27, 17, 23, 23, 20, 26, 28, 23, 26, 31, 25, 31, 19, 32, 26, 18, 19, 29, 20, 21, 15, 25, 27, 29,
22, 22, 22, 26, 23, 22, 23, 29, 28, 20, 21, 22, 20, 22, 27, 25, 23, 32, 23, 20, 31, 20, 27, 26,
34, 20, 22, 36, 21, 29, 25, 20, 21, 22, 29, 29, 25, 22, 24, 22
};
vtkm::Float32 fieldC[nVerts] = {
3, 1, 4, 6, 5, 4, 8, 7, 2, 9, 2, 0, 0, 4, 3, 2, 5, 2, 3, 6, 3, 8, 3, 4,
3, 3, 2, 7, 2, 10, 9, 6, 1, 1, 4, 7, 3, 3, 1, 4, 4, 3, 9, 4, 4, 7, 3, 2,
4, 7, 3, 3, 2, 10, 1, 6, 2, 2, 3, 8, 3, 3, 6, 9, 4, 1, 4, 3, 16, 7, 0, 1,
8, 7, 13, 3, 5, 0, 3, 8, 10, 3, 5, 5, 1, 5, 2, 1, 3, 2, 5, 3, 4, 3, 3, 3,
3, 1, 13, 2, 3, 1, 2, 7, 3, 4, 1, 2, 5, 4, 4, 4, 2, 6, 3, 2, 7, 8, 1, 3,
4, 1, 2, 0, 1, 6, 1, 8, 8, 1, 1, 4, 2, 1, 4, 3, 5, 4, 6, 4, 2, 3, 8, 8,
3, 3, 3, 4, 5, 8, 8, 16, 7, 12, 4, 3, 14, 8, 3, 12, 5, 0, 5, 3, 5, 2, 9, 2,
9, 4, 1, 0, 0, 4, 4, 6, 3, 4, 11, 2, 4, 7, 4, 2, 1, 9, 4, 3, 2, 5, 1, 5,
3, 8, 2, 8, 1, 8, 0, 4, 1, 3, 2, 1, 2, 3, 2, 1, 8, 5, 4, 1, 9, 9, 1, 3,
5, 0, 1, 6, 10, 8, 3, 12, 3, 4, 4, 7, 1, 3, 6, 4, 4, 6, 1, 4, 7, 5, 6, 11,
6, 5, 2, 7, 2, 5, 3, 5, 6, 3, 6, 2, 1, 10, 8, 3, 7, 0, 2, 6, 9, 3, 11, 3,
2, 5, 1, 4, 6, 10, 9, 1, 4, 3, 7, 12, 3, 10, 0, 2, 11, 2, 1, 0, 4, 1, 2, 16,
5, 17, 7, 8, 2, 10, 10, 3, 1, 3, 2, 2, 4, 8, 4, 3, 2, 4, 4, 6, 8, 6, 2, 3,
2, 4, 2, 4, 7, 10, 5, 3, 5, 2, 4, 6, 9, 3, 1, 1, 1, 1, 4, 2, 2, 7, 4, 9,
2, 3, 5, 6, 2, 5, 1, 6, 5, 7, 8, 3, 7, 2, 2, 8, 6, 2, 10, 2, 1, 4, 5, 1,
1, 1, 5, 6, 1, 1, 4, 5, 4, 2, 4, 3, 2, 7, 19, 4, 7, 2, 7, 5, 2, 5, 3, 8,
4, 6, 7, 2, 0, 0, 2, 12, 6, 2, 2, 3, 5, 9, 4, 9, 2, 2, 7, 8, 3, 3, 10, 6,
3, 2, 1, 6, 2, 4, 6, 3, 5, 8, 2, 3, 6, 14, 0, 3, 6, 5, 2, 7, 0, 3, 8, 5,
3, 2, 2, 5, 1, 3, 12, 11, 16, 2, 1, 3, 7, 3, 1, 6, 4, 3, 12, 5, 1, 3, 1, 4,
9, 1, 3, 3, 4, 4, 6, 7, 7, 5, 2, 4, 2, 3, 2, 2, 6, 4, 2, 2, 3, 5, 1, 4,
9, 1, 0, 7, 6, 4, 3, 3, 7, 3, 3, 6, 2, 7, 9, 3, 1, 16, 5, 4, 3, 6, 3, 2,
5, 2, 2, 4, 3, 1, 3, 3, 6, 3, 5, 9, 1, 10, 1, 7, 2, 2, 6, 7, 3, 5, 3, 7,
2, 2, 2, 2, 6, 4, 3, 2, 5, 5, 3, 15, 4, 2, 7, 7, 4, 3, 3, 5, 1, 2, 9, 0,
5, 7, 12, 2, 4, 8, 5, 7, 8, 3, 2, 2, 18, 1, 7, 2, 2, 1, 3, 3, 3, 7, 1, 9,
8, 4, 3, 7, 6, 4, 5, 2, 0, 5, 1, 5, 10, 4, 2, 8, 2, 2, 0, 5, 6, 4, 5, 0,
1, 5, 11, 3, 3, 4, 4, 2, 3, 5, 1, 6, 5, 7, 2, 2, 5, 7, 4, 8, 4, 1, 1, 7,
2, 3, 9, 6, 13, 1, 5, 4, 6, 2, 4, 11, 2, 5, 5, 1, 4, 1, 4, 7, 1, 5, 8, 3,
1, 10, 9, 13, 1, 7, 2, 9, 4, 3, 3, 10, 12, 2, 0, 4, 6, 5, 5, 1, 4, 7, 2, 12,
7, 6, 5, 0, 6, 4, 4, 12, 1, 3, 10, 1, 9, 2, 2, 2, 1, 5, 5, 6, 9, 6, 4, 1,
11, 6, 9, 3, 2, 7, 1, 7, 4, 3, 0, 3, 1, 12, 17, 2, 1, 6, 4, 4, 2, 1, 5, 5,
3, 2, 2, 4, 6, 5, 4, 6, 11, 3, 12, 6, 3, 6, 3, 0, 6, 3, 7, 4, 8, 5, 14, 5,
1, 9, 4, 6, 5, 3, 9, 3, 1, 1, 0, 3, 7, 3, 5, 1, 6, 2, 2, 6, 2, 12, 1, 0,
6, 3, 3, 5, 4, 7, 2, 2, 15, 7, 3, 10, 4, 2, 6, 3, 4, 8, 3, 1, 5, 5, 5, 4,
3, 7, 3, 4, 5, 5, 2, 4, 2, 5, 1, 12, 5, 6, 3, 2, 8, 5, 2, 3, 11, 11, 6, 5,
0, 3, 3, 9, 4, 2, 11, 1, 5, 3, 5, 6, 3, 6, 4, 2, 4, 10, 11, 3, 3, 4, 1, 1,
1, 3, 5, 5, 1, 1, 4, 1, 5, 1, 6, 8, 6, 4, 6, 7, 6, 3, 5, 3, 6, 6, 6, 4,
0, 6, 3, 1, 2, 4, 2, 6, 1, 1, 1, 2, 2, 4, 7, 2, 6, 2, 5, 7, 6, 4, 6, 3,
1, 4, 5, 1, 4, 6, 2, 3, 0, 6, 11, 2, 9, 2, 6, 4, 5, 6, 2, 19, 2, 10, 4, 2,
3, 3, 11, 7, 3, 3, 1, 5, 3, 6, 4, 3, 0, 6, 6, 6, 4, 2, 5, 2, 2, 2, 6, 10,
4, 9, 3, 7, 7, 0, 6, 8, 5, 2, 3, 2, 3, 3, 3, 1, 6, 1, 8, 2, 5, 3, 6, 11,
5, 7, 2, 6, 7, 3, 4, 1, 0, 5, 8, 3, 2, 9, 3, 1, 2, 3, 3, 9, 5, 6, 5, 1,
4, 5, 6, 7, 6, 1, 5, 1, 6, 6, 2, 6, 7, 2, 4, 6
};
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vtkm::Id3(xVerts, yVerts, 1));
dataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", coordinates));
// Set point scalars
dataSet.AddField(vtkm::cont::Field("fieldA", vtkm::cont::Field::ASSOC_POINTS, fieldA, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldB", vtkm::cont::Field::ASSOC_POINTS, fieldB, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldC", vtkm::cont::Field::ASSOC_POINTS, fieldC, nVerts));
return dataSet;
}
//
// Create a dataset with known point data
// Extract the three field
// Run NDimsEntropy to calculate the entropy
//
void TestNDimsEntropy()
{
// Data attached is the poisson distribution
vtkm::cont::DataSet ds = MakeTestDataSet();
vtkm::worklet::NDimsEntropy ndEntropy;
ndEntropy.SetNumOfDataPoints(ds.GetField(0).GetData().GetNumberOfValues(),
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Add field one by one
ndEntropy.AddField(ds.GetField("fieldA").GetData(), 10, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
ndEntropy.AddField(ds.GetField("fieldB").GetData(), 10, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
ndEntropy.AddField(ds.GetField("fieldC").GetData(), 10, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Run worklet to calculate multi-variate entropy
vtkm::Float64 entropy = ndEntropy.Run(VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
VTKM_TEST_ASSERT(fabs(entropy - 7.457857) < 0.001,
"N-Dimentional entropy calculation is incorrect");
} // TestNDimsEntropy
}
int UnitTestNDimsEntropy(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestNDimsEntropy);
}

@ -0,0 +1,317 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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/worklet/NDimsHistMarginalization.h>
#include <vtkm/worklet/NDimsHistogram.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
// Make testing dataset with three fields(variables), each one has 1000 values
vtkm::cont::DataSet MakeTestDataSet()
{
vtkm::cont::DataSet dataSet;
const int xVerts = 20;
const int yVerts = 50;
const int nVerts = xVerts * yVerts;
vtkm::Float32 fieldA[nVerts] = {
8, 10, 9, 8, 14, 11, 12, 9, 19, 7, 8, 11, 7, 10, 11, 11, 11, 6, 8, 8, 7, 15, 9, 7,
8, 10, 9, 10, 10, 12, 7, 6, 14, 10, 14, 10, 7, 11, 13, 9, 13, 11, 10, 10, 12, 12, 7, 12,
10, 11, 12, 8, 13, 9, 5, 12, 11, 9, 5, 9, 12, 9, 6, 10, 11, 9, 9, 11, 9, 7, 7, 18,
16, 13, 12, 8, 10, 11, 9, 8, 17, 3, 15, 15, 9, 10, 10, 8, 10, 9, 7, 9, 8, 10, 13, 9,
7, 11, 7, 10, 13, 10, 11, 9, 10, 7, 10, 6, 12, 6, 9, 7, 6, 12, 12, 9, 12, 12, 11, 6,
1, 12, 8, 13, 14, 8, 8, 10, 7, 7, 6, 7, 5, 11, 6, 11, 13, 8, 13, 5, 9, 12, 7, 11,
10, 15, 11, 9, 7, 12, 15, 7, 8, 7, 12, 8, 21, 16, 13, 11, 10, 14, 12, 11, 12, 14, 7, 11,
7, 12, 16, 8, 10, 8, 9, 7, 8, 7, 13, 13, 11, 15, 7, 7, 6, 11, 7, 12, 12, 13, 14, 11,
13, 11, 11, 9, 15, 8, 6, 11, 12, 10, 11, 7, 6, 14, 11, 10, 12, 5, 8, 9, 11, 15, 11, 10,
17, 14, 9, 10, 10, 12, 11, 13, 13, 12, 11, 7, 8, 10, 7, 11, 10, 5, 8, 10, 13, 13, 12, 6,
10, 7, 13, 8, 11, 7, 10, 7, 8, 7, 14, 16, 9, 11, 8, 11, 9, 15, 11, 10, 10, 12, 7, 7,
11, 7, 5, 17, 9, 11, 11, 11, 10, 17, 10, 15, 7, 11, 12, 16, 9, 8, 11, 14, 9, 22, 8, 8,
8, 13, 12, 12, 1, 14, 15, 6, 15, 8, 11, 16, 14, 8, 6, 9, 8, 9, 9, 10, 8, 6, 13, 8,
6, 12, 11, 12, 13, 8, 6, 6, 5, 6, 10, 9, 11, 12, 14, 12, 10, 11, 10, 10, 8, 13, 8, 11,
7, 13, 13, 12, 12, 13, 15, 4, 9, 16, 7, 9, 8, 10, 6, 9, 11, 12, 6, 7, 14, 6, 4, 15,
5, 18, 9, 9, 11, 12, 9, 5, 6, 7, 15, 6, 11, 14, 8, 12, 6, 9, 5, 9, 14, 9, 12, 6,
9, 14, 11, 12, 12, 13, 15, 9, 8, 7, 13, 12, 7, 13, 6, 9, 10, 10, 10, 9, 11, 5, 9, 13,
16, 9, 10, 8, 9, 6, 13, 12, 8, 12, 9, 12, 17, 8, 11, 10, 8, 7, 11, 7, 13, 13, 10, 14,
11, 9, 6, 6, 14, 16, 5, 9, 13, 11, 12, 7, 4, 6, 9, 11, 11, 10, 12, 9, 7, 13, 8, 8,
12, 5, 10, 7, 11, 11, 10, 10, 14, 6, 8, 8, 3, 12, 16, 11, 11, 7, 6, 12, 11, 5, 9, 12,
9, 13, 7, 8, 9, 9, 12, 7, 9, 8, 12, 11, 6, 10, 6, 7, 6, 11, 10, 8, 9, 8, 4, 19,
12, 6, 10, 9, 6, 12, 9, 14, 7, 8, 11, 7, 7, 12, 13, 9, 13, 12, 8, 6, 10, 17, 19, 10,
10, 13, 5, 11, 8, 10, 8, 16, 12, 6, 6, 7, 10, 9, 12, 8, 5, 10, 7, 18, 9, 12, 10, 4,
9, 9, 15, 15, 6, 7, 7, 11, 12, 4, 8, 18, 5, 12, 12, 11, 10, 14, 9, 9, 10, 8, 10, 8,
10, 9, 9, 4, 10, 12, 5, 13, 6, 9, 7, 5, 12, 8, 11, 10, 9, 17, 9, 9, 8, 11, 18, 11,
10, 9, 4, 13, 10, 15, 5, 10, 9, 7, 7, 8, 10, 6, 6, 19, 10, 16, 7, 7, 9, 10, 10, 13,
10, 10, 14, 13, 12, 8, 7, 13, 12, 11, 13, 12, 9, 8, 6, 8, 10, 3, 8, 8, 12, 12, 13, 13,
10, 5, 10, 7, 13, 7, 9, 5, 13, 7, 10, 8, 13, 11, 17, 9, 6, 14, 10, 10, 13, 9, 15, 8,
15, 9, 12, 11, 12, 8, 3, 9, 8, 10, 12, 8, 14, 13, 12, 11, 12, 9, 18, 10, 13, 7, 4, 4,
11, 8, 3, 7, 9, 10, 12, 7, 11, 21, 9, 7, 8, 9, 10, 10, 11, 9, 15, 13, 21, 12, 8, 11,
9, 10, 11, 9, 17, 8, 9, 8, 14, 6, 13, 9, 8, 11, 12, 12, 12, 11, 6, 13, 7, 9, 11, 15,
17, 17, 11, 10, 7, 8, 11, 8, 6, 9, 13, 7, 9, 6, 5, 10, 7, 16, 16, 9, 7, 6, 14, 8,
13, 16, 7, 7, 10, 11, 6, 10, 9, 9, 8, 14, 11, 9, 11, 9, 10, 11, 9, 8, 14, 11, 7, 12,
11, 8, 9, 9, 10, 11, 11, 10, 9, 6, 6, 11, 16, 10, 7, 6, 6, 13, 18, 8, 12, 11, 14, 13,
8, 8, 10, 17, 17, 6, 6, 10, 18, 5, 8, 11, 6, 6, 14, 10, 9, 6, 11, 6, 13, 12, 10, 6,
9, 9, 9, 13, 7, 17, 10, 14, 10, 9, 10, 10, 11, 10, 11, 15, 13, 6, 12, 19, 10, 12, 12, 15,
13, 10, 10, 13, 11, 13, 13, 17, 6, 5, 6, 7, 6, 9, 13, 11, 8, 12, 9, 6, 10, 16, 11, 12,
5, 12, 14, 13, 13, 16, 11, 6, 12, 12, 15, 8, 7, 11, 8, 5, 10, 8, 9, 11, 9, 12, 10, 5,
12, 11, 9, 6, 14, 12, 10, 11, 9, 6, 7, 12, 8, 12, 8, 15, 9, 8, 7, 9, 3, 6, 14, 7,
8, 11, 9, 10, 12, 9, 10, 9, 8, 6, 12, 11, 6, 8, 9, 8, 15, 11, 7, 18, 12, 11, 10, 13,
11, 11, 10, 7, 9, 8, 8, 11, 11, 13, 6, 12, 13, 16, 11, 11, 5, 12, 14, 15, 9, 14, 15, 6,
8, 7, 6, 8, 9, 19, 7, 12, 11, 8, 14, 12, 10, 9, 3, 7
};
vtkm::Float32 fieldB[nVerts] = {
24, 19, 28, 19, 25, 28, 25, 22, 27, 26, 35, 26, 30, 28, 24, 23, 21, 31, 20, 11, 21, 22, 14, 25,
20, 24, 24, 21, 24, 29, 26, 21, 32, 29, 23, 28, 31, 25, 23, 30, 18, 24, 22, 25, 33, 24, 22, 23,
21, 17, 20, 28, 30, 18, 20, 32, 25, 24, 32, 15, 27, 24, 27, 19, 30, 27, 17, 24, 29, 23, 22, 19,
24, 19, 28, 24, 25, 24, 25, 30, 24, 31, 30, 27, 25, 25, 25, 15, 29, 23, 29, 29, 21, 25, 35, 24,
28, 10, 31, 23, 22, 22, 22, 33, 29, 27, 18, 27, 27, 24, 20, 20, 21, 29, 23, 31, 23, 23, 22, 23,
30, 27, 28, 31, 16, 29, 25, 19, 33, 28, 25, 24, 15, 27, 37, 29, 15, 19, 14, 19, 24, 23, 30, 29,
35, 22, 19, 26, 26, 14, 24, 30, 32, 23, 30, 29, 26, 27, 25, 23, 17, 26, 32, 29, 20, 17, 21, 23,
22, 20, 36, 12, 26, 23, 15, 29, 24, 22, 26, 33, 24, 23, 20, 26, 22, 17, 26, 26, 34, 22, 26, 17,
23, 18, 29, 27, 21, 29, 28, 29, 24, 25, 28, 19, 18, 21, 23, 23, 27, 25, 24, 25, 24, 25, 21, 25,
21, 27, 23, 20, 29, 15, 28, 30, 24, 27, 17, 23, 16, 21, 25, 17, 27, 28, 21, 13, 19, 27, 16, 30,
31, 25, 30, 17, 17, 25, 26, 22, 21, 17, 24, 17, 25, 22, 27, 14, 27, 24, 27, 25, 26, 31, 21, 23,
30, 30, 22, 19, 23, 22, 23, 25, 24, 25, 24, 28, 26, 30, 18, 25, 30, 37, 27, 34, 28, 34, 25, 10,
25, 22, 35, 30, 24, 32, 24, 34, 19, 29, 26, 16, 27, 17, 26, 23, 27, 25, 26, 21, 31, 21, 28, 15,
32, 24, 23, 23, 18, 15, 22, 25, 16, 25, 31, 26, 25, 28, 24, 26, 23, 25, 33, 20, 27, 28, 24, 29,
32, 20, 24, 20, 19, 32, 24, 6, 24, 21, 26, 18, 15, 30, 19, 26, 22, 30, 35, 23, 22, 30, 20, 22,
18, 30, 28, 25, 16, 25, 27, 30, 18, 24, 30, 28, 20, 19, 20, 28, 21, 24, 15, 33, 20, 18, 20, 36,
30, 26, 25, 18, 28, 27, 31, 31, 15, 26, 16, 22, 27, 14, 17, 27, 27, 22, 32, 30, 22, 34, 22, 25,
20, 22, 26, 29, 28, 33, 18, 23, 20, 20, 27, 24, 28, 21, 25, 27, 25, 19, 19, 25, 19, 32, 29, 27,
23, 21, 28, 33, 23, 23, 28, 26, 31, 19, 21, 29, 21, 27, 23, 32, 24, 26, 21, 28, 28, 24, 17, 31,
27, 21, 19, 32, 28, 23, 30, 23, 29, 15, 26, 26, 15, 20, 25, 26, 27, 31, 21, 23, 23, 33, 28, 19,
23, 22, 22, 25, 27, 17, 23, 17, 25, 28, 26, 30, 32, 31, 19, 25, 25, 19, 23, 29, 27, 23, 34, 22,
13, 21, 32, 10, 20, 33, 21, 17, 29, 31, 14, 24, 23, 19, 19, 22, 17, 26, 37, 26, 22, 26, 38, 29,
29, 27, 30, 20, 31, 14, 32, 32, 24, 23, 23, 18, 21, 31, 24, 20, 28, 15, 21, 25, 25, 20, 30, 25,
22, 21, 21, 25, 24, 25, 18, 23, 28, 30, 20, 27, 27, 19, 10, 32, 24, 20, 29, 26, 25, 20, 25, 29,
28, 24, 32, 26, 22, 19, 23, 27, 27, 29, 20, 25, 21, 30, 28, 31, 24, 19, 23, 19, 19, 18, 30, 18,
16, 24, 20, 20, 30, 25, 29, 25, 31, 21, 28, 31, 24, 26, 27, 21, 24, 23, 26, 18, 32, 26, 28, 26,
24, 26, 29, 30, 22, 20, 24, 28, 25, 29, 20, 21, 22, 15, 30, 27, 33, 26, 22, 32, 30, 31, 20, 19,
24, 26, 27, 31, 17, 17, 33, 27, 16, 27, 27, 22, 27, 19, 24, 21, 17, 24, 28, 23, 26, 24, 19, 26,
20, 24, 22, 19, 22, 21, 21, 28, 29, 39, 19, 16, 25, 29, 31, 22, 22, 29, 26, 22, 22, 22, 26, 23,
23, 23, 30, 25, 25, 25, 27, 29, 18, 33, 21, 12, 22, 29, 12, 20, 35, 22, 34, 28, 18, 29, 21, 20,
24, 33, 24, 26, 23, 34, 31, 25, 31, 22, 35, 21, 20, 29, 27, 22, 30, 22, 27, 23, 22, 32, 16, 19,
27, 22, 24, 27, 21, 33, 25, 25, 19, 28, 20, 27, 21, 25, 28, 20, 27, 22, 21, 20, 26, 30, 33, 23,
20, 24, 17, 23, 28, 35, 14, 23, 22, 28, 28, 26, 25, 18, 20, 28, 28, 22, 13, 24, 22, 20, 30, 26,
26, 18, 22, 20, 23, 24, 20, 27, 34, 28, 18, 24, 34, 33, 25, 33, 37, 21, 20, 31, 19, 23, 29, 22,
21, 24, 19, 27, 19, 32, 25, 23, 33, 26, 33, 27, 29, 30, 19, 22, 30, 19, 18, 24, 25, 17, 31, 19,
31, 26, 22, 23, 28, 28, 25, 24, 19, 19, 27, 28, 23, 21, 29, 26, 31, 22, 22, 25, 16, 29, 21, 22,
23, 25, 22, 21, 22, 19, 27, 26, 28, 30, 22, 21, 24, 22, 23, 26, 28, 22, 18, 25, 23, 27, 31, 19,
15, 29, 20, 19, 27, 25, 21, 29, 22, 24, 25, 17, 36, 29, 22, 22, 24, 28, 27, 22, 26, 31, 29, 31,
18, 25, 23, 16, 37, 27, 21, 31, 25, 24, 20, 23, 28, 33, 24, 21, 26, 20, 18, 31, 20, 24, 23, 19,
27, 17, 23, 23, 20, 26, 28, 23, 26, 31, 25, 31, 19, 32, 26, 18, 19, 29, 20, 21, 15, 25, 27, 29,
22, 22, 22, 26, 23, 22, 23, 29, 28, 20, 21, 22, 20, 22, 27, 25, 23, 32, 23, 20, 31, 20, 27, 26,
34, 20, 22, 36, 21, 29, 25, 20, 21, 22, 29, 29, 25, 22, 24, 22
};
vtkm::Float32 fieldC[nVerts] = {
3, 1, 4, 6, 5, 4, 8, 7, 2, 9, 2, 0, 0, 4, 3, 2, 5, 2, 3, 6, 3, 8, 3, 4,
3, 3, 2, 7, 2, 10, 9, 6, 1, 1, 4, 7, 3, 3, 1, 4, 4, 3, 9, 4, 4, 7, 3, 2,
4, 7, 3, 3, 2, 10, 1, 6, 2, 2, 3, 8, 3, 3, 6, 9, 4, 1, 4, 3, 16, 7, 0, 1,
8, 7, 13, 3, 5, 0, 3, 8, 10, 3, 5, 5, 1, 5, 2, 1, 3, 2, 5, 3, 4, 3, 3, 3,
3, 1, 13, 2, 3, 1, 2, 7, 3, 4, 1, 2, 5, 4, 4, 4, 2, 6, 3, 2, 7, 8, 1, 3,
4, 1, 2, 0, 1, 6, 1, 8, 8, 1, 1, 4, 2, 1, 4, 3, 5, 4, 6, 4, 2, 3, 8, 8,
3, 3, 3, 4, 5, 8, 8, 16, 7, 12, 4, 3, 14, 8, 3, 12, 5, 0, 5, 3, 5, 2, 9, 2,
9, 4, 1, 0, 0, 4, 4, 6, 3, 4, 11, 2, 4, 7, 4, 2, 1, 9, 4, 3, 2, 5, 1, 5,
3, 8, 2, 8, 1, 8, 0, 4, 1, 3, 2, 1, 2, 3, 2, 1, 8, 5, 4, 1, 9, 9, 1, 3,
5, 0, 1, 6, 10, 8, 3, 12, 3, 4, 4, 7, 1, 3, 6, 4, 4, 6, 1, 4, 7, 5, 6, 11,
6, 5, 2, 7, 2, 5, 3, 5, 6, 3, 6, 2, 1, 10, 8, 3, 7, 0, 2, 6, 9, 3, 11, 3,
2, 5, 1, 4, 6, 10, 9, 1, 4, 3, 7, 12, 3, 10, 0, 2, 11, 2, 1, 0, 4, 1, 2, 16,
5, 17, 7, 8, 2, 10, 10, 3, 1, 3, 2, 2, 4, 8, 4, 3, 2, 4, 4, 6, 8, 6, 2, 3,
2, 4, 2, 4, 7, 10, 5, 3, 5, 2, 4, 6, 9, 3, 1, 1, 1, 1, 4, 2, 2, 7, 4, 9,
2, 3, 5, 6, 2, 5, 1, 6, 5, 7, 8, 3, 7, 2, 2, 8, 6, 2, 10, 2, 1, 4, 5, 1,
1, 1, 5, 6, 1, 1, 4, 5, 4, 2, 4, 3, 2, 7, 19, 4, 7, 2, 7, 5, 2, 5, 3, 8,
4, 6, 7, 2, 0, 0, 2, 12, 6, 2, 2, 3, 5, 9, 4, 9, 2, 2, 7, 8, 3, 3, 10, 6,
3, 2, 1, 6, 2, 4, 6, 3, 5, 8, 2, 3, 6, 14, 0, 3, 6, 5, 2, 7, 0, 3, 8, 5,
3, 2, 2, 5, 1, 3, 12, 11, 16, 2, 1, 3, 7, 3, 1, 6, 4, 3, 12, 5, 1, 3, 1, 4,
9, 1, 3, 3, 4, 4, 6, 7, 7, 5, 2, 4, 2, 3, 2, 2, 6, 4, 2, 2, 3, 5, 1, 4,
9, 1, 0, 7, 6, 4, 3, 3, 7, 3, 3, 6, 2, 7, 9, 3, 1, 16, 5, 4, 3, 6, 3, 2,
5, 2, 2, 4, 3, 1, 3, 3, 6, 3, 5, 9, 1, 10, 1, 7, 2, 2, 6, 7, 3, 5, 3, 7,
2, 2, 2, 2, 6, 4, 3, 2, 5, 5, 3, 15, 4, 2, 7, 7, 4, 3, 3, 5, 1, 2, 9, 0,
5, 7, 12, 2, 4, 8, 5, 7, 8, 3, 2, 2, 18, 1, 7, 2, 2, 1, 3, 3, 3, 7, 1, 9,
8, 4, 3, 7, 6, 4, 5, 2, 0, 5, 1, 5, 10, 4, 2, 8, 2, 2, 0, 5, 6, 4, 5, 0,
1, 5, 11, 3, 3, 4, 4, 2, 3, 5, 1, 6, 5, 7, 2, 2, 5, 7, 4, 8, 4, 1, 1, 7,
2, 3, 9, 6, 13, 1, 5, 4, 6, 2, 4, 11, 2, 5, 5, 1, 4, 1, 4, 7, 1, 5, 8, 3,
1, 10, 9, 13, 1, 7, 2, 9, 4, 3, 3, 10, 12, 2, 0, 4, 6, 5, 5, 1, 4, 7, 2, 12,
7, 6, 5, 0, 6, 4, 4, 12, 1, 3, 10, 1, 9, 2, 2, 2, 1, 5, 5, 6, 9, 6, 4, 1,
11, 6, 9, 3, 2, 7, 1, 7, 4, 3, 0, 3, 1, 12, 17, 2, 1, 6, 4, 4, 2, 1, 5, 5,
3, 2, 2, 4, 6, 5, 4, 6, 11, 3, 12, 6, 3, 6, 3, 0, 6, 3, 7, 4, 8, 5, 14, 5,
1, 9, 4, 6, 5, 3, 9, 3, 1, 1, 0, 3, 7, 3, 5, 1, 6, 2, 2, 6, 2, 12, 1, 0,
6, 3, 3, 5, 4, 7, 2, 2, 15, 7, 3, 10, 4, 2, 6, 3, 4, 8, 3, 1, 5, 5, 5, 4,
3, 7, 3, 4, 5, 5, 2, 4, 2, 5, 1, 12, 5, 6, 3, 2, 8, 5, 2, 3, 11, 11, 6, 5,
0, 3, 3, 9, 4, 2, 11, 1, 5, 3, 5, 6, 3, 6, 4, 2, 4, 10, 11, 3, 3, 4, 1, 1,
1, 3, 5, 5, 1, 1, 4, 1, 5, 1, 6, 8, 6, 4, 6, 7, 6, 3, 5, 3, 6, 6, 6, 4,
0, 6, 3, 1, 2, 4, 2, 6, 1, 1, 1, 2, 2, 4, 7, 2, 6, 2, 5, 7, 6, 4, 6, 3,
1, 4, 5, 1, 4, 6, 2, 3, 0, 6, 11, 2, 9, 2, 6, 4, 5, 6, 2, 19, 2, 10, 4, 2,
3, 3, 11, 7, 3, 3, 1, 5, 3, 6, 4, 3, 0, 6, 6, 6, 4, 2, 5, 2, 2, 2, 6, 10,
4, 9, 3, 7, 7, 0, 6, 8, 5, 2, 3, 2, 3, 3, 3, 1, 6, 1, 8, 2, 5, 3, 6, 11,
5, 7, 2, 6, 7, 3, 4, 1, 0, 5, 8, 3, 2, 9, 3, 1, 2, 3, 3, 9, 5, 6, 5, 1,
4, 5, 6, 7, 6, 1, 5, 1, 6, 6, 2, 6, 7, 2, 4, 6
};
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vtkm::Id3(xVerts, yVerts, 1));
dataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", coordinates));
// Set point scalars
dataSet.AddField(vtkm::cont::Field("fieldA", vtkm::cont::Field::ASSOC_POINTS, fieldA, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldB", vtkm::cont::Field::ASSOC_POINTS, fieldB, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldC", vtkm::cont::Field::ASSOC_POINTS, fieldC, nVerts));
return dataSet;
}
// The Condition function for non-marginal variable
// var is index of variable (data array)
// binId is bin index
// return true indicates considering this bin into final
// marginal histogram
// In this example, we have three variable var0, var1, var2
// the condition is P(Var0, Var2 | 1<Var1<4)
// because var0 and var2 are the marginal varaible, we do not care the case var==0 or var==2
// it supposes that we do not have input as var ==0 or 2
struct VariableCondition
{
VTKM_EXEC
bool operator()(vtkm::Id var, vtkm::Id binId) const
{
if (var == 1)
{
if (1 < binId && binId < 4)
return true;
}
return false;
}
};
//
// Create a dataset with known point data and cell data (statistical distributions)
// Extract arrays of point and cell fields
// Create output structure to hold histogram bins
// Run FieldHistogram filter
//
void TestNDimsHistMarginalization()
{
// Setup and calculate a ND-histogram without marginalization first
// Create the output bin array
vtkm::cont::ArrayHandle<vtkm::Range> range;
vtkm::cont::ArrayHandle<vtkm::Float64> delta;
// Data attached is the poisson distribution
vtkm::cont::DataSet ds = MakeTestDataSet();
// Compute Nd histogram first
vtkm::worklet::NDimsHistogram ndHistogram;
// Set the number of data points
ndHistogram.SetNumOfDataPoints(ds.GetField(0).GetData().GetNumberOfValues(),
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Add field one by one
vtkm::Range rangeFieldA;
vtkm::Float64 deltaFieldA;
ndHistogram.AddField(ds.GetField("fieldA").GetData(),
10,
rangeFieldA,
deltaFieldA,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::Range rangeFieldB;
vtkm::Float64 deltaFieldB;
ndHistogram.AddField(ds.GetField("fieldB").GetData(),
10,
rangeFieldB,
deltaFieldB,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::Range rangeFieldC;
vtkm::Float64 deltaFieldC;
ndHistogram.AddField(ds.GetField("fieldC").GetData(),
10,
rangeFieldC,
deltaFieldC,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// the return binIds and freqs is sparse distribution representation
// (we do not keep the 0 frequency entities)
// e.g. we have three variable(data arrays) in this example
// binIds[0, 1, 2][j] is a combination of bin ID of three variable,
// freqs[j] is the freqncy of this bin IDs combination
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> binIds;
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
ndHistogram.Run(binIds, freqs, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// setup for histogram marginalization
// use a bool array to indicate the marginal variable (true -> marginal variable)
// the length of this array has to be equal to number of input variables
bool marginalVariableAry[3] = { true, false, true };
vtkm::cont::ArrayHandle<bool> marginalVariable =
vtkm::cont::make_ArrayHandle(marginalVariableAry, 3);
std::vector<vtkm::Id> numberOfBinsVec(3, 10);
vtkm::cont::ArrayHandle<vtkm::Id> numberOfBins = vtkm::cont::make_ArrayHandle(numberOfBinsVec);
// calculate marginal histogram by the setup (return sparse represetnation)
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> marginalBinIds;
vtkm::cont::ArrayHandle<vtkm::Id> marginalFreqs;
vtkm::worklet::NDimsHistMarginalization ndHistMarginalization;
ndHistMarginalization.Run(binIds,
freqs,
numberOfBins,
marginalVariable,
VariableCondition(),
marginalBinIds,
marginalFreqs,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Ground truth ND histogram
vtkm::Id gtNonSparseBins = 40;
vtkm::Id gtIdx0[40] = { 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4,
4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 8, 9 };
vtkm::Id gtIdx1[40] = { 1, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 7, 0, 1, 2, 3, 4, 5, 0, 1,
2, 3, 4, 5, 7, 8, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 1, 2, 0, 1 };
vtkm::Id gtFreq[40] = { 1, 2, 1, 2, 1, 4, 7, 6, 3, 2, 2, 1, 6, 6, 8, 6, 2, 2, 6, 9,
10, 2, 5, 1, 1, 1, 6, 7, 9, 6, 3, 3, 2, 3, 2, 2, 3, 2, 1, 1 };
// Check result
vtkm::Id nonSparseBins = marginalBinIds[0].GetPortalControl().GetNumberOfValues();
VTKM_TEST_ASSERT(nonSparseBins == gtNonSparseBins,
"Incorrect ND-histogram marginalization results");
for (int i = 0; i < nonSparseBins; i++)
{
vtkm::Id idx0 = marginalBinIds[0].GetPortalControl().Get(i);
vtkm::Id idx1 = marginalBinIds[1].GetPortalControl().Get(i);
vtkm::Id f = marginalFreqs.GetPortalControl().Get(i);
VTKM_TEST_ASSERT(idx0 == gtIdx0[i] && idx1 == gtIdx1[i] && f == gtFreq[i],
"Incorrect ND-histogram marginalization results");
}
} // TestNDimsHistMarginalization
}
int UnitTestNDimsHistMarginalization(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestNDimsHistMarginalization);
}

@ -0,0 +1,141 @@
//============================================================================
// 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 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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/worklet/NDimsHistogram.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
// Make testing dataset with three fields(variables), each one has 100 values
vtkm::cont::DataSet MakeTestDataSet()
{
vtkm::cont::DataSet dataSet;
const int nVerts = 100;
vtkm::Float32 fieldA[nVerts] = { 8, 10, 9, 8, 14, 11, 12, 9, 19, 7, 8, 11, 7, 10, 11,
11, 11, 6, 8, 8, 7, 15, 9, 7, 8, 10, 9, 10, 10, 12,
7, 6, 14, 10, 14, 10, 7, 11, 13, 9, 13, 11, 10, 10, 12,
12, 7, 12, 10, 11, 12, 8, 13, 9, 5, 12, 11, 9, 5, 9,
12, 9, 6, 10, 11, 9, 9, 11, 9, 7, 7, 18, 16, 13, 12,
8, 10, 11, 9, 8, 17, 3, 15, 15, 9, 10, 10, 8, 10, 9,
7, 9, 8, 10, 13, 9, 7, 11, 7, 10 };
vtkm::Float32 fieldB[nVerts] = { 24, 19, 28, 19, 25, 28, 25, 22, 27, 26, 35, 26, 30, 28, 24,
23, 21, 31, 20, 11, 21, 22, 14, 25, 20, 24, 24, 21, 24, 29,
26, 21, 32, 29, 23, 28, 31, 25, 23, 30, 18, 24, 22, 25, 33,
24, 22, 23, 21, 17, 20, 28, 30, 18, 20, 32, 25, 24, 32, 15,
27, 24, 27, 19, 30, 27, 17, 24, 29, 23, 22, 19, 24, 19, 28,
24, 25, 24, 25, 30, 24, 31, 30, 27, 25, 25, 25, 15, 29, 23,
29, 29, 21, 25, 35, 24, 28, 10, 31, 23 };
vtkm::Float32 fieldC[nVerts] = {
3, 1, 4, 6, 5, 4, 8, 7, 2, 9, 2, 0, 0, 4, 3, 2, 5, 2, 3, 6, 3, 8, 3, 4, 3,
3, 2, 7, 2, 10, 9, 6, 1, 1, 4, 7, 3, 3, 1, 4, 4, 3, 9, 4, 4, 7, 3, 2, 4, 7,
3, 3, 2, 10, 1, 6, 2, 2, 3, 8, 3, 3, 6, 9, 4, 1, 4, 3, 16, 7, 0, 1, 8, 7, 13,
3, 5, 0, 3, 8, 10, 3, 5, 5, 1, 5, 2, 1, 3, 2, 5, 3, 4, 3, 3, 3, 3, 1, 13, 2
};
// Set point scalars
dataSet.AddField(vtkm::cont::Field("fieldA", vtkm::cont::Field::ASSOC_POINTS, fieldA, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldB", vtkm::cont::Field::ASSOC_POINTS, fieldB, nVerts));
dataSet.AddField(vtkm::cont::Field("fieldC", vtkm::cont::Field::ASSOC_POINTS, fieldC, nVerts));
return dataSet;
}
void TestNDimsHistogram()
{
// Create a dataset
vtkm::cont::DataSet ds = MakeTestDataSet();
vtkm::worklet::NDimsHistogram ndHistogram;
// Set the number of data points
ndHistogram.SetNumOfDataPoints(ds.GetField(0).GetData().GetNumberOfValues(),
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Add field one by one
vtkm::Range rangeFieldA;
vtkm::Float64 deltaFieldA;
ndHistogram.AddField(ds.GetField("fieldA").GetData(),
4,
rangeFieldA,
deltaFieldA,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::Range rangeFieldB;
vtkm::Float64 deltaFieldB;
ndHistogram.AddField(ds.GetField("fieldB").GetData(),
4,
rangeFieldB,
deltaFieldB,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::Range rangeFieldC;
vtkm::Float64 deltaFieldC;
ndHistogram.AddField(ds.GetField("fieldC").GetData(),
4,
rangeFieldC,
deltaFieldC,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// the return binIds and freqs is sparse distribution representation
// (we do not keep the 0 frequency entities)
// e.g. we have three variable(data arrays) in this example
// binIds[0, 1, 2][j] is a combination of bin ID of three variable,
// freqs[j] is the freqncy of this bin IDs combination
std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> binIds;
vtkm::cont::ArrayHandle<vtkm::Id> freqs;
ndHistogram.Run(binIds, freqs, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// Ground truth ND histogram
vtkm::Id gtNonSparseBins = 33;
vtkm::Id gtIdx0[33] = { 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3 };
vtkm::Id gtIdx1[33] = { 1, 1, 2, 3, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3,
0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 1, 1, 2, 2, 2, 3 };
vtkm::Id gtIdx2[33] = { 0, 1, 1, 0, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 3,
0, 0, 1, 0, 1, 2, 3, 0, 1, 2, 0, 2, 0, 1, 2, 1 };
vtkm::Id gtFreq[33] = { 1, 1, 1, 3, 2, 1, 1, 6, 6, 3, 17, 8, 2, 6, 2, 1, 2,
1, 1, 4, 11, 4, 1, 1, 3, 3, 1, 1, 1, 1, 1, 2, 1 };
// Check result
vtkm::Id nonSparseBins = binIds[0].GetPortalControl().GetNumberOfValues();
VTKM_TEST_ASSERT(nonSparseBins == gtNonSparseBins, "Incorrect ND-histogram results");
for (int i = 0; i < nonSparseBins; i++)
{
vtkm::Id idx0 = binIds[0].GetPortalControl().Get(i);
vtkm::Id idx1 = binIds[1].GetPortalControl().Get(i);
vtkm::Id idx2 = binIds[2].GetPortalControl().Get(i);
vtkm::Id f = freqs.GetPortalControl().Get(i);
VTKM_TEST_ASSERT(idx0 == gtIdx0[i] && idx1 == gtIdx1[i] && idx2 == gtIdx2[i] && f == gtFreq[i],
"Incorrect ND-histogram results");
}
} // TestNDHistogram
}
int UnitTestNDimsHistogram(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestNDimsHistogram);
}