Mantid
Loading...
Searching...
No Matches
DampedGaussNewtonMinimizer.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
13
16#include "MantidAPI/IFunction.h"
17
18#include "MantidKernel/Logger.h"
19
20#include <cmath>
21#include <gsl/gsl_blas.h>
22
24
25namespace {
27Kernel::Logger g_log("DampedGaussNewtonMinimizer");
28} // namespace
29
30DECLARE_FUNCMINIMIZER(DampedGaussNewtonMinimizer, Damped GaussNewton)
31
32
34 declareProperty("Damping", 0.0, "The damping parameter.");
35 declareProperty("Verbose", false, "Make output more verbose.");
36}
37
39void DampedGaussNewtonMinimizer::initialize(API::ICostFunction_sptr function, size_t /*maxIterations*/) {
40 m_leastSquares = std::dynamic_pointer_cast<CostFunctions::CostFuncLeastSquares>(function);
41 if (!m_leastSquares) {
42 throw std::invalid_argument("Damped Gauss-Newton minimizer works only with least "
43 "squares. Different function was given.");
44 }
45}
46
48bool DampedGaussNewtonMinimizer::iterate(size_t /*iteration*/) {
49 const bool verbose = getProperty("Verbose");
50 const double damping = getProperty("Damping");
51
52 if (!m_leastSquares) {
53 throw std::runtime_error("Cost function isn't set up.");
54 }
55 size_t n = m_leastSquares->nParams();
56
57 if (n == 0) {
58 m_errorString = "No parameters to fit";
59 return false;
60 }
61
62 // calculate the first and second derivatives of the cost function.
63 m_leastSquares->valDerivHessian();
64
65 // copy the hessian
66 EigenMatrix H(m_leastSquares->getHessian());
67 EigenVector dd(m_leastSquares->getDeriv());
68
69 for (size_t i = 0; i < n; ++i) {
70 double tmp = H.get(i, i) + damping;
71 if (tmp == 0.0) {
72 m_errorString = "Function doesn't depend on parameter " + m_leastSquares->parameterName(i);
73 return false;
74 }
75 H.set(i, i, tmp);
76 }
77
78 if (verbose) {
79 g_log.warning() << "H:\n" << H;
80 g_log.warning() << "-----------------------------\n";
81 for (size_t j = 0; j < n; ++j) {
82 g_log.warning() << dd.get(j) << ' ';
83 }
84 g_log.warning() << '\n';
85 }
86
88 EigenVector dx(n);
89 // To find dx solve the system of linear equations H * dx == -m_der
90 dd *= -1.0;
91 try {
92 H.solve(dd, dx);
93 } catch (std::runtime_error &e) {
94 m_errorString = e.what();
95 return false;
96 }
97
98 if (verbose) {
99 for (size_t j = 0; j < n; ++j) {
100 g_log.warning() << dx.get(j) << ' ';
101 }
102 g_log.warning() << "\n\n";
103 }
104
105 // Update the parameters of the cost function.
106 for (size_t i = 0; i < n; ++i) {
107 if (!std::isfinite(dx[i])) {
108 m_errorString = "Encountered an infinite number or NaN.";
109 return false;
110 }
111 double d = m_leastSquares->getParameter(i) + dx.get(i);
112 m_leastSquares->setParameter(i, d);
113 if (verbose) {
114 g_log.warning() << i << " Parameter " << m_leastSquares->parameterName(i) << ' ' << d << '\n';
115 }
116 }
117 m_leastSquares->getFittingFunction()->applyTies();
118
119 // --- prepare for the next iteration --- //
120
121 // Try the stop condition
122 EigenVector p(n);
123 m_leastSquares->getParameters(p);
124 gsl_vector_view dx_gsl = getGSLVectorView(dx.mutator());
125 double dx_norm = gsl_blas_dnrm2(&dx_gsl.vector);
126 return dx_norm >= m_relTol;
127}
128
131 if (!m_leastSquares) {
132 throw std::runtime_error("Cost function isn't set up.");
133 }
134 return m_leastSquares->val();
135}
136
137} // namespace Mantid::CurveFitting::FuncMinimisers
gsl_vector * tmp
#define DECLARE_FUNCMINIMIZER(classname, username)
Macro for declaring a new type of minimizers to be used with the FuncMinimizerFactory.
An interface for function minimizers.
std::string m_errorString
Error string.
A wrapper around Eigen::Matrix.
Definition: EigenMatrix.h:33
A wrapper around Eigen::Vector.
Definition: EigenVector.h:27
double get(const size_t i) const
Get an element.
vec_map_type & mutator()
Get the map of the eigen vector.
Definition: EigenVector.h:51
Implements a Gauss-Newton minimization algorithm with damping for use with least squares cost functio...
std::shared_ptr< CostFunctions::CostFuncLeastSquares > m_leastSquares
Pointer to the cost function. Must be the least squares.
void initialize(API::ICostFunction_sptr function, size_t maxIterations=0) override
Initialize minimizer, i.e. pass a function to minimize.
double costFunctionVal() override
Return current value of the cost function.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< ICostFunction > ICostFunction_sptr
define a shared pointer to a cost function
Definition: ICostFunction.h:60
gsl_vector_view getGSLVectorView(vec_map_type &v)
take data from Eigen Vector and take a gsl view
Definition: GSLFunctions.h:54