Merge pull request #475 from mantidproject/11384_use_estimate_peak_errors_in_poldi_2d_fir

Refactor PoldiFitPeaks2D to use EstimatePeakErrors

Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h

@@ -5,6 +5,10 @@
 #include "MantidSINQ/DllConfig.h"
 #include "MantidAPI/Algorithm.h"
 #include "MantidAPI/IFunction.h"
+#include "MantidAPI/IPeakFunction.h"
+
+#include "MantidKernel/Matrix.h"
+
 #include "MantidDataObjects/TableWorkspace.h"
 #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
 #include "MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h"
@@ -58,6 +62,10 @@ public:
 
   virtual const std::string summary() const;
 
+  boost::shared_ptr<Kernel::DblMatrix> getLocalCovarianceMatrix(
+      const boost::shared_ptr<const Kernel::DblMatrix> &covarianceMatrix,
+      size_t parameterOffset, size_t nParams) const;
+
 protected:
   PoldiPeakCollection_sptr
   getPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const;
@@ -72,7 +80,7 @@ protected:
                            PoldiPeakCollection_sptr &to) const;
 
   PoldiPeakCollection_sptr
-  getPeakCollectionFromFunction(const API::IFunction_sptr &fitFunction) const;
+  getPeakCollectionFromFunction(const API::IFunction_sptr &fitFunction);
   Poldi2DFunction_sptr getFunctionFromPeakCollection(
       const PoldiPeakCollection_sptr &peakCollection) const;
   void addBackgroundTerms(Poldi2DFunction_sptr poldi2DFunction) const;

Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp

@@ -50,9 +50,7 @@ const std::string PoldiFitPeaks2D::name() const { return "PoldiFitPeaks2D"; }
 int PoldiFitPeaks2D::version() const { return 1; }
 
 /// Algorithm's category for identification. @see Algorithm::category
-const std::string PoldiFitPeaks2D::category() const {
-  return "SINQ\\Poldi";
-}
+const std::string PoldiFitPeaks2D::category() const { return "SINQ\\Poldi"; }
 
 /// Very short algorithm summary. @see Algorith::summary
 const std::string PoldiFitPeaks2D::summary() const {
@@ -102,6 +100,37 @@ void PoldiFitPeaks2D::init() {
                   "Table workspace with fitted peaks.");
 }
 
+/**
+ * Extracts the covariance matrix for the supplied function
+ *
+ * This method extracts the covariance matrix for a sub-function from
+ * the global covariance matrix. If no matrix is given, a zero-matrix is
+ * returned.
+ *
+ * @param covarianceMatrix :: Global covariance matrix.
+ * @param parameterOffset :: Offset for the parameters of profileFunction.
+ * @param nParams :: Number of parameters of the local function.
+ * @return
+ */
+boost::shared_ptr<DblMatrix> PoldiFitPeaks2D::getLocalCovarianceMatrix(
+    const boost::shared_ptr<const Kernel::DblMatrix> &covarianceMatrix,
+    size_t parameterOffset, size_t nParams) const {
+  if (covarianceMatrix) {
+    boost::shared_ptr<Kernel::DblMatrix> localCov =
+        boost::make_shared<Kernel::DblMatrix>(nParams, nParams, false);
+    for (size_t j = 0; j < nParams; ++j) {
+      for (size_t k = 0; k < nParams; ++k) {
+        (*localCov)[j][k] =
+            (*covarianceMatrix)[parameterOffset + j][parameterOffset + k];
+      }
+    }
+
+    return localCov;
+  }
+
+  return boost::make_shared<Kernel::DblMatrix>(nParams, nParams, false);
+}
+
 /**
  * Construct a PoldiPeakCollection from a Poldi2DFunction
  *
@@ -114,7 +143,7 @@ void PoldiFitPeaks2D::init() {
  * @return PoldiPeakCollection containing peaks with normalized intensities
  */
 PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
