Fixing bug in the matlab implementation.

Forcing to use the interior-point-convex algorithm which
produce more stable results. Still, it is not in the bounds.

Commented the normalization in the rational function class,
it seems to introduce some errors (weird).

sources/core/rational_function.cpp

@@ -34,7 +34,7 @@ void rational_function_1d::update(const vec& in_a,
 	a.resize(in_a.size()) ;
 	b.resize(in_b.size()) ;
 
-    const double b0 = in_b[0];
+    const double b0 = 1.0;//in_b[0];
 
     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; }
@@ -201,9 +201,9 @@ double rational_function_1d::p(const vec& x, int i) const
 	std::vector<int> deg = index2degree(i);
 	double res = 1.0;
 	for(int k=0; k<dimX(); ++k)
-	{
-//		res *= pow(x[k], deg[k]) ;
+    {
         res *= pow(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]) ;
+        //res *= pow(x[k], deg[k]) ;
 	}
 
 	return res ;
@@ -388,19 +388,6 @@ void rational_function::load(std::istream& in)
 	for(int k=0; k<dimX(); ++k) {in >> max[k]; }
 	setMax(max);
 
-
-	// Check for the polynomial basis type
-	in >> token;
-   if(token.compare("#BASIS") != 0) 
-	{
-		std::cerr << "<<ERROR>> the file is not specifying the polynomial basis." << std::endl; 
-	}
-	in >> token;
-   if(token.compare("LEGENDRE") != 0) 
-	{
-		std::cerr << "<<ERROR>> the basis is different than LEGENDRE." << std::endl; 
-	}
-
 	vec a(_np), b(_nq);
 	for(int i=0; i<_nY; ++i)
 	{
@@ -432,7 +419,6 @@ void rational_function::save_call(std::ostream& out, const arguments&) const
 	out << "#NQ " << nq << std::endl ;
 	out << "#MIN "; for(int k=0; k<_nX; ++k) { out << _min[k] << " "; } out << std::endl;
 	out << "#MAX "; for(int k=0; k<_nX; ++k) { out << _max[k] << " "; } out << std::endl; 
-	out << "#BASIS LEGENDRE" << std::endl ;
 
 	for(int k=0; k<_nY; ++k)
 	{

sources/core/rational_function.h

@@ -74,7 +74,7 @@ class rational_function_1d : public function
 
 		// Store the coefficients for the moment, I assume
 		// the functions to be polynomials.
-		vec a, b ;
+        vec a, b ;
 } ;
 
 #ifdef TODO

sources/core/vertical_segment.cpp

@@ -100,26 +100,29 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
 		for(int i=0; i<dimY(); ++i)
 			linestream >> v[dimX() + i] ;
 
-		for(int i=0; i<dimY(); ++i)
+        // Check if the data containt a vertical segment around the mean
+        // value.
+        for(int i=0; i<dimY(); ++i)
         {
             double min_dt = 0.0;
             double max_dt = 0.0;
 
-			if(vs[i] == 2) 
-			{
+
+            if(vs[i] == 2)
+            {
                 linestream >> min_dt ;
                 linestream >> max_dt ;
                 min_dt = min_dt-v[dimX()+i];
                 max_dt = max_dt-v[dimX()+i];
-			} 
-			else if(vs[i] == 1)
-			{
-				double dt ;
-				linestream >> dt ;
+            }
+            else if(vs[i] == 1)
+            {
+                double dt ;
+                linestream >> dt ;
                 min_dt = -dt;
                 max_dt =  dt;
-			}
-			else 
+            }
+            else
             {
                 double dt = args.get_float("dt", 0.1f);
                 min_dt = -dt;
@@ -127,16 +130,16 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
             }
 
             if(args.is_defined("dt-relative"))
-               {
-                    v[dimX() + dimY()+i] = v[dimX() + i] * (1.0 + min_dt) ;
-                    v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + max_dt) ;
-                }
-                else
-                {
-                    v[dimX() + dimY()+i] = v[dimX() + i] + min_dt ;
-                    v[dimX() + 2*dimY()+i] = v[dimX() + i] + max_dt ;
-                }
-		}
+            {
+                v[dimX() + dimY()+i]   = v[dimX() + i] * (1.0 + min_dt) ;
+                v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + max_dt) ;
+            }
+            else
+            {
+                v[dimX() + dimY()+i]   = v[dimX() + i] + min_dt ;
+                v[dimX() + 2*dimY()+i] = v[dimX() + i] + max_dt ;
+            }
+        }
 
 		// If data is not in the interval of fit
 		bool is_in = true ;

