diff --git a/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h b/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h index 16623ba3089dbc0c0b39bafbb3cc5b57ae1c561c..63100182aa7c00d7c65084f54dc1ef25b5971b87 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h @@ -48,7 +48,7 @@ public: return {"Stitch1DMany"}; } /// Does the x-axis have non-zero errors - bool hasNonzeroErrors(Mantid::API::MatrixWorkspace_sptr ws); + bool hasNonzeroErrors(Mantid::API::MatrixWorkspace_sptr &ws); /// Cross-check properties with each other @see IAlgorithm::validateInputs std::map<std::string, std::string> validateInputs() override; @@ -68,9 +68,10 @@ private: const double &intesectionMax) const; /// Get the rebin parameters - std::vector<double> getRebinParams(Mantid::API::MatrixWorkspace_sptr &lhsWS, - Mantid::API::MatrixWorkspace_sptr &rhsWS, - const bool scaleRHS) const; + std::vector<double> + getRebinParams(Mantid::API::MatrixWorkspace_const_sptr &lhsWS, + Mantid::API::MatrixWorkspace_const_sptr &rhsWS, + const bool scaleRHS) const; /// Perform rebin Mantid::API::MatrixWorkspace_sptr rebin(Mantid::API::MatrixWorkspace_sptr &input, @@ -79,7 +80,7 @@ private: Mantid::API::MatrixWorkspace_sptr integration(Mantid::API::MatrixWorkspace_sptr &input, const double &start, const double &stop); - Mantid::API::MatrixWorkspace_sptr singleValueWS(double val); + Mantid::API::MatrixWorkspace_sptr singleValueWS(const double &val); /// Calculate the weighted mean Mantid::API::MatrixWorkspace_sptr weightedMean(Mantid::API::MatrixWorkspace_sptr &inOne, @@ -89,8 +90,7 @@ private: conjoinXAxis(Mantid::API::MatrixWorkspace_sptr &inOne, Mantid::API::MatrixWorkspace_sptr &inTwo); /// Sort x axis - Mantid::API::MatrixWorkspace_sptr - sortXAxis(Mantid::API::MatrixWorkspace_sptr &inWS); + void sortXAxis(Mantid::API::MatrixWorkspace_sptr &ws); /// Find the start and end indexes boost::tuple<int, int> findStartEndIndexes(double startOverlap, double endOverlap, @@ -99,7 +99,7 @@ private: Mantid::API::MatrixWorkspace_sptr maskAllBut(int a1, int a2, Mantid::API::MatrixWorkspace_sptr &source); /// Mask out everything but the data in the ranges, but do it inplace. - void maskInPlace(int a1, int a2, Mantid::API::MatrixWorkspace_sptr source); + void maskInPlace(int a1, int a2, Mantid::API::MatrixWorkspace_sptr &source); /// Add back in any special values void reinsertSpecialValues(Mantid::API::MatrixWorkspace_sptr ws); /// Range tolerance @@ -108,10 +108,10 @@ private: double m_scaleFactor; double m_errorScaleFactor; /// Scale workspace (left hand side or right hand side) - void scaleWorkspace(Mantid::API::MatrixWorkspace_sptr ws, - API::MatrixWorkspace_sptr divident, - API::MatrixWorkspace_sptr divisor, - Mantid::API::MatrixWorkspace_sptr dxWS); + void scaleWorkspace(API::MatrixWorkspace_sptr &ws, + API::MatrixWorkspace_sptr ÷nt, + API::MatrixWorkspace_sptr &divisor, + API::MatrixWorkspace_const_sptr &dxWS); /// Index per workspace spectra of Nans SpecialTypeIndexes m_nanEIndexes; SpecialTypeIndexes m_nanYIndexes; diff --git a/Framework/Algorithms/src/ConjoinXRuns.cpp b/Framework/Algorithms/src/ConjoinXRuns.cpp index b9fc326800d0c28120bc9f78085fd9a4785ba6d1..fdfa226750fa19e8478216330bcdaf3763f450fa 100644 --- a/Framework/Algorithms/src/ConjoinXRuns.cpp +++ b/Framework/Algorithms/src/ConjoinXRuns.cpp @@ -3,6 +3,7 @@ #include "MantidAPI/ADSValidator.h" #include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/Axis.h" +#include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/Run.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/WorkspaceGroup.h" diff --git a/Framework/Algorithms/src/Stitch1D.cpp b/Framework/Algorithms/src/Stitch1D.cpp index a18f5bbbdf110eac87151980c30add666272c131..10416dce15f41b2cc06c60287b2454185ac9cf57 100644 --- a/Framework/Algorithms/src/Stitch1D.cpp +++ b/Framework/Algorithms/src/Stitch1D.cpp @@ -1,6 +1,7 @@ #include "MantidAlgorithms/Stitch1D.h" #include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/Workspace.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/WorkspaceProperty.h" #include "MantidHistogramData/HistogramE.h" @@ -30,8 +31,8 @@ using Mantid::HistogramData::Points; namespace { using MinMaxTuple = boost::tuple<double, double>; -MinMaxTuple calculateXIntersection(MatrixWorkspace_sptr lhsWS, - MatrixWorkspace_sptr rhsWS) { +MinMaxTuple calculateXIntersection(MatrixWorkspace_const_sptr &lhsWS, + MatrixWorkspace_const_sptr &rhsWS) { return MinMaxTuple(rhsWS->x(0).front(), lhsWS->x(0).back()); } @@ -93,8 +94,9 @@ MatrixWorkspace_sptr Stitch1D::maskAllBut(int a1, int a2, * @param a1 : start position in X * @param a2 : end position in X * @param source : Workspace to mask. + * @return Masked workspace. */ -void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr source) { +void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr &source) { const int histogramCount = static_cast<int>(source->getNumberHistograms()); PARALLEL_FOR_IF(Kernel::threadSafe(*source)) for (int i = 0; i < histogramCount; ++i) { @@ -164,8 +166,8 @@ void Stitch1D::init() { */ std::map<std::string, std::string> Stitch1D::validateInputs(void) { std::map<std::string, std::string> issues; - MatrixWorkspace_sptr lhs = getProperty("LHSWorkspace"); - MatrixWorkspace_sptr rhs = getProperty("RHSWorkspace"); + MatrixWorkspace_const_sptr lhs = getProperty("LHSWorkspace"); + MatrixWorkspace_const_sptr rhs = getProperty("RHSWorkspace"); if (lhs->isHistogramData() && !rhs->isHistogramData()) issues["RHSWorkspace"] = "Must be a histogram like LHSWorkspace."; if (!lhs->isHistogramData() && rhs->isHistogramData()) @@ -241,8 +243,8 @@ double Stitch1D::getEndOverlap(const double &intesectionMin, @param scaleRHS :: Scale the right hand side workspace @return a vector<double> contianing the rebinning parameters */ -std::vector<double> Stitch1D::getRebinParams(MatrixWorkspace_sptr &lhsWS, - MatrixWorkspace_sptr &rhsWS, +std::vector<double> Stitch1D::getRebinParams(MatrixWorkspace_const_sptr &lhsWS, + MatrixWorkspace_const_sptr &rhsWS, const bool scaleRHS) const { std::vector<double> inputParams = this->getProperty("Params"); Property *prop = this->getProperty("Params"); @@ -392,8 +394,8 @@ MatrixWorkspace_sptr Stitch1D::weightedMean(MatrixWorkspace_sptr &inOne, */ MatrixWorkspace_sptr Stitch1D::conjoinXAxis(MatrixWorkspace_sptr &inOne, MatrixWorkspace_sptr &inTwo) { - std::string in1 = "__Stitch1D_intermediate_workspace_1__"; - std::string in2 = "__Stitch1D_intermediate_workspace_2__"; + const std::string in1 = "__Stitch1D_intermediate_workspace_1__"; + const std::string in2 = "__Stitch1D_intermediate_workspace_2__"; Mantid::API::AnalysisDataService::Instance().addOrReplace(in1, inOne); Mantid::API::AnalysisDataService::Instance().addOrReplace(in2, inTwo); auto conjoinX = this->createChildAlgorithm("ConjoinXRuns"); @@ -407,51 +409,57 @@ MatrixWorkspace_sptr Stitch1D::conjoinXAxis(MatrixWorkspace_sptr &inOne, } /** Sort x axis - @param inWS :: Input workspace + @param ws :: Input workspace @return A shared pointer to the resulting MatrixWorkspace */ -MatrixWorkspace_sptr Stitch1D::sortXAxis(MatrixWorkspace_sptr &inWS) { +void Stitch1D::sortXAxis(MatrixWorkspace_sptr &ws) { // using a std::multimap (keys are sorted) std::multimap<double, double> sorter; - for (size_t i = 0; i < inWS->getNumberHistograms(); ++i) { - for (size_t k = 0; k < inWS->size(); ++k) { - sorter.insert(std::pair<double, double>(inWS->x(i)[k], inWS->y(i)[k])); - sorter.insert(std::pair<double, double>(inWS->x(i)[k], inWS->e(i)[k])); - if (inWS->hasDx(i)) - sorter.insert(std::pair<double, double>(inWS->x(i)[k], inWS->dx(i)[k])); + double x_value; + PARALLEL_FOR_IF(Kernel::threadSafe(*ws)) + for (size_t i = 0; i < ws->getNumberHistograms(); ++i) { + PARALLEL_START_INTERUPT_REGION + for (size_t k = 0; k < ws->size(); ++k) { + x_value = ws->x(i)[k]; + sorter.insert(std::pair<double, double>(x_value, ws->y(i)[k])); + sorter.insert(std::pair<double, double>(x_value, ws->e(i)[k])); + if (ws->hasDx(i)) + sorter.insert(std::pair<double, double>(x_value, ws->dx(i)[k])); } - Mantid::MantidVec vecx(inWS->size()); - Mantid::MantidVec vecy(inWS->size()); - Mantid::MantidVec vece(inWS->size()); - Mantid::MantidVec vecdx(inWS->size()); + size_t ws_size = ws->size(); + Mantid::MantidVec vecx(ws_size); + Mantid::MantidVec vecy(ws_size); + Mantid::MantidVec vece(ws_size); + Mantid::MantidVec vecdx(ws_size); int l = 0; for (auto it = sorter.cbegin(); it != sorter.cend(); ++it) { vecx[l] = it->first; vecy[l] = it->second; vece[l] = (++it)->second; - if (inWS->hasDx(i)) + if (ws->hasDx(i)) vecdx[l] = (++it)->second; ++l; } auto x = make_cow<HistogramX>(std::move(vecx)); auto y = make_cow<HistogramY>(std::move(vecy)); auto e = make_cow<HistogramE>(std::move(vece)); - inWS->setSharedX(i, x); - inWS->setSharedY(i, y); - inWS->setSharedE(i, e); - if (inWS->hasDx(i)) { + ws->setSharedX(i, x); + ws->setSharedY(i, y); + ws->setSharedE(i, e); + if (ws->hasDx(i)) { auto dx = make_cow<HistogramDx>(std::move(vecdx)); - inWS->setSharedDx(i, dx); + ws->setSharedDx(i, dx); } + PARALLEL_END_INTERUPT_REGION } - return inWS; + PARALLEL_CHECK_INTERUPT_REGION } /** Runs the CreateSingleValuedWorkspace Algorithm as a child @param val :: The double to convert to a single value workspace @return A shared pointer to the resulting MatrixWorkspace */ -MatrixWorkspace_sptr Stitch1D::singleValueWS(double val) { +MatrixWorkspace_sptr Stitch1D::singleValueWS(const double &val) { auto singleValueWS = this->createChildAlgorithm("CreateSingleValuedWorkspace"); singleValueWS->initialize(); @@ -484,7 +492,7 @@ Stitch1D::findStartEndIndexes(double startOverlap, double endOverlap, @param ws :: The input workspace @return True if there are any non-zero errors in the workspace */ -bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws) { +bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr &ws) { int64_t ws_size = static_cast<int64_t>(ws->getNumberHistograms()); bool hasNonZeroErrors = false; PARALLEL_FOR_IF(Kernel::threadSafe(*ws)) @@ -515,10 +523,10 @@ bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws) { * @param divisor * @param dxWS :: A MatrixWorkspace (size of ws) containing Dx values */ -void Stitch1D::scaleWorkspace(MatrixWorkspace_sptr ws, - MatrixWorkspace_sptr divident, - MatrixWorkspace_sptr divisor, - MatrixWorkspace_sptr dxWS) { +void Stitch1D::scaleWorkspace(MatrixWorkspace_sptr &ws, + MatrixWorkspace_sptr ÷nt, + MatrixWorkspace_sptr &divisor, + MatrixWorkspace_const_sptr &dxWS) { const auto ratio = divident / divisor; ws *= ratio; // We lost Dx values (Multiply) and need to get them back for point data @@ -545,8 +553,8 @@ void Stitch1D::scaleWorkspace(MatrixWorkspace_sptr ws, /** Execute the algorithm. */ void Stitch1D::exec() { - MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace"); - MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace"); + MatrixWorkspace_const_sptr lhsWS = this->getProperty("LHSWorkspace"); + MatrixWorkspace_const_sptr rhsWS = this->getProperty("RHSWorkspace"); const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS); const size_t histogramCount = rhsWS->getNumberHistograms(); @@ -577,8 +585,8 @@ void Stitch1D::exec() { const bool scaleRHS = this->getProperty("ScaleRHSWorkspace"); - MatrixWorkspace_sptr lhs; - MatrixWorkspace_sptr rhs; + MatrixWorkspace_sptr lhs = boost::const_pointer_cast<MatrixWorkspace>(lhsWS); + MatrixWorkspace_sptr rhs = boost::const_pointer_cast<MatrixWorkspace>(rhsWS); if (lhsWS->isHistogramData()) { MantidVec params = getRebinParams(lhsWS, rhsWS, scaleRHS); const double &xMin = params.front(); @@ -605,11 +613,8 @@ void Stitch1D::exec() { endOverlap % xMax); throw std::runtime_error(message); } - lhs = rebin(lhsWS, params); - rhs = rebin(rhsWS, params); - } else { - lhs = lhsWS; - rhs = rhsWS; + lhs = rebin(lhs, params); + rhs = rebin(rhs, params); } m_scaleFactor = this->getProperty("ManualScaleFactor"); @@ -622,10 +627,8 @@ void Stitch1D::exec() { else lhs *= manualScaleFactorWS; } else { - const auto rhsOverlapIntegrated = - integration(rhs, startOverlap, endOverlap); - const auto lhsOverlapIntegrated = - integration(lhs, startOverlap, endOverlap); + auto rhsOverlapIntegrated = integration(rhs, startOverlap, endOverlap); + auto lhsOverlapIntegrated = integration(lhs, startOverlap, endOverlap); if (scaleRHS) scaleWorkspace(rhs, lhsOverlapIntegrated, rhsOverlapIntegrated, rhsWS); else @@ -664,10 +667,10 @@ void Stitch1D::exec() { result = lhs + overlapave + rhs; reinsertSpecialValues(result); } else { // The input workspaces are point data ... join & sort - MatrixWorkspace_sptr ws = conjoinXAxis(lhs, rhs); - if (!ws) + result = conjoinXAxis(lhs, rhs); + if (!result) g_log.error("Could not retrieve joined workspace."); - result = sortXAxis(ws); + sortXAxis(result); } setProperty("OutputWorkspace", result); setProperty("OutScaleFactor", m_scaleFactor); diff --git a/Framework/Algorithms/test/Stitch1DTest.h b/Framework/Algorithms/test/Stitch1DTest.h index b098ebf63b3d6af762b855be2f0b07ea94b4687a..cd1b2a6dce1f5ab280d58da1f48ade7f35ce49cc 100644 --- a/Framework/Algorithms/test/Stitch1DTest.h +++ b/Framework/Algorithms/test/Stitch1DTest.h @@ -3,6 +3,7 @@ #include <cxxtest/TestSuite.h> +#include "MantidAlgorithms/CompareWorkspaces.h" #include "MantidAlgorithms/Stitch1D.h" #include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/Axis.h" @@ -279,6 +280,38 @@ public: TS_ASSERT_EQUALS(scaleFactor, 1.); // Default scale factor } + void test_point_data_input_workspace_not_modified() { + const auto &x1 = HistogramX(3, LinearGenerator(1., 1.)); + const auto &y1 = HistogramY(3, LinearGenerator(1., 1.)); + const auto &e = HistogramE(3, LinearGenerator(7., -1.)); + const auto &dx1 = HistogramDx(3, LinearGenerator(-3., 1.)); + auto ws1 = createWorkspace(x1, y1, e, dx1); + + const auto &x2 = HistogramX(3, LinearGenerator(2.1, 1.)); + const auto &y2 = HistogramY(3, LinearGenerator(5., 1.)); + const auto &dx2 = HistogramDx(3, LinearGenerator(-9., 0.)); + auto ws2 = createWorkspace(x2, y2, e, dx2); + + Stitch1D alg; + alg.setChild(true); + alg.setRethrows(true); + alg.initialize(); + alg.setProperty("LHSWorkspace", ws1); + alg.setProperty("RHSWorkspace", ws2); + alg.setProperty("UseManualScaleFactor", true); + alg.setPropertyValue("OutputWorkspace", "dummy_value"); + alg.execute(); + TS_ASSERT(alg.isExecuted()); + Mantid::Algorithms::CompareWorkspaces compare; + compare.initialize(); + compare.setRethrows(true); + TS_ASSERT_THROWS_NOTHING(compare.setProperty("Workspace1", ws1)); + TS_ASSERT_THROWS_NOTHING(compare.setProperty("Workspace2", ws2)); + TS_ASSERT(compare.execute()); + TS_ASSERT(compare.isExecuted()); + TS_ASSERT_EQUALS(compare.getPropertyValue("Result"), "1"); + } + void test_point_data_with_dx() { const auto &x1 = HistogramX(3, LinearGenerator(1., 1.)); const auto &y1 = HistogramY(3, LinearGenerator(1., 1.)); @@ -301,7 +334,7 @@ public: alg.execute(); TS_ASSERT(alg.isExecuted()); double scaleFactor = alg.getProperty("OutScaleFactor"); - TS_ASSERT_EQUALS(scaleFactor, 1./3.); + TS_ASSERT_EQUALS(scaleFactor, 1. / 3.); } void test_point_data_without_dx() { @@ -325,18 +358,26 @@ public: alg.execute(); TS_ASSERT(alg.isExecuted()); double scaleFactor = alg.getProperty("OutScaleFactor"); - TS_ASSERT_EQUALS(scaleFactor, 1./3.); + TS_ASSERT_EQUALS(scaleFactor, 1. / 3.); } void test_histogram_workspaces_pass() { - auto histo_ws = make_arbitrary_histogram_ws(); + auto ws1 = make_arbitrary_histogram_ws(); const auto &x = HistogramX(3, LinearGenerator(-1.2, 0.2)); const auto &y = HistogramY(2, LinearGenerator(1., 1.0)); const auto &e = HistogramE(2, 1.); const auto &dx = HistogramDx(2, LinearGenerator(-3., 0.1)); - auto histo_ws_2 = createWorkspace(x, y, e, dx); + auto ws2 = createWorkspace(x, y, e, dx); TSM_ASSERT_THROWS_NOTHING("Histogram workspaces should pass", - do_stitch1D(histo_ws_2, histo_ws)); + do_stitch1D(ws2, ws1)); + Mantid::Algorithms::CompareWorkspaces compare; + compare.initialize(); + compare.setRethrows(true); + TS_ASSERT_THROWS_NOTHING(compare.setProperty("Workspace1", ws1)); + TS_ASSERT_THROWS_NOTHING(compare.setProperty("Workspace2", ws2)); + TS_ASSERT(compare.execute()); + TS_ASSERT(compare.isExecuted()); + TS_ASSERT_EQUALS(compare.getPropertyValue("Result"), "1"); } void test_input_validation() { @@ -497,7 +538,8 @@ public: TS_ASSERT_EQUALS(temp, 0.); } // Check that the output X-Values are correct. - // truncate the input and oputput x values to 6 decimal places to eliminate + // truncate the input and oputput x values to 6 decimal places to + // eliminate // insignificant error auto xCopy = this->x; std::transform(stitched_x.begin(), stitched_x.end(), stitched_x.begin(), @@ -526,7 +568,8 @@ public: TS_ASSERT_EQUALS(temp, 0.); } // Check that the output X-Values are correct. - // truncate the input and oputput x values to 6 decimal places to eliminate + // truncate the input and oputput x values to 6 decimal places to + // eliminate // insignificant error auto xCopy = this->x; std::transform(stitched_x.begin(), stitched_x.end(), stitched_x.begin(), @@ -556,7 +599,8 @@ public: TS_ASSERT_EQUALS(temp, 0.); } // Check that the output X-Values are correct. - // truncate the input and oputput x values to 6 decimal places to eliminate + // truncate the input and oputput x values to 6 decimal places to + // eliminate // insignificant error auto xCopy = this->x; std::transform(stitched_x.begin(), stitched_x.end(), stitched_x.begin(), @@ -586,7 +630,8 @@ public: TS_ASSERT_EQUALS(temp, 0); } // Check that the output X-Values are correct. - // truncate the input and oputput x values to 6 decimal places to eliminate + // truncate the input and oputput x values to 6 decimal places to + // eliminate // insignificant error auto xCopy = this->x; std::transform(stitched_x.begin(), stitched_x.end(), stitched_x.begin(), @@ -732,7 +777,8 @@ public: class Stitch1DTestPerformance : public CxxTest::TestSuite { public: - // This pair of boilerplate methods prevent the suite being created statically + // This pair of boilerplate methods prevent the suite being created + // statically // This means the constructor isn't called when running other tests static Stitch1DTestPerformance *createSuite() { return new Stitch1DTestPerformance();