Commit d9ea8068 authored by PACANOWSKI Romain's avatar PACANOWSKI Romain

Merging

parents 5c6bac5b d89b5c68
......@@ -56,7 +56,8 @@ NLOPT_DIR = ['#external/build/lib']
NLOPT_LIBS = ['nlopt']
NLOPT_OPT_LIBS = []
# MATLAB library and Engine
## MATLAB library and Engine
##
MATLAB_INC = ['/home/pac/MATLAB/R2013b/extern/include/']
MATLAB_DIR = ['/home/pac/MATLAB/R2013b/bin/glnxa64/']
MATLAB_LIBS = ['']
MATLAB_LIBS = ['engine']
import os
Import('env')
env = env.Clone()
if env.GetOption('clean'):
print "Nothing to do in external"
......
......@@ -214,7 +214,8 @@ function* plugins_manager::get_function(const std::string& filename)
// Parameters of the function object
int nX, nY;
params::input param_in; params::output param_out;
params::input param_in = params::UNKNOWN_INPUT;
params::output param_out = params::UNKNOWN_OUTPUT;
arguments args;
// Test for the first line of the file. Should be a ALTA FUNC HEADER
......@@ -331,7 +332,7 @@ function* plugins_manager::get_function(const arguments& args)
}
else
{
std::string filename = args["func"];
const std::string filename = args["func"];
FunctionPrototype myFunction = open_library<FunctionPrototype>(filename, "provide_function");
if(myFunction != NULL)
{
......
Import('env')
env = env.Clone()
env.AppendUnique(CPPPATH = env['MATLAB_INC'])
env.AppendUnique(LIBPATH = env['MATLAB_DIR'])
env.AppendUnique(LIBS = env['MATLAB_LIBS'])
env.AppendUnique(LIBS = ['core'])
conf = Configure(env)
sources = ['rational_fitter.cpp']
if conf.CheckLib(env['MATLAB_LIBS']):
env.SharedLibrary('../../build/rational_fitter_matlab', sources)
#end
conf.Finish()
......@@ -4,15 +4,19 @@ env = env.Clone()
env.AppendUnique(CPPPATH = env['QUADPROG_INC'])
env.AppendUnique(LIBPATH = env['QUADPROG_DIR'])
env.AppendUnique(LIBS = env['QUADPROG_LIBS'])
env.AppendUnique(CCFLAGS = env['OPENMP_FLAGS'])
env.AppendUnique(LIBS = env['OPENMP_LIBS'])
env.AppendUnique(LIBS = ['core'])
conf = Configure(env)
sources = ['rational_fitter.cpp']
if conf.CheckLib(env['QUADPROG_LIBS']) and conf.CheckLib(env['OPENMP_LIBS']):
if conf.CheckLib(env['QUADPROG_LIBS']):
env.AppendUnique(CCFLAGS = env['OPENMP_FLAGS'])
if conf.CheckLib(env['OPENMP_LIBS']):
env.AppendUnique(LIBS = env['OPENMP_LIBS'])
#end
env.SharedLibrary('../../build/rational_fitter_parallel', sources)
#end
conf.Finish()
......@@ -122,14 +122,16 @@ class quadratic_program
QuadProgPP::Vector<double> ci(0.0, n) ;
QuadProgPP::Matrix<double> CE(0.0, m, 0) ;
QuadProgPP::Vector<double> ce(0.0, 0) ;
/*
// 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 ;
/*/
delta = 1.0;
//*/
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<m; ++i)
......
......@@ -14,8 +14,9 @@
#include <cmath>
#include <string>
#include <list>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "quadratic_program.h"
......@@ -64,20 +65,21 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
timer time ;
time.start() ;
int nb_cores = args.get_int("nb-cores", omp_get_num_procs());
#ifdef _OPENMP
const int nb_cores = args.get_int("nb-cores", omp_get_num_procs());
#ifdef DEBUG
std::cout << "<<DEBUG>> will use " << nb_cores << " threads to compute the quadratic programs" << std::endl ;
#endif
omp_set_num_threads(nb_cores) ;
#endif
double min_delta = std::numeric_limits<double>::max();
double mean_delta = 0.0;
int nb_sol_found = 0;
int nb_sol_tested = 0;
#pragma omp parallel for shared(nb_sol_found, nb_sol_tested, min_delta, mean_delta), schedule(dynamic,1)
#pragma omp parallel for shared(args, nb_sol_found, nb_sol_tested, min_delta, mean_delta), schedule(dynamic,1)
for(int j=1; j<i; ++j)
{
// Compute the number of coefficients in the numerator and in the denominator
......@@ -90,18 +92,22 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
// Allocate a rational function and set it to the correct size, dimensions
// and parametrizations.
rational_function* rk = dynamic_cast<rational_function*>(plugins_manager::get_function(args));
rational_function* rk = NULL;
#pragma omp critical (args)
{
rk = dynamic_cast<rational_function*>(plugins_manager::get_function(args));
}
if(rk == NULL)
{
std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
throw;
}
rk->setParametrization(r->input_parametrization());
rk->setParametrization(r->output_parametrization());
rk->setDimX(r->dimX()) ;
rk->setDimY(r->dimY()) ;
rk->setMin(r->min()) ;
rk->setMax(r->max()) ;
if(rk == NULL)
{
std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
throw;
}
// Set the rational function size
rk->setSize(temp_np, temp_nq);
......@@ -135,7 +141,8 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
}
}
delete rk; // memory clean
if(rk != NULL)
delete rk; // memory clean
#pragma omp critical (nb_sol_tested)
{
......@@ -229,8 +236,10 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
do
{
#ifdef _OPENMP
#ifdef DEBUG
std::cout << "<<DEBUG>> thread " << omp_get_thread_num() << ", number of intervals tested = " << qp.nb_constraints()/2 << std::endl ;
#endif
#endif
QuadProgPP::Vector<double> x(n);
bool solves_qp = qp.solve_program(x, delta, p, q);
......
......@@ -217,10 +217,11 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int n
std::cout << "]" << std::endl ;
}
#endif
/*
// Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12)
Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(eCI, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(eCI);
const double sigma_m = svd.singularValues()(std::min(2*M, N)-1) ;
const double sigma_M = svd.singularValues()(0) ;
......@@ -235,6 +236,10 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int n
double delta = sigma_m / sigma_M ;
#ifndef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
if(isnan(delta) || (std::abs(delta) == std::numeric_limits<double>::infinity()))
{
std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
......@@ -244,10 +249,9 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int n
{
delta = 1.0 ;
}
#ifndef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
/*/
const double delta = 1.0;
//*/
for(int i=0; i<2*M; ++i)
{
ci[i] = ci[i] * delta ;
......
......@@ -29,11 +29,11 @@ rational_function_chebychev_1d::rational_function_chebychev_1d()
setDimY(0);
}
rational_function_chebychev_1d::rational_function_chebychev_1d(int np, int nq) :
rational_function_chebychev_1d::rational_function_chebychev_1d(int nX, int nY, int np, int nq) :
rational_function_1d(np, nq)
{
setDimX(0);
setDimY(0);
setDimX(nX);
setDimY(nY);
this->resize(np, nq);
}
......@@ -115,9 +115,7 @@ rational_function_1d* rational_function_chebychev::get(int i)
{
if(rs[i] == NULL)
{
rs[i] = new rational_function_chebychev_1d(np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
rs[i] = new rational_function_chebychev_1d(dimX(), dimY(), np, nq);
// Test if the input domain is not empty. If one dimension of
// the input domain is a point, I manually inflate this dimension
......
......@@ -17,7 +17,7 @@ class rational_function_chebychev_1d : public rational_function_1d
public: // methods
rational_function_chebychev_1d() ;
rational_function_chebychev_1d(int np, int nq) ;
rational_function_chebychev_1d(int nX, int nY, int np, int nq) ;
rational_function_chebychev_1d(const vec& a, const vec& b) ;
virtual ~rational_function_chebychev_1d() {}
......@@ -52,17 +52,19 @@ public: // methods
// Update the function
virtual void update(const vec& in_a,
const vec& in_b)
const vec& in_b)
{
rational_function_1d::update(in_a, in_b);
// Get the size of the input vector
const int np = in_a.size();
const int nq = in_b.size();
for(int k=0; k<np; ++k)
{
_p_coeffs[k].a = in_a[k];
}
const int nq = in_b.size();
for(int k=0; k<nq; ++k)
{
_q_coeffs[k].a = in_b[k];
......
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