Commit fc6add90 authored by Laurent Belcour's avatar Laurent Belcour

Fixing the quadratic program to support streaming of constraints. It works now

but the streaming algorithm is dumb.
parent bdf57e7e
...@@ -9,195 +9,233 @@ ...@@ -9,195 +9,233 @@
class quadratic_program class quadratic_program
{ {
public: public:
//! \brief Constructor need to specify the number of coefficients //! \brief Constructor need to specify the number of coefficients
quadratic_program(int np, int nq) : _np(np), _nq(nq), CI(0.0, _np+_nq, 0) { } quadratic_program(int np, int nq) : _np(np), _nq(nq), CI(0.0, _np+_nq, 0) { }
//! \brief Remove the already defined constraints //! \brief Remove the already defined constraints
void clear_constraints() void clear_constraints()
{ {
CI.resize(_np+_nq, 0); CI.resize(_np+_nq, 0);
} }
//! \brief Add a constraint by specifying the vector //! \brief Add a constraint by specifying the vector
void add_constraints(const vec& c) void add_constraints(const vec& c)
{ {
const int m = CI.nrows(); const int m = CI.nrows();
const int n = CI.ncols(); const int n = CI.ncols();
if(n > 0) if(n > 0)
{ {
// Construct temp buffer // Construct temp buffer
double* temp = new double[n*m]; double* temp = new double[n*m];
for(int u=0; u<n; ++u) for(int u=0; u<n; ++u)
{ {
for(int v=0; v<m; ++v) for(int v=0; v<m; ++v)
{ {
temp[u*m + v] = CI[v][u]; temp[u*m + v] = CI[v][u];
} }
} }
// Resize matrix CI // Resize matrix CI
CI.resize(m, n+1); CI.resize(m, n+1);
// Recopy data // Recopy data
for(int u=0; u<n+1; ++u) for(int u=0; u<n+1; ++u)
{ {
if(u==n) if(u==n)
{ {
for(int v=0; v<m; ++v) for(int v=0; v<m; ++v)
CI[v][u] = c[v]; CI[v][u] = c[v];
} }
else else
{ {
for(int v=0; v<m; ++v) for(int v=0; v<m; ++v)
CI[v][u] = temp[u*m + v]; CI[v][u] = temp[u*m + v];
} }
} }
delete[] temp; delete[] temp;
} }
else else
{ {
// Resize matrix CI // Resize matrix CI
CI.resize(m, 1); CI.resize(m, 1);
// Recopy data // Recopy data
for(int u=0; u<m; ++u) for(int u=0; u<m; ++u)
CI[n][u] = c[u]; CI[n][u] = c[u];
} }
} }
//! \brief Provide the number of constraints //! \brief Provide the number of constraints
int nb_constraints() const int nb_constraints() const
{ {
return CI.ncols(); return CI.ncols();
} }
//! \brief Solves the quadratic program and update the p and //! \brief Solves the quadratic program and update the p and
//! q vector if necessary. //! q vector if necessary.
inline bool solve_program(QuadProgPP::Vector<double>& x, double& delta, vec& p, vec& q) inline bool solve_program(QuadProgPP::Vector<double>& x, double& delta, vec& p, vec& q)
{ {
bool solves_qp = solve_program(x, delta) ; bool solves_qp = solve_program(x, delta) ;
if(solves_qp) if(solves_qp)
{ {
double norm = 0.0; double norm = 0.0;
for(int i=0; i<_np+_nq; ++i) for(int i=0; i<_np+_nq; ++i)
{ {
const double v = x[i]; const double v = x[i];
norm += v*v ; norm += v*v ;
if(i < _np) if(i < _np)
{ {
p[i] = v ; p[i] = v ;
} }
else else
{ {
q[i-_np] = v ; q[i-_np] = v ;
} }
} }
return norm > 0.0; return norm > 0.0;
} }
else else
{ {
return false ; return false ;
} }
} }
//! \brief Solves the quadratic program //! \brief Solves the quadratic program
inline bool solve_program(QuadProgPP::Vector<double>& v, double& delta) inline bool solve_program(QuadProgPP::Vector<double>& v, double& delta)
{ {
const int m = CI.nrows(); const int m = CI.nrows();
const int n = CI.ncols(); const int n = CI.ncols();
QuadProgPP::Matrix<double> G (0.0, m, m) ; QuadProgPP::Matrix<double> G (0.0, m, m) ;
QuadProgPP::Vector<double> g (0.0, m) ; QuadProgPP::Vector<double> g (0.0, m) ;
QuadProgPP::Vector<double> ci(0.0, n) ; QuadProgPP::Vector<double> ci(0.0, n) ;
QuadProgPP::Matrix<double> CE(0.0, m, 0) ; QuadProgPP::Matrix<double> CE(0.0, m, 0) ;
QuadProgPP::Vector<double> ce(0.0, 0) ; QuadProgPP::Vector<double> ce(0.0, 0) ;
// Update the ci column with the delta parameter // Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12) // (See Celis et al. 2007 p.12)
Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(Eigen::MatrixXd::Map(&CI[0][0], m, n)); Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(Eigen::MatrixXd::Map(&CI[0][0], m, n));
const double sigma_m = svd.singularValues()(std::min(m, n)-1) ; const double sigma_m = svd.singularValues()(std::min(m, n)-1) ;
const double sigma_M = svd.singularValues()(0) ; const double sigma_M = svd.singularValues()(0) ;
delta = sigma_M / sigma_m ; delta = sigma_M / sigma_m ;
// Select the size of the result vector to // Select the size of the result vector to
// be equal to the dimension of p + q // be equal to the dimension of p + q
for(int i=0; i<m; ++i) for(int i=0; i<m; ++i)
{ {
G[i][i] = 1.0 ; G[i][i] = 1.0 ;
} }
// Each constraint (fitting interval or point // Each constraint (fitting interval or point
// add another dimension to the constraint // add another dimension to the constraint
// matrix // matrix
for(int i=0; i<n; ++i) for(int i=0; i<n; ++i)
{ {
// Norm of the row vector // Norm of the row vector
double norm = 0.0 ; double norm = 0.0 ;
for(int j=0; j<m; ++j) for(int j=0; j<m; ++j)
{ {
norm += CI[j][i]*CI[j][i] ; ; norm += CI[j][i]*CI[j][i] ; ;
} }
// Set the c vector, will later be updated using the // Set the c vector, will later be updated using the
// delta parameter. // delta parameter.
ci[i] = -delta * sqrt(norm) ; ci[i] = -delta * sqrt(norm) ;
} }
// Compute the solution // Compute the solution
const double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, v); const double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, v);
bool solves_qp = !(cost == std::numeric_limits<double>::infinity()); bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
return solves_qp; return solves_qp;
} }
//! \brief Test all the constraints of the data //! \brief Test all the constraints of the data
bool test_constraints(const rational_function* r, const vertical_segment* data) bool test_constraints(int ny, const rational_function* r, const vertical_segment* data)
{ {
int nb_failed = 0; int nb_failed = 0;
for(int n=0; n<data->size(); ++n) for(int n=0; n<data->size(); ++n)
{ {
vec x, yl, yu; vec x, yl, yu;
data->get(n, x, yl, yu); data->get(n, x, yl, yu);
vec y = r->value(x); vec y = r->value(x);
if(y < yl || y > yu) bool fail_upper = y > yu;
{ bool fail_lower = y < yl;
nb_failed++; if(fail_lower || fail_upper)
} {
} nb_failed++;
vec cu, cl;
get_constraint(x, yl, yu, ny, r, cu, cl);
add_constraints(cu);
add_constraints(cl);
}
}
#ifdef DEBUG #ifdef DEBUG
std::cout << "<<TRACE>> " << nb_failed << " constraints where not satified." << std::endl; std::cout << "<<TRACE>> " << nb_failed << " constraints where not satified." << std::endl;
#endif #endif
return nb_failed == 0;
return nb_failed == 0; }
}
//! \brief Generate two constraint vectors from a vertical segment and a
//! \brief Give the next position in the data that is not satisfied. //! ration function type.
//! This method works only for a single color channel ny ! inline void get_constraint(const vec& xi, const vec& yl, const vec& yu, int ny,
static int next_unmatching_constraint(int i, int ny, const rational_function* r, const rational_function* func,
const vertical_segment* data) vec& cu, vec& cl)
{ {
for(int n=i; n<data->size(); ++n) cu.resize(_np+_nq);
{ cl.resize(_np+_nq);
vec x, yl, yu;
data->get(n, x, yl, yu); // Create two vector of constraints
for(int j=0; j<_np+_nq; ++j)
vec y = r->value(x); {
if(y[ny] < yl[ny] || y[ny] > yu[ny]) // Filling the p part
{ if(j<_np)
{
const double pi = func->p(xi, j) ;
cu[j] = pi ;
cl[j] = -pi ;
}
// Filling the q part
else
{
const double qi = func->q(xi, j-_np) ;
cu[j] = -yu[ny] * qi ;
cl[j] = yl[ny] * qi ;
}
}
}
//! \brief Give the next position in the data that is not satisfied.
//! This method works only for a single color channel ny !
static int next_unmatching_constraint(int i, int ny, const rational_function* r,
const vertical_segment* data)
{
for(int n=i; n<data->size(); ++n)
{
vec x, yl, yu;
data->get(n, x, yl, yu);
vec y = r->value(x);
if(y[ny] < yl[ny] || y[ny] > yu[ny])
{
return n; return n;
} }
} }
return data->size(); return data->size();
} }
protected: protected:
int _np, _nq ; int _np, _nq ;
QuadProgPP::Matrix<double> CI; QuadProgPP::Matrix<double> CI;
}; };
...@@ -52,13 +52,14 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar ...@@ -52,13 +52,14 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
r->setMin(d->min()) ; r->setMin(d->min()) ;
r->setMax(d->max()) ; r->setMax(d->max()) ;
const int _min_np = args.get_int("min-np", 10); const int _min_np = args.get_int("min-np", 10);
const int _max_np = args.get_int("np", _min_np); const int _max_np = args.get_int("np", _min_np);
std::cout << "<<INFO>> N in [" << _min_np << ", " << _max_np << "]" << std::endl ; std::cout << "<<INFO>> N in [" << _min_np << ", " << _max_np << "]" << std::endl ;
for(int i=_min_np; i<=_max_np; ++i) for(int i=_min_np; i<=_max_np; ++i)
{ {
std::cout << "<<INFO>> fit using np+nq = " << i << "\r" ; std::cout << "<<INFO>> fit using np+nq = " << i << std::endl ;
std::cout.flush() ; std::cout.flush() ;
QTime time ; QTime time ;
time.start() ; time.start() ;
...@@ -82,31 +83,32 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar ...@@ -82,31 +83,32 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
#ifdef DEBUG #ifdef DEBUG
std::cout << "<<DEBUG>> will use " << nb_cores << " threads to compute the quadratic programs" << std::endl ; std::cout << "<<DEBUG>> will use " << nb_cores << " threads to compute the quadratic programs" << std::endl ;
#endif #endif
omp_set_num_threads(nb_cores) ; omp_set_num_threads(nb_cores) ;
std::vector<rational_function*> rs; std::vector<rational_function*> rs;
for(int j=0; j<nb_cores; ++j) for(int j=0; j<nb_cores; ++j)
{ {
rational_function* rj = dynamic_cast<rational_function*>(plugins_manager::get_function(args["func"])); rational_function* rj = dynamic_cast<rational_function*>(plugins_manager::get_function(args["func"]));
rj->setDimX(d->dimX()) ; rj->setDimX(d->dimX()) ;
rj->setDimY(d->dimY()) ; rj->setDimY(1) ;
rj->setMin(d->min()) ; rj->setMin(d->min()) ;
rj->setMax(d->max()) ; rj->setMax(d->max()) ;
if(rj == NULL) if(rj == NULL)
{ {
std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl; std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
return false; return false;
} }
rs.push_back(rj); rs.push_back(rj);
} }
double min_delta = std::numeric_limits<double>::max(); double min_delta = std::numeric_limits<double>::max();
int nb_sol_found = 0; int nb_sol_found = 0;
int np, nq ; int np, nq ;
#pragma omp parallel for // #pragma omp parallel for
for(int j=1; j<i; ++j) for(int j=1; j<i; ++j)
{ {
int temp_np = i - j; int temp_np = i - j;
...@@ -115,10 +117,10 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar ...@@ -115,10 +117,10 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
vec p(temp_np*r->dimY()), q(temp_nq*r->dimY()); vec p(temp_np*r->dimY()), q(temp_nq*r->dimY());
double delta; double delta;
bool is_fitted = fit_data(d, temp_np, temp_nq, rs[omp_get_thread_num()], p, q, delta); bool is_fitted = fit_data(d, temp_np, temp_nq, rs[omp_get_thread_num()], p, q, delta);
if(is_fitted) if(is_fitted)
{ {
#pragma omp critical #pragma omp critical
{ {
++nb_sol_found ; ++nb_sol_found ;
if(delta < min_delta) if(delta < min_delta)
...@@ -136,11 +138,11 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar ...@@ -136,11 +138,11 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
} }
} }
// Clean memory // Clean memory
for(int j=0; j<nb_cores; ++j) for(int j=0; j<nb_cores; ++j)
{ {
delete rs[j]; delete rs[j];
} }
if(min_delta < std::numeric_limits<double>::max()) if(min_delta < std::numeric_limits<double>::max())
{ {
...@@ -165,8 +167,8 @@ void rational_fitter_parallel::set_parameters(const arguments& args) ...@@ -165,8 +167,8 @@ void rational_fitter_parallel::set_parameters(const arguments& args)
bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq,
rational_function* r, rational_function* r,
vec& P, vec& Q, double& delta) vec& P, vec& Q, double& delta)
{ {
for(int j=0; j<d->dimY(); ++j) for(int j=0; j<d->dimY(); ++j)
{ {
...@@ -186,71 +188,81 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n ...@@ -186,71 +188,81 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals) // y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
// the function return a ration BRDF function and a boolean // the function return a ration BRDF function and a boolean
bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, int ny, bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, int ny,
rational_function* r, rational_function* r,
vec& p, vec& q, double& delta) vec& p, vec& q, double& delta)
{ {
const int m = d->size(); // 2*m = number of constraints const int m = d->size(); // 2*m = number of constraints
const int n = np+nq; // n = np+nq const int n = np+nq; // n = np+nq
quadratic_program qp(np, nq); quadratic_program qp(np, nq);
#ifndef TODO_PUT_IN_METHOD #ifndef TODO_PUT_IN_METHOD
for(int i=0; i<d->size()/10; ++i) for(int i=0; i<d->size()/100; ++i)
{ {
// Create two vector of constraints // Create two vector of constraints
vec c1(n), c2(n); vec c1(n), c2(n);
get_constraint(10*i, np, nq, ny, d, r, c1, c2); get_constraint(100*i, np, nq, ny, d, r, c1, c2);
qp.add_constraints(c1); qp.add_constraints(c1);
qp.add_constraints(c2); qp.add_constraints(c2);
} }
#endif #endif
while(true) while(qp.nb_constraints() < 2*m)
{ {
#ifdef DEBUG
std::cout << "<<DEBUG>> number of constraints = " << qp.nb_constraints() << std::endl ;
#endif
QuadProgPP::Vector<double> x(n);