Add parsec directory to solve multiple problem in parallel

sources/plugins/rational_fitter_parsec_multi/quadratic_program.h

+#pragma once
+
+#include <Eigen/SVD>
+#include <Array.hh>
+#include <QuadProg++.hh>
+
+#include <core/rational_function.h>
+#include <core/vertical_segment.h>
+
+class quadratic_program
+{
+ public:
+    //! \brief Constructor need to specify the number of coefficients
+ quadratic_program(int np, int nq, bool compute_delta = false) :
+    _np(np), _nq(nq), _compute_delta(compute_delta), CI(0.0, _np+_nq, 0)
+        { }
+
+    //! \brief Remove the already defined constraints
+    void clear_constraints()
+    {
+        CI.resize(_np+_nq, 0);
+    }
+
+    //! \brief Add a constraint by specifying the vector
+    void add_constraints(const vec& c)
+    {
+        const int m = CI.nrows();
+        const int n = CI.ncols();
+
+        if(n > 0)
+            {
+                // Construct temp buffer
+                double* temp = new double[n*m];
+                for(int u=0; u<n; ++u)
+                    {
+                        for(int v=0; v<m; ++v)
+                            {
+                                temp[u*m + v] = CI[v][u];
+                            }
+                    }
+
+                // Resize matrix CI
+                CI.resize(m, n+1);
+
+                // Recopy data
+                for(int u=0; u<n+1; ++u)
+                    {
+                        if(u==n)
+                            {
+                                for(int v=0; v<m; ++v)
+                                    CI[v][u] = c[v];
+                            }
+                        else
+                            {
+                                for(int v=0; v<m; ++v)
+                                    CI[v][u] = temp[u*m + v];
+                            }
+                    }
+                delete[] temp;
+            }
+        else
+            {
+                // Resize matrix CI
+                CI.resize(m, 1);
+
+                // Recopy data
+                for(int u=0; u<m; ++u)
+                    CI[n][u] = c[u];
+            }
+    }
+
+    //! \brief Provide the number of constraints
+    int nb_constraints() const
+    {
+        return CI.ncols();
+    }
+
+    //! Set the indices of the remaining data
+    void set_training_set(const std::list<unsigned int>& ts)
+    {
+        this->training_set = ts;
+    }
+
+    //! \brief Solves the quadratic program and update the p and
+    //! q vector if necessary.
+    inline bool solve_program(QuadProgPP::Vector<double>& x, double& delta, vec& p, vec& q)
+    {
+        bool solves_qp = solve_program(x, delta) ;
+
+        if(solves_qp)
+            {
+                double norm = 0.0;
+                for(int i=0; i<_np+_nq; ++i)
+                    {
+                        const double v = x[i];
+                        norm += v*v ;
+                        if(i < _np)
+                            {
+                                p[i] = v ;
+                            }
+                        else
+                            {
+                                q[i-_np] = v ;
+                            }
+                    }
+                return norm > 0.0;
+            }
+        else
+            {
+                return false ;
+            }
+    }
+
+    //! \brief Solves the quadratic program
+    inline bool solve_program(QuadProgPP::Vector<double>& v, double& delta)
+    {
+        const int m = CI.nrows();
+        const int n = CI.ncols();
+
+        QuadProgPP::Matrix<double> G (0.0, m, m) ;
+        QuadProgPP::Vector<double> g (0.0, m) ;
+        QuadProgPP::Vector<double> ci(0.0, n) ;
+        QuadProgPP::Matrix<double> CE(0.0, m, 0) ;
+        QuadProgPP::Vector<double> ce(0.0, 0) ;
+
+	if(_compute_delta)
+	    {
+		// Update the ci column with the delta parameter
+		// (See Celis et al. 2007 p.12)
+		Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(Eigen::MatrixXd::Map(&CI[0][0], m, n));
+		const double sigma_m = svd.singularValues()(std::min(m, n)-1) ;
+		const double sigma_M = svd.singularValues()(0) ;
+		delta = sigma_M / sigma_m ;
+	    }
+
+        // Select the size of the result vector to
+        // be equal to the dimension of p + q
+        for(int i=0; i<m; ++i)
+            {
+                G[i][i] = 1.0 ;
+            }
+
+        // Each constraint (fitting interval or point
+        // add another dimension to the constraint
+        // matrix
+        for(int i=0; i<n; ++i)
+            {
+                // Norm of the row vector
+                double norm = 0.0 ;
+
+                for(int j=0; j<m; ++j)
+                    {
+                        norm += CI[j][i]*CI[j][i] ; ;
+                    }
+
+                // Set the c vector, will later be updated using the
+                // delta parameter.
+                ci[i] = -delta * sqrt(norm) ;
+            }
+
+        // Compute the solution
+        const double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, v);
+
+        bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
+        return solves_qp;
+    }
+
+#define PACANOWSKI2012
+
+    //! \brief Test all the constraints of the data.
+    //! Add the sample that is farest away from the function.
+    bool test_constraints(int ny, const rational_function_1d* r, const vertical_segment* data)
+    {
+#ifdef PACANOWSKI2012
+        int nb_failed = 0;
+        double max_dev = 0.0; // Maximum absolute distance of the current
+        // solution to the data.
+        std::list<unsigned int>::iterator max_ind;
+        vec cu, cl;
+
+        std::list<unsigned int>::iterator it;
+        for(it = training_set.begin(); it != training_set.end(); it++)
+            {
+                vec x, yl, yu;
+                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];
+                if(fail_lower || fail_upper)
+                    {
+                        const double dev = std::abs(0.5*(yu[ny]+yl[ny]) - y[ny]);
+
+                        nb_failed++;
+
+                        if(max_dev < dev)
+                            {
+                                get_constraint(x, yl, yu, ny, r, cu, cl);
+                                max_dev = dev;
+                                max_ind = it;
+                            }
+                    }
+            }
+
+#ifdef DEBUG
+        std::cout << "<<TRACE>> " << nb_failed << " constraints where not satified." << std::endl;
+        std::cout << "<<TRACE>> an interval failed the test with distance = " << max_dev << std::endl;
+#endif
+
+        if(nb_failed > 0)
+            {
+                add_constraints(cu);
+                add_constraints(cl);
+                training_set.erase(max_ind);
+#ifdef DEBUG
+                std::cout << "<<DEBUG>> number of remaining training elements: " << training_set.size() << std::endl;
+#endif
+
+                return false;
+            }
+        else
+            {
+                return true;
+            }
+#else
+        int n = next_unmatching_constraint(0, ny, r, data);
+        if(n < data->size())
+            {
+                vec x, yl, yu;
+                data->get(n, x, yl, yu);
+
+                vec cu, cl;
+                get_constraint(x, yl, yu, ny, r, cu, cl);
+
+                add_constraints(cu);
+                add_constraints(cl);
+
+                return false;
+            }
+        else
+            {
+                return true;
+            }
+#endif
+    }
+
+    //! \brief Generate two constraint vectors from a vertical segment and a
+    //! ration function type.
+    inline void get_constraint(const vec& xi, const vec& yl, const vec& yu, int ny,
+                               const rational_function_1d* func,
+                               vec& cu, vec& cl);
+
+    //! \brief Give the next position in the data that is not satisfied.
+    //! This method works only for a single color channel ny !
+    static int next_unmatching_constraint(int i, int ny, const rational_function_1d* r,
+                                          const vertical_segment* data);
+
+ protected:
+    int _np, _nq;
+    bool _compute_delta;
+    QuadProgPP::Matrix<double> CI;
+
+    //! Contains the indices of the vertical segment unused during the
+    //! rational interpolation.
+    std::list<unsigned int> training_set;
+};
+
+
+inline void quadratic_program::get_constraint(const vec& xi, const vec& yl, const vec& yu,
+                                              int ny, const rational_function_1d* func,
+                                              vec& cu, vec& cl)
+{
+    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
+		{
+		    const double qi = func->q(xi, j-_np) ;
+
+		    cu[j] = -yu[ny] * qi ;
+		    cl[j] =  yl[ny] * qi ;
+		}
+	}
+}
+
+int quadratic_program::next_unmatching_constraint(int i, int ny, const rational_function_1d* r,
+                                                  const vertical_segment* data)
+{
+    for(int n=i; n<data->size(); ++n)
+        {
+            vec x, yl, yu;
+            data->get(n, x, yl, yu);
+
+            vec y = r->value(x);
+            if(y[ny] < yl[ny] || y[ny] > yu[ny])
+                {
+                    return n;
+                }
+        }
+    return data->size();
+}

