Updating the interface for rational function

sources/core/function.cpp

@@ -168,8 +168,9 @@ double function::Linf_distance(const ptr<data>& d) const
 	vec mean = vec::Zero(dimY());
 	vec var  = vec::Zero(dimY());
 
-
+#ifdef DEBUG
 	std::cout << "<<DEBUG>> input param here = " << params::get_name(input_parametrization()) << std::endl;
+#endif
 
 	double linf_dist = 0.0;
 	for(int i=0; i<d->size(); ++i)

sources/core/rational_function.cpp

@@ -114,6 +114,25 @@ void rational_function_1d::update(const vec& in_a,
 		_q_coeffs[k].a = in_b[k] / b0;
 	}
 }
+		
+void rational_function_1d::update(const rational_function_1d* r)
+{
+	// Get the size of the input vector
+	const unsigned int np = r->_p_coeffs.size();
+	const unsigned int nq = r->_q_coeffs.size();
+
+	// Resize the rational function
+	resize(np, nq);
+
+	for(unsigned int k=0; k<np; ++k) {
+		_p_coeffs[k].a = r->_p_coeffs[k].a;
+	}
+
+	for(unsigned int k=0; k<nq; ++k) {
+		_q_coeffs[k].a = r->_q_coeffs[k].a;
+	}
+
+}
 
 void rational_function_1d::resize(unsigned int np, unsigned int nq)
 {
@@ -388,7 +407,7 @@ rational_function::~rational_function()
 
 void rational_function::update(int i, rational_function_1d* r)
 {
-	rs[i] = r;
+	rs[i]->update(r);
 }
 
 

sources/core/rational_function.h

@@ -71,6 +71,11 @@ class rational_function_1d : public function
 		//! coefficients.
 		virtual void update(const vec& in_a, 
 		                    const vec& in_b) ;
+		
+		//! Update the 1D rational function with another one. Note that
+		//! this function can change the dimensions of the coefficients
+		//! vectors.
+		virtual void update(const rational_function_1d* r) ;
 
 		//! Resize the polynomial.
 		virtual void resize(unsigned int np, unsigned int nq);
@@ -159,6 +164,16 @@ class rational_function : public function
 		virtual bool load(std::istream& in) ;
 
 		// Update the function
+		virtual void update(const ptr<rational_function>& r) 
+		{
+			assert(r->dimX() == dimX());
+			assert(r->dimY() == dimY());
+
+			for(int k=0; k<dimY(); ++k) 
+			{
+				get(k)->update(r->get(k));
+			}
+		}
 		virtual void update(int i, rational_function_1d* r) ;
 
 		//! Get the 1D function associated with color channel i. If no one exist, 

sources/plugins/plugins.pro

@@ -11,7 +11,6 @@ SUBDIRS  = \
 	rational_fitter_parsec_one					\
 	rational_fitter_parsec_multi					\
 	rational_function_chebychev					\
-	rational_function_chebychev_opt					\
 	rational_function_legendre					\
 	rational_function_cosine					\
 	nonlinear_fitter_eigen                 				\

sources/plugins/rational_fitter_parallel/quadratic_program.h

@@ -185,11 +185,11 @@ class quadratic_program
                 data->get(*it, x, yl, yu);
 
 				vec y = r->value(x);
-				bool fail_upper = y[ny] > yu[ny];
-				bool fail_lower = y[ny] < yl[ny];
+				bool fail_upper = y[0] > yu[ny];
+				bool fail_lower = y[0] < yl[ny];
 				if(fail_lower || fail_upper)
 				{
-                    const double dev = std::abs(0.5*(yu[ny]+yl[ny]) - y[ny]);
+                    const double dev = std::abs(0.5*(yu[ny]+yl[ny]) - y[0]);
 
 					nb_failed++;
 
@@ -303,7 +303,7 @@ int quadratic_program::next_unmatching_constraint(int i, int ny, const rational_
 		data->get(n, x, yl, yu);
 
 		vec y = r->value(x);
-        if(y[ny] < yl[ny] || y[ny] > yu[ny])
+        if(y[0] < yl[ny] || y[0] > yu[ny])
 		{
 			return n;
 		}

sources/plugins/rational_fitter_parallel/rational_fitter.cpp

@@ -106,7 +106,7 @@ bool rational_fitter_parallel::fit_data(const ptr<data>& dat, ptr<function>& fit
 		int nb_sol_found  = 0;
 		int nb_sol_tested = 0;
 
-        #pragma omp parallel for shared(args, nb_sol_found, nb_sol_tested, min_delta, mean_delta), schedule(dynamic,1)
+        #pragma omp parallel for shared(r, args, nb_sol_found, nb_sol_tested, min_delta, mean_delta), schedule(dynamic,1)
         for(int j=1; j<i; ++j)
         {
             // Compute the number of coefficients in the numerator and in the denominator
@@ -115,7 +115,7 @@ bool rational_fitter_parallel::fit_data(const ptr<data>& dat, ptr<function>& fit
             int temp_np = i - j;
             int temp_nq = j;
 
-            vec p(temp_np*r->dimY()), q(temp_nq*r->dimY());
+            //vec p(temp_np*r->dimY()), q(temp_nq*r->dimY());
 
             // Allocate a rational function and set it to the correct size, dimensions
             // and parametrizations.
@@ -144,12 +144,13 @@ bool rational_fitter_parallel::fit_data(const ptr<data>& dat, ptr<function>& fit
             bool is_fitted = fit_data(d, temp_np, temp_nq, rk, args, delta, linf_dist, l2_dist);
             if(is_fitted)
             {
-                #pragma omp critical (nb_sol_found)
+                #pragma omp critical (r)
                 {
                     ++nb_sol_found ;
                     mean_delta += delta ;
 
-                    std::cout << "<<INFO>> found a solution with np=" << temp_np << ", nq = " << temp_nq << std::endl;
+                    std::cout << "<<INFO>> found a solution with np=" << temp_np 
+							         << ", nq = " << temp_nq << std::endl;
                     std::cout << "<<INFO>> Linf error = " << linf_dist << std::endl;
                     std::cout << "<<INFO>> L2   error = " << l2_dist << std::endl;
                     std::cout << "<<INFO>>      delta = " << delta << std::endl;
@@ -162,10 +163,7 @@ bool rational_fitter_parallel::fit_data(const ptr<data>& dat, ptr<function>& fit
                         min_delta   = delta ;
                         min_l2_dist = l2_dist ;
                         r->setSize(temp_np, temp_nq);
-                        for(int y=0; y<r->dimY(); ++y)
-                        {
-                            r->update(y, rk->get(y));
-                        }
+								r->update(rk);
                     }
                 }
             }

sources/plugins/rational_function_chebychev_opt/SConscript

-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)

sources/plugins/rational_function_chebychev_opt/rational_function.cpp

-#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 nX, int np, int nq) :
-	rational_function_1d(nX, np, nq)
-{
-	this->resize(np, nq);
-}
-
-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(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
-			// 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;
-	}
-}

sources/plugins/rational_function_chebychev_opt/rational_function.h

-#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 nX, int np, int nq) ;
-    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 = _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 += _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)
-	{
-		// Get the size of the input vector
-		const int np = in_a.size();
-		const int nq = in_b.size();
-
-		for(int k=0; k<np; ++k)
-		{
-			_p_coeffs[k].a = in_a[k];
-		}
-		
-		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
-			}
-		}
-} ;
-

sources/plugins/rational_function_chebychev_opt/rational_function_chebychev_opt.pro

-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
-