Commit 1d1303b2 authored by Laurent Belcour's avatar Laurent Belcour

Merge

parents 3eb82600 745bf8aa
......@@ -29,10 +29,11 @@ directory.
Dependencies:
ALTA core: Eigen
Plugin rational_eigen: Eigen 3.x
Plugin rational_quadprog: Quadprog++
Plugin rational_quadprog: Quadprog++
Plugin rational_cgal: The CGAL library
Plugin rational_parallel: The OpenMP library, Quadprog++ library and Eigen
Plugin rational_matlab: Matlab engine (matlab.prf required)
Plugin rational_parsec_*: PLASMA coreblas, PaRSEC runtime
Plugin nonlinear_eigen: Eigen
Plugin nonlinear_ceres: CERES library and its dependencies
Plugin nonlinear_nlopt: NLOpt library and its dependencies
......
TEMPLATE = subdirs
SUBDIRS = sources \
external/quadprog++
external/quadprog++-v2
sources.depends = external/quadprog++
sources.depends = external/quadprog++-v2
TEMPLATE = subdirs
SUBDIRS = quadprog++
SUBDIRS = quadprog++-v2
// This file is part of QuadProg++: a C++ library implementing
// the algorithm of Goldfarb and Idnani for the solution of a (convex)
// Quadratic Programming problem by means of an active-set dual method.
// Copyright (C) 2007-2009 Luca Di Gaspero.
// Copyright (C) 2009 Eric Moyer.
//
// QuadProg++ is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// QuadProg++ is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with QuadProg++. If not, see <http://www.gnu.org/licenses/>.
#include "Array.hh"
/**
Index utilities
*/
namespace QuadProgPP{
std::set<unsigned int> seq(unsigned int s, unsigned int e)
{
std::set<unsigned int> tmp;
for (unsigned int i = s; i <= e; i++)
tmp.insert(i);
return tmp;
}
std::set<unsigned int> singleton(unsigned int i)
{
std::set<unsigned int> tmp;
tmp.insert(i);
return tmp;
}
}
This diff is collapsed.
This diff is collapsed.
/*
The quadprog_solve() function implements the algorithm of Goldfarb and Idnani
for the solution of a (convex) Quadratic Programming problem
by means of an active-set dual method.
The problem is in the form:
min 0.5 * x G x + g0 x
s.t.
CE^T x + ce0 = 0
CI^T x + ci0 >= 0
The matrix and vectors dimensions are as follows:
G: n * n
g0: n
CE: n * p
ce0: p
CI: n * m
ci0: m
x: n
The function will return the cost of the solution written in the x vector or
std::numeric_limits::infinity() if the problem is infeasible. In the latter case
the value of the x vector is not correct.
References: D. Goldfarb, A. Idnani. A numerically stable dual method for solving
strictly convex quadratic programs. Mathematical Programming 27 (1983) pp. 1-33.
Notes:
1. pay attention in setting up the vectors ce0 and ci0.
If the constraints of your problem are specified in the form
A^T x = b and C^T x >= d, then you should set ce0 = -b and ci0 = -d.
2. The matrix G is modified within the function since it is used to compute
the G = L^T L cholesky factorization for further computations inside the function.
If you need the original matrix G you should make a copy of it and pass the copy
to the function.
Author: Luca Di Gaspero
DIEGM - University of Udine, Italy
l.digaspero@uniud.it
http://www.diegm.uniud.it/digaspero/
The author will be grateful if the researchers using this software will
acknowledge the contribution of this function in their research papers.
LICENSE
This file is part of QuadProg++: a C++ library implementing
the algorithm of Goldfarb and Idnani for the solution of a (convex)
Quadratic Programming problem by means of an active-set dual method.
Copyright (C) 2007-2009 Luca Di Gaspero.
Copyright (C) 2009 Eric Moyer.
QuadProg++ is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
QuadProg++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with QuadProg++. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _QUADPROGPP
#define _QUADPROGPP
#include "Array.hh"
namespace QuadProgPP{
double solve_quadprog(Matrix<double>& G, Vector<double>& g0,
const Matrix<double>& CE, const Vector<double>& ce0,
const Matrix<double>& CI, const Vector<double>& ci0,
Vector<double>& x);
}
#endif // #define _QUADPROGPP
TARGET = quadprog++
TEMPLATE = lib
CONFIG *= static \
qt \
plugin \
eigen
DESTDIR = ../build/lib
INCLUDEPATH += ../..
HEADERS = QuadProg++.hh \
Array.hh
SOURCES = QuadProg++.cc \
Array.cc
This diff is collapsed.
TEMPLATE = subdirs
SUBDIRS = \
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_parallel \
rational_fitter_eigen \
rational_fitter_leastsquare \
rational_fitter_matlab \
rational_fitter_dca \
rational_function_chebychev \
rational_function_chebychev_opt \
rational_function_legendre \
rational_function_cosine \
nonlinear_fitter_eigen \
nonlinear_fitter_ceres \
nonlinear_fitter_ipopt \
nonlinear_fitter_nlopt \
nonlinear_fresnel_schlick \
nonlinear_fresnel_normalized_schlick \
nonlinear_fresnel_retroschlick \
nonlinear_shadowing_smith \
nonlinear_shadowing_schlick \
nonlinear_function_diffuse \
# nonlinear_function_microfacets \
nonlinear_function_abc \
nonlinear_function_beckmann \
nonlinear_function_retrobeckmann \
nonlinear_function_blinn \
nonlinear_function_retroblinn \
nonlinear_function_retroyoo \
nonlinear_function_ward \
nonlinear_function_spherical_gaussian \
nonlinear_function_lafortune \
nonlinear_function_isotropic_lafortune \
data_merl \
data_brdf_slice \
data_interpolant \
# data_interpolant_matlab \
# data_astm
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_parallel \
rational_fitter_eigen \
rational_fitter_leastsquare \
rational_fitter_matlab \
rational_fitter_dca \
rational_fitter_parsec_one \
rational_fitter_parsec_multi \
rational_function_chebychev \
rational_function_chebychev_opt \
rational_function_legendre \
rational_function_cosine \
nonlinear_fitter_eigen \
nonlinear_fitter_ceres \
nonlinear_fitter_ipopt \
nonlinear_fitter_nlopt \
nonlinear_fresnel_schlick \
nonlinear_fresnel_normalized_schlick \
nonlinear_fresnel_retroschlick \
nonlinear_shadowing_smith \
nonlinear_shadowing_schlick \
nonlinear_function_diffuse \
# nonlinear_function_microfacets \
nonlinear_function_abc \
nonlinear_function_beckmann \
nonlinear_function_retrobeckmann \
nonlinear_function_blinn \
nonlinear_function_retroblinn \
nonlinear_function_retroyoo \
nonlinear_function_ward \
nonlinear_function_spherical_gaussian \
nonlinear_function_lafortune \
nonlinear_function_isotropic_lafortune \
data_merl \
data_brdf_slice \
data_interpolant \
# data_interpolant_matlab \
# data_astm
#pragma once
#include <Eigen/SVD>
#include <Array.hh>
#include <QuadProg++.hh>
#include <core/rational_function.h>
#include <core/vertical_segment.h>
class quadratic_program
{
public:
//! \brief Constructor need to specify the number of coefficients
quadratic_program(int np, int nq, bool compute_delta = false) :
_np(np), _nq(nq), _compute_delta(compute_delta), CI(0.0, _np+_nq, 0)
{ }
//! \brief Remove the already defined constraints
void clear_constraints()
{
CI.resize(_np+_nq, 0);
}
//! \brief Add a constraint by specifying the vector
void add_constraints(const vec& c)
{
const int m = CI.nrows();
const int n = CI.ncols();
if(n > 0)
{
// Construct temp buffer
double* temp = new double[n*m];
for(int u=0; u<n; ++u)
{
for(int v=0; v<m; ++v)
{
temp[u*m + v] = CI[v][u];
}
}
// Resize matrix CI
CI.resize(m, n+1);
// Recopy data
for(int u=0; u<n+1; ++u)
{
if(u==n)
{
for(int v=0; v<m; ++v)
CI[v][u] = c[v];
}
else
{
for(int v=0; v<m; ++v)
CI[v][u] = temp[u*m + v];
}
}
delete[] temp;
}
else
{
// Resize matrix CI
CI.resize(m, 1);
// Recopy data
for(int u=0; u<m; ++u)
CI[n][u] = c[u];
}
}
//! \brief Provide the number of constraints
int nb_constraints() const
{
return CI.ncols();
}
//! Set the indices of the remaining data
void set_training_set(const std::list<unsigned int>& ts)
{
this->training_set = ts;
}
//! \brief Solves the quadratic program and update the p and
//! q vector if necessary.
inline bool solve_program(QuadProgPP::Vector<double>& x, double& delta, vec& p, vec& q)
{
bool solves_qp = solve_program(x, delta) ;
if(solves_qp)
{
double norm = 0.0;
for(int i=0; i<_np+_nq; ++i)
{
const double v = x[i];
norm += v*v ;
if(i < _np)
{
p[i] = v ;
}
else
{
q[i-_np] = v ;
}
}
return norm > 0.0;
}
else
{
return false ;
}
}
//! \brief Solves the quadratic program
inline bool solve_program(QuadProgPP::Vector<double>& v, double& delta)
{
const int m = CI.nrows();
const int n = CI.ncols();
QuadProgPP::Matrix<double> G (0.0, m, m) ;
QuadProgPP::Vector<double> g (0.0, m) ;
QuadProgPP::Vector<double> ci(0.0, n) ;
QuadProgPP::Matrix<double> CE(0.0, m, 0) ;
QuadProgPP::Vector<double> ce(0.0, 0) ;
if(_compute_delta)
{
// Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12)
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()(0) ;
delta = sigma_M / sigma_m ;
}
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<m; ++i)
{
G[i][i] = 1.0 ;
}
// Each constraint (fitting interval or point
// add another dimension to the constraint
// matrix
for(int i=0; i<n; ++i)
{
// Norm of the row vector
double norm = 0.0 ;
for(int j=0; j<m; ++j)
{
norm += CI[j][i]*CI[j][i] ; ;
}
// Set the c vector, will later be updated using the
// delta parameter.
ci[i] = -delta * sqrt(norm) ;
}
// Compute the solution
const double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, v);
bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
return solves_qp;
}
#define PACANOWSKI2012
//! \brief Test all the constraints of the data.
//! Add the sample that is farest away from the function.
bool test_constraints(int ny, const rational_function_1d* r, const vertical_segment* data)
{
#ifdef PACANOWSKI2012
int nb_failed = 0;
double max_dev = 0.0; // Maximum absolute distance of the current
// solution to the data.
std::list<unsigned int>::iterator max_ind;
vec cu, cl;
std::list<unsigned int>::iterator it;
for(it = training_set.begin(); it != training_set.end(); it++)
{
vec x, yl, yu;
data->get(*it, x, yl, yu);
vec y = r->value(x);
bool fail_upper = y[ny] > yu[ny];
bool fail_lower = y[ny] < yl[ny];
if(fail_lower || fail_upper)
{
const double dev = std::abs(0.5*(yu[ny]+yl[ny]) - y[ny]);
nb_failed++;
if(max_dev < dev)
{
get_constraint(x, yl, yu, ny, r, cu, cl);
max_dev = dev;
max_ind = it;
}
}
}
#ifdef DEBUG
std::cout << "<<TRACE>> " << nb_failed << " constraints where not satified." << std::endl;
std::cout << "<<TRACE>> an interval failed the test with distance = " << max_dev << std::endl;
#endif
if(nb_failed > 0)
{
add_constraints(cu);
add_constraints(cl);
training_set.erase(max_ind);
#ifdef DEBUG
std::cout << "<<DEBUG>> number of remaining training elements: " << training_set.size() << std::endl;
#endif
return false;
}
else
{
return true;
}
#else
int n = next_unmatching_constraint(0, ny, r, data);
if(n < data->size())
{
vec x, yl, yu;
data->get(n, x, yl, yu);
vec cu, cl;
get_constraint(x, yl, yu, ny, r, cu, cl);
add_constraints(cu);
add_constraints(cl);
return false;
}
else
{
return true;
}
#endif
}
//! \brief Generate two constraint vectors from a vertical segment and a
//! ration function type.
inline void get_constraint(const vec& xi, const vec& yl, const vec& yu, int ny,
const rational_function_1d* func,
vec& cu, vec& cl);
//! \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_1d* r,
const vertical_segment* data);
protected:
int _np, _nq;
bool _compute_delta;
QuadProgPP::Matrix<double> CI;
//! Contains the indices of the vertical segment unused during the
//! rational interpolation.
std::list<unsigned int> training_set;
};
inline void quadratic_program::get_constraint(const vec& xi, const vec& yl, const vec& yu,
int ny, const rational_function_1d* func,
vec& cu, vec& cl)
{
cu.resize(_np+_nq);
cl.resize(_np+_nq);
// Create two vector of constraints
for(int j=0; j<_np+_nq; ++j)
{
// 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 ;
}
}
}
int quadratic_program::next_unmatching_constraint(int i, int ny, const rational_function_1d* 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 data->size();
}
#include "rational_fitter.h"
#include <core/plugins_manager.h>
#include <Eigen/SVD>
#include <Array.hh>
#include <QuadProg++.hh>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <string>
#include <list>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "quadratic_program.h"