Refs #9782. Added gaussian with hand-coded derivatives

Code/Mantid/Framework/CurveFitting/CMakeLists.txt

@@ -48,6 +48,7 @@ set ( SRC_FILES
 	src/Gaussian.cpp
 	src/GaussianAutoDiff.cpp
 	src/GaussianComptonProfile.cpp
+	src/GaussianHandCoded.cpp
 	src/GaussianNumDiff.cpp
 	src/GramCharlierComptonProfile.cpp
 	src/IkedaCarpenterPV.cpp
@@ -155,6 +156,7 @@ set ( INC_FILES
 	inc/MantidCurveFitting/Gaussian.h
 	inc/MantidCurveFitting/GaussianAutoDiff.h
 	inc/MantidCurveFitting/GaussianComptonProfile.h
+	inc/MantidCurveFitting/GaussianHandCoded.h
 	inc/MantidCurveFitting/GaussianNumDiff.h
 	inc/MantidCurveFitting/GramCharlierComptonProfile.h
 	inc/MantidCurveFitting/IkedaCarpenterPV.h
@@ -251,6 +253,7 @@ set ( TEST_FILES
 	GausOscTest.h
 	GaussianAutoDiffTest.h
 	GaussianComptonProfileTest.h
+	GaussianHandCodedTest.h
 	GaussianNumDiffTest.h
 	GaussianTest.h
 	GramCharlierComptonProfileTest.h

Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/GaussianHandCoded.h

@@ -0,0 +1,55 @@
+#ifndef MANTID_CURVEFITTING_GAUSSIANHANDCODED_H_
+#define MANTID_CURVEFITTING_GAUSSIANHANDCODED_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/IFunction1D.h"
+#include "MantidAPI/ParamFunction.h"
+
+
+namespace Mantid
+{
+namespace CurveFitting
+{
+
+  /** GaussianHandCoded : TODO: DESCRIPTION
+    
+    Copyright © 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+    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 DLLExport GaussianHandCoded : public API::IFunction1D, public API::ParamFunction
+  {
+  public:
+    GaussianHandCoded();
+    virtual ~GaussianHandCoded();
+    std::string name() const { return "GaussianHandCoded"; }
+
+    void function1D(double *out, const double *xValues, const size_t nData) const;
+    void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData);
+
+  protected:
+    void init();
+
+};
+
+
+} // namespace CurveFitting
+} // namespace Mantid
+
+#endif  /* MANTID_CURVEFITTING_GAUSSIANHANDCODED_H_ */

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

@@ -1,7 +1,10 @@
 #include "MantidCurveFitting/AutoDiffTestAlg.h"
 #include "MantidCurveFitting/GaussianNumDiff.h"
-
 #include "MantidCurveFitting/GaussianAutoDiff.h"
+#include "MantidCurveFitting/GaussianHandCoded.h"
+#include "MantidAPI/FunctionDomain1D.h"
+#include "MantidAPI/FunctionValues.h"
+#include "MantidCurveFitting/Jacobian.h"
 #include "MantidAPI/CompositeFunction.h"
 #include "MantidKernel/Timer.h"
 #include "MantidKernel/Logger.h"
@@ -53,8 +56,7 @@ const std::string AutoDiffTestAlg::summary() const { return "Test"; }
 void AutoDiffTestAlg::init()
 {
     declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input), "Data with 20 gaussian peaks.");
-    declareProperty(new WorkspaceProperty<ITableWorkspace>("OutputWorkspace","",Direction::Output), "Data with 20 gaussian peaks.");
-    declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspacePlot","",Direction::Output), "Data with 20 gaussian peaks.");
+    declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "Data with 20 gaussian peaks.");
     declareProperty("DerivativeType","adept","How to calculate derivatives.");
 
 }
@@ -71,7 +73,7 @@ namespace
     if ( type == "adept" )
     {
       return IFunction_sptr(new GaussianAutoDiff);
-    }
+    } else if ( type == "num") {
     /*
     API::CompositeFunction_sptr comp(new API::CompositeFunction);
     auto gauss = IFunction_sptr(new GaussianNumDiff);
@@ -81,6 +83,9 @@ namespace
     return comp;
     */
     return IFunction_sptr(new GaussianNumDiff);
+    }
+
+    return IFunction_sptr(new GaussianHandCoded);
   }
 }
 
@@ -160,6 +165,7 @@ void AutoDiffTestAlg::exec()
     sigma.push_back(5.657079363007425);
 
     // make composite function
+    /*
     CompositeFunction_sptr bigFunction(new CompositeFunction);
     for(size_t i = 0; i < 20; ++i) {
         IFunction_sptr g = getFunction(m_derType);
@@ -172,26 +178,76 @@ void AutoDiffTestAlg::exec()
 
         bigFunction->addFunction(g);
     }
+    */
 
     Kernel::Timer timer;
 
-    IAlgorithm_sptr fitAlgorithm = createChildAlgorithm("Fit", -1, -1, true);
-    fitAlgorithm->setProperty("CreateOutput", true);
-    fitAlgorithm->setProperty("Output", "FitPeaks1D");
-    fitAlgorithm->setProperty("CalcErrors", true);
-    fitAlgorithm->setProperty("Function", boost::static_pointer_cast<IFunction>(bigFunction));
-    fitAlgorithm->setProperty("InputWorkspace", fitData);
-    fitAlgorithm->setProperty("WorkspaceIndex", 0);
-
-    fitAlgorithm->execute();
+    for(size_t i = 0; i < 10; ++i) {
+        IFunction_sptr g = getFunction(m_derType);
+        g->initialize();
+        g->setParameter(0, 4.0); // h
+        g->setParameter(1, 0.0); // x0
+        g->setParameter(2, 3.5); // s
+
+        IAlgorithm_sptr fitAlgorithm = createChildAlgorithm("Fit", -1, -1, false);
+        fitAlgorithm->setProperty("CreateOutput", true);
+        fitAlgorithm->setProperty("Output", "FitPeaks1D");
+        fitAlgorithm->setProperty("CalcErrors", true);
+        fitAlgorithm->setProperty("Function", g);
+        fitAlgorithm->setProperty("InputWorkspace", fitData);
+        fitAlgorithm->setProperty("WorkspaceIndex", 0);
+
+        fitAlgorithm->execute();
+    }
 
     g_log.warning() << "Fit took " << timer.elapsed() << " seconds to complete" << std::endl;
 
