2015-09-10 17:13:08 +00:00
|
|
|
//============================================================================
|
|
|
|
// Copyright (c) Kitware, Inc.
|
|
|
|
// All rights reserved.
|
|
|
|
// See LICENSE.txt for details.
|
2019-04-15 23:24:21 +00:00
|
|
|
//
|
2015-09-10 17:13:08 +00:00
|
|
|
// 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.
|
|
|
|
//============================================================================
|
|
|
|
|
|
|
|
#ifndef vtk_m_worklet_FieldStatistics_h
|
|
|
|
#define vtk_m_worklet_FieldStatistics_h
|
|
|
|
|
|
|
|
#include <vtkm/Math.h>
|
2018-10-17 19:17:29 +00:00
|
|
|
#include <vtkm/cont/Algorithm.h>
|
2019-08-22 15:12:59 +00:00
|
|
|
#include <vtkm/cont/ArrayGetValues.h>
|
2015-09-10 17:13:08 +00:00
|
|
|
#include <vtkm/cont/ArrayHandle.h>
|
2017-05-18 14:51:24 +00:00
|
|
|
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
2015-09-10 17:13:08 +00:00
|
|
|
#include <vtkm/worklet/DispatcherMapField.h>
|
|
|
|
#include <vtkm/worklet/WorkletMapField.h>
|
|
|
|
|
|
|
|
#include <vtkm/cont/Field.h>
|
|
|
|
|
|
|
|
#include <stdio.h>
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
namespace vtkm
|
|
|
|
{
|
|
|
|
namespace worklet
|
|
|
|
{
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
//simple functor that prints basic statistics
|
2018-10-17 19:17:29 +00:00
|
|
|
template <typename FieldType>
|
2015-09-10 17:13:08 +00:00
|
|
|
class FieldStatistics
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
// For moments readability
|
2018-02-27 14:25:25 +00:00
|
|
|
static constexpr vtkm::Id FIRST = 0;
|
|
|
|
static constexpr vtkm::Id SECOND = 1;
|
|
|
|
static constexpr vtkm::Id THIRD = 2;
|
|
|
|
static constexpr vtkm::Id FOURTH = 3;
|
|
|
|
static constexpr vtkm::Id NUM_POWERS = 4;
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
struct StatInfo
|
|
|
|
{
|
|
|
|
FieldType minimum;
|
|
|
|
FieldType maximum;
|
|
|
|
FieldType median;
|
|
|
|
FieldType mean;
|
|
|
|
FieldType variance;
|
|
|
|
FieldType stddev;
|
|
|
|
FieldType skewness;
|
|
|
|
FieldType kurtosis;
|
|
|
|
FieldType rawMoment[4];
|
|
|
|
FieldType centralMoment[4];
|
|
|
|
};
|
|
|
|
|
|
|
|
class CalculatePowers : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
public:
|
2019-01-10 18:59:25 +00:00
|
|
|
using ControlSignature = void(FieldIn value,
|
|
|
|
FieldOut pow1Array,
|
|
|
|
FieldOut pow2Array,
|
|
|
|
FieldOut pow3Array,
|
|
|
|
FieldOut pow4Array);
|
2018-05-25 21:18:41 +00:00
|
|
|
using ExecutionSignature = void(_1, _2, _3, _4, _5);
|
2018-02-22 13:29:13 +00:00
|
|
|
using InputDomain = _1;
|
2016-11-23 14:10:57 +00:00
|
|
|
|
2015-09-10 17:13:08 +00:00
|
|
|
vtkm::Id numPowers;
|
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2017-05-18 14:29:41 +00:00
|
|
|
CalculatePowers(vtkm::Id num)
|
|
|
|
: numPowers(num)
|
|
|
|
{
|
|
|
|
}
|
2016-11-23 14:10:57 +00:00
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2017-05-26 17:53:28 +00:00
|
|
|
void operator()(const FieldType& value,
|
|
|
|
FieldType& pow1,
|
|
|
|
FieldType& pow2,
|
|
|
|
FieldType& pow3,
|
2017-05-18 14:29:41 +00:00
|
|
|
FieldType& pow4) const
|
2015-09-10 17:13:08 +00:00
|
|
|
{
|
|
|
|
pow1 = value;
|
|
|
|
pow2 = pow1 * value;
|
|
|
|
pow3 = pow2 * value;
|
|
|
|
pow4 = pow3 * value;
|
|
|
|
}
|
|
|
|
};
|
|
|
|
|
|
|
|
class SubtractConst : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
public:
|
2019-01-10 18:59:25 +00:00
|
|
|
using ControlSignature = void(FieldIn value, FieldOut diff);
|
2018-05-25 21:18:41 +00:00
|
|
|
using ExecutionSignature = _2(_1);
|
2018-02-22 13:29:13 +00:00
|
|
|
using InputDomain = _1;
|
2016-11-23 14:10:57 +00:00
|
|
|
|
2015-09-10 17:13:08 +00:00
|
|
|
FieldType constant;
|
2016-11-23 14:10:57 +00:00
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2017-05-18 14:29:41 +00:00
|
|
|
SubtractConst(const FieldType& constant0)
|
|
|
|
: constant(constant0)
|
2015-09-10 17:13:08 +00:00
|
|
|
{
|
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
|
|
|
|
VTKM_EXEC
|
|
|
|
FieldType operator()(const FieldType& value) const { return (value - constant); }
|
2015-09-10 17:13:08 +00:00
|
|
|
};
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename Storage>
|
|
|
|
void Run(vtkm::cont::ArrayHandle<FieldType, Storage> fieldArray, StatInfo& statinfo)
|
2015-09-10 17:13:08 +00:00
|
|
|
{
|
2018-10-17 19:17:29 +00:00
|
|
|
using DeviceAlgorithms = vtkm::cont::Algorithm;
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
// Copy original data to array for sorting
|
|
|
|
vtkm::cont::ArrayHandle<FieldType> tempArray;
|
|
|
|
DeviceAlgorithms::Copy(fieldArray, tempArray);
|
|
|
|
DeviceAlgorithms::Sort(tempArray);
|
|
|
|
|
2019-08-22 15:12:59 +00:00
|
|
|
vtkm::Id dataSize = tempArray.GetNumberOfValues();
|
2015-09-10 17:13:08 +00:00
|
|
|
FieldType numValues = static_cast<FieldType>(dataSize);
|
2019-08-22 15:12:59 +00:00
|
|
|
const auto firstAndMedian = vtkm::cont::ArrayGetValues({ 0, dataSize / 2 }, tempArray);
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
// Median
|
2019-08-22 15:12:59 +00:00
|
|
|
statinfo.median = firstAndMedian[1];
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
// Minimum and maximum
|
2019-08-22 15:12:59 +00:00
|
|
|
const vtkm::Vec<FieldType, 2> initValue(firstAndMedian[0]);
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Vec<FieldType, 2> result =
|
|
|
|
DeviceAlgorithms::Reduce(fieldArray, initValue, vtkm::MinAndMax<FieldType>());
|
2016-11-25 18:04:26 +00:00
|
|
|
statinfo.minimum = result[0];
|
|
|
|
statinfo.maximum = result[1];
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
// Mean
|
|
|
|
FieldType sum = DeviceAlgorithms::ScanInclusive(fieldArray, tempArray);
|
|
|
|
statinfo.mean = sum / numValues;
|
|
|
|
statinfo.rawMoment[FIRST] = sum / numValues;
|
|
|
|
|
|
|
|
// Create the power sum vector for each value
|
|
|
|
vtkm::cont::ArrayHandle<FieldType> pow1Array, pow2Array, pow3Array, pow4Array;
|
|
|
|
pow1Array.Allocate(dataSize);
|
|
|
|
pow2Array.Allocate(dataSize);
|
|
|
|
pow3Array.Allocate(dataSize);
|
|
|
|
pow4Array.Allocate(dataSize);
|
|
|
|
|
|
|
|
// Raw moments via Worklet
|
2018-08-28 20:36:50 +00:00
|
|
|
vtkm::worklet::DispatcherMapField<CalculatePowers> calculatePowersDispatcher(
|
2017-05-18 14:29:41 +00:00
|
|
|
CalculatePowers(4));
|
2015-09-10 17:13:08 +00:00
|
|
|
calculatePowersDispatcher.Invoke(fieldArray, pow1Array, pow2Array, pow3Array, pow4Array);
|
|
|
|
|
|
|
|
// Accumulate the results using ScanInclusive
|
|
|
|
statinfo.rawMoment[FIRST] = DeviceAlgorithms::ScanInclusive(pow1Array, pow1Array) / numValues;
|
|
|
|
statinfo.rawMoment[SECOND] = DeviceAlgorithms::ScanInclusive(pow2Array, pow2Array) / numValues;
|
|
|
|
statinfo.rawMoment[THIRD] = DeviceAlgorithms::ScanInclusive(pow3Array, pow3Array) / numValues;
|
|
|
|
statinfo.rawMoment[FOURTH] = DeviceAlgorithms::ScanInclusive(pow4Array, pow4Array) / numValues;
|
|
|
|
|
|
|
|
// Subtract the mean from every value and leave in tempArray
|
2018-08-28 20:36:50 +00:00
|
|
|
vtkm::worklet::DispatcherMapField<SubtractConst> subtractConstDispatcher(
|
2017-05-18 14:29:41 +00:00
|
|
|
SubtractConst(statinfo.mean));
|
2015-09-10 17:13:08 +00:00
|
|
|
subtractConstDispatcher.Invoke(fieldArray, tempArray);
|
|
|
|
|
|
|
|
// Calculate sums of powers on the (value - mean) array
|
|
|
|
calculatePowersDispatcher.Invoke(tempArray, pow1Array, pow2Array, pow3Array, pow4Array);
|
|
|
|
|
|
|
|
// Accumulate the results using ScanInclusive
|
2017-05-18 14:29:41 +00:00
|
|
|
statinfo.centralMoment[FIRST] =
|
|
|
|
DeviceAlgorithms::ScanInclusive(pow1Array, pow1Array) / numValues;
|
|
|
|
statinfo.centralMoment[SECOND] =
|
|
|
|
DeviceAlgorithms::ScanInclusive(pow2Array, pow2Array) / numValues;
|
|
|
|
statinfo.centralMoment[THIRD] =
|
|
|
|
DeviceAlgorithms::ScanInclusive(pow3Array, pow3Array) / numValues;
|
|
|
|
statinfo.centralMoment[FOURTH] =
|
|
|
|
DeviceAlgorithms::ScanInclusive(pow4Array, pow4Array) / numValues;
|
2015-09-10 17:13:08 +00:00
|
|
|
|
|
|
|
// Statistics from the moments
|
|
|
|
statinfo.variance = statinfo.centralMoment[SECOND];
|
|
|
|
statinfo.stddev = Sqrt(statinfo.variance);
|
2017-05-18 14:29:41 +00:00
|
|
|
statinfo.skewness =
|
|
|
|
statinfo.centralMoment[THIRD] / Pow(statinfo.stddev, static_cast<FieldType>(3.0));
|
|
|
|
statinfo.kurtosis =
|
|
|
|
statinfo.centralMoment[FOURTH] / Pow(statinfo.stddev, static_cast<FieldType>(4.0));
|
2016-06-27 19:39:03 +00:00
|
|
|
}
|
2015-09-10 17:13:08 +00:00
|
|
|
};
|
|
|
|
}
|
|
|
|
} // namespace vtkm::worklet
|
|
|
|
|
|
|
|
#endif // vtk_m_worklet_FieldStatistics_h
|