Multidimensional fitting on the output space

sources/core/common.h

@@ -2,26 +2,31 @@
 
 #include <vector>
 #include <iostream>
+#include <cassert>
 
 class vec : public std::vector<double>
 {
 	public:
-/*		// Constructor & Destructors
+//*		// Constructor & Destructors
 		//
-		vec(int dim)
+		vec() : std::vector<double>()
+		{
+		}
+		vec(int dim) : std::vector<double>(dim)
 		{
 			assign(dim, 0.0) ;
 		} ;
 		virtual ~vec() 
 		{
 		} ;
-
+//*/
+		
 		// Mathematical operators
 		//
 		friend vec operator-(const vec& a)
 		{
-			vec b(a.size) ;
-			for(int i=0; i<a.size(); ++i)
+			vec b(a.size()) ;
+			for(unsigned int i=0; i<a.size(); ++i)
 			{
 				b[i] = -a[i] ;
 			}
@@ -32,8 +37,8 @@ class vec : public std::vector<double>
 #ifdef DEBUG
 			assert(a.size() == b.size()) ;
 #endif
-			vec c(a.size) ;
-			for(int i=0; i<a.size(); ++i)
+			vec c(a.size()) ;
+			for(unsigned int i=0; i<a.size(); ++i)
 			{
 				c[i] = a[i] - b[i];
 			}
@@ -44,8 +49,8 @@ class vec : public std::vector<double>
 #ifdef DEBUG
 			assert(a.size() == b.size()) ;
 #endif
-			vec c(a.size) ;
-			for(int i=0; i<a.size(); ++i)
+			vec c(a.size()) ;
+			for(unsigned int i=0; i<a.size(); ++i)
 			{
 				c[i] = a[i] + b[i];
 			}
@@ -56,14 +61,14 @@ class vec : public std::vector<double>
 #ifdef DEBUG
 			assert(a.size() == b.size()) ;
 #endif
-			vec c(a.size) ;
-			for(int i=0; i<a.size(); ++i)
+			vec c(a.size()) ;
+			for(unsigned int i=0; i<a.size(); ++i)
 			{
 				c[i] = a[i] * b[i];
 			}
 			return c ;
 		}
-
+/*
 		friend double norm(const vec& a)
 		{
 

sources/core/fitter.h

@@ -19,7 +19,7 @@ class fitter
 		// underling function class. Return the best
 		// fit (along with fitting information ?)
 		//
-		virtual bool fit_data(const data* d, function*& f) = 0 ;
+		virtual bool fit_data(const data* d, function* f) = 0 ;
 
 
 		virtual void set_parameters(const arguments& args) = 0 ;

sources/plugins/rational_fitter_cgal/rational_fitter_cgal.cpp

@@ -6,7 +6,6 @@
 #include <CGAL/MP_Float.h>
 #include <Eigen/SVD>
 
-#include <boost/regex.hpp>
 #include <string>
 #include <iostream>
 #include <fstream>
@@ -14,9 +13,6 @@
 #include <algorithm>
 #include <cmath>
 
-#include <rational_function.h>
-#include <rational_data.h>
-
 typedef CGAL::MP_Float ET ;
 typedef CGAL::Quadratic_program<ET> Program ;
 typedef CGAL::Quadratic_program_solution<ET> Solution ;
@@ -29,18 +25,46 @@ rational_fitter_cgal::~rational_fitter_cgal()
 {
 }
 
-bool rational_fitter_cgal::fit_data(const data* d, function*& fit)
+bool rational_fitter_cgal::fit_data(const data* dat, function* fit)
 {
-	std::cout << "<<INFO>> np in  [" << _min_np << ", " << _max_np << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
+	rational_function* r = dynamic_cast<rational_function*>(fit) ;
+	const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
+	if(r == NULL || d == NULL)
+	{
+		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
+		return false ;
+	}
+
+	// I need to set the dimension of the resulting function to be equal
+	// to the dimension of my fitting problem
+	r->setDimX(d->dimX()) ;
+	r->setDimY(d->dimY()) ;
+
+#ifdef DEBUG_DEGREES
+	for(int i=0; i<100; ++i)
+	{
+		std::vector<int> v = r->index2degree(i) ;
+		std::cout << i << " = " ;
+		for(int j=0; j<v.size(); ++j)
+			std::cout << v[j] << "\t" ;
+		std::cout << std::endl ;
+	}
+	throw ;
+#endif
+
+	std::cout << "<<INFO>> np in  [" << _min_np << ", " << _max_np 
+	          << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
+
 	int temp_np = _min_np, temp_nq = _min_nq ;
 	while(temp_np <= _max_np || temp_nq <= _max_nq)
 	{
-		if(fit_data(d, temp_np, temp_nq, fit))
+		if(fit_data(d, temp_np, temp_nq, r))
 		{
+			std::cout << "<<INFO>> got a fit using np = " << temp_np << " & nq =  " << temp_nq << std::endl ;
 			return true ;
 		}
 
-		std::cout << "<<INFO>> fitt using np = " << temp_np << " & nq =  " << temp_nq << " failed\r"  ;
+		std::cout << "<<INFO>> fit using np = " << temp_np << " & nq =  " << temp_nq << " failed\r"  ;
 		std::cout.flush() ;
 
 		if(temp_np < _max_np)
@@ -64,25 +88,35 @@ void rational_fitter_cgal::set_parameters(const arguments& args)
 	_min_nq = args.get_float("min-nq", _max_nq) ;
 }
 		
-
-bool rational_fitter_cgal::fit_data(const data* dat, int np, int nq, function*& rf) 
+// TODO Finish
+bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, rational_function* r) 
 {
-	// by default, we have a nonnegative QP with Ax - b >= 0
-	Program qp (CGAL::LARGER, false, 0, false, 0) ; 
 
-	rational_function* r = dynamic_cast<rational_function*>(rf) ;
-	const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
-	if(r == NULL || d == NULL)
+	// Multidimensional coefficients
+	std::vector<double> Pn ; Pn.reserve(d->dimY()*np) ;
+	std::vector<double> Qn ; Qn.reserve(d->dimY()*nq) ;
+
+	for(int j=0; j<d->dimY(); ++j)
 	{
-		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
-		return false ;
+		if(!fit_data(d, np, nq, j, r))
+			return false ;
+
+		for(int i=0; i<np; ++i) { Pn.push_back(r->getP(i)) ; }
+		for(int i=0; i<nq; ++i) { Qn.push_back(r->getQ(i)) ; }
 	}
 
-	// I need to set the dimension of the resulting function to be equal
-	// to the dimension of my fitting problem
-	r->setDimX(d->dimX()) ;
-	r->setDimY(d->dimY()) ;
+	r->update(Pn, Qn) ;
+	return true ;
+}
 
+// dat is the data object, it contains all the points to fit
+// np and nq are the degree of the RP to fit to the data
+// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
+// the function return a ration BRDF function and a boolean
+bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, int ny, rational_function* r) 
+{
+	// by default, we have a nonnegative QP with Ax - b >= 0
+	Program qp (CGAL::LARGER, false, 0, false, 0) ; 
 
 	// Select the size of the result vector to
 	// be equal to the dimension of p + q
@@ -126,13 +160,13 @@ bool rational_fitter_cgal::fit_data(const data* dat, int np, int nq, function*&
 				d->get(i, yl, yu) ;
 
 				const double qi = r->q(d->get(i), j-np) ;
-				a0_norm += qi*qi * (yl[0]*yl[0]) ;
-				qp.set_a(j, i, ET(-yl[0] * qi)) ;
-				CI(j, i) = -yl[0] * qi ;
+				a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
+				qp.set_a(j, i, ET(-yl[ny] * qi)) ;
+				CI(j, i) = -yl[ny] * qi ;
 				
-				a1_norm += qi*qi * (yu[0]*yu[0]) ;
-				qp.set_a(j, i+d->size(),  ET(yu[0] * qi)) ;
-				CI(j, i+d->size()) = yu[0] * qi ;
+				a1_norm += qi*qi * (yu[ny]*yu[ny]) ;
+				qp.set_a(j, i+d->size(),  ET(yu[ny] * qi)) ;
+				CI(j, i+d->size()) = yu[ny] * qi ;
 			}
 		}
 	
@@ -257,26 +291,25 @@ bool rational_fitter_cgal::fit_data(const data* dat, int np, int nq, function*&
 	{
 		// Recopy the vector d
 		std::vector<double> p, q;
+		p.assign(np, 0.0) ; q.assign(nq, 0.0) ;
 		for(int i=0; i<np+nq; ++i)
 		{
 			const double v = CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
 
 			if(i < np)
 			{
-				p.push_back(v) ;
+				p[i] = v ;
 			}
 			else
 			{
-				q.push_back(v) ;
+				q[i-np] = v ;
 			}
 		}
-
-		if(r != NULL)
-		{
-			delete r ;
-		}
+/*
 		r = new rational_function(p, q);
-			
+/*/
+ 		r->update(p, q) ;
+//*/			
 		std::cout << "<<INFO>> got solution " << *r << std::endl ;
 		return true;
 	}

