Merging

sources/core/data.h

@@ -64,20 +64,7 @@ class data : public parametrized
 		// Get data size, e.g. the number of samples to fit
 		virtual int size() const = 0 ;
 
-
-		// Get min and max input space values
-		virtual vec min() const = 0 ;
-		virtual vec max() const = 0 ;
-
-
-		// Get the size of the input X domain and output Y domain
-		virtual int dimX() const { return _nX ; }
-		virtual int dimY() const { return _nY ; }
-
 	protected: // data
-
-		// Dimensions of the data
-		int _nX, _nY ;
 } ;
 
 /*! \brief Change the parametrization of data to fit the parametrization of the

sources/core/function.h

@@ -37,32 +37,6 @@ class function : public parametrized
 		//! example. The default behaviour is to load a function file.
 		virtual void bootstrap(const data*, const arguments& args);
 
-		//! Provide the dimension of the input space of the function
-		virtual int dimX() const { return _nX ; }
-		//! Provide the dimension of the output space of the function
-		virtual int dimY() const { return _nY ; }
-
-		//! Set the dimension of the input space of the function
-		virtual void setDimX(int nX) { _nX = nX ; }
-		//! Set the dimension of the output space of the function
-		virtual void setDimY(int nY) { _nY = nY ; }
-
-
-		/* DEFINITION DOMAIN OF THE FUNCTION */
-
-		//! \brief Set the minimum value the input can take
-		virtual void setMin(const vec& min) ;
-
-		//! \brief Set the maximum value the input can take
-		virtual void setMax(const vec& max) ;
-
-		//! \brief Get the minimum value the input can take
-		virtual vec getMin() const ;
-
-		//! \brief Get the maximum value the input can take
-		virtual vec getMax() const ;
-
-
 		/* IMPORT/EXPORT FUNCTIONS */
 
 		//! Load function specific files
@@ -109,10 +83,6 @@ class function : public parametrized
 
 
 	protected: // data
-
-		// Dimension of the function & domain of definition.
-		int _nX, _nY ;
-		vec _min, _max ;
 };
 
 /*! \brief Non-linear function interface

sources/core/params.h

@@ -335,8 +335,41 @@ class parametrized
 			}
 		}
 
+
+		/* DIMENSION OF THE INPUT AND OUTPUT DOMAIN */
+
+		//! Provide the dimension of the input space of the function
+		virtual int dimX() const { return _nX ; }
+		//! Provide the dimension of the output space of the function
+		virtual int dimY() const { return _nY ; }
+
+		//! Set the dimension of the input space of the function
+		virtual void setDimX(int nX) { _nX = nX ; }
+		//! Set the dimension of the output space of the function
+		virtual void setDimY(int nY) { _nY = nY ; }
+
+
+		/* DEFINITION DOMAIN OF THE FUNCTION */
+
+		//! \brief Set the minimum value the input can take
+		virtual void setMin(const vec& min) ;
+
+		//! \brief Set the maximum value the input can take
+		virtual void setMax(const vec& max) ;
+
+		//! \brief Get the minimum value the input can take
+		virtual vec getMin() const ;
+
+		//! \brief Get the maximum value the input can take
+		virtual vec getMax() const ;
+
+
 	protected:
 		// Input and output parametrization
 		params::input  _in_param ;
 		params::output _out_param ;
+
+		// Dimension of the function & domain of definition.
+		int _nX, _nY ;
+		vec _min, _max ;
 };

sources/core/rational_function.h

@@ -21,7 +21,7 @@ class rational_function_1d : public function
 	public: // methods
 
 		rational_function_1d() ;
-		rational_function_1d(int np, int nq, bool separable = true) ;
+		rational_function_1d(int np, int nq, bool separable = false) ;
         rational_function_1d(const vec& a, const vec& b) ;
 		virtual ~rational_function_1d() {}
 

sources/plugins/data_interpolant/data.cpp

