Commit 6c8b917e authored by Laurent Belcour's avatar Laurent Belcour

Add the cosine formula of the Chebychev poylnomial. And the quadprog fitter

constraints on the boundary.
parent 9ae03658
......@@ -53,6 +53,10 @@ class rational_function : public QObject, public function
static int estimate_dk(int k, int d);
static void populate(std::vector<int>& vec, int N, int M, int j);
//! \brief Output the rational function as a gnuplot file. It requires
//! the data object to output the function at the input location only.
virtual void save_gnuplot(const std::string& filename, const data* d, const arguments& args) const ;
protected: // functions
//! Convert a 1D index into a vector of degree for a
......@@ -61,11 +65,7 @@ class rational_function : public QObject, public function
std::vector<int> index2degree(int i) const ;
//! \brief Save the rational function to the rational format (see \ref formating).
virtual void save(const std::string& filename) const ;
//! \brief Output the rational function as a gnuplot file. It requires
//! the data object to output the function at the input location only.
virtual void save_gnuplot(const std::string& filename, const data* d, const arguments& args) const ;
virtual void save(const std::string& filename) const ;
//! \brief Output the rational function using a C++ function formating.
virtual void save_cpp(const std::string& filename, const arguments& args) const ;
......
......@@ -119,26 +119,27 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
// TODO Specify the delta in case
// Handle multiple dim
double dt = args.get_float("dt", 0.1);
if(std::abs(v[dimX() + i]) > 0.0)
{
#ifdef RELATIVE_ERROR
v[dimX() + dimY()+i] = v[dimX() + i] * (1.0 - dt) ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + dt) ;
#endif
v[dimX() + dimY()+i] = std::max(v[dimX() + i] - dt, 0.0) ;
v[dimX() + dimY()+i] = v[dimX() + i] * (1.0 - dt) ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + dt) ;
#else
if(v[dimX() + i] > dt)
{
v[dimX() + dimY()+i] = v[dimX() + i] - dt ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] + dt ;
}
else
{
v[dimX() + dimY()+i] = 0.0 ;
v[dimX() + 2*dimY()+i] = dt ;
}
#endif
}
}
// If data is not in the interval of fit
// TODO: Update to more dims
bool is_in = true ;
for(int i=0; i<dimX(); ++i)
{
......
......@@ -201,6 +201,7 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
delete rs[j];
}
std::cout << " \r";
if(min_l2 < std::numeric_limits<double>::max())
{
std::cout << "<<INFO>> min L2 = " << min_l2 << std::endl;
......
......@@ -34,7 +34,7 @@ function* rational_fitter_quadprog::provide_function() const
return new rational_function() ;
}
rational_fitter_quadprog::rational_fitter_quadprog() : QObject()
rational_fitter_quadprog::rational_fitter_quadprog() : QObject(), _boundary(1.0)
{
}
rational_fitter_quadprog::~rational_fitter_quadprog()
......@@ -100,7 +100,8 @@ void rational_fitter_quadprog::set_parameters(const arguments& args)
_max_np = args.get_float("np", 10) ;
_max_nq = args.get_float("nq", 10) ;
_min_np = args.get_float("min-np", _max_np) ;
_min_nq = args.get_float("min-nq", _max_nq) ;
_min_nq = args.get_float("min-nq", _max_nq) ;
_boundary = args.get_float("boundary-constraint", 1.0f);
}
......@@ -168,7 +169,13 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
double a0_norm = 0.0 ;
double a1_norm = 0.0 ;
vec xi = d->get(i) ;
vec xi = d->get(i) ;
bool is_boundary = false;
for(int i=0; i<d->dimX(); ++i)
{
is_boundary = is_boundary || (xi[i] <= d->min()[i]) || (xi[i] >= d->max()[i]);
}
// A row of the constraint matrix has this
// form: [p_{0}(x_i), .., p_{np}(x_i), -f(x_i) q_{0}(x_i), .., -f(x_i) q_{nq}(x_i)]
......@@ -195,6 +202,14 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
vec yl, yu ;
d->get(i, yl, yu) ;
// Add a constraints for boundary conditions
if(is_boundary)
{
vec mean = 0.5*(yl+yu);
yl = mean + _boundary * (yl - mean);
yu = mean + _boundary * (yu - mean);
}
const double qi = r->q(xi, j-np) ;
a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
CI[j][i] = -yl[ny] * qi ;
......
......@@ -39,6 +39,8 @@ class rational_fitter_quadprog : public QObject, public fitter
virtual data* provide_data() const ;
virtual function* provide_function() const ;
private: // methods
protected: // data
// Fitting a data object using np elements in the numerator and nq
......@@ -49,5 +51,9 @@ class rational_fitter_quadprog : public QObject, public fitter
// min and Max usable np and nq values for the fitting
int _max_np, _max_nq ;
int _min_np, _min_nq ;
// Add constraints to the boundary of the domain. You can shrink it of
// the parameter --boundary-constraint *double*
double _boundary;
} ;
......@@ -27,6 +27,7 @@ rational_function_chebychev::~rational_function_chebychev()
double chebychev(double x, int i)
{
#ifdef RECURSIVE_FORM
if(i == 0)
{
return 1;
......@@ -39,6 +40,9 @@ double chebychev(double x, int i)
{
return 2.0*x*chebychev(x, i-1) - chebychev(x, i-2) ;
}
#else
return cos(i * acos(x));
#endif
}
// Get the p_i and q_j function
......@@ -67,5 +71,224 @@ double rational_function_chebychev::q(const vec& x, int i) const
return res ;
}
//! \todo it should handle parametrization
void rational_function_chebychev::save_matlab(const std::string& filename, const arguments& args) const
{
unsigned int np = a.size() / _nY ;
unsigned int nq = b.size() / _nY ;
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "function y = brdf(x)" << std::endl;
file << std::endl;
file << "\ts = [";
for(int i=0; i<dimX(); ++i)
{
file << 1.0 / (_max[i]-_min[i]);
if(i < dimX()-1)
{
file << ", ";
}
}
file << "];" << std::endl;
file << "\tc = [";
for(int i=0; i<dimX(); ++i)
{
file << _min[i];
if(i < dimX()-1)
{
file << ", ";
}
}
file << "];" << std::endl;
file << std::endl ;
// Export each color channel independantly
for(int j=0; j<dimY(); ++j)
{
// Export the numerator of the jth color channel
file << "\tp(" << j+1 << ",:) = ";
for(unsigned int i=0; i<np; ++i)
{
if(i > 0 && a[np*j + i] >= 0.0)
file << " + ";
else if(a[np*j + i] < 0.0)
file << " " ;
file << a[np*j + i];
std::vector<int> degree = index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << ".*chebychevpoly(1," << degree[k] << ", 2.0*((x(" << k+1 << ",:)"
<< "-c(" << k+1 << "))*s(" << k+1 << ") - 0.5))" ;
}
}
file << ";" << std::endl;
// Export the denominator of the jth color channel
file << "\tq(" << j+1 << ",:) = ";
for(unsigned int i=0; i<nq; ++i)
{
if(i > 0 && b[np*j + i] >= 0.0)
file << " + ";
else if(b[np*j + i] < 0.0)
file << " " ;
file << b[np*j + i] ;
std::vector<int> degree = index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << ".*chebychevpoly(1," << degree[k] << ", 2.0*((x(" << k+1 << ",:)"
<< "-c(" << k+1 << "))*s(" << k+1 << ") - 0.5))" ;
}
}
file << ";" << std::endl;
file << "\ty(" << j+1 << ",:) = p./q;" << std::endl;
if(j < dimY()-1)
{
file << std::endl;
}
}
file << "endfunction" << std::endl;
file.close() ;
}
//! \todo it should handle parametrization
void rational_function_chebychev::save_cpp(const std::string& filename, const arguments& args) const
{
unsigned int np = a.size() / _nY ;
unsigned int nq = b.size() / _nY ;
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "double s[" << dimX() << "] = {";
for(int i=0; i<dimX(); ++i)
{
file << 1.0 / (_max[i]-_min[i]);
if(i < dimX()-1)
{
file << ", ";
}
}
file << "};" << std::endl;
file << "double c[" << dimX() << "] = {";
for(int i=0; i<dimX(); ++i)
{
file << _min[i];
if(i < dimX()-1)
{
file << ", ";
}
}
file << "};" << std::endl;
file << std::endl ;
file << "// The Chebychev polynomial of order i evaluated in x" << std::endl;
file << "double l(double x, int i)" << std::endl;
file << "{" << std::endl;
file << "\treturn cos(i * acos(x));" << std::endl;
file << "}" << std::endl;
file << std::endl;
file << "void brdf(double* x, double* y)" << std::endl;
file << "{" << std::endl;
file << "\tdouble p, q;" << std::endl;
// Export each color channel independantly
for(int j=0; j<dimY(); ++j)
{
// Export the numerator of the jth color channel
file << "\tp = ";
for(unsigned int i=0; i<np; ++i)
{
if(i > 0 && a[np*j + i] >= 0.0)
{
file << " + ";
}
file << a[np*j + i];
std::vector<int> degree = index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << "*l(2.0*((x\[" << k << "\]-c[" << k << "])*s[" << k << "] - 0.5), " << degree[k] << ")" ;
}
}
file << ";" << std::endl;
// Export the denominator of the jth color channel
file << "\tq = ";
for(unsigned int i=0; i<nq; ++i)
{
if(i > 0 && b[np*j + i] >= 0.0)
file << " + ";
else if(b[np*j + i] < 0.0)
file << " " ;
file << b[np*j + i] ;
std::vector<int> degree = index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << "*l(2.0*((x\[" << k << "\]-c[" << k << "])*s[" << k << "] - 0.5), " << degree[k] << ")" ;
}
}
file << ";" << std::endl;
file << "\ty[" << j << "] = p/q;" << std::endl;
if(j < dimY()-1)
{
file << std::endl;
}
}
file << "}" << std::endl;
file.close() ;
}
void rational_function_chebychev::save(const std::string& filename) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl ;
file << "#NP " << a.size() / _nY << std::endl ;
file << "#NQ " << b.size() / _nY << std::endl ;
file << "#BASIS CHEBYCHEV" << std::endl ;
unsigned int np = a.size() / _nY ;
unsigned int nq = b.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
for(unsigned int i=0; i<np; ++i)
{
std::vector<int> index = index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
file << index[j] << "\t" ;
}
file << a[i+np*k] << std::endl ;
}
for(unsigned int i=0; i<nq; ++i)
{
std::vector<int> index = index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
file << index[j] << "\t" ;
}
file << b[i+nq*k] << std::endl ;
}
}
}
Q_EXPORT_PLUGIN2(rational_function_chebychev, rational_function_chebychev)
......@@ -24,5 +24,16 @@ class rational_function_chebychev : public rational_function
// Get the p_i and q_j function
virtual double p(const vec& x, int i) const ;
virtual double q(const vec& x, int j) const ;
protected: // methods
//! \brief Save the rational function to the rational format (see \ref formating).
virtual void save(const std::string& filename) const ;
//! \brief Output the rational function using a C++ function formating.
virtual void save_cpp(const std::string& filename, const arguments& args) const ;
//! \brief Output the rational function using a C++ function formating.
virtual void save_matlab(const std::string& filename, const arguments& args) const ;
} ;
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