SourceXtractorPlusPlus  0.11
Please provide a description of the project.
Example_SimpleFit.cpp
Go to the documentation of this file.
1 
23 #include <iostream>
24 #include <vector>
25 #include <chrono>
31 #include "utils.h"
32 
33 using namespace ModelFitting;
35 
37 private:
40 
41 public:
44  : m_A{a}, m_lambda{lambda}, m_b{b}, m_y{y}, m_t{t} {
45  }
46 
47  virtual ~ExpResidualProvider() = default;
48 
49  std::size_t numberOfResiduals() const override {
50  return m_y.size();
51  }
52 
53  void populateResidualBlock(IterType iter) override {
54  std::cout << std::fixed << std::setprecision(4) << "A=" << m_A->getValue() << "\tlambda=" << m_lambda->getValue() << "\tb=" << m_b->getValue();
55  double e = 0;
56  for (size_t i = 0; i < m_y.size(); ++i) {
57  double Yi = m_b->getValue() + m_A->getValue() * exp(-m_lambda->getValue() * m_t[i]);
58  *iter = Yi - m_y[i];
59  e += *iter * *iter;
60  ++iter;
61  }
62  std::cout << std::scientific << "\tsum(e^2)=" << e << std::endl;
63  }
64 };
65 
66 /*
67  * Adapted from
68  * https://www.gnu.org/software/gsl/doc/html/nls.html#exponential-fitting-example
69  * "Simple" exponential function, so it is easier to understand and visualize what's going on
70  */
71 int main(int argc, char **argv) {
72  std::string engine_impl("levmar");
73  if (argc > 1) {
74  engine_impl = argv[1];
75  }
76 
77  // Generate data points, using the exponential function
78  int n = 100;
79  double tmax = 40.;
80  std::vector<double> y(n), t(n);
81 
82  for (int i = 0; i < n; ++i) {
83  double ti = i * tmax / (n - 1.0);
84  double yi = 1. + 5. * exp(-0.1 * ti);
85  t[i] = ti;
86  y[i] = yi;
87  }
88 
89  // Initialize and run engine
90  auto A = std::make_shared<EngineParameter>(2., make_unique<NeutralConverter>());
91  auto lambda = std::make_shared<EngineParameter>(1., make_unique<NeutralConverter>());
92  auto b = std::make_shared<EngineParameter>(0., make_unique<NeutralConverter>());
93 
94  EngineParameterManager manager {};
95  manager.registerParameter(A);
96  manager.registerParameter(lambda);
97  manager.registerParameter(b);
98 
99  ResidualEstimator res_estimator {};
100  res_estimator.registerBlockProvider(make_unique<ExpResidualProvider>(A, lambda, b, y, t));
101 
102  std::cout << "Using engine " << engine_impl << std::endl;
103  auto engine = LeastSquareEngineManager::create(engine_impl);
104  auto t1 = std::chrono::steady_clock::now();
105  auto solution = engine->solveProblem(manager, res_estimator);
106  auto t2 = std::chrono::steady_clock::now();
107 
109  std::cout << "Time of fitting: " << std::chrono::duration<double, std::milli>(t2 - t1).count() << " ms" << std::endl;
110  std::cout << "A (5.) : " << A->getValue() << std::endl;
111  std::cout << "lambda (0.1): " << lambda->getValue() << std::endl;
112  std::cout << "b (1.) : " << b->getValue() << std::endl;
113  printLevmarInfo(boost::any_cast<std::array<double, 10>>(solution.underlying_framework_info));
114 }
ExpResidualProvider::populateResidualBlock
void populateResidualBlock(IterType iter) override
Provides the residual values.
Definition: Example_SimpleFit.cpp:53
ModelFitting::ResidualEstimator
Provides to the LeastSquareEngine the residual values.
Definition: ResidualEstimator.h:50
std::setprecision
T setprecision(T... args)
ModelFitting::LeastSquareEngineManager::create
static std::shared_ptr< LeastSquareEngine > create(const std::string &name, unsigned max_iterations=1000)
Definition: LeastSquareEngineManager.cpp:65
std::string
STL class.
std::shared_ptr< BasicParameter >
ModelFitting::ResidualBlockProvider
Interface of a class which can provide a block of residuals for least square minimization solving.
Definition: ResidualBlockProvider.h:35
ExpResidualProvider::ExpResidualProvider
ExpResidualProvider(std::shared_ptr< BasicParameter > a, std::shared_ptr< BasicParameter > lambda, std::shared_ptr< BasicParameter > b, const std::vector< double > &y, const std::vector< double > &t)
Definition: Example_SimpleFit.cpp:42
std::vector< double >
std::vector::size
T size(T... args)
std::chrono::duration
ExpResidualProvider::numberOfResiduals
std::size_t numberOfResiduals() const override
Returns the number of residuals provided by this provider.
Definition: Example_SimpleFit.cpp:49
ModelFitting::BasicParameter::getValue
virtual double getValue() const
Definition: BasicParameter.h:62
ModelFitting::ResidualBlockProvider::IterType
double * IterType
Definition: ResidualBlockProvider.h:45
ExpResidualProvider::m_lambda
std::shared_ptr< BasicParameter > m_lambda
Definition: Example_SimpleFit.cpp:38
utils.h
ExpResidualProvider::~ExpResidualProvider
virtual ~ExpResidualProvider()=default
NeutralConverter.h
NormalizedConverter.h
std::cout
LeastSquareEngineManager.h
ModelFitting::ResidualEstimator::registerBlockProvider
void registerBlockProvider(std::unique_ptr< ResidualBlockProvider > provider)
Registers a ResidualBlockProvider to the ResidualEstimator.
Definition: ResidualEstimator.cpp:29
std::array
STL class.
ModelFitting::EngineParameterManager
Class responsible for managing the parameters the least square engine minimizes.
Definition: EngineParameterManager.h:61
ModelFitting::EngineParameterManager::registerParameter
void registerParameter(std::shared_ptr< EngineParameter > parameter)
Registers an EngineParameter to the EngineParameterManager.
Definition: EngineParameterManager.cpp:29
Euclid::make_unique
std::unique_ptr< T > make_unique(Args &&... args)
e
constexpr double e
main
int main(int argc, char **argv)
Definition: Example_SimpleFit.cpp:71
std::endl
T endl(T... args)
std::exp
T exp(T... args)
printLevmarInfo
void printLevmarInfo(std::array< double, 10 > info)
Definition: utils.h:118
std::chrono::duration::count
T count(T... args)
std::fixed
T fixed(T... args)
std::size_t
y
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
Definition: MoffatModelFittingTask.cpp:93
ExpResidualProvider::m_y
std::vector< double > m_y
Definition: Example_SimpleFit.cpp:39
ModelFitting
Definition: AsinhChiSquareComparator.h:30
memory_tools.h
ExpResidualProvider
Definition: Example_SimpleFit.cpp:36
EngineParameterManager.h
std::chrono::steady_clock::now
T now(T... args)