rational_fitter.cpp 7.44 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
#include "rational_fitter.h"

#include <Eigen/SVD>
#include <Array.hh>
#include <QuadProg++.hh>

#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>

14 15
#include <QTime>

16

17 18
data* rational_fitter_quadprog::provide_data() const
{
Laurent Belcour's avatar
Laurent Belcour committed
19
	return new vertical_segment() ;
20 21 22 23 24 25 26
}

function* rational_fitter_quadprog::provide_function() const 
{
	return new rational_function() ;
}

27
rational_fitter_quadprog::rational_fitter_quadprog() : QObject()
28 29
{
}
30
rational_fitter_quadprog::~rational_fitter_quadprog() 
31 32 33
{
}

34
bool rational_fitter_quadprog::fit_data(const data* dat, function* fit)
35
{
36
	rational_function* r = dynamic_cast<rational_function*>(fit) ;
Laurent Belcour's avatar
Laurent Belcour committed
37
	const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
38 39 40 41 42 43 44 45 46 47
	if(r == NULL || d == NULL)
	{
		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
		return false ;
	}

	// 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()) ;
48 49
	r->setMin(d->min()) ;
	r->setMax(d->max()) ;
50 51 52 53 54

	std::cout << "<<INFO>> np in  [" << _min_np << ", " << _max_np 
	          << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;


55
	int temp_np = _min_np, temp_nq = _min_nq ;
56
	while(temp_np <= _max_np || temp_nq <= _max_nq)
57
	{
58 59
		QTime time ;
		time.start() ;
60
		
61
		if(fit_data(d, temp_np, temp_nq, r))
62
		{
63 64 65 66 67 68 69
			int msec = time.elapsed() ;
			int sec  = (msec / 1000) % 60 ;
			int min  = (msec / 60000) % 60 ;
			int hour = (msec / 3600000) ;
			std::cout << "<<INFO>> got a fit using np = " << temp_np << " & nq =  " << temp_nq << "      " << std::endl ;
			std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;

70 71 72
			return true ;
		}

73
		std::cout << "<<INFO>> fit using np = " << temp_np << " & nq =  " << temp_nq << " failed\r"  ;
74 75
		std::cout.flush() ;

76
		if(temp_np <= _max_np)
77 78 79
		{
			++temp_np ;
		}
80
		if(temp_nq <= _max_nq)
81 82 83 84 85 86 87 88
		{
			++temp_nq ;
		}
	}

	return false ;
}

89
void rational_fitter_quadprog::set_parameters(const arguments& args)
90 91 92 93 94 95 96 97
{
	_max_np = args.get_float("np", 10) ;
	_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) ;
}
		

Laurent Belcour's avatar
Laurent Belcour committed
98
bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int nq, rational_function* r) 
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
{

	// Multidimensional coefficients
	std::vector<double> Pn ; Pn.reserve(d->dimY()*np) ;
	std::vector<double> Qn ; Qn.reserve(d->dimY()*nq) ;

	for(int j=0; j<d->dimY(); ++j)
	{
		if(!fit_data(d, np, nq, j, r))
			return false ;

		for(int i=0; i<np; ++i) { Pn.push_back(r->getP(i)) ; }
		for(int i=0; i<nq; ++i) { Qn.push_back(r->getQ(i)) ; }
	}

	r->update(Pn, Qn) ;
	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 return a ration BRDF function and a boolean
Laurent Belcour's avatar
Laurent Belcour committed
122
bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int nq, int ny, rational_function* rf) 
123 124
{
	rational_function* r = dynamic_cast<rational_function*>(rf) ;
Laurent Belcour's avatar
Laurent Belcour committed
125
	const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
126 127 128 129 130 131
	if(r == NULL || d == NULL)
	{
		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
		return false ;
	}

132 133 134 135
	// 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() ;

136 137
	// Matrices of the problem
	QuadProgPP::Matrix<double> G (0.0, np+nq, np+nq) ;
138
	QuadProgPP::Vector<double> g (0.0, np+nq) ;
139 140
	QuadProgPP::Matrix<double> CI(0.0, np+nq, 2*d->size()) ;
	QuadProgPP::Vector<double> ci(0.0, 2*d->size()) ;
141 142
	QuadProgPP::Matrix<double> CE(0.0, np+nq, 0) ;
	QuadProgPP::Vector<double> ce(0.0, 0) ;
143

144 145
	Eigen::MatrixXd eCI(np+nq, 2*d->size()) ;

146 147 148 149
	// Select the size of the result vector to
	// be equal to the dimension of p + q
	for(int i=0; i<np+nq; ++i)
	{
150
		G[i][i] = 1.0 ; 
151 152 153 154 155 156 157 158 159 160 161
	}
	
	// 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 ;

162
		vec xi = d->get(i) ;
163
/*		for(int k=0; k<d->dimX(); ++k)
164 165 166
		{
			xi[k] /= dmax[k] ;
		}
167
*/
168 169 170 171 172 173 174 175 176
		// 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)
			{
177
				const double pi = r->p(xi, j) ;
178 179
				a0_norm += pi*pi ;
				a1_norm += pi*pi ;
180 181
				CI[j][i] =  pi ;
				CI[j][i+d->size()] = -pi ;
182 183 184 185

				// Updating Eigen matrix
				eCI(j,i) = pi ;
				eCI(j,i+d->size()) = -pi ;
186 187 188 189 190 191 192
			}
			// Filling the q part
			else
			{
				vec yl, yu ; 
				d->get(i, yl, yu) ;

193
				const double qi = r->q(xi, j-np) ;
194 195
				a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
				CI[j][i] = -yl[ny] * qi ;
196
				
197 198
				a1_norm += qi*qi * (yu[ny]*yu[ny]) ;
				CI[j][i+d->size()] = yu[ny] * qi ;
199 200
				
				// Updating Eigen matrix
201 202
				eCI(j,i) = -yl[ny] * qi ;
				eCI(j,i+d->size()) = yu[ny] * qi ;
203 204 205 206 207
			}
		}
	
		// Set the c vector, will later be updated using the
		// delta parameter.
