This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
| Other format: | [Raw text] | |
On Tue, 4 Jan 2005 18:59:58 +0100
Andrej Prsa <andrej.prsa@guest.arnes.si> wrote:
> is being minimized. I never have a fully constrained problem, but at
> least several dimensions may be bounded and that's what I'm after
> here.
It is quite possible that I did not understand your problem, so forgive
me if what I'm going to say will result totally off topic.
Judging from what you wrote and the example you gave to the mailing
list, it seems that you want to change the value of the function
argument in order to mimic some kind of constraint on the function
domain. This things are dangerous.
As I mentioned before, there is a quite straightforward, and definitely
safer way in which a bounded minimization problem can be transformed in
an unbounded one. It is by means of change of variables. Let me give you
an example. Assume that you want x* solution of
x* = argmin{ f(x) } with x \in [-1,1]
Notice that the domain here is a closed interval. Of course this is a
constrained problem. But consider instead the problem
y* = argmin{ f(cos(y) } with y in (-\infty,+\infty)
inspecting the first and second order conditions, you can immediately
see that the extremant y* of this second problem are of two kinds:
1) points y^* such that x^*=arcsin(y^*) is interior to [-1,1] and
is an interior solution of the first problem
2) point y^* such that arcsin(y^*) is on the boundary of the interval of
the first problem, i.e. -1 or 1, and is a boundary solution of the
first problem
You can use an unconstrained optimization algorithm to find out the
set of solutions {y^*}, and then transform back this set to {x^*}.
Now, quite probably what I said above sound trivial to you. Then, let me
conclude by attaching a piece of software that is, if possible,
still more trivial. It is an interface to gsl unbounded optimization
routines to compute constrained optimization on closed, open, bounded
and semi-bounded intervals. It implements different changes of variables
for the different cases. I use it in a package for the maximum
likelihood estimation of the power exponential function (subbotools)
which you can find on my web page (if you want to see an application
of the function). Don't expect anything sophisticated. It's a simple,
minimalistic piece of code. A single function with NAG-like bloated
list of arguments.
I thought that, maybe, someone else could be interested, so I decided to
post it to the mailing list. Forgive me for the size of the resulting
message.
Best,
G.
--
Giulio Bottazzi PGP Key ID:BAB0A33F
/*
multimin.c (ver. 0.9.5) -- Interface to GSL multidim. minimization
Copyright (C) 2002-2003 Giulio Bottazzi
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
(version 2) as published by the Free Software Foundation;
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/*
multimin is an interface to the various GSL minimization
routines. When invoked, all the information necessary to perform the
minimization are passed as formal parameters. This generate a pretty
long, FORTRAN-like, list of parameters. This approach allows, however,
to black-box, as far as possible, the interior functioning of the
routine from the rest of the program.
Let's analyse the calling convention in details:
multimin(size_t n,double *x,double *fun,
unsigned *type,double *xmin,double *xmax,
void (*f) (const size_t,const double *,void *,double *),
void (* df) (const size_t,const double *, void *,double *),
void (* fdf) (const size_t,const double *, void *,double *,double *),
void *fparams,
const struct multimin_params oparams)
where
--------------------------------------------------------------------
n
INPUT: dimension of the problem, number of independent variables of
the function.
--------------------------------------------------------------------
x
INPUT: pointer to an array of n values x[0],...x[n-1] containing the
initial estimate of the minimum point
OUTPUT: contains the final estimation of the minimum position
--------------------------------------------------------------------
type
a pointer to an array of integer type[1],...,type[n-1] describing the
boundary conditions for the different variables. The problem is solved
as an uncostrained one on a suitably trasformed variable y. Possible
values are:
Interval: Trasformation:
0 unconstrained x=y
1 semi-closed right half line [ xmin,+infty ) x=xmin+y^2
2 semi-closed left half line ( -infty,xmax ] x=xmax-y^2
3 closed interval [ xmin,xmax ] x=SS+SD*sin(y)
4 open right half line ( xmin,+infty ) x=xmin+exp(y)
5 open left half line ( -infty,xmax ) x=xmax-exp(y)
6 open interval ( xmin,xmax ) x=SS+SD*tanh(y)
where SS=.5(xmin+xmax) SD=.5(xmax-xmin)
There are also other UNSUPPORTED trasformations used in various test
7 open interval ( xmin,xmax ) x=SS+SD*(1+y/sqrt(1+y^2))
8 open right half line ( xmin,+infty ) x=xmin+.5*(y+sqrt(1+y^2))
--------------------------------------------------------------------
xmin
xmax
pointers to arrays of double containing respectively the lower and
upper boundaries of the different variables. For a given variable,
only the values that are implied by the type of constraints, defined
as in *type, are actually inspected.
--------------------------------------------------------------------
f
f calculates the objective function at a specified point x. Its
specification is
void (*f) (const size_t n, const double *x,void *fparams,double *fval)
n
INPUT: the number of variables
x
INPUT:the point at which the function is required
fparams
pointer to a structure containing parameters required by the
function. If no external parameter are required it can be set to
NULL.
fval
OUTPUT: the value of the objective function at the current point
x.
--------------------------------------------------------------------
df
df calculates the gradient of the objective function at a specified
point x. Its specification is
void (*df) (const size_t n, const double *x,void *fparams,double *grad)
n
INPUT: the number of variables
x
INPUT:the point at which the function is required
fparams
pointer to a structure containing parameters required by the
function. If no external parameter are required it can be set to
NULL.
grad
OUTPUT: the values of the gradient of the objective function at
the current point x are stored in grad[0],...,grad[n-1].
--------------------------------------------------------------------
fdf
fdf calculates the value and the gradient of the objective function at
a specified point x. Its specification is
void (*df) (const size_t n, const double *x,void *fparams,double *fval,double *grad)
n
INPUT: the number of variables
x
INPUT:the point at which the function is required
fparams
pointer to a structure containing parameters required by the
function. If no external parameter are required it can be set to
NULL.
fval
OUTPUT: the value of the objective function at the current point
x.
grad
OUTPUT: the values of the gradient of the objective function at
the current point x are stored in grad[0],...,grad[n-1].
--------------------------------------------------------------------
fparams
pointer to a structure containing parameters required by the
function. If no external parameter are required it can be set to NULL.
--------------------------------------------------------------------
oparams
structure of the type "multimin_params" containing the optimization
parameters. The members are
double step_size
size of the first trial step
double tol
accuracy of the line minimization
unsigned maxiter
maximum number of iterations
double epsabs
accuracy of the minimization
unsigned method
method to use. Possible values are:
0: Fletcher-Reeves conjugate gradient
1: Polak-Ribiere conjugate gradient
2: Vector Broyden-Fletcher-Goldfarb-Shanno method
3: Steepest descent algorithm
unsigned verbosity
if greater then 0 print info on intermediate steps
*/
#include "common.h"
struct g_params {
size_t n;
unsigned *type;
double *xmin;
double *xmax;
void (*f) (const size_t,const double *,void *,double *);
void (* df) (const size_t,const double *, void *,double *);
void (* fdf) (const size_t,const double *, void *,double *,double *);
void *fparams;
};
static double g(const gsl_vector *y,void *gparams){
struct g_params *p= (struct g_params *) gparams;
size_t i;
double dtmp1;
double res=GSL_NAN;/* the function is forced to return a value */
/* dereference usefull stuff */
const size_t n = p->n;
const unsigned *type = p->type;
const double * xmin = p->xmin;
const double * xmax = p->xmax;
double *x = (double *) my_alloc(sizeof(double)*n);
/* compute values of x and of dx/dy */
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
x[i]= GET(y,i);
break;
case 1:/* [a,+inf) */
x[i]= xmin[i]+GET(y,i)*GET(y,i);
break;
case 2:/* (-inf,a] */
x[i]= xmax[i]-GET(y,i)*GET(y,i);
break;
case 3:/* [a,b] */
dtmp1 = sin( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 4:/* (a,+inf) */
dtmp1 = exp( GET(y,i) );
x[i]= xmin[i]+dtmp1;
break;
case 5:/* (-inf,a) */
dtmp1 = -exp( GET(y,i) );
x[i]= xmax[i]+dtmp1;
break;
case 6:/* (a,b) */
dtmp1 = tanh( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 7:/* (a,b) second approach */
dtmp1 = GET(y,i)/sqrt(1.+GET(y,i)*GET(y,i));
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 8:/* (a,+inf) second approach */
dtmp1 = sqrt(1.+GET(y,i)*GET(y,i));
x[i]= xmin[i] + .5*(dtmp1+GET(y,i));
break;
}
}
p->f(n,x,p->fparams,&res) ;
free (x);
return(res);
}
static void dg(const gsl_vector *y,void *gparams,gsl_vector *dg){
struct g_params *p= (struct g_params *) gparams;
size_t i;
double dtmp1,dtmp2;
/* dereference usefull stuff */
const size_t n = p->n;
const unsigned *type = p->type;
const double * xmin = p->xmin;
const double * xmax = p->xmax;
double *x = (double *) my_alloc(sizeof(double)*n);
double *dx = (double *) my_alloc(sizeof(double)*n);
double *df = (double *) my_alloc(sizeof(double)*n);
/* compute values of x and of dx/dy */
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
x[i]= GET(y,i);
dx[i]= 1;
break;
case 1:/* [a,+inf) */
x[i]= xmin[i]+GET(y,i)*GET(y,i);
dx[i]= 2.*GET(y,i);
break;
case 2:/* (-inf,a] */
x[i]= xmax[i]-GET(y,i)*GET(y,i);
dx[i]= -2.*GET(y,i);
break;
case 3:/* [a,b] */
dtmp1 = sin( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])*cos(GET(y,i));
break;
case 4:/* (a,+inf) */
dtmp1 = exp( GET(y,i) );
x[i]= xmin[i]+dtmp1;
dx[i]= dtmp1;
break;
case 5:/* (-inf,a) */
dtmp1 = -exp( GET(y,i) );
x[i]= xmax[i]+dtmp1;
dx[i]= dtmp1;
break;
case 6:/* (a,b) */
dtmp1 = tanh( GET(y,i) );
dtmp2 = cosh( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])/(dtmp2*dtmp2);
break;
case 7:/* (a,b) second approach */
dtmp1 = GET(y,i)/sqrt(1.+GET(y,i)*GET(y,i));
dtmp2 = (1.+GET(y,i)*GET(y,i))*sqrt(1.+GET(y,i)*GET(y,i)) ;
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])/dtmp2;
break;
case 8:/* (a,+inf) second approach */
dtmp1 = GET(y,i);
dtmp2 = sqrt(1.+dtmp1*dtmp1);
x[i]= xmin[i] + .5*(dtmp1+dtmp2);
dx[i]= .5*(dtmp1+dtmp2)/dtmp2;
break;
}
}
p->df(n,x,p->fparams,df);
/* fprintf(stderr,"#dg: x=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(x,i)); */
/* fprintf(stderr,") dx=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(dx,i)); */
/* fprintf(stderr,") df=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(dg,i)); */
for(i=0;i<n;i++){
SET(dg,i,df[i]*dx[i]);
}
/* fprintf(stderr,") dg=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(dg,i)); */
/* fprintf(stderr,")\n"); */
free (x);
free (dx);
free (df);
}
static void gdg(const gsl_vector *y,void *gparams,double *g,gsl_vector *dg){
struct g_params *p= (struct g_params *) gparams;
size_t i;
double dtmp1,dtmp2;
/* dereference usefull stuff */
const size_t n = p->n;
const unsigned *type = p->type;
const double * xmin = p->xmin;
const double * xmax = p->xmax;
double *x = (double *) my_alloc(sizeof(double)*n);
double *dx = (double *) my_alloc(sizeof(double)*n);
double *df = (double *) my_alloc(sizeof(double)*n);
/* compute values of x and of dx/dy */
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
x[i]= GET(y,i);
dx[i]= 1;
break;
case 1:/* [a,+inf) */
x[i]= xmin[i]+GET(y,i)*GET(y,i);
dx[i]= 2.*GET(y,i);
break;
case 2:/* (-inf,a] */
x[i]= xmax[i]-GET(y,i)*GET(y,i);
dx[i]= -2.*GET(y,i);
break;
case 3:/* [a,b] */
dtmp1 = sin( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])*cos(GET(y,i));
break;
case 4:/* (a,+inf) */
dtmp1 = exp( GET(y,i) );
x[i]= xmin[i]+dtmp1;
dx[i]= dtmp1;
break;
case 5:/* (-inf,a) */
dtmp1 = -exp( GET(y,i) );
x[i]= xmax[i]+dtmp1;
dx[i]= dtmp1;
break;
case 6:/* (a,b) */
dtmp1 = tanh( GET(y,i) );
dtmp2 = cosh( GET(y,i) );
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])/(dtmp2*dtmp2);
break;
case 7:/* (a,b) second approach */
dtmp1 = GET(y,i)/sqrt(1.+GET(y,i)*GET(y,i));
dtmp2 = (1.+GET(y,i)*GET(y,i))*sqrt(1.+GET(y,i)*GET(y,i)) ;
x[i]= .5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
dx[i]= .5*(xmax[i]-xmin[i])/dtmp2;
break;
case 8:/* (a,+inf) second approach */
dtmp1 = GET(y,i);
dtmp2 = sqrt(1.+dtmp1*dtmp1);
x[i]= xmin[i] + .5*(dtmp1+dtmp2);
dx[i]= .5*(dtmp1+dtmp2)/dtmp2;
break;
}
}
p->fdf(n,x,p->fparams,g,df);
/* fprintf(stderr,"#gdg: f=%f x=( ",g); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(x,i)); */
/* fprintf(stderr,") dx=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(dx,i)); */
/* fprintf(stderr,") df=( "); */
/* for(i=0;i<n;i++) */
/* fprintf(stderr,"%f ",GET(dg,i)); */
/* fprintf(stderr,")\n"); */
for(i=0;i<n;i++){
SET(dg,i,df[i]*dx[i]);
}
free (x);
free (dx);
free (df);
}
/*
n the dimension of the problem
x INPUT: initial guess OUTPUT: minimum point
f INPUT: ------------- OUTPUT: minimum value
type the types of the boundaries
xmin the minimum values
xmax the maximum values
f the structure of a function
fparams the parameters of the provided function
oparams parameters of the optimization
*/
void
multimin(size_t n,double *x,double *fun,
unsigned *type,double *xmin,double *xmax,
void (*f) (const size_t,const double *,void *,double *),
void (* df) (const size_t,const double *, void *,double *),
void (* fdf) (const size_t,const double *, void *,double *,double *),
void *fparams,
const struct multimin_params oparams)
{
unsigned iter=0;
int status;
size_t i;
double dtmp1;
const gsl_multimin_fdfminimizer_type *T;
gsl_multimin_fdfminimizer *s;
gsl_vector * y = gsl_vector_alloc (n);
gsl_multimin_function_fdf GdG;
struct g_params gparams;
/* set the algorithm */
switch(oparams.method){
case 0:/* Fletcher-Reeves conjugate gradient */
T = gsl_multimin_fdfminimizer_conjugate_fr;
break;
case 1:/* Polak-Ribiere conjugate gradient */
T = gsl_multimin_fdfminimizer_conjugate_pr;
break;
case 2:/* Vector Broyden-Fletcher-Goldfarb-Shanno method */
T = gsl_multimin_fdfminimizer_vector_bfgs;
break;
case 3:/* Steepest descent algorithm */
T =gsl_multimin_fdfminimizer_steepest_descent;
break;
default:
fprintf(stderr,"Optimization method not recognized. Try -h\n");
exit(EXIT_FAILURE);
};
s = gsl_multimin_fdfminimizer_alloc(T, n);
/* --- OUPUT ---------------------------------- */
if(oparams.verbosity>0){
fprintf(stderr,"#--- MULTIMIN START [method %s]\n",T->name);
fprintf(stderr,"# [step_size=%.3e tol=%.3e maxiter=%u epsabs=%.3e]\n",
oparams.step_size,oparams.tol,oparams.maxiter,oparams.epsabs);
}
/* -------------------------------------------- */
/* set the parameters of the new function */
gparams.n = n;
gparams.type = type;
gparams.xmin = xmin;
gparams.xmax = xmax;
gparams.f = f;
gparams.df = df;
gparams.fdf = fdf;
gparams.fparams = fparams;
/* set the function to solve */
GdG.f=g;
GdG.df=dg;
GdG.fdf=gdg;
GdG.n=n;
GdG.params=(void *) &gparams;
/* compute values of y for initial condition */
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
SET(y,i,x[i]);
break;
case 1:/* [a,+inf) */
SET(y,i,sqrt( x[i]-xmin[i] ));
break;
case 2:/* (-inf,a] */
SET(y,i,sqrt( xmax[i]-x[i] ));
break;
case 3:/* [a,b] */
dtmp1 = (xmax[i]>xmin[i]?
(2.*x[i]-xmax[i]-xmin[i])/(xmax[i]-xmin[i]) : 0);
/* dtmp1 = (2.*x[i]-xmax[i]-xmin[i])/(xmax[i]-xmin[i]); */
SET(y,i,asin( dtmp1 ));
break;
case 4:/* (a,+inf) */
SET(y,i,log( x[i]-xmin[i] ));
break;
case 5:/* (-inf,a) */
SET(y,i,log( xmax[i]-x[i] ));
break;
case 6:/* (a,b) */
dtmp1 = (2.*x[i]-xmax[i]-xmin[i])/(xmax[i]-xmin[i]);
SET(y,i,gsl_atanh ( dtmp1 ));
break;
case 7:/* (a,b) second approach */
dtmp1 = (2.*x[i]-xmax[i]-xmin[i])/(xmax[i]-xmin[i]);
SET(y,i, dtmp1/sqrt(1-dtmp1*dtmp1));
break;
case 8:/* (a,+inf) second approach */
dtmp1 = x[i]-xmin[i];
SET(y,i, dtmp1-1./(4.*dtmp1));
break;
}
}
/* --- OUPUT ---------------------------------- */
if(oparams.verbosity>1){
fprintf(stderr,"# - initial values and boundaries\n");
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
fprintf(stderr,"# x[%d]=%e (-inf,+inf) -> %e\n",(int) i,x[i],GET(y,i));
break;
case 1:/* [a,+inf) */
fprintf(stderr,"# x[%d]=%e [%f,+inf) -> %e\n",(int) i,x[i],xmin[i],GET(y,i));
break;
case 2:/* (-inf,a] */
fprintf(stderr,"# x[%d]=%e (-inf,%f] -> %e\n",(int) i,x[i],xmax[i],GET(y,i));
break;
case 3:/* [a,b] */
fprintf(stderr,"# x[%d]=%e [%f,%f] -> %e\n",(int) i,x[i],xmin[i],xmax[i],GET(y,i));
break;
case 4:/* (a,+inf) */
fprintf(stderr,"# x[%d]=%e (%f,+inf) -> %e\n",(int) i,x[i],xmin[i],GET(y,i));
break;
case 5:/* (-inf,a) */
fprintf(stderr,"# x[%d]=%e (-inf,%f) -> %e\n",(int) i,x[i],xmax[i],GET(y,i));
break;
case 6:/* (a,b) */
case 7:
fprintf(stderr,"# x[%d]=%e (%f,%f) -> %e\n",(int) i,x[i],xmin[i],xmax[i],GET(y,i));
break;
case 8:/* [a,+inf) */
fprintf(stderr,"# x[%d]=%e (%f,+inf) -> %e\n",(int) i,x[i],xmin[i],GET(y,i));
break;
}
}
}
/* -------------------------------------------- */
/* initialize minimizer */
status=gsl_multimin_fdfminimizer_set(s,&GdG,y,oparams.step_size,oparams.tol);
if(status)
{
fprintf(stderr,"#ERROR: %s\n",gsl_strerror (status));
exit(EXIT_FAILURE);
}
if(oparams.verbosity>2)
fprintf(stderr,"# - start minimization \n");
do
{
iter++;
status = gsl_multimin_fdfminimizer_iterate (s);
/* --- OUPUT ---------------------------------- */
if(oparams.verbosity>2){
fprintf(stderr,"# g=%f y=( ",s->f);
for(i=0;i<n;i++)
fprintf(stderr,"%f ",GET(s->x,i));
fprintf(stderr,") dg=( ");
for(i=0;i<n;i++)
fprintf(stderr,"%f ",GET(s->gradient,i));
fprintf(stderr,")\n");
}
/* -------------------------------------------- */
if(status){
fprintf(stderr,"#WARNING: %s\n", gsl_strerror (status));
break;
}
status = gsl_multimin_test_gradient (s->gradient,oparams.epsabs);
if(iter>=oparams.maxiter) break;
}
while (status == GSL_CONTINUE);
/* compute values of x */
for(i=0;i<n;i++){
switch(type[i]){
case 0:/* (-inf,+inf) */
x[i]=GET(s->x,i);
break;
case 1:/* [a,+inf) */
x[i]=xmin[i]+GET(s->x,i)*GET(s->x,i);
break;
case 2:/* (-inf,a] */
x[i]=xmax[i]-GET(s->x,i)*GET(s->x,i);
break;
case 3:/* [a,b] */
dtmp1 = sin( GET(s->x,i) );
x[i]=.5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 4:/* (a,+inf) */
dtmp1 = exp( GET(s->x,i) );
x[i]=xmin[i]+dtmp1;
break;
case 5:/* (-inf,a) */
dtmp1 = -exp( GET(s->x,i) );
x[i]=xmax[i]+dtmp1;
break;
case 6:/* (a,b) */
dtmp1 = tanh( GET(s->x,i) );
x[i]=.5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 7:/* (a,b) second approach */
dtmp1 = GET(s->x,i) ;
dtmp1 = dtmp1/sqrt(1.+dtmp1*dtmp1);
x[i]=.5*(xmin[i]*(1-dtmp1) +xmax[i]*(1+dtmp1));
break;
case 8:/* (a,+inf) second approach */
dtmp1 = sqrt(1.+GET(s->x,i)*GET(s->x,i));
x[i]= xmin[i] + .5*(dtmp1+GET(s->x,i));
break;
}
}
*fun=s->f;
/* --- OUPUT ---------------------------------- */
if(oparams.verbosity>0){
fprintf(stderr,"#--- MULTIMIN END --- ");
fprintf(stderr,"iterations %u\n",iter);
}
/* -------------------------------------------- */
gsl_multimin_fdfminimizer_free (s);
gsl_vector_free (y);
}
/*insert GNU extensions*/
#define _GNU_SOURCE
/*in particular, use of NAN extension*/
/* used by <errno.h> */
extern int errno;
/* GSL ----------- */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>
/* --------------- */
#define GET(x,i) gsl_vector_get(x,i)
#define SET(x,i,y) gsl_vector_set(x,i,y)
struct multimin_params {
double step_size;
double tol;
unsigned maxiter;
double epsabs;
unsigned method;
unsigned verbosity;
};
void
multimin(size_t,double *,double *,
unsigned *,double *,double *,
void (*) (const size_t,const double *,void *,double *),
void (*) (const size_t,const double *, void *,double *),
void (*) (const size_t,const double *, void *,double *,double *),
void *,
const struct multimin_params);
Attachment:
pgp00000.pgp
Description: PGP signature
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |