fitter.cpp 5.57 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#include "fitter.h"

#include <nlopt.h>

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

#include <QTime>

#include <core/common.h>

ALTA_DLL_EXPORT fitter* provide_fitter()
{
18
	return new nonlinear_fitter_nlopt();
19 20
}

21
void print_nlopt_error(nlopt_result res, const std::string& string)
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
{
	if(res == NLOPT_FAILURE)
	{
		std::cerr << "<<ERROR>> generic failure for \"" << string << "\"" << std::endl;
	}
	else if(res == NLOPT_INVALID_ARGS)
	{
		std::cerr << "<<ERROR>> invalid arguments for \"" << string << "\"" << std::endl;
	}
	else if(res == NLOPT_OUT_OF_MEMORY)
	{
		std::cerr << "<<ERROR>> not enough memory for \"" << string << "\"" << std::endl;
	}
}

// The parameter of the function _f should be set prior to this function
// call. If not it will produce undesirable results.
void df(double* fjac, const nonlinear_function* f, const data* d)
{
41 42 43
	// Clean memory
	memset(fjac, 0.0, f->nbParameters()*sizeof(double));

44 45 46 47 48 49 50 51
	// Each constraint is of the form data point * color channel
	for(int s=0; s<d->size(); ++s)
	{
		vec xi = d->get(s);

		// Get the jacobian of the function at position x_i for the current
		// set of parameters (set prior to function call)
		vec _jac = f->parametersJacobian(xi);
52

53 54 55 56 57 58 59 60 61
		vec _di = vec(f->dimY());
		for(int i=0; i<f->dimY(); ++i)
		{
			_di[i] = xi[f->dimX() + i];
		}

		// Should add the resulting vector completely
		vec _y = (*f)(xi) - _di;

62 63
		// Fill the columns of the matrix
		for(int j=0; j<f->nbParameters(); ++j)
64
		{
65 66 67
			// For each output channel, update the subpart of the
			// vector row
			for(int i=0; i<f->dimY(); ++i)
68
			{
69
				fjac[j] += 2 * _y[i] * _jac[i*f->nbParameters() + j];
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
			}
		}
	}
}

double f(unsigned n, const double* x, double* dy, void* dat)
{
	nonlinear_function* _f = (nonlinear_function*)(((void**)dat)[0]);
	const data* _d = (const data*)(((void**)dat)[1]);

	// Update the parameters vector
	vec _p(_f->nbParameters());
	for(int i=0; i<_f->nbParameters(); ++i) { _p[i] = x[i]; }
	_f->setParameters(_p);

	// Create the result vector
	double y = 0.0;

	// Each constraint is of the form data point * color channel
	for(int s=0; s<_d->size(); ++s)
	{
		vec xi = _d->get(s);
		vec _di = vec(_f->dimY());
		for(int i=0; i<_f->dimY(); ++i)
		{
			_di[i] = xi[_f->dimX() + i];
		}

		// Should add the resulting vector completely
		vec _y = (*_f)(xi) - _di;
		for(int i=0; i<_f->dimY(); ++i)
		{
			y += pow(_y[i], 2);
		}
	}

	if(dy != NULL)
	{
		df(dy, _f, _d);
	}

	return y;
}

nonlinear_fitter_nlopt::nonlinear_fitter_nlopt() 
{
}
nonlinear_fitter_nlopt::~nonlinear_fitter_nlopt() 
{
}

