Commit 2385eb6f authored by Laurent Belcour's avatar Laurent Belcour

Working DCA, without Papamarkos boostrap

parent 39cc8cb1
......@@ -35,19 +35,19 @@ bool rational_fitter_dca::fit_data(const data* dat, function* fit)
{
rational_function* r = dynamic_cast<rational_function*>(fit) ;
const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
if(r == NULL || d == NULL || ep == NULL)
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// Create matlab engine
if (!(ep = engOpen("")))
{
std::cerr << "<ERROR>> can't start MATLAB engine" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
......@@ -57,20 +57,20 @@ bool rational_fitter_dca::fit_data(const data* dat, function* fit)
r->setMax(d->max()) ;
QTime time ;
time.start() ;
if(fit_data(d, _max_np, _max_nq, r))
{
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 = " << _max_np << " & nq = " << _max_nq << " " << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
QTime time ;
time.start() ;
return true ;
}
if(fit_data(d, _max_np, _max_nq, r))
{
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 = " << _max_np << " & nq = " << _max_nq << " " << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
return true ;
}
engClose(ep);
return false ;
......@@ -80,30 +80,28 @@ void rational_fitter_dca::set_parameters(const arguments& args)
{
_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) ;
}
double distance(const rational_function* f, const data* d)
{
double distance = 0.0;
for(int i=0; i<d->size(); ++i)
{
vec xi = d->get(i) ;
vec y = f->value(xi);
double current_d = 0.0;
for(int j=0; j<d->dimY(); ++j)
{
double diff = y[j] - xi[d->dimX()+j];
current_d += diff;
}
current_d = std::abs(current_d);
if(current_d > distance)
distance = current_d;
}
return distance;
double distance = 0.0;
for(int i=0; i<d->size(); ++i)
{
vec xi = d->get(i) ;
vec y = f->value(xi);
double current_d = 0.0;
for(int j=0; j<d->dimY(); ++j)
{
double diff = y[j] - xi[d->dimX()+j];
current_d += diff;
}
current_d = std::abs(current_d);
if(current_d > distance)
distance = current_d;
}
return distance;
}
// Bootstrap the DCA algorithm with the Papamarkos fitting
......@@ -111,12 +109,12 @@ double distance(const rational_function* f, const data* d)
// \todo Finish the Papamarkos implementation
void rational_fitter_dca::bootstrap(const data* d, int np, int nq, rational_function* fit, double& delta)
{
vec p(np*d->dimY());
vec q(nq*d->dimY());
q[0] = 1.0;
vec p(np*d->dimY());
vec q(nq*d->dimY());
q[0] = 1.0;
fit->update(p, q);
delta = distance(fit, d);
fit->update(p, q);
delta = distance(fit, d);
}
// dat is the data object, it contains all the points to fit
......@@ -130,29 +128,59 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
int M = d->size() ;
int nY = d->dimY();
// Bootstrap the delta and rational function using the Papamarkos
// algorithm.
double delta = 0.0;
bootstrap(d, np, nq, r, delta);
// Bootstrap the delta and rational function using the Papamarkos
// algorithm.
double delta = 0.0;
bootstrap(d, np, nq, r, delta);
#ifndef DEBUG
std::cout << "<<DEBUG>> delta value after boostrap: " << delta << std::endl;
std::cout << "<<DEBUG>> delta value after boostrap: " << delta << std::endl;
std::cout << "<<DEBUG>> r: " << *r << std::endl;
#endif
// Create the MATLAB defintion of objects
// MATLAB defines a linear prog as
// min f' x with A x <= b
//
mxArray *f, *A, *b, *x, *flag;
f = mxCreateDoubleMatrix(N*nY, 1, mxREAL);
A = mxCreateDoubleMatrix( 2*M*nY, N*nY, mxREAL);
b = mxCreateDoubleMatrix( 2*M*nY, 1, mxREAL);
// Matrices of the problem in Eigen format
mxArray *f, *A, *b, *x, *flag, *u, *l;
f = mxCreateDoubleMatrix( N*nY, 1, mxREAL);
A = mxCreateDoubleMatrix( 4*M*nY, N*nY, mxREAL);
b = mxCreateDoubleMatrix( 4*M*nY, 1, mxREAL);
u = mxCreateDoubleMatrix( N*nY, 1, mxREAL);
l = mxCreateDoubleMatrix( N*nY, 1, mxREAL);
// Matrices of the problem in Eigen format
Eigen::VectorXd g (nY*N) ;
Eigen::MatrixXd CI(2*M*nY, N*nY) ;
Eigen::VectorXd ci(2*M*nY) ;
Eigen::MatrixXd CI(4*M*nY, N*nY) ;
Eigen::VectorXd ci(4*M*nY) ;
Eigen::VectorXd li(N*nY) ;
Eigen::VectorXd ui(N*nY) ;
// Updating the upper and lower values for the solution
for(int i=0; i<N*nY; ++i)
{
if(i<np*nY)
{
ui(i) = 100.0; // std::numeric_limits<double>::max();
li(i) = -100.0; //-std::numeric_limits<double>::max();
}
else if(i<(np+nq)*nY)
{
ui(i) = 1.0;
li(i) = -1.0;
}
else
{
ui(i) = 100.0; // std::numeric_limits<double>::max();
li(i) = -100.0; //-std::numeric_limits<double>::max();
}
}
memcpy((void *)mxGetPr(u), (void *)ui.data(), N*nY*sizeof(double));
memcpy((void *)mxGetPr(l), (void *)li.data(), N*nY*sizeof(double));
engPutVariable(ep, "u", u);
engPutVariable(ep, "l", l);
double delta_k = delta;
double delta_k = delta;
// Loop until you get a converge solution \delta > \delta_k
// \todo add the correct looping condition
......@@ -173,21 +201,21 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
// matrix
for(int i=0; i<M; ++i)
{
//
//
vec xi = d->get(i) ;
// For each input data i \in M, the following constraints have to
// be fulfilled:
// [ f_i + \delta_k] * q_i - p_i + qk_i \delta >= 0
// [-f_i + \delta_k] * q_i + p_i + qk_i \delta >= 0
// pi >= 0
// qi >= 0
// p >= 0
// q > 0
//
// The constraints matrix has the following form
// [-p_j(x_i) ..-p_j(x_i), [ f_i+\delta_k]*q_j(x_i) ..[ fi+\delta_k]*q_j(x_i), qk(x_i)]
// [ p_j(x_i) .. p_j(x_i), [-f_i+\delta_k]*q_j(x_i) ..[-fi+\delta_k]*q_j(x_i), qk(x_i)]
//
for(int j=0; j<N; ++j)
for(int j=0; j<N; ++j)
{
// Filling the p part
if(j<np)
......@@ -197,8 +225,10 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
// Updating Eigen matrix
for(int y=0; y<nY; ++y)
{
CI(2*(nY*i + y)+0, nY*j + y) = pi ;
CI(2*(nY*i + y)+1, nY*j + y) = -pi ;
CI(2*(nY*i + y)+0, nY*j + y) = pi ;
CI(2*(nY*i + y)+1, nY*j + y) = -pi ;
CI(2*M*nY + 2*(nY*i+y) + 0, nY*j+y) = -pi ;
CI(2*M*nY + 2*(nY*i+y) + 1, nY*j+y) = 0.0 ;
}
}
// Filling the q part
......@@ -209,18 +239,22 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
// Updating Eigen matrix
for(int y=0; y<nY; ++y)
{
CI(2*(nY*i + y)+0, nY*j + y) = -(delta_k+xi[d->dimX()+y]) * qi ;
CI(2*(nY*i + y)+1, nY*j + y) = -(delta_k-xi[d->dimX()+y]) * qi ;
CI(2*(nY*i + y)+0, nY*j + y) = -(delta_k+xi[d->dimX()+y]) * qi ;
CI(2*(nY*i + y)+1, nY*j + y) = -(delta_k-xi[d->dimX()+y]) * qi ;
CI(2*M*nY + 2*(nY*i+y) + 0, nY*j + y) = 0.0 ;
CI(2*M*nY + 2*(nY*i+y) + 1, nY*j + y) = - qi ;
}
}
else
{
// Last column of the constraint matrix
vec qk = r->q(xi) ;
vec qk = r->q(xi) ;
for(int y=0; y<nY; ++y)
{
CI(2*(nY*i + y)+0, nY*j + y) = -qk[y] ;
CI(2*(nY*i + y)+1, nY*j + y) = -qk[y] ;
CI(2*(nY*i + y)+0, nY*j + y) = -qk[y] ;
CI(2*(nY*i + y)+1, nY*j + y) = -qk[y] ;
CI(2*M*nY + 2*(nY*i+y) + 0, nY*j + y) = 0.0 ;
CI(2*M*nY + 2*(nY*i+y) + 1, nY*j + y) = 0.0 ;
}
}
}
......@@ -229,22 +263,22 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
// Set the c vector
for(int y=0; y<nY; ++y)
{
ci(2*(nY*i+y)+0) = 0.0 ;
ci(2*(nY*i+y)+1) = 0.0 ;
ci(2*(nY*i+y)+0) = 0.0 ;
ci(2*(nY*i+y)+1) = 0.0 ;
ci(2*M*nY + 2*(nY*i) + 0) = 0.0 ;
ci(2*M*nY + 2*(nY*i) + 1) = 0.0 ;
}
}
std::cout << CI << std::endl << std::endl;
// Copy the data to matlab and execute the linear program
//
memcpy((void *)mxGetPr(f), (void *) g.data(), N*nY*sizeof(double));
memcpy((void *)mxGetPr(A), (void *)CI.data(), (2*M*nY)*N*nY*sizeof(double));
memcpy((void *)mxGetPr(b), (void *)ci.data(), (2*M*nY)*sizeof(double));
memcpy((void *)mxGetPr(f), (void *) g.data(), N*nY*sizeof(double));
memcpy((void *)mxGetPr(A), (void *)CI.data(), 4*M*nY*N*nY*sizeof(double));
memcpy((void *)mxGetPr(b), (void *)ci.data(), 4*M*nY*sizeof(double));
engPutVariable(ep, "f", f);
engPutVariable(ep, "A", A);
engPutVariable(ep, "b", b);
engPutVariable(ep, "f", f);
engPutVariable(ep, "A", A);
engPutVariable(ep, "b", b);
char* output = new char[BUFFER_SIZE+1];
output[BUFFER_SIZE] = '\0';
......@@ -256,9 +290,13 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
std::cout << output << std::endl ;
engEvalString(ep, "display(b)");
std::cout << output << std::endl ;
engEvalString(ep, "display(u)");
std::cout << output << std::endl ;
engEvalString(ep, "display(l)");
std::cout << output << std::endl ;
#endif
engEvalString(ep, "[x, fval, flag] = linprog(f,A,b);");
engEvalString(ep, "[x, fval, flag] = linprog(f,A,b,[],[], l, u);");
#ifndef DEBUG
std::cout << output << std::endl ;
#endif
......@@ -281,7 +319,7 @@ bool rational_fitter_dca::fit_data(const data* d, int np, int nq, rational_funct
{
mxDestroyArray(x);
mxDestroyArray(flag);
#ifdef DEBUG
std::cerr << "<<ERROR>> flag is not equal to 1" << std::endl ;
#endif
......
......@@ -52,7 +52,6 @@ class rational_fitter_dca : public QObject, public fitter
// min and Max usable np and nq values for the fitting
//
int _max_np, _max_nq ;
int _min_np, _min_nq ;
Engine *ep;
} ;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment