Merge pull request #28302 from mantidproject/set_uncertainties_inplace_2

inplace specialization for SetUncertainties

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);
 }
 

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 {

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
 ------------