SourceXtractorPlusPlus  0.15
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 
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 
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  int info = 0;
185  int ret = gsl_multifit_nlinear_driver(
186  m_itmax, // Maximum number of iterations
187  m_xtol, // Tolerance for the step
188  m_gtol, // Tolerance for the gradient
189  m_ftol, // Tolerance for chi^2 (may be unused)
190  nullptr, // Callback
191  nullptr, // Callback parameters
192  &info, // Convergence information if GSL_SUCCESS
193  workspace // The solver workspace
194  );
195 
196  // Final cost
197  double chisq;
198  gsl_blas_ddot(residual, residual, &chisq);
199 
200  // Build run summary
201  std::vector<double> covariance_matrix(
202  parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
203 
204  LeastSquareSummary summary;
205  summary.status_flag = getStatusFlag(ret);
206  summary.engine_stop_reason = ret;
207  summary.iteration_no = gsl_multifit_nlinear_niter(workspace);
208  summary.parameter_sigmas = {};
209 
210  // Covariance matrix. Note: Do not free J! It is owned by the workspace.
211  gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
212  gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.data(), parameter_manager.numberOfParameters(),
213  parameter_manager.numberOfParameters());
214  gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
215 
216  // We have done an unweighted minimization, so, from the documentation, we need
217  // to multiply the covariance matrix by the variance of the residuals
218  // See: https://www.gnu.org/software/gsl/doc/html/nls.html#covariance-matrix-of-best-fit-parameters
219  double sigma2 = 0;
220  for (size_t i = 0; i < residual->size; ++i) {
221  auto v = gsl_vector_get(residual, i);
222  sigma2 += v * v;
223  }
224  sigma2 /= (fdf.n - fdf.p);
225 
226  for (auto ci = covariance_matrix.begin(); ci != covariance_matrix.end(); ++ci) {
227  *ci *= sigma2;
228  }
229 
230  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
231  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
232  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
233  }
234 
235  // Levmar-compatible engine info
236  int levmar_reason = 0;
237  if (ret == GSL_SUCCESS) {
238  levmar_reason = (info == 1) ? 2 : 1;
239  }
240  else if (ret == GSL_EMAXITER) {
241  levmar_reason = 3;
242  }
243 
245  chisq0, // ||e||_2 at initial p
246  chisq, // ||e||_2
247  gsl_blas_dnrm2(workspace->g), // ||J^T e||_inf
248  gsl_blas_dnrm2(workspace->dx), // ||Dp||_2
249  0., // mu/max[J^T J]_ii
250  static_cast<double>(summary.iteration_no), // Number of iterations
251  static_cast<double>(levmar_reason), // Reason for finishing
252  static_cast<double>(fdf.nevalf), // Function evaluations
253  static_cast<double>(fdf.nevaldf), // Jacobian evaluations
254  0. // Linear systems solved
255  };
256 
257  // Free resources and return
258  gsl_multifit_nlinear_free(workspace);
259  return summary;
260 }
261 
262 } // end namespace ModelFitting
ModelFitting::GslVectorIterator::GslVectorIterator
GslVectorIterator(gsl_vector *v)
Definition: GSLEngine.cpp:54
ModelFitting::ResidualEstimator
Provides to the LeastSquareEngine the residual values.
Definition: ResidualEstimator.h:50
ModelFitting::LeastSquareEngineManager::StaticEngine
Definition: LeastSquareEngineManager.h:88
ModelFitting::getStatusFlag
static LeastSquareSummary::StatusFlag getStatusFlag(int ret)
Definition: GSLEngine.cpp:106
ModelFitting::GslVectorIterator::operator*
double operator*() const
Definition: GSLEngine.cpp:69
ModelFitting::EngineParameterManager::convertCovarianceMatrixToWorldSpace
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
Definition: EngineParameterManager.cpp:37
ModelFitting::GslVectorIterator::operator++
GslVectorIterator operator++(int)
Definition: GSLEngine.cpp:63
std::shared_ptr
STL class.
ModelFitting::LeastSquareSummary::engine_stop_reason
int engine_stop_reason
Engine-specific reason for stopping the fitting.
Definition: LeastSquareSummary.h:54
ModelFitting::GslVectorIterator
Definition: GSLEngine.cpp:47
GSLEngine.h
ModelFitting::LeastSquareSummary::MAX_ITER
@ MAX_ITER
Definition: LeastSquareSummary.h:41
std::vector< double >
std::vector::size
T size(T... args)
std::iterator
ModelFitting::GslVectorConstIterator::GslVectorConstIterator
GslVectorConstIterator(const gsl_vector *v)
Definition: GSLEngine.cpp:86
ModelFitting::GslVectorIterator::operator*
double & operator*()
Definition: GSLEngine.cpp:73
ModelFitting::EngineParameterManager::getEngineValues
void getEngineValues(DoubleIter output_iter) const
Returns the engine values of the managed parameters.
ModelFitting::ResidualEstimator::numberOfResiduals
std::size_t numberOfResiduals() const
Definition: ResidualEstimator.cpp:34
ModelFitting::LeastSquareSummary
Class containing the summary information of solving a least square minimization problem.
Definition: LeastSquareSummary.h:38
std::sqrt
T sqrt(T... args)
ModelFitting::EngineParameterManager::updateEngineValues
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
std::tie
T tie(T... args)
std::vector::push_back
T push_back(T... args)
ModelFitting::GSLEngine::GSLEngine
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
ModelFitting::GslVectorConstIterator::operator++
GslVectorConstIterator & operator++()
Definition: GSLEngine.cpp:90
ModelFitting::GSLEngine::solveProblem
LeastSquareSummary solveProblem(EngineParameterManager &parameter_manager, ResidualEstimator &residual_estimator) override
Definition: GSLEngine.cpp:117
ModelFitting::LeastSquareSummary::status_flag
StatusFlag status_flag
Flag indicating if the minimization was successful.
Definition: LeastSquareSummary.h:45
ModelFitting::LeastSquareSummary::iteration_no
size_t iteration_no
The number of iterations.
Definition: LeastSquareSummary.h:48
Exception.h
ModelFitting::LeastSquareSummary::ERROR
@ ERROR
Definition: LeastSquareSummary.h:41
LeastSquareEngineManager.h
ModelFitting::GSLEngine::m_ftol
double m_ftol
Definition: GSLEngine.h:71
ModelFitting::ResidualEstimator::populateResiduals
void populateResiduals(DoubleIter output_iter) const
std::array
STL class.
ModelFitting::EngineParameterManager
Class responsible for managing the parameters the least square engine minimizes.
Definition: EngineParameterManager.h:61
ModelFitting::GslVectorConstIterator::operator*
double operator*() const
Definition: GSLEngine.cpp:101
Elements::Exception
ModelFitting::levmar_engine
static LeastSquareEngineManager::StaticEngine levmar_engine
Definition: GSLEngine.cpp:38
ModelFitting::LeastSquareSummary::SUCCESS
@ SUCCESS
Definition: LeastSquareSummary.h:41
ModelFitting::GSLEngine::m_gtol
double m_gtol
Definition: GSLEngine.h:71
ModelFitting::createLevmarEngine
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:34
ModelFitting::GslVectorConstIterator::m_v
const gsl_vector * m_v
Definition: GSLEngine.cpp:81
ModelFitting::GSLEngine::m_delta
double m_delta
Definition: GSLEngine.h:71
ModelFitting::GSLEngine::m_xtol
double m_xtol
Definition: GSLEngine.h:71
ModelFitting::GSLEngine::m_itmax
int m_itmax
Definition: GSLEngine.h:70
ModelFitting::GslVectorIterator::GslVectorIterator
GslVectorIterator(const GslVectorIterator &)=default
ModelFitting::EngineParameterManager::numberOfParameters
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.
Definition: EngineParameterManager.cpp:33
x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
Definition: MoffatModelFittingTask.cpp:94
std::vector::begin
T begin(T... args)
ModelFitting::GslVectorConstIterator
Definition: GSLEngine.cpp:79
ModelFitting::GslVectorIterator::m_v
gsl_vector * m_v
Definition: GSLEngine.cpp:49
ModelFitting::GslVectorConstIterator::m_i
size_t m_i
Definition: GSLEngine.cpp:82
ModelFitting::LeastSquareSummary::underlying_framework_info
boost::any underlying_framework_info
Definition: LeastSquareSummary.h:63
ModelFitting::GslVectorConstIterator::operator++
GslVectorConstIterator operator++(int)
Definition: GSLEngine.cpp:95
std::vector::end
T end(T... args)
ModelFitting::GslVectorIterator::m_i
size_t m_i
Definition: GSLEngine.cpp:50
ModelFitting::GslVectorIterator::operator++
GslVectorIterator & operator++()
Definition: GSLEngine.cpp:58
ModelFitting
Definition: AsinhChiSquareComparator.h:30
ModelFitting::GslVectorConstIterator::GslVectorConstIterator
GslVectorConstIterator(const GslVectorConstIterator &)=default
std::vector::data
T data(T... args)
ModelFitting::LeastSquareSummary::StatusFlag
StatusFlag
Definition: LeastSquareSummary.h:40
ModelFitting::LeastSquareSummary::parameter_sigmas
std::vector< double > parameter_sigmas
1-sigma margin of error for all the parameters
Definition: LeastSquareSummary.h:51