Commit 0dfeeee6 authored by Laurent Belcour's avatar Laurent Belcour
Browse files

Correction of the Rational Polynomial Chebychev plugin

parent 94a6cfff
......@@ -285,6 +285,7 @@ std::ostream& operator<< (std::ostream& out, const rational_function_1d& r)
return out ;
}
#ifndef TODO
rational_function::rational_function() : np(0), nq(0)
{
}
......@@ -442,235 +443,6 @@ void rational_function::load(std::istream& in)
}
}
//! \todo it should handle parametrization
void rational_function::save_matlab(const std::string& filename, const arguments&) const
{
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)
{
rational_function_1d* rf = get(j);
vec a = rf->getP();
vec b = rf->getQ();
// Export the numerator of the jth color channel
file << "\tp(" << j+1 << ",:) = ";
for(int i=0; i<np; ++i)
{
if(i > 0 && a[i] >= 0.0)
{
file << " + ";
}
else if(a[i] < 0.0)
{
file << " " ;
}
file << a[i];
std::vector<int> degree = rf->index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << ".*legendrepoly(" << 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(int i=0; i<nq; ++i)
{
if(i > 0 && b[i] >= 0.0)
{
file << " + ";
}
else if(b[i] < 0.0)
{
file << " " ;
}
file << b[i] ;
std::vector<int> degree = rf->index2degree(i);
for(unsigned int k=0; k<degree.size(); ++k)
{
file << ".*legendrepoly(" << 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::save_cpp(const std::string& filename, const arguments&) const
{
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 Legendre polynomial of order i evaluated in x" << std::endl;
file << "double l(double x, int i)" << std::endl;
file << "{" << std::endl;
file << " if(i == 0)" << std::endl;
file << " {" << std::endl;
file << " return 1;" << std::endl;
file << " }" << std::endl;
file << " else if(i == 1)" << std::endl;
file << " {" << std::endl;
file << " return x;" << std::endl;
file << " }" << std::endl;
file << " else" << std::endl;
file << " {" << std::endl;
file << " return ((2*i-1)*x*l(x, i-1) - (i-1)*l(x, i-2)) / (double)i ;" << std::endl;
file << " }" << 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)
{
rational_function_1d* rf = get(j);
vec a = rf->getP();
vec b = rf->getQ();
// Export the numerator of the jth color channel
file << "\tp = ";
for(int i=0; i<np; ++i)
{
if(i > 0 && a[i] >= 0.0)
{
file << " + ";
}
else if(a[i] < 0.0)
{
file << " " ;
}
file << a[i];
std::vector<int> degree = rf->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(int i=0; i<nq; ++i)
{
if(i > 0 && b[i] >= 0.0)
{
file << " + ";
}
else if(b[i] < 0.0)
{
file << " " ;
}
file << b[i] ;
std::vector<int> degree = rf->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::save_gnuplot(const std::string& filename, const data* d, const arguments&) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
for(int i=0; i<d->size(); ++i)
{
vec v = d->get(i) ;
// vec y1 ; y1.assign(d->dimY(), 0.0) ;
// for(int k=0; k<d->dimY(); ++k) { y1[k] = v[d->dimX() + k] ; }
vec y2 = value(v) ;
for(int u=0; u<d->dimX(); ++u)
file << v[u] << "\t" ;
for(int u=0; u<d->dimY(); ++u)
file << y2[u] << "\t" ;
file << std::endl ;
}
file.close();
}
void rational_function::save_call(std::ostream& out, const arguments&) const
{
......@@ -709,45 +481,6 @@ void rational_function::save_call(std::ostream& out, const arguments&) const
}
}
void rational_function::save(const std::string& filename) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl ;
file << "#NP " << np << std::endl ;
file << "#NQ " << nq << std::endl ;
file << "#BASIS LEGENDRE" << std::endl ;
file << "#INPUT_PARAM " << params::get_name(this->parametrization()) << std::endl;
for(int k=0; k<_nY; ++k)
{
rational_function_1d* rf = get(k);
vec a = rf->getP();
vec b = rf->getQ();
for(int i=0; i<np; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
file << index[j] << "\t" ;
}
file << a[i] << std::endl ;
}
for(int i=0; i<nq; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
file << index[j] << "\t" ;
}
file << b[i] << std::endl ;
}
}
}
std::ostream& operator<< (std::ostream& out, rational_function& r)
{
for(int i=0; i<r.dimY(); ++i)
......@@ -768,4 +501,4 @@ std::ostream& operator<< (std::ostream& out, rational_function& r)
}
#endif
......@@ -36,7 +36,7 @@ class rational_function_1d : public function
virtual double q(const vec& x, int j) const ;
// IO function to text files
virtual void load(std::istream& in);
virtual void load(std::istream& in);
// Update the function
virtual void update(const vec& in_a,
......@@ -75,6 +75,279 @@ class rational_function_1d : public function
vec a, b ;
} ;
#ifdef TODO
template<class T> class rational_function_t : public function
{
public: // methods
rational_function_t() : np(0), nq(0) {}
rational_function_t(int np, int nq) : np(np), nq(nq) {}
virtual ~rational_function_t() {}
// Overload the function operator
virtual vec value(const vec& x) const
{
vec res(_nY) ;
for(int k=0; k<_nY; ++k)
{
res[k] = rs[k]->value(x)[0] ;
}
return res ;
}
virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files
virtual void load(std::istream& in)
{
// Parse line until the next comment
while(in.peek() != '#')
{
char line[256];
in.getline(line, 256);
}
// Checking for the comment line #FUNC nonlinear_function_lafortune
std::string token;
in >> token;
if(token.compare("#FUNC") != 0)
{
std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl;
}
in >> token;
if(token.compare("rational_function") != 0)
{
std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl;
}
int _np, _nq;
// Shoudl have the #NP [int]
in >> token >> _np;
// Shoudl have the #NQ [int]
in >> token >> _nq;
setSize(_np, _nq);
// Check for the MIN and MAX vector
vec min(dimX()), max(dimX());
in >> token;
if(token.compare("#MIN") != 0)
{
std::cerr << "<<ERROR>> the min value for the input space is not defined." << std::endl;
}
for(int k=0; k<dimX(); ++k) {in >> min[k];}
setMin(min);
in >> token;
if(token.compare("#MAX") != 0)
{
std::cerr << "<<ERROR>> the max value for the input space is not defined." << std::endl;
}
for(int k=0; k<dimX(); ++k) {in >> max[k]; }
setMax(max);
// Check for the polynomial basis type
in >> token;
if(token.compare("#BASIS") != 0)
{
std::cerr << "<<ERROR>> the file is not specifying the polynomial basis." << std::endl;
}
in >> token;
if(token.compare("LEGENDRE") != 0)
{
std::cerr << "<<ERROR>> the basis is different than LEGENDRE." << std::endl;
}
vec a(_np), b(_nq);
for(int i=0; i<_nY; ++i)
{
// Parse the p_i coefficients
for(int j=0; j<_np; ++j)
{
in >> token >> a[j];
}
// Parse the q_i coefficients
for(int j=0; j<_nq; ++j)
{
in >> token >> b[j];
}
std::cout << a << std::endl;
std::cout << b << std::endl;
// Update the i_th color channel
get(i)->update(a, b);
}
}
#ifdef OLD
// Update the function
virtual void update(int i, T* r) ;
#endif
//! Get the 1D function associated with color channel i. If no one exist,
//! this function allocates a new element. If i > nY, it returns NULL.
T* get(int i)
{
// Check for consistency in the index of color channel
if(i < _nY)
{
if(rs[i] == NULL)
{
rs[i] = new T(np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
rs[i]->setMin(getMin()) ;
rs[i]->setMax(getMax()) ;
}
return rs[i];
}
else
{
std::cout << "<<ERROR>> tried to access out of bound 1D RF"
<< std::endl;
return NULL;
}
}
T* get(int i) const
{
// Check for consistency in the index of color channel
if(i < _nY)
{
return rs[i];
}
else
{
std::cout << "<<ERROR>> tried to access out of bound 1D RF"
<< std::endl;
return NULL;
}
}
#ifdef OLD
// STL stream ouput
friend std::ostream& operator<<(std::ostream& out, rational_function_t& r)
{
for(int i=0; i<r.dimY(); ++i)
{
T* rf = r.get(i);
out << "dimension " << i << ": ";
if(rf != NULL)
{
out << *rf << std::endl;
}
else
{
out << "[NULL]" << std::endl;
}
}
return out ;
}
#endif
//! Set the dimension of the output space of the function. This function
//! will update the size of the rs vector size.
virtual void setDimY(int nY)
{
_nY = nY ;
rs.resize(nY);
}
//! \brief Set the size of the rational function. Any newly created 1D
//! function will have np and nq fixed as defined.
virtual void setSize(int np, int nq)
{
this->np = np;
this->nq = nq;
clear();
}
//! \brief Clear the vector of 1D rational functions.
virtual void clear()
{
rs.clear();
rs.resize(_nY);
}
virtual void setMin(const vec& min)
{
function::setMin(min);
for(int i=0; i<dimY(); ++i)
{
get(i)->setMin(min);
}
}
virtual void setMax(const vec& max)
{
function::setMax(max);
for(int i=0; i<dimY(); ++i)
{
get(i)->setMax(max);
}
}
//! \brief Save the rational function to the rational format (see \ref
//! formating).
virtual void save_call(std::ostream& out, const arguments& args) const
{
out << "#FUNC rational_function" << std::endl;
out << "#NP " << np << std::endl ;
out << "#NQ " << nq << std::endl ;
out << "#MIN "; for(int k=0; k<_nX; ++k) { out << _min[k] << " "; }
out << std::endl;
out << "#MAX "; for(int k=0; k<_nX; ++k) { out << _max[k] << " "; }
out << std::endl;
out << "#BASIS LEGENDRE" << std::endl ;
for(int k=0; k<_nY; ++k)
{
T* rf = get(k);
vec a = rf->getP();
vec b = rf->getQ();
for(int i=0; i<np; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << a[i] << std::endl ;
}
for(int i=0; i<nq; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << b[i] << std::endl ;
}
}
}
protected: // data
//! Store the y \in R rational functions. Each channel is a distinct
//! polynomial and should be fitted separately.
std::vector<T*> rs ;
//! Size of the polynomials
//! \todo Change it by a more adaptive scheme, with different np, nq per
//! color channel?
int np, nq;
}
typedef rational_function rational_function_t<rational_function_1d>;
#else
class rational_function : public function
{
public: // methods
......@@ -88,34 +361,29 @@ class rational_function : public function
virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files
virtual void load(std::istream& in) ;
virtual void load(std::istream& in) ;
// Update the function
virtual void update(int i, rational_function_1d* r) ;
//! Get the 1D function associated with color channel i. If no one exist, this
//! function allocates a new element. If i > nY, it returns NULL.
//! Get the 1D function associated with color channel i. If no one exist,
//! this function allocates a new element. If i > nY, it returns NULL.
virtual rational_function_1d* get(int i) ;
virtual rational_function_1d* get(int i) const ;
// STL stream ouput
friend std::ostream& operator<<(std::ostream& out, rational_function& r) ;