+#include "data.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+
+data_interpolant::data_interpolant()
+{
+	_kdtree = new flann::Index< flann::L2<double> >(flann::KDTreeIndexParams(4));
+
+	_knn = 10;
+}
+
+data_interpolant::~data_interpolant()
+{
+	delete _data;
+	delete _kdtree;
+}
+
+// Load data from a file
+void data_interpolant::load(const std::string& filename) 
+{
+	// Load the data
+	_data->load(filename);
+
+	// Copy the informations
+	setDimX(_data->dimX());
+	setDimY(_data->dimY());
+	setMin(_data->getMin());
+	setMax(_data->getMax());
+
+	// Update the KDtreee by inserting all points
+	for(int i=0; i<_data->size(); ++i)
+	{
+		vec x = _data->get(i);
+		flann::Matrix<double> pts(&x[0], dimX(), 1);
+
+		_kdtree->addPoints(pts);
+	}
+	_kdtree->buildIndex();
+}
+void data_interpolant::load(const std::string& filename, const arguments&)
+{
+	load(filename);
+}
+
+void data_interpolant::save(const std::string& filename) const 
+{
+}
+
+// Acces to data
+vec data_interpolant::get(int id) const 
+{
+	vec res(dimX() + dimY()) ;
+	return res ;
+}
+vec data_interpolant::operator[](int i) const 
+{
+	return get(i) ;
+}
+
+//! \todo Test this function
+void data_interpolant::set(vec x)
+{
+	assert(x.size() == dimX());
+}
+
+vec data_interpolant::value(vec, vec) const
+{
+	vec res(dimY());
+	std::cerr << "<<ERROR>> Deprecated function: " << __func__ << std::endl;
+	return res;
+}
+vec data_interpolant::value(vec x) const
+{
+	vec res = vec::Zero(dimY());
+
+	// Query point
+	flann::Matrix<double> pts(&x[0], dimX(), 1);
+	std::vector< std::vector<int> >    indices;
+	std::vector< std::vector<double> > dists;
+
+	_kdtree->knnSearch(pts, indices, dists, _knn, flann::SearchParams());
+
+	// Interpolate the value using the indices
+	double cum_dist = 0.0;
+	for(int i=0; i<indices[0].size(); ++i)
+	{
+		int indice = indices[0][i];
+		vec y = _data->get(indice);
+
+		for(int j=0; j<dimY(); ++j)
+		{
+			res[j] += dists[0][i] * y[dimX() + j];
+		}
+		cum_dist += dists[0][i];
+	}
+	for(int j=0; j<dimY(); ++j)
+	{
+		res[j] /= cum_dist;
+	}
+
+   return res;
+}
+
+// Get data size, e.g. the number of samples to fit
+int data_interpolant::size() const 
+{
+	return _data->size();
+}
+
+ALTA_DLL_EXPORT data* provide_data()
+{
+    return new data_interpolant();
+}
+
+

sources/plugins/data_interpolant/data.h

+#pragma once
+
+#include <core/data.h>
+#include <core/common.h>
+#include <core/args.h>
+
+#include <flann/flann.hpp>
+
+/*! \brief Load a data file, but provide access to an interpolated version of
+ *  the data points.
+ *
+ *  <h2>Requirements:</h2>
+ *  The FLANN library.
+ *  On linux plateforms it can be obtained using the package manager: 
+ *  \verbatim 
+ *		sudo apt-get install libflann-dev
+ *	\endverbatim
+ */
+class data_interpolant : public data
+{
+	public: // methods
+
+		data_interpolant();
+		~data_interpolant();
+
+		// Load data from a file
+		virtual void load(const std::string& filename) ;
+		virtual void load(const std::string& filename, const arguments& args) ;
+
+		virtual void save(const std::string& filename) const ;
+
+		// Acces to data
+		virtual vec get(int i) const ;
+		virtual vec operator[](int i) const ;
+
+		virtual vec value(vec in, vec out) const;
+		virtual vec value(vec x) const;
+
+		// Set data
+		virtual void set(vec x);
+
+		// Get data size, e.g. the number of samples to fit
+		virtual int size() const ;
+
+	private: // data
+		
+		// The data object used to load sparse points sets
+		data* _data;
+
+		// Interpolation 
+		flann::Index< flann::L2<double> >* _kdtree;
+		flann::Matrix<double>* _points;
+		int _knn;
+} ;

