Cleaning quadprog and matlab code to obtain the same thing

sources/core/function.cpp

@@ -48,10 +48,11 @@ vec function::getMax() const
 
 void function::save(const std::string& filename, const arguments& args) const
 {
-    bool is_cpp      = args["export"] == "C++";
-    bool is_explorer = args["export"] == "explorer";
-    bool is_shader   = args["export"] == "shader" || is_explorer;
-    bool is_matlab   = args["export"] == "matlab";
+	bool is_alta     = !args.is_defined("export") || args["export"] == "alta";
+	bool is_cpp      = args["export"] == "C++";
+	bool is_explorer = args["export"] == "explorer";
+	bool is_shader   = args["export"] == "shader" || is_explorer;
+	bool is_matlab   = args["export"] == "matlab";
 
 	// Open the file
 	std::ofstream file(filename.c_str());
@@ -60,60 +61,66 @@ void function::save(const std::string& filename, const arguments& args) const
 		std::cerr << "<<ERROR>> unable to open output file for writing" << std::endl;
 	}
 
-    if(is_explorer)
-    {
-        file << "analytic" << std::endl;
-        file << std::endl;
-        file << "::begin shader" << std::endl;
-    }
+	// If the export is the alta format, use the maximum precision formatting
+	if(is_alta)
+	{
+		file.precision(10);
+	}
+
+	if(is_explorer)
+	{
+		file << "analytic" << std::endl;
+		file << std::endl;
+		file << "::begin shader" << std::endl;
+	}
 
 	// Save common header
-    save_header(file, args);
+	save_header(file, args);
 
 	// Save function definition
-    save_body(file, args);
+	save_body(file, args);
 
-    if(is_cpp)
-    {
-        file << "vec brdf(const vec& in, const vec& file)" << std::endl;
-        file << "{" << std::endl;
-        file << "\tvec res(" << dimY() << ");" << std::endl;
-        file << "\t";
-    }
-    else if(is_matlab)
-    {
-        file << "function res = brdf(in, file)" << std::endl;
-        file << "{" << std::endl;
-        file << "\tres = zeros(" << dimY() << ");" << std::endl;
-        file << "\t";
-    }
-    else if(is_shader)
-    {
-        file << "vec3 BRDF(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y)" << std::endl;
-        file << "{" << std::endl;
-        file << "\tvec3 res = ";
-    }
+	if(is_cpp)
+	{
+		file << "vec brdf(const vec& in, const vec& file)" << std::endl;
+		file << "{" << std::endl;
+		file << "\tvec res(" << dimY() << ");" << std::endl;
+		file << "\t";
+	}
+	else if(is_matlab)
+	{
+		file << "function res = brdf(in, file)" << std::endl;
+		file << "{" << std::endl;
+		file << "\tres = zeros(" << dimY() << ");" << std::endl;
+		file << "\t";
+	}
+	else if(is_shader)
+	{
+		file << "vec3 BRDF(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y)" << std::endl;
+		file << "{" << std::endl;
+		file << "\tvec3 res = ";
+	}
 
 	// Save fit data
-    save_call(file, args);
-    if(is_cpp || is_shader)
-    {
-        file << ";" << std::endl;
-        file << "\treturn res;" << std::endl;
-        file << "}" << std::endl;
-    }
-    else if(is_matlab)
-    {
-        file << ";" << std::endl;
-        file << "\treturn res;" << std::endl;
-        file << "endfunction" << std::endl;
-    }
-    file << std::endl;
+	save_call(file, args);
+	if(is_cpp || is_shader)
+	{
+		file << ";" << std::endl;
+		file << "\treturn res;" << std::endl;
+		file << "}" << std::endl;
+	}
+	else if(is_matlab)
+	{
+		file << ";" << std::endl;
+		file << "\treturn res;" << std::endl;
+		file << "endfunction" << std::endl;
+	}
+	file << std::endl;
 
-    if(is_explorer)
-    {
-        file << "::end shader" << std::endl;
-    }
+	if(is_explorer)
+	{
+		file << "::end shader" << std::endl;
+	}
 }
 		
 //! \brief save the header of the output function file. The header should
@@ -128,9 +135,10 @@ void function::save_header(std::ostream& out, const arguments& args) const
 		out << "#DIM " << _nX << " " << _nY << std::endl;
 		out << "#PARAM_IN  " << params::get_name(input_parametrization()) << std::endl;
 		//out << "#PARAM_OUT " << params::get_name(output_parametrization()) << std::endl;*
-        if(args.is_defined("export-append")) {
-            out << args["export-append"] << std::endl;
-        }
+		if(args.is_defined("export-append")) 
+		{
+			out << args["export-append"] << std::endl;
+		}
 		out << "#ALTA HEADER END" << std::endl;
 		out << std::endl;
 	}

