Re #4158. Optimized LeastSquares for use in MD fitting

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

@@ -42,7 +42,7 @@ namespace CurveFitting
 class DLLExport CostFuncFitting : public API::ICostFunction 
 {
 public:
-
+  CostFuncFitting();
   /// Get i-th parameter
   /// @param i :: Index of a parameter
   /// @return :: Value of the parameter
@@ -65,7 +65,7 @@ class DLLExport CostFuncFitting : public API::ICostFunction
   /// @param covar :: Returned covariance matrix, here as 
   /// @param epsrel :: Is used to remove linear-dependent columns
   ///
-  virtual void calCovarianceMatrix(GSLMatrix& covar, double epsrel = 0.0001);
+  virtual void calCovarianceMatrix(GSLMatrix& covar, double epsrel = 1e-8);
 
   /// Calculate fitting errors
   virtual void calFittingErrors(const GSLMatrix& covar);
@@ -74,14 +74,24 @@ class DLLExport CostFuncFitting : public API::ICostFunction
 
 protected:
 
+  /**
+   * Calculates covariance matrix for fitting function's active parameters. 
+   */
+  virtual void calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel = 1e-8);
+
   bool isValid() const;
   void checkValidity() const;
   void calTransformationMatrixNumerically(GSLMatrix& tm);
+  void setDirty();
 
   API::IFunction_sptr m_function;
   API::FunctionDomain_sptr m_domain;
   API::FunctionValues_sptr m_values;
   std::vector<size_t> m_indexMap;
+
+  mutable bool m_dirtyVal;
+  mutable bool m_dirtyDeriv;
+  mutable bool m_dirtyHessian;
 };
 
 } // namespace CurveFitting

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

@@ -40,14 +40,13 @@ namespace CurveFitting
 class DLLExport CostFuncLeastSquares : public CostFuncFitting
 {
 public:
+  /// Constructor
+  CostFuncLeastSquares():CostFuncFitting(),m_value(0),m_pushed(false){}
   /// Virtual destructor
   virtual ~CostFuncLeastSquares() {}
 
-  /// Constructor
-  CostFuncLeastSquares() : m_name("Least squares") { }
-
   /// Get name of minimizer
-  virtual std::string name() const { return m_name;}
+  virtual std::string name() const { return "Least squares";}
 
   /// Get short name of minimizer - useful for say labels in guis
   virtual std::string shortName() const {return "Chi-sq";};
@@ -64,14 +63,30 @@ class DLLExport CostFuncLeastSquares : public CostFuncFitting
   /// @return :: The value of the function
   virtual double valAndDeriv(std::vector<double>& der) const;
 
-  virtual double valDerivHessian(GSLVector& der, GSLMatrix& hessian, bool evalFunction = true) const;
+  virtual double valDerivHessian(bool evalFunction = true) const;
+  const GSLVector& getDeriv() const;
+  const GSLMatrix& getHessian() const;
+  void push();
+  void pop();
+  void drop();
 
   void setParameters(const GSLVector& params);
   void getParameters(GSLVector& params) const;
 
+protected:
+
+  virtual void calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel = 1e-8);
+
 private:
-  /// name of this minimizer
-  const std::string m_name;
+
+  mutable double m_value;
+  mutable GSLVector m_der;
+  mutable GSLMatrix m_hessian;
+
+  mutable bool m_pushed;
+  mutable double m_pushedValue;
+  mutable GSLVector m_pushedParams;
+
 };
 
 } // namespace CurveFitting

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

@@ -2,12 +2,16 @@
 #define MANTID_CURVEFITTING_GSLMATRIX_H_
 
 #include "MantidCurveFitting/DllConfig.h"
+#include "MantidCurveFitting/GSLVector.h"
 
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
 
 #include <vector>
 #include <stdexcept>
+#include <iostream>
+#include <iomanip>
 
 namespace Mantid
 {
@@ -250,6 +254,43 @@ namespace Mantid
 
       return *this;
     }
+
+    /// Solve system of linear equations M*x == rhs, M is this matrix
+    /// This matrix is destroyed.
+    void solve(const GSLVector& rhs, GSLVector& x)
+    {
+      if (size1() != size2())
+      {
+        throw std::runtime_error("System of linear equations: the matrix must be square.");
+      }
+      size_t n = size1();
+      if (rhs.size() != n)
+      {
+        throw std::runtime_error("System of linear equations: right-hand side vector has wrong size.");
+      }
+      x.resize(n);
+      int s;
+      gsl_permutation * p = gsl_permutation_alloc( n );
+      gsl_linalg_LU_decomp( gsl(), p, &s ); // matrix is modified at this moment
+      gsl_linalg_LU_solve( gsl(), p, rhs.gsl(), x.gsl() );
+      gsl_permutation_free( p );
+    }
+
+    /// Invert this matrix
+    void invert()
+    {
+      if (size1() != size2())
+      {
+        throw std::runtime_error("Matrix inverse: the matrix must be square.");
+      }
+      size_t n = size1();
+      int s;
+      GSLMatrix LU(*this);
+      gsl_permutation * p = gsl_permutation_alloc( n );
+      gsl_linalg_LU_decomp( LU.gsl(), p, &s );
+      gsl_linalg_LU_invert( LU.gsl(), p, gsl() );
+      gsl_permutation_free( p );
+    }
   };
 
   inline GSLMatrixMult2 operator*(const GSLMatrix& m1, const GSLMatrix& m2)