208 209
		ci[i] = -sqrt(a0_norm) ;
		ci[i+d->size()] = -sqrt(a1_norm) ;
210 211
	}
#ifdef DEBUG
212
	std::cout << "CI = [" ;
213 214 215
	for(int j=0; j<d->size()*2; ++j)
	{
		for(int i=0; i<np+nq; ++i)
216 217 218 219 220 221 222 223
		{
			std::cout << CI[i][j] ;
			if(i != np+nq-1) std::cout << ", ";
		}
		if(j != d->size()*2-1)
			std::cout << ";" << std::endl; 
		else
			std::cout << "]" << std::endl ;
224 225 226 227
	}
#endif
	// Update the ci column with the delta parameter
	// (See Celis et al. 2007 p.12)
228
	Eigen::JacobiSVD<Eigen::MatrixXd> svd(eCI);
229 230
	const double sigma_m = svd.singularValues()(std::min(2*d->size(), np+nq)-1) ;
	const double sigma_M = svd.singularValues()(0) ;
231

232 233 234 235 236 237 238 239
#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
240
	
241 242 243 244 245 246 247 248 249 250 251 252 253
	double 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 ;
	}

254

255
#ifndef DEBUG
256 257 258 259
	std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
	for(int i=0; i<2*d->size(); ++i)	
	{		
260
		ci[i] = ci[i] * delta ; 
261 262 263
#ifdef DEBUG
		std::cout << ci[i] << "\t" ;
#endif
264
	}
265 266
#ifdef DEBUG
	std::cout << std::endl << std::endl ;
267

268 269
	std::cout << eCI << std::endl << std::endl ;
#endif
270 271 272
	// Compute the solution
	QuadProgPP::Vector<double> x;
	double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, x);
273 274


275
	bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
276 277
	for(int i=0; i<np+nq; ++i)
	{
278
		const double v = x[i];
279 280 281 282 283 284 285
		solves_qp = solves_qp && !std::isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
	}

	if(solves_qp)
	{
		// Recopy the vector d
		std::vector<double> p, q;
286
		double norm = 0.0 ;
287 288
		for(int i=0; i<np+nq; ++i)
		{
289
			const double v = x[i];
290
			norm += v*v ;
291 292
			if(i < np)
			{
293
				p.push_back(v) ;
294 295 296
			}
			else
			{
297
				q.push_back(v) ;
298 299 300
			}
		}

301
		r->update(p, q);
302
		std::cout << "<<INFO>> got solution " << *r << std::endl ;
303
		
304
		return norm > 0.0;
305 306
	}
	else
307

308 309 310 311 312
	{
		return false; 
	}
}

313
Q_EXPORT_PLUGIN2(rational_fitter_quadprog, rational_fitter_quadprog)