SourceXtractorPlusPlus  0.15
Please provide a description of the project.
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 namespace {
33 
34 __attribute__((unused))
35 void printLevmarInfo(std::array<double, 10> info) {
36  std::cerr << "\nMinimization info:\n";
37  std::cerr << " ||e||_2 at initial p: " << info[0] << '\n';
38  std::cerr << " ||e||_2: " << info[1] << '\n';
39  std::cerr << " ||J^T e||_inf: " << info[2] << '\n';
40  std::cerr << " ||Dp||_2: " << info[3] << '\n';
41  std::cerr << " mu/max[J^T J]_ii: " << info[4] << '\n';
42  std::cerr << " # iterations: " << info[5] << '\n';
43  switch ((int) info[6]) {
44  case 1:
45  std::cerr << " stopped by small gradient J^T e\n";
46  break;
47  case 2:
48  std::cerr << " stopped by small Dp\n";
49  break;
50  case 3:
51  std::cerr << " stopped by itmax\n";
52  break;
53  case 4:
54  std::cerr << " singular matrix. Restart from current p with increased mu\n";
55  break;
56  case 5:
57  std::cerr << " no further error reduction is possible. Restart with increased mu\n";
58  break;
59  case 6:
60  std::cerr << " stopped by small ||e||_2\n";
61  break;
62  case 7:
63  std::cerr << " stopped by invalid (i.e. NaN or Inf) func values; a user error\n";
64  break;
65  }
66  std::cerr << " # function evaluations: " << info[7] << '\n';
67  std::cerr << " # Jacobian evaluations: " << info[8] << '\n';
68  std::cerr << " # linear systems solved: " << info[9] << "\n\n";
69 }
70 
71 }
72 
73 namespace ModelFitting {
74 
76 
77 static std::shared_ptr<LeastSquareEngine> createLevmarEngine(unsigned max_iterations) {
78  return std::make_shared<LevmarEngine>(max_iterations);
79 }
80 
82 
84  if (res == -1) {
86  }
87  switch (static_cast<int>(info[6])) {
88  case 7:
89  case 4:
91  case 3:
93  default:
95  }
96 }
97 
98 LevmarEngine::LevmarEngine(size_t itmax, double tau, double epsilon1,
99  double epsilon2, double epsilon3, double delta)
100  : m_itmax{itmax}, m_opts{tau, epsilon1, epsilon2, epsilon3, delta} {
101 #ifdef LINSOLVERS_RETAIN_MEMORY
102  logger.warn() << "Using a non thread safe levmar! Parallelism will be reduced.";
103 #endif
104 }
105 
106 LevmarEngine::~LevmarEngine() = default;
107 
108 
109 #ifdef LINSOLVERS_RETAIN_MEMORY
110 // If the Levmar library is not configured for multithreading, this mutex is used to ensure only one thread
111 // at a time can enter levmar
112 namespace {
113  std::mutex levmar_mutex;
114 }
115 #endif
116 
118  ResidualEstimator& residual_estimator) {
119  // Create a tuple which keeps the references to the given manager and estimator
120  auto adata = std::tie(parameter_manager, residual_estimator);
121 
122  // The function which is called by the levmar loop
123  auto levmar_res_func = [](double *p, double *hx, int, int, void *extra) {
124 #ifdef LINSOLVERS_RETAIN_MEMORY
125  levmar_mutex.unlock();
126 #endif
127  auto* extra_ptr = (decltype(adata)*)extra;
128  EngineParameterManager& pm = std::get<0>(*extra_ptr);
129  pm.updateEngineValues(p);
130  ResidualEstimator& re = std::get<1>(*extra_ptr);
131  re.populateResiduals(hx);
132 
133 #ifdef LINSOLVERS_RETAIN_MEMORY
134  levmar_mutex.lock();
135 #endif
136  };
137 
138  // Create the vector which will be used for keeping the parameter values
139  // and initialize it to the current values of the parameters
140  std::vector<double> param_values (parameter_manager.numberOfParameters());
141  parameter_manager.getEngineValues(param_values.begin());
142 
143  // Create a vector for getting the information of the minimization
145 
146  std::vector<double> covariance_matrix (parameter_manager.numberOfParameters() * parameter_manager.numberOfParameters());
147 
148 #ifdef LINSOLVERS_RETAIN_MEMORY
149  levmar_mutex.lock();
150 #endif
151  // Call the levmar library
152  auto res = dlevmar_dif(levmar_res_func, // The function called from the levmar algorithm
153  param_values.data(), // The pointer where the parameter values are
154  NULL, // We don't use any measurement vector
155  parameter_manager.numberOfParameters(), // The number of free parameters
156  residual_estimator.numberOfResiduals(), // The number of residuals
157  m_itmax, // The maximum number of iterations
158  m_opts.data(), // The minimization options
159  info.data(), // Where the information of the minimization is stored
160  NULL, // Working memory is allocated internally
161  covariance_matrix.data(),
162  &adata // No additional data needed
163  );
164 #ifdef LINSOLVERS_RETAIN_MEMORY
165  levmar_mutex.unlock();
166 #endif
167 
168  // Create and return the summary object
169  LeastSquareSummary summary {};
170 
171  auto converted_covariance_matrix = parameter_manager.convertCovarianceMatrixToWorldSpace(covariance_matrix);
172  for (unsigned int i=0; i<parameter_manager.numberOfParameters(); i++) {
173  summary.parameter_sigmas.push_back(sqrt(converted_covariance_matrix[i*(parameter_manager.numberOfParameters()+1)]));
174  }
175 
176  summary.status_flag = getStatusFlag(info, res);
177  summary.engine_stop_reason = info[6];
178  summary.iteration_no = info[5];
179  summary.underlying_framework_info = info;
180  return summary;
181 }
182 
183 } // end of namespace ModelFitting
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::EngineParameterManager::convertCovarianceMatrixToWorldSpace
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
Definition: EngineParameterManager.cpp:37
std::shared_ptr
STL class.
__attribute__
#define __attribute__(x)
ModelFitting::LeastSquareSummary::MAX_ITER
@ MAX_ITER
Definition: LeastSquareSummary.h:41
ModelFitting::logger
static Elements::Logging logger
Definition: LevmarEngine.cpp:75
Elements::Logging
std::vector< double >
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::LevmarEngine::LevmarEngine
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.
Definition: LevmarEngine.cpp:98
std::cerr
Exception.h
ModelFitting::LeastSquareSummary::ERROR
@ ERROR
Definition: LeastSquareSummary.h:41
LeastSquareEngineManager.h
ModelFitting::LevmarEngine::m_opts
std::vector< double > m_opts
Definition: LevmarEngine.h:72
ModelFitting::ResidualEstimator::populateResiduals
void populateResiduals(DoubleIter output_iter) const
std::array
STL class.
Elements::Logging::warn
void warn(const std::string &logMessage)
ModelFitting::EngineParameterManager
Class responsible for managing the parameters the least square engine minimizes.
Definition: EngineParameterManager.h:61
ModelFitting::levmar_engine
static LeastSquareEngineManager::StaticEngine levmar_engine
Definition: GSLEngine.cpp:38
ModelFitting::LeastSquareSummary::SUCCESS
@ SUCCESS
Definition: LeastSquareSummary.h:41
ModelFitting::createLevmarEngine
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
Definition: GSLEngine.cpp:34
ModelFitting::LevmarEngine::solveProblem
LeastSquareSummary solveProblem(EngineParameterManager &parameter_manager, ResidualEstimator &residual_estimator) override
Definition: LevmarEngine.cpp:117
Elements::Logging::getLogger
static Logging getLogger(const std::string &name="")
ModelFitting::EngineParameterManager::numberOfParameters
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.
Definition: EngineParameterManager.cpp:33
std::vector::begin
T begin(T... args)
LevmarEngine.h
std::mutex
STL class.
Logging.h
ModelFitting
Definition: AsinhChiSquareComparator.h:30
std::vector::data
T data(T... args)
ModelFitting::LevmarEngine::m_itmax
size_t m_itmax
Definition: LevmarEngine.h:71
ModelFitting::LeastSquareSummary::StatusFlag
StatusFlag
Definition: LeastSquareSummary.h:40
ModelFitting::LevmarEngine::~LevmarEngine
virtual ~LevmarEngine()
Destructor.