sources/plugins/rational_fitter_cgal/rational_fitter_cgal.h

@@ -11,6 +11,9 @@
 #include <core/fitter.h>
 #include <core/args.h>
 
+#include <rational_function.h>
+#include <rational_data.h>
+
 class rational_fitter_cgal : public QObject, public fitter
 {
 	Q_OBJECT
@@ -22,11 +25,12 @@ class rational_fitter_cgal : public QObject, public fitter
 		virtual ~rational_fitter_cgal() ;
 
 		// Fitting a data object
-		virtual bool fit_data(const data* d, function*& fit) ;
+		virtual bool fit_data(const data* d, function* fit) ;
 
 		// Fitting a data object using np elements in the numerator and nq 
 		// elements in the denominator
-		virtual bool fit_data(const data* d, int np, int nq, function*& fit) ;
+		virtual bool fit_data(const rational_data* d, int np, int nq, rational_function* fit) ;
+		virtual bool fit_data(const rational_data* dat, int np, int nq, int ny, rational_function* fit) ;
 
 		virtual void set_parameters(const arguments& args) ;
 

sources/plugins/rational_fitter_cgal/rational_fitter_cgal.pro

@@ -2,8 +2,7 @@ TEMPLATE        = lib
 CONFIG         *= qt      \
                   plugin  \
 						eigen   \
-						cgal    \
-#						debug
+						cgal    
 
 DESTDIR         = ../../build
  
@@ -14,8 +13,7 @@ INCLUDEPATH    += ../rational_function \
 HEADERS         = rational_fitter_cgal.h
 SOURCES         = rational_fitter_cgal.cpp
 
-LIBS           += -lCGAL				   \
-                  -L../build           \
+LIBS           += -L../build           \
 						-lrational_function	\
 						-lrational_data
 

sources/plugins/rational_fitter_quadprog++/rational_fitter.cpp

+#include "rational_fitter.h"
+
+#include <Eigen/SVD>
+#include <Array.hh>
+#include <QuadProg++.hh>
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <limits>
+#include <algorithm>
+#include <cmath>
+
+#include <rational_function.h>
+#include <rational_data.h>
+
+rational_fitter_quadprog++::rational_fitter_quadprog++() : QObject()
+{
+}
+rational_fitter_quadprog++::~rational_fitter_quadprog++() 
+{
+}
+
+bool rational_fitter_quadprog++::fit_data(const data* d, function*& fit)
+{
+	std::cout << "<<INFO>> np in  [" << _min_np << ", " << _max_np << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
+	int temp_np = _min_np, temp_nq = _min_nq ;
+	while(temp_np <= _max_np || temp_nq <= _max_nq)
+	{
+		if(fit_data(d, temp_np, temp_nq, fit))
+		{
+			return true ;
+		}
+
+		std::cout << "<<INFO>> fitt using np = " << temp_np << " & nq =  " << temp_nq << " failed\r"  ;
+		std::cout.flush() ;
+
+		if(temp_np < _max_np)
+		{
+			++temp_np ;
+		}
+		if(temp_nq < _max_nq)
+		{
+			++temp_nq ;
+		}
+	}
+
+	return false ;
+}
+
+void rational_fitter_quadprog++::set_parameters(const arguments& args)
+{
+	_max_np = args.get_float("np", 10) ;
+	_max_nq = args.get_float("nq", 10) ;
+	_min_np = args.get_float("min-np", _max_np) ;
+	_min_nq = args.get_float("min-nq", _max_nq) ;
+}
+		
+
+bool rational_fitter_quadprog++::fit_data(const data* dat, int np, int nq, function*& rf) 
+{
+	// by default, we have a nonnegative QP with Ax - b >= 0
+	Program qp (CGAL::LARGER, false, 0, false, 0) ; 
+
+	rational_function* r = dynamic_cast<rational_function*>(rf) ;
+	const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
+	if(r == NULL || d == NULL)
+	{
+		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
+		return false ;
+	}
+
+	// I need to set the dimension of the resulting function to be equal
+	// to the dimension of my fitting problem
+	r->setDimX(d->dimX()) ;
+	r->setDimY(d->dimY()) ;
+
+	// Matrices of the problem
+	QuadProgPP::Matrix<double> G (0.0, np+nq, np+nq) ;
+	QuadProgPP::Matrix<double> CI(0.0, np+nq, 2*d->size()) ;
+	QuadProgPP::Vector<double> ci(0.0, 2*d->size()) ;
+	QuadProgPP::Matrix<double> CE(0.0, np+nq, 2*d->size()) ;
+	QuadProgPP::Vector<double> ce(0.0, 2*d->size()) ;
+
+	// Select the size of the result vector to
+	// be equal to the dimension of p + q
+	for(int i=0; i<np+nq; ++i)
+	{
+		G.set(1.0, i, i) ; 
+	}
+	
+	// Each constraint (fitting interval or point
+	// add another dimension to the constraint
+	// matrix
+	for(int i=0; i<d->size(); ++i)	
+	{		
+		// Norm of the row vector
+		double a0_norm = 0.0 ;
+		double a1_norm = 0.0 ;
+
+		// A row of the constraint matrix has this 
+		// form: [p_{0}(x_i), .., p_{np}(x_i), -f(x_i) q_{0}(x_i), .., -f(x_i) q_{nq}(x_i)]
+		// For the lower constraint and negated for 
+		// the upper constraint
+		for(int j=0; j<np+nq; ++j)
+		{
+			// Filling the p part
+			if(j<np)
+			{
+				const double pi = r->p(d->get(i), j) ;
+				a0_norm += pi*pi ;
+				a1_norm += pi*pi ;
+				qp.set_a(j, i,  ET(pi)) ;
+				qp.set_a(j, i+d->size(), ET(-pi)) ;
+				CI(j, i) =  pi ;
+				CI(j, i+d->size()) = -pi ;
+			}
+			// Filling the q part
+			else
+			{
+				vec yl, yu ; 
+				d->get(i, yl, yu) ;
+
+				const double qi = r->q(d->get(i), j-np) ;
+				a0_norm += qi*qi * (yl[0]*yl[0]) ;
+				qp.set_a(j, i, ET(-yl[0] * qi)) ;
+				CI(j, i) = -yl[0] * qi ;
+				
+				a1_norm += qi*qi * (yu[0]*yu[0]) ;
+				qp.set_a(j, i+d->size(),  ET(yu[0] * qi)) ;
+				CI(j, i+d->size()) = yu[0] * qi ;
+			}
+		}
+	
+		// Set the c vector, will later be updated using the
+		// delta parameter.
+		ci(i) = sqrt(a0_norm) ;
+		ci(i+d->size()) = sqrt(a1_norm) ;
+	}
+#ifdef DEBUG
+	for(int j=0; j<d->size()*2; ++j)
+	{
+		for(int i=0; i<np+nq; ++i)
+			std::cout << CI(i,j) << "\t";
+		std::cout << std::endl; 
+	}
+	std::cout << std::endl ;
+#endif
+	// Update the ci column with the delta parameter
+	// (See Celis et al. 2007 p.12)
+	Eigen::JacobiSVD<Eigen::MatrixXd> svd(CI);
+	const double sigma_m = svd.singularValues()(std::min(2*d->size(), np+nq)-1) ;
+	const double sigma_M = svd.singularValues()(0) ;
+
+#ifdef DEBUG
+	std::cout << "<<DEBUG>> SVD = [ " ;
+	for(int i=0; i<std::min(2*d->size(), np+nq); ++i)
+	{
+		std::cout << svd.singularValues()(i) << ", " ;
+	}
+	std::cout << " ]" << std::endl ;
+#endif
+	
+	double delta = sigma_m / sigma_M ;
+	if(std::isnan(delta) || (std::abs(delta) == std::numeric_limits<double>::infinity()))
+	{
+#ifdef DEBUG
+		std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
+#endif
+		return false ;
+	}
+	else if(delta == 0.0)
+	{
+		delta = 1.0 ;
+	}
+
+#ifdef DEBUG
+	std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
+#endif
+	for(int i=0; i<2*d->size(); ++i)	
+	{		
+		qp.set_b(i, ET(delta * ci(i))) ;
+	}
+
+#ifdef DEBUG
+	// Export some informations on the problem to solve
+	std::cout << "<<DEBUG>> " << qp.get_n() << " variables" << std::endl ;
+	std::cout << "<<DEBUG>> " << qp.get_m() << " constraints" << std::endl ;
+#endif
+
+	if(qp.is_linear() && !qp.is_nonnegative())
+	{
+		std::cerr << "<<ERROR>> the problem should not be linear or negative!" << std::endl ;
+		return false ;
+	}
+
+#ifdef EXTREM_DEBUG
+	// Iterate over the rows
+	std::cout << std::endl ;
+	std::cout << "A = [" ;
+	for(int j=0; j<qp.get_m(); ++j)
+	{
+		if(j == 0) std::cout << "   " ;
+
+		// Iterate over the columns
+		for(int i=0; i<qp.get_n(); ++i)
+		{
+			if(i > 0) std::cout << ",\t" ;
+			std::cout << *(*(qp.get_a()+i) +j) ;
+		}
+		std::cout << ";" ;
+		if(j < qp.get_n()-1) std::cout << std::endl ;
+	}
+	std::cout << "]" << std::endl << std::endl ;
+	
+	std::cout << std::endl ;
+	std::cout << "D = [" ;
+	for(int j=0; j<np+nq; ++j)
+	{
+		if(j == 0) std::cout << "   " ;
+
+		// Iterate over the columns
+		for(int i=0; i<np+nq; ++i)
+		{
+			if(i > 0) std::cout << ",\t" ;
+			std::cout << *(*(qp.get_d()+i) +j) ;
+		}
+		std::cout << ";" ;
+	}
+	std::cout << "]" << std::endl << std::endl ;
+#endif
+
+	// solve the program, using ET as the exact type
+	Options  o ;
+	o.set_auto_validation(true) ;
+	Solution s = CGAL::solve_quadratic_program(qp, ET(), o) ;
+
+#ifdef DEBUG
+	if(s.is_infeasible())
+	{
+		std::cout << "<<DEBUG>> the current program is infeasible" << std::endl ;
+	}
+#endif
+
+	bool solves_qp = !s.is_infeasible() && s.solves_quadratic_program(qp)  ;
+	for(int i=0; i<np+nq; ++i)
+	{
+		const double v = (double)CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
+		solves_qp = solves_qp && !std::isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
+	}
+
+	if(solves_qp)
+	{
+		// Recopy the vector d
+		std::vector<double> p, q;
+		for(int i=0; i<np+nq; ++i)
+		{
+			const double v = CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
+
+			if(i < np)
+			{
+				p.push_back(v) ;
+			}
+			else
+			{
+				q.push_back(v) ;
+			}
+		}
+
+		if(r != NULL)
+		{
+			delete r ;
+		}
+		r = new rational_function(p, q);
+			
+		std::cout << "<<INFO>> got solution " << *r << std::endl ;
+		return true;
+	}
+	else
+	{
+		return false; 
+	}