@@ -292,6 +333,19 @@ namespace Mantid
     return GSLMatrixMult3(mm,m);
   }
 
+  inline std::ostream& operator<<(std::ostream& ostr, const GSLMatrix& m)
+  {
+    ostr << std::scientific << std::setprecision(6);
+    for(size_t i = 0; i < m.size1(); ++i)
+    {
+      for(size_t j = 0; j < m.size2(); ++j)
+      {
+        ostr << std::setw(13) << m.get(i,j) << ' ';
+      }
+      ostr << std::endl;
+    }
+    return ostr;
+  }
 
   } // namespace CurveFitting
 } // namespace Mantid

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

@@ -71,11 +71,6 @@ class DLLExport LevenbergMarquardtMDMinimizer : public IFuncMinimizer
   double m_rho;
   /// To keep function value
   double m_F;
-  /// To keep first derivatives
-  GSLVector m_der;
-  /// To keep second derivatives
-  GSLMatrix m_hessian;
-  double m_oldDder;
   std::vector<double> m_D;
 	/// Static reference to the logger class
 	static Kernel::Logger& g_log;

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

@@ -12,6 +12,20 @@ namespace Mantid
 namespace CurveFitting
 {
 
+CostFuncFitting::CostFuncFitting():
+m_dirtyVal(true),
+m_dirtyDeriv(true),
+m_dirtyHessian(true)
+{
+}
+
+void CostFuncFitting::setDirty()
+{
+  m_dirtyVal = true;
+  m_dirtyDeriv = true;
+  m_dirtyHessian = true;
+}
+
 /// Get i-th parameter
 /// @param i :: Index of a parameter
 /// @return :: Value of the parameter
@@ -28,6 +42,7 @@ void CostFuncFitting::setParameter(size_t i, const double& value)
 {
   checkValidity();
   m_function->setActiveParameter(m_indexMap[i],value);
+  setDirty();
 }
 
 /// Number of parameters
@@ -87,24 +102,45 @@ void CostFuncFitting::checkValidity() const
   }
 }
 
-/** Calculates covariance matrix
-  *
-  * @param covar :: Returned covariance matrix
-  * @param epsrel :: Is used to remove linear-dependent columns
+/**
+  * Calculates covariance matrix for fitting function's active parameters. 
   */
-void CostFuncFitting::calCovarianceMatrix( GSLMatrix& covar, double epsrel )
+void CostFuncFitting::calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel)
 {
   // construct the jacobian
   GSLJacobian J( m_function, m_values->size() );
   size_t na = this->nParams(); // number of active parameters
   assert( J.getJ()->size2 == na );
+  covar.resize(na,na);
 
   // calculate the derivatives
   m_function->functionDeriv( *m_domain, J );
 
+  //std::cerr << "J=\n";
+  //for(size_t i = 0; i < J.getJ()->size1; ++i)
+  //{
+  //  for(size_t j = 0; j < J.getJ()->size2; ++j)
+  //  {
+  //    std::cerr << std::scientific << std::setprecision(6) << std::setw(13);
+  //    std::cerr << gsl_matrix_get(J.getJ(),i,j) << ' ';
+  //  }
+  //  std::cerr << std::endl;
+  //}
+
   // let the GSL to compute the covariance matrix
-  GSLMatrix c( na, na );
-  gsl_multifit_covar( J.getJ(), epsrel, c.gsl() );
+  gsl_multifit_covar( J.getJ(), epsrel, covar.gsl() );
+
+}
+
+/** Calculates covariance matrix
+  *
+  * @param covar :: Returned covariance matrix
+  * @param epsrel :: Is used to remove linear-dependent columns
+  */
+void CostFuncFitting::calCovarianceMatrix( GSLMatrix& covar, double epsrel )
+{
+  GSLMatrix c;
+  calActiveCovarianceMatrix( c, epsrel );
 
   size_t np = m_function->nParams();
 
@@ -118,6 +154,7 @@ void CostFuncFitting::calCovarianceMatrix( GSLMatrix& covar, double epsrel )
     ++ii;
   }
 
+  //std::cerr << "c=\n" << c << std::endl;
   if (isTransformationIdentity)
   {
     // if the transformation is identity simply copy the matrix
@@ -129,8 +166,11 @@ void CostFuncFitting::calCovarianceMatrix( GSLMatrix& covar, double epsrel )
     GSLMatrix tm;
     calTransformationMatrixNumerically(tm);
     covar = Tr(tm) * c * tm;
+    //std::cerr << "tm:\n" << tm << std::endl;
   }
 
+  //std::cerr << "Covar:\n" << covar << std::endl;
+
 }
 
 /**