Commit 65a9c4a8 authored by Laurent Belcour's avatar Laurent Belcour

Got a working class QuadraticPrograming. Next step is the iterative scheme.

parent 7af7ad29
......@@ -17,22 +17,59 @@ public:
}
//! \brief Add a constraint by specifying the vector
void add_constraints(const QuadProgPP::Vector<double> v)
void add_constraints(const QuadProgPP::Vector<double> vec)
{
const int m = CI.nrows();
const int n = CI.ncols();
CI.resize(m, n+1);
for(int i=0; i<m; ++i)
if(n > 0)
{
CI[i][n] = v[i];
// 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] = vec[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] = vec[u];
}
}
//! \brief Solves the quadratic program
bool solve_program(QuadProgPP::Vector<double>& v)
{
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) ;
......@@ -44,7 +81,7 @@ public:
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) ;
const double delta = sigma_M / sigma_m ;
delta = sigma_M / sigma_m ;
// Select the size of the result vector to
// be equal to the dimension of p + q
......
......@@ -56,7 +56,7 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit)
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 << "\r" ;
std::cout.flush() ;
QTime time ;
time.start() ;
......@@ -171,11 +171,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
}
#ifdef OLD
// Get the maximum value in data to scale the input parameter space
// so that it reduces the values of the polynomial
vec dmax = d->max() ;
// Matrices of the problem
// Matrices of the problem
QuadProgPP::Matrix<double> G (0.0, np+nq, np+nq) ;
QuadProgPP::Vector<double> g (0.0, np+nq) ;
QuadProgPP::Matrix<double> CI(0.0, np+nq, 2*d->size()) ;
......@@ -247,20 +243,10 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
ci[i+d->size()] = -sqrt(a1_norm) ;
}
#ifdef DEBUG
std::cout << "CI = [" ;
for(int j=0; j<d->size()*2; ++j)
{
for(int i=0; i<np+nq; ++i)
{
std::cout << CI[i][j] ;
if(i != np+nq-1) std::cout << ", ";
}
if(j != d->size()*2-1)
std::cout << ";" << std::endl;
else
std::cout << "]" << std::endl ;
}
#ifndef DEBUG
std::cout << CI << std::endl ;
#endif
// Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12)
......@@ -352,7 +338,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
}
QuadProgPP::Vector<double> x(np+nq);
bool solves_qp = qp.solve_program(x);
bool solves_qp = qp.solve_program(x, delta);
#endif
if(solves_qp)
{
......@@ -372,7 +358,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
}
}
#ifdef DEBUG
#ifndef DEBUG
std::cout << "<<INFO>> got solution " << *r << std::endl ;
#endif
return norm > 0.0;
......
......@@ -24,12 +24,23 @@ int main(int argc, char** argv)
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
const float y = 100.0f * exp(-10.0 * x*x) * x*x - 0.01 *x*x*x ;
const float y = 100.0f * exp(-10.0 * x*x) * x*x - 0.01 *x*x*x + 0.1 ;
f << x << "\t" << y << "\t" << 0.1f << std::endl ;
}
}
else if(k == 2)
else if(k == 2)
{
f << "#DIM 1 1" << std::endl ;
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
const float y = (1.0 + 7.0*x - 10.5*x*x) / (1.0 + 7.0 * x) ;
f << x << "\t" << y << "\t" << 0.1f << std::endl ;
}
}
else if(k == 3)
{
f << "#DIM 2 1" << std::endl ;
for(int i=0; i<nbx; ++i)
......@@ -42,7 +53,7 @@ int main(int argc, char** argv)
f << x << "\t" << y << "\t" << z << "\t" << 0.1f << std::endl ;
}
}
else if(k == 3)
else if(k == 4)
{
f << "#DIM 2 1" << std::endl ;
for(int i=0; i<nbx; ++i)
......
......@@ -92,7 +92,7 @@ int main(int argc, char** argv)
size_t n = args["output"].find('.') ;
std::string gnuplot_filename = args["output"].substr(0,n);
gnuplot_filename.append(".gnuplot") ;
/*
/*
f->save_gnuplot(gnuplot_filename, d, args);
/*/
std::ofstream file(gnuplot_filename.c_str(), std::ios_base::trunc);
......@@ -110,7 +110,8 @@ int main(int argc, char** argv)
file << y2[u] << "\t" ;
file << std::endl ;
}
}
file.close();
//*/
std::string error_filename = args["output"].substr(0,n);
error_filename.append(".errorplot") ;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment