Merge pull request #17985 from mantidproject/17801_montecarlo_abs_cspline

Add CSpline interpolation option to MonteCarloAbsoprtion

Framework/Algorithms/CMakeLists.txt

@@ -18,9 +18,9 @@ set ( SRC_FILES
 	src/AverageLogData.cpp
 	src/BinaryOperateMasks.cpp
 	src/BinaryOperation.cpp
-	src/CalculateCountRate.cpp
 	src/CalMuonDeadTime.cpp
 	src/CalMuonDetectorPhases.cpp
+	src/CalculateCountRate.cpp
 	src/CalculateDIFC.cpp
 	src/CalculateEfficiency.cpp
 	src/CalculateFlatBackground.cpp
@@ -152,6 +152,7 @@ set ( SRC_FILES
 	src/IntegrateByComponent.cpp
 	src/Integration.cpp
 	src/InterpolatingRebin.cpp
+	src/InterpolationOption.cpp
 	src/InvertMask.cpp
 	src/Logarithm.cpp
 	src/LorentzCorrection.cpp
@@ -326,9 +327,9 @@ set ( INC_FILES
 	inc/MantidAlgorithms/BinaryOperateMasks.h
 	inc/MantidAlgorithms/BinaryOperation.h
 	inc/MantidAlgorithms/BoostOptionalToAlgorithmProperty.h
-	inc/MantidAlgorithms/CalculateCountRate.h
 	inc/MantidAlgorithms/CalMuonDeadTime.h
 	inc/MantidAlgorithms/CalMuonDetectorPhases.h
+	inc/MantidAlgorithms/CalculateCountRate.h
 	inc/MantidAlgorithms/CalculateDIFC.h
 	inc/MantidAlgorithms/CalculateEfficiency.h
 	inc/MantidAlgorithms/CalculateFlatBackground.h
@@ -461,6 +462,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/IntegrateByComponent.h
 	inc/MantidAlgorithms/Integration.h
 	inc/MantidAlgorithms/InterpolatingRebin.h
+	inc/MantidAlgorithms/InterpolationOption.h
 	inc/MantidAlgorithms/InvertMask.h
 	inc/MantidAlgorithms/Logarithm.h
 	inc/MantidAlgorithms/LorentzCorrection.h
@@ -646,9 +648,9 @@ set ( TEST_FILES
 	AverageLogDataTest.h
 	BinaryOperateMasksTest.h
 	BinaryOperationTest.h
-	CalculateCountRateTest.h
 	CalMuonDeadTimeTest.h
 	CalMuonDetectorPhasesTest.h
+	CalculateCountRateTest.h
 	CalculateDIFCTest.h
 	CalculateEfficiencyTest.h
 	CalculateFlatBackgroundTest.h
@@ -772,6 +774,7 @@ set ( TEST_FILES
 	IntegrateByComponentTest.h
 	IntegrationTest.h
 	InterpolatingRebinTest.h
+	InterpolationOptionTest.h
 	InvertMaskTest.h
 	LogarithmTest.h
 	LorentzCorrectionTest.h

Framework/Algorithms/inc/MantidAlgorithms/InterpolationOption.h

+#ifndef MANTID_ALGORITHMS_INTERPOLATIONOPTION_H_
+#define MANTID_ALGORITHMS_INTERPOLATIONOPTION_H_
+
+#include "MantidAlgorithms/DllConfig.h"
+#include <memory>
+
+namespace Mantid {
+namespace HistogramData {
+class Histogram;
+}
+namespace Kernel {
+class Property;
+}
+namespace Algorithms {
+
+/**
+  Class to provide a consistent interface to an interpolation option on
+  algorithms.
+
+  Copyright &copy; 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+  National Laboratory & European Spallation Source
+
+  This file is part of Mantid.
+
+  Mantid is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  Mantid is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+  File change history is stored at: <https://github.com/mantidproject/mantid>
+  Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+class MANTID_ALGORITHMS_DLL InterpolationOption {
+public:
+  // Indices must match the order in static string array
+  enum class Value : uint8_t { Linear, CSpline };
+
+  InterpolationOption();
+  void set(Value kind);
+  void set(const std::string &kind);
+
+  std::unique_ptr<Kernel::Property> property() const;
+  std::string propertyDoc() const;
+
+  void applyInplace(HistogramData::Histogram &inOut, size_t stepSize) const;
+
+private:
+  using MethodImpl = void(HistogramData::Histogram &inOut,
+                          const size_t stepSize);
+  MethodImpl *m_impl = nullptr;
+};
+
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif /* MANTID_ALGORITHMS_INTERPOLATIONOPTION_H_ */

Framework/Algorithms/inc/MantidAlgorithms/MonteCarloAbsorption.h

@@ -6,6 +6,7 @@
 //------------------------------------------------------------------------------
 #include "MantidAPI/Algorithm.h"
 #include "MantidAlgorithms/SampleCorrections/IBeamProfile.h"
+#include "MantidAlgorithms/InterpolationOption.h"
 
 namespace Mantid {
 namespace API {
@@ -61,8 +62,9 @@ private:
   void init() override;
   void exec() override;
 
-  API::MatrixWorkspace_sptr doSimulation(const API::MatrixWorkspace &inputWS,
-                                         size_t nevents, int nlambda, int seed);
+  API::MatrixWorkspace_sptr
+  doSimulation(const API::MatrixWorkspace &inputWS, size_t nevents, int nlambda,
+               int seed, const InterpolationOption &interpolateOpt);
   API::MatrixWorkspace_sptr
   createOutputWorkspace(const API::MatrixWorkspace &inputWS) const;
   std::unique_ptr<IBeamProfile>

Framework/Algorithms/src/InterpolationOption.cpp

+#include "MantidAlgorithms/InterpolationOption.h"
+#include "MantidHistogramData/Interpolate.h"
+#include "MantidKernel/ListValidator.h"
+#include "MantidKernel/make_unique.h"
+#include "MantidKernel/PropertyWithValue.h"
+
+#include <cassert>
+
+namespace {
+// The name of the interpolation property
+std::string PROP_NAME("Interpolation");
+std::string LINEAR_OPT("Linear");
+std::string CSPLINE_OPT("CSpline");
+std::vector<std::string> OPTIONS{LINEAR_OPT, CSPLINE_OPT};
+}
+
+namespace Mantid {
+using HistogramData::interpolateLinearInplace;
+using HistogramData::interpolateCSplineInplace;
+using Kernel::Property;
+
+namespace Algorithms {
+
+/**
+ * Constructs an object with the default interpolation method
+ */
+InterpolationOption::InterpolationOption() { set(Value::Linear); }
+
+/**
+ * Set the interpolation option
+ * @param kind Set the type of interpolation on the call to apply
+ */
+void InterpolationOption::set(InterpolationOption::Value kind) {
+  if (kind == Value::Linear) {
+    m_impl = interpolateLinearInplace;
+  } else {
+    m_impl = interpolateCSplineInplace;
+  }
+}
+
+/**
+ * Set the interpolation option
+ * @param kind Set the type of interpolation on the call to apply
+ */
+void InterpolationOption::set(const std::string &kind) {
+  if (kind == LINEAR_OPT) {
+    m_impl = interpolateLinearInplace;
+  } else if (kind == CSPLINE_OPT) {
+    m_impl = interpolateCSplineInplace;
+  } else {
+    throw std::invalid_argument(
+        "InterpolationOption::set() - Unknown interpolation method '" + kind +
+        "'");
+  }
+}
+
+/**
+ * Create a property suitable to attach to an algorithm to support interpolation
+ * @return A new Property containing the valid list of interpolation methods
+ */
+
+std::unique_ptr<Property> InterpolationOption::property() const {
+  using Kernel::StringListValidator;
+  using StringProperty = Kernel::PropertyWithValue<std::string>;
+
+  return Kernel::make_unique<StringProperty>(
+      PROP_NAME, LINEAR_OPT, boost::make_shared<StringListValidator>(OPTIONS));
+}
+
+/**
+ * @return The documentation string for the property
+ */
+std::string InterpolationOption::propertyDoc() const {
+  return "Method of interpolation used to compute unsimulated values.";
+}
+
+/**
+ * Apply the interpolation method to the given histogram
+ * @param inOut A reference to a histogram to interpolate
+ * @param stepSize The step size of calculated points
+ */
+void InterpolationOption::applyInplace(HistogramData::Histogram &inOut,
+                                       size_t stepSize) const {
+  assert(m_impl);
+  (*m_impl)(inOut, stepSize);
+}
+
+} // namespace Algorithms
+} // namespace Mantid

Framework/Algorithms/src/MonteCarloAbsorption.cpp

@@ -2,6 +2,7 @@
 // Includes
 //------------------------------------------------------------------------------
 #include "MantidAlgorithms/MonteCarloAbsorption.h"
+#include "MantidAlgorithms/InterpolationOption.h"
 #include "MantidAPI/ExperimentInfo.h"
 #include "MantidAPI/InstrumentValidator.h"
 #include "MantidAPI/Sample.h"
@@ -103,6 +104,9 @@ void MonteCarloAbsorption::init() {
       "The number of \"neutron\" events to generate per simulated point");
   declareProperty("SeedValue", DEFAULT_SEED, positiveInt,
                   "Seed the random number generator with this value");
+
+  InterpolationOption interpolateOpt;
+  declareProperty(interpolateOpt.property(), interpolateOpt.propertyDoc());
 }
 
 /**
@@ -113,9 +117,11 @@ void MonteCarloAbsorption::exec() {
   const int nevents = getProperty("EventsPerPoint");
   const int nlambda = getProperty("NumberOfWavelengthPoints");
   const int seed = getProperty("SeedValue");
+  InterpolationOption interpolateOpt;
+  interpolateOpt.set(getPropertyValue("Interpolation"));
 
-  auto outputWS =
-      doSimulation(*inputWS, static_cast<size_t>(nevents), nlambda, seed);
+  auto outputWS = doSimulation(*inputWS, static_cast<size_t>(nevents), nlambda,
+                               seed, interpolateOpt);
 
   setProperty("OutputWorkspace", outputWS);
 }
@@ -127,11 +133,13 @@ void MonteCarloAbsorption::exec() {
  * @param nlambda Number of wavelength points to simulate. The remainder
  * are computed using interpolation
  * @param seed Seed value for the random number generator
+ * @param interpolateOpt Method of interpolation to compute unsimulated points
  * @return A new workspace containing the correction factors & errors
  */
 MatrixWorkspace_sptr
 MonteCarloAbsorption::doSimulation(const MatrixWorkspace &inputWS,
-                                   size_t nevents, int nlambda, int seed) {
+                                   size_t nevents, int nlambda, int seed,
+                                   const InterpolationOption &interpolateOpt) {
   auto outputWS = createOutputWorkspace(inputWS);
   // Cache information about the workspace that will be used repeatedly
   auto instrument = inputWS.getInstrument();
@@ -203,7 +211,7 @@ MonteCarloAbsorption::doSimulation(const MatrixWorkspace &inputWS,
     // Interpolate through points not simulated
     if (lambdaStepSize > 1) {
       auto histnew = outputWS->histogram(i);
-      interpolateLinearInplace(histnew, lambdaStepSize);
+      interpolateOpt.applyInplace(histnew, lambdaStepSize);
       outputWS->setHistogram(i, histnew);
     }
 

Framework/Algorithms/test/InterpolationOptionTest.h

+#ifndef MANTID_ALGORITHMS_INTERPOLATIONOPTIONTEST_H_
+#define MANTID_ALGORITHMS_INTERPOLATIONOPTIONTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/InterpolationOption.h"
+#include "MantidHistogramData/Histogram.h"
+#include "MantidHistogramData/Points.h"
+#include "MantidHistogramData/LinearGenerator.h"
+#include "MantidKernel/PropertyWithValue.h"
+#include "MantidTestHelpers/HistogramDataTestHelper.h"
+
+using Mantid::Algorithms::InterpolationOption;
+using namespace Mantid::HistogramData;
+
+class InterpolationOptionTest : public CxxTest::TestSuite {
+public:
+  // This pair of boilerplate methods prevent the suite being created statically
+  // This means the constructor isn't called when running other tests
+  static InterpolationOptionTest *createSuite() {
+    return new InterpolationOptionTest();
+  }
+  static void destroySuite(InterpolationOptionTest *suite) { delete suite; }
+
+  //----------------------------------------------------------------------------
+  // Success tests
+  //----------------------------------------------------------------------------
+  void test_Property_Defaults_To_Linear_Interpolation() {
+    InterpolationOption interpolateOpt;
+    auto prop = interpolateOpt.property();
+
+    TS_ASSERT(prop);
+    TS_ASSERT_EQUALS("Interpolation", prop->name());
+    TS_ASSERT_EQUALS("Linear", prop->getDefault());
+  }
+
+  void test_Documentation_Is_Not_Empty() {
+    InterpolationOption interpolateOpt;
+
+    TS_ASSERT(!interpolateOpt.propertyDoc().empty())
+  }
+
+  void test_Apply_With_Linear_Succeeds() {
+    using namespace Mantid::HistogramData;
+    InterpolationOption interpolateOpt;
+
+    Histogram inOut(Points(7, LinearGenerator(0, 0.5)),
+                    Counts({-3, 0, -4, 0, 4, 0, 3}));
+    Histogram input(inOut);
+    interpolateOpt.applyInplace(inOut, 2);
+
+    const std::vector<double> expectedY = {-3, -3.5, -4, 0, 4, 3.5, 3};
+    checkData(input, inOut, expectedY);
+  }
+
+  void test_Apply_With_CSpline_Succeeds() {
+    InterpolationOption interpolateOptEnum;
+    // Set by enum
+    interpolateOptEnum.set(InterpolationOption::Value::CSpline);
+
+    Histogram inOut(Points(7, LinearGenerator(0, 0.5)),
+                    Counts({-3, 0, -4, 0, 4, 0, 3}));
+    const Histogram input(inOut);
+    interpolateOptEnum.applyInplace(inOut, 2);
+
+    const std::vector<double> expectedY = {-3, -4.625, -4, 0., 4, 4.625, 3};
+    checkData(input, inOut, expectedY);
+
+    // Set by string
+    InterpolationOption interpolateOptStr;
+    interpolateOptStr.set("CSpline");
+
+    Histogram inOutStr(input);
+    interpolateOptStr.applyInplace(inOutStr, 2);
+
+    checkData(input, inOutStr, expectedY);
+  }
+
+  //----------------------------------------------------------------------------
+  // Failure tests
+  //----------------------------------------------------------------------------
+  void test_set_From_String_Throws_With_Unknown_Type() {
+    InterpolationOption interpolateOpt;
+    TS_ASSERT_THROWS(interpolateOpt.set("Unknown"), std::invalid_argument);
+  }
+
+  void test_set_From_String_Throws_With_Empty_String() {
+    InterpolationOption interpolateOpt;
+    TS_ASSERT_THROWS(interpolateOpt.set(""), std::invalid_argument);
+  }
+
+private:
+  void checkData(const Histogram &input, const Histogram &output,
+                 const std::vector<double> &expectedY) {
+    TS_ASSERT_EQUALS(input.x(), output.x());
+    TS_ASSERT_EQUALS(input.xMode(), output.xMode());
+    TS_ASSERT_EQUALS(input.yMode(), output.yMode());
+    const auto &outY = output.y();
+    for (size_t i = 0; i < expectedY.size(); ++i) {
+      TS_ASSERT_DELTA(expectedY[i], outY[i], 1e-14);
+    }
+  }
+};
+
+#endif /* MANTID_ALGORITHMS_INTERPOLATIONOPTIONTEST_H_ */

Framework/Algorithms/test/MonteCarloAbsorptionTest.h

@@ -174,14 +174,30 @@ public:
     TestWorkspaceDescriptor wsProps = {1, 10, Environment::SampleOnly,
                                        DeltaEMode::Elastic, -1, -1};
     const int nlambda(5);
-    auto outputWS = runAlgorithm(wsProps, nlambda);
+    const std::string interpolation("Linear");
+    auto outputWS = runAlgorithm(wsProps, nlambda, interpolation);
 
     verifyDimensions(wsProps, outputWS);
     const double delta(1e-05);
-    const size_t middle_index(4);
+    TS_ASSERT_DELTA(0.019012, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.0044149, outputWS->y(0)[3], delta);
+    TS_ASSERT_DELTA(0.0026940, outputWS->y(0)[4], delta);
+    TS_ASSERT_DELTA(0.00042886, outputWS->y(0).back(), delta);
+  }
 
+  void test_CSpline_Interpolation() {
+    using Mantid::Kernel::DeltaEMode;
+    TestWorkspaceDescriptor wsProps = {1, 10, Environment::SampleOnly,
+                                       DeltaEMode::Elastic, -1, -1};
+    const int nlambda(5);
+    const std::string interpolation("CSpline");
+    auto outputWS = runAlgorithm(wsProps, nlambda, interpolation);
+
+    verifyDimensions(wsProps, outputWS);
+    const double delta(1e-05);
     TS_ASSERT_DELTA(0.019012, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.0026940, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.0036653, outputWS->y(0)[3], delta);
+    TS_ASSERT_DELTA(0.0026940, outputWS->y(0)[4], delta);
     TS_ASSERT_DELTA(0.00042886, outputWS->y(0).back(), delta);
   }
 
@@ -214,7 +230,8 @@ public:
 
 private:
   Mantid::API::MatrixWorkspace_const_sptr
-  runAlgorithm(const TestWorkspaceDescriptor &wsProps, int nlambda = -1) {
+  runAlgorithm(const TestWorkspaceDescriptor &wsProps, int nlambda = -1,
+               const std::string &interpolate = "") {
     auto inputWS = setUpWS(wsProps);
     auto mcabs = createAlgorithm();
     TS_ASSERT_THROWS_NOTHING(mcabs->setProperty("InputWorkspace", inputWS));
@@ -222,6 +239,10 @@ private:
       TS_ASSERT_THROWS_NOTHING(
           mcabs->setProperty("NumberOfWavelengthPoints", nlambda));
     }
+    if (!interpolate.empty()) {
+      TS_ASSERT_THROWS_NOTHING(
+          mcabs->setProperty("Interpolation", interpolate));
+    }
     mcabs->execute();
     return getOutputWorkspace(mcabs);
   }

docs/source/algorithms/MonteCarloAbsorption-v1.rst

@@ -77,7 +77,12 @@ The algorithm proceeds as follows. For each spectrum:
    * average the accumulated attentuation factors over `NEvents` and assign this as the correction factor for
      this :math:`\lambda_{step}`.
 
-#. finally, perform an interpolation through the unsimulated wavelength points
+#. finally, interpolate through the unsimulated wavelength points using the selected method
+
+Interpolation
+#############
+
+The default linear interpolation method will produce an absorption curve that is not smooth. CSpline interpolation will produce a smoother result by using a 3rd-order polynomial to approximate the original points. 
 
 Usage
 -----
@@ -96,6 +101,22 @@ Usage
    abscor = MonteCarloAbsorption(data, NumberOfWavelengthPoints=50)
    corrected = data/abscor
 
+**Example: A cylindrical sample with no container, interpolating with a CSpline**
+
+.. testcode:: ExCylinderSampleOnlyAndSpline
+
+   data = CreateSampleWorkspace(WorkspaceType='Histogram', NumBanks=1)
+   data = ConvertUnits(data, Target="Wavelength")
+   # Default up axis is Y
+   SetSample(data, Geometry={'Shape': 'Cylinder', 'Height': 5.0, 'Radius': 1.0,
+                     'Center': [0.0,0.0,0.0]},
+                   Material={'ChemicalFormula': '(Li7)2-C-H4-N-Cl6', 'SampleNumberDensity': 0.07})
+   # Simulating every data point can be slow. Use a smaller set and interpolate
+   abscor = MonteCarloAbsorption(data, NumberOfWavelengthPoints=50,
+                                 Interpolation='CSpline')
+   corrected = data/abscor
+
+
 **Example: A cylindrical sample setting a beam size**
 
 .. testcode:: ExCylinderSampleAndBeamSize

docs/source/release/v3.9.0/framework.rst

@@ -18,10 +18,12 @@ New
 Improved
 ########
 
-- :ref:`CalculateFlatBackground <algm-CalculateFlatBackground>` has now a new mode 'Moving Average' which takes the minimum of a moving window average as the flat background.
+- :ref:`CalculateFlatBackground <algm-CalculateFlatBackground>` has a new mode 'Moving Average' which takes the minimum of a moving window average as the flat background.
 - :ref:`StartLiveData <algm-StartLiveData>` and its dialog now support dynamic listener properties, based on the specific LiveListener being used.
 - :ref: All algorithms using `AsciiPointBase` now have a new property 'Separator' which allows the delimiter to be set to either comma, space or tab. This affects `SaveReflCustomAscii <algm-SaveReflCustomAscii>`, `SaveReflThreeColumnAscii <algm-SaveReflThreeColumnAscii>`, `SaveANSTOAscii <algm-SaveANSTOAscii>` and `SaveILLCosmosAscii <algm-SaveILLCosmosAscii>`.
 - :ref:`ReplaceSpecialValues <algm_ReplaceSpecialValues>` now can replace 'small' values below a user specified threshold.
+- :ref:`MonteCarloAbsorption <algm-MonteCarloAbsorption>` gained a new option: `Interpolation`.
+  This controls the method used for interpolation. Availabile options are: `Linear` & `CSpline`.
 
 Deprecated
 ##########