sources/plugins/data_interpolant/data_interpolant.pro

+TEMPLATE        = lib
+CONFIG         *=  plugin \
+						 eigen
+
+DESTDIR         = ../../build
+
+INCLUDEPATH    += ../..
+HEADERS         = data.h
+SOURCES         = data.cpp
+
+
+LIBS           += -L../../build \
+						-lcore
+
+

sources/plugins/nonlinear_fitter_ceres/fitter.cpp

@@ -176,8 +176,8 @@ bool nonlinear_fitter_ceres::fit_data(const data* d, function* fit, const argume
     // to the dimension of my fitting problem
     fit->setDimX(d->dimX()) ;
     fit->setDimY(d->dimY()) ;
-    fit->setMin(d->min()) ;
-    fit->setMax(d->max()) ;
+    fit->setMin(d->getMin()) ;
+    fit->setMax(d->getMax()) ;
 
 	 // Convert the function and bootstrap it with the data
     if(dynamic_cast<nonlinear_function*>(fit) == NULL)

sources/plugins/nonlinear_fitter_eigen/fitter.cpp

@@ -273,8 +273,8 @@ bool nonlinear_fitter_eigen::fit_data(const data* d, function* fit, const argume
     // to the dimension of my fitting problem
     fit->setDimX(d->dimX()) ;
     fit->setDimY(d->dimY()) ;
-    fit->setMin(d->min()) ;
-    fit->setMax(d->max()) ;
+    fit->setMin(d->getMin()) ;
+    fit->setMax(d->getMax()) ;
 
 	 // Convert the function and bootstrap it with the data
     if(dynamic_cast<nonlinear_function*>(fit) == NULL)

sources/plugins/nonlinear_fitter_ipopt/fitter.cpp

@@ -218,8 +218,8 @@ bool nonlinear_fitter_ipopt::fit_data(const data* d, function* fit, const argume
 	// to the dimension of my fitting problem
 	fit->setDimX(d->dimX()) ;
 	fit->setDimY(d->dimY()) ;
-	fit->setMin(d->min()) ;
-	fit->setMax(d->max()) ;
+	fit->setMin(d->getMin()) ;
+	fit->setMax(d->getMax()) ;
 
 	// Convert the function and bootstrap it with the data
 	if(dynamic_cast<nonlinear_function*>(fit) == NULL)

sources/plugins/nonlinear_fitter_nlopt/fitter.cpp

@@ -122,8 +122,8 @@ bool nonlinear_fitter_nlopt::fit_data(const data* d, function* fit, const argume
 	// to the dimension of my fitting problem
 	fit->setDimX(d->dimX()) ;
 	fit->setDimY(d->dimY()) ;
-	fit->setMin(d->min()) ;
-	fit->setMax(d->max()) ;
+	fit->setMin(d->getMin()) ;
+	fit->setMax(d->getMax()) ;
 
 	// Convert the function and bootstrap it with the data
 	if(dynamic_cast<nonlinear_function*>(fit) == NULL)

sources/plugins/plugins.pro

@@ -28,5 +28,6 @@ SUBDIRS  = \
 				nonlinear_function_isotropic_lafortune	\
 				data_merl										\
 				data_brdf_slice								\
+				data_interpolant								\
 #				data_astm
 

sources/softs/data2moments/main.cpp

@@ -15,189 +15,277 @@
 
 int main(int argc, char** argv)
 {
-    arguments args(argc, argv) ;
-
-    if(args.is_defined("help")) {
-        std::cout << "<<HELP>> data2moments --input data.file --output gnuplot.file --data loader.so" << std::endl ;
-        std::cout << " - input, output and data are mandatory parameters" << std::endl ;
-        return 0 ;
-    }
-
-    if(! args.is_defined("input")) {
-        std::cerr << "<<ERROR>> the input filename is not defined" << std::endl ;
-        return 1 ;
-    }
-    if(! args.is_defined("output")) {
-        std::cerr << "<<ERROR>> the output filename is not defined" << std::endl ;
-        return 1 ;
-    }
-    if(! args.is_defined("data")) {
-        std::cerr << "<<ERROR>> the data loader is not defined" << std::endl ;
-        return 1 ;
-    }
-
-    // Import data
-    data* d = NULL ;
-    d = plugins_manager::get_data(args["data"]) ;
-    d->load(args["input"]);
-
-    // Create output file
-    std::ofstream file(args["output"].c_str(), std::ios_base::trunc);
-
-    if(d != NULL)
-    {
-        // Data parametrization
-        params::input data_param = d->parametrization();
-        const int nX = d->dimX();
-        const int nY = d->dimY();
-
-        // Options
-        bool with_cosine = args.is_defined("cosine");
-        bool use_angular = args.is_defined("angular-moments");
-
-        // Normalisation factor for the integration
-        const int nb_theta_int = 90;
-        const int nb_phi_int   = (use_angular) ? 2 : 360;
-        const double normalization = ((use_angular) ? M_PI : M_PI*M_PI) / (double)(nb_phi_int*nb_theta_int);
-
-        double in_angle[4] = {0.0, 0.0, 0.0, 0.0} ;
-
-
-        // Raw moments
-        vec m_0(nY);
-        vec m_x(nY);
-        vec m_y(nY);
-        vec m_xx(nY);
-        vec m_xy(nY);
-        vec m_yy(nY);
-
-        // Cumulants
-        vec k_x(nY);
-        vec k_y(nY);
-        vec k_xx(nY);
-        vec k_xy(nY);
-        vec k_yy(nY);
-
-        // The X and Y directions to compute the moments. This is an argument of the command
-        // line. The different directions can be: using stereographic coordinates, using the
-        // theta of the classical parametrization (the second coordinate is then 0).
-        vec xy(2);
-
-        for(int theta_in=0; theta_in<90; theta_in++)
-        {
-            in_angle[0] = theta_in * 0.5*M_PI / 90.0;
-
-            m_0  = vec::Zero(nY);
-            m_x  = vec::Zero(nY);
-            m_y  = vec::Zero(nY);
-            m_xx = vec::Zero(nY);
-            m_xy = vec::Zero(nY);
-            m_yy = vec::Zero(nY);
-
-            // Theta angle in [0 .. PI / 2]
-            for(int theta_out=0; theta_out<nb_theta_int; theta_out++)
-            {
-                in_angle[2] = theta_out * 0.5*M_PI / (double)nb_theta_int;
-
-                // Azimutal angle in [-PI .. PI]
-                for(int phi_out=0; phi_out<nb_phi_int; phi_out++)
-                {
-                    in_angle[3] = phi_out * 2.0*M_PI / (double)nb_phi_int - M_PI;
-
-                    vec in(nX), stereographics(4);
-                    params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, data_param, &in[0]);
-
-                    if(use_angular)
-                    {
-                        xy[0] = (std::abs(in_angle[3]) < 0.5*M_PI) ? in_angle[2] : -in_angle[2];
-                        xy[1] = 0.0;
-                    }
-                    else
-                    {
-                        params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, params::STEREOGRAPHIC, &stereographics[0]);
-                        xy[0] = stereographics[2];
-                        xy[1] = stereographics[3];
-                    }
-
-                    // Copy the input vector
-                    vec x = d->value(in);
-
-
-                    for(int i=0; i<nY; ++i)
-                    {
-                        double val = x[i] * normalization;
-                        if(!use_angular) {
-                            val *= sin(in_angle[2]);
-                        }
-                        if(with_cosine)
-                        {
-                            val *= cos(in_angle[2]);
-                        }
-
-                        m_0[i]  += val ;
-                        m_x[i]  += val * xy[0];
-                        m_y[i]  += val * xy[1];
-                        m_xx[i] += val * xy[0] * xy[0];
-                        m_xy[i] += val * xy[0] * xy[1];
-                        m_yy[i] += val * xy[1] * xy[1];
-                    }
-                }
-            }
-
-
-            for(int i=0; i<nY; ++i)
-            {
-                m_x[i]  /= m_0[i];
-                m_y[i]  /= m_0[i];
-                m_xx[i] /= m_0[i];
-                m_xy[i] /= m_0[i];
-                m_yy[i] /= m_0[i];
-            }
-
-            // compute cumulants
-            k_x  = m_x;
-            k_y  = m_y;
-            k_xx = m_xx - product(m_x,m_x);
-            k_xy = m_xy - product(m_x,m_y);
-            k_yy = m_yy - product(m_y,m_y);
-
-            // Output the value into the file
-            file << in_angle[0] << "\t";
-
-            for(int i=0; i<nY; ++i)
-                file << m_0[i] << "\t";
-
-            for(int i=0; i<nY; ++i)
-                file << k_x[i] << "\t";
-
-            if(!use_angular) {
-                for(int i=0; i<nY; ++i) {
-                    file << k_y[i] << "\t";
-                }
-            }
-
-            for(int i=0; i<nY; ++i)
-                file << k_xx[i] << "\t";
-
-            if(!use_angular) {
-                for(int i=0; i<nY; ++i) {
-                    file << k_xy[i] << "\t";
-                }
-
-                for(int i=0; i<nY; ++i) {
-                    file << k_yy[i] << "\t";
-                }
-            }
-            file << std::endl;
-
-        }
-
-        file.close();
-    }
-    else
-    {
-        std::cerr << "<<ERROR>> load file \"" << args["input"] << "\"" << std::endl ;
-    }
-
-    return 0 ;
+	arguments args(argc, argv) ;
+
+	if(args.is_defined("help")) {
+		std::cout << "<<HELP>> data2moments --input data.file --output gnuplot.file --data loader.so" << std::endl ;
+		std::cout << " - input, output and data are mandatory parameters" << std::endl ;
+		return 0 ;
+	}
+
+	if(! args.is_defined("input")) {
+		std::cerr << "<<ERROR>> the input filename is not defined" << std::endl ;
+		return 1 ;
+	}
+	if(! args.is_defined("output")) {
+		std::cerr << "<<ERROR>> the output filename is not defined" << std::endl ;
+		return 1 ;
+	}
+	if(! args.is_defined("data")) {
+		std::cerr << "<<ERROR>> the data loader is not defined" << std::endl ;
+		return 1 ;
+	}
+
+	// Import data
+	data* d = NULL ;
+	d = plugins_manager::get_data(args["data"]) ;
+	d->load(args["input"]);
+
+	// Create output file
+	std::ofstream file(args["output"].c_str(), std::ios_base::trunc);
+
+	if(d != NULL)
+	{
+		// Data parametrization
+		params::input data_param = d->parametrization();
+		const int nX = d->dimX();
+		const int nY = d->dimY();
+
+		// Raw moments
+		vec m_0(nY);
+		vec m_x(nY);
+		vec m_y(nY);
+		vec m_xx(nY);
+		vec m_xy(nY);
+		vec m_yy(nY);
+
+		// Cumulants
+		vec k_x(nY);
+		vec k_y(nY);
+		vec k_xx(nY);
+		vec k_xy(nY);
+		vec k_yy(nY);
+
+#ifndef OLD
+		// Number of elements per dimension
+		const vec dim = args.get_vec("dim", nX, 100.0);
+		assert(dim.size() == nX);
+
+		// Compute the volume in which we integrate and then compute the
+		// dt to apply for each element.
+		double dt = 1.0;
+		int nb = 0;
+		for(int k=0; k<nX; ++k)
+		{
+			dt *= d->max()[k] - d->min()[k];
+			nb += int(dim[k]);
+		}
+		dt /= double(nb);
+		
+		// Set all values to zero
+		m_0  = vec::Zero(nY);
+		m_x  = vec::Zero(nY);
+		m_y  = vec::Zero(nY);
+		m_xx = vec::Zero(nY);
+		m_xy = vec::Zero(nY);
+		m_yy = vec::Zero(nY);