sources/plugins/rational_fitter_matlab/rational_fitter.cpp

@@ -20,7 +20,11 @@ ALTA_DLL_EXPORT fitter* provide_fitter()
 rational_fitter_matlab::rational_fitter_matlab() 
 {
 	// Create matlab engine
-	if (!(ep = engOpen(""))) 
+#ifdef WIN32
+    if (!(ep = engOpen(NULL)))
+#else
+    if (!(ep = engOpen("matlab -nosplash")))
+#endif
 	{
 		std::cerr << "<ERROR>> can't start MATLAB engine" << std::endl ;
 	}
@@ -89,6 +93,8 @@ void rational_fitter_matlab::set_parameters(const arguments& args)
 	_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) ;	
+
+    _use_matlab = !args.is_defined("use-qpas");
 }
 		
 bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq, rational_function* r) 
@@ -207,7 +213,7 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 
 	if(std::isnan(delta) || (std::abs(delta) == std::numeric_limits<double>::infinity()))
 	{
-#ifdef DEBUG
+#ifdef DEBUG_MATRICES
 		std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
 #endif
 		return false ;
@@ -217,7 +223,7 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 		delta = 1.0 ;
 	}
 
-#ifdef DEBUG
+#ifdef DEBUG_MATRICES
 	std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
 #endif
 	
@@ -226,7 +232,7 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 		ci(i) = ci(i) * delta ; 
 	}
 	
-#ifdef DEBUG
+#ifdef DEBUG_MATRICES
 	std::cout << "CI = " << CI << std::endl << std::endl ;
 #endif
 
@@ -253,7 +259,7 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 	char* output = new char[BUFFER_SIZE+1];
 	output[BUFFER_SIZE] = '\0';
 	engOutputBuffer(ep, output, BUFFER_SIZE) ;
-#ifdef DEBUG
+#ifdef DEBUG_MATRICES
 	engEvalString(ep, "display(H)");
 	std::cout << output << std::endl ;
 	engEvalString(ep, "display(f)");
@@ -264,21 +270,27 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 	std::cout << output << std::endl ;
 #endif
 
-#ifdef USE_MATLAB
-	engEvalString(ep, "[x, fval, flag] = quadprog(H,f,A,b);");
+    // Use Matlab quadratic programming solver
+    if(_use_matlab)
+    {
+        engEvalString(ep, "options = optimset('Algorithm', 'interior-point-convex')");
+        engEvalString(ep, "[x, fval, flag] = quadprog(H, f, A, b, [], [], [], [], [], options);");
 #ifdef DEBUG
-	std::cout << output << std::endl ;
+        std::cout << output << std::endl ;
 #endif
-#else
-    engEvalString(ep, "cd matlab;");
-    engEvalString(ep, "[x, err] = qpas(H,f,A,b);");
+    }
+    else // Use QPAS solver
+    {
+        engEvalString(ep, "cd matlab;");
+        engEvalString(ep, "[x, err] = qpas(H,f,A,b);");
 #ifdef DEBUG
-	std::cout << output << std::endl ;
+        std::cout << output << std::endl ;
 #endif
-    engEvalString(ep, "flag = err == 0.0;");
-    engEvalString(ep, "cd ..;");
-#endif
-#ifdef DEBUG
+        engEvalString(ep, "flag = err == 0.0;");
+        engEvalString(ep, "cd ..;");
+    }
+
+#ifdef DEBUG_MATRICES
 	engEvalString(ep, "display(x)");
 	std::cout << output << std::endl ;
 	engEvalString(ep, "display(flag)");