bool nonlinear_fitter_nlopt::fit_data(const data* d, function* fit, const arguments &args)
{
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
	// I need to set the dimension of the resulting function to be equal
	// to the dimension of my fitting problem
	fit->setDimX(d->dimX()) ;
	fit->setDimY(d->dimY()) ;
	fit->setMin(d->min()) ;
	fit->setMax(d->max()) ;

	// Convert the function and bootstrap it with the data
	if(dynamic_cast<nonlinear_function*>(fit) == NULL)
	{
		std::cerr << "<<ERROR>> the function is not a non-linear function" << std::endl;
		return false;
	}
	nonlinear_function* nf = dynamic_cast<nonlinear_function*>(fit);
	nf->bootstrap(d, args);
138 139

#ifndef DEBUG
140
	std::cout << "<<DEBUG>> number of parameters: " << nf->nbParameters() << std::endl;
141
#endif
142 143 144 145 146 147 148 149 150 151 152
	if(nf->nbParameters() == 0)
	{
		return true;
	}

	// the following starting values provide a rough fit is the bootstrap flag is
	// enabled
	vec p = nf->parameters();

	// Check the optimizer to use 
	nlopt_algorithm algorithm;
153
	nlopt_result res;
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
	std::string optimizerName = args["nlopt-optimizer"];
	if(optimizerName == "COBYLA")
	{
		algorithm = NLOPT_LN_COBYLA;
	}
	else if(optimizerName == "BOBYQA")
	{
		algorithm = NLOPT_LN_BOBYQA;
	}
	else if(optimizerName == "NEWUOA")
	{
		algorithm = NLOPT_LN_NEWUOA_BOUND;
	}
	else if(optimizerName == "PRAXIS")
	{
		algorithm = NLOPT_LN_PRAXIS;
	}
171
    else if(optimizerName == "Nelder-Mead")
172 173 174 175 176 177 178
	{
		algorithm = NLOPT_LN_NELDERMEAD;
	}
	else if(optimizerName == "Sbplx")
	{
		algorithm = NLOPT_LN_SBPLX;
	}
179 180 181 182
	else if(optimizerName == "controlled-random-search")
	{
		algorithm = NLOPT_GN_CRS2_LM;
	}
183 184 185 186
	else if(optimizerName == "MMA")
	{
		algorithm = NLOPT_LD_MMA; 
	}
187 188 189 190
	else if(optimizerName == "SQP")
	{
		algorithm = NLOPT_LD_SLSQP;
	}
191 192
	else
	{
193
		algorithm = NLOPT_LD_SLSQP;
194 195 196 197 198 199 200 201 202 203
	}

	// Create the optimizer
	nlopt_opt opt = nlopt_create(algorithm, nf->nbParameters());
	if(opt == NULL)
	{
		std::cerr << "<<ERROR>> unable to create the optimizer" << std::endl;
		return false;
	}

204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

	// Set the bounds to the parameters
	vec p_min = nf->getParametersMin();
	vec p_max = nf->getParametersMax();
	nlopt_set_lower_bounds(opt, &p_min[0]);
	nlopt_set_upper_bounds(opt, &p_max[0]);

	// Set the stopping criterion to a maximum number of evaluation
	if(args.is_defined("nlop-max-num-iterations"))
	{
		res = nlopt_set_maxeval(opt, args.get_int("nlop-max-num-iterations", 10));
		if(res < 0) { print_nlopt_error(res, "nlopt_set_maxeval"); }
	}
	nlopt_set_xtol_rel(opt, 1e-4);

219 220 221 222
	// Create the problem
	void* dat[2];
	dat[0] = (void*)nf;
	dat[1] = (void*)d;
223
	res = nlopt_set_min_objective(opt, f, dat);
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
	if(res < 0)
	{
		print_nlopt_error(res, "nlopt_set_min_objective");
		return false;
	}


	// Launch minimization
	double f_end;
	res = nlopt_optimize(opt, &p[0], &f_end);

	if(res > 0)
	{
		std::cout << "<<INFO>> optimized distance: " << f_end << std::endl;
		std::cout << "<<INFO>> found parameters: " << p << std::endl;
		nf->setParameters(p);
	}
	else
	{
		print_nlopt_error(res, "nlopt_optimize");
	}

	nlopt_destroy(opt);
	return res > 0;
248 249 250 251 252
}

void nonlinear_fitter_nlopt::set_parameters(const arguments& args)
{
}