Refs #9992. Checkpointing work

Results look promising already, need some fine tuning for input/output etc.

Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp

@@ -267,9 +267,8 @@ DECLARE_FUNCTION(PawleyFunction)
 
 /// Constructor
 PawleyFunction::PawleyFunction()
-    : IPawleyFunction(), m_compositeFunction(),
-      m_pawleyParameterFunction(), m_peakProfileComposite(), m_hkls(),
-      m_dUnit(), m_wsUnit() {}
+    : IPawleyFunction(), m_compositeFunction(), m_pawleyParameterFunction(),
+      m_peakProfileComposite(), m_hkls(), m_dUnit(), m_wsUnit() {}
 
 void PawleyFunction::setMatrixWorkspace(
     boost::shared_ptr<const MatrixWorkspace> workspace, size_t wi,
@@ -366,16 +365,24 @@ void PawleyFunction::function(const FunctionDomain &domain,
                               FunctionValues &values) const {
   UnitCell cell = m_pawleyParameterFunction->getUnitCellFromParameters();
   double zeroShift = m_pawleyParameterFunction->getParameter("ZeroShift");
+  std::string centreName =
+      m_pawleyParameterFunction->getProfileFunctionCenterParameterName();
 
   for (size_t i = 0; i < m_hkls.size(); ++i) {
     double centre = getTransformedCenter(cell.d(m_hkls[i]));
 
-    m_peakProfileComposite->getFunction(i)->setParameter(
-        m_pawleyParameterFunction->getProfileFunctionCenterParameterName(),
-        centre + zeroShift);
+    m_peakProfileComposite->getFunction(i)
+        ->setParameter(centreName, centre + zeroShift);
   }
 
   m_peakProfileComposite->function(domain, values);
+
+  for (size_t i = 0; i < m_hkls.size(); ++i) {
+    m_peakProfileComposite->getFunction(i)
+        ->setParameter(centreName, m_peakProfileComposite->getFunction(i)
+                                           ->getParameter(centreName) -
+                                       zeroShift);
+  }
 }
 
 /// Removes all peaks from the function.

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

@@ -5,6 +5,7 @@
 #include "MantidSINQ/DllConfig.h"
 #include "MantidAPI/Algorithm.h"
 #include "MantidAPI/IFunction.h"
+#include "MantidAPI/IPeakFunction.h"
 #include "MantidDataObjects/TableWorkspace.h"
 #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
 #include "MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h"
@@ -58,6 +59,18 @@ class MANTID_SINQ_DLL PoldiFitPeaks2D : public API::Algorithm {
 
   virtual const std::string summary() const;
 
+  Poldi2DFunction_sptr getFunctionIndividualPeaks(
+      std::string profileFunctionName,
+      const PoldiPeakCollection_sptr &peakCollection) const;
+
+  Poldi2DFunction_sptr
+  getFunctionPawley(std::string profileFunctionName,
+                    const PoldiPeakCollection_sptr &peakCollection) const;
+
+  PoldiPeak_sptr
+  getPeakFromPeakFunction(API::IPeakFunction_sptr profileFunction,
+                          const Kernel::V3D &hkl) const;
+
 protected:
   PoldiPeakCollection_sptr
   getPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const;

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

@@ -12,13 +12,15 @@ use the Build/wiki_maker.py script to generate your full wiki page.
 #include "MantidAPI/MultiDomainFunction.h"
 #include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
 #include "MantidSINQ/PoldiUtilities/PoldiSpectrumLinearBackground.h"
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h"
 #include "MantidAPI/FunctionDomain1D.h"
 
 #include "MantidSINQ/PoldiUtilities/IPoldiFunction1D.h"
 #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
 #include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
 #include "MantidSINQ/PoldiUtilities/PoldiDeadWireDecorator.h"
 #include "MantidAPI/IPeakFunction.h"
+#include "MantidAPI/IPawleyFunction.h"
 
 #include "MantidSINQ/PoldiUtilities/Poldi2DFunction.h"
 
@@ -50,9 +52,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 {
@@ -71,6 +71,17 @@ void PoldiFitPeaks2D::init() {
                   "Profile function to use for integrating the peak profiles "
                   "before calculating the spectrum.");
 
+  declareProperty("PawleyFit", false, "Instead of refining individual peaks, "
+                                      "refine a unit cell. Peaks must be "
+                                      "indexed, an initial cell must be "
+                                      "provided, as well as a crystal system.");
+  declareProperty("InitialCell", "", "Initial unit cell parameters as 6 "
+                                     "space-separated numbers. Only used when "
+                                     "PawleyFit is checked.");
+  declareProperty("CrystalSystem", "", "Crystal system to use for constraining "
+                                       "unit cell parameters. Only used when "
+                                       "PawleyFit is checked.");
+
   declareProperty("FitConstantBackground", true,
                   "Add a constant background term to the fit.");
   declareProperty("ConstantBackgroundParameter", 0.0,
@@ -102,6 +113,47 @@ void PoldiFitPeaks2D::init() {
                   "Table workspace with fitted peaks.");
 }
 
+PoldiPeak_sptr
+PoldiFitPeaks2D::getPeakFromPeakFunction(IPeakFunction_sptr profileFunction,
+                                         const V3D &hkl) const {
+  double centre = profileFunction->centre();
+  double height = profileFunction->height();
+
+  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(profileFunction->getParameter(dIndex),
+                   profileFunction->getError(dIndex));
+
+  // size_t iIndex = peakFunction->parameterIndex("Area");
+  UncertainValue intensity(profileFunction->getParameter(iIndex),
+                           profileFunction->getError(iIndex));
+
+  // size_t fIndex = peakFunction->parameterIndex("Sigma");
+  double fwhmValue = profileFunction->fwhm();
+  UncertainValue fwhm(fwhmValue, fwhmValue /
+                                     profileFunction->getParameter(fIndex) *
+                                     profileFunction->getError(fIndex));
+
+  PoldiPeak_sptr peak =
+      PoldiPeak::create(MillerIndices(hkl), d, intensity, UncertainValue(1.0));
+  peak->setFwhm(fwhm, PoldiPeak::FwhmRelation::AbsoluteD);
+
+  return peak;
+}
+
 /**
  * Construct a PoldiPeakCollection from a Poldi2DFunction
  *
@@ -127,77 +179,72 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
       boost::make_shared<PoldiPeakCollection>(PoldiPeakCollection::Integral);
 
   for (size_t i = 0; i < poldi2DFunction->nFunctions(); ++i) {
-    boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
-        boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
+    boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction =
+        boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>(
             poldi2DFunction->getFunction(i));
 
-    if (peakFunction) {
-      IPeakFunction_sptr profileFunction =
-          boost::dynamic_pointer_cast<IPeakFunction>(
-              peakFunction->getProfileFunction());
+    if (poldiPawleyFunction) {
+      IPawleyFunction_sptr pawleyFunction =
+          poldiPawleyFunction->getPawleyFunction();
+      if (pawleyFunction) {
+        CompositeFunction_sptr wrappedFn =
+            boost::dynamic_pointer_cast<CompositeFunction>(
+                pawleyFunction->getDecoratedFunction());
+        IFunction_sptr pawleyParamFn = wrappedFn->getFunction(0);
+        for (size_t j = 0; j < pawleyParamFn->nParams(); ++j) {
+          std::cout << j << " " << pawleyParamFn->parameterName(j) << " "
+                    << pawleyParamFn->getParameter(j) << " "
+                    << pawleyParamFn->getError(j) << std::endl;
+        }
+
+        for (size_t j = 0; j < pawleyFunction->getPeakCount(); ++j) {
+          IPeakFunction_sptr profileFunction =
+              pawleyFunction->getPeakFunction(j);
+          V3D peakHKL = pawleyFunction->getPeakHKL(j);
 
-      double centre = profileFunction->centre();
-      double height = profileFunction->height();
+          PoldiPeak_sptr peak =
+              getPeakFromPeakFunction(profileFunction, peakHKL);
 
-      size_t dIndex = 0;
-      size_t iIndex = 0;
-      size_t fIndex = 0;
+          normalizedPeaks->addPeak(peak);
 
-      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;
+          for (size_t p = 0; p < profileFunction->nParams(); ++p) {
+            std::cout << j << " " << p << " "
+                      << profileFunction->parameterName(p) << " "
+                      << profileFunction->getParameter(p) << std::endl;
+          }
         }
       }
+      break;
+    }
 
-      // size_t dIndex = peakFunction->parameterIndex("Centre");
-      UncertainValue d(peakFunction->getParameter(dIndex),
-                       peakFunction->getError(dIndex));
-
-      // size_t iIndex = peakFunction->parameterIndex("Area");
-      UncertainValue intensity(peakFunction->getParameter(iIndex),
-                               peakFunction->getError(iIndex));
+    boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
+        boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
+            poldi2DFunction->getFunction(i));
 
-      // size_t fIndex = peakFunction->parameterIndex("Sigma");
-      double fwhmValue = profileFunction->fwhm();
-      UncertainValue fwhm(fwhmValue, fwhmValue /
-                                         peakFunction->getParameter(fIndex) *
-                                         peakFunction->getError(fIndex));
+    if (peakFunction) {
+      IPeakFunction_sptr profileFunction =
+          boost::dynamic_pointer_cast<IPeakFunction>(
+              peakFunction->getProfileFunction());
 
       PoldiPeak_sptr peak =
-          PoldiPeak::create(MillerIndices(), d, intensity, UncertainValue(1.0));
-      peak->setFwhm(fwhm, PoldiPeak::FwhmRelation::AbsoluteD);
+          getPeakFromPeakFunction(profileFunction, V3D(0, 0, 0));
 
       normalizedPeaks->addPeak(peak);
+      continue;
     }
   }
 
   return normalizedPeaks;
 }
 
-/**
- * Constructs a proper function from a peak collection
- *
- * This method constructs a Poldi2DFunction and assigns one
- *PoldiSpectrumDomainFunction to it for
- * each peak contained in the peak collection.
- *
- * @param peakCollection :: PoldiPeakCollection containing peaks with integral
- *intensities
- * @return Poldi2DFunction with one PoldiSpectrumDomainFunction per peak
- */
-Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionIndividualPeaks(
+    std::string profileFunctionName,
     const PoldiPeakCollection_sptr &peakCollection) const {
   Poldi2DFunction_sptr mdFunction(new Poldi2DFunction);
 
   for (size_t i = 0; i < peakCollection->peakCount(); ++i) {
     PoldiPeak_sptr peak = peakCollection->peak(i);
 
-    std::string profileFunctionName = getProperty("PeakProfileFunction");
-
     boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
         boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
             FunctionFactory::Instance().createFunction(
@@ -226,6 +273,86 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
   return mdFunction;
 }
 
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley(
+    std::string profileFunctionName,
+    const PoldiPeakCollection_sptr &peakCollection) const {
+  Poldi2DFunction_sptr mdFunction(new Poldi2DFunction);
+
+  boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction =
+      boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>(
+          FunctionFactory::Instance().createFunction(
+              "PoldiSpectrumPawleyFunction"));
+
+  if (!poldiPawleyFunction) {
+    throw std::invalid_argument("Could not create pawley function.");
+  }
+
+  poldiPawleyFunction->setDecoratedFunction("PawleyFunction");
+
+  IPawleyFunction_sptr pawleyFunction =
+      poldiPawleyFunction->getPawleyFunction();
+  pawleyFunction->setProfileFunction(profileFunctionName);
+
+  std::string crystalSystem = getProperty("CrystalSystem");
+  pawleyFunction->setCrystalSystem(crystalSystem);
+
+  std::string initialCell = getProperty("InitialCell");
+  pawleyFunction->setUnitCell(initialCell);
+
+  IPeakFunction_sptr pFun = boost::dynamic_pointer_cast<IPeakFunction>(
+      FunctionFactory::Instance().createFunction(profileFunctionName));
+
+  for (size_t i = 0; i < peakCollection->peakCount(); ++i) {
+    PoldiPeak_sptr peak = peakCollection->peak(i);
+
+    pFun->setCentre(peak->d());
+    pFun->setFwhm(peak->fwhm(PoldiPeak::AbsoluteD));
+    pFun->setIntensity(peak->intensity());
+
+    pawleyFunction->addPeak(peak->hkl().asV3D(),
+                            peak->fwhm(PoldiPeak::AbsoluteD),
+                            pFun->height());
+
+    IPeakFunction_sptr p = pawleyFunction->getPeakFunction(i);
+    V3D h = pawleyFunction->getPeakHKL(i);
+    std::cout << p->centre() << " " << p->fwhm() << " " << h << std::endl;
+  }
+
+  pawleyFunction->fix(pawleyFunction->parameterIndex("f0.ZeroShift"));
+
+  mdFunction->addFunction(poldiPawleyFunction);
+
+  for (size_t i = 0; i < mdFunction->nParams(); ++i) {
+    std::cout << i << " " << mdFunction->parameterName(i) << " "
+              << mdFunction->getParameter(i) << std::endl;
+  }
+
+  return mdFunction;
+}
+
+/**
+ * Constructs a proper function from a peak collection
+ *
+ * This method constructs a Poldi2DFunction and assigns one
+ *PoldiSpectrumDomainFunction to it for
+ * each peak contained in the peak collection.
+ *
+ * @param peakCollection :: PoldiPeakCollection containing peaks with integral
+ *intensities
+ * @return Poldi2DFunction with one PoldiSpectrumDomainFunction per peak
+ */
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
+    const PoldiPeakCollection_sptr &peakCollection) const {
+  std::string profileFunctionName = getProperty("PeakProfileFunction");
+
+  bool pawleyFit = getProperty("PawleyFit");
+  if (pawleyFit) {
+    return getFunctionPawley(profileFunctionName, peakCollection);
+  }
+
+  return getFunctionIndividualPeaks(profileFunctionName, peakCollection);
+}
+
 /// Executes the algorithm
 void PoldiFitPeaks2D::exec() {
   TableWorkspace_sptr peakTable = getProperty("PoldiPeakWorkspace");