diff --git a/Framework/Algorithms/src/SetUncertainties.cpp b/Framework/Algorithms/src/SetUncertainties.cpp index a7aefc2ae5e3ed1903ee4146c376371aa29e472e..c652e851177251d3b3866023630a50a512c2cea8 100644 --- a/Framework/Algorithms/src/SetUncertainties.cpp +++ b/Framework/Algorithms/src/SetUncertainties.cpp @@ -8,6 +8,7 @@ #include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/SpectrumInfo.h" #include "MantidAPI/WorkspaceFactory.h" +#include "MantidDataObjects/EventWorkspace.h" #include "MantidGeometry/IDetector.h" #include "MantidGeometry/Instrument.h" #include "MantidKernel/BoundedValidator.h" @@ -111,6 +112,10 @@ void SetUncertainties::init() { void SetUncertainties::exec() { MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace"); + auto inputEventWorkspace = + boost::dynamic_pointer_cast<const DataObjects::EventWorkspace>( + inputWorkspace); + MatrixWorkspace_sptr outputWorkspace = getProperty("OutputWorkspace"); std::string errorType = getProperty("SetError"); bool zeroError = (errorType == ZERO); bool takeSqrt = ((errorType == SQRT) || (errorType == SQRT_OR_ONE)); @@ -120,32 +125,45 @@ void SetUncertainties::exec() { double valueToSet = resetOne ? 1.0 : getProperty("SetErrorTo"); double valueToCompare = resetOne ? 0.0 : getProperty("IfEqualTo"); int precision = getProperty("Precision"); - double tolerance = resetOne ? 1E-10 : std::pow(10.0, precision * (-1)); + double tolerance = resetOne ? 1E-10 : std::pow(10.0, -1. * precision); // Create the output workspace. This will copy many aspects from the input // one. - MatrixWorkspace_sptr outputWorkspace = - WorkspaceFactory::Instance().create(inputWorkspace); - // ...but not the data, so do that here. + const int64_t numHists = + static_cast<int64_t>(inputWorkspace->getNumberHistograms()); + if (inputEventWorkspace) { + if (inputWorkspace == outputWorkspace && errorType == SQRT && + inputEventWorkspace->getEventType() == Mantid::API::EventType::TOF) { + // Uncertainty is already square root of the y value + return; + } + // Copy the histogram representation over to a Workspace2D + outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace); + PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace)) + for (int64_t i = 0; i < numHists; ++i) { + PARALLEL_START_INTERUPT_REGION + // copy X/Y/E + outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i)); + outputWorkspace->setSharedY(i, inputWorkspace->sharedY(i)); + outputWorkspace->setSharedE(i, inputWorkspace->sharedE(i)); + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + } else if (inputWorkspace != outputWorkspace) { + outputWorkspace = inputWorkspace->clone(); + } + const auto &spectrumInfo = inputWorkspace->spectrumInfo(); - const size_t numHists = inputWorkspace->getNumberHistograms(); Progress prog(this, 0.0, 1.0, numHists); - PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace)) - for (int64_t i = 0; i < int64_t(numHists); ++i) { + for (int64_t i = 0; i < numHists; ++i) { PARALLEL_START_INTERUPT_REGION - - // copy the X/Y - outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i)); - outputWorkspace->setSharedY(i, inputWorkspace->sharedY(i)); // copy the E or set to zero depending on the mode if (errorType == ONE_IF_ZERO || customError) { - outputWorkspace->setSharedE(i, inputWorkspace->sharedE(i)); } else { outputWorkspace->mutableE(i) = 0.0; } - // ZERO mode doesn't calculate anything further if ((!zeroError) && (!(spectrumInfo.hasDetectors(i) && spectrumInfo.isMasked(i)))) { @@ -159,13 +177,10 @@ void SetUncertainties::exec() { SetError(valueToSet, valueToCompare, tolerance)); } } - prog.report(); - PARALLEL_END_INTERUPT_REGION } PARALLEL_CHECK_INTERUPT_REGION - setProperty("OutputWorkspace", outputWorkspace); } diff --git a/Framework/Algorithms/test/SetUncertaintiesTest.h b/Framework/Algorithms/test/SetUncertaintiesTest.h index c9f1cbab774b9923156a2f35c3f4e55c35138600..02ac7f9fbac25d88d2d3f24a0727cb857d9c0b99 100644 --- a/Framework/Algorithms/test/SetUncertaintiesTest.h +++ b/Framework/Algorithms/test/SetUncertaintiesTest.h @@ -55,6 +55,38 @@ public: return outWS; } + /** + * Create and execute the algorithm in the specified mode. + * @param mode The name of the SetError property SetUncertainties + * @return The name that the output workspace will be registered as. + */ + API::MatrixWorkspace_sptr runAlgInPlace(const std::string &mode) { + // random data mostly works + auto inWksp = WorkspaceCreationHelper::create1DWorkspaceRand(30, true); + // Ensure first elements of random workspace are zero so test don't + // pass randomly + auto &E = inWksp->mutableE(0); + E[0] = 0.; + auto &Y = inWksp->mutableY(0); + Y[1] = 0.; // stress sqrtOrOne + + std::string outWSname = "SetUncertainties_" + mode; + WorkspaceCreationHelper::storeWS(outWSname, inWksp); + SetUncertainties alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()); + alg.setProperty("InputWorkspace", inWksp); + alg.setProperty("SetError", mode); + alg.setProperty("OutputWorkspace", outWSname); + TS_ASSERT_THROWS_NOTHING(alg.execute()); + TS_ASSERT(alg.isExecuted()); + + const auto outWS = + API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>( + outWSname); + TS_ASSERT(bool(outWS)); // non-null pointer + return outWS; + } + API::MatrixWorkspace_sptr runAlgCustom(const double toSet, const double toReplace, const double errorVal, const int precision, const int position) { @@ -97,6 +129,17 @@ public: API::AnalysisDataService::Instance().remove(outWS->getName()); } + void test_zero_in_place() { + const auto outWS = runAlgInPlace("zero"); + + const auto &E = outWS->e(0); + for (const auto item : E) { + TS_ASSERT_EQUALS(item, 0.); + } + + API::AnalysisDataService::Instance().remove(outWS->getName()); + } + void test_sqrt() { const auto outWS = runAlg("sqrt"); @@ -109,6 +152,18 @@ public: API::AnalysisDataService::Instance().remove(outWS->getName()); } + void test_sqrt_in_place() { + const auto outWS = runAlgInPlace("sqrt"); + + const auto &E = outWS->e(0); + const auto &Y = outWS->y(0); + for (size_t i = 0; i < E.size(); ++i) { + TS_ASSERT_DELTA(Y[i], E[i] * E[i], .001); + } + + API::AnalysisDataService::Instance().remove(outWS->getName()); + } + void test_oneIfZero() { const auto outWS = runAlg("oneIfZero"); @@ -119,6 +174,16 @@ public: API::AnalysisDataService::Instance().remove(outWS->getName()); } + void test_oneIfZero_in_place() { + const auto outWS = runAlgInPlace("oneIfZero"); + + const auto &E = outWS->e(0); + for (const auto item : E) { + TS_ASSERT(item > 0.); + } + API::AnalysisDataService::Instance().remove(outWS->getName()); + } + void test_sqrtOrOne() { const auto outWS = runAlg("sqrtOrOne"); @@ -131,7 +196,21 @@ public: TS_ASSERT_DELTA(Y[i], E[i] * E[i], .001); } } + API::AnalysisDataService::Instance().remove(outWS->getName()); + } + + void test_sqrtOrOne_in_place() { + const auto outWS = runAlgInPlace("sqrtOrOne"); + const auto &E = outWS->e(0); + const auto &Y = outWS->y(0); + for (size_t i = 0; i < E.size(); ++i) { + if (Y[i] == 0.) { + TS_ASSERT_EQUALS(E[i], 1.); + } else { + TS_ASSERT_DELTA(Y[i], E[i] * E[i], .001); + } + } API::AnalysisDataService::Instance().remove(outWS->getName()); } @@ -201,6 +280,42 @@ public: } API::AnalysisDataService::Instance().remove(outWS->getName()); } + + /** + * Create and execute the algorithm in the specified mode. + * @param mode The name of the SetError property SetUncertainties + * @return The name that the output workspace will be registered as. + */ + void test_EventWorkspacePlace() { + // random data mostly works + auto inWksp = + WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument2(10, + 10); + std::string outWSname = "SetUncertainties_oneIfZero"; + WorkspaceCreationHelper::storeWS(outWSname, inWksp); + SetUncertainties alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()); + alg.setProperty("InputWorkspace", inWksp); + alg.setProperty("SetError", "oneIfZero"); + alg.setProperty("OutputWorkspace", outWSname); + TS_ASSERT_THROWS_NOTHING(alg.execute()); + TS_ASSERT(alg.isExecuted()); + + const auto outWS = + API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>( + outWSname); + TS_ASSERT(bool(outWS)); // non-null pointer + + const auto &E = outWS->e(0); + const auto &Y = outWS->y(0); + for (size_t i = 0; i < E.size(); ++i) { + if (Y[i] == 0.) { + TS_ASSERT_EQUALS(E[i], 1.); + } else { + TS_ASSERT_DELTA(Y[i], E[i] * E[i], .001); + } + } + } }; class SetUncertaintiesTestPerformance : public CxxTest::TestSuite { diff --git a/docs/source/release/v5.1.0/framework.rst b/docs/source/release/v5.1.0/framework.rst index a893818aae15fa7ccbaa5cd3a782ec694d43d5f4..02d64bb3567e0b6ba1ea322fbcb7eca5d6cf95a1 100644 --- a/docs/source/release/v5.1.0/framework.rst +++ b/docs/source/release/v5.1.0/framework.rst @@ -15,6 +15,10 @@ Concepts Algorithms ---------- +- Add specialization to :ref:`SetUncertainties <algm-SetUncertainties>` for the + case where InputWorkspace == OutputWorkspace. Where possible, avoid the + cost of cloning the inputWorkspace. + Data Objects ------------