SourceXtractorPlusPlus  0.16
Please provide a description of the project.
GSLEngine.cpp
Go to the documentation of this file.
1 
26 #include <chrono>
27 #include <gsl/gsl_blas.h>
28 #include <gsl/gsl_multifit_nlinear.h>
29 #include <iostream>
30 
31 namespace ModelFitting {
32 
33 
34 static std::shared_ptr<LeastSquareEngine> createGslEngine(unsigned max_iterations) {
35  return std::make_shared<GSLEngine>(max_iterations);
36 }
37 
39 
40 GSLEngine::GSLEngine(int itmax, double xtol, double gtol, double ftol, double delta):
41  m_itmax{itmax}, m_xtol{xtol}, m_gtol{gtol}, m_ftol{ftol}, m_delta{delta} {
42  // Prevent an abort() call from GSL, we handle the return code ourselves
43  gsl_set_error_handler_off();
44 }
45 
46 // Provide an iterator for gsl_vector
47 class GslVectorIterator: public std::iterator<std::output_iterator_tag, double> {
48 private:
49  gsl_vector *m_v;
50  size_t m_i;
51 
52 public:
53 
54  explicit GslVectorIterator(gsl_vector *v): m_v{v}, m_i{0} {}
55 
57 
59  ++m_i;
60  return *this;
61  }
62 
64  auto c = GslVectorIterator(*this);
65  ++m_i;
66  return c;
67  }
68 
69  double operator*() const {
70  return gsl_vector_get(m_v, m_i);
71  }
72 
73  double &operator*() {
74  return *gsl_vector_ptr(m_v, m_i);
75  }
76 };
77 
78 // Provide a const iterator for gsl_vector
79 class GslVectorConstIterator: public std::iterator<std::input_iterator_tag, double> {
80 private:
81  const gsl_vector *m_v;
82  size_t m_i;
83 
84 public:
85 
86  explicit GslVectorConstIterator(const gsl_vector *v): m_v{v}, m_i{0} {}
87 
89 
91  ++m_i;
92  return *this;
93  }
94 
96  auto c = GslVectorConstIterator(*this);
97  ++m_i;
98  return c;
99  }
100 
101  double operator*() const {
102  return gsl_vector_get(m_v, m_i);
103  }
104 };
105 
107  switch (ret) {
108  case GSL_SUCCESS:
110  case GSL_EMAXITER:
112  default:
114  }
115 }
116 
118  ModelFitting::ResidualEstimator& residual_estimator) {
119  // Create a tuple which keeps the references to the given manager and estimator
120  // If we capture, we can not use the lambda for the function pointer
121  auto adata = std::tie(parameter_manager, residual_estimator);
122 
123  // Only type supported by GSL
124  const gsl_multifit_nlinear_type *type = gsl_multifit_nlinear_trust;
125 
126  // Initialize parameters
127  // NOTE: These settings are set from experimenting with the examples/runs, and are those
128  // that match closer Levmar residuals, without increasing runtime too much
129  gsl_multifit_nlinear_parameters params = gsl_multifit_nlinear_default_parameters();
130  // gsl_multifit_nlinear_trs_lm is the only one that converges reasonably fast
131  // gsl_multifit_nlinear_trs_lmaccel *seems* to be faster when fitting a single source, but turns out to be slower
132  // when fitting a whole image
133  params.trs = gsl_multifit_nlinear_trs_lm;
134  // This is the only scale method that converges reasonably in SExtractor++
135  // When using gsl_multifit_nlinear_scale_more or gsl_multifit_nlinear_scale_marquardt the results are *very* bad
136  params.scale = gsl_multifit_nlinear_scale_levenberg;
137  // Others work too, but this one is faster
138  params.solver = gsl_multifit_nlinear_solver_cholesky;
139  // This is the default, and requires half the operations than GSL_MULTIFIT_NLINEAR_CTRDIFF
140  params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
141  // Passed via constructor
142  params.h_df = m_delta;
143 
144  // Allocate the workspace for the resolver. It may return null if there is no memory.
145  gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(
146  type, &params,
147  residual_estimator.numberOfResiduals(), parameter_manager.numberOfParameters()
148  );
149  if (workspace == nullptr) {
150  throw Elements::Exception() << "Insufficient memory for initializing the GSL solver";
151  }
152 
153  // Allocate space for the parameters and initialize with the guesses
154  std::vector<double> param_values(parameter_manager.numberOfParameters());
155  parameter_manager.getEngineValues(param_values.begin());
156  gsl_vector_view gsl_param_view = gsl_vector_view_array(param_values.data(), param_values.size());
157 
158  // Function to be minimized
159  auto function = [](const gsl_vector *x, void *extra, gsl_vector *f) -> int {
160  auto *extra_ptr = (decltype(adata) *) extra;
161  EngineParameterManager& pm = std::get<0>(*extra_ptr);
163  ResidualEstimator& re = std::get<1>(*extra_ptr);
165  return GSL_SUCCESS;
166  };
167  gsl_multifit_nlinear_fdf fdf;
168  fdf.f = function;
169  fdf.df = nullptr;
170  fdf.fvv = nullptr;
171  fdf.n = residual_estimator.numberOfResiduals();
172  fdf.p = parameter_manager.numberOfParameters();
173  fdf.params = &adata;
174 
175  // Setup the solver
176  gsl_multifit_nlinear_init(&gsl_param_view.vector, &fdf, workspace);
177 
178  // Initial cost
179  double chisq0;
180  gsl_vector *residual = gsl_multifit_nlinear_residual(workspace);
181  gsl_blas_ddot(residual, residual, &chisq0);
182 
183  // Solve
184  auto start = std::chrono::steady_clock::now();
185  int info = 0;
186  int ret = gsl_multifit_nlinear_driver(
187  m_itmax, // Maximum number of iterations
188  m_xtol, // Tolerance for the step
189  m_gtol, // Tolerance for the gradient
190  m_ftol, // Tolerance for chi^2 (may be unused)
191  nullptr, // Callback
192  nullptr, // Callback parameters
193  &info, // Convergence information if GSL_SUCCESS
194  workspace // The solver workspace
195  );
197  std::chrono::duration<float> elapsed = end - start;
198 
199  // Final cost
200  double chisq;
201  gsl_blas_ddot(residual, residual, &chisq);
202 
203  // Build run summary
204  std::vector<double> covariance_matrix(
205  parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
206 
207  LeastSquareSummary summary;
208  summary.status_flag = getStatusFlag(ret);
209  summary.engine_stop_reason = ret;
210  summary.iteration_no = gsl_multifit_nlinear_niter(workspace);
211  summary.parameter_sigmas = {};
212  summary.duration = elapsed.count();
213 
214  // Covariance matrix. Note: Do not free J! It is owned by the workspace.
215  gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
216  gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.data(), parameter_manager.numberOfParameters(),
217  parameter_manager.numberOfParameters());
218  gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
219 
220  // We have done an unweighted minimization, so, from the documentation, we need
221  // to multiply the covariance matrix by the variance of the residuals
222  // See: https://www.gnu.org/software/gsl/doc/html/nls.html#covariance-matrix-of-best-fit-parameters
223  double sigma2 = 0;
224  for (size_t i = 0; i < residual->size; ++i) {
225  auto v = gsl_vector_get(residual, i);
226  sigma2 += v * v;
227  }
228  sigma2 /= (fdf.n - fdf.p);
229 
230  for (auto ci = covariance_matrix.begin(); ci != covariance_matrix.end(); ++ci) {
231  *ci *= sigma2;
232  }
233 
234  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
235  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
236  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
237  }
238 
239  // Levmar-compatible engine info
240  int levmar_reason = 0;
241  if (ret == GSL_SUCCESS) {
242  levmar_reason = (info == 1) ? 2 : 1;
243  }
244  else if (ret == GSL_EMAXITER) {
245  levmar_reason = 3;
246  }
247 
249  chisq0, // ||e||_2 at initial p
250  chisq, // ||e||_2
251  gsl_blas_dnrm2(workspace->g), // ||J^T e||_inf
252  gsl_blas_dnrm2(workspace->dx), // ||Dp||_2
253  0., // mu/max[J^T J]_ii
254  static_cast<double>(summary.iteration_no), // Number of iterations
255  static_cast<double>(levmar_reason), // Reason for finishing
256  static_cast<double>(fdf.nevalf), // Function evaluations
257  static_cast<double>(fdf.nevaldf), // Jacobian evaluations
258  0. // Linear systems solved
259  };
260 
261  // Free resources and return
262  gsl_multifit_nlinear_free(workspace);
263  return summary;
264 }
265 
266 } // end namespace ModelFitting
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
T begin(T... args)
Class responsible for managing the parameters the least square engine minimizes.
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.
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:40
LeastSquareSummary solveProblem(EngineParameterManager &parameter_manager, ResidualEstimator &residual_estimator) override
Definition: GSLEngine.cpp:117
GslVectorConstIterator & operator++()
Definition: GSLEngine.cpp:90
GslVectorConstIterator(const GslVectorConstIterator &)=default
GslVectorConstIterator(const gsl_vector *v)
Definition: GSLEngine.cpp:86
GslVectorConstIterator operator++(int)
Definition: GSLEngine.cpp:95
GslVectorIterator operator++(int)
Definition: GSLEngine.cpp:63
GslVectorIterator & operator++()
Definition: GSLEngine.cpp:58
GslVectorIterator(gsl_vector *v)
Definition: GSLEngine.cpp:54
GslVectorIterator(const GslVectorIterator &)=default
Provides to the LeastSquareEngine the residual values.
void populateResiduals(DoubleIter output_iter) const
T data(T... args)
T end(T... args)
static LeastSquareSummary::StatusFlag getStatusFlag(int ret)
Definition: GSLEngine.cpp:106
static std::shared_ptr< LeastSquareEngine > createGslEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:34
static LeastSquareEngineManager::StaticEngine gsl_engine
Definition: GSLEngine.cpp:38
T push_back(T... args)
T size(T... args)
T sqrt(T... args)
Class containing the summary information of solving a least square minimization problem.
StatusFlag status_flag
Flag indicating if the minimization was successful.
size_t iteration_no
The number of iterations.
float duration
Runtime (in seconds)
std::vector< double > parameter_sigmas
1-sigma margin of error for all the parameters
int engine_stop_reason
Engine-specific reason for stopping the fitting.
T tie(T... args)