@@ -307,26 +319,24 @@ bool rational_fitter_matlab::fit_data(const vertical_segment* d, int np, int nq,
 				return false ;
 			}
 
-			double  total = 0.0;
 			double* val = (double*)mxGetData(x) ;
 			vec a(np), b(nq);
 			for(int i=0; i<N; ++i)
 			{
-				total += val[i]*val[i] ;
 				if(i < np)
 				{
-					a[i] =val[i] ;
+                    a[i] = val[i] ;
 				}
 				else
 				{
-					b[i] = val[i] ;
+                    b[i-np] = val[i] ;
 				}
 			}
 			r->update(a, b) ;
 
 			mxDestroyArray(x);
 			mxDestroyArray(flag);
-			return total > 0.0 ;
+            return true ;
 		}
 		else
 		{

sources/plugins/rational_fitter_matlab/rational_fitter.h

@@ -13,6 +13,9 @@
 #include <core/rational_function.h>
 #include <core/vertical_segment.h>
 
+//#define USE_MATLAB
+#define DEBUG
+
 class rational_fitter_matlab : public fitter
 {
 	public: // methods
@@ -42,5 +45,8 @@ class rational_fitter_matlab : public fitter
 		int _max_np, _max_nq ;
 		int _min_np, _min_nq ;
 		Engine *ep;
+
+        // Decide which quadratic solver to use. Matlab's or the QPAS solver
+        bool _use_matlab;
 } ;
 

sources/softs/brdf2gnuplot/main.cpp

@@ -101,7 +101,7 @@ int main(int argc, char** argv)
 			{
 				if(plot_error)
 				{
-					file << (v[d->dimX() + u] - y2[u])/v[d->dimX()+u] << "\t" ;
+                    file << (v[d->dimX() + u] - y2[u]) << "\t" ;
 				}
 				else if(linear_plot)
 				{

sources/softs/generate_data/main.cpp

@@ -13,7 +13,7 @@ int main(int argc, char** argv)
 
 	std::cout.precision(10);
 
-	int nbx = 1000;
+    int nbx = 1000;
 	int nby = 10;
 	int nbz = 10;
 	if(args.is_defined("nbx"))
@@ -27,7 +27,7 @@ int main(int argc, char** argv)
 	{
 		f << "#DIM 1 1" << std::endl ;
 		f << "#PARAM_IN UNKNOWN" << std::endl;
-		//f << "#VS 2" << std::endl;
+        f << "#VS 2" << std::endl;
 		for(int i=0; i<nbx; ++i)
 		{
 			const float x  = i / (float)nbx ;
@@ -35,7 +35,7 @@ int main(int argc, char** argv)
 			//const float y  = 1000.0f * exp(- x*x) * x*x + 00.1 * exp(-100.0 * xp*xp)  *x*x*x + 0.1 ;
 			const float y = (1.0) / (1.0E-10 + x);
 
-			f << x << "\t" << y << /*"\t" << y*0.9f << "\t" << y*1.1f <<*/ std::endl ;
+            f << x << "\t" << y << "\t" << y*0.9f << "\t" << y*1.1f << std::endl ;
 		}
 	}
 	else if(k == K++)
@@ -43,13 +43,14 @@ int main(int argc, char** argv)
 		std::cout << "<<INFO>> output 1.0E4 x^2 e^{- x^2) + 0.1 x^3 e^{-100.0 (x - 9.0)^2} + 0.1" << std::endl;
 		f << "#DIM 1 1" << std::endl ;
 		f << "#PARAM_IN UNKNOWN" << std::endl;
+        f << "#VS 2" << std::endl;
 		for(int i=0; i<nbx; ++i)
 		{
 			const float x = 10.0 * i / (float)nbx ;
 			const float xp = (x - 9);
 			const float y  = 1000.0f * exp(- x*x) * x*x + 00.1 * exp(-100.0 * xp*xp)  *x*x*x + 0.1 ;
 
-			f << x << "\t" << y << std::endl ;
+            f << x << "\t" << y << "\t" << y*(0.9) << "\t" << y*(1.1)<< std::endl ;
 		}
 	}
 	else if(k == K++)