Merge remote-tracking branch 'origin/20898_unweighted_CalculatePolynomialBackground'

Framework/Algorithms/src/CalculatePolynomialBackground.cpp

@@ -9,18 +9,26 @@
 #include "MantidKernel/ArrayOrderedPairsValidator.h"
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/ListValidator.h"
 
 #include <utility>
 
 namespace {
 /// String constants for algorithm's properties.
 namespace Prop {
+const static std::string COST_FUNCTION = "CostFunction";
 const static std::string INPUT_WS = "InputWorkspace";
 const static std::string OUTPUT_WS = "OutputWorkspace";
 const static std::string POLY_DEGREE = "Degree";
 const static std::string XRANGES = "XRanges";
 }
 
+/// String constants for cost function options.
+namespace CostFunc {
+const static std::string UNWEIGHTED_LEAST_SQUARES = "Unweighted least squares";
+const static std::string WEIGHTED_LEAST_SQUARES = "Least squares";
+}
+
 /** Filters ranges completely outside the histogram X values.
  *  @param ranges a vector of start-end pairs to filter
  *  @param ws a workspace
@@ -155,11 +163,10 @@ std::vector<double> invertRanges(const std::vector<double> &ranges) {
  *  @param ranges a vector defining the fitting intervals
  *  @return a vector of final fitted parameters
  */
-std::vector<double> executeFit(Mantid::API::Algorithm &fit,
-                               const std::string &function,
-                               Mantid::API::MatrixWorkspace_sptr &ws,
-                               const size_t wsIndex,
-                               const std::vector<double> &ranges) {
+std::vector<double>
+executeFit(Mantid::API::Algorithm &fit, const std::string &function,
+           Mantid::API::MatrixWorkspace_sptr &ws, const size_t wsIndex,
+           const std::vector<double> &ranges, const std::string &costFunction) {
   const auto fitRanges = histogramRanges(ranges, *ws, wsIndex);
   const auto excludedRanges = invertRanges(fitRanges);
   fit.setProperty("Function", function);
@@ -168,6 +175,8 @@ std::vector<double> executeFit(Mantid::API::Algorithm &fit,
   fit.setProperty("StartX", fitRanges.front());
   fit.setProperty("EndX", fitRanges.back());
   fit.setProperty("Exclude", excludedRanges);
+  fit.setProperty("Minimizer", "Levenberg-MarquardtMD");
+  fit.setProperty(Prop::COST_FUNCTION, costFunction);
   fit.setProperty("CreateOutput", true);
   fit.executeAsChildAlg();
   Mantid::API::ITableWorkspace_sptr fitResult =
@@ -268,6 +277,12 @@ void CalculatePolynomialBackground::init() {
   declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
                       Prop::XRANGES, std::vector<double>(), orderedPairs),
                   "A list of fitting ranges given as pairs of X values.");
+  std::array<std::string, 2> costFuncOpts{
+      {CostFunc::WEIGHTED_LEAST_SQUARES, CostFunc::UNWEIGHTED_LEAST_SQUARES}};
+  declareProperty(
+      Prop::COST_FUNCTION, CostFunc::WEIGHTED_LEAST_SQUARES.c_str(),
+      boost::make_shared<Kernel::ListValidator<std::string>>(costFuncOpts),
+      "The cost function to be passed to the Fit algorithm.");
 }
 
 //----------------------------------------------------------------------------------------------
@@ -279,6 +294,7 @@ void CalculatePolynomialBackground::exec() {
   API::MatrixWorkspace_sptr outWS{
       DataObjects::create<DataObjects::Workspace2D>(*inWS)};
   const std::vector<double> inputRanges = getProperty(Prop::XRANGES);
+  const std::string costFunction = getProperty(Prop::COST_FUNCTION);
   const auto polyDegree =
       static_cast<size_t>(static_cast<int>(getProperty(Prop::POLY_DEGREE)));
   const std::vector<double> initialParams(polyDegree + 1, 0.1);
@@ -290,7 +306,8 @@ void CalculatePolynomialBackground::exec() {
     PARALLEL_START_INTERUPT_REGION
     const bool logging{false};
     auto fit = createChildAlgorithm("Fit", 0, 0, logging);
-    const auto parameters = executeFit(*fit, fitFunction, inWS, i, inputRanges);
+    const auto parameters =
+        executeFit(*fit, fitFunction, inWS, i, inputRanges, costFunction);
     const auto bkgFunction = makeFunctionString(parameters);
     evaluateInPlace(bkgFunction, *outWS, i);
     progress.report();

Framework/Algorithms/test/CalculatePolynomialBackgroundTest.h

@@ -98,7 +98,73 @@ public:
       const auto &bkgEs = outWS->e(histI);
       const auto &bkgXs = outWS->x(histI);
       for (size_t binI = 0; binI < nBin; ++binI) {
-        TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-12)
+        TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-10)
+        TS_ASSERT_EQUALS(bkgEs[binI], 0)
+        TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
+      }
+    }
+  }
+
+  void test_costFuctionLeastSquares() {
+    using namespace WorkspaceCreationHelper;
+    const size_t nHist{2};
+    const HistogramData::Counts counts{0, 4, 0, 0};
+    const HistogramData::CountStandardDeviations stdDevs{0, 0.001, 0, 0};
+    const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
+    API::MatrixWorkspace_sptr ws(
+        DataObjects::create<DataObjects::Workspace2D>(
+            nHist, HistogramData::Histogram(edges, counts, stdDevs)).release());
+    auto alg = makeAlgorithm();
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("CostFunction", "Least squares"))
+    TS_ASSERT_THROWS_NOTHING(alg->execute())
+    TS_ASSERT(alg->isExecuted())
+    API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
+    TS_ASSERT(outWS)
+    for (size_t histI = 0; histI < nHist; ++histI) {
+      const auto &xs = ws->x(histI);
+      const auto &bkgYs = outWS->y(histI);
+      const auto &bkgEs = outWS->e(histI);
+      const auto &bkgXs = outWS->x(histI);
+      for (size_t binI = 0; binI < counts.size(); ++binI) {
+        // Number 4 in counts is heavily weighted by the small error.
+        TS_ASSERT_DELTA(bkgYs[binI], 4, 1e-5)
+        TS_ASSERT_EQUALS(bkgEs[binI], 0)
+        TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
+      }
+    }
+  }
+
+  void test_costFuctionUnweightedLeastSquares() {
+    using namespace WorkspaceCreationHelper;
+    const size_t nHist{2};
+    const HistogramData::Counts counts{0, 4, 0, 0};
+    const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
+    API::MatrixWorkspace_sptr ws(
+        DataObjects::create<DataObjects::Workspace2D>(
+            nHist, HistogramData::Histogram(edges, counts)).release());
+    auto alg = makeAlgorithm();
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
+    TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
+    TS_ASSERT_THROWS_NOTHING(
+        alg->setProperty("CostFunction", "Unweighted least squares"))
+    TS_ASSERT_THROWS_NOTHING(alg->execute())
+    TS_ASSERT(alg->isExecuted())
+    API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
+    TS_ASSERT(outWS)
+    // Unweighted fitting of zeroth order polynomial is equivalent to the mean.
+    const double result = std::accumulate(counts.cbegin(), counts.cend(), 0.0) /
+                          static_cast<double>(counts.size());
+    for (size_t histI = 0; histI < nHist; ++histI) {
+      const auto &xs = ws->x(histI);
+      const auto &bkgYs = outWS->y(histI);
+      const auto &bkgEs = outWS->e(histI);
+      const auto &bkgXs = outWS->x(histI);
+      for (size_t binI = 0; binI < counts.size(); ++binI) {
+        TS_ASSERT_DELTA(bkgYs[binI], result, 1e-5)
         TS_ASSERT_EQUALS(bkgEs[binI], 0)
         TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
       }
@@ -176,7 +242,7 @@ public:
     const auto &bkgXs = outWS->x(0);
     const std::vector<double> expected{1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
     for (size_t binI = 0; binI < nBin; ++binI) {
-      TS_ASSERT_DELTA(bkgYs[binI], expected[binI], 1e-12)
+      TS_ASSERT_DELTA(bkgYs[binI], expected[binI], 1e-10)
       TS_ASSERT_EQUALS(bkgEs[binI], 0)
       TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
     }

docs/source/algorithms/CalculatePolynomialBackground-v1.rst

@@ -12,6 +12,8 @@ Description
 
 This algorithm calculates backgrounds for the histograms in a workspace by fitting a polynomial to ranges given by *XRanges*. The fitting is done using :ref:`algm-Fit`. The backgrounds are returned as the *OutputWorkspace* which can be subsequently subtracted from the original workspace. The degree of the polynomial can be given by *Degree*. A zeroth degree polynomial would correspond to fitting a constant value. The *XRanges* property is a list of range start and range end value pairs in the X axis units. Overlapping ranges are merged. If no *XRanges* are given, full histograms are used.
 
+The value of *CostFunction* is passed to :ref:`algm-Fit` as-is. The default option is 'Least squares' which uses the histogram errors as weights. This might not be desirable, e.g. when there are bins with zero counts and zero errors. An 'Unweighted least squares' option is available to deal with such cases.
+
 Usage
 -----
 

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

@@ -23,6 +23,7 @@ Algorithms
 ----------
 
 :ref:`NormaliseToMonitor <algm-NormaliseToMonitor>` now supports workspaces with detector scans and workspaces with single-count point data.
+- It is now possible to choose between weighted and unweighted fitting in :ref:`CalculatePolynomialBackground <algm-CalculatePolynomialBackground>`.
 - :ref:`CreateWorkspace <algm-CreateWorkspace>` will no longer create a default (and potentially wrong) mapping from spectra to detectors, unless a parent workspace is given. This change ensures that accidental bad mappings that could lead to corrupted data are not created silently anymore. This change does *not* affect the use of this algorithm if: (1) a parent workspace is given, or (2) no instrument is loaded into to workspace at a later point, or (3) an instrument is loaded at a later point but ``LoadInstrument`` is used with ``RewriteSpectraMapping=True``. See also the algorithm documentation for details.
 - :ref:`ConjoinWorkspaces <algm-ConjoinWorkspaces>` now supports non-constant bins.
 - :ref:`Fit <algm-Fit>` will now respect excluded ranges when ``CostFunction = 'Unweighted least squares'``.