Commit fe624fdf authored by Utkarsh Ayachit's avatar Utkarsh Ayachit
Browse files

Add ability to generate histograms for non-generic arrays VTK.

parent 3f925775
......@@ -8,16 +8,12 @@ set (sources
ConfigurableAnalysis.h
DataAdaptor.cxx
DataAdaptor.h
Histogram.cxx
Histogram.h
PosthocIO.cxx
PosthocIO.h
)
if(VTK_HAS_GENERIC_ARRAYS)
list(APPEND sources
Histogram.cxx
Histogram.h)
endif()
if(ENABLE_CATALYST)
list(APPEND sources
catalyst/AnalysisAdaptor.cxx
......@@ -64,7 +60,7 @@ else()
endif()
if(VTK_HAS_GENERIC_ARRAYS)
target_compile_definitions(sensei PRIVATE ENABLE_HISTOGRAM)
target_compile_definitions(sensei PRIVATE VTK_HAS_GENERIC_ARRAYS)
endif()
if(ENABLE_ADIOS)
......
......@@ -7,9 +7,7 @@
#include "Autocorrelation.h"
#include "PosthocIO.h"
#ifdef ENABLE_HISTOGRAM
# include "Histogram.h"
#endif
#include "Histogram.h"
#ifdef ENABLE_ADIOS
# include "adios/AnalysisAdaptor.h"
#endif
......@@ -112,30 +110,25 @@ class ConfigurableAnalysis::vtkInternals
public:
std::vector<vtkSmartPointer<AnalysisAdaptor> > Analyses;
#ifdef ENABLE_HISTOGRAM
void AddHistogram(MPI_Comm comm, pugi::xml_node node)
int AddHistogram(MPI_Comm comm, pugi::xml_node node)
{
if (node.attribute("enabled") && !node.attribute("enabled").as_int())
return -1;
pugi::xml_attribute array = node.attribute("array");
if (array)
{
int association = GetAssociation(node.attribute("association"));
int bins = node.attribute("bins")? node.attribute("bins").as_int() : 10;
vtkNew<Histogram> histogram;
histogram->Initialize(comm, bins, association, array.value());
this->Analyses.push_back(histogram.GetPointer());
}
else
if (!array)
{
ConfigurableAnalysisError(<< "'histogram' missing required attribute 'array'. Skipping.");
return -1;
}
int association = GetAssociation(node.attribute("association"));
int bins = node.attribute("bins")? node.attribute("bins").as_int() : 10;
vtkNew<Histogram> histogram;
histogram->Initialize(comm, bins, association, array.value());
this->Analyses.push_back(histogram.GetPointer());
return 0;
}
#endif
#ifdef ENABLE_ADIOS
int AddAdios(MPI_Comm comm, pugi::xml_node node)
......@@ -308,11 +301,9 @@ bool ConfigurableAnalysis::Initialize(MPI_Comm comm, const std::string& filename
analysis; analysis = analysis.next_sibling("analysis"))
{
std::string type = analysis.attribute("type").value();
#ifdef ENABLE_HISTOGRAM
if ((type == "histogram") &&
!this->Internals->AddHistogram(comm, analysis))
continue;
#endif
#ifdef ENABLE_ADIOS
if ((type == "adios") &&
!this->Internals->AddAdios(comm, analysis))
......
#include "Histogram.h"
#include <vtkAOSDataArrayTemplate.h>
#include <vtkArrayDispatch.h>
#include <vtkCompositeDataIterator.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataArray.h>
......@@ -10,123 +8,18 @@
#include <vtkObjectFactory.h>
#include <vtkSmartPointer.h>
#ifdef VTK_HAS_GENERIC_ARRAYS
# include "HistogramInternals-GenericArrays.h"
#else
# include "HistogramInternals-NoGenericArrays.h"
#endif
#include "DataAdaptor.h"
#include <timer/Timer.h>
#include <algorithm>
#include <vector>
namespace
{
// Private worker for histogram method. Computes the local histogram on
// array (passed to operator()). To be used with vtkArrayDispatch.
//
// Inputs:
// range: Global range of data
// bins: Number of histogram bins
// array: Local data.
//
// Outputs:
// Histogram: The histogram of the local data.
struct HistogramWorker
{
const double *Range;
int Bins;
std::vector<unsigned int> Histogram;
HistogramWorker(const double *range, int bins) : Range(range), Bins(bins) {}
template <typename ArrayT>
void operator()(ArrayT *array)
{
assert(array);
assert(array->GetNumberOfComponents() == 1);
typedef typename ArrayT::ValueType ValueType;
ValueType width = static_cast<ValueType>((this->Range[1] - this->Range[0]) / this->Bins);
ValueType min = static_cast<ValueType>(this->Range[0]);
vtkIdType numTuples = array->GetNumberOfTuples();
// + 1 to store val == max. These will be moved to this last bin before
// returning (Avoids having to branch in the loop below);
this->Histogram.resize(this->Bins + 1, 0);
for (vtkIdType tIdx = 0; tIdx < numTuples; ++tIdx)
{
int bin = static_cast<int>((array->GetComponent(tIdx, 0) - min) / width);
++this->Histogram[bin];
}
// Merge the last two bins (the last is only when val == max)
this->Histogram[this->Bins-1] += this->Histogram[this->Bins];
this->Histogram.resize(this->Bins);
}
};
class vtkHistogram
{
double Range[2];
HistogramWorker *Worker;
public:
vtkHistogram()
{
this->Range[0] = VTK_DOUBLE_MAX;
this->Range[1] = VTK_DOUBLE_MIN;
this->Worker = NULL;
}
~vtkHistogram()
{
delete this->Worker;
}
void AddRange(vtkDataArray* da)
{
if (da)
{
double crange[2];
da->GetRange(crange);
this->Range[0] = std::min(this->Range[0], crange[0]);
this->Range[1] = std::max(this->Range[1], crange[1]);
}
}
void PreCompute(MPI_Comm comm, int bins)
{
double g_range[2];
// Find the global max/min
MPI_Allreduce(&this->Range[0], &g_range[0], 1, MPI_DOUBLE, MPI_MIN, comm);
MPI_Allreduce(&this->Range[1], &g_range[1], 1, MPI_DOUBLE, MPI_MAX, comm);
this->Range[0] = g_range[0];
this->Range[1] = g_range[1];
this->Worker = new HistogramWorker(this->Range, bins);
}
void Compute(vtkDataArray* da)
{
if (da)
{
vtkArrayDispatch::Dispatch::Execute(da, *this->Worker);
}
}
void PostCompute(MPI_Comm comm, int bins, const std::string& name)
{
std::vector<unsigned int> gHist(bins, 0);
MPI_Reduce(&this->Worker->Histogram[0], &gHist[0], bins, MPI_UNSIGNED, MPI_SUM, 0, comm);
int rank;
MPI_Comm_rank(comm, &rank);
if (rank == 0)
{
std::cout << "Histogram '" << name << "' (VTK):\n";
double width = (this->Range[1] - this->Range[0]) / bins;
for (int i = 0; i < bins; ++i)
{
printf(" %f-%f: %d\n", this->Range[0] + i*width, this->Range[0] + (i+1)*width, gHist[i]);
}
}
}
};
}
namespace sensei
{
......
#include <vtkAOSDataArrayTemplate.h>
#include <vtkArrayDispatch.h>
#include <vtkCompositeDataIterator.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataArray.h>
#include <vtkDataObject.h>
#include <vtkFieldData.h>
#include <vtkObjectFactory.h>
#include <vtkSmartPointer.h>
#include "DataAdaptor.h"
#include <timer/Timer.h>
#include <algorithm>
#include <vector>
namespace
{
// Private worker for histogram method. Computes the local histogram on
// array (passed to operator()). To be used with vtkArrayDispatch.
//
// Inputs:
// range: Global range of data
// bins: Number of histogram bins
// array: Local data.
//
// Outputs:
// Histogram: The histogram of the local data.
struct HistogramWorker
{
const double *Range;
int Bins;
std::vector<unsigned int> Histogram;
HistogramWorker(const double *range, int bins) : Range(range), Bins(bins) {}
template <typename ArrayT>
void operator()(ArrayT *array)
{
assert(array);
assert(array->GetNumberOfComponents() == 1);
typedef typename ArrayT::ValueType ValueType;
ValueType width = static_cast<ValueType>((this->Range[1] - this->Range[0]) / this->Bins);
ValueType min = static_cast<ValueType>(this->Range[0]);
vtkIdType numTuples = array->GetNumberOfTuples();
// + 1 to store val == max. These will be moved to this last bin before
// returning (Avoids having to branch in the loop below);
this->Histogram.resize(this->Bins + 1, 0);
for (vtkIdType tIdx = 0; tIdx < numTuples; ++tIdx)
{
int bin = static_cast<int>((array->GetComponent(tIdx, 0) - min) / width);
++this->Histogram[bin];
}
// Merge the last two bins (the last is only when val == max)
this->Histogram[this->Bins-1] += this->Histogram[this->Bins];
this->Histogram.resize(this->Bins);
}
};
class vtkHistogram
{
double Range[2];
HistogramWorker *Worker;
public:
vtkHistogram()
{
this->Range[0] = VTK_DOUBLE_MAX;
this->Range[1] = VTK_DOUBLE_MIN;
this->Worker = NULL;
}
~vtkHistogram()
{
delete this->Worker;
}
void AddRange(vtkDataArray* da)
{
if (da)
{
double crange[2];
da->GetRange(crange);
this->Range[0] = std::min(this->Range[0], crange[0]);
this->Range[1] = std::max(this->Range[1], crange[1]);
}
}
void PreCompute(MPI_Comm comm, int bins)
{
double g_range[2];
// Find the global max/min
MPI_Allreduce(&this->Range[0], &g_range[0], 1, MPI_DOUBLE, MPI_MIN, comm);
MPI_Allreduce(&this->Range[1], &g_range[1], 1, MPI_DOUBLE, MPI_MAX, comm);
this->Range[0] = g_range[0];
this->Range[1] = g_range[1];
this->Worker = new HistogramWorker(this->Range, bins);
}
void Compute(vtkDataArray* da)
{
if (da)
{
vtkArrayDispatch::Dispatch::Execute(da, *this->Worker);
}
}
void PostCompute(MPI_Comm comm, int bins, const std::string& name)
{
std::vector<unsigned int> gHist(bins, 0);
MPI_Reduce(&this->Worker->Histogram[0], &gHist[0], bins, MPI_UNSIGNED, MPI_SUM, 0, comm);
int rank;
MPI_Comm_rank(comm, &rank);
if (rank == 0)
{
std::cout << "Histogram '" << name << "' (VTK):\n";
double width = (this->Range[1] - this->Range[0]) / bins;
for (int i = 0; i < bins; ++i)
{
printf(" %f-%f: %d\n", this->Range[0] + i*width, this->Range[0] + (i+1)*width, gHist[i]);
}
}
}
};
}
#include <vtkDataArrayDispatcher.h>
#include <vtkDataArrayTemplate.h>
#include <vtkSmartPointer.h>
#include "DataAdaptor.h"
#include <timer/Timer.h>
#include <algorithm>
#include <vector>
namespace
{
// Private worker for histogram method. Computes the local histogram on
// array (passed to operator()). To be used with vtkArrayDispatch.
//
// Inputs:
// range: Global range of data
// bins: Number of histogram bins
// array: Local data.
//
// Outputs:
// Histogram: The histogram of the local data.
struct HistogramWorker
{
const double *Range;
int Bins;
std::vector<unsigned int> Histogram;
HistogramWorker(const double *range, int bins) : Range(range), Bins(bins) {}
template <typename T>
void operator()(const vtkDataArrayDispatcherPointer<T>& array)
{
assert(array.NumberOfComponents == 1);
typedef T ValueType;
ValueType width = static_cast<ValueType>((this->Range[1] - this->Range[0]) / this->Bins);
ValueType min = static_cast<ValueType>(this->Range[0]);
vtkIdType numTuples = array.NumberOfTuples;
// + 1 to store val == max. These will be moved to this last bin before
// returning (Avoids having to branch in the loop below);
this->Histogram.resize(this->Bins + 1, 0);
for (vtkIdType tIdx = 0; tIdx < numTuples; ++tIdx)
{
int bin = static_cast<int>((array.RawPointer[tIdx] - min) / width);
++this->Histogram[bin];
}
// Merge the last two bins (the last is only when val == max)
this->Histogram[this->Bins-1] += this->Histogram[this->Bins];
this->Histogram.resize(this->Bins);
}
};
class vtkHistogram
{
double Range[2];
HistogramWorker *Worker;
public:
vtkHistogram()
{
this->Range[0] = VTK_DOUBLE_MAX;
this->Range[1] = VTK_DOUBLE_MIN;
this->Worker = NULL;
}
~vtkHistogram()
{
delete this->Worker;
}
void AddRange(vtkDataArray* da)
{
if (da)
{
double crange[2];
da->GetRange(crange);
this->Range[0] = std::min(this->Range[0], crange[0]);
this->Range[1] = std::max(this->Range[1], crange[1]);
}
}
void PreCompute(MPI_Comm comm, int bins)
{
double g_range[2];
// Find the global max/min
MPI_Allreduce(&this->Range[0], &g_range[0], 1, MPI_DOUBLE, MPI_MIN, comm);
MPI_Allreduce(&this->Range[1], &g_range[1], 1, MPI_DOUBLE, MPI_MAX, comm);
this->Range[0] = g_range[0];
this->Range[1] = g_range[1];
this->Worker = new HistogramWorker(this->Range, bins);
}
void Compute(vtkDataArray* da)
{
if (da)
{
vtkDataArrayDispatcher<HistogramWorker> dispatcher(*this->Worker);
dispatcher.Go(da);
}
}
void PostCompute(MPI_Comm comm, int bins, const std::string& name)
{
std::vector<unsigned int> gHist(bins, 0);
MPI_Reduce(&this->Worker->Histogram[0], &gHist[0], bins, MPI_UNSIGNED, MPI_SUM, 0, comm);
int rank;
MPI_Comm_rank(comm, &rank);
if (rank == 0)
{
std::cout << "Histogram '" << name << "' (VTK):\n";
double width = (this->Range[1] - this->Range[0]) / bins;
for (int i = 0; i < bins; ++i)
{
printf(" %f-%f: %d\n", this->Range[0] + i*width, this->Range[0] + (i+1)*width, gHist[i]);
}
}
}
};
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment