From 79183ebb79025d068076e849ad2cb0bd12db49d4 Mon Sep 17 00:00:00 2001
From: Martyn Gigg <martyn.gigg@stfc.ac.uk>
Date: Fri, 24 Jun 2011 12:08:35 +0000
Subject: [PATCH] The performance can be improved but Rebin2D should give the
correct values with the exception of masked bins at the moment. Re #3197
---
.../Framework/Algorithms/CMakeLists.txt | 9 +-
.../Algorithms/inc/MantidAlgorithms/Rebin2D.h | 42 +++++
.../Framework/Algorithms/src/Rebin2D.cpp | 155 ++++++++++++++--
.../Framework/Algorithms/test/Rebin2DTest.h | 174 ++++++++++++++++++
4 files changed, 358 insertions(+), 22 deletions(-)
create mode 100644 Code/Mantid/Framework/Algorithms/test/Rebin2DTest.h
diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
index 6f69c8afade..6ad803204f3 100644
--- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt
+++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
@@ -1,4 +1,5 @@
set ( SRC_FILES
+ # src/Q1DTOF.cpp
src/AbsorptionCorrection.cpp
src/AddSampleLog.cpp
src/AlignDetectors.cpp
@@ -111,12 +112,12 @@ set ( SRC_FILES
src/PowerLawCorrection.cpp
src/Q1D.cpp
src/Q1D2.cpp
-# src/Q1DTOF.cpp
src/Q1DWeighted.cpp
src/Qxy.cpp
src/ReadGroupsFromFile.cpp
src/RealFFT.cpp
src/Rebin.cpp
+ src/Rebin2D.cpp
src/RebinToWorkspace.cpp
src/Rebunch.cpp
src/Regroup.cpp
@@ -152,6 +153,7 @@ set ( SRC_FILES
)
set ( INC_FILES
+ # inc/MantidAlgorithms/Q1DTOF.h
inc/MantidAlgorithms/AbsorptionCorrection.h
inc/MantidAlgorithms/AddSampleLog.h
inc/MantidAlgorithms/AlignDetectors.h
@@ -265,12 +267,12 @@ set ( INC_FILES
inc/MantidAlgorithms/PowerLawCorrection.h
inc/MantidAlgorithms/Q1D.h
inc/MantidAlgorithms/Q1D2.h
-# inc/MantidAlgorithms/Q1DTOF.h
inc/MantidAlgorithms/Q1DWeighted.h
inc/MantidAlgorithms/Qxy.h
inc/MantidAlgorithms/ReadGroupsFromFile.h
inc/MantidAlgorithms/RealFFT.h
inc/MantidAlgorithms/Rebin.h
+ inc/MantidAlgorithms/Rebin2D.h
inc/MantidAlgorithms/RebinToWorkspace.h
inc/MantidAlgorithms/Rebunch.h
inc/MantidAlgorithms/Regroup.h
@@ -405,12 +407,13 @@ set ( TEST_FILES
test/PolynomialCorrectionTest.h
test/PowerLawCorrectionTest.h
test/PowerTest.h
- test/Q1DTest.h
test/Q1D2Test.h
+ test/Q1DTest.h
test/Q1DWeightedTest.h
test/QxyTest.h
test/ReadGroupsFromFileTest.h
test/RealFFTTest.h
+ test/Rebin2DTest.h
test/RebinTest.h
test/RebinToWorkspaceTest.h
test/RebunchTest.h
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin2D.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin2D.h
index 6c92efdf7c2..767353dea12 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin2D.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin2D.h
@@ -9,6 +9,14 @@
namespace Mantid
{
+ //------------------------------------------------------------------------------
+ // Forward declarations
+ //------------------------------------------------------------------------------
+ namespace Geometry
+ {
+ class ConvexPolygon;
+ }
+
namespace Algorithms
{
@@ -48,12 +56,46 @@ namespace Mantid
virtual const std::string category() const { return "Rebin";}
private:
+ /// A struct to store information about an intersection
+ struct BinWithWeight
+ {
+ /** Constructor
+ * @param i :: The index in the Y direction of the data bin
+ * @param j :: The index in the X direction of the data bin
+ * @param pointWeight :: The weight this point carries
+ */
+ BinWithWeight(const size_t i, const size_t j, const double pointWeight)
+ : yIndex(i), xIndex(j), weight(pointWeight) {}
+ /// The index in the Y direction of the data bin
+ size_t yIndex;
+ /// The index in the X direction of the data bin
+ size_t xIndex;
+ /// The weight this point carries
+ double weight;
+ };
+
/// Sets documentation strings for this algorithm
virtual void initDocs();
/// Initialise the properties
void init();
/// Run the algorithm
void exec();
+ /// Create a polygon
+ Geometry::ConvexPolygon createPolygon(const double x_i, const double x_ip1,
+ const double y_i, const double y_ip1) const;
+ /// Calculate the Y and E values for the given possible overlap
+ std::pair<double,double> calculateYE(API::MatrixWorkspace_const_sptr inputWS,
+ const MantidVec & oldYBins,
+ const Geometry::ConvexPolygon & outputPoly) const;
+ /// Find the overlap of the inputWS with the given polygon
+ std::vector<BinWithWeight> findIntersections(const MantidVec & oldAxis1,
+ const MantidVec & oldAxis2,
+ const Geometry::ConvexPolygon & poly) const;
+ /// Setup the output workspace
+ API::MatrixWorkspace_sptr createOutputWorkspace(API::MatrixWorkspace_const_sptr parent,
+ MantidVec &newXBins,
+ MantidVec &newYBins) const;
+
};
} // namespace Algorithms
diff --git a/Code/Mantid/Framework/Algorithms/src/Rebin2D.cpp b/Code/Mantid/Framework/Algorithms/src/Rebin2D.cpp
index 9b1da309b4d..884e702c8c6 100644
--- a/Code/Mantid/Framework/Algorithms/src/Rebin2D.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Rebin2D.cpp
@@ -7,6 +7,8 @@
#include "MantidKernel/VectorHelper.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/NumericAxis.h"
+#include "MantidGeometry/Math/ConvexPolygon.h"
+#include "MantidGeometry/Math/PolygonIntersection.h"
namespace Mantid
{
@@ -14,10 +16,13 @@ namespace Mantid
{
// Register the algorithm into the AlgorithmFactory
- //DECLARE_ALGORITHM(Rebin2D)
+ DECLARE_ALGORITHM(Rebin2D)
+
+ using namespace API;
+ using Geometry::ConvexPolygon;
+ using Geometry::Vertex2DList;
+ using Kernel::V2D;
- using API::MatrixWorkspace_const_sptr;
-
//--------------------------------------------------------------------------
// Public methods
//--------------------------------------------------------------------------
@@ -33,6 +38,9 @@ namespace Mantid
"new bin on the workspace");
}
+ //--------------------------------------------------------------------------
+ // Private methods
+ //--------------------------------------------------------------------------
/**
* Initialize the algorithm's properties.
*/
@@ -52,29 +60,138 @@ namespace Mantid
declareProperty(new ArrayProperty<double>("Axis2Binning", new RebinParamsValidator), docString);
}
- //----------------------------------------------------------------------------------------------
- /** Execute the algorithm.
+ /**
+ * Execute the algorithm.
*/
void Rebin2D::exec()
{
+ // Information to form input grid
MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
+ NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis*>(inputWS->getAxis(1));
+ if( !oldAxis2 )
+ {
+ throw std::invalid_argument("Vertical axis is not a numeric axis, cannot rebin. "
+ "If it is a spectra axis try running ConvertSpectrumAxis first.");
+ }
+ const std::vector<double> oldYBins = oldAxis2->createBinBoundaries();
- // We require two grids. The input workspace has axis 1 as non-histogram data so we
- // need to make bins out of axis1.
- const MantidVec oldXBins = inputWS->readX(0);
- const std::vector<double> oldYBins =
- dynamic_cast<API::NumericAxis*>(inputWS->getAxis(1))->createBinBoundaries();
- // Output bins
+ // Output grid and workspace. Fills in the new X and Y bin vectors
MantidVecPtr newXBins;
- const int nXValues =
- Kernel::VectorHelper::createAxisFromRebinParams(getProperty("Axis1Binning"), newXBins.access());
- MantidVecPtr newYBins;
- const int nYValues =
- Kernel::VectorHelper::createAxisFromRebinParams(getProperty("Axis2Binning"), newYBins.access());
- assert(nXValues > 0);
- assert(nYValues > 0);
-
+ MantidVec newYBins;
+ MatrixWorkspace_sptr outputWS = createOutputWorkspace(inputWS, newXBins.access(), newYBins);
+ const int64_t numHist = static_cast<int64_t>(outputWS->getNumberHistograms());
+ const int64_t numSig = static_cast<int64_t>(outputWS->blocksize());
+ Axis* newAxis2 = outputWS->getAxis(1);
+
+ Progress progress(this,0.0,1.0, numHist);
+ // For each box on the output workspace find the overlapping boxes on the input workspac
+ PARALLEL_FOR2(inputWS,outputWS)
+ for(int64_t i = 0; i < numHist; ++i)
+ {
+ PARALLEL_START_INTERUPT_REGION
+
+ outputWS->setX(i, newXBins);
+ // References to the Y and E data
+ MantidVec & outY = outputWS->dataY(i);
+ MantidVec & outE = outputWS->dataE(i);
+ // The Y axis value is the bin centre of the new Y bins
+ const double y_i = newYBins[i];
+ const double y_ip1 = newYBins[i+1];
+ newAxis2->setValue(i, 0.5*(y_i + y_ip1));
+
+ // For each bin on the new workspace, test compute overlap
+ for(int64_t j = 0; j < numSig; ++j)
+ {
+ ConvexPolygon newBin = ConvexPolygon((*newXBins)[j],(*newXBins)[j+1], y_i, y_ip1);
+ std::pair<double,double> contrib = calculateYE(inputWS, oldYBins, newBin);
+ outY[j] = contrib.first;
+ outE[j] = contrib.second;
+ }
+ progress.report();
+ PARALLEL_END_INTERUPT_REGION
+ }
+ PARALLEL_CHECK_INTERUPT_REGION
+
+ setProperty("OutputWorkspace", outputWS);
+ }
+
+ /**
+ * Calculate the Y and E values for the given possible overlap
+ * @param inputWS :: A pointer to the inputWS
+ * @param oldYBins :: A reference to the set of Y bin boundaries on the input WS
+ * @param newBin :: A reference to a polygon to test for overlap
+ */
+ std::pair<double,double>
+ Rebin2D::calculateYE(API::MatrixWorkspace_const_sptr inputWS,
+ const MantidVec & oldYBins, const ConvexPolygon & newBin) const
+ {
+ const MantidVec & oldXBins = inputWS->readX(0);
+ // Build a list intersection locations in terms of workspace indices
+ // along with corresponding weights from that location
+ std::vector<BinWithWeight> overlaps = findIntersections(oldXBins, oldYBins, newBin);
+ double totalY(0.0), totalE(0.0);
+ std::vector<BinWithWeight>::const_iterator iend = overlaps.end();
+ for(std::vector<BinWithWeight>::const_iterator itr = overlaps.begin();
+ itr != iend; ++itr)
+ {
+ totalY += inputWS->readY(itr->yIndex)[itr->xIndex] * itr->weight;
+ totalE += std::pow(inputWS->readE(itr->yIndex)[itr->xIndex] * itr->weight, 2);
+ }
+ return std::make_pair(totalY,std::sqrt(totalE));
+ }
+
+
+ /**
+ * Find the overlap of the inputWS with the given polygon
+ * @param oldAxis1 :: Axis 1 bin boundaries from the input grid
+ * @param oldAxis2 :: Axis 2 bin boundaries from the input grid
+ * @param newPoly :: The new polygon to test
+ * @returns A list of intersection locations with weights of the overlap
+ */
+ std::vector<Rebin2D::BinWithWeight>
+ Rebin2D::findIntersections(const MantidVec & oldAxis1, const MantidVec & oldAxis2,
+ const Geometry::ConvexPolygon & newPoly) const
+ {
+ std::vector<BinWithWeight> overlaps;
+ overlaps.reserve(5); // Have a guess at a posible limit
+
+ const size_t nxpoints(oldAxis1.size()-1), nypoints(oldAxis2.size()-1);
+ Vertex2DList vertices(4);
+ for(size_t i = 0; i < nypoints; ++i)
+ {
+ for(size_t j = 0; j < nxpoints; ++j)
+ {
+ ConvexPolygon oldPoly(oldAxis1[j], oldAxis1[j+1], oldAxis2[i], oldAxis2[i+1]);
+ try
+ {
+ ConvexPolygon overlap = chasingEdgeIntersect(newPoly, oldPoly);
+ overlaps.push_back(BinWithWeight(i,j,overlap.area()/oldPoly.area()));
+ }
+ catch(std::runtime_error &)
+ {}
+ }
+ }
+ return overlaps;
+ }
+
+ /**
+ * Setup the output workspace
+ * @param parent :: A pointer to the input workspace
+ * @param newXBins[out] :: An output vector to be filled with the new X bin boundaries
+ * @param newYBins[out] :: An output vector to be filled with the new Y bin boundaries
+ * @return A pointer to the output workspace
+ */
+ MatrixWorkspace_sptr Rebin2D::createOutputWorkspace(MatrixWorkspace_const_sptr parent,
+ MantidVec & newXBins,
+ MantidVec & newYBins) const
+ {
+ using Kernel::VectorHelper::createAxisFromRebinParams;
+ // First create the two sets of bin boundaries
+ const int newXSize = createAxisFromRebinParams(getProperty("Axis1Binning"), newXBins);
+ const int newYSize = createAxisFromRebinParams(getProperty("Axis2Binning"), newYBins);
+ // and now the workspace
+ return WorkspaceFactory::Instance().create(parent,newYSize-1,newXSize,newXSize-1);
}
} // namespace Algorithms
diff --git a/Code/Mantid/Framework/Algorithms/test/Rebin2DTest.h b/Code/Mantid/Framework/Algorithms/test/Rebin2DTest.h
new file mode 100644
index 00000000000..8619276dc69
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/test/Rebin2DTest.h
@@ -0,0 +1,174 @@
+#ifndef MANTID_ALGORITHMS_REBIN2DTEST_H_
+#define MANTID_ALGORITHMS_REBIN2DTEST_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidAlgorithms/Rebin2D.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidAPI/NumericAxis.h"
+
+#include "MantidKernel/Timer.h"
+
+using Mantid::Algorithms::Rebin2D;
+using namespace Mantid::API;
+
+namespace
+{
+ /**
+ * Shared methods between the unit test and performance test
+ */
+
+ /// Return the input workspace. All Y values are 2 and E values sqrt(2)
+ MatrixWorkspace_sptr makeInputWS(const bool large = false)
+ {
+ size_t nhist(0), nbins(0);
+ double x0(0.0), deltax(0.0);
+
+ if( large )
+ {
+ nhist = 200;
+ nbins = 200;
+ x0 = 100.;
+ deltax = 200.;
+ }
+ else
+ {
+ nhist = 10;
+ nbins = 10;
+ x0 = 5.0;
+ deltax = 1.0;
+ }
+
+ MatrixWorkspace_sptr ws = WorkspaceCreationHelper::Create2DWorkspaceBinned(int(nhist), int(nbins), x0, deltax);
+
+ // We need something other than a spectrum axis, call this one theta
+ NumericAxis* const thetaAxis = new NumericAxis(nhist);
+ for(size_t i = 0; i < nhist; ++i)
+ {
+ thetaAxis->setValue(i, static_cast<double>(i));
+ }
+ ws->replaceAxis(1, thetaAxis);
+ return ws;
+ }
+
+ MatrixWorkspace_sptr runAlgorithm(MatrixWorkspace_sptr inputWS,
+ const std::string & axis1Params,
+ const std::string & axis2Params)
+ {
+ // Name of the output workspace.
+ std::string outWSName("Rebin2DTest_OutputWS");
+
+ Rebin2D alg;
+ TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+ TS_ASSERT( alg.isInitialized() )
+ TS_ASSERT_THROWS_NOTHING( alg.setProperty("InputWorkspace", inputWS) );
+ TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) );
+ TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Axis1Binning", axis1Params) );
+ TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Axis2Binning", axis2Params) );
+ TS_ASSERT_THROWS_NOTHING( alg.execute(); );
+ TS_ASSERT( alg.isExecuted() );
+
+ MatrixWorkspace_sptr outputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(
+ AnalysisDataService::Instance().retrieve(outWSName));
+ TS_ASSERT(outputWS);
+ return outputWS;
+ }
+
+}
+
+class Rebin2DTest : public CxxTest::TestSuite
+{
+public:
+
+ void test_Init()
+ {
+ Rebin2D alg;
+ TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+ TS_ASSERT( alg.isInitialized() )
+ }
+
+ void test_Rebin2D_With_Axis2_Unchanged()
+ {
+ MatrixWorkspace_sptr inputWS = makeInputWS(); //10 histograms, 10 bins
+ MatrixWorkspace_sptr outputWS = runAlgorithm(inputWS, "5.,2.,15.", "-0.5,1,9.5");
+
+ // Check values
+ const size_t nxvalues(6), nhist(10);
+ TS_ASSERT_EQUALS(outputWS->getNumberHistograms(), nhist);
+ // Axis sizes
+ TS_ASSERT_EQUALS(outputWS->getAxis(0)->length(), nxvalues);
+ TS_ASSERT_EQUALS(outputWS->getAxis(1)->length(), nhist);
+ TS_ASSERT_EQUALS(outputWS->readX(0).size(), nxvalues);
+ TS_ASSERT_EQUALS(outputWS->readY(0).size(), nxvalues - 1);
+
+ Mantid::API::Axis *newYAxis = outputWS->getAxis(1);
+ for(size_t i = 0; i < nhist; ++i)
+ {
+ for(size_t j = 0; j < nxvalues - 1; ++j)
+ {
+ const double expected(5.0 + 2.0*static_cast<double>(j));
+ TS_ASSERT_DELTA(outputWS->readX(i)[j], expected, DBL_EPSILON);
+ TS_ASSERT_DELTA(outputWS->readY(i)[j], 4.0, DBL_EPSILON);
+ TS_ASSERT_DELTA(outputWS->readE(i)[j], 2.0, DBL_EPSILON);
+ }
+ // Final X boundary
+ TS_ASSERT_DELTA(outputWS->readX(i)[nxvalues-1], 15.0, DBL_EPSILON);
+ // The new Y axis value should be the centre point bin values
+ TS_ASSERT_DELTA((*newYAxis)(i), static_cast<double>(i), DBL_EPSILON);
+ }
+
+ // Clean up
+ AnalysisDataService::Instance().remove(outputWS->getName());
+ }
+
+ void test_Rebin2D_With_Axis1_Unchanged()
+ {
+ MatrixWorkspace_sptr inputWS = makeInputWS(); //10 histograms, 10 bins
+ MatrixWorkspace_sptr outputWS = runAlgorithm(inputWS, "5.,1.,15.", "-0.5,2,9.5");
+ // Check values
+ const size_t nxvalues(11), nhist(5);
+ TS_ASSERT_EQUALS(outputWS->getNumberHistograms(), nhist);
+ // Axis sizes
+ TS_ASSERT_EQUALS(outputWS->getAxis(0)->length(), nxvalues);
+ TS_ASSERT_EQUALS(outputWS->getAxis(1)->length(), nhist);
+ TS_ASSERT_EQUALS(outputWS->readX(0).size(), nxvalues);
+ TS_ASSERT_EQUALS(outputWS->readY(0).size(), nxvalues - 1);
+
+ Mantid::API::Axis *newYAxis = outputWS->getAxis(1);
+ for(size_t i = 0; i < nhist; ++i)
+ {
+ for(size_t j = 0; j < nxvalues - 1; ++j)
+ {
+ TS_ASSERT_DELTA(outputWS->readX(i)[j], 5.0 + static_cast<double>(j), DBL_EPSILON);
+ TS_ASSERT_DELTA(outputWS->readY(i)[j], 4.0, DBL_EPSILON);
+ TS_ASSERT_DELTA(outputWS->readE(i)[j], 2.0, DBL_EPSILON);
+ }
+ // Final X boundary
+ TS_ASSERT_DELTA(outputWS->readX(i)[nxvalues-1], 15.0, DBL_EPSILON);
+ // The new Y axis value should be the centre point bin values
+ TS_ASSERT_DELTA((*newYAxis)(i), 0.5 + 2.0*static_cast<double>(i), DBL_EPSILON);
+ }
+ }
+};
+
+
+//------------------------------------------------------------------------------
+// Performance test
+//------------------------------------------------------------------------------
+
+class Rebin2DTestPerformance : public CxxTest::TestSuite
+{
+
+public:
+
+ void test_On_Large_Workspace()
+ {
+ MatrixWorkspace_sptr inputWS = makeInputWS();
+ runAlgorithm(inputWS, "200,250,40000", "-0.5,5,199.5");
+ }
+
+
+};
+
+
+#endif /* MANTID_ALGORITHMS_REBIN2DTEST_H_ */
+
--
GitLab