diff --git a/docs/changelog/statistics-filter.md b/docs/changelog/statistics-filter.md new file mode 100644 index 000000000..d20cfa74d --- /dev/null +++ b/docs/changelog/statistics-filter.md @@ -0,0 +1,3 @@ +# New Statistics filter + +The statistics filter computes the descriptive statistics of the fields specified by users based on `DescriptiveStatistics`. Users can set `RequiredStatsList` to specify which statistics will be stored in the output data set. diff --git a/vtkm/filter/density_estimate/CMakeLists.txt b/vtkm/filter/density_estimate/CMakeLists.txt index bac3f2a06..2fb862448 100644 --- a/vtkm/filter/density_estimate/CMakeLists.txt +++ b/vtkm/filter/density_estimate/CMakeLists.txt @@ -15,6 +15,7 @@ set(density_estimate_headers ParticleDensityBase.h ParticleDensityCloudInCell.h ParticleDensityNearestGridPoint.h + Statistics.h ) set(density_estimate_sources_device @@ -25,6 +26,7 @@ set(density_estimate_sources_device ParticleDensityBase.cxx ParticleDensityCloudInCell.cxx ParticleDensityNearestGridPoint.cxx + Statistics.cxx ) vtkm_library( diff --git a/vtkm/filter/density_estimate/Statistics.cxx b/vtkm/filter/density_estimate/Statistics.cxx new file mode 100644 index 000000000..2b9a459f8 --- /dev/null +++ b/vtkm/filter/density_estimate/Statistics.cxx @@ -0,0 +1,114 @@ +//============================================================================ +// 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. +//============================================================================ +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace density_estimate +{ + +VTKM_CONT vtkm::cont::DataSet Statistics::DoExecute(const vtkm::cont::DataSet& inData) +{ + vtkm::worklet::DescriptiveStatistics worklet; + vtkm::cont::DataSet output; + + auto resolveType = [&](const auto& concrete) { + auto result = worklet.Run(concrete); + + for (size_t i = 0; i < RequiredStatsList.size(); i++) + { + vtkm::cont::ArrayHandle stat; + stat.Allocate(1); + Stats statEnum = RequiredStatsList[i]; + + switch (statEnum) + { + case Stats::N: + { + stat.WritePortal().Set(0, static_cast(result.N())); + break; + } + case Stats::Min: + { + stat.WritePortal().Set(0, static_cast(result.Min())); + break; + } + case Stats::Max: + { + stat.WritePortal().Set(0, static_cast(result.Max())); + break; + } + case Stats::Sum: + { + stat.WritePortal().Set(0, static_cast(result.Sum())); + break; + } + case Stats::Mean: + { + stat.WritePortal().Set(0, static_cast(result.Mean())); + break; + } + case Stats::SampleStdDev: + { + stat.WritePortal().Set(0, static_cast(result.SampleStddev())); + break; + } + case Stats::PopulationStdDev: + { + stat.WritePortal().Set(0, static_cast(result.PopulationStddev())); + break; + } + case Stats::SampleVariance: + { + stat.WritePortal().Set(0, static_cast(result.SampleVariance())); + break; + } + case Stats::PopulationVariance: + { + stat.WritePortal().Set(0, static_cast(result.PopulationVariance())); + break; + } + case Stats::Skewness: + { + stat.WritePortal().Set(0, static_cast(result.Skewness())); + break; + } + case Stats::Kurtosis: + { + stat.WritePortal().Set(0, static_cast(result.Kurtosis())); + break; + } + default: + { + throw vtkm::cont::ErrorFilterExecution( + "Unsupported statistics variable in statistics filter."); + } + } + + output.AddField({ this->StatsName[static_cast(statEnum)], + vtkm::cont::Field::Association::WholeDataSet, + stat }); + } + }; + const auto& fieldArray = this->GetFieldFromDataSet(inData).GetData(); + fieldArray + .CastAndCallForTypesWithFloatFallback( + resolveType); + + return output; +} + +} // namespace density_estimate +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/density_estimate/Statistics.h b/vtkm/filter/density_estimate/Statistics.h new file mode 100644 index 000000000..769fcb00f --- /dev/null +++ b/vtkm/filter/density_estimate/Statistics.h @@ -0,0 +1,81 @@ +//============================================================================ +// 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. +//============================================================================ +#ifndef vtk_m_filter_density_estimate_Statistics_h +#define vtk_m_filter_density_estimate_Statistics_h + +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace density_estimate +{ +/// \brief Calculating the statistics of input fields +/// +/// This filter calculates the statistics of input fields. +/// +class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT Statistics : public vtkm::filter::FilterField +{ +public: + enum struct Stats + { + N = 0, + Min, + Max, + Sum, + Mean, + SampleStdDev, + PopulationStdDev, + SampleVariance, + PopulationVariance, + Skewness, + Kurtosis + }; + + + /// \{ + /// \brief The output statistical variables for executing the statistics filter. + /// + void SetRequiredStats(const std::vector StatsList) { RequiredStatsList = StatsList; } + const std::vector& GetRequiredStats() const { return this->RequiredStatsList; } + /// \} +private: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override; + std::vector RequiredStatsList{ Stats::N, + Stats::Min, + Stats::Max, + Stats::Sum, + Stats::Mean, + Stats::SampleStdDev, + Stats::PopulationStdDev, + Stats::SampleVariance, + Stats::PopulationVariance, + Stats::Skewness, + Stats::Kurtosis }; + // This string vector stores variables names stored in the output dataset + std::vector StatsName{ "N", + "Min", + "Max", + "Sum", + "Mean", + "SampleStddev", + "PopulationStdDev", + "SampleVariance", + "PopulationVariance", + "Skewness", + "Kurtosis" }; +}; +} // namespace density_estimate +} // namespace filter +} // namespace vtkm + +#endif //vtk_m_filter_density_estimate_Statistics_h diff --git a/vtkm/filter/density_estimate/testing/CMakeLists.txt b/vtkm/filter/density_estimate/testing/CMakeLists.txt index 7d28e5409..d36470d03 100644 --- a/vtkm/filter/density_estimate/testing/CMakeLists.txt +++ b/vtkm/filter/density_estimate/testing/CMakeLists.txt @@ -14,6 +14,7 @@ set(unit_tests UnitTestNDEntropyFilter.cxx UnitTestNDHistogramFilter.cxx UnitTestPartitionedDataSetHistogramFilter.cxx + UnitTestStatisticsFilter.cxx ) set(unit_tests_device diff --git a/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx b/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx new file mode 100644 index 000000000..15dade562 --- /dev/null +++ b/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx @@ -0,0 +1,86 @@ +//============================================================================ +// 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. +//============================================================================ + +#include + +#include +#include + + +namespace +{ + +vtkm::FloatDefault getStatsFromArray(vtkm::cont::DataSet dataset, const std::string statName) +{ + vtkm::cont::ArrayHandle array; + dataset.GetField(statName).GetData().AsArrayHandle(array); + vtkm::cont::ArrayHandle::ReadPortalType portal = array.ReadPortal(); + vtkm::FloatDefault value = portal.Get(0); + return value; +} + +void TestStatisticsPartial() +{ + vtkm::cont::DataSet dataSet; + constexpr vtkm::FloatDefault N = 1000; + // the number is from 0 to 999 + auto scalaArray = + vtkm::cont::ArrayHandleCounting(0.0f, 1.0f, static_cast(N)); + dataSet.AddPointField("scalarField", scalaArray); + + using STATS = vtkm::filter::density_estimate::Statistics; + STATS statisticsFilter; + + //set required states + std::vector RequiredStatsList{ STATS::Stats::N, STATS::Stats::Sum, + STATS::Stats::Mean, STATS::Stats::SampleVariance, + STATS::Stats::Skewness, STATS::Stats::Kurtosis }; + + //default RequiredStatsList contains all statistics variables + statisticsFilter.SetRequiredStats(RequiredStatsList); + using AsscoType = vtkm::cont::Field::Association; + + statisticsFilter.SetActiveField("scalarField", AsscoType::Points); + + // We use the same test cases with the UnitTestDescriptiveStatistics.h + vtkm::cont::DataSet resultDataSet = statisticsFilter.Execute(dataSet); + + vtkm::FloatDefault NValueFromFilter = getStatsFromArray(resultDataSet, "N"); + VTKM_TEST_ASSERT(test_equal(NValueFromFilter, N)); + + vtkm::FloatDefault SumFromFilter = getStatsFromArray(resultDataSet, "Sum"); + VTKM_TEST_ASSERT(test_equal(SumFromFilter, N * (N - 1) / 2)); + + vtkm::FloatDefault MeanFromFilter = getStatsFromArray(resultDataSet, "Mean"); + VTKM_TEST_ASSERT(test_equal(MeanFromFilter, (N - 1) / 2)); + + vtkm::FloatDefault SVFromFilter = getStatsFromArray(resultDataSet, "SampleVariance"); + VTKM_TEST_ASSERT(test_equal(SVFromFilter, 83416.66)); + + vtkm::FloatDefault SkewnessFromFilter = getStatsFromArray(resultDataSet, "Skewness"); + VTKM_TEST_ASSERT(test_equal(SkewnessFromFilter, 0)); + + // we use fisher=False when computing the Kurtosis value + vtkm::FloatDefault KurtosisFromFilter = getStatsFromArray(resultDataSet, "Kurtosis"); + VTKM_TEST_ASSERT(test_equal(KurtosisFromFilter, 1.8)); +} + + +void TestStatistics() +{ + TestStatisticsPartial(); +} // TestFieldStatistics +} + +//More deatiled tests can be found in the UnitTestStatisticsFilter +int UnitTestStatisticsFilter(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestStatistics, argc, argv); +} diff --git a/vtkm/worklet/FieldStatistics.h b/vtkm/worklet/FieldStatistics.h index 8e02f37ac..701e7bcb0 100644 --- a/vtkm/worklet/FieldStatistics.h +++ b/vtkm/worklet/FieldStatistics.h @@ -27,10 +27,9 @@ namespace vtkm { namespace worklet { - //simple functor that prints basic statistics template -class FieldStatistics +class VTKM_DEPRECATED(2.1, "Use DescriptiveStatistics or the statistics filter.") FieldStatistics { public: // For moments readability diff --git a/vtkm/worklet/testing/UnitTestFieldStatistics.cxx b/vtkm/worklet/testing/UnitTestFieldStatistics.cxx index 0b238dce2..078cdb993 100644 --- a/vtkm/worklet/testing/UnitTestFieldStatistics.cxx +++ b/vtkm/worklet/testing/UnitTestFieldStatistics.cxx @@ -283,6 +283,7 @@ vtkm::cont::DataSet Make2DUniformStatDataSet1() // // Create a dataset with known point data and cell data (statistical distributions) // +VTKM_DEPRECATED_SUPPRESS_BEGIN void PrintStatInfo(vtkm::worklet::FieldStatistics::StatInfo statinfo) { std::cout << " Median " << statinfo.median << std::endl; @@ -302,12 +303,14 @@ void PrintStatInfo(vtkm::worklet::FieldStatistics::StatInfo stati std::cout << statinfo.centralMoment[i] << " "; std::cout << "]" << std::endl; } +VTKM_DEPRECATED_SUPPRESS_END // // Test simple dataset 2D Uniform with 10 cells // void TestFieldSimple() { + VTKM_DEPRECATED_SUPPRESS_BEGIN // Create the output structure vtkm::worklet::FieldStatistics::StatInfo statinfo; @@ -332,6 +335,7 @@ void TestFieldSimple() VTKM_TEST_ASSERT(test_equal(statinfo.stddev, expected[5]), "Error in stddev"); VTKM_TEST_ASSERT(test_equal(statinfo.skewness, expected[6]), "Error in skewness"); VTKM_TEST_ASSERT(test_equal(statinfo.kurtosis, expected[7]), "Error in kurtosis"); + VTKM_DEPRECATED_SUPPRESS_END } // @@ -339,6 +343,7 @@ void TestFieldSimple() // void TestFieldStandardDistributions() { + VTKM_DEPRECATED_SUPPRESS_BEGIN // Create the output structure vtkm::worklet::FieldStatistics::StatInfo statinfo; @@ -374,6 +379,8 @@ void TestFieldStandardDistributions() vtkm::worklet::FieldStatistics().Run(p_uniform, statinfo); std::cout << "Uniform distributed POINT data:" << std::endl; PrintStatInfo(statinfo); + VTKM_DEPRECATED_SUPPRESS_END + } // TestFieldStatistics void TestFieldStatistics()