Updating the parallel rational fitter and the associated quadratic

program class. Now one can check if all the constraints of datas
are satisfied by a solution. I need to break down the code a little
more to implement progressive estimation. I will now work on a branch.

sources/core/common.h

@@ -152,6 +152,21 @@ class vec : public std::vector<double>
 
 			return res;
 		}
+        friend bool operator<(const vec& a, const vec& b)
+        {
+            bool lessthan = true ;
+            for(int i=0; i<a.size(); ++i)
+                lessthan &= a[i] < b[i];
+            return lessthan;
+        }
+        friend bool operator>(const vec& a, const vec& b)
+        {
+            bool greatthan = true ;
+            for(int i=0; i<a.size(); ++i)
+                greatthan &= a[i] > b[i];
+            return greatthan;
+        }
+
 /*
 		friend double norm(const vec& a)
 		{

sources/core/vertical_segment.cpp

@@ -7,6 +7,7 @@
 #include <limits>
 #include <algorithm>
 #include <cmath>
+#include <cassert>
 
 void vertical_segment::load(const std::string& filename) 
 {
@@ -143,18 +144,22 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
 	std::cout << "<<INFO>> " << _data.size() << " elements to fit" << std::endl ;
 }
 
-bool vertical_segment::get(int i, double& x, double& yl, double& yu) const
+void vertical_segment::get(int i, vec& x, vec& yl, vec& yu) const
 {
-	if(i >= (int)_data.size())
-	{
-		return false ;
-	}
-
-	x  = _data[i][0] ;
-	yl = _data[i][1] ;
-	yu = _data[i][2] ;
+#ifdef DEBUG
+    assert(i >= 0 && i < _data.size());
+#endif
+    x.resize(dimX()); yl.resize(dimY()) ; yu.resize(dimY()) ;
+    for(int j=0; j<dimX(); ++j)
+    {
+        x[j] = _data[i][j] ;
+    }
+    for(int j=0; j<dimY(); ++j)
+    {
+        yl[j] = _data[i][dimX() + dimY() + j] ;
+        yu[j] = _data[i][dimX() + 2*dimY() + j] ;
+    }
 
-	return true ;
 }
 		
 void vertical_segment::get(int i, vec& yl, vec& yu) const

sources/core/vertical_segment.h

@@ -19,7 +19,7 @@ class vertical_segment : public data
 		virtual void load(const std::string& filename, const arguments& args) ;
 
 		// Acces to data
-		virtual bool get(int i, double& x, double& yl, double& yu) const ;
+        virtual void get(int i, vec &x, vec &yl, vec &yu) const ;
 		virtual vec get(int i) const ;		
 		virtual void get(int i, vec& yl, vec& yu) const ;		
 		virtual vec operator[](int i) const ;

sources/plugins/rational_fitter_parallel/quadratic_program.h

@@ -4,6 +4,9 @@
 #include <Array.hh>
 #include <QuadProg++.hh>
 
+#include <core/rational_function.h>
+#include <core/vertical_segment.h>
+
 class quadratic_program
 {
 public:
@@ -115,6 +118,29 @@ public:
         return solves_qp;
     }
 
+    //! \brief Test all the constraints of the data
+    bool test_constraints(const rational_function* r, const vertical_segment* data)
+    {
+        int nb_failed = 0;
+        for(int n=0; n<data->size(); ++n)
+        {
+            vec x, yl, yu;
+            data->get(n, x, yl, yu);
+
+            vec y = r->value(x);
+            if(y < yl || y > yu)
+            {
+                nb_failed++;
+            }
+        }
+
+#ifdef DEBUG
+        std::cout << "<<TRACE>> " << nb_failed << " constraints where not satified." << std::endl;
+#endif
+
+        return nb_failed == 0;
+    }
+
 protected:
     int _np, _nq ;
     QuadProgPP::Matrix<double> CI;

sources/plugins/rational_fitter_parallel/rational_fitter.cpp

@@ -85,7 +85,7 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit)
 		double min_delta = std::numeric_limits<double>::max();
 		int nb_sol_found = 0;
 		int np, nq ;
-//		#pragma omp parallel for
+        #pragma omp parallel for
 		for(int j=1; j<i; ++j)
 		{
 			int temp_np = i - j;
@@ -170,140 +170,6 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
 		return false ;
 	}
 
-#ifdef OLD
-    // Matrices of the problem
-	QuadProgPP::Matrix<double> G (0.0, np+nq, np+nq) ;
-	QuadProgPP::Vector<double> g (0.0, 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, 0) ;
-	QuadProgPP::Vector<double> ce(0.0, 0) ;
-
-	Eigen::MatrixXd eCI(np+nq, 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[i][i] = 1.0 ; 
-	}
-	
-	// 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 ;
-
-		vec xi = d->get(i) ;
-
-		// 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(xi, j) ;
-				a0_norm += pi*pi ;
-				a1_norm += pi*pi ;
-				CI[j][i] =  pi ;
-				CI[j][i+d->size()] = -pi ;
-
-				// Updating Eigen matrix
-				eCI(j,i) = pi ;
-				eCI(j,i+d->size()) = -pi ;
-			}
-			// Filling the q part
-			else
-			{
-				vec yl, yu ; 
-				d->get(i, yl, yu) ;
-
-				const double qi = r->q(xi, j-np) ;
-				a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
-				CI[j][i] = -yl[ny] * qi ;
-				
-				a1_norm += qi*qi * (yu[ny]*yu[ny]) ;
-				CI[j][i+d->size()] = yu[ny] * qi ;
-				
-				// Updating Eigen matrix
-				eCI(j,i) = -yl[ny] * qi ;
-				eCI(j,i+d->size()) = yu[ny] * 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) ;
-	}
-
-
-
-#ifndef DEBUG
-    std::cout << CI << std::endl ;
-#endif
-	// Update the ci column with the delta parameter
-	// (See Celis et al. 2007 p.12)
-	Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(eCI);
-	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
-	
-	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)	
-	{		
-		ci[i] = ci[i] / delta ; 
-#ifdef DEBUG
-		std::cout << ci[i] << "\t" ;
-#endif
-	}
-#ifdef DEBUG
-	std::cout << std::endl << std::endl ;
-
-	std::cout << eCI << std::endl << std::endl ;
-#endif
-	// Compute the solution
-	QuadProgPP::Vector<double> x;
-	double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, x);
-
-
-	bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
-	for(int i=0; i<np+nq; ++i)
-	{
-		const double v = x[i];
-		solves_qp = solves_qp && !std::isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
-	}
-#else
     quadratic_program qp(np, nq);
     for(int i=0; i<d->size(); ++i)
     {
@@ -339,7 +205,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
 
     QuadProgPP::Vector<double> x(np+nq);
     bool solves_qp = qp.solve_program(x, delta);
-#endif
+
 	if(solves_qp)
 	{
 		// Recopy the vector d
@@ -358,7 +224,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
 			}
 		}
 
-#ifndef DEBUG
+#ifdef DEBUG
 		std::cout << "<<INFO>> got solution " << *r << std::endl ;
 #endif
 		return norm > 0.0;
@@ -370,4 +236,35 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
 	}
 }
 
+void rational_fitter_parallel::get_constraint(int i, int np, int nq, int ny, const vertical_segment* data, const rational_function* func, vec& cu, vec& cl)
+{
+    const vec xi = data->get(i) ;
+    cu.resize(np+nq);
+    cl.resize(np+nq);
+
+    // Create two vector of constraints
+    for(int j=0; j<np+nq; ++j)
+    {
+        // Filling the p part
+        if(j<np)
+        {
+            const double pi = func->p(xi, j) ;
+            cu[j] =  pi ;
+            cl[j] = -pi ;
+
+        }
+        // Filling the q part
+        else
+        {
+            vec yl, yu ;
+            data->get(i, yl, yu) ;
+            const double qi = func->q(xi, j-np) ;
+
+            cl[j] = -yl[ny] * qi ;
+            cu[j] = yu[ny] * qi ;
+        }
+    }
+
+}
+
 Q_EXPORT_PLUGIN2(rational_fitter_parallel, rational_fitter_parallel)

sources/plugins/rational_fitter_parallel/rational_fitter.h

@@ -50,6 +50,10 @@ class rational_fitter_parallel : public QObject, public fitter
 		virtual bool fit_data(const vertical_segment* d, int np, int nq, rational_function* fit, vec& p, vec& q, double& delta) ;
 		virtual bool fit_data(const vertical_segment* dat, int np, int nq, int ny, rational_function* fit, vec& p, vec& q, double& delta) ;
 
+        //! \brief Create a constraint vector given its index i in the data
+        // object and the rational function object to fit.
+        virtual void get_constraint(int i, int np, int nq, int ny, const vertical_segment* data, const rational_function* func, vec& cu, vec& cl);
+
 		// min and Max usable np and nq values for the fitting
 		int _max_np, _max_nq ;
 		int _min_np, _min_nq ;