-    ITableWorkspace_sptr t = fitAlgorithm->getProperty("OutputParameters");
-    MatrixWorkspace_sptr mw = fitAlgorithm->getProperty("OutputWorkspace");
+    IFunction_sptr g = getFunction(m_derType);
+    g->initialize();
+    g->setParameter(0, 4.0); // h
+    g->setParameter(1, 0.0); // x0
+    g->setParameter(2, 3.5); // s
+
+    FunctionDomain1DVector x(fitData->readX(0));
+    FunctionValues y(x);
+    CurveFitting::Jacobian J(500, 3);
+
+    g->function(x, y);
+    g->functionDeriv(x, J);
+
+    MatrixWorkspace_sptr t = WorkspaceFactory::Instance().create(fitData);
+
+    MantidVec &yd = t->dataY(0);
+    const MantidVec &yr = fitData->readY(0);
+    for(size_t i = 0; i < 500; ++i) {
+        yd[i] = yr[i] - y[i];
+    }
+
+    MantidVec &dfdx0 = t->dataY(1);
+    const MantidVec &dfdx0r = fitData->readY(1);
+    for(size_t i = 0; i < 500; ++i) {
+        dfdx0[i] = dfdx0r[i] - J.get(i, 1);
+    }
+
+    MantidVec &dfdh = t->dataY(2);
+    const MantidVec &dfdhr = fitData->readY(2);
+    for(size_t i = 0; i < 500; ++i) {
+        dfdh[i] = dfdhr[i] - J.get(i, 0);
+    }
+
+    MantidVec &dfds = t->dataY(3);
+    const MantidVec &dfdsr = fitData->readY(3);
+    for(size_t i = 0; i < 500; ++i) {
+        dfds[i] = dfdsr[i] - J.get(i, 2);
+    }
+
+    for(size_t i = 0; i < 4; ++i) {
+        MantidVec &xd = t->dataX(i);
+        xd = fitData->readX(i);
+    }
 
+    // test accuracy
     setProperty("OutputWorkspace", t);
-    setProperty("OutputWorkspacePlot", mw);
 }
 
 

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

@@ -0,0 +1,59 @@
+#include "MantidCurveFitting/GaussianHandCoded.h"
+
+namespace Mantid
+{
+namespace CurveFitting
+{
+
+using namespace API;
+  //----------------------------------------------------------------------------------------------
+  /** Constructor
+   */
+  GaussianHandCoded::GaussianHandCoded()
+  {
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Destructor
+   */
+  GaussianHandCoded::~GaussianHandCoded()
+  {
+  }
+
+  void GaussianHandCoded::function1D(double *out, const double *xValues, const size_t nData) const
+  {
+      double height = getParameter("Height");
+      double centre = getParameter("PeakCentre");
+      double sigma = getParameter("Sigma");
+
+      for(size_t i = 0; i < nData; ++i) {
+          out[i] = height * exp(-0.5 * pow( (xValues[i] - centre)/sigma, 2.0) );
+      }
+  }
+
+  void GaussianHandCoded::functionDeriv1D(Jacobian *out, const double *xValues, const size_t nData)
+  {
+      double height = getParameter("Height");
+      double centre = getParameter("PeakCentre");
+      double sigma = getParameter("Sigma");
+
+      for(size_t i = 0; i < nData; ++i) {
+          double xDiff = xValues[i] - centre;
+          double expTerm = exp(-0.5 * pow( (xDiff)/sigma, 2.0));
+
+          out->set(i, 0, expTerm);
+          out->set(i, 1, (height*expTerm*(2 * xDiff))/(2*pow(sigma, 2)));
+          out->set(i, 2, (height*expTerm*pow(xDiff, 2))/pow(sigma, 3));
+      }
+  }
+
+  void GaussianHandCoded::init() {
+      declareParameter("Height");
+      declareParameter("PeakCentre");
+      declareParameter("Sigma");
+  }
+
+
+
+} // namespace CurveFitting
+} // namespace Mantid

Code/Mantid/Framework/CurveFitting/test/GaussianHandCodedTest.h

@@ -0,0 +1,29 @@
+#ifndef MANTID_CURVEFITTING_GAUSSIANHANDCODEDTEST_H_
+#define MANTID_CURVEFITTING_GAUSSIANHANDCODEDTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidCurveFitting/GaussianHandCoded.h"
+
+using Mantid::CurveFitting::GaussianHandCoded;
+using namespace Mantid::API;
+
+class GaussianHandCodedTest : 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 GaussianHandCodedTest *createSuite() { return new GaussianHandCodedTest(); }
+  static void destroySuite( GaussianHandCodedTest *suite ) { delete suite; }
+
+
+  void test_Something()
+  {
+    TSM_ASSERT( "You forgot to write a test!", 0);
+  }
+
+
+};
+
+
+#endif /* MANTID_CURVEFITTING_GAUSSIANHANDCODEDTEST_H_ */