Mise à jour terminée. Pour connaître les apports de la version 13.8.4 par rapport à notre ancienne version vous pouvez lire les "Release Notes" suivantes :
https://about.gitlab.com/releases/2021/02/11/security-release-gitlab-13-8-4-released/
https://about.gitlab.com/releases/2021/02/05/gitlab-13-8-3-released/

Commit b83c9185 authored by Laurent Belcour's avatar Laurent Belcour

Towards multidimensional fitting

parent a7864faa
......@@ -5,6 +5,8 @@
#include <QtPlugin>
#include "common.h"
class data
{
public: // methods
......@@ -14,22 +16,23 @@ class data
virtual void load(const std::string& filename, double min, double max) = 0 ;
// Acces to data
virtual bool get(int i, double& x, double& y, double& t) const = 0 ;
virtual const std::vector<double>& get(int i) const = 0 ;
virtual const std::vector<double>& operator[](int i) const = 0 ;
// virtual bool get(int i, double& x, double& y, double& t) const = 0 ;
virtual const vec& get(int i) const = 0 ;
virtual const vec& operator[](int i) const = 0 ;
// Get data size, e.g. the number of samples to fit
virtual int size() const = 0 ;
// Get the dimension of the input and output space
// WARNING: this dimension is defined after loading
// the data!
virtual int input_dimension() const = 0;
virtual int output_dimension() const = 0;
// Get min and max input space values
virtual double min() const = 0 ;
virtual double max() const = 0 ;
virtual vec min() const = 0 ;
virtual vec max() const = 0 ;
virtual int dimX() const { return _nX ; } ;
virtual int dimY() const { return _nY ; } ;
protected:
// Dimension of the function
int _nX, _nY ;
} ;
Q_DECLARE_INTERFACE(data, "Fitter.Data")
......@@ -3,6 +3,7 @@
#include "function.h"
#include "data.h"
#include "args.h"
#include "common.h"
#include <QtPlugin>
......
......@@ -5,17 +5,28 @@
#include <QtPlugin>
class function : public std::function<double(double)>
#include "common.h"
class function //: public std::function<double(double)>
{
public: // methods
// Overload the function operator
virtual double operator()(double x) const = 0 ;
virtual vec operator()(const vec& x) const = 0 ;
// IO function to text files
virtual void load(const std::string& filename) = 0 ;
virtual void save() const = 0 ;
virtual int dimX() const { return _nX ; } ;
virtual int dimY() const { return _nY ; } ;
virtual int setDimX(int nX) { _nX = nX ; } ;
virtual int setDimY(int nY) { _nY = nY ; } ;
protected:
// Dimension of the function
int _nX, _nY ;
} ;
Q_DECLARE_INTERFACE(function, "Fitter.Function")
......@@ -17,16 +17,16 @@ void rational_data::load(const std::string& filename)
void rational_data::load(const std::string& filename, double min, double max)
{
std::ifstream file(filename) ;
_min = std::numeric_limits<double>::max() ;
_max = -std::numeric_limits<double>::max() ;
if(!file.is_open())
{
std::cerr << "<<ERROR>> unable to open file \"" << filename << "\"" << std::endl ;
}
// N-Floats regexp
boost::regex e ("^([0-9]*\.?[0-9]+[\\t ]?)+");
// boost::regex e ("^([0-9]*\.?[0-9]+[\\t ]?)");
_nX = 0 ; _nY = 0 ;
std::vector<int> vs ;
double x, y, dy ;
while(file.good())
......@@ -36,30 +36,78 @@ void rational_data::load(const std::string& filename, double min, double max)
std::stringstream linestream(line) ;
// Discard incorrect lines
if(!boost::regex_match(line,e))
if(linestream.peek() == '#')
{
linestream.ignore(1) ;
std::string comment ;
linestream >> comment ;
std::cout << comment << std::endl ;
if(comment == std::string("DIM"))
{
linestream >> _nX >> _nY ;
vs.resize(dimY()) ;
for(int k=0; k<dimY(); ++k)
{
vs[k] = 0 ;
}
_min.resize(dimX()) ;
_max.resize(dimX()) ;
for(int k=0; k<dimY(); ++k)
{
_min[k] = std::numeric_limits<double>::max() ;
_max[k] = -std::numeric_limits<double>::max() ;
}
}
continue ;
}
else if(line.empty()/*!boost::regex_match(line,e)*/)
{
continue ;
}
linestream >> x >> y ;
if(linestream.good()) {
linestream >> dy ;
} else {
// TODO Specify the delta in case
dy = 0.01f ;
vec v ;
v.resize(dimX() + 3*dimY()) ;
for(int i=0; i<dimX(); ++i)
linestream >> v[i] ;
for(int i=0; i<dimY(); ++i)
linestream >> v[dimX() + i] ;
for(int i=0; i<dimY(); ++i)
{
// TODO, the firts case does not account for the
// dimension of the ouput vector
if(linestream.good())
{
linestream >> v[dimX() + dimY()+i] ;
}
else
{
// TODO Specify the delta in case
// Handle multiple dim
v[dimX() + dimY()+i] = v[dimX() + i] - 0.01f ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] + 0.01f ;
}
}
// If data is not in the interval of fit
// TODO: Update to more dims
if(v[0] < min || v[0] > max)
{
continue ;
}
if(x <= max && x >= min)
_data.push_back(v) ;
// Update min and max
for(int k=0; k<dimX(); ++k)
{
std::vector<double> v ;
v.push_back(x) ;
v.push_back(y-dy) ;
v.push_back(y+dy) ;
_data.push_back(v) ;
// Update min and max
_min = std::min(_min, x) ;
_max = std::max(_max, x) ;
_min[k] = std::min(_min[k], v[k]) ;
_max[k] = std::max(_max[k], v[k]) ;
}
}
......@@ -75,10 +123,11 @@ void rational_data::load(const std::string& filename, double min, double max)
_data = temp ;
*/
// Sort the vector
std::sort(_data.begin(), _data.end(), [](const std::vector<double>& a, const std::vector<double>& b){return (a[0]<b[0]);});
// std::sort(_data.begin(), _data.end(), [](const std::vector<double>& a, const std::vector<double>& b){return (a[0]<b[0]);});
std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ;
std::cout << "<<INFO>> data inside [" << _min << ", " << _max << "]" << std::endl ;
std::cout << "<<INFO>> data inside " << _min << " ... " << _max << std::endl ;
std::cout << "<<INFO>> loading data file of R^" << dimX() << " -> R^" << dimY() << std::endl ;
std::cout << "<<INFO>> " << _data.size() << " elements to fit" << std::endl ;
}
......@@ -96,20 +145,21 @@ bool rational_data::get(int i, double& x, double& yl, double& yu) const
return true ;
}
int rational_data::input_dimension() const
{
return 1 ;
}
int rational_data::output_dimension() const
void rational_data::get(int i, vec& yl, vec& yu) const
{
return 1 ;
yl.resize(dimY()) ; yu.resize(dimY()) ;
for(int j=0; j<dimY(); ++j)
{
yl[j] = _data[i][dimX() + dimY() + j] ;
yu[j] = _data[i][dimX() + 2*dimY() + j] ;
}
}
const std::vector<double>& rational_data::operator[](int i) const
const vec& rational_data::operator[](int i) const
{
return _data[i] ;
}
const std::vector<double>& rational_data::get(int i) const
const vec& rational_data::get(int i) const
{
return _data[i] ;
}
......@@ -119,12 +169,12 @@ int rational_data::size() const
return _data.size() ;
}
double rational_data::min() const
vec rational_data::min() const
{
return _min ;
}
double rational_data::max() const
vec rational_data::max() const
{
return _max ;
}
......
......@@ -26,30 +26,25 @@ class rational_data : public QObject, public data
// Acces to data
virtual bool get(int i, double& x, double& yl, double& yu) const ;
virtual const std::vector<double>& get(int i) const ;
virtual const std::vector<double>& operator[](int i) const ;
virtual const vec& get(int i) const ;
virtual void get(int i, vec& yl, vec& yu) const ;
virtual const vec& operator[](int i) const ;
// Get data size
virtual int size() const ;
// Get min and max input parameters
virtual double min() const ;
virtual double max() const ;
// Get the dimension of the input and output space
// WARNING: this dimension is defined after loading
// the data!
virtual int input_dimension() const ;
virtual int output_dimension() const ;
virtual vec min() const ;
virtual vec max() const ;
private: // data
// Store for each point of data, the upper
// and lower value
std::vector<std::vector<double> > _data ;
std::vector<vec> _data ;
// Store the min and max value on the input
// domain
double _min, _max ;
vec _min, _max ;
} ;
......@@ -10,4 +10,4 @@ SOURCES = rational_data.cpp
LIBS += -lboost_regex
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC -g
......@@ -64,18 +64,25 @@ void rational_fitter_cgal::set_parameters(const arguments& args)
}
bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf)
bool rational_fitter_cgal::fit_data(const data* dat, int np, int nq, function*& rf)
{
// by default, we have a nonnegative QP with Ax - b >= 0
Program qp (CGAL::LARGER, false, 0, false, 0) ;
rational_function* r = dynamic_cast<rational_function*>(rf) ;
if(r == nullptr)
const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
if(r == nullptr || d == nullptr)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<np+nq; ++i)
......@@ -103,7 +110,7 @@ bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf
// Filling the p part
if(j<np)
{
const double pi = r->p((*d)[i][0], j) ;
const double pi = r->p(d->get(i), j) ;
a0_norm += pi*pi ;
a1_norm += pi*pi ;
qp.set_a(j, i, ET(pi)) ;
......@@ -114,14 +121,17 @@ bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf
// Filling the q part
else
{
const double qi = r->q((*d)[i][0], j-np) ;
a0_norm += qi*qi * ((*d)[i][1]*(*d)[i][1]) ;
qp.set_a(j, i, ET(-(*d)[i][1] * qi)) ;
CI(j, i) = -(*d)[i][1] * qi ;
vec yl, yu ;
d->get(i, yl, yu) ;
const double qi = r->q(d->get(i), j-np) ;
a0_norm += qi*qi * (yl[0]*yl[0]) ;
qp.set_a(j, i, ET(-yl[0] * qi)) ;
CI(j, i) = -yl[0] * qi ;
a1_norm += qi*qi * ((*d)[i][2]*(*d)[i][2]) ;
qp.set_a(j, i+d->size(), ET((*d)[i][2] * qi)) ;
CI(j, i+d->size()) = (*d)[i][2] * qi ;
a1_norm += qi*qi * (yu[0]*yu[0]) ;
qp.set_a(j, i+d->size(), ET(yu[0] * qi)) ;
CI(j, i+d->size()) = yu[0] * qi ;
}
}
......@@ -146,7 +156,7 @@ bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf
std::cout << " ]" << std::endl ;
#endif
const double delta = sigma_m / sigma_M ;
double delta = sigma_m / sigma_M ;
if(isnan(delta) || (abs(delta) == std::numeric_limits<double>::infinity()))
{
#ifdef DEBUG
......@@ -154,6 +164,10 @@ bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf
#endif
return false ;
}
else if(delta == 0.0)
{
delta = 1.0 ;
}
#ifdef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
......@@ -253,6 +267,8 @@ bool rational_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf
delete r ;
}
r = new rational_function(p, q);
std::cout << "<<INFO>> got solution " << *r << std::endl ;
return true;
}
else
......
......@@ -13,5 +13,5 @@ LIBS += -lCGAL -lboost_regex \
-L../build \
-lrational_function -lrational_data
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC -g
......@@ -23,32 +23,86 @@ rational_function::~rational_function()
}
// Overload the function operator
double rational_function::operator()(double x) const
vec rational_function::operator()(const vec& x) const
{
double p = 0.0f ;
double q = 0.0f ;
vec res ;
res.resize(dimY()) ;
for(int i=a.size()-1; i>=0; --i)
for(int k=0; k<_nY; ++k)
{
p = x*p + a[i] ;
}
double p = 0.0f ;
double q = 0.0f ;
for(int l=0; l<_nX; ++l)
{
for(int i=b.size()-1; i>=0; --i)
{
q = x*q + b[i] ;
}
for(int i=a.size()-1; i>=0; --i)
{
p = x[l]*p + a[i] ;
}
return p/q ;
for(int i=b.size()-1; i>=0; --i)
{
q = x[l]*q + b[i] ;
}
}
res[k] = p/q ;
}
return res ;
}
// Get the p_i and q_j function
double rational_function::p(double x, int i) const
double rational_function::p(const vec& x, int i) const
{
return pow(x, i) ;
std::vector<int> deg ; deg.assign(dimY(), 0) ;
double res = 1.0 ;
if(i == 0)
return res ;
int temp_i = i ;
int temp_c ;
while(temp_i != 0)
{
temp_c = (temp_i-1) % dimX() ;
temp_i = (temp_i - temp_c) / dimX() ;
deg[temp_c] += 1 ;
}
for(int k=0; k<dimX(); ++k)
{
res *= pow(x[k], deg[k]) ;
}
return res ;
}
double rational_function::q(double x, int j) const
double rational_function::q(const vec& x, int i) const
{
return pow(x, j) ;
std::vector<int> deg ; deg.assign(dimY(), 0) ;
double res = 1.0 ;
if(i == 0)
return res ;
int temp_i = i ;
int temp_c ;
while(temp_i != 0)
{
temp_c = (temp_i-1) % dimX() ;
temp_i = (temp_i - temp_c) / dimX() ;
deg[temp_c] += 1 ;
}
for(int k=0; k<dimX(); ++k)
{
res *= pow(x[k], deg[k]) ;
}
return res ;
}
// IO function to text files
......
......@@ -12,6 +12,7 @@
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
#include <core/common.h>
class rational_function : public QObject, public function
{
......@@ -25,11 +26,11 @@ class rational_function : public QObject, public function
virtual ~rational_function() ;
// Overload the function operator
virtual double operator()(double x) const ;
virtual vec operator()(const vec& x) const ;
// Get the p_i and q_j function
virtual double p(double x, int i) const ;
virtual double q(double x, int j) const ;
virtual double p(const vec& x, int i) const ;
virtual double q(const vec& x, int j) const ;
// IO function to text files
void load(const std::string& filename) ;
......
......@@ -10,5 +10,5 @@ SOURCES = rational_function.cpp
LIBS += -lboost_regex
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC -rdynamic
QMAKE_CXXFLAGS += -std=c++11 -frounding-math -fPIC -rdynamic -g
......@@ -17,14 +17,23 @@
int main(int argc, char** argv)
{
QApplication app(argc, argv, false);
// QCoreApplication::addLibraryPath() ;
QCoreApplication::addLibraryPath(QString("/home/belcour/Projects/alta/sources/tests/plugin_loader/")) ;
QApplication app(argc, argv, false);
arguments args(argc, argv) ;
std::vector<function*> functions ;
std::vector<data*> datas ;
std::vector<fitter*> fitters ;
QDir pluginsDir = QDir(app.applicationDirPath());
if(args.is_defined("plugins"))
{
pluginsDir = QDir(args["plugins"].c_str()) ;
}
foreach (QString fileName, pluginsDir.entryList(QDir::Files))
{
QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
......@@ -62,7 +71,6 @@ int main(int argc, char** argv)
}
}
arguments args(argc, argv) ;
if(! args.is_defined("input")) {
std::cerr << "<<ERROR>> the input filename is not defined" << std::endl ;
return 1 ;
......@@ -100,16 +108,13 @@ int main(int argc, char** argv)
// Display the result
if(is_fitted)
{
// std::cout << *f << std::endl ;
//*
std::ofstream file(args["output"], std::ios_base::trunc);
const double dt = (d->max() - d->min()) / 100.0f ;
for(double x=d->min(); x<=d->max(); x+=dt)
const double dt = (d->max()[0] - d->min()[0]) / 100.0f ;
for(double x=d->min()[0]; x<=d->max()[0]; x+=dt)
{
file << x << "\t" << (*f)(x) << std::endl ;
vec vx ; vx.push_back(x) ;
file << x << "\t" << ((*f)(vx))[0] << std::endl ;
}
//*/
}