sources/core/plugins_manager.cpp

@@ -208,6 +208,9 @@ function* plugins_manager::get_function(const std::string& filename)
 		throw ;
 	}
 
+	// Set the precision of the input
+	file.precision(10);
+
 	// Parameters of the function object
 	int nX, nY;
 	arguments args;

sources/core/rational_function.cpp

@@ -35,7 +35,13 @@ void rational_function_1d::update(const vec& in_a,
 	a.resize(in_a.size()) ;
 	b.resize(in_b.size()) ;
 
-    const double b0 = 1.0;//in_b[0];
+#define NORMALIZE
+
+#ifdef NORMALIZE
+    const double b0 = (std::abs(in_b[0]) > 1.0E-10) ? in_b[0] : 1.0;
+#else
+    const double b0 = 1.0;
+#endif
 
     for(int i=0; i<a.size(); ++i) { a[i] = in_a[i] / b0; }
     for(int i=0; i<b.size(); ++i) { b[i] = in_b[i] / b0; }

sources/plugins/rational_fitter_matlab/rational_fitter.cpp

@@ -124,11 +124,7 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 	// Size of the problem
 	int N = np+nq ;
 	int M = d->size() ;
-/*
-	// Get the maximum value in data to scale the input parameter space
-	// so that it reduces the values of the polynomial
-	vec dmax = d->max() ;
-*/
+	
 	// Matrices of the problem
 	Eigen::MatrixXd G (N, N) ;
 	Eigen::VectorXd g (N) ;
@@ -171,26 +167,27 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 			if(j<np)
 			{
 				const double pi = r->p(xi, j) ;
-				a0_norm += pi*pi ;
-				a1_norm += pi*pi ;
+				
 				// Updating Eigen matrix
-				CI(2*i+0, j) =  pi ;
-				CI(2*i+1, j) = -pi ;
+				CI(2*i+0, j) =  pi;
+				CI(2*i+1, j) = -pi;
 			}
 			// Filling the q part
 			else
 			{
-				vec yl, yu ; 
-				d->get(i, yl, yu) ;
+				vec yl, yu; 
+				d->get(i, yl, yu);
 
-				const double qi = r->q(xi, j-np) ;
-				a0_norm += qi*qi * (yu[ny]*yu[ny]) ;
-				a1_norm += qi*qi * (yl[ny]*yl[ny]) ;
+				const double qi = r->q(xi, j-np);
 				
 				// Updating Eigen matrix
-				CI(2*i+0, j) = -yu[ny] * qi ;
-				CI(2*i+1, j) =  yl[ny] * qi ;
+				CI(2*i+0, j) = -yu[ny] * qi;
+				CI(2*i+1, j) =  yl[ny] * qi;
 			}
+
+			// Update the norm of the row
+			a0_norm += CI(2*i+0, j)*CI(2*i+0, j);
+			a1_norm += CI(2*i+1, j)*CI(2*i+1, j);
 		}
 	
 		// Set the c vector, will later be updated using the

sources/plugins/rational_fitter_quadprog/rational_fitter.cpp

@@ -122,32 +122,25 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int n
 // 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_quadprog::fit_data(const vertical_segment* dat, int np, int nq, int ny, rational_function_1d* r) 
+bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int nq, int ny, rational_function_1d* r) 
 {
-	const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
-	if(r == NULL || d == NULL)
-	{
-		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
-		return false ;
-	}
-
-	// Get the maximum value in data to scale the input parameter space
-	// so that it reduces the values of the polynomial
-	vec dmax = d->max() ;
+	// Size of the problem
+	const int N = np+nq ;
+	const int M = d->size() ;
 
 	// 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::Matrix<double> G (0.0, N, N) ;
+	QuadProgPP::Vector<double> g (0.0, N) ;
+	QuadProgPP::Matrix<double> CI(0.0, N, 2*M) ;
+	QuadProgPP::Vector<double> ci(0.0, 2*M) ;
+	QuadProgPP::Matrix<double> CE(0.0, N, 0) ;
 	QuadProgPP::Vector<double> ce(0.0, 0) ;
 
-	Eigen::MatrixXd eCI(np+nq, 2*d->size()) ;
+	Eigen::MatrixXd eCI(2*M, N) ;
 
 	// Select the size of the result vector to
 	// be equal to the dimension of p + q
-	for(int i=0; i<np+nq; ++i)
+	for(int i=0; i<N; ++i)
 	{
 		G[i][i] = 1.0 ; 
 	}
@@ -155,7 +148,7 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
 	// Each constraint (fitting interval or point
 	// add another dimension to the constraint
 	// matrix
-	for(int i=0; i<d->size(); ++i)	
+	for(int i=0; i<M; ++i)	
 	{		
 		// Norm of the row vector
 		double a0_norm = 0.0 ;
@@ -163,30 +156,23 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
 
         vec xi = d->get(i) ;
 
-        bool is_boundary = false;
-        for(int k=0; k<d->dimX(); ++k)
-        {
-            is_boundary = is_boundary || (xi[k] <= d->min()[k]) || (xi[k] >= d->max()[k]);
-        }
-
 		// 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)
+		for(int j=0; j<N; ++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 ;
+				const double pi = r->p(xi, j);
+
+				CI[j][2*i+0] =  pi;
+				CI[j][2*i+1] = -pi;
 
 				// Updating Eigen matrix
-				eCI(j,i) = pi ;
-				eCI(j,i+d->size()) = -pi ;
+				eCI(2*i+0, j) = CI[j][2*i+0];
+				eCI(2*i+1, j) = CI[j][2*i+1];
 			}
 			// Filling the q part
 			else
@@ -194,67 +180,51 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
 				vec yl, yu ; 
 				d->get(i, yl, yu) ;
 
-                // Add a constraints for boundary conditions
-                if(is_boundary)
-                {
-                    vec mean = 0.5*(yl+yu);
-                    yl = mean + _boundary * (yl - mean);
-                    yu = mean + _boundary * (yu - mean);
-                }
-
-				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 ;
+				const double qi = r->q(xi, j-np);
+
+				CI[j][2*i+0] = -yu[ny] * qi;
+				CI[j][2*i+1] =  yl[ny] * qi;
 				
 				// Updating Eigen matrix
-				eCI(j,i) = -yl[ny] * qi ;
-				eCI(j,i+d->size()) = yu[ny] * qi ;
+				eCI(2*i+0, j) = CI[j][2*i+0];
+				eCI(2*i+1, j) = CI[j][2*i+1];
 			}
+
+			// Update the norm of the row
+			a0_norm += CI[j][2*i+0]*CI[j][2*i+0];
+			a1_norm += CI[j][2*i+1]*CI[j][2*i+1];
 		}
 	
 		// Set the c vector, will later be updated using the
 		// delta parameter.
-		ci[i] = -sqrt(a0_norm) ;
-		ci[i+d->size()] = -sqrt(a1_norm) ;
+		ci[2*i+0] = -sqrt(a0_norm) ;
+		ci[2*i+1] = -sqrt(a1_norm) ;
 	}
 #ifdef DEBUG
 	std::cout << "CI = [" ;
-	for(int j=0; j<d->size()*2; ++j)
+	for(int j=0; j<2*M; ++j)
 	{
-		for(int i=0; i<np+nq; ++i)
+		for(int i=0; i<N; ++i)
 		{
 			std::cout << CI[i][j] ;
-			if(i != np+nq-1) std::cout << ", ";
+			if(i != N-1) std::cout << ", ";
 		}
-		if(j != d->size()*2-1)
+		if(j != 2*M-1)
 			std::cout << ";" << std::endl; 
 		else
 			std::cout << "]" << std::endl ;
 	}
 #endif
-/*		
-	QTime time ;
-	time.start() ;
-*/
+	
 	// 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()(std::min(2*M, N)-1) ;
 	const double sigma_M = svd.singularValues()(0) ;
 
-/*
-	int msec = time.elapsed() ;
-	int sec  = (msec / 1000) % 60 ;
-	int min  = (msec / 60000) % 60 ;
-	int hour = (msec / 3600000) ;
-	std::cout << "<<INFO>> delta took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
-*/
 #ifdef DEBUG
 	std::cout << "<<DEBUG>> SVD = [ " ;
-	for(int i=0; i<std::min(2*d->size(), np+nq); ++i)
+	for(int i=0; i<std::min(2*M, N); ++i)
 	{
 		std::cout << svd.singularValues()(i) << ", " ;
 	}
@@ -270,16 +240,16 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
 #endif
 		return false ;
 	}
-    else if(delta <= 0.0)
+   else if(delta <= 0.0)
 	{
 		delta = 1.0 ;
 	}
 
 
-#ifdef DEBUG
+#ifndef DEBUG
 	std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
 #endif
-	for(int i=0; i<2*d->size(); ++i)	
+	for(int i=0; i<2*M; ++i)	
 	{		
 		ci[i] = ci[i] * delta ; 
 #ifdef DEBUG
@@ -308,7 +278,7 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
 		// Recopy the vector d
 		vec p(np), q(nq);
 		double norm = 0.0 ;
-		for(int i=0; i<np+nq; ++i)
+		for(int i=0; i<N; ++i)
 		{
 			const double v = x[i];
 			norm += v*v ;