SourceXtractorPlusPlus  0.8
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GSLEngine.cpp
Go to the documentation of this file.
1 
7 #include <gsl/gsl_multifit_nlinear.h>
8 #include <gsl/gsl_blas.h>
10 #include <iostream>
13 
14 
15 namespace ModelFitting {
16 
17 
18 static std::shared_ptr<LeastSquareEngine> createLevmarEngine(unsigned max_iterations) {
19  return std::make_shared<GSLEngine>(max_iterations);
20 }
21 
23 
24 GSLEngine::GSLEngine(int itmax, double xtol, double gtol, double ftol, double delta):
25  m_itmax{itmax}, m_xtol{xtol}, m_gtol{gtol}, m_ftol{ftol}, m_delta{delta} {
26  // Prevent an abort() call from GSL, we handle the return code ourselves
27  gsl_set_error_handler_off();
28 }
29 
30 // Provide an iterator for gsl_vector
31 class GslVectorIterator: public std::iterator<std::output_iterator_tag, double> {
32 private:
33  gsl_vector *m_v;
34  size_t m_i;
35 
36 public:
37 
38  GslVectorIterator(gsl_vector *v): m_v{v}, m_i{0} {}
39 
40  GslVectorIterator(const GslVectorIterator&) = default;
41 
43  ++m_i;
44  return *this;
45  }
46 
48  auto c = GslVectorIterator(*this);
49  ++m_i;
50  return c;
51  }
52 
53  double operator*() const {
54  return gsl_vector_get(m_v, m_i);
55  }
56 
57  double &operator*() {
58  return *gsl_vector_ptr(m_v, m_i);
59  }
60 };
61 
62 // Provide a const iterator for gsl_vector
63 class GslVectorConstIterator: public std::iterator<std::input_iterator_tag, double> {
64 private:
65  const gsl_vector *m_v;
66  size_t m_i;
67 
68 public:
69 
70  GslVectorConstIterator(const gsl_vector *v): m_v{v}, m_i{0} {}
71 
73 
75  ++m_i;
76  return *this;
77  }
78 
80  auto c = GslVectorConstIterator(*this);
81  ++m_i;
82  return c;
83  }
84 
85  double operator*() const {
86  return gsl_vector_get(m_v, m_i);
87  }
88 };
89 
90 
92  ModelFitting::ResidualEstimator& residual_estimator) {
93  // Create a tuple which keeps the references to the given manager and estimator
94  // If we capture, we can not use the lambda for the function pointer
95  auto adata = std::tie(parameter_manager, residual_estimator);
96 
97  // Only type supported by GSL
98  const gsl_multifit_nlinear_type *type = gsl_multifit_nlinear_trust;
99 
100  // Initialize parameters
101  // NOTE: These settings are set from experimenting with the examples/runs, and are those
102  // that match closer Levmar residuals, without increasing runtime too much
103  gsl_multifit_nlinear_parameters params = gsl_multifit_nlinear_default_parameters();
104  // gsl_multifit_nlinear_trs_lm is the only one that converges reasonably fast
105  // gsl_multifit_nlinear_trs_lmaccel *seems* to be faster when fitting a single source, but turns out to be slower
106  // when fitting a whole image
107  params.trs = gsl_multifit_nlinear_trs_lm;
108  // This is the only scale method that converges reasonably in SExtractor++
109  // When using gsl_multifit_nlinear_scale_more or gsl_multifit_nlinear_scale_marquardt the results are *very* bad
110  params.scale = gsl_multifit_nlinear_scale_levenberg;
111  // Others work too, but this one is faster
112  params.solver = gsl_multifit_nlinear_solver_cholesky;
113  // This is the default, and requires half the operations than GSL_MULTIFIT_NLINEAR_CTRDIFF
114  params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
115  // Passed via constructor
116  params.h_df = m_delta;
117 
118  // Allocate the workspace for the resolver. It may return null if there is no memory.
119  gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(
120  type, &params,
121  residual_estimator.numberOfResiduals(), parameter_manager.numberOfParameters()
122  );
123  if (workspace == nullptr) {
124  throw Elements::Exception() << "Insufficient memory for initializing the GSL solver";
125  }
126 
127  // Allocate space for the parameters and initialize with the guesses
128  std::vector<double> param_values (parameter_manager.numberOfParameters());
129  parameter_manager.getEngineValues(param_values.begin());
130  gsl_vector_view gsl_param_view = gsl_vector_view_array(param_values.data(), param_values.size());
131 
132  // Function to be minimized
133  auto function = [](const gsl_vector* x, void *extra, gsl_vector *f) -> int {
134  auto *extra_ptr = (decltype(adata) *) extra;
135  EngineParameterManager& pm = std::get<0>(*extra_ptr);
137  ResidualEstimator& re = std::get<1>(*extra_ptr);
139  return GSL_SUCCESS;
140  };
141  gsl_multifit_nlinear_fdf fdf;
142  fdf.f = function;
143  fdf.df = nullptr;
144  fdf.fvv = nullptr;
145  fdf.n = residual_estimator.numberOfResiduals();
146  fdf.p = parameter_manager.numberOfParameters();
147  fdf.params = &adata;
148 
149  // Setup the solver
150  gsl_multifit_nlinear_init(&gsl_param_view.vector, &fdf, workspace);
151 
152  // Initial cost
153  double chisq0;
154  gsl_vector *residual = gsl_multifit_nlinear_residual(workspace);
155  gsl_blas_ddot(residual, residual, &chisq0);
156 
157  // Solve
158  int info = 0;
159  int ret = gsl_multifit_nlinear_driver(
160  m_itmax, // Maximum number of iterations
161  m_xtol, // Tolerance for the step
162  m_gtol, // Tolerance for the gradient
163  m_ftol, // Tolerance for chi^2 (may be unused)
164  nullptr, // Callback
165  nullptr, // Callback parameters
166  &info, // Convergence information if GSL_SUCCESS
167  workspace // The solver workspace
168  );
169 
170  // Final cost
171  double chisq;
172  gsl_blas_ddot(residual, residual, &chisq);
173 
174  // Build run summary
175  std::vector<double> covariance_matrix (parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
176 
177  LeastSquareSummary summary;
178  summary.success_flag = (ret == GSL_SUCCESS);
179  summary.iteration_no = gsl_multifit_nlinear_niter(workspace);
180  summary.parameter_sigmas = {};
181 
182  // Covariance matrix. Note: Do not free J! It is owned by the workspace.
183  gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
184  gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.data(), parameter_manager.numberOfParameters(),
185  parameter_manager.numberOfParameters());
186  gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
187 
188  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
189  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
190  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
191  }
192 
193  // Levmar-compatible engine info
194  int levmar_reason = 0;
195  if (ret == GSL_SUCCESS) {
196  levmar_reason = (info == 1) ? 2 : 1;
197  }
198  else if (ret == GSL_EMAXITER) {
199  levmar_reason = 3;
200  }
201 
202  summary.underlying_framework_info = std::array<double, 10>{
203  chisq0, // ||e||_2 at initial p
204  chisq, // ||e||_2
205  gsl_blas_dnrm2(workspace->g), // ||J^T e||_inf
206  gsl_blas_dnrm2(workspace->dx), // ||Dp||_2
207  0., // mu/max[J^T J]_ii
208  static_cast<double>(summary.iteration_no), // Number of iterations
209  static_cast<double>(levmar_reason), // Reason for finishing
210  static_cast<double>(fdf.nevalf), // Function evaluations
211  static_cast<double>(fdf.nevaldf), // Jacobian evaluations
212  0. // Linear systems solved
213  };
214 
215  // Free resources and return
216  gsl_multifit_nlinear_free(workspace);
217  return summary;
218 }
219 
220 } // end namespace ModelFitting
GslVectorConstIterator(const gsl_vector *v)
Definition: GSLEngine.cpp:70
T tie(T...args)
GslVectorIterator & operator++()
Definition: GSLEngine.cpp:42
Class containing the summary information of solving a least square minimization problem.
void populateResiduals(DoubleIter output_iter) const
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void getEngineValues(DoubleIter output_iter) const
Returns the engine values of the managed parameters.
GSLEngine(int itmax=1000, double xtol=1e-8, double gtol=1e-8, double ftol=1e-8, double delta=1e-4)
Constructs a new instance of the engine.
Definition: GSLEngine.cpp:24
LeastSquareSummary solveProblem(EngineParameterManager &parameter_manager, ResidualEstimator &residual_estimator) override
Definition: GSLEngine.cpp:91
GslVectorConstIterator & operator++()
Definition: GSLEngine.cpp:74
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
T push_back(T...args)
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:18
std::size_t numberOfResiduals() const
GslVectorIterator(gsl_vector *v)
Definition: GSLEngine.cpp:38
static LeastSquareEngineManager::StaticEngine levmar_engine
Definition: GSLEngine.cpp:22
GslVectorIterator operator++(int)
Definition: GSLEngine.cpp:47
STL class.
Class responsible for managing the parameters the least square engine minimizes.
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
T sqrt(T...args)
GslVectorConstIterator operator++(int)
Definition: GSLEngine.cpp:79
Provides to the LeastSquareEngine the residual values.
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.