sources/plugins/rational_fitter_parsec_multi/rational_fitter.cpp

+#include "rational_fitter.h"
+
+#include <core/plugins_manager.h>
+
+#include <Eigen/SVD>
+#include <Array.hh>
+#include <QuadProg++.hh>
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <limits>
+#include <algorithm>
+#include <cmath>
+#include <string>
+#include <list>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "quadratic_program.h"
+
+ALTA_DLL_EXPORT fitter* provide_fitter()
+{
+    return new rational_fitter_parallel();
+}
+
+rational_fitter_parallel::rational_fitter_parallel() : nb_starting_points(100)
+{
+}
+rational_fitter_parallel::~rational_fitter_parallel()
+{
+}
+
+bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const arguments &args)
+{
+    rational_function* r = dynamic_cast<rational_function*>(fit) ;
+    if(r == NULL)
+	{
+	    std::cerr << "<<ERROR>> not passing the correct function class to the fitter: must be a rational_function" << std::endl ;
+	    return false ;
+	}
+
+    const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
+    if(d == NULL)
+	{
+	    std::cerr << "<<WARNING>> automatic convertion of the data object to vertical_segment," << std::endl;
+	    std::cerr << "<<WARNING>> we advise you to perform convertion with a separate command." << std::endl;
+
+	    vertical_segment* vs = new vertical_segment();
+	    for(int i=0; i<dat->size(); ++i) 
+		{
+		    const vec x = dat->get(i);
+		    vec y(dat->dimX() + 3*dat->dimY());
+
+		    for(int k=0; k<x.size()   ; ++k) { y[k]                          = x[k]; }
+		    for(int k=0; k<dat->dimY(); ++k) { y[k + x.size()]               = (1.0 - args.get_float("dt", 0.1)) * x[k + dat->dimX() +   dat->dimY()]; }
+		    for(int k=0; k<dat->dimY(); ++k) { y[k + x.size() + dat->dimY()] = (1.0 + args.get_float("dt", 0.1)) * x[k + dat->dimX() + 2*dat->dimY()]; }
+
+		    vs->set(y);
+		}
+
+	    d = vs;
+	}
+
+    // 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->setMin(d->min()) ;
+    r->setMax(d->max()) ;
+
+    const int _min_np = args.get_int("min-np", 10);
+    const int _max_np = args.get_int("np", _min_np);
+    std::cout << "<<INFO>> N in  [" << _min_np << ", " << _max_np << "]"  << std::endl ;
+
+    const int nb_starting_points = args.get_int("nb-starting-points", 100);
+    std::cout << "<<INFO>> number of data point used in start: " << nb_starting_points << std::endl;
+
+    const int step = args.get_int("np-step", 1);
+    const bool use_delta = args.is_defined("use_delta");
+
+    for(int i=_min_np; i<=_max_np; i+=step)
+        {
+            std::cout << "<<INFO>> fit using np+nq = " << i << std::endl ;
+            std::cout.flush() ;
+            timer time ;
+            time.start() ;
+
+            int nb_cores = args.get_int("nb-cores", omp_get_num_procs());
+
+#ifdef DEBUG
+            std::cout << "<<DEBUG>> will use " << nb_cores << " threads to compute the quadratic programs" << std::endl ;
+#endif
+
+            omp_set_num_threads(nb_cores) ;
+
+	    double min_delta   = std::numeric_limits<double>::max();
+	    double min_l2_dist = std::numeric_limits<double>::max();
+	    double mean_delta = 0.0;
+	    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)
+	    for(int j=1; j<i; ++j)
+		{
+		    // Compute the number of coefficients in the numerator and in the denominator
+		    // from the current number of coefficients i and the current index in the
+		    // loop j.
+		    int temp_np = i - j;
+		    int temp_nq = j;
+
+                    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.
+		    rational_function* rk = NULL;
+#pragma omp critical (args)
+		    {
+			rk = dynamic_cast<rational_function*>(plugins_manager::get_function(args));
+		    }
+		    if(rk == NULL)
+			{
+			    std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
+			    throw;
+			}
+		    rk->setParametrization(r->input_parametrization());
+		    rk->setParametrization(r->output_parametrization());
+		    rk->setDimX(r->dimX()) ;
+		    rk->setDimY(r->dimY()) ;
+		    rk->setMin(r->min()) ;
+		    rk->setMax(r->max()) ;
+
+                    // Set the rational function size
+                    rk->setSize(temp_np, temp_nq);
+
+		    double delta = 1.0;
+		    double linf_dist, l2_dist;
+		    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)
+			    {
+				++nb_sol_found ;
+				mean_delta += delta ;
+
+                                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;
+                                std::cout << std::endl;
+
+				// Get the solution with the minimum delta or the minimum L2 distance, 
+				// and update the main rational function r.
+				if((use_delta && delta < min_delta) || (!use_delta && l2_dist < min_l2_dist))
+				    {
+					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));
+					    }
+				    }
+			    }
+
+			    if(rk != NULL)
+				delete rk; // memory clean
+
+			    time.stop();
+			    std::cout << "<<INFO>> got a fit using N = " << i << std::endl ;
+			    std::cout << "<<INFO>> it took " << time << std::endl ;
+			    std::cout << "<<INFO>> I got " << nb_sol_found << " solutions to the QP" << std::endl ;
+			    return true ;
+			}
+		}
+
+	    return false ;
+	}
+
+void rational_fitter_parallel::set_parameters(const arguments&)
+{
+}
+
+
+bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, 
+                                        rational_function* r, const arguments &args,
+                                        double& delta, double& linf_dist, double& l2_dist)
+{
+	// Fit the different output dimension independantly
+	for(int j=0; j<d->dimY(); ++j)
+	{
+		vec p(np), q(nq);
+		rational_function_1d* rf = r->get(j);
+		rf->resize(np, nq);
+
+		if(!fit_data(d, np, nq, j, rf, args, p, q, delta))
+		{
+			return false ;
+		}
+            rf->update(p, q);
+        }
+
+    linf_dist = r->Linf_distance(d);
+    l2_dist   = r->L2_distance(d);
+
+    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 returns a rational BRDF function and a boolean
+bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, int ny,
+                                        rational_function_1d* r, const arguments& args,
+                                        vec& p, vec& q, double& delta)
+{
+    const int m = d->size(); // 2*m = number of constraints
+    const int n = np+nq;     // n = np+nq
+
+    quadratic_program qp(np, nq, args.is_defined("use_delta"));
+
+    // Starting with only a nb_starting_points vertical segments
+    std::list<unsigned int> training_set;
+    const int di = std::max((m-1) / (nb_starting_points-1), 1);
+    for(int i=0; i<m; ++i)
+        {
+            if(i % di == 0)
+                {
+                    // Create two vector of constraints
+                    vec c1(n), c2(n);
+                    get_constraint(i, np, nq, ny, d, r, c1, c2);
+
+                    qp.add_constraints(c1);
+                    qp.add_constraints(c2);
+
+                }
+            else
+                {
+                    training_set.push_back(i);
+                }
+        }
+    qp.set_training_set(training_set);
+
+    do
+	{
+#ifdef _OPENMP
+#ifdef DEBUG
+            std::cout << "<<DEBUG>> thread " << omp_get_thread_num() << ", number of intervals tested = " << qp.nb_constraints()/2 << std::endl ;
+#endif
+#endif
+		QuadProgPP::Vector<double> x(n);
+		bool solves_qp = qp.solve_program(x, delta, p, q);
+		r->update(p, q);
+
+            if(solves_qp)
+                {
+                    if(qp.test_constraints(ny, r, d))
+                        {
+#ifdef DEBUG
+                            std::cout << "<<INFO>> got solution " << *r << std::endl ;
+#endif
+                            return true;
+                        }
+                }
+            else
+                {
+#ifdef DEBUG
+                    std::cout << "<<DEB