Commit 03d5939d authored by Laurent Belcour's avatar Laurent Belcour

Reformating the rational function class to be more efficient.

parent b63dbe6d
......@@ -12,17 +12,14 @@ rational_function_1d::rational_function_1d()
{
}
rational_function_1d::rational_function_1d(int np, int nq, bool separable)
rational_function_1d::rational_function_1d(int nX, unsigned int np, unsigned int nq,
bool separable)
{
resize(np, nq);
_separable = separable;
}
setDimX(nX);
setDimY(1);
rational_function_1d::rational_function_1d(const vec& a,
const vec& b) : a(a), b(b)
{
//update(a, b);
_separable = false;
resize(np, nq);
_separable = separable;
}
bool rational_function_1d::load(std::istream&)
......@@ -34,70 +31,70 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c
{
bool is_matlab = args["export"] == "matlab";
bool is_alta = !args.is_defined("export") || args["export"] == "alta";
if(is_alta)
{
for(int i=0; i<a.size(); ++i)
{
std::vector<int> index = this->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<b.size(); ++i)
{
std::vector<int> index = this->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << b[i] << std::endl ;
}
}
else if(is_matlab)
{
out << "(";
for(int i=0; i<a.size(); ++i)
{
out << a[i];
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != a.size()-1) { out << " + "; }
}
out << ") / (";
for(int i=0; i<b.size(); ++i)
{
out << b[i] << "x.^" << i;
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != b.size()-1) { out << " + "; }
}
out << ")";
}
else
{
NOT_IMPLEMENTED();
}
const unsigned int np = _p_coeffs.size();
const unsigned int nq = _q_coeffs.size();
if(is_alta)
{
for(unsigned int i=0; i<np; ++i)
{
for(int j=0; j<dimX(); ++j)
{
out << _p_coeffs[i].deg[j] << "\t" ;
}
out << _p_coeffs[i].a << std::endl ;
}
for(unsigned int i=0; i<nq; ++i)
{
for(int j=0; j<dimX(); ++j)
{
out << _q_coeffs[i].deg[j] << "\t" ;
}
out << _q_coeffs[i].a << std::endl ;
}
}
else if(is_matlab)
{
out << "(";
for(unsigned int i=0; i<np; ++i)
{
out << _p_coeffs[i].a;
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << _q_coeffs[i].deg[k];
}
if(i != np-1) { out << " + "; }
}
out << ") / (";
for(unsigned int i=0; i<nq; ++i)
{
out << _p_coeffs[i].a << "x.^" << i;
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << _q_coeffs[i].deg[k];
}
if(i != nq-1) { out << " + "; }
}
out << ")";
}
else
{
NOT_IMPLEMENTED();
}
}
void rational_function_1d::update(const vec& in_a,
const vec& in_b)
{
a.resize(in_a.size()) ;
b.resize(in_b.size()) ;
// Get the size of the input vector
const int np = in_a.size();
const int nq = in_b.size();
#define NORMALIZE
......@@ -107,22 +104,45 @@ void rational_function_1d::update(const vec& in_a,
const double b0 = 1.0;
#endif
for(int i=0; i<a.size(); ++i) { a[i] = in_a[i] / b0; }
for(int i=0; i<b.size(); ++i) { b[i] = in_b[i] / b0; }
for(int k=0; k<np; ++k)
{
_p_coeffs[k].a = in_a[k] / b0;
}
for(int k=0; k<nq; ++k)
{
_q_coeffs[k].a = in_b[k] / b0;
}
}
void rational_function_1d::resize(int np, int nq)
void rational_function_1d::resize(unsigned int np, unsigned int nq)
{
const int old_np = a.size();
const int old_nq = b.size();
// It is not possible to resize the multi-dimensional vector
// if the input dimensions size is not defined. This can
// happen at the creation of the rational function object.
if(dimX() == 0) { return; }
// Resize the vector
a.resize(np);
b.resize(nq);
// Resize the numerator
if(_p_coeffs.size() != np)
{
_p_coeffs.resize(np);
for(unsigned int k=0; k<np; ++k)
{
std::vector<int> deg = index2degree(k);
_p_coeffs[k].deg = deg;
}
}
// Set the new coeffs to zero
for(int i=old_np; i<np; ++i) { a[i] = 0.0; }
for(int i=old_nq; i<nq; ++i) { b[i] = 0.0; }
// Resize the denominator
if(_q_coeffs.size() != nq)
{
_q_coeffs.resize(nq);
for(unsigned int k=0; k<nq; ++k)
{
std::vector<int> deg = index2degree(k);
_q_coeffs[k].deg = deg;
}
}
}
......@@ -130,40 +150,26 @@ void rational_function_1d::resize(int np, int nq)
// Get the p_i and q_j function
vec rational_function_1d::p(const vec& x) const
{
vec res(_nY) ;
unsigned int const np = a.size() / _nY ;
vec res(1) ;
for(int k=0; k<_nY; ++k)
const unsigned int np = _p_coeffs.size() ;
for(unsigned int i=0; i<np; ++i)
{
double p = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[k*_nY + i]*this->p(x, i) ;
}
res[k] = p ;
res[0] += _p_coeffs[i].a * this->p(x, i) ;
}
return res ;
}
vec rational_function_1d::q(const vec& x) const
{
vec res(_nY) ;
unsigned int const nq = b.size() / _nY ;
vec res(1) ;
for(int k=0; k<_nY; ++k)
const unsigned int np = _q_coeffs.size() ;
for(unsigned int i=0; i<np; ++i)
{
double q = 0.0f ;
for(unsigned int i=0; i<nq; ++i)
{
q += b[k*_nY + i]*this->q(x, i) ;
}
res[k] = q ;
res[0] += _q_coeffs[i].a * this->p(x, i) ;
}
return res ;
}
......@@ -287,19 +293,25 @@ std::vector<int> rational_function_1d::index2degree(int i) const
// Get the p_i and q_j function
double rational_function_1d::p(const vec& x, int i) const
{
std::vector<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
res *= pow(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]) ;
//res *= pow(x[k], deg[k]) ;
{
const double xp = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= pow(xp, _p_coeffs[i].deg[k]) ;
}
return res ;
}
double rational_function_1d::q(const vec& x, int i) const
{
return p(x, i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
const double xp = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= pow(xp, _q_coeffs[i].deg[k]) ;
}
return res ;
}
// Overload the function operator
......@@ -307,20 +319,20 @@ vec rational_function_1d::value(const vec& x) const
{
vec res(1) ;
unsigned int const np = a.size() / _nY ;
unsigned int const nq = b.size() / _nY ;
unsigned int const np = _p_coeffs.size();
unsigned int const nq = _q_coeffs.size();
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[i]*this->p(x, i) ;
p += _p_coeffs[i].a*this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += b[i]*this->q(x, i) ;
q += _q_coeffs[i].a*this->q(x, i) ;
}
res[0] = p/q ;
......@@ -330,25 +342,28 @@ vec rational_function_1d::value(const vec& x) const
std::ostream& operator<< (std::ostream& out, const rational_function_1d& r)
{
const unsigned int np = r._p_coeffs.size();
const unsigned int nq = r._q_coeffs.size();
std::cout << "p = [" ;
for(int i=0; i<r.a.size(); ++i)
for(unsigned int i=0; i<np; ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.a[i] ;
std::cout << r._p_coeffs[i].a ;
}
std::cout << "]" << std::endl ;
std::cout << "q = [" ;
for(int i=0; i<r.b.size(); ++i)
for(unsigned int i=0; i<nq; ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.b[i] ;
std::cout << r._q_coeffs[i].a ;
}
std::cout << "]" << std::endl ;
......@@ -384,9 +399,7 @@ rational_function_1d* rational_function::get(int i)
{
if(rs[i] == NULL)
{
rs[i] = new rational_function_1d(np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
rs[i] = new rational_function_1d(dimX(), 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
......
......@@ -21,28 +21,28 @@ class rational_function_1d : public function
public: // methods
rational_function_1d() ;
rational_function_1d(int np, int nq, bool separable = false) ;
rational_function_1d(const vec& a, const vec& b) ;
rational_function_1d(int nX, unsigned int np, unsigned int nq,
bool separable = false) ;
virtual ~rational_function_1d() {}
/* FUNCTION INHERITANCE */
/* FUNCTION INHERITANCE */
// Overload the function operator
//! Overload the function operator
virtual vec value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files
virtual bool load(std::istream& in);
//! IO function to text files
virtual bool load(std::istream& in);
//! \brief Save the rational function expansion. It should
//! not be store in any variable (e.g. "y = rf(x);") as the
//! nD rational function can append factor to the 1D rational
//! function.
virtual void save_body(std::ostream&, const arguments&) const;
//! \brief Save the rational function expansion. It should
//! not be store in any variable (e.g. "y = rf(x);") as the
//! nD rational function can append factor to the 1D rational
//! function.
virtual void save_body(std::ostream&, const arguments&) const;
/* RATIONAL FUNCTION SPECIFIC */
/* RATIONAL FUNCTION SPECIFIC */
//! Evaluate the numerator \f$p(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast
......@@ -50,6 +50,7 @@ class rational_function_1d : public function
//! algorithm to efficiently evaluate recursively defined
//! polynomials.
virtual vec p(const vec& x) const ;
//! Evaluate the denominator \f$q(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast
//! implementation. For example one can use the Clenshaw
......@@ -64,19 +65,41 @@ class rational_function_1d : public function
//! denominator of the rational function.
virtual double q(const vec& x, int j) const ;
// Update the function
//! Update the coefficient vectors with new values. The new values
//! are normalized by the first element of the denominator
//! coefficients.
virtual void update(const vec& in_a,
const vec& in_b) ;
// Resize the polynomial
virtual void resize(int np, int nq);
//! Resize the polynomial.
virtual void resize(unsigned int np, unsigned int nq);
// Get the coefficients
virtual double getP(int i) const { return a[i]; }
virtual double getQ(int i) const { return b[i]; }
virtual vec getP() const { return a; }
virtual vec getQ() const { return b; }
//! Get the i-th coefficient of the numerator.
virtual double getP(int i) const { return _p_coeffs[i].a; }
//! Get the i-th coefficient of the denominator.
virtual double getQ(int i) const { return _q_coeffs[i].a; }
//! Get the vector of coefficient for the numerator.
virtual vec getP() const
{
const int np = _p_coeffs.size();
vec t(np);
for(int i=0; i<np; ++i) {t[i] = _p_coeffs[i].a; }
return t;
}
//! Get the vector of coefficient for the denominator.
virtual vec getQ() const
{
const int nq = _q_coeffs.size();
vec t(nq);
for(int i=0; i<nq; ++i) {t[i] = _q_coeffs[i].a; }
return t;
}
// STL stream ouput
friend std::ostream& operator<< (std::ostream& out,
......@@ -96,14 +119,28 @@ class rational_function_1d : public function
protected: // data
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
vec a, b ;
//! Is the function separable with respect to its input dimensions?
//! \todo Make possible to have only part of the dimensions
//! separable.
bool _separable;
// Structure to store a multi-dimensional coefficient. The
// coeffcient, a, is associated with the vector of degree
// for each dimension.
struct coeff
{
coeff() {}
coeff(double a, std::vector<int> deg) :
a(a), deg(deg) { }
double a;
std::vector<int> deg;
};
// Table of coefficients and indices, sorted with respect
// to the indices.
std::vector<coeff> _p_coeffs;
std::vector<coeff> _q_coeffs;
//! Is the function separable with respect to its input dimensions?
//! \todo Make possible to have only part of the dimensions
//! separable.
bool _separable;
} ;
class rational_function : public function
......
......@@ -27,14 +27,8 @@ rational_function_chebychev_1d::rational_function_chebychev_1d()
{
}
rational_function_chebychev_1d::rational_function_chebychev_1d(int np, int nq) :
rational_function_1d(np, nq)
{
}
rational_function_chebychev_1d::rational_function_chebychev_1d(const vec& a,
const vec& b) :
rational_function_1d(a, b)
rational_function_chebychev_1d::rational_function_chebychev_1d(int nX, int np, int nq) :
rational_function_1d(nX, np, nq)
{
}
......@@ -92,9 +86,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(), 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,8 +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(const vec& a, const vec& b) ;
rational_function_chebychev_1d(int nX, int np, int nq) ;
virtual ~rational_function_chebychev_1d() {}
// Get the p_i and q_j function
......
......@@ -29,41 +29,12 @@ rational_function_chebychev_1d::rational_function_chebychev_1d()
setDimY(0);
}
rational_function_chebychev_1d::rational_function_chebychev_1d(int nX, int nY, int np, int nq) :
rational_function_1d(np, nq)
rational_function_chebychev_1d::rational_function_chebychev_1d(int nX, int np, int nq) :
rational_function_1d(nX, np, nq)
{
setDimX(nX);
setDimY(nY);
this->resize(np, nq);
}
rational_function_chebychev_1d::rational_function_chebychev_1d(const vec& a,
const vec& b) :
rational_function_1d(a, b)
{
setDimX(0);
setDimY(0);
const int np = a.size();
const int nq = b.size();
this->resize(np, nq);
// Update the numerator coefficient array
for(int k=0; k<np; ++k)
{
_p_coeffs[k].a = a[k];
}
// Update the denominator coefficient array
for(int k=0; k<nq; ++k)
{
_q_coeffs[k].a = b[k];
}
}
double chebychev(double x, int i)
{
#ifdef RECURSIVE_FORM
......@@ -115,7 +86,7 @@ rational_function_1d* rational_function_chebychev::get(int i)
{
if(rs[i] == NULL)
{
rs[i] = new rational_function_chebychev_1d(dimX(), dimY(), np, nq);
rs[i] = new rational_function_chebychev_1d(dimX(), 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,8 +17,7 @@ class rational_function_chebychev_1d : public rational_function_1d
public: // methods
rational_function_chebychev_1d() ;
rational_function_chebychev_1d(int nX, int nY, int np, int nq) ;
rational_function_chebychev_1d(const vec& a, const vec& b) ;
rational_function_chebychev_1d(int nX, int np, int nq) ;
virtual ~rational_function_chebychev_1d() {}
// Get the p_i and q_j function
......@@ -30,8 +29,8 @@ public: // methods