Commit 5c73705f authored by Utkarsh Ayachit's avatar Utkarsh Ayachit
Browse files

Make histogram robust for composite blocks distributed across ranks.

parent 48bfa437
......@@ -6,8 +6,6 @@ set (sources
if(ENABLE_SENSEI)
list(APPEND sources
vtk_histogram.cxx
vtk_histogram.h
HistogramAnalysisAdaptor.cxx
HistogramAnalysisAdaptor.h
)
......
#include "HistogramAnalysisAdaptor.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataSet.h"
#include "vtkDataObject.h"
#include "vtkFieldData.h"
#include "vtkInsituDataAdaptor.h"
#include "vtkObjectFactory.h"
#include "vtkSmartPointer.h"
#include <vtkAOSDataArrayTemplate.h>
#include <vtkArrayDispatch.h>
#include <vtkCompositeDataIterator.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataArray.h>
#include <vtkDataObject.h>
#include <vtkFieldData.h>
#include <vtkInsituDataAdaptor.h>
#include <vtkObjectFactory.h>
#include <vtkSmartPointer.h>
#include "vtk_histogram.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(g_range, bins);
}
void Compute(vtkDataArray* da)
{
if (da)
{
vtkArrayDispatch::Dispatch::Execute(da, *this->Worker);
}
}
void PostCompute(MPI_Comm comm, int bins, const char* 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]);
}
}
}
};
}
vtkStandardNewMacro(HistogramAnalysisAdaptor);
//-----------------------------------------------------------------------------
......@@ -37,38 +152,50 @@ void HistogramAnalysisAdaptor::Initialize(
//-----------------------------------------------------------------------------
bool HistogramAnalysisAdaptor::Execute(vtkInsituDataAdaptor* data)
{
bool retval = false;
vtkHistogram histogram;
vtkDataObject* mesh = data->GetMesh(/*structure_only*/true);
if (mesh && data->AddArray(mesh, this->Association, this->ArrayName.c_str()))
if (mesh == NULL || !data->AddArray(mesh, this->Association, this->ArrayName.c_str()))
{
if (vtkCompositeDataSet* cd = vtkCompositeDataSet::SafeDownCast(mesh))
histogram.PreCompute(this->Communicator, this->Bins);
histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
return true;
}
if (vtkCompositeDataSet* cd = vtkCompositeDataSet::SafeDownCast(mesh))
{
vtkSmartPointer<vtkCompositeDataIterator> iter;
iter.TakeReference(cd->NewIterator());
for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
{
vtkSmartPointer<vtkCompositeDataIterator> iter;
iter.TakeReference(cd->NewIterator());
for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
{
retval = this->Execute(iter->GetCurrentDataObject()) || retval;
}
vtkDataArray* array = this->GetArray(iter->GetCurrentDataObject());
histogram.AddRange(array);
}
else
histogram.PreCompute(this->Communicator, this->Bins);
for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
{
retval = this->Execute(mesh);
vtkDataArray* array = this->GetArray(iter->GetCurrentDataObject());
histogram.Compute(array);
}
histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
}
return retval;
else
{
vtkDataArray* array = this->GetArray(mesh);
histogram.AddRange(array);
histogram.PreCompute(this->Communicator, this->Bins);
histogram.Compute(array);
histogram.PostCompute(this->Communicator, this->Bins, this->ArrayName.c_str());
}
return true;
}
//-----------------------------------------------------------------------------
bool HistogramAnalysisAdaptor::Execute(vtkDataObject* mesh)
vtkDataArray* HistogramAnalysisAdaptor::GetArray(vtkDataObject* dobj)
{
if (vtkFieldData* fd = mesh->GetAttributesAsFieldData(this->Association))
assert(dobj != NULL && vtkCompositeDataSet::SafeDownCast(dobj) == NULL);
if (vtkFieldData* fd = dobj->GetAttributesAsFieldData(this->Association))
{
if (vtkDataArray* da = fd->GetArray(this->ArrayName.c_str()))
{
vtk_histogram(this->Communicator, da, this->Bins);
}
return true;
return fd->GetArray(this->ArrayName.c_str());
}
return false;
return NULL;
}
......@@ -6,6 +6,7 @@
#include <mpi.h>
class vtkDataObject;
class vtkDataArray;
/// HistogramAnalysisAdaptor is a vtkInsituAnalysisAdaptor specialization for
/// histogram analysis.
class HistogramAnalysisAdaptor : public vtkInsituAnalysisAdaptor
......@@ -23,7 +24,8 @@ protected:
HistogramAnalysisAdaptor();
virtual ~HistogramAnalysisAdaptor();
bool Execute(vtkDataObject* mesh);
vtkDataArray* GetArray(vtkDataObject* dobj);
MPI_Comm Communicator;
int Bins;
......
......@@ -27,9 +27,11 @@ target_link_libraries(oscillator PRIVATE util ${CMAKE_THREAD_LIBS_INIT})
if(ENABLE_SENSEI)
target_link_libraries(oscillator PRIVATE core)
endif()
if(ENABLE_HISTOGRAM AND FALSE)
if(ENABLE_HISTOGRAM)
target_link_libraries(oscillator PRIVATE histogram)
endif()
if(ENABLE_AUTOCORRELATION)
target_link_libraries(oscillator PRIVATE autocorrelation)
endif()
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