freefem++-doc 3.19.1-1 / usr / share / doc / freefem++ / examples / examples++-load / ffnewuoa.edp

/*  README NEWUOA:
     This is the Fortran version of NEWUOA. Its purpose is to seek
the least value of a function F of several variables, when derivatives
are not available, where F is specified by the user through a subroutine
called CALFUN. The algorithm is intended to change the variables to values
that are close to a local minimum of F. The user, however, should assume
responsibility for finding out if the calculations are satisfactory, by
considering carefully the values of F that occur. The method is described
in the report "The NEWUOA software for unconstrained optimization without
derivatives", which is available on the web at www.damtp.cam.ac.uk, where
you have to click on Numerical Analysis and then on Reports, the number
of the report being NA2004/08. Let N be the number of variables. The main
new feature of the method is that quadratic models are updated using only
about NPT=2N+1 interpolation conditions, the remaining freedom being taken
up by minimizing the Frobenius norm of the change to the second derivative
matrix of the model.

     The new software was developed from UOBYQA, which also forms quadratic
models from interpolation conditions. That method requires NPT=(N+1)(N+2)/2
conditions, however, because they have to define all the parameters of the
model. The least Frobenius norm updating procedure with NPT=2N+1 is usually
much more efficient when N is large, because the work of each iteration is
much less than before, and in some experiments the number of calculations
of the objective function seems to be only of magnitude N.

     The attachments in sequence are a suitable Makefile, followed by a main
program and a CALFUN routine for the Chebyquad problems, in order to provide
an example for testing. Then NEWUOA and its five auxiliary routines, namely
NEWUOB, BIGDEN, BIGLAG, TRSAPP and UPDATE, are given. Finally, the computed
output that the author obtained for the Chebyquad problems is listed.

     The way of calling NEWUOA should be clear from the Chebyquad example
and from the comments of that subroutine. It is hoped that the software will
be helpful to much future research and to many applications. There are no
restrictions on or charges for its use. If you wish to refer to it, please
cite the DAMTP report that is mentioned above, which has been submitted for
publication in the proceedings of the 40th Workshop on Large Scale Nonlinear
Optimization (Erice, Italy, 2004).

December 16th, 2004                    M.J.D. Powell (mjdp@cam.ac.uk)
*/

load "ffnewuoa"
/*
      SUBROUTINE CALFUN (N,X,F,IWF)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(*),Y(10,10),IWF(*)
      DO 10 J=1,N
      Y(1,J)=1.0D0
   10 Y(2,J)=2.0D0*X(J)-1.0D0
      DO 20 I=2,N
      DO 20 J=1,N
   20 Y(I+1,J)=2.0D0*Y(2,J)*Y(I,J)-Y(I-1,J)
      F=0.0D0
      NP=N+1
      IW=1
      DO 40 I=1,NP
      SUM=0.0D0
      DO 30 J=1,N
   30 SUM=SUM+Y(I,J)
      SUM=SUM/DFLOAT(N)
      IF (IW .GT. 0) SUM=SUM+1.0D0/DFLOAT(I*I-2*I)
      IW=-IW
   40 F=F+SUM*SUM
      RETURN
      END 

*/
int N=2;
real[int,int] Y(N+1,N);
real[int] X(N);
func real J(real[int] & X)
{
  Y(0,:)=1.;
  Y(1,:)=2.*X;
  for(int i=2;i<=N;++i)
    for(int j=0;j<N;++j)
      Y(i,j) = Y(1,j) * Y(i-1,j)- Y(i-2,j);

  real f=0;
  int np=N+1,iw=1;
  for(int i=0;i<=N;i++)
    {
      real s=Y(i,:).sum/N;
      int i1=i+1;
      if(iw>0) s+= 1./(i1*i1-2.*i1);
      iw=-iw;
      f += s*s;
    }
  return 2.14+f;
}
real delta=1;
mesh Th=square(30,30,[(x-0.5)*delta,(y-0.5)*delta]);
fespace Vh(Th,P1);
Vh u;
for(int i=0;i<Th.nv;++i)
  {X[0]=Th(i).x;X[1]=Th(i).y;
    u[][i]=J(X);
  }
plot(u,wait=1);

for(int i=0;i<N;++i)
  X[i]=(i+1)/(N+1.);
real mincost=newuoa(J,X,rhobeg=2*X[0],rhoend=1e-6,npt=2*N+1);
cout << " min " << mincost << "  at :" << X << endl;