SourceXtractorPlusPlus  0.12
Please provide a description of the project.
GSLEngine.cpp
Go to the documentation of this file.
1 
23 #include <gsl/gsl_multifit_nlinear.h>
24 #include <gsl/gsl_blas.h>
26 #include <iostream>
29 
30 
31 namespace ModelFitting {
32 
33 
34 static std::shared_ptr<LeastSquareEngine> createLevmarEngine(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  GslVectorIterator(gsl_vector *v): m_v{v}, m_i{0} {}
55 
56  GslVectorIterator(const GslVectorIterator&) = default;
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  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 
106 
108  ModelFitting::ResidualEstimator& residual_estimator) {
109  // Create a tuple which keeps the references to the given manager and estimator
110  // If we capture, we can not use the lambda for the function pointer
111  auto adata = std::tie(parameter_manager, residual_estimator);
112 
113  // Only type supported by GSL
114  const gsl_multifit_nlinear_type *type = gsl_multifit_nlinear_trust;
115 
116  // Initialize parameters
117  // NOTE: These settings are set from experimenting with the examples/runs, and are those
118  // that match closer Levmar residuals, without increasing runtime too much
119  gsl_multifit_nlinear_parameters params = gsl_multifit_nlinear_default_parameters();
120  // gsl_multifit_nlinear_trs_lm is the only one that converges reasonably fast
121  // gsl_multifit_nlinear_trs_lmaccel *seems* to be faster when fitting a single source, but turns out to be slower
122  // when fitting a whole image
123  params.trs = gsl_multifit_nlinear_trs_lm;
124  // This is the only scale method that converges reasonably in SExtractor++
125  // When using gsl_multifit_nlinear_scale_more or gsl_multifit_nlinear_scale_marquardt the results are *very* bad
126  params.scale = gsl_multifit_nlinear_scale_levenberg;
127  // Others work too, but this one is faster
128  params.solver = gsl_multifit_nlinear_solver_cholesky;
129  // This is the default, and requires half the operations than GSL_MULTIFIT_NLINEAR_CTRDIFF
130  params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
131  // Passed via constructor
132  params.h_df = m_delta;
133 
134  // Allocate the workspace for the resolver. It may return null if there is no memory.
135  gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(
136  type, &params,
137  residual_estimator.numberOfResiduals(), parameter_manager.numberOfParameters()
138  );
139  if (workspace == nullptr) {
140  throw Elements::Exception() << "Insufficient memory for initializing the GSL solver";
141  }
142 
143  // Allocate space for the parameters and initialize with the guesses
144  std::vector<double> param_values(parameter_manager.numberOfParameters());
145  parameter_manager.getEngineValues(param_values.begin());
146  gsl_vector_view gsl_param_view = gsl_vector_view_array(param_values.data(), param_values.size());
147 
148  // Function to be minimized
149  auto function = [](const gsl_vector *x, void *extra, gsl_vector *f) -> int {
150  auto *extra_ptr = (decltype(adata) *) extra;
151  EngineParameterManager& pm = std::get<0>(*extra_ptr);
153  ResidualEstimator& re = std::get<1>(*extra_ptr);
155  return GSL_SUCCESS;
156  };
157  gsl_multifit_nlinear_fdf fdf;
158  fdf.f = function;
159  fdf.df = nullptr;
160  fdf.fvv = nullptr;
161  fdf.n = residual_estimator.numberOfResiduals();
162  fdf.p = parameter_manager.numberOfParameters();
163  fdf.params = &adata;
164 
165  // Setup the solver
166  gsl_multifit_nlinear_init(&gsl_param_view.vector, &fdf, workspace);
167 
168  // Initial cost
169  double chisq0;
170  gsl_vector *residual = gsl_multifit_nlinear_residual(workspace);
171  gsl_blas_ddot(residual, residual, &chisq0);
172 
173  // Solve
174  int info = 0;
175  int ret = gsl_multifit_nlinear_driver(
176  m_itmax, // Maximum number of iterations
177  m_xtol, // Tolerance for the step
178  m_gtol, // Tolerance for the gradient
179  m_ftol, // Tolerance for chi^2 (may be unused)
180  nullptr, // Callback
181  nullptr, // Callback parameters
182  &info, // Convergence information if GSL_SUCCESS
183  workspace // The solver workspace
184  );
185 
186  // Final cost
187  double chisq;
188  gsl_blas_ddot(residual, residual, &chisq);
189 
190  // Build run summary
191  std::vector<double> covariance_matrix(
192  parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
193 
194  LeastSquareSummary summary;
195  summary.success_flag = (ret == GSL_SUCCESS);
196  summary.iteration_no = gsl_multifit_nlinear_niter(workspace);
197  summary.parameter_sigmas = {};
198 
199  // Covariance matrix. Note: Do not free J! It is owned by the workspace.
200  gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
201  gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.data(), parameter_manager.numberOfParameters(),
202  parameter_manager.numberOfParameters());
203  gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
204 
205  // We have done an unweighted minimization, so, from the documentation, we need
206  // to multiply the covariance matrix by the variance of the residuals
207  // See: https://www.gnu.org/software/gsl/doc/html/nls.html#covariance-matrix-of-best-fit-parameters
208  double sigma2 = 0;
209  for (size_t i = 0; i < residual->size; ++i) {
210  auto v = gsl_vector_get(residual, i);
211  sigma2 += v * v;
212  }
213  sigma2 /= (fdf.n - fdf.p);
214 
215  for (auto ci = covariance_matrix.begin(); ci != covariance_matrix.end(); ++ci) {
216  *ci *= sigma2;
217  }
218 
219  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
220  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
221  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
222  }
223 
224  // Levmar-compatible engine info
225  int levmar_reason = 0;
226  if (ret == GSL_SUCCESS) {
227  levmar_reason = (info == 1) ? 2 : 1;
228  }
229  else if (ret == GSL_EMAXITER) {
230  levmar_reason = 3;
231  }
232 
233  summary.underlying_framework_info = std::array<double, 10>{
234  chisq0, // ||e||_2 at initial p
235  chisq, // ||e||_2
236  gsl_blas_dnrm2(workspace->g), // ||J^T e||_inf
237  gsl_blas_dnrm2(workspace->dx), // ||Dp||_2
238  0., // mu/max[J^T J]_ii
239  static_cast<double>(summary.iteration_no), // Number of iterations
240  static_cast<double>(levmar_reason), // Reason for finishing
241  static_cast<double>(fdf.nevalf), // Function evaluations
242  static_cast<double>(fdf.nevaldf), // Jacobian evaluations
243  0. // Linear systems solved
244  };
245 
246  // Free resources and return
247  gsl_multifit_nlinear_free(workspace);
248  return summary;
249 }
250 
251 } // end namespace ModelFitting
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
GslVectorConstIterator(const gsl_vector *v)
Definition: GSLEngine.cpp:86
T tie(T... args)
GslVectorIterator & operator++()
Definition: GSLEngine.cpp:58
Class containing the summary information of solving a least square minimization problem.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
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:107
GslVectorConstIterator & operator++()
Definition: GSLEngine.cpp:90
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
T push_back(T... args)
void populateResiduals(DoubleIter output_iter) const
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:34
GslVectorIterator(gsl_vector *v)
Definition: GSLEngine.cpp:54
static LeastSquareEngineManager::StaticEngine levmar_engine
Definition: GSLEngine.cpp:38
GslVectorIterator operator++(int)
Definition: GSLEngine.cpp:63
STL class.
Class responsible for managing the parameters the least square engine minimizes.
T sqrt(T... args)
GslVectorConstIterator operator++(int)
Definition: GSLEngine.cpp:95
Provides to the LeastSquareEngine the residual values.
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.