GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GOptimizerLM.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GOptimizerLM.cpp - Levenberg Marquardt optimizer *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GOptimizerLM.cpp
23  * @brief Levenberg-Marquardt optimizer class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GOptimizerLM.hpp"
32 #include "GTools.hpp"
33 #include "GException.hpp"
34 
35 /* __ Method name definitions ____________________________________________ */
36 
37 /* __ Macros _____________________________________________________________ */
38 
39 /* __ Coding definitions _________________________________________________ */
40 
41 /* __ Debug definitions __________________________________________________ */
42 //#define G_DEBUG_OPT //!< Define to debug optimize() method
43 //#define G_DEBUG_ITER //!< Define to debug iteration() method
44 //#define G_DEBUG_SHOW_GRAD_COVAR //!< Define to show grad and curvature
45 
46 
47 /*==========================================================================
48  = =
49  = Constructors/destructors =
50  = =
51  ==========================================================================*/
52 
53 /***********************************************************************//**
54  * @brief Void constructor
55  ***************************************************************************/
57 {
58  // Initialise private members for clean destruction
59  init_members();
60 
61  // Return
62  return;
63 }
64 
65 
66 /***********************************************************************//**
67  * @brief Constructor with logger
68  *
69  * @param[in] log Logger to use in optimizer.
70  ***************************************************************************/
72 {
73  // Initialise private members for clean destruction
74  init_members();
75 
76  // Set pointer to logger
77  m_logger = log;
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Copy constructor
86  *
87  * @param[in] opt Optimizer from which the instance should be built.
88  ***************************************************************************/
90 {
91  // Initialise private members for clean destruction
92  init_members();
93 
94  // Copy members
95  copy_members(opt);
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Destructor
104  ***************************************************************************/
106 {
107  // Free members
108  free_members();
109 
110  // Return
111  return;
112 }
113 
114 
115 /*==========================================================================
116  = =
117  = Operators =
118  = =
119  ==========================================================================*/
120 
121 /***********************************************************************//**
122  * @brief Assignment operator
123  *
124  * @param[in] opt Optimizer to be assigned.
125  ***************************************************************************/
127 {
128  // Execute only if object is not identical
129  if (this != &opt) {
130 
131  // Copy base class members
132  this->GOptimizer::operator=(opt);
133 
134  // Free members
135  free_members();
136 
137  // Initialise private members for clean destruction
138  init_members();
139 
140  // Copy members
141  copy_members(opt);
142 
143  } // endif: object was not identical
144 
145  // Return
146  return *this;
147 }
148 
149 
150 /*==========================================================================
151  = =
152  = Public methods =
153  = =
154  ==========================================================================*/
155 
156 /***********************************************************************//**
157  * @brief Clear object
158  *
159  * This method properly resets the object to an initial state.
160  ***************************************************************************/
162 {
163  // Free class members (base and derived classes, derived class first)
164  free_members();
165  this->GOptimizer::free_members();
166 
167  // Initialise members
168  this->GOptimizer::init_members();
169  init_members();
170 
171  // Return
172  return;
173 }
174 
175 
176 /***********************************************************************//**
177  * @brief Clone object
178 ***************************************************************************/
180 {
181  return new GOptimizerLM(*this);
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief Optimize function parameters
187  *
188  * @param[in] fct Optimization function.
189  * @param[in] pars Function parameters.
190  ***************************************************************************/
192 {
193  // Initialize optimizer parameters
195  m_value = 0.0;
196  m_delta = 0.0;
198  m_iter = 0;
199  m_num_dec = 0;
200 
201  // Get number of parameters. Continue only if there are free parameters
202  m_npars = pars.size();
203  m_nfree = pars.nfree();
204  if (m_nfree > 0) {
205 
206  // Initialise bookkeeping arrays
207  m_hit_boundary.clear();
208  m_hit_minimum.clear();
209  m_hit_maximum.clear();
210  m_par_freeze.clear();
211  m_par_remove.clear();
212  m_hit_boundary.reserve(m_npars);
213  m_hit_minimum.reserve(m_npars);
214  m_hit_maximum.reserve(m_npars);
215  m_par_freeze.reserve(m_npars);
216  m_par_remove.reserve(m_npars);
217  for (int i = 0; i < m_npars; ++i) {
218  m_hit_boundary.push_back(false);
219  m_hit_minimum.push_back(0);
220  m_hit_maximum.push_back(0);
221  m_par_freeze.push_back(false);
222  m_par_remove.push_back(false);
223  }
224 
225  // Initial function evaluation
226  fct.eval(pars);
227 
228  // If a free parameter has a zero diagonal element in the curvature
229  // matrix then remove this parameter definitely from the fit as it
230  // otherwise will block the fit. The problem appears in the unbinned
231  // fitting where parameter gradients may be zero (due to the truncation
232  // of the PSF), but the Npred gradient is not zero. In principle we
233  // could use the Npred gradient for fitting (I guess), but I still
234  // have to figure out how ... (the diagonal loading was not so
235  // successful as it faked early convergence)
236  for (int ipar = 0; ipar < m_npars; ++ipar) {
237  if (pars[ipar]->is_free()) {
238  if ((*fct.curvature())(ipar,ipar) == 0.0) {
239  if (m_logger != NULL) {
240  *m_logger << " Parameter \"" << pars[ipar]->name();
241  *m_logger << "\" has zero curvature.";
242  *m_logger << " Fix parameter." << std::endl;
243  }
244  m_par_remove[ipar] = true;
245  pars[ipar]->fix();
246  }
247  }
248  }
249 
250  // Save function value
251  m_value = fct.value();
252 
253  // Save initial statistics and lambda values
254  double lambda_old = m_lambda;
255  int lambda_inc = 0;
256 
257  // Optionally write initial iteration into logger
258  if (m_logger != NULL) {
259  (*m_logger)(">Iteration %3d: -logL=%.3f, Lambda=%.1e",
260  0, m_value, m_lambda);
261  }
262  #if defined(G_DEBUG_OPT)
263  std::cout << "Initial iteration: func=" << m_value << ", Lambda="
264  << m_lambda << std::endl;
265  #endif
266  #if defined(G_DEBUG_SHOW_GRAD_COVAR)
267  if (m_logger != NULL) {
268  *m_logger << *fct.gradient() << std::endl;
269  *m_logger << *fct.curvature() << std::endl;
270  }
271  #endif
272 
273  // Initialise iteration flags
274  bool check_for_freeze = true;
275 
276  // Iterative fitting
277  for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) {
278 
279  // Perform one iteration
280  double step = iteration(fct, pars);
281 
282  // Determine maximum (scaled) gradient
283  double grad_max = 0.0;
284  int grad_imax = -1;
285  for (int ipar = 0; ipar < m_npars; ++ipar) {
286  if (pars[ipar]->is_free()) {
287  double grad = pars[ipar]->factor_gradient();
288  if (std::abs(grad) > std::abs(grad_max)) {
289  grad_max = grad;
290  grad_imax = ipar;
291  }
292  if (grad == 0.0) {
293  if (m_logger != NULL) {
294  *m_logger << "Parameter " << ipar;
295  *m_logger << " (" << pars[ipar]->name() << ")";
296  *m_logger << " has a zero gradient." << std::endl;
297  }
298  }
299  }
300  }
301 
302  // Optionally write iteration results into logger
303  if (m_logger != NULL) {
304  std::string stalled = "";
305  std::string status = "";
306  if (m_lambda > lambda_old) {
307  status = " ";
308  stalled = " (stalled)";
309  }
310  else {
311  status = ">";
312  }
313  std::string parname = "";
314  if (grad_imax != -1) {
315  parname = " [" + pars[grad_imax]->name() + ":" +
316  gammalib::str(grad_imax) + "]";
317  }
318  (*m_logger)("%sIteration %3d: -logL=%.3f, Lambda=%.1e,"
319  " delta=%.3f, step=%.1e, max(|grad|)=%f%s%s",
320  status.c_str(), m_iter, m_value, lambda_old,
321  m_delta, step, grad_max,
322  parname.c_str(), stalled.c_str());
323  }
324  #if defined(G_DEBUG_OPT)
325  std::cout << "Iteration " << m_iter << ": func="
326  << m_value << ", Lambda=" << lambda_old
327  << ", delta=" << m_delta << std::endl;
328  #endif
329  #if defined(G_DEBUG_SHOW_GRAD_COVAR)
330  if (m_logger != NULL) {
331  *m_logger << *fct.gradient() << std::endl;
332  *m_logger << *fct.curvature() << std::endl;
333  }
334  #endif
335 
336  // Reset lambda increment if we had success
337  if (m_lambda < lambda_old) {
338  lambda_inc = 0;
339  }
340 
341  // Stop if convergence was reached. Before stopping, check
342  // if some parameters were frozen, and if this was the case,
343  // free them now and continue. We do this only once, i.e.
344  // the next time a parameter is frozen and convergence is
345  // reached we really stop.
346  if ((m_lambda <= lambda_old) &&
347  (m_delta >= 0.0) &&
348  (m_delta < step*m_eps)) {
349 
350  // Check for frozen parameters, and if some exist, free
351  // them and start over
352  if (check_for_freeze) {
353 
354  // Signal that we won't check frozen parameters
355  // anymore
356  check_for_freeze = false;
357 
358  // Free frozen parameters and determine how many
359  // have been frozen
360  int nfrozen = 0;
361  for (int ipar = 0; ipar < m_npars; ++ipar) {
362  if (m_par_freeze[ipar]) {
363  nfrozen++;
364  pars[ipar]->free();
365  if (m_logger != NULL) {
366  *m_logger << " Free parameter \""
367  << pars[ipar]->name()
368  << "\" after convergence was"
369  " reached with frozen"
370  " parameter." << std::endl;
371  }
372  }
373  }
374 
375  // If there were frozen parameters then start over
376  // again (initialise optimizer)
377  if (nfrozen > 0) {
379  lambda_old = m_lambda;
380  lambda_inc = 0;
381  continue;
382  }
383  }
384 
385  // ... otherwise, convergence was reached and we can stop
386  // now
387  break;
388 
389  } // endif: convergence check
390 
391  // Monitor the number of subsequent increases of lambda and
392  // stop if the number of increases exceeds threshold
393  lambda_inc = (m_lambda > lambda_old) ? lambda_inc + 1 : 0;
394  if (lambda_inc >= m_max_stall) {
396  break;
397  }
398 
399  // Bookkeeping of actual result (we always store the last
400  // lambda to detect turn arounds in the lambda tendency; however
401  // we always keep the best function value)
402  lambda_old = m_lambda;
403 
404  } // endfor: iterations
405 
406  // Free now all temporarily frozen parameters so that the resulting
407  // model has the same attributes as the initial model
408  for (int ipar = 0; ipar < m_npars; ++ipar) {
409  if (m_par_freeze[ipar] || m_par_remove[ipar]) {
410  pars[ipar]->free();
411  if (m_logger != NULL) {
412  *m_logger << " Free parameter \""
413  << pars[ipar]->name()
414  << "\" after convergence was"
415  " reached with frozen"
416  " parameter." << std::endl;
417  }
418  }
419  }
420 
421  } // endif: there were free parameters to fit
422 
423  // ... otherwise just evaluate function and store the parameters
424  else {
425 
426  // Evaluate function
427  fct.eval(pars);
428 
429  // Save function value
430  m_value = fct.value();
431 
432  // Optionally write initial iteration into logger
433  if (m_logger != NULL) {
434  (*m_logger)(">Iteration %3d: -logL=%.3f, Lambda=%.1e (no free parameters)",
435  0, m_value, m_lambda);
436  }
437 
438  }
439 
440  // Return
441  return;
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Compute parameter uncertainties
447  *
448  * @param[in] fct Optimizer function.
449  * @param[in] pars Function parameters.
450  *
451  * Compute parameter uncertainties from the diagonal elements of the
452  * curvature matrix.
453  ***************************************************************************/
455 {
456  // Get number of parameters
457  int npars = pars.size();
458 
459  // Perform final parameter evaluation
460  fct.eval(pars);
461 
462  // Fetch sparse matrix pointer. We have to do this after the eval()
463  // method since eval() will allocate new memory for the curvature
464  // matrix!
465  GMatrixSparse* curvature = fct.curvature();
466 
467  // Save best fitting value
468  m_value = fct.value();
469 
470  // Save curvature matrix
471  GMatrixSparse save_curvature = GMatrixSparse(*curvature);
472 
473  // Signal no diagonal element loading
474  bool diag_loaded = false;
475 
476  // Loop over error computation (maximum 2 turns)
477  for (int i = 0; i < 2; ++i) {
478 
479  // Solve: curvature * X = unit
480  try {
481  GMatrixSparse decomposition = curvature->cholesky_decompose(true);
482  GVector unit(npars);
483  for (int ipar = 0; ipar < npars; ++ipar) {
484  unit[ipar] = 1.0;
485  GVector x = decomposition.cholesky_solver(unit, true);
486  if (x[ipar] >= 0.0) {
487  pars[ipar]->factor_error(sqrt(x[ipar]));
488  }
489  else {
490  pars[ipar]->factor_error(0.0);
492  }
493  unit[ipar] = 0.0;
494  }
495  }
496  catch (GException::matrix_zero &e) {
498  if (m_logger != NULL) {
499  *m_logger << "GOptimizerLM::terminate: "
500  << "All curvature matrix elements are zero."
501  << std::endl;
502  }
503  break;
504  }
506 
507  // Load diagonal if this has not yet been tried
508  if (!diag_loaded) {
509 
510  // Flag errors as inaccurate
512  if (m_logger != NULL) {
513  *m_logger << "Non-Positive definite curvature matrix encountered."
514  << std::endl;
515  *m_logger << "Load diagonal elements with 1e-10."
516  << " Fit errors may be inaccurate."
517  << std::endl;
518  }
519 
520  // Try now with diagonal loaded matrix
521  *curvature = save_curvature;
522  for (int ipar = 0; ipar < npars; ++ipar) {
523  (*curvature)(ipar,ipar) += 1.0e-10;
524  }
525 
526  // Signal loading
527  diag_loaded = true;
528 
529  // Try again
530  continue;
531 
532  } // endif: diagonal has not yet been loaded
533 
534  // ... otherwise signal an error
535  else {
537  if (m_logger != NULL) {
538  *m_logger << "Non-Positive definite curvature matrix encountered,"
539  << " even after diagonal loading." << std::endl;
540  }
541  break;
542  }
543  }
544  catch (std::exception &e) {
545  throw;
546  }
547 
548  // If no error occured then break now
549  break;
550 
551  } // endfor: looped over error computation
552 
553  // Return
554  return;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Print optimizer information
560  *
561  * @param[in] chatter Chattiness.
562  * @return String containing optimizer information.
563  ***************************************************************************/
564 std::string GOptimizerLM::print(const GChatter& chatter) const
565 {
566  // Initialise result string
567  std::string result;
568 
569  // Continue only if chatter is not silent
570  if (chatter != SILENT) {
571 
572  // Append header
573  result.append("=== GOptimizerLM ===");
574 
575  // Append information
576  result.append("\n"+gammalib::parformat("Optimized function value"));
577  result.append(gammalib::str(m_value, 3));
578  result.append("\n"+gammalib::parformat("Absolute precision"));
579  result.append(gammalib::str(m_eps));
580  result.append("\n"+gammalib::parformat("Acceptable value decrease"));
581  result.append(gammalib::str(std::abs(m_accept_dec)));
582 
583  // Append status
584  result.append("\n"+gammalib::parformat("Optimization status"));
585  switch (m_status) {
586  case G_LM_CONVERGED:
587  result.append("converged");
588  break;
589  case G_LM_STALLED:
590  result.append("stalled");
591  break;
592  case G_LM_SINGULAR:
593  result.append("singular curvature matrix encountered");
594  break;
596  result.append("curvature matrix not positive definite");
597  break;
598  case G_LM_BAD_ERRORS:
599  result.append("errors are inaccurate");
600  break;
601  default:
602  result.append("unknown");
603  break;
604  }
605 
606  // Append further information
607  result.append("\n"+gammalib::parformat("Number of parameters"));
608  result.append(gammalib::str(m_npars));
609  result.append("\n"+gammalib::parformat("Number of free parameters"));
610  result.append(gammalib::str(m_nfree));
611  result.append("\n"+gammalib::parformat("Number of iterations"));
612  result.append(gammalib::str(m_iter));
613  result.append("\n"+gammalib::parformat("Lambda"));
614  result.append(gammalib::str(m_lambda));
615 
616  } // endif: chatter was not silent
617 
618  // Return result
619  return result;
620 }
621 
622 
623 /*==========================================================================
624  = =
625  = Private methods =
626  = =
627  ==========================================================================*/
628 
629 /***********************************************************************//**
630  * @brief Initialise class members
631  ***************************************************************************/
633 {
634  // Initialise optimizer parameters
635  m_npars = 0;
636  m_nfree = 0;
637  m_lambda_start = 1.0e-3;
638  m_lambda_inc = 10.0;
639  m_lambda_dec = 0.1;
640  m_eps = 5.0e-3; //!< Changed on 30/10/2014 from 1.0e-6
641  m_accept_dec = 2.0; //!< Allow to decrease by 2.0
642  m_max_iter = 100; //!< Changed on 30/10/2014 from 1000
643  m_max_stall = 10;
644  m_max_hit = 3; //!< Maximum successive boundary hits before freeze
645  m_max_dec = 3; //!< Maximum number of function decreases
646  m_step_adjust = true;
647 
648  // Initialise bookkeeping arrays
649  m_hit_boundary.clear();
650  m_hit_minimum.clear();
651  m_hit_maximum.clear();
652  m_par_freeze.clear();
653  m_par_remove.clear();
654 
655  // Initialise optimizer values
657  m_value = 0.0;
658  m_delta = 0.0;
659  m_status = 0;
660  m_iter = 0;
661  m_num_dec = 0;
662 
663  // Initialise pointer to logger
664  m_logger = NULL;
665 
666  // Return
667  return;
668 }
669 
670 
671 /***********************************************************************//**
672  * @brief Copy class members
673  *
674  * @param[in] opt GOptimizerLM members to be copied.
675  ***************************************************************************/
677 {
678  // Copy attributes
682  m_eps = opt.m_eps;
684  m_max_iter = opt.m_max_iter;
685  m_max_stall = opt.m_max_stall;
686  m_max_dec = opt.m_max_dec;
693  m_lambda = opt.m_lambda;
694  m_value = opt.m_value;
695  m_delta = opt.m_delta;
696  m_status = opt.m_status;
697  m_iter = opt.m_iter;
698  m_logger = opt.m_logger;
699 
700  // Return
701  return;
702 }
703 
704 
705 /***********************************************************************//**
706  * @brief Delete class members
707  ***************************************************************************/
709 {
710  // Return
711  return;
712 }
713 
714 
715 /***********************************************************************//**
716  * @brief Perform one LM iteration
717  *
718  * @param[in] fct Optimizer function.
719  * @param[in] pars Function parameters.
720  * @return Step size for iteration.
721  *
722  * This method performs one LM iteration. Note that the method only acts on
723  * the parameter value, i.e. it does not worry about the true scale of the
724  * parameter. It calls the eval() method of the optimizer function which
725  * is assumed to return the gradient and curvature matrix with respect to the
726  * parameter values (and not the scaled true values).
727  ***************************************************************************/
729 {
730  // Debug option: dump header
731  #if defined(G_DEBUG_ITER)
732  std::cout << "GOptimizerLM::iteration: enter" << std::endl;
733  #endif
734 
735  // Initialise step size
736  double step = 0.0;
737 
738  // Single loop for common exit point
739  do {
740 
741  // Initialise iteration parameters
742  GVector* grad = fct.gradient();
743  GMatrixSparse* curvature = fct.curvature();
744 
745  // Save function value, gradient and curvature matrix
746  double save_value = m_value;
747  GVector save_grad = GVector(*grad);
748  GMatrixSparse save_curvature = GMatrixSparse(*curvature);
749 
750  // Save parameter values in vector
751  GVector save_pars(m_npars);
752  for (int ipar = 0; ipar < m_npars; ++ipar) {
753  save_pars[ipar] = pars[ipar]->factor_value();
754  }
755 
756  // Setup matrix and vector for curvature computation
757  for (int ipar = 0; ipar < m_npars; ++ipar) {
758  (*curvature)(ipar,ipar) *= (1.0 + m_lambda);
759  (*grad)[ipar] = -(*grad)[ipar];
760  }
761 
762  // Debug option: dump gradient and curvature matrix
763  #if defined(G_DEBUG_ITER)
764  std::cout << "Gradient : ";
765  for (int ipar = 0; ipar < m_npars; ++ipar) {
766  std::cout << (*grad)[ipar] << " ";
767  }
768  std::cout << std::endl;
769  std::cout << "Curvature: ";
770  for (int ipar = 0; ipar < m_npars; ++ipar) {
771  for (int jpar = 0; jpar < m_npars; ++jpar) {
772  std::cout << (*curvature)(ipar,jpar) << " ";
773  }
774  }
775  std::cout << std::endl;
776  #endif
777 
778  // Solve: curvature * X = grad. Handle matrix problems
779  try {
780  *grad = curvature->solve(*grad);
781  }
782  catch (GException::matrix_zero &e) {
784  if (m_logger != NULL) {
785  *m_logger << "GOptimizerLM::iteration: "
786  << "All curvature matrix elements are zero."
787  << std::endl;
788  }
789  continue;
790  }
793  if (m_logger != NULL) {
794  *m_logger << "GOptimizerLM::iteration: "
795  << "Curvature matrix not positive definite."
796  << std::endl;
797  }
798  continue;
799  }
800  catch (std::exception &e) {
801  throw;
802  }
803 
804  // Get LM step size
805  step = step_size(*grad, pars);
806 
807  // Debug option: dump step size
808  #if defined(G_DEBUG_ITER)
809  std::cout << "Step size = " << step << std::endl;
810  #endif
811 
812  // Derive new parameter vector
813  for (int ipar = 0; ipar < m_npars; ++ipar) {
814 
815  // Consider only free parameters
816  if (pars[ipar]->is_free()) {
817 
818  // Get actual parameter value and limits
819  double p = pars[ipar]->factor_value();
820  double p_min = pars[ipar]->factor_min();
821  double p_max = pars[ipar]->factor_max();
822 
823  // Compute new parameter value
824  p += (*grad)[ipar] * step;
825 
826  // Debug option: dump new parameter value
827  #if defined(G_DEBUG_ITER)
828  std::cout << "Trial factor value of ";
829  std::cout << pars[ipar]->name() << " = ";
830  std::cout << p << std::endl;
831  #endif
832 
833  // Constrain parameter to within the valid range
834  if (pars[ipar]->has_min() && p < p_min) {
835  if (m_hit_minimum[ipar] >= m_max_hit) {
836  if (m_logger != NULL) {
837  *m_logger << " Parameter \"" << pars[ipar]->name();
838  *m_logger << "\" hits minimum " << p_min << " more than ";
839  *m_logger << m_max_hit << " times.";
840  *m_logger << " Fix parameter at minimum for now." << std::endl;
841  }
842  m_par_freeze[ipar] = true;
843  pars[ipar]->fix();
844  }
845  else {
846  if (m_logger != NULL) {
847  *m_logger << " Parameter \"" << pars[ipar]->name();
848  *m_logger << "\" hits minimum: ";
849  *m_logger << p << " < " << p_min;
850  *m_logger << " (" << m_hit_minimum[ipar]+1 << ")" << std::endl;
851  }
852  }
853  m_hit_minimum[ipar]++;
854  p = p_min;
855  }
856  else if (pars[ipar]->has_max() && p > p_max) {
857  if (m_hit_maximum[ipar] >= m_max_hit) {
858  if (m_logger != NULL) {
859  *m_logger << " Parameter \"" << pars[ipar]->name();
860  *m_logger << "\" hits maximum " << p_max << " more than ";
861  *m_logger << m_max_hit << " times.";
862  *m_logger << " Fix parameter at maximum for now." << std::endl;
863  }
864  m_par_freeze[ipar] = true;
865  pars[ipar]->fix();
866  }
867  else {
868  if (m_logger != NULL) {
869  *m_logger << " Parameter \"" << pars[ipar]->name();
870  *m_logger << "\" hits maximum: ";
871  *m_logger << p << " > " << p_max;
872  *m_logger << " (" << m_hit_maximum[ipar]+1 << ")" << std::endl;
873  }
874  }
875  m_hit_maximum[ipar]++;
876  p = p_max;
877  }
878  else {
879  m_hit_minimum[ipar] = 0;
880  m_hit_maximum[ipar] = 0;
881  }
882 
883  // Set new parameter value
884  pars[ipar]->factor_value(p);
885 
886  // Debug option: dump new parameter value
887  #if defined(G_DEBUG_ITER)
888  std::cout << "New value of ";
889  std::cout << pars[ipar]->name() << " = ";
890  std::cout << pars[ipar]->value() << std::endl;
891  #endif
892 
893  } // endif: Parameter was free
894 
895  } // endfor: computed new parameter vector
896 
897  // Evaluate function at new parameters
898  fct.eval(pars);
899 
900  // If a free parameter has a zero diagonal element in the curvature
901  // matrix then remove this parameter definitely from the fit as it
902  // otherwise will block the fit. The problem appears in the unbinned
903  // fitting where parameter gradients may be zero (due to the truncation
904  // of the PSF), but the Npred gradient is not zero. In principle we
905  // could use the Npred gradient for fitting (I guess), but I still
906  // have to figure out how ... (the diagonal loading was not so
907  // successful as it faked early convergence)
908  for (int ipar = 0; ipar < m_npars; ++ipar) {
909  if (pars[ipar]->is_free()) {
910  if ((*fct.curvature())(ipar,ipar) == 0.0) {
911  if (m_logger != NULL) {
912  *m_logger << " Parameter \"" << pars[ipar]->name();
913  *m_logger << "\" has zero curvature.";
914  *m_logger << " Fix parameter." << std::endl;
915  }
916  m_par_remove[ipar] = true;
917  pars[ipar]->fix();
918  }
919  }
920  }
921 
922  // Fetch new pointers since eval will allocate new memory
923  grad = fct.gradient();
924  curvature = fct.curvature();
925 
926  // Retrieve new function value
927  m_value = fct.value();
928 
929  // Debug option: dump new function value
930  #if defined(G_DEBUG_ITER)
931  std::cout << "Function value " << m_value;
932  std::cout << " (was " << save_value << ")" << std::endl;
933  #endif
934 
935  // Determine how many parameters have changed
936  int par_change = 0;
937  for (int ipar = 0; ipar < m_npars; ++ipar) {
938  if (pars[ipar]->factor_value() != save_pars[ipar]) {
939  par_change++;
940  }
941  }
942 
943  // If the function has decreased then accept the new solution
944  // and decrease lambda ...
945  m_delta = save_value - m_value;
946  if (m_delta > 0.0) {
948  }
949 
950  // ... if function is identical then accept new solution. If the
951  // parameters have changed then increase lambda, otherwise decrease
952  // lambda
953  else if (m_delta == 0.0) {
954  if (par_change) {
956  }
957  else {
959  }
960  }
961 
962  // ... if function worsened by an amount smaller than m_accept_dec
963  // then accept new solution and increase lambda for the next
964  // iteration. This kluge which may help to overcome some local
965  // minimum will only for m_max_dec times at maximum to avoid
966  // keeping oscillating in a loop.
967  else if ((m_num_dec < m_max_dec) && (m_delta >= -std::abs(m_accept_dec))) {
969  m_num_dec++;
970  }
971 
972  // ... otherwise, if the statistics did not improve then use old
973  // parameters and increase lamdba. Restore also the best statistics
974  // value that was reached so far, the gradient vector and the curve
975  // matrix.
976  else {
978  m_value = save_value;
979  *grad = save_grad;
980  *curvature = save_curvature;
981  for (int ipar = 0; ipar < m_npars; ++ipar) {
982  pars[ipar]->factor_value(save_pars[ipar]);
983  }
984  }
985 
986  } while (0); // endwhile: main loop
987 
988  // Debug option: dump trailer
989  #if defined(G_DEBUG_ITER)
990  std::cout << "GOptimizerLM::iteration: exit" << std::endl;
991  #endif
992 
993  // Return step size
994  return step;
995 
996 }
997 
998 
999 /***********************************************************************//**
1000  * @brief Return LM step size
1001  *
1002  * @param[in] grad Function gradient.
1003  * @param[in] pars Function parameters.
1004  *
1005  * Determine the size of the LM step. By default a step size of 1 is taken.
1006  * If m_step_adjust=true then the step size will be estimated so that the
1007  * next parameter vector should stay within the parameter boundaries
1008  * (provided that boundaries exist).
1009  ***************************************************************************/
1010 double GOptimizerLM::step_size(const GVector& grad, const GOptimizerPars& pars)
1011 {
1012  // Initialise step size
1013  double step = 1.0;
1014 
1015  // Check if we should reduce the step size
1016  if (m_step_adjust) {
1017 
1018  // Initialise the parameter index that constrains most the fit
1019  int ipar_bnd = -1;
1020 
1021  // Loop over all parameters
1022  for (int ipar = 0; ipar < pars.size(); ++ipar) {
1023 
1024  // Get parameter attributes
1025  double p = pars[ipar]->factor_value();
1026  double p_min = pars[ipar]->factor_min();
1027  double p_max = pars[ipar]->factor_max();
1028  double delta = grad[ipar];
1029 
1030  // Check if a parameter minimum requires a reduced step size
1031  if (pars[ipar]->has_min()) {
1032  double step_min = (delta < 0.0) ? (p_min - p)/delta : 1.0;
1033  if (step_min > 0.0) {
1034  if (step_min < step) {
1035  if (!m_hit_boundary[ipar]) {
1036  ipar_bnd = ipar;
1037  step = step_min;
1038  }
1039  }
1040  else if (m_hit_boundary[ipar]) {
1041  m_hit_boundary[ipar] = false;
1042  if (m_logger != NULL) {
1043  *m_logger << " Parameter \"";
1044  *m_logger << pars[ipar]->name();
1045  *m_logger << "\" does not drive optimization step anymore.";
1046  *m_logger << std::endl;
1047  }
1048  }
1049  }
1050  }
1051 
1052  // Check if a parameter maximum requires a reduced step size
1053  if (pars[ipar]->has_max()) {
1054  double step_max = (delta > 0.0) ? (p_max - p)/delta : 1.0;
1055  if (step_max > 0.0) {
1056  if (step_max < step) {
1057  if (!m_hit_boundary[ipar]) {
1058  ipar_bnd = ipar;
1059  step = step_max;
1060  }
1061  }
1062  else if (m_hit_boundary[ipar]) {
1063  m_hit_boundary[ipar] = false;
1064  if (m_logger != NULL) {
1065  *m_logger << " Parameter \"";
1066  *m_logger << pars[ipar]->name();
1067  *m_logger << "\" does not drive optimization step anymore.";
1068  *m_logger << std::endl;
1069  }
1070  }
1071  }
1072  }
1073 
1074  } // endfor: looped over all parameters
1075 
1076  // Signal if a parameter is driving the optimization step
1077  if (ipar_bnd != -1) {
1078  m_hit_boundary[ipar_bnd] = true;
1079  if (m_logger != NULL) {
1080  *m_logger << " Parameter \"";
1081  *m_logger << pars[ipar_bnd]->name();
1082  *m_logger << "\" drives optimization step (step=";
1083  *m_logger << step << ")" << std::endl;
1084  }
1085  }
1086 
1087  } // endif: automatic step size adjustment requested
1088 
1089  // Return step size
1090  return step;
1091 }
#define G_LM_SINGULAR
GOptimizerLM & operator=(const GOptimizerLM &opt)
Assignment operator.
const double & lambda_inc(void) const
Return lambda increment value.
double step_size(const GVector &grad, const GOptimizerPars &pars)
Return LM step size.
Optimizer function abstract base class.
void init_members(void)
Initialise class members.
Sparse matrix class interface definition.
#define G_LM_STALLED
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
virtual GMatrixSparse * curvature(void)=0
Optimizer parameter container class.
virtual void optimize(GOptimizerFunction &fct, GOptimizerPars &pars)
Optimize function parameters.
double m_lambda_start
Initial start value.
double m_lambda
Actual lambda.
void copy_members(const GOptimizerLM &opt)
Copy class members.
int m_iter
Iteration.
double m_eps
Absolute precision.
std::vector< bool > m_hit_boundary
Bookkeeping array for boundary hits.
virtual void eval(const GOptimizerPars &pars)=0
int m_npars
Number of parameters.
Gammalib tools definition.
#define G_LM_CONVERGED
GOptimizerLM(void)
Void constructor.
virtual GOptimizer & operator=(const GOptimizer &opt)
Assignment operator.
Definition: GOptimizer.cpp:101
std::vector< bool > m_par_remove
Bookkeeping of parameter removal.
int m_max_dec
Maximum number of function decrease.
#define G_LM_BAD_ERRORS
double m_value
Actual function value.
GLog * m_logger
Pointer to optional logger.
Information logger interface definition.
Definition: GLog.hpp:62
virtual GOptimizerLM * clone(void) const
Clone object.
int m_nfree
Number of free parameters.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
GVector solve(const GVector &vector) const
Solves linear matrix equation.
bool m_step_adjust
Adjust step size to boundaries.
int m_status
Fit status.
virtual void clear(void)
Clear object.
double m_delta
Function improvement.
double m_accept_dec
Acceptable function decrease.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
#define G_LM_NOT_POSTIVE_DEFINITE
virtual ~GOptimizerLM(void)
Destructor.
double m_lambda_dec
Lambda decrease.
Abstract optimizer abstract base class.
Definition: GOptimizer.hpp:54
std::vector< bool > m_par_freeze
Bookkeeping of parameter freeze.
GChatter
Definition: GTypemaps.hpp:33
double m_lambda_inc
Lambda increase.
Levenberg Marquardt optimizer class interface definition.
int m_num_dec
Number of function decreases.
void free_members(void)
Delete class members.
int m_max_iter
Maximum number of iterations.
int m_max_hit
Maximum number of successive hits.
std::vector< int > m_hit_maximum
Bookkeeping of successive maximum hits.
void init_members(void)
Initialise class members.
Definition: GOptimizer.cpp:137
virtual GVector * gradient(void)=0
GVector cholesky_solver(const GVector &vector, const bool &compress=true) const
Cholesky solver.
int nfree(void) const
Return number of free parameters.
virtual double value(void) const =0
virtual void errors(GOptimizerFunction &fct, GOptimizerPars &pars)
Compute parameter uncertainties.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print optimizer information.
Exception handler interface definition.
Vector class.
Definition: GVector.hpp:46
int size(void) const
Return number of parameters in container.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
std::vector< int > m_hit_minimum
Bookkeeping of successive minimum hits.
GMatrixSparse cholesky_decompose(const bool &compress=true) const
Return Cholesky decomposition.
void free_members(void)
Delete class members.
Definition: GOptimizer.cpp:159
Levenberg Marquardt optimizer class.
int m_max_stall
Maximum number of stalls.
virtual int status(void) const
Return optimizer status.
double iteration(GOptimizerFunction &fct, GOptimizerPars &pars)
Perform one LM iteration.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413