Commit 8825a87b authored by Laurent Belcour's avatar Laurent Belcour
Browse files

Merge

parents c836654a 03d5939d
...@@ -12,19 +12,16 @@ rational_function_1d::rational_function_1d() ...@@ -12,19 +12,16 @@ 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)
{ {
setDimX(nX);
setDimY(1);
resize(np, nq); resize(np, nq);
_separable = separable; _separable = separable;
} }
rational_function_1d::rational_function_1d(const vec& a,
const vec& b) : a(a), b(b)
{
//update(a, b);
_separable = false;
}
bool rational_function_1d::load(std::istream&) bool rational_function_1d::load(std::istream&)
{ {
return true; return true;
...@@ -35,55 +32,54 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c ...@@ -35,55 +32,54 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c
bool is_matlab = args["export"] == "matlab"; bool is_matlab = args["export"] == "matlab";
bool is_alta = !args.is_defined("export") || args["export"] == "alta"; bool is_alta = !args.is_defined("export") || args["export"] == "alta";
const unsigned int np = _p_coeffs.size();
const unsigned int nq = _q_coeffs.size();
if(is_alta) if(is_alta)
{ {
for(int i=0; i<a.size(); ++i) for(unsigned int i=0; i<np; ++i)
{ {
std::vector<int> index = this->index2degree(i) ; for(int j=0; j<dimX(); ++j)
for(unsigned int j=0; j<index.size(); ++j)
{ {
out << index[j] << "\t" ; out << _p_coeffs[i].deg[j] << "\t" ;
} }
out << a[i] << std::endl ; out << _p_coeffs[i].a << std::endl ;
} }
for(int i=0; i<b.size(); ++i) for(unsigned int i=0; i<nq; ++i)
{ {
std::vector<int> index = this->index2degree(i) ; for(int j=0; j<dimX(); ++j)
for(unsigned int j=0; j<index.size(); ++j)
{ {
out << index[j] << "\t" ; out << _q_coeffs[i].deg[j] << "\t" ;
} }
out << b[i] << std::endl ; out << _q_coeffs[i].a << std::endl ;
} }
} }
else if(is_matlab) else if(is_matlab)
{ {
out << "("; out << "(";
for(int i=0; i<a.size(); ++i) for(unsigned int i=0; i<np; ++i)
{ {
out << a[i]; out << _p_coeffs[i].a;
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k) for(int k=0; k<dimX(); ++k)
{ {
if(k != dimX()-1) { out << ".*"; } if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k]; out << "x(k).^" << _q_coeffs[i].deg[k];
} }
if(i != a.size()-1) { out << " + "; } if(i != np-1) { out << " + "; }
} }
out << ") / ("; out << ") / (";
for(int i=0; i<b.size(); ++i) for(unsigned int i=0; i<nq; ++i)
{ {
out << b[i] << "x.^" << i; out << _p_coeffs[i].a << "x.^" << i;
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k) for(int k=0; k<dimX(); ++k)
{ {
if(k != dimX()-1) { out << ".*"; } if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k]; out << "x(k).^" << _q_coeffs[i].deg[k];
} }
if(i != b.size()-1) { out << " + "; } if(i != nq-1) { out << " + "; }
} }
out << ")"; out << ")";
} }
...@@ -96,8 +92,9 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c ...@@ -96,8 +92,9 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c
void rational_function_1d::update(const vec& in_a, void rational_function_1d::update(const vec& in_a,
const vec& in_b) const vec& in_b)
{ {
a.resize(in_a.size()) ; // Get the size of the input vector
b.resize(in_b.size()) ; const int np = in_a.size();
const int nq = in_b.size();
#define NORMALIZE #define NORMALIZE
...@@ -107,22 +104,45 @@ void rational_function_1d::update(const vec& in_a, ...@@ -107,22 +104,45 @@ void rational_function_1d::update(const vec& in_a,
const double b0 = 1.0; const double b0 = 1.0;
#endif #endif
for(int i=0; i<a.size(); ++i) { a[i] = in_a[i] / b0; } for(int k=0; k<np; ++k)
for(int i=0; i<b.size(); ++i) { b[i] = in_b[i] / b0; } {
_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(); // It is not possible to resize the multi-dimensional vector
const int old_nq = b.size(); // 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 // Resize the numerator
a.resize(np); if(_p_coeffs.size() != np)
b.resize(nq); {
_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 // Resize the denominator
for(int i=old_np; i<np; ++i) { a[i] = 0.0; } if(_q_coeffs.size() != nq)
for(int i=old_nq; i<nq; ++i) { b[i] = 0.0; } {
_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) ...@@ -130,40 +150,26 @@ void rational_function_1d::resize(int np, int nq)
// Get the p_i and q_j function // Get the p_i and q_j function
vec rational_function_1d::p(const vec& x) const vec rational_function_1d::p(const vec& x) const
{ {
vec res(_nY) ; vec res(1) ;
unsigned int const np = a.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
double p = 0.0f ;
const unsigned int np = _p_coeffs.size() ;
for(unsigned int i=0; i<np; ++i) for(unsigned int i=0; i<np; ++i)
{ {
p += a[k*_nY + i]*this->p(x, i) ; res[0] += _p_coeffs[i].a * this->p(x, i) ;
} }
res[k] = p ;
}
return res ; return res ;
} }
vec rational_function_1d::q(const vec& x) const vec rational_function_1d::q(const vec& x) const
{ {
vec res(_nY) ; vec res(1) ;
unsigned int const nq = b.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
double q = 0.0f ;
for(unsigned int i=0; i<nq; ++i) const unsigned int np = _q_coeffs.size() ;
for(unsigned int i=0; i<np; ++i)
{ {
q += b[k*_nY + i]*this->q(x, i) ; res[0] += _q_coeffs[i].a * this->p(x, i) ;
} }
res[k] = q ;
}
return res ; return res ;
} }
...@@ -287,19 +293,25 @@ std::vector<int> rational_function_1d::index2degree(int i) const ...@@ -287,19 +293,25 @@ std::vector<int> rational_function_1d::index2degree(int i) const
// Get the p_i and q_j function // Get the p_i and q_j function
double rational_function_1d::p(const vec& x, int i) const double rational_function_1d::p(const vec& x, int i) const
{ {
std::vector<int> deg = index2degree(i);
double res = 1.0; double res = 1.0;
for(int k=0; k<dimX(); ++k) for(int k=0; k<dimX(); ++k)
{ {
res *= pow(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]) ; const double xp = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
//res *= pow(x[k], deg[k]) ; res *= pow(xp, _p_coeffs[i].deg[k]) ;
} }
return res ; return res ;
} }
double rational_function_1d::q(const vec& x, int i) const 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 // Overload the function operator
...@@ -307,20 +319,20 @@ vec rational_function_1d::value(const vec& x) const ...@@ -307,20 +319,20 @@ vec rational_function_1d::value(const vec& x) const
{ {
vec res(1) ; vec res(1) ;
unsigned int const np = a.size() / _nY ; unsigned int const np = _p_coeffs.size();
unsigned int const nq = b.size() / _nY ; unsigned int const nq = _q_coeffs.size();
double p = 0.0f ; double p = 0.0f ;
double q = 0.0f ; double q = 0.0f ;
for(unsigned int i=0; i<np; ++i) 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) 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 ; res[0] = p/q ;
...@@ -330,25 +342,28 @@ vec rational_function_1d::value(const vec& x) const ...@@ -330,25 +342,28 @@ vec rational_function_1d::value(const vec& x) const
std::ostream& operator<< (std::ostream& out, const rational_function_1d& r) 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 = [" ; std::cout << "p = [" ;
for(int i=0; i<r.a.size(); ++i) for(unsigned int i=0; i<np; ++i)
{ {
if(i != 0) if(i != 0)
{ {
std::cout << ", " ; std::cout << ", " ;
} }
std::cout << r.a[i] ; std::cout << r._p_coeffs[i].a ;
} }
std::cout << "]" << std::endl ; std::cout << "]" << std::endl ;
std::cout << "q = [" ; std::cout << "q = [" ;
for(int i=0; i<r.b.size(); ++i) for(unsigned int i=0; i<nq; ++i)
{ {
if(i != 0) if(i != 0)
{ {
std::cout << ", " ; std::cout << ", " ;
} }
std::cout << r.b[i] ; std::cout << r._q_coeffs[i].a ;
} }
std::cout << "]" << std::endl ; std::cout << "]" << std::endl ;
...@@ -384,9 +399,7 @@ rational_function_1d* rational_function::get(int i) ...@@ -384,9 +399,7 @@ rational_function_1d* rational_function::get(int i)
{ {
if(rs[i] == NULL) if(rs[i] == NULL)
{ {
rs[i] = new rational_function_1d(np, nq); rs[i] = new rational_function_1d(dimX(), np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
// Test if the input domain is not empty. If one dimension of // Test if the input domain is not empty. If one dimension of
// the input domain is a point, I manually inflate this dimension // the input domain is a point, I manually inflate this dimension
......
...@@ -21,18 +21,18 @@ class rational_function_1d : public function ...@@ -21,18 +21,18 @@ class rational_function_1d : public function
public: // methods public: // methods
rational_function_1d() ; rational_function_1d() ;
rational_function_1d(int np, int nq, bool separable = false) ; rational_function_1d(int nX, unsigned int np, unsigned int nq,
rational_function_1d(const vec& a, const vec& b) ; bool separable = false) ;
virtual ~rational_function_1d() {} 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 value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; } virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files //! IO function to text files
virtual bool load(std::istream& in); virtual bool load(std::istream& in);
//! \brief Save the rational function expansion. It should //! \brief Save the rational function expansion. It should
...@@ -50,6 +50,7 @@ class rational_function_1d : public function ...@@ -50,6 +50,7 @@ class rational_function_1d : public function
//! algorithm to efficiently evaluate recursively defined //! algorithm to efficiently evaluate recursively defined
//! polynomials. //! polynomials.
virtual vec p(const vec& x) const ; virtual vec p(const vec& x) const ;
//! Evaluate the denominator \f$q(\mathbf{x})\f$ of the rational //! Evaluate the denominator \f$q(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast //! function. This function is provided to allow fast
//! implementation. For example one can use the Clenshaw //! implementation. For example one can use the Clenshaw
...@@ -64,19 +65,41 @@ class rational_function_1d : public function ...@@ -64,19 +65,41 @@ class rational_function_1d : public function
//! denominator of the rational function. //! denominator of the rational function.
virtual double q(const vec& x, int j) const ; 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, virtual void update(const vec& in_a,
const vec& in_b) ; const vec& in_b) ;
// Resize the polynomial //! Resize the polynomial.
virtual void resize(int np, int nq); virtual void resize(unsigned int np, unsigned int nq);
//! Get the i-th coefficient of the numerator.
virtual double getP(int i) const { return _p_coeffs[i].a; }
// Get the coefficients //! Get the i-th coefficient of the denominator.
virtual double getP(int i) const { return a[i]; } virtual double getQ(int i) const { return _q_coeffs[i].a; }
virtual double getQ(int i) const { return b[i]; }
//! 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;
}
virtual vec getP() const { return a; }
virtual vec getQ() const { return b; }
// STL stream ouput // STL stream ouput
friend std::ostream& operator<< (std::ostream& out, friend std::ostream& operator<< (std::ostream& out,
...@@ -96,9 +119,23 @@ class rational_function_1d : public function ...@@ -96,9 +119,23 @@ class rational_function_1d : public function
protected: // data protected: // data
// Store the coefficients for the moment, I assume // Structure to store a multi-dimensional coefficient. The
// the functions to be polynomials. // coeffcient, a, is associated with the vector of degree
vec a, b ; // 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? //! Is the function separable with respect to its input dimensions?
//! \todo Make possible to have only part of the dimensions //! \todo Make possible to have only part of the dimensions
......
...@@ -27,14 +27,8 @@ rational_function_chebychev_1d::rational_function_chebychev_1d() ...@@ -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_chebychev_1d::rational_function_chebychev_1d(int nX, int np, int nq) :
rational_function_1d(np, nq) rational_function_1d(nX, np, nq)
{
}
rational_function_chebychev_1d::rational_function_chebychev_1d(const vec& a,
const vec& b) :
rational_function_1d(a, b)
{ {
} }
...@@ -92,9 +86,7 @@ rational_function_1d* rational_function_chebychev::get(int i) ...@@ -92,9 +86,7 @@ rational_function_1d* rational_function_chebychev::get(int i)
{ {
if(rs[i] == NULL) if(rs[i] == NULL)
{ {
rs[i] = new rational_function_chebychev_1d(np, nq); rs[i] = new rational_function_chebychev_1d(dimX(), np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
// Test if the input domain is not empty. If one dimension of // Test if the input domain is not empty. If one dimension of
// the input domain is a point, I manually inflate this dimension // the input domain is a point, I manually inflate this dimension
......
...@@ -17,8 +17,7 @@ class rational_function_chebychev_1d : public rational_function_1d ...@@ -17,8 +17,7 @@ class rational_function_chebychev_1d : public rational_function_1d
public: // methods public: // methods
rational_function_chebychev_1d() ; rational_function_chebychev_1d() ;