Commit 9e836626 authored by Laurent Belcour's avatar Laurent Belcour
Browse files

Adding a first attempt to optimize the rational functions

parent e23239c2
......@@ -44,12 +44,24 @@ class rational_function_1d : public function
/* RATIONAL FUNCTION SPECIFIC */
// Get the numerator (p) and denominator (q) functions
//! Evaluate the numerator \f$p(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast
//! implementation. For example one can use the Clenshaw
//! 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
//! algorithm to efficiently evaluate recursively defined
//! polynomials.
virtual vec q(const vec& x) const ;
// Get the p_i and q_j function
//! Evaluate the basis function \f$p_i(\mathbf{x})\f$ for the
//! numerator of the rational function.
virtual double p(const vec& x, int i) const ;
//! Evaluate the basis function \f$q_i(\mathbf{x})\f$ for the
//! denominator of the rational function.
virtual double q(const vec& x, int j) const ;
// Update the function
......
......@@ -26,6 +26,7 @@ SConscript('rational_fitter_parallel/SConscript')
# Building rational functions
SConscript('rational_function_legendre/SConscript')
SConscript('rational_function_chebychev/SConscript')
SConscript('rational_function_chebychev_opt/SConscript')
SConscript('rational_function_cosine/SConscript')
# Building retro functions
......
env = Environment()
env.Append(CPPPATH = ['../../../external/build/include', '../../'])
env.Append(LIBPATH = ['../../../external/build/lib', '../../build'])
sources = ['rational_function.cpp']
libs = ['core']
env.SharedLibrary('../../build/rational_function_chebychev_opt', sources, LIBS=libs)
#include "rational_function.h"
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <core/common.h>
ALTA_DLL_EXPORT function* provide_function()
{
return new rational_function_chebychev();
}
rational_function_chebychev::rational_function_chebychev() : rational_function()
{
}
rational_function_chebychev::~rational_function_chebychev()
{
}
rational_function_chebychev_1d::rational_function_chebychev_1d()
{
setDimX(0);
setDimY(0);
}
rational_function_chebychev_1d::rational_function_chebychev_1d(int np, int nq) :
rational_function_1d(np, nq)
{
setDimX(0);
setDimY(0);
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
if(i == 0)
{
return 1;
}
else if(i == 1)
{
return x;
}
else
{
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
double rational_function_chebychev_1d::p(const vec& x, int i) const
{
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
double xk = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= chebychev(xk, _p_coeffs[i].deg[k]);
}
return res ;
}
double rational_function_chebychev_1d::q(const vec& x, int i) const
{
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
double xk = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= chebychev(xk, _q_coeffs[i].deg[k]);
}
return res ;
}
rational_function_1d* rational_function_chebychev::get(int i)
{
// Check for consistency in the index of color channel
if(i < _nY)
{
if(rs[i] == NULL)
{
rs[i] = new rational_function_chebychev_1d(np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
// Test if the input domain is not empty. If one dimension of
// the input domain is a point, I manually inflate this dimension
// to avoid numerical issues.
vec _min = min();
vec _max = max();
for(int k=0; k<dimX(); ++k)
{
if(_min[k] == _max[k])
{
_min[k] -= 1.0;
_max[k] += 1.0;
}
}
rs[i]->setMin(_min) ;
rs[i]->setMax(_max) ;
}
return rs[i];
}
else
{
std::cout << "<<ERROR>> tried to access out of bound 1D RF" << std::endl;
return NULL;
}
}
#pragma once
// Include STL
#include <vector>
#include <string>
// Interface
#include <core/function.h>
#include <core/rational_function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
#include <core/common.h>
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) ;
virtual ~rational_function_chebychev_1d() {}
// 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 ;
// Overload the function operator
virtual vec 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 += _p_coeffs[i].a * this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += _q_coeffs[i].a * this->q(x, i) ;
}
res[0] = p/q ;
return res ;
}
// Update the function
virtual void update(const vec& in_a,
const vec& in_b)
{
rational_function_1d::update(in_a, in_b);
const int np = in_a.size();
for(int k=0; k<np; ++k)
{
_p_coeffs[k].a = in_a[k];
}
const int nq = in_b.size();
for(int k=0; k<nq; ++k)
{
_q_coeffs[k].a = in_b[k];
}
}
// Resize the polynomial
virtual void resize(int np, int nq)
{
if(dimX() == 0) return;
if(_p_coeffs.size() != np)
{
_p_coeffs.resize(np);
for(int k=0; k<np; ++k)
{
std::vector<int> deg = index2degree(k);
_p_coeffs[k].deg = deg;
}
}
if(_q_coeffs.size() != nq)
{
_q_coeffs.resize(nq);
for(int k=0; k<nq; ++k)
{
std::vector<int> deg = index2degree(k);
_q_coeffs[k].deg = deg;
}
}
}
// Get the coefficients
virtual double getP(int i) const { return _p_coeffs[i].a; }
virtual double getQ(int i) const { return _q_coeffs[i].a; }
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;
}
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;
}
protected: // data
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;
} ;
class rational_function_chebychev : public rational_function
{
public: // methods
rational_function_chebychev() ;
virtual ~rational_function_chebychev() ;
//! 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) ;
//! Update the y-1D function for the ith dimension.
//! \note It will test if the 1D function provided is of the dynamic type
//! \name rational_function_chebychev_1d
virtual void update(int i, rational_function_1d* r)
{
if(dynamic_cast<rational_function_chebychev_1d*>(r) != NULL)
{
rational_function::update(i, r);
}
else
{
#ifdef DEBUG
std::cerr << "<<ERROR>> the function provided is not of type \"rational_function_chebychev\"" << std::endl;
#endif
}
}
} ;
TEMPLATE = lib
CONFIG *= qt \
plugin \
eigen
DESTDIR = ../../build
INCLUDEPATH += ../..
HEADERS = rational_function.h
SOURCES = rational_function.cpp
LIBS += -L../../build \
-lcore
#QMAKE_CXXFLAGS += -frounding-math \
# -fPIC \
# -g
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