vtkh stats to vtkm stats (#1150)

replace vtkh stats worklet with vtkm stats

src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp

@@ -3059,16 +3059,19 @@ 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());
     MPI_Comm_rank(mpi_comm, &rank);
 #endif
     if(rank == 0)
     {
-      res.Print(std::cout);
+      res->PrintSummary(std::cout);
     }
 }
 //-----------------------------------------------------------------------------

src/libs/vtkh/filters/Statistics.cpp

@@ -1,11 +1,11 @@
 #include <vtkh/filters/Statistics.hpp>
+#include <vtkh/vtkm_filters/vtkmStatistics.hpp>
 #include <vtkh/Error.hpp>
 #include <vtkh/Logger.hpp>
 #include <vtkm/cont/Algorithm.h>
 #include <vtkm/cont/ArrayHandleCast.h>
 #include <vtkm/cont/Invoker.h>
-#include <vtkm/worklet/DispatcherMapField.h>
-#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/cont/PartitionedDataSet.h>
 #include <vector>
 
 #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 <typename T>
-  VTKM_EXEC void operator()(const T &in, vtkm::Float32 &out) const
-  {
-    out = static_cast<vtkm::Float32>(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 <typename T>
-  VTKM_EXEC vtkm::Float32 operator()(T value) const
-  {
-    return static_cast<vtkm::Float32>(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 <typename T>
-  VTKM_EXEC vtkm::Float32 operator()(T value) const
-  {
-    return static_cast<vtkm::Float32>(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 <typename T>
-  VTKM_EXEC vtkm::Float32 operator()(T value) const
-  {
-    return static_cast<vtkm::Float32>(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<vtkm::cont::ArrayHandle<vtkm::Float32>> fields;
-  int total_values = 0;
 
+  std::vector<vtkm::cont::DataSet> 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<vtkm::Float32> float_field;
-      vtkm::worklet::DispatcherMapField<detail::CopyToFloat>(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<vtkm::cont::ArrayHandle<vtkm::Float32>> unbiased_fields(fields.size());
-  std::vector<vtkm::cont::ArrayHandle<vtkm::Float32>> 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

src/libs/vtkh/filters/Statistics.hpp

@@ -4,32 +4,26 @@
 #include <vtkh/vtkh_exports.h>
 #include <vtkh/vtkh.hpp>
 #include <vtkh/DataSet.hpp>
+#include <vtkh/filters/Filter.hpp>
 
 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    : "<<mean<<"\n";
-      out<<"Variance: "<<variance<<"\n";
-      out<<"Skewness: "<<skewness<<"\n";
-      out<<"Kurtosis: "<<kurtosis<<"\n";
-    }
-  };
+  void SetField(const std::string &field_name);
+  std::string GetField() const;
+protected:
+  void PreExecute() override;
+  void PostExecute() override;
+  void DoExecute() override;
 
-  Statistics();
-  ~Statistics();
-  Statistics::Result Run(vtkh::DataSet &data_set, const std::string field_name);
+  std::string m_field_name;
 
 };
 

src/libs/vtkh/vtkm_filters/CMakeLists.txt

@@ -18,6 +18,7 @@ set(vtkm_filters_headers
     vtkmMeshQuality.hpp
     vtkmPointAverage.hpp
     vtkmPointTransform.hpp
+    vtkmStatistics.hpp
     vtkmTetrahedralize.hpp
     vtkmTriangulate.hpp
     vtkmThreshold.hpp
@@ -41,6 +42,7 @@ set(vtkm_filters_sources
     vtkmMeshQuality.cpp
     vtkmPointAverage.cpp
     vtkmPointTransform.cpp
+    vtkmStatistics.cpp
     vtkmTetrahedralize.cpp
     vtkmTriangulate.cpp
     vtkmThreshold.cpp

src/libs/vtkh/vtkm_filters/vtkmStatistics.cpp

@@ -0,0 +1,19 @@
+#include "vtkmStatistics.hpp"
+#include <vtkm/filter/density_estimate/Statistics.h>
+
+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

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 <vtkm/cont/PartitionedDataSet.h>
 #include <vtkm/cont/DataSet.h>
-#include <vtkm/filter/density_estimation/Statistics.h>
+#include <vtkm/filter/density_estimate/Statistics.h>
 
 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