From 8baa78c9f2f32f7e170d3281ee7168226bc3cb41 Mon Sep 17 00:00:00 2001 From: Nicole Date: Fri, 3 May 2024 07:54:13 -0700 Subject: [PATCH] vtkh stats to vtkm stats (#1150) replace vtkh stats worklet with vtkm stats --- .../ascent_runtime_vtkh_filters.cpp | 7 +- src/libs/vtkh/filters/Statistics.cpp | 211 +++++------------- src/libs/vtkh/filters/Statistics.hpp | 30 +-- src/libs/vtkh/vtkm_filters/CMakeLists.txt | 2 + src/libs/vtkh/vtkm_filters/vtkmStatistics.cpp | 19 ++ src/libs/vtkh/vtkm_filters/vtkmStatistics.hpp | 15 +- src/tests/vtkh/t_vtk-h_statistics_par.cpp | 10 +- 7 files changed, 105 insertions(+), 189 deletions(-) create mode 100644 src/libs/vtkh/vtkm_filters/vtkmStatistics.cpp diff --git a/src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp b/src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp index a9ff6e95e..5c8976673 100644 --- a/src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp +++ b/src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp @@ -3059,8 +3059,11 @@ VTKHStats::execute() vtkh::DataSet &data = collection->dataset_by_topology(topo_name); vtkh::Statistics stats; + stats.SetField(field_name); + stats.SetInput(&data); + stats.Update(); - vtkh::Statistics::Result res = stats.Run(data, field_name); + vtkh::DataSet* res = stats.GetOutput(); int rank = 0; #ifdef ASCENT_MPI_ENABLED MPI_Comm mpi_comm = MPI_Comm_f2c(Workspace::default_mpi_comm()); @@ -3068,7 +3071,7 @@ VTKHStats::execute() #endif if(rank == 0) { - res.Print(std::cout); + res->PrintSummary(std::cout); } } //----------------------------------------------------------------------------- diff --git a/src/libs/vtkh/filters/Statistics.cpp b/src/libs/vtkh/filters/Statistics.cpp index 1bcb93df4..64efe0bf5 100644 --- a/src/libs/vtkh/filters/Statistics.cpp +++ b/src/libs/vtkh/filters/Statistics.cpp @@ -1,11 +1,11 @@ #include +#include #include #include #include #include #include -#include -#include +#include #include #ifdef VTKH_PARALLEL @@ -18,199 +18,94 @@ namespace vtkh namespace detail { -class CopyToFloat : public vtkm::worklet::WorkletMapField -{ -public: - typedef void ControlSignature(FieldIn, FieldOut); - typedef void ExecutionSignature(_1, _2); - - VTKM_CONT CopyToFloat() - { } - - template - VTKM_EXEC void operator()(const T &in, vtkm::Float32 &out) const - { - out = static_cast(in); - } -}; +} // namespace detail -class SubtractMean : public vtkm::worklet::WorkletMapField +Statistics::Statistics() { -public: - // typedef void ControlSignature() - using ControlSignature = void(FieldIn, FieldOut); - using ExecutionSignature = _2(_1); - using InputDomain = _1; - - vtkm::Float32 mean; - // define constructor to set mean as the input mean - VTKM_CONT SubtractMean(vtkm::Float32 mean) : mean(mean) { } +} - template - VTKM_EXEC vtkm::Float32 operator()(T value) const - { - return static_cast(value - this->mean); - } -}; -// create worklet class to exponentiate vector elements by power -// Create a worklet -class VecPow : public vtkm::worklet::WorkletMapField +Statistics::~Statistics() { -public: - // take each element of the vector to the power p - // typedef void ControlSignature() - using ControlSignature = void(FieldIn, FieldOut); - using ExecutionSignature = _2(_1); - using InputDomain = _1; - int power; - - // define constructor to set mean as the input mean - VTKM_CONT VecPow(int p) : power(p) { } - - template - VTKM_EXEC vtkm::Float32 operator()(T value) const - { - return static_cast(vtkm::Pow(value,this->power)); - } -}; +} -// worklet class to divide by scalar -class DivideByScalar : public vtkm::worklet::WorkletMapField +void +Statistics::SetField(const std::string &field_name) { -public: - // typedef void ControlSignature() - using ControlSignature = void(FieldIn, FieldOut); - using ExecutionSignature = _2(_1); - using InputDomain = _1; - - vtkm::Float32 scale; - - // define constructor to set scale as the input scale - VTKM_CONT DivideByScalar(vtkm::Float32 s) : scale(s) { } - - template - VTKM_EXEC vtkm::Float32 operator()(T value) const - { - return static_cast(value/this->scale); - } -}; - -} // namespace detail + m_field_name = field_name; +} -Statistics::Statistics() +std::string +Statistics::GetField() const { - + return m_field_name; } -Statistics::~Statistics() +void +Statistics::PreExecute() { + Filter::PreExecute(); +} +void +Statistics::PostExecute() +{ + Filter::PostExecute(); } -Statistics::Result Statistics::Run(vtkh::DataSet &data_set, const std::string field_name) +void Statistics::DoExecute() { VTKH_DATA_OPEN("statistics"); VTKH_DATA_ADD("device", GetCurrentDevice()); - VTKH_DATA_ADD("input_cells", data_set.GetNumberOfCells()); - VTKH_DATA_ADD("input_domains", data_set.GetNumberOfDomains()); - const int num_domains = data_set.GetNumberOfDomains(); + VTKH_DATA_ADD("input_cells", this->m_input->GetNumberOfCells()); + VTKH_DATA_ADD("input_domains", this->m_input->GetNumberOfDomains()); + const int num_domains = this->m_input->GetNumberOfDomains(); + this->m_output = new DataSet(); - if(!data_set.GlobalFieldExists(field_name)) + if(!this->m_input->GlobalFieldExists(m_field_name)) { - throw Error("Statistics: field : '"+field_name+"' does not exist'"); + throw Error("Statistics: field : '"+m_field_name+"' does not exist'"); } - std::vector> fields; - int total_values = 0; + std::vector vtkm_ds; + + for(int i = 0; i < num_domains; ++i) { vtkm::Id domain_id; vtkm::cont::DataSet dom; - data_set.GetDomain(i, dom, domain_id); - if(dom.HasField(field_name)) + this->m_input->GetDomain(i, dom, domain_id); + if(dom.HasField(m_field_name)) { - vtkm::cont::Field field = dom.GetField(field_name); - vtkm::cont::ArrayHandle float_field; - vtkm::worklet::DispatcherMapField(detail::CopyToFloat()) - .Invoke(field.GetData().ResetTypes(vtkm::TypeListFieldScalar(),VTKM_DEFAULT_STORAGE_LIST{}), float_field); - fields.push_back(float_field); - total_values += float_field.GetNumberOfValues(); + vtkm_ds.push_back(dom); } } - Statistics::Result res; -#ifdef VTKH_PARALLEL - MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle()); - MPI_Allreduce(MPI_IN_PLACE, &total_values, 1, MPI_INT, MPI_SUM, mpi_comm); -#endif + vtkm::cont::PartitionedDataSet data_pds(vtkm_ds); + vtkmStatistics stats; + auto result = stats.Run(data_pds, m_field_name); - // mean - vtkm::Float32 sum = 0; - for(size_t i = 0; i < fields.size(); ++i) + int size = result.GetNumberOfFields(); + vtkm::cont::DataSet dom; + + for(int i = 0; i < size; i++) { - sum += vtkm::cont::Algorithm::Reduce(fields[i],0.0f,vtkm::Sum()); + //g_field will have assoc=Global which only goes with vtkm::PDS + //convert to new field with assoc=WholeDataSet to put in vtkm::DS + vtkm::cont::Field g_field = result.GetField(i); + vtkm::cont::Field field(g_field.GetName(),vtkm::cont::Field::Association::WholeDataSet,g_field.GetData()); + dom.AddField(field); } - -#ifdef VTKH_PARALLEL - MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_FLOAT, MPI_SUM, mpi_comm); -#endif - vtkm::Float32 mean = sum / vtkm::Float32(total_values); - res.mean = mean; - - std::vector> unbiased_fields(fields.size()); - std::vector> temp_fields(fields.size()); - - vtkm::cont::Invoker invoke; - - // variance - vtkm::Float32 sum2 = 0; - for(size_t i = 0; i < fields.size(); ++i) - { - invoke(detail::SubtractMean(mean), fields[i], unbiased_fields[i]); - invoke(detail::VecPow(2),unbiased_fields[i],temp_fields[i]); - sum2 += vtkm::cont::Algorithm::Reduce(temp_fields[i],0.0f,vtkm::Sum()); - } - -#ifdef VTKH_PARALLEL - MPI_Allreduce(MPI_IN_PLACE, &sum2, 1, MPI_FLOAT, MPI_SUM, mpi_comm); -#endif - - vtkm::Float32 variance = sum2 / vtkm::Float32(total_values - 1); - res.variance = variance; - - // skewness - vtkm::Float32 sum3 = 0; - for(size_t i = 0; i < fields.size(); ++i) - { - invoke(detail::VecPow(3),unbiased_fields[i],temp_fields[i]); - sum3 += vtkm::cont::Algorithm::Reduce(temp_fields[i],0.0f,vtkm::Sum()); - } - -#ifdef VTKH_PARALLEL - MPI_Allreduce(MPI_IN_PLACE, &sum3, 1, MPI_FLOAT, MPI_SUM, mpi_comm); -#endif - vtkm::Float32 skewness = (sum3/vtkm::Float32(total_values))/(vtkm::Pow(variance,1.5f)); - res.skewness = skewness; - - // kertosis - vtkm::Float32 sum4 = 0; - for(size_t i = 0; i < fields.size(); ++i) - { - invoke(detail::VecPow(4),unbiased_fields[i],temp_fields[i]); - sum4 += vtkm::cont::Algorithm::Reduce(temp_fields[i],0.0f,vtkm::Sum()); - } - -#ifdef VTKH_PARALLEL - MPI_Allreduce(MPI_IN_PLACE, &sum4, 1, MPI_FLOAT, MPI_SUM, mpi_comm); -#endif - vtkm::Float32 kurtosis = (sum4/vtkm::Float32(total_values))/vtkm::Pow(variance,2.f) - 3.0f; - res.kurtosis = kurtosis; + this->m_output->AddDomain(dom,0); VTKH_DATA_CLOSE(); - return res; } +std::string +Statistics::GetName() const +{ + return "vtkh::Statistics"; +} } // namespace vtkh diff --git a/src/libs/vtkh/filters/Statistics.hpp b/src/libs/vtkh/filters/Statistics.hpp index 80b63c36d..a6a91d164 100644 --- a/src/libs/vtkh/filters/Statistics.hpp +++ b/src/libs/vtkh/filters/Statistics.hpp @@ -4,32 +4,26 @@ #include #include #include +#include namespace vtkh { -class VTKH_API Statistics +class VTKH_API Statistics: public Filter { public: + Statistics(); + virtual ~Statistics(); + std::string GetName() const override; - struct Result - { - vtkm::Float32 mean; - vtkm::Float32 variance; - vtkm::Float32 skewness; - vtkm::Float32 kurtosis; - void Print(std::ostream &out) - { - out<<"Mean : "< + +namespace vtkh +{ +vtkm::cont::PartitionedDataSet +vtkmStatistics::Run(vtkm::cont::PartitionedDataSet &p_input, + std::string field_name) +{ + vtkm::filter::density_estimate::Statistics stats; + + stats.SetActiveField(field_name); + + auto output = stats.Execute(p_input); + //output.PrintSummary(std::cerr); + return output; +} + +} // namespace vtkh diff --git a/src/libs/vtkh/vtkm_filters/vtkmStatistics.hpp b/src/libs/vtkh/vtkm_filters/vtkmStatistics.hpp index ffd490bd3..8eb6d5f9a 100644 --- a/src/libs/vtkh/vtkm_filters/vtkmStatistics.hpp +++ b/src/libs/vtkh/vtkm_filters/vtkmStatistics.hpp @@ -1,19 +1,18 @@ -#ifndef VTK_H_VTKM_POINT_AVERAGE_HPP -#define VTK_H_VTKM_POINT_AVERAGE_HPP +#ifndef VTK_H_VTKM_STATISTICS_HPP +#define VTK_H_VTKM_STATISTICS_HPP +#include #include -#include +#include namespace vtkh { -class vtkmPointAverage +class vtkmStatistics { public: - vtkm::cont::DataSet Run(vtkm::cont::DataSet &input, - std::string field_name, - std::string output_field_name, - vtkm::filter::FieldSelection map_fields); + vtkm::cont::PartitionedDataSet Run(vtkm::cont::PartitionedDataSet &p_input, + std::string field_name); }; } #endif diff --git a/src/tests/vtkh/t_vtk-h_statistics_par.cpp b/src/tests/vtkh/t_vtk-h_statistics_par.cpp index 70540454d..7a139e1af 100644 --- a/src/tests/vtkh/t_vtk-h_statistics_par.cpp +++ b/src/tests/vtkh/t_vtk-h_statistics_par.cpp @@ -39,12 +39,16 @@ TEST(vtkh_statistics_par, vtkh_stats) data_set.AddDomain(CreateTestData(domain_id, num_blocks, base_size), domain_id); } - vtkh::Statistics::Result res; + vtkh::DataSet* res; vtkh::Statistics stats; - res = stats.Run(data_set,"point_data_Float64"); + stats.SetField("point_data_Float64"); + stats.SetInput(&data_set); + stats.Update(); + res = stats.GetOutput(); - if(rank == 0) res.Print(std::cout); + + if(rank == 0) res->PrintSummary(std::cout); MPI_Finalize(); }