-    const IFunction_sptr &fitFunction) const {
+    const IFunction_sptr &fitFunction) {
   Poldi2DFunction_sptr poldi2DFunction =
       boost::dynamic_pointer_cast<Poldi2DFunction>(fitFunction);
 
@@ -126,6 +155,11 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
   PoldiPeakCollection_sptr normalizedPeaks =
       boost::make_shared<PoldiPeakCollection>(PoldiPeakCollection::Integral);
 
+  boost::shared_ptr<const Kernel::DblMatrix> covarianceMatrix =
+      poldi2DFunction->getCovarianceMatrix();
+
+  size_t offset = 0;
+
   for (size_t i = 0; i < poldi2DFunction->nFunctions(); ++i) {
     boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
         boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
@@ -136,36 +170,35 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
           boost::dynamic_pointer_cast<IPeakFunction>(
               peakFunction->getProfileFunction());
 
+      // Get local covariance matrix
+      size_t nLocalParams = profileFunction->nParams();
+      boost::shared_ptr<Kernel::DblMatrix> localCov =
+          getLocalCovarianceMatrix(covarianceMatrix, offset, nLocalParams);
+      profileFunction->setCovarianceMatrix(localCov);
+
+      // Increment offset for next function
+      offset += nLocalParams;
+
+      IAlgorithm_sptr errorAlg = createChildAlgorithm("EstimatePeakErrors");
+      errorAlg->setProperty(
+          "Function", boost::dynamic_pointer_cast<IFunction>(profileFunction));
+      errorAlg->setPropertyValue("OutputWorkspace", "Errors");
+      errorAlg->execute();
+
       double centre = profileFunction->centre();
       double height = profileFunction->height();
+      double fwhmValue = profileFunction->fwhm();
 
-      size_t dIndex = 0;
-      size_t iIndex = 0;
-      size_t fIndex = 0;
-
-      for (size_t j = 0; j < profileFunction->nParams(); ++j) {
-        if (profileFunction->getParameter(j) == centre) {
-          dIndex = j;
-        } else if (profileFunction->getParameter(j) == height) {
-          iIndex = j;
-        } else {
-          fIndex = j;
-        }
-      }
-
-      // size_t dIndex = peakFunction->parameterIndex("Centre");
-      UncertainValue d(peakFunction->getParameter(dIndex),
-                       peakFunction->getError(dIndex));
+      ITableWorkspace_sptr errorTable =
+          errorAlg->getProperty("OutputWorkspace");
 
-      // size_t iIndex = peakFunction->parameterIndex("Area");
-      UncertainValue intensity(peakFunction->getParameter(iIndex),
-                               peakFunction->getError(iIndex));
+      double centreError = errorTable->cell<double>(0, 2);
+      double heightError = errorTable->cell<double>(1, 2);
+      double fwhmError = errorTable->cell<double>(2, 2);
 
-      // size_t fIndex = peakFunction->parameterIndex("Sigma");
-      double fwhmValue = profileFunction->fwhm();
-      UncertainValue fwhm(fwhmValue, fwhmValue /
-                                         peakFunction->getParameter(fIndex) *
-                                         peakFunction->getError(fIndex));
+      UncertainValue d(centre, centreError);
+      UncertainValue intensity(height, heightError);
+      UncertainValue fwhm(fwhmValue, fwhmError);
 
       PoldiPeak_sptr peak =
           PoldiPeak::create(MillerIndices(), d, intensity, UncertainValue(1.0));

Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h

@@ -6,6 +6,8 @@
 #include "MantidAPI/FrameworkManager.h"
 #include "MantidAPI/AlgorithmManager.h"
 
+#include "MantidKernel/Matrix.h"
+
 #include "MantidSINQ/PoldiFitPeaks2D.h"
 #include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
 #include "MantidSINQ/PoldiUtilities/PoldiMockInstrumentHelpers.h"
@@ -14,6 +16,7 @@
 using namespace Mantid::Poldi;
 using namespace Mantid::API;
 using namespace Mantid::DataObjects;
+using namespace Mantid::Kernel;
 
 class PoldiFitPeaks2DTest : public CxxTest::TestSuite
 {
@@ -234,6 +237,18 @@ public:
         PoldiPeakCollection_sptr peaks = PoldiPeakCollectionHelpers::createPoldiPeakCollectionNormalized();
         IFunction_sptr poldi2DFunction = spectrumCalculator.getFunctionFromPeakCollection(peaks);
 
+        size_t nParams = poldi2DFunction->nParams();
+
+        // Make a matrix with diagonal elements = 0.05 and set as covariance matrix
+        boost::shared_ptr<DblMatrix> matrix = boost::make_shared<DblMatrix>(nParams, nParams, true);
+        matrix->operator *=(0.05);
+        poldi2DFunction->setCovarianceMatrix(matrix);
+
+        // Also set errors for old behavior
+        for(size_t i = 0; i < nParams; ++i) {
+            poldi2DFunction->setError(i, sqrt(0.05));
+        }
+
         PoldiPeakCollection_sptr peaksFromFunction = spectrumCalculator.getPeakCollectionFromFunction(poldi2DFunction);
 
         TS_ASSERT_EQUALS(peaksFromFunction->peakCount(), peaks->peakCount());
@@ -241,8 +256,11 @@ public:
             PoldiPeak_sptr functionPeak = peaksFromFunction->peak(i);
             PoldiPeak_sptr referencePeak = peaks->peak(i);
 
-            TS_ASSERT_EQUALS(functionPeak->d(), referencePeak->d());
-            TS_ASSERT_EQUALS(functionPeak->fwhm(), referencePeak->fwhm());
+            TS_ASSERT_EQUALS(functionPeak->d().value(), referencePeak->d().value());
+            TS_ASSERT_EQUALS(functionPeak->fwhm().value(), referencePeak->fwhm().value());
+            TS_ASSERT_DELTA(functionPeak->d().error(), sqrt(0.05), 1e-6);
+            TS_ASSERT_DELTA(functionPeak->fwhm(PoldiPeak::AbsoluteD).error(), sqrt(0.05) * (2.0 * sqrt(2.0 * log(2.0))), 1e-6);
+
         }
     }