Attention une mise à jour du serveur va être effectuée le vendredi 16 avril entre 12h et 12h30. Cette mise à jour va générer une interruption du service de quelques minutes.

Commit 709fc8f3 authored by Laurent Belcour's avatar Laurent Belcour

Updating the test suite to include mutidimensional tests.

The rational function class is being refactorized to handle better multidimensional
functions. The base class is now a template that takes 1D RF as input. To use it
I still need to remove the dependancy to Qt since we cannot use templates with QObjects.
parent 9fb2259a
......@@ -26,3 +26,4 @@ SOURCES = plugins_manager.cpp \
params.cpp \
function.cpp \
# clustering.cpp
rational_function.inl
......@@ -151,18 +151,18 @@ void params::from_cartesian(const double* invec, params::input outtype,
// 2D Parametrizations
case params::COS_TH_TD:
outvec[0] = half[2];
outvec[1] = half[0]*outvec[0] + half[1]*outvec[1] + half[2]*outvec[2];
outvec[1] = half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2];
break;
case params::RUSIN_TH_TD:
outvec[0] = acos(half[2]);
outvec[2] = acos(half[0]*outvec[0] + half[1]*outvec[1] + half[2]*outvec[2]);
outvec[2] = acos(half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2]);
break;
// 3D Parametrization
case params::RUSIN_TH_PH_TD:
outvec[0] = acos(half[2]);
outvec[1] = atan2(half[0], half[1]);
outvec[2] = acos(half[0]*outvec[0] + half[1]*outvec[1] + half[2]*outvec[2]);
outvec[2] = acos(half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2]);
break;
case params::RUSIN_TH_TD_PD:
outvec[0] = acos(half[2]);
......
......@@ -8,61 +8,320 @@
#include <algorithm>
#include <cmath>
rational_function::rational_function() : a(), b()
#ifdef NEW
rational_function_1d::rational_function_1d()
{
}
rational_function::rational_function(int np, int nq) : a(np), b(nq)
rational_function_1d::rational_function_1d(int np, int nq)
{
a.resize(np);
b.resize(nq);
}
rational_function::rational_function(const std::vector<double>& a,
const std::vector<double>& b) :
rational_function_1d::rational_function_1d(const std::vector<double>& a,
const std::vector<double>& b) :
a(a), b(b)
{
}
rational_function::~rational_function()
void rational_function_1d::load(const std::string& filename)
{
}
void rational_function::update(const std::vector<double>& in_a,
const std::vector<double>& in_b)
void rational_function_1d::update(const std::vector<double>& in_a,
const std::vector<double>& in_b)
{
a.reserve(in_a.size()) ;
b.reserve(in_b.size()) ;
a = in_a ;
b = in_b ;
}
void rational_function_1d::resize(int np, int nq)
{
const int old_np = a.size();
const int old_nq = b.size();
// Overload the function operator
vec rational_function::value(const vec& x) const
// Resize the vector
a.resize(np);
b.resize(nq);
// 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; }
}
// 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 ;
unsigned int const nq = b.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[k*_nY + i]*this->p(x, i) ;
}
res[k] = p ;
}
return res ;
}
vec rational_function_1d::q(const vec& x) const
{
vec res(_nY) ;
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)
{
q += b[k*_nY + i]*this->q(x, i) ;
}
res[k] = p/q ;
res[k] = q ;
}
return res ;
}
// Estimate the number of configuration for an indice
// vector of dimension d with maximum element value
// being k.
int rational_function_1d::estimate_dk(int k, int d)
{
if(d == 1)
{
return 1;
}
else if(d ==2)
{
return k+1;
}
else
{
int res = 0;
for(int i=0; i<=k; ++i)
{
res += estimate_dk(k-i, d-1);
}
return res;
}
}
// Populate a vector of degrees of dimension N using a
// maximum degree of M. The index at the current level
// is j
void rational_function_1d::populate(std::vector<int>& vec, int N, int M, int j)
{
// For each dimension, estimate the current level
// based on the number of configurations in the
// other dimensions
int current_M = M ;
int nb_conf = 0;
for(int d=0; d<N-1; ++d)
{
int k;
for(k=0; k<=current_M; ++k)
{
int oracle = estimate_dk(current_M-k, N-(d+1));
if(nb_conf <= j && j < nb_conf+oracle)
{
break;
}
nb_conf += oracle;
}
vec[N-1 - d] = k ;
current_M -= k ;
}
vec[0] = current_M;
}
std::vector<int> rational_function_1d::index2degree(int i) const
{
std::vector<int> deg ; deg.assign(dimX(), 0) ;
if(i == 0)
return deg ;
if(dimX() == 1)
{
deg[0] = i;
}
else if(dimX() == 2)
{
int Nk = 1 ;
int k = 1 ;
while(!(i >= Nk && i < Nk+k+1))
{
Nk += k+1 ;
++k ;
}
int r = i-Nk ;
deg[0] = k-r;
deg[1] = r;
}
else
{
int Nk = 1 ;
int k = 1 ;
int dk = estimate_dk(k, dimX()) ;
while(!(i >= Nk && i < Nk+dk))
{
Nk += dk ;
++k ;
dk = estimate_dk(k, dimX()) ;
}
// Populate the vector from front to back
int j = i-Nk ;
populate(deg, dimX(), k, j) ;
}
return deg ;
}
double legendre(double x, int i)
{
if(i == 0)
{
return 1;
}
else if(i == 1)
{
return x;
}
else
{
return ((2*i-1)*x*legendre(x, i-1) - (i-1)*legendre(x, i-2)) / (double)i ;
}
}
//#define POLYNOMIALS
// 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)
{
#ifdef POLYNOMIALS
res *= pow(x[k], deg[k]) ;
#else // LEGENDRE
res *= legendre(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]);
#endif
}
return res ;
}
double rational_function_1d::q(const vec& x, int i) const
{
std::vector<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
#ifdef POLYNOMIALS
res *= pow(x[k], deg[k]) ;
#else // LEGENDRE
res *= legendre(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]);
#endif
}
return res ;
}
// Overload the function operator
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 ;
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[i]*this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += b[i]*this->q(x, i) ;
}
res[0] = p/q ;
return res ;
}
std::ostream& operator<< (std::ostream& out, const rational_function_1d& r)
{
std::cout << "p = [" ;
for(unsigned int i=0; i<r.a.size(); ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.a[i] ;
}
std::cout << "]" << std::endl ;
std::cout << "q = [" ;
for(unsigned int i=0; i<r.b.size(); ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.b[i] ;
}
std::cout << "]" << std::endl ;
return out ;
}
#else
rational_function::rational_function() : a(), b()
{
}
rational_function::rational_function(int np, int nq) : a(np), b(nq)
{
}
rational_function::rational_function(const std::vector<double>& a,
const std::vector<double>& b) :
a(a), b(b)
{
}
void rational_function::update(const std::vector<double>& in_a,
const std::vector<double>& in_b)
{
a.reserve(in_a.size()) ;
b.reserve(in_b.size()) ;
a = in_a ;
b = in_b ;
}
// Get the p_i and q_j function
vec rational_function::p(const vec& x) const
{
......@@ -252,6 +511,34 @@ double rational_function::q(const vec& x, int i) const
return res ;
}
// Overload the function operator
vec rational_function::value(const vec& x) const
{
vec res(_nY) ;
unsigned int const np = a.size() / _nY ;
unsigned int const nq = b.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[k*_nY + i]*this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += b[k*_nY + i]*this->q(x, i) ;
}
res[k] = p/q ;
}
return res ;
}
// IO function to text files
void rational_function::load(const std::string& filename)
{
......@@ -635,5 +922,11 @@ std::ostream& operator<< (std::ostream& out, const rational_function& r)
return out ;
}
#endif
//Q_EXPORT_PLUGIN2(rational_function, rational_function)
......@@ -3,6 +3,7 @@
// Include STL
#include <vector>
#include <string>
#include <sstream>
// Interface
#include <QObject>
......@@ -12,6 +13,142 @@
#include "args.h"
#include "common.h"
//! \todo template this
#ifdef NEW
class rational_function_1d : public function
{
public: // methods
rational_function_1d() ;
rational_function_1d(int np, int nq) ;
rational_function_1d(const std::vector<double>& a, const std::vector<double>& b) ;
virtual ~rational_function_1d() {}
// Overload the function operator
virtual vec value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; }
// Get the numerator (p) and denominator (q) functions
virtual vec p(const vec& x) const ;
virtual vec q(const vec& x) const ;
// 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 ;
// IO function to text files
virtual void load(const std::string& filename) ;
// Update the function
virtual void update(const std::vector<double>& in_a,
const std::vector<double>& in_b) ;
// Resize the polynomial
virtual void resize(int np, int nq);
// Get the coefficients
virtual double getP(int i) const { return a[i]; }
virtual double getQ(int i) const { return b[i]; }
virtual std::vector<double> getP() const { return a; }
virtual std::vector<double> getQ() const { return b; }
// STL stream ouput
friend std::ostream& operator<< (std::ostream& out, const rational_function_1d& r) ;
protected: // functions
//! Convert a 1D index into a vector of degree for a
//! multinomial coeffcient. The resulting vector v should
//! be used as prod_k x[k]^v[k] for the monomial basis
std::vector<int> index2degree(int i) const ;
static int estimate_dk(int k, int d);
static void populate(std::vector<int>& vec, int N, int M, int j);
protected: // data
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
std::vector<double> a ;
std::vector<double> b ;
} ;
// Prior definition for the standard output declaration
template<class RF1D> class rational_function_t ;
template<class RF1D> std::ostream& operator<< (std::ostream& out, rational_function_t<RF1D>& r) ;
template<class RF1D> class rational_function_t : public QObject, public function
{
Q_OBJECT
Q_INTERFACES(function)
public: // methods
rational_function_t() ;
rational_function_t(int np, int nq) ;
rational_function_t(const std::vector<double>& a, const std::vector<double>& b) ;
virtual ~rational_function_t() ;
// 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 void load(const std::string& filename) ;
// Update the function
virtual void update(int i, RF1D* 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.
virtual RF1D* get(int i) ;
// STL stream ouput
friend std::ostream& operator<< <> (std::ostream& out, rational_function_t<RF1D>& r) ;
//! \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 ;
//! 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);
}
protected: // functions
//! \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 ;
protected: // data
// Store the y \in R rational functions. Each channel is a distinct polynomial
// and should be fitted separately.
std::vector<RF1D*> rs ;
// Size of the polynomials
//! \todo Change it by a more adaptive scheme, with different np, nq per color
//!channel?
int np, nq;
} ;
#include "rational_function.inl"
typedef rational_function_t<rational_function_1d> rational_function;
#else
class rational_function : public QObject, public function
{
Q_OBJECT
......@@ -22,7 +159,7 @@ class rational_function : public QObject, public function
rational_function() ;
rational_function(int np, int nq) ;
rational_function(const std::vector<double>& a, const std::vector<double>& b) ;
virtual ~rational_function() ;
virtual ~rational_function() { } ;
// Overload the function operator
virtual vec value(const vec& x) const ;
......@@ -83,3 +220,4 @@ class rational_function : public QObject, public function
std::vector<double> a ;
std::vector<double> b ;
} ;
#endif
template <class RF1D>
rational_function_t<RF1D>::rational_function_t() : np(0), nq(0)
{
}
template <class RF1D>
rational_function_t<RF1D>::rational_function_t(int np, int nq) : np(np), nq(nq)
{
}
//! \todo clean memory here
template <class RF1D>
rational_function_t<RF1D>::~rational_function_t()
{
}
template <class RF1D>
void rational_function_t<RF1D>::update(int i,
RF1D* r)
{
rs[i] = r;
}
template <class RF1D>
RF1D* rational_function_t<RF1D>::get(int i)
{
// Check for consistency in the index of color channel
if(i < _nY)
{
if(rs[i] == NULL)