SourceXtractorPlusPlus  0.14
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevmarEngine.cpp
Go to the documentation of this file.
1 
23 #include <cmath>
24 #include <mutex>
25 
26 #include <levmar.h>
28 #include <ElementsKernel/Logging.h>
31 
32 
33 namespace ModelFitting {
34 
36 
37 static std::shared_ptr<LeastSquareEngine> createLevmarEngine(unsigned max_iterations) {
38  return std::make_shared<LevmarEngine>(max_iterations);
39 }
40 
42 
43 LevmarEngine::LevmarEngine(size_t itmax, double tau, double epsilon1,
44  double epsilon2, double epsilon3, double delta)
45  : m_itmax{itmax}, m_opts{tau, epsilon1, epsilon2, epsilon3, delta} {
46 #ifdef LINSOLVERS_RETAIN_MEMORY
47  logger.warn() << "Using a non thread safe levmar! Parallelism will be reduced.";
48 #endif
49 }
50 
51 LevmarEngine::~LevmarEngine() = default;
52 
53 
54 #ifdef LINSOLVERS_RETAIN_MEMORY
55 // If the Levmar library is not configured for multithreading, this mutex is used to ensure only one thread
56 // at a time can enter levmar
57 namespace {
58  std::mutex levmar_mutex;
59 }
60 #endif
61 
63  ResidualEstimator& residual_estimator) {
64  // Create a tuple which keeps the references to the given manager and estimator
65  auto adata = std::tie(parameter_manager, residual_estimator);
66 
67  // The function which is called by the levmar loop
68  auto levmar_res_func = [](double *p, double *hx, int, int, void *extra) {
69 #ifdef LINSOLVERS_RETAIN_MEMORY
70  levmar_mutex.unlock();
71 #endif
72  auto* extra_ptr = (decltype(adata)*)extra;
73  EngineParameterManager& pm = std::get<0>(*extra_ptr);
74  pm.updateEngineValues(p);
75  ResidualEstimator& re = std::get<1>(*extra_ptr);
76  re.populateResiduals(hx);
77 
78 #ifdef LINSOLVERS_RETAIN_MEMORY
79  levmar_mutex.lock();
80 #endif
81  };
82 
83  // Create the vector which will be used for keeping the parameter values
84  // and initialize it to the current values of the parameters
85  std::vector<double> param_values (parameter_manager.numberOfParameters());
86  parameter_manager.getEngineValues(param_values.begin());
87 
88  // Create a vector for getting the information of the minimization
90 
91  std::vector<double> covariance_matrix (parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
92 
93 #ifdef LINSOLVERS_RETAIN_MEMORY
94  levmar_mutex.lock();
95 #endif
96  // Call the levmar library
97  auto res = dlevmar_dif(levmar_res_func, // The function called from the levmar algorithm
98  param_values.data(), // The pointer where the parameter values are
99  NULL, // We don't use any measurement vector
100  parameter_manager.numberOfParameters(), // The number of free parameters
101  residual_estimator.numberOfResiduals(), // The number of residuals
102  m_itmax, // The maximum number of iterations
103  m_opts.data(), // The minimization options
104  info.data(), // Where the information of the minimization is stored
105  NULL, // Working memory is allocated internally
106  covariance_matrix.data(),
107  &adata // No additional data needed
108  );
109 #ifdef LINSOLVERS_RETAIN_MEMORY
110  levmar_mutex.unlock();
111 #endif
112 
113  // Create and return the summary object
114  LeastSquareSummary summary {};
115 
116  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
117  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
118  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
119  }
120 
121  summary.success_flag = (res != -1);
122  summary.iteration_no = info[5];
123  summary.underlying_framework_info = info;
124  return summary;
125 }
126 
127 } // end of namespace ModelFitting
virtual ~LevmarEngine()
Destructor.
T tie(T...args)
Class containing the summary information of solving a least square minimization problem.
void populateResiduals(DoubleIter output_iter) const
void getEngineValues(DoubleIter output_iter) const
Returns the engine values of the managed parameters.
static Elements::Logging logger
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
T push_back(T...args)
T data(T...args)
void warn(const std::string &logMessage)
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:34
std::vector< double > m_opts
Definition: LevmarEngine.h:72
LevmarEngine(size_t itmax=1000, double tau=1E-3, double epsilon1=1E-8, double epsilon2=1E-8, double epsilon3=1E-8, double delta=1E-4)
Constructs a new instance of the engine.
std::size_t numberOfResiduals() const
static LeastSquareEngineManager::StaticEngine levmar_engine
Definition: GSLEngine.cpp:38
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)
LeastSquareSummary solveProblem(EngineParameterManager &parameter_manager, ResidualEstimator &residual_estimator) override
Provides to the LeastSquareEngine the residual values.
static Logging getLogger(const std::string &name="")
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.