Nelder-Mead method [algorithm] by 바죠

http://en.wikipedia.org/wiki/Nelder-Mead_method
http://www.cs.ucsb.edu/~kyleklein/publications/neldermead.pdf

다차원 독립변수들에 의한 목적 함수 최적화에 사용되는 방법이다. 아메바 방법이라고도 한다. 명백히 stochastic 접근 방법이다.  도함수(gradient) 를 활용하지 않는다. 따라서, 매우 일반적인 목적 함수 최적화에 응용될 수 있다. 문제에 따라서는 목적 함수의 도함수가 잘 정의되지 않는 경우가 많다. 이러할 경우 목적 함수 최적화가 어려워진다. 왜냐하면 목적 함수의 변화에 대한 정보를 최적화 과정에서 사용하지 않기 때문이다. 효율성은 떨어질 수 있다. 하지만, 실제 응용문제 풀이에서는 그러한 경우의 문제들이 매우 많이 산재한다.

Nelder-Mead 방법은 simulated annealing 처럼 계속해서 새로운 해를 찾아나서는 시도를 반복한다. 함수 피팅같은 간단한 계산에서 유용하게 사용될 수 있다. 몇 개의 독립변수들을 최적화 시킬 때 매우 유용하다고 알려져 있다. 국소 최적화 방법으로 사용될 수 있다. 특히 도함수가 알려지지 않을 경우에 유용하다. 도함수가 알려진 경우에는 BFGS 알고리듬이 유용하다. 도함수 정보를 반드시 이용해야만 효율적인 계산이 가능하다. 물론 도함수가 알려진 경우에 한해서 그렇다. 문제에 따라서는 도함수가 알려지지 않은 경우도 많다. 

아래의 그림에서처럼, 2차원의 경우 3개의 점을 이용한다. 서로 다른 2차원 상의 3점에서 각각 함수값을 계산한다. 그렇게 하면 2차원에서 3각형을 만들 수 있다. 이 때, 우리는 어느 방향으로 가야만 더 낮은 함수값을 주는 영역인가를 어느 정도 알 수 있다. 얼마만큼 크게 움직여야 적당한지는 새로운 시도를 통해서 대략적으로 알 수 있다. 너무 많이 가면 뒤로 되돌아와야만 할 것이다. 2차원에서 3점은 결국 아메바처럼 움직인다.  2차원 면을 기어다니는 꼴이다. 국소점을 찾는 것은 결국 아메바의 크기가 작아진다는 것을 의미한다. 결국에는 조금씩 움직여야만 할 것이다. 처음에는 다소 크게 움직이는 것이 가능할 것이다.  초반에 크게 움직일 때, 우리는 국소점을 놓칠 확률이 있다. 당연한 것이다. 하지만, 아메바의 크기가 궁극적으로는 작아진다. 이렇게 될 경우, 국소점은 찾아지게 된다. 3점이 만드는 삼각형은 일단 크기가 변한다. 회전할 수도 있다. 꼭지점이 반대편 모서리 쪽으로 진출하여 새로운 형태의 삼각형을 만들 수 있다. 삼각형 내부에 국소점이 놓이게하는 것이 알고리듬의 핵심 추구사항이다.

Nelder-Mead 방법을 실제로 사용할 때는 매우 간단하게 적용하는 방법을 택한다. 즉, 사용자는 단순하게 fn(x)라는 함수를 설계함으로써 기존의 Nelder-Mead 방법을 그대로 이용할 수 있다. 사용자 함수 fn(x)는 입력과 출력이 각각 x, fn으로 구별지어진다. 따라서 사용자는 fn(x)를 설계할 때, module 을 사용하여 사용자가 제공해야만 하는 도메인 정보를 fn 함수에 제공할 수 있다. fn함수는 x 이외에는 정보를 제공받을 수 없다. 출력은 fn이다. 입력은 x 벡터이다. 물론, x가 벡터량이기 때문에 벡터의 길이를 변화시키면서 필요한 정보를 제공할 수도 있다.  적절한 방법은 아닐 것이다.  그렇게 하지 않고 단순히 사용자 module 을 정의 하고, fn 함수에서 USE 를 사용함으로써, 변수들 또는 부함수들을 불러서 도메인에서 필요한 계산을 수행할 수 있다. 이렇게 하면, Nelder-Mead 방법을 전혀 손대지 않고 그대로 사용할 수 있다.  아래의 예제들을 보자. Nelder-Mead 방법이 간단하기 때문에 이것을 변형시켜서 메인 프로그램을 만들려고 시도하려는 유혹을 느낀다. 하지만, 그렇게 할 필요가 없다. 계산과 프로그램 관리에 전혀 도움이 되지 않는다.

기억하자. 새로운 수학이 나오지 않는 한, 대부분의 수치 프로그램은 결국에는 궁극에는 잘 알려진 기본 수학 라이버러리들의 사용으로 귀결되게 된다. 사용자는 결국, 가장 적절한 라이버러리를 찾고 문제에 맞게 사용함으로써 최고의 계산 효율을 얻어낼 수 있게 되는 것이다.  사용자는 가장 확실한 계산 방법들을 연속적으로 찾는 것이다. 그것이 프로그래밍이다. 그 이상은 알고리듬을 만들어 내는 것이다. 그 경우에도 마찬가지이다. 기존의 알고리듬과 비교는 불가피하다.


http://people.sc.fsu.edu/~jburkardt/f_src/asa047/asa047_prb.f90

http://people.sc.fsu.edu/~jburkardt/f_src/asa047/asa047.f90






-------------------------------------------------------------------------------------------------------------------

!      Written by In-Ho Lee, KRISS, October 14, 2017.
       subroutine nmbfgsbody(start0,rslt)
       implicit none
       real*8 start0(nn),rslt
       integer icount,ifault,numres
       real*8 ynewlo
       real*8, external :: fn

       if(lgradient)then
       call lbfgsmini(rslt,start0,nn)
                    else
       start=start0
       call nelmin(fn,nn,start,xmin,ynewlo,reqmin,step,konvge,kcount,icount,numres,ifault)
       start0=xmin ; rslt=ynewlo
                    endif
       end subroutine nmbfgsbody
!234567890
!      Written by In-Ho Lee, KRISS, March 3, 2003.
       subroutine lbfgsmini(object,xxxx,nn1)
       implicit none
       integer nn1
       real*8 object,xxxx(nn1)
       integer m,msave
       real*8, allocatable :: diag(:),w(:)
       real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax
       integer iprint(2),iflag,icall,mp,lp,nwork
       logical diagco
       real*8, external :: fn
!
!     The driver for LBFGS must always declare LB2 as EXTERNAL
       external lb2
       common /lb3/mp,lp,gtol,stpmin,stpmax
!
       m=5 ; msave=7 ; nwork=nn1*(2*msave+1)+2*msave
       allocate(diag(nn1),w(nwork))
       iprint(1)= 1 ; iprint(2)= 0
!
!     We do not wish to provide the diagonal matrices Hk0, and
!     therefore set DIAGCO to FALSE.
       diagco= .false. ; eps= 1.0d-4 ; xtol= 1.0d-16 ; icall=0 ; iflag=0
!
 20   continue
      object=fn(xxxx)
!
      call lbfgsac(nn1,m,xxxx,object,gggg,diagco,diag,iprint,eps,xtol,w,iflag)
!
      if(iflag <= 0) go to 50
      icall=icall + 1
!     We allow at most 2000 evaluations of Function and Gradient
      if(icall >  kcount) go to 50
      go to 20
  50  continue
       deallocate(diag,w)
       end subroutine lbfgsmini


-------------------------------------------------------------------------------------------------------------------




subroutine nelmin ( fn, n, start, xmin, ynewlo, reqmin, step, konvge, kcount, &
  icount, numres, ifault )

!*****************************************************************************80
!
!! NELMIN minimizes a function using the Nelder-Mead algorithm.
!
!  Discussion:
!
!    This routine seeks the minimum value of a user-specified function.
!
!    Simplex function minimisation procedure due to Nelder and Mead (1965),
!    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
!    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
!    25, 97) and Hill(1978, 27, 380-2)
!
!    The function to be minimized must be defined by a function of
!    the form
!
!      function fn ( x)
!      real ( kind = 8 ) fn
!      real ( kind = 8 ) x(*)
!
!    and the name of this subroutine must be declared EXTERNAL in the
!    calling routine and passed as the argument FN.
!
!    This routine does not include a termination test using the
!    fitting of a quadratic surface.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    27 February 2008
!
!  Author:
!
!    Original FORTRAN77 version by R ONeill.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    John Nelder, Roger Mead,
!    A simplex method for function minimization,
!    Computer Journal,
!    Volume 7, 1965, pages 308-313.
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, external FN, the name of the function which evaluates
!    the function to be minimized.
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!    0 < N is required.
!
!    Input/output, real ( kind = 8 ) START(N).  On input, a starting point
!    for the iteration.  On output, this data may have been overwritten.
!
!    Output, real ( kind = 8 ) XMIN(N), the coordinates of the point which
!    is estimated to minimize the function.
!
!    Output, real ( kind = 8 ) YNEWLO, the minimum value of the function.
!
!    Input, real ( kind = 8 ) REQMIN, the terminating limit for the variance
!    of the function values.  0 < REQMIN is required.
!
!    Input, real ( kind = 8 ) STEP(N), determines the size and shape of the
!    initial simplex.  The relative magnitudes of its elements should reflect
!    the units of the variables.
!
!    Input, integer ( kind = 4 ) KONVGE, the convergence check is carried out
!    every KONVGE iterations. 0 < KONVGE is required.
!
!    Input, integer ( kind = 4 ) KCOUNT, the maximum number of function
!    evaluations.
!
!    Output, integer ( kind = 4 ) ICOUNT, the number of function evaluations
!    used.
!
!    Output, integer ( kind = 4 ) NUMRES, the number of restarts.
!
!    Output, integer ( kind = 4 ) IFAULT, error indicator.
!    0, no errors detected.
!    1, REQMIN, N, or KONVGE has an illegal value.
!    2, iteration terminated because KCOUNT was exceeded without convergence.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ), parameter :: ccoeff = 0.5D+00
  real ( kind = 8 ) del
  real ( kind = 8 ), parameter :: ecoeff = 2.0D+00
  real ( kind = 8 ), parameter :: eps = 0.001D+00
  real ( kind = 8 ), external :: fn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jcount
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) l
  integer ( kind = 4 ) numres
  real ( kind = 8 ) p(n,n+1)
  real ( kind = 8 ) p2star(n)
  real ( kind = 8 ) pbar(n)
  real ( kind = 8 ) pstar(n)
  real ( kind = 8 ), parameter :: rcoeff = 1.0D+00
  real ( kind = 8 ) reqmin
  real ( kind = 8 ) rq
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) x
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) y(n+1)
  real ( kind = 8 ) y2star
  real ( kind = 8 ) ylo
  real ( kind = 8 ) ynewlo
  real ( kind = 8 ) ystar
  real ( kind = 8 ) z
!
!  Check the input parameters.
!
  if ( reqmin <= 0.0D+00 ) then
    ifault = 1
    return
  end if

  if ( n < 1 ) then
    ifault = 1
    return
  end if

  if ( konvge < 1 ) then
    ifault = 1
    return
  end if
!
!  Initialization.
!
  icount = 0
  numres = 0
  jcount = konvge
  del = 1.0D+00
  rq = reqmin * real ( n, kind = 8 )
!
!  Initial or restarted loop.
!
  do

    p(1:n,n+1) = start(1:n)
    y(n+1) = fn ( start )
    icount = icount + 1
!
!  Define the initial simplex.
!
    do j = 1, n
      x = start(j)
      start(j) = start(j) + step(j) * del
      p(1:n,j) = start(1:n)
      y(j) = fn ( start )
      icount = icount + 1
      start(j) = x
    end do
!
!  Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
!  the vertex of the simplex to be replaced.
!
    ilo = minloc ( y(1:n+1), 1 )
    ylo = y(ilo)
!
!  Inner loop.
!
    do while ( icount < kcount )
!
!  YNEWLO is, of course, the HIGHEST value???
!
      ihi = maxloc ( y(1:n+1), 1 )
      ynewlo = y(ihi)
!
!  Calculate PBAR, the centroid of the simplex vertices
!  excepting the vertex with Y value YNEWLO.
!
      do i = 1, n
        pbar(i) = ( sum ( p(i,1:n+1) ) - p(i,ihi) ) / real ( n, kind = 8 )
      end do
!
!  Reflection through the centroid.
!
      pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
      ystar = fn ( pstar )
      icount = icount + 1
!
!  Successful reflection, so extension.
!
      if ( ystar < ylo ) then

        p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
        y2star = fn ( p2star )
        icount = icount + 1
!
!  Retain extension or contraction.
!
        if ( ystar < y2star ) then
          p(1:n,ihi) = pstar(1:n)
          y(ihi) = ystar
        else
          p(1:n,ihi) = p2star(1:n)
          y(ihi) = y2star
        end if
!
!  No extension.
!
      else

        l = 0
        do i = 1, n + 1
          if ( ystar < y(i) ) then
            l = l + 1
          end if
        end do

        if ( 1 < l ) then

          p(1:n,ihi) = pstar(1:n)
          y(ihi) = ystar
!
!  Contraction on the Y(IHI) side of the centroid.
!
        else if ( l == 0 ) then

          p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
          y2star = fn ( p2star )
          icount = icount + 1
!
!  Contract the whole simplex.
!
          if ( y(ihi) < y2star ) then

            do j = 1, n + 1
              p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5D+00
              xmin(1:n) = p(1:n,j)
              y(j) = fn ( xmin )
              icount = icount + 1
            end do

            ilo = minloc ( y(1:n+1), 1 )
            ylo = y(ilo)

            cycle
!
!  Retain contraction.
!
          else
            p(1:n,ihi) = p2star(1:n)
            y(ihi) = y2star
          end if
!
!  Contraction on the reflection side of the centroid.
!
        else if ( l == 1 ) then

          p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
          y2star = fn ( p2star )
          icount = icount + 1
!
!  Retain reflection?
!
          if ( y2star <= ystar ) then
            p(1:n,ihi) = p2star(1:n)
            y(ihi) = y2star
          else
            p(1:n,ihi) = pstar(1:n)
            y(ihi) = ystar
          end if

        end if

      end if
!
!  Check if YLO improved.
!
      if ( y(ihi) < ylo ) then
        ylo = y(ihi)
        ilo = ihi
      end if

      jcount = jcount - 1

      if ( 0 < jcount ) then
        cycle
      end if
!
!  Check to see if minimum reached.
!
      if ( icount <= kcount ) then

        jcount = konvge

        x = sum ( y(1:n+1) ) / real ( n + 1, kind = 8 )
        z = sum ( ( y(1:n+1) - x )**2 )

        if ( z <= rq ) then
          exit
        end if

      end if

    end do
!
!  Factorial tests to check that YNEWLO is a local minimum.
!
    xmin(1:n) = p(1:n,ilo)
    ynewlo = y(ilo)

    if ( kcount < icount ) then
      ifault = 2
      exit
    end if

    ifault = 0

    do i = 1, n
      del = step(i) * eps
      xmin(i) = xmin(i) + del
      z = fn ( xmin )
      icount = icount + 1
      if ( z < ynewlo ) then
        ifault = 2
        exit
      end if
      xmin(i) = xmin(i) - del - del
      z = fn ( xmin )
      icount = icount + 1
      if ( z < ynewlo ) then
        ifault = 2
        exit
      end if
      xmin(i) = xmin(i) + del
    end do

    if ( ifault == 0 ) then
      exit
    end if
!
!  Restart the procedure.
!
    start(1:n) = xmin(1:n)
    del = eps
    numres = numres + 1

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end


program main

!*****************************************************************************80
!
!! MAIN is the main program for ASA047_PRB.
!
!  Discussion:
!
!    ASA047_PRB calls the ASA047 routines.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ASA047_PRB:'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test the ASA047 library.'

  call test01 ( )
  call test02 ( )
  call test03 ( )
  call test04 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ASA047_PRB:'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test01 ( )

!*****************************************************************************80
!
!! TEST01 demonstrates the use of NELMIN on ROSENBROCK.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 2

  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) numres
  real ( kind = 8 ) reqmin
  real ( kind = 8 ), external :: rosenbrock
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Apply NELMIN to the ROSENBROCK function.'

  start(1:n) = (/ -1.2D+00, 1.0D+00 /)

  reqmin = 1.0D-08

  step(1:n) = (/ 1.0D+00, 1.0D+00 /)

  konvge = 10
  kcount = 500

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Starting point X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) start(i)
  end do

  ynewlo = rosenbrock ( start )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X) = ', ynewlo

  call nelmin ( rosenbrock, n, start, xmin, ynewlo, reqmin, step, &
    konvge, kcount, icount, numres, ifault )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Return code IFAULT = ', ifault
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Estimate of minimizing value X*:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) xmin(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X*) = ', ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations = ', icount
  write ( *, '(a,i8)' ) '  Number of restarts =   ', numres

  return
end
function rosenbrock ( x )

!*****************************************************************************80
!
!! ROSENBROCK evaluates the Rosenbrock parabolic value function.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X(2), the argument.
!
!    Output, real ( kind = 8 ) ROSENBROCK, the value of the function.
!
  implicit none

  real ( kind = 8 ) fx
  real ( kind = 8 ) fx1
  real ( kind = 8 ) fx2
  real ( kind = 8 ) rosenbrock
  real ( kind = 8 ) x(3)

  fx1 = x(2) - x(1) * x(1)
  fx2 = 1.0D+00 - x(1)

  fx = 100.0D+00 * fx1 * fx1 &
     +             fx2 * fx2

  rosenbrock = fx

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! TEST02 demonstrates the use of NELMIN on POWELL.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 4

  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) numres
  real ( kind = 8 ), external :: powell
  real ( kind = 8 ) reqmin
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Apply NELMIN to POWELL quartic function.'

  start(1:n) = (/  3.0D+00, - 1.0D+00,   0.0D+00,   1.0D+00 /)

  reqmin = 1.0D-08

  step(1:n) = (/ 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)

  konvge = 10
  kcount = 500

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Starting point X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) start(i)
  end do

  ynewlo = powell ( start )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X) = ', ynewlo

  call nelmin ( powell, n, start, xmin, ynewlo, reqmin, step, &
    konvge, kcount, icount, numres, ifault )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Return code IFAULT = ', ifault
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Estimate of minimizing value X*:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) xmin(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X*) = ', ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations = ', icount
  write ( *, '(a,i8)' ) '  Number of restarts =   ', numres

  return
end
function powell ( x )

!*****************************************************************************80
!
!! POWELL evaluates the Powell quartic function.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X(4), the argument.
!
!    Output, real ( kind = 8 ) POWELL, the value of the function.
!
  implicit none

  real ( kind = 8 ) fx
  real ( kind = 8 ) fx1
  real ( kind = 8 ) fx2
  real ( kind = 8 ) fx3
  real ( kind = 8 ) fx4
  real ( kind = 8 ) powell
  real ( kind = 8 ) x(4)

  fx1 = x(1) + 10.0D+00 * x(2)
  fx2 = x(3) - x(4)
  fx3 = x(2) - 2.0D+00 * x(3)
  fx4 = x(1) - x(4)

  fx =            fx1 * fx1 &
     +  5.0D+00 * fx2 * fx2 &
     +            fx3 * fx3 * fx3 * fx3 &
     + 10.0D+00 * fx4 * fx4 * fx4 * fx4

  powell = fx

  return
end
subroutine test03 ( )

!*****************************************************************************80
!
!! TEST03 demonstrates the use of NELMIN on HELICAL.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 3

  real ( kind = 8 ), external :: helical
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) numres
  real ( kind = 8 ) reqmin
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Apply NELMIN to the HELICAL function.'

  start(1:n) = (/ - 1.0D+00,   0.0D+00,   0.0D+00 /)

  reqmin = 1.0D-08

  step(1:n) = (/ 1.0D+00, 1.0D+00, 1.0D+00 /)

  konvge = 10
  kcount = 500

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Starting point X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) start(i)
  end do

  ynewlo = helical ( start )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X) = ', ynewlo

  call nelmin ( helical, n, start, xmin, ynewlo, reqmin, step, &
    konvge, kcount, icount, numres, ifault )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Return code IFAULT = ', ifault
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Estimate of minimizing value X*:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) xmin(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X*) = ', ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations = ', icount
  write ( *, '(a,i8)' ) '  Number of restarts =   ', numres

  return
end
function helical ( x )

!*****************************************************************************80
!
!! HELICAL evaluates the Fletcher-Powell helical valley function.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X(3), the argument.
!
!    Output, real ( kind = 8 ) HELICAL, the value of the function.
!
  implicit none

  real ( kind = 8 ) fx
  real ( kind = 8 ) fx1
  real ( kind = 8 ) fx2
  real ( kind = 8 ) fx3
  real ( kind = 8 ) helical
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) theta
  real ( kind = 8 ) x(3)

  if ( 0.0D+00 < x(1) ) then
    theta = atan2 ( x(2), x(1) ) / 2.0D+00 / pi
  else if ( x(1) < 0.0D+00 ) then
    theta = 0.5D+00 + atan2 ( x(2), x(1) ) / 2.0D+00 / pi
  else if ( x(1) == 0.0D+00 ) then
    theta = 0.25D+00
  end if

  fx1 = x(3) - 10.0D+00 * theta
  fx2 = sqrt ( x(1) * x(1) + x(2) * x(2) )
  fx3 = x(3)

  fx = 100.0D+00 * fx1 * fx1 &
     +             fx2 * fx2 &
     +             fx3 * fx3

  helical = fx

  return
end
subroutine test04 ( )

!*****************************************************************************80
!
!! TEST04 demonstrates the use of NELMIN on QUARTIC.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: n = 10

  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) numres
  real ( kind = 8 ), external :: quartic
  real ( kind = 8 ) reqmin
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  Apply NELMIN to the QUARTIC function.'

  start(1:n) = 1.0D+00

  reqmin = 1.0D-08

  step(1:n) = 1.0D+00

  konvge = 10
  kcount = 2000

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Starting point X:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) start(i)
  end do

  ynewlo = quartic ( start )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X) = ', ynewlo

  call nelmin ( quartic, n, start, xmin, ynewlo, reqmin, step, &
    konvge, kcount, icount, numres, ifault )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Return code IFAULT = ', ifault
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Estimate of minimizing value X*:'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,g14.6)' ) xmin(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  F(X*) = ', ynewlo

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations = ', icount
  write ( *, '(a,i8)' ) '  Number of restarts =   ', numres

  return
end
function quartic ( x )

!*****************************************************************************80
!
!! QUARTIC evaluates a function defined by a sum of fourth powers.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    19 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X(10), the argument.
!
!    Output, real ( kind = 8 ) QUARTIC, the value of the function.
!
  implicit none

  real ( kind = 8 ) quartic
  real ( kind = 8 ) x(10)

  quartic = sum ( x(1:10)**4 )

  return
end



http://incredible.egloos.com/4856511
http://incredible.egloos.com/4784590


!23456789
       implicit none
       real ( kind = 8 ), external :: fn
       integer ( kind = 4 ) n,konvge,kcount,icount,numres,ifault
       real ( kind = 8 ) ynewlo,reqmin
       real ( kind = 8 ), allocatable :: start(:),step(:),xmin(:)
       integer i
       real*8 tmp

       n=3
       allocate(start(n),step(n),xmin(n))
       do i=1,n
       call random_number(tmp)
       start(i)=tmp-0.5d0
       step(i)=(tmp-0.5d0)*0.1d0
       enddo
       reqmin=1.d-3
       konvge=10
       kcount=100000
       write(6,*) kcount, ' the maximum number of function evaluations'
       write(6,*) konvge, ' the convergence check is carried out every', konvge,' iterations'
       write(6,*) fn(start)
       write(6,*) start
       call nelmin(fn,n,start,xmin,ynewlo,reqmin,step,konvge,kcount,icount,numres,ifault)
       write(6,*) numres,' the number of restarts'
       write(6,*) icount,' the number of function evaluations'
       write(6,*) ifault,' ifault'
       if(ifault == 0) write(6,*) ' no errors detected'
       if(ifault == 1) write(6,*) ' REQMIN, N, or KONVGE has an illegal value'
       if(ifault == 2) write(6,*) ' iteration terminated because KCOUNT was exceeded without convergence'
       write(6,*) fn(xmin)
       write(6,*) xmin
       deallocate(start,step,xmin)
       stop
       end
!23456789
       function fn (x)
       real ( kind = 8 ) fn
       real ( kind = 8 ) x(*)

       fn=(x(1)-0.d0)**2 +(x(2)-1.d0)**2 +(x(3)-2.d0)**2
       end
!23456789




!234567890
       implicit none
       integer i,j,k,nn
       real*8 xx,yy,zz,w0,xc,pi,amp
       real*8 xi,xf,dx,fac1,fac2,tmp,tmq,sm1,sm2
       real*8 facsig,facgam,sig,gam,fwhmg,fwhml

       xc=0.d0
       amp=1.d0
       w0=1.d0
       pi=4.d0*atan(1.d0)
       fac1=amp/(sqrt(pi/2.d0)*w0)
       fac2=2.d0*amp/pi

       sig=1.d0
       gam=1.d0
       facsig=1.d0/(sig*sqrt(2.d0*pi))
       facgam=1.d0/(pi*gam)
       fwhmg=2.d0*sqrt(2.d0*log(2.d0))*sig
       fwhml=2.d0*gam

       nn=100000
       xi=-10.d0
       xf=10.d0
       dx=(xf-xi)/float(nn-1)
       sm1=0.d0 ; sm2=0.d0
       do i=1,nn
       xx=xi+float(i-1)*dx
       tmp=-2.d0*((xx-xc)/w0)**2 ; if(tmp < -50.d0) tmp=-50.d0 ; if(tmp >  50.d0) tmp=50.d0
       yy=fac1*exp(tmp)
       tmp=-(xx-xc)**2/(2.d0*sig*sig) ; if(tmp < -50.d0) tmp=-50.d0 ; if(tmp >  50.d0) tmp=50.d0
       yy=facsig*exp(tmp)
       tmq=4.d0*(xx-xc)**2 ; zz=fac2*((w0)/(w0**2+tmq))
       zz=fac2*w0/(w0**2+4.d0*(xx-xc)**2)
       zz=facgam/(1.d0+((xx-xc)/gam)**2)
       sm1=sm1+dx*yy
       sm2=sm2+dx*zz
!      write(6,*) xx,yy
       enddo
       write(6,*) sm1,sm2


       stop
       end

!234567890
       function gauss(xx,xc,sig)
       implicit none
       real*8 xx,xc,sig,gauss
       real*8 tmp,pi,facsig,fwhmg

       pi=4.d0*atan(1.d0)
       if(sig < 0.d0)then
       fwhmg=-2.d0*sqrt(2.d0*log(2.d0))*sig
       sig=fwhmg/(2.d0*sqrt(2.d0*log(2.d0)))
                     endif
       facsig=1.d0/(sig*sqrt(2.d0*pi))
       tmp=-(xx-xc)**2/(2.d0*sig*sig) ; if(tmp < -50.d0) tmp=-50.d0 ; if(tmp >  50.d0) tmp=50.d0
       gauss=facsig*exp(tmp)
       end
!234567890
       function alorentz(xx,xc,gam)
       implicit none
       real*8 xx,xc,gam,alorentz
       real*8 pi,facgam,fwhml

       pi=4.d0*atan(1.d0)
       if(gam < 0.d0)then
       fwhml=-2.d0*gam
       gam=fwhml/2.d0
                     endif
       facgam=1.d0/(pi*gam)
       alorentz=facgam/(1.d0+((xx-xc)/gam)**2)
       end

!234567890
       function gauss(xx,xc,sig0)
       implicit none
       real*8 xx,xc,sig0,gauss
       real*8 tmp,pi,facsig,fwhmg,sig

       sig=sig0
       pi=4.d0*atan(1.d0)
       if(sig < 0.d0)then
       fwhmg=-2.d0*sqrt(2.d0*log(2.d0))*sig
       sig=fwhmg/(2.d0*sqrt(2.d0*log(2.d0)))
                     endif
       facsig=1.d0/(sig*sqrt(2.d0*pi))
       tmp=-(xx-xc)**2/(2.d0*sig*sig) ; if(tmp < -50.d0) tmp=-50.d0 ; if(tmp >  50.d0) tmp=50.d0
       gauss=facsig*exp(tmp)
       end
!234567890
       function alorentz(xx,xc,gam0)
       implicit none
       real*8 xx,xc,gam0,alorentz
       real*8 pi,facgam,fwhml,gam

       gam=gam0
       pi=4.d0*atan(1.d0)
       if(gam < 0.d0)then
       fwhml=-2.d0*gam
       gam=fwhml/2.d0
                     endif
       facgam=1.d0/(pi*gam)
       alorentz=facgam/(1.d0+((xx-xc)/gam)**2)
       end


--------------------------------------------------------------

Nelder-Mead algorithm
함수 f(x)를 최소화하는 프로그램. f(x) 이 최소가 되는 x를 찾는것이다. x는 다차원의 벡터이다. 목적함수 f(x)의 미분이 정의되지 않은 경우, Nelder-Mead 알고리듬을 활용한다.

미분이 정의된 경우는 BFGS 알고리듬을 활용한다. 일반적인 문제 풀이에서 최고의 성능을 제공한다고 알려진 알고리듬이다.
https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
http://incredible.egloos.com/4838059
http://people.sc.fsu.edu/~jburkardt/f_src/asa047/asa047.f90
http://people.sc.fsu.edu/~jburkardt/f_src/asa047/asa047.html
http://people.sc.fsu.edu/~jburkardt/f_src/f_src.html

curse of dimensionality (차원의 저주)

차원이 높은 문제를 풀수록 해를 찾기가 어려워진다. 일차 미분을 사용하지 않기 때문에 해를 찾는데 더욱 어려움이 있다고 할 수 있다.

driver 설계

함수 정의, module을 사용하여 필요한 변수들을 확보한다.
Nelder-Mead 방법을 아래와 같이 복사된 서브루틴 형식으로 사용하면 편리하다. 이러할 경우, 반드시 함수를 정의해야 한다.
정의 될 함수는 반드시 f(x)형식으로 정의어야만 한다. 즉, 벡터 x 가 입력이고 반환되는 목적 함수는 f가 된다.
아래의 예제에서와 같이 응용문제 풀이에 동원된 모든 변수, 서브루틴들을 아래의 함수에서 이용할 수 있다. 이렇게 함으로써 모든 일반적인 응용문제 목적함수를 계산할 수 있게된다. 벡터 x로 표현된 독립변수들을 응용문제에 맞게 재정의하여 사용할 수 있다.
즉, 아래의 함수에서 서브루틴, 기타 타함수를 불러서 계산할 수 있다.

!234567890
       real*8 function fn(pvec)
       USE csa, ONLY : nvar
       implicit none
       real*8 pvec(nvar)
       real*8 tmp,ftn,dftn,pi
       integer i

       pi=4.d0*atan(1.d0)
       tmp=10.d0*nvar
       do i=1,nvar
       tmp=tmp+pvec(i)**2-10.d0*cos(2.d0*pi*pvec(i))
       enddo
       fn=tmp
       end

!234567890
!      Written by In-Ho Lee, KRISS, October 14, 2017.
       module neldermead
       implicit none
       private
       save
       integer n,kcount,konvge
       real*8 reqmin
       real*8, allocatable :: start(:),step(:),xmin(:)
       public :: neldermeadbody,initialneldermead,finalneldermead

       contains
!234567890
!      Written by In-Ho Lee, KRISS, October 14, 2017.
       subroutine initialneldermead(n0)
       implicit none
       integer n0
       integer i

       n=n0
       reqmin = 1.0D-08 ; konvge = 10 ; kcount = 500
       kcount=200*n0
       allocate(start(n),step(n),xmin(n))
       do i=1,n
       step(i)=1.d0
       enddo
       end subroutine initialneldermead
!234567890
!      Written by In-Ho Lee, KRISS, October 14, 2017.
       subroutine finalneldermead()
       implicit none

       deallocate(start,step,xmin)
       end subroutine finalneldermead
!234567890
!      Written by In-Ho Lee, KRISS, October 14, 2017.
       subroutine neldermeadbody(start0,rslt)
       implicit none
       real*8 start0(n),rslt
       integer icount,ifault,numres
       real*8 ynewlo
       real*8, external :: fn

       start=start0
       call nelmin(fn,n,start,xmin,ynewlo,reqmin,step,konvge,kcount,icount,numres,ifault)
       start0=xmin ; rslt=ynewlo
       end subroutine neldermeadbody
       end module neldermead
!234567890
subroutine nelmin ( fn, n, start, xmin, ynewlo, reqmin, step, konvge, kcount, &
  icount, numres, ifault )

!*****************************************************************************80
!
!! NELMIN minimizes a function using the Nelder-Mead algorithm.
!
!  Discussion:
!
!    This routine seeks the minimum value of a user-specified function.
!
!    Simplex function minimisation procedure due to Nelder and Mead (1965),
!    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
!    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
!    25, 97) and Hill(1978, 27, 380-2)
!
!    The function to be minimized must be defined by a function of
!    the form
!
!      function fn ( x )
!      real ( kind = 8 ) fn
!      real ( kind = 8 ) x(*)
!
!    and the name of this subroutine must be declared EXTERNAL in the
!    calling routine and passed as the argument FN.
!
!    This routine does not include a termination test using the
!    fitting of a quadratic surface.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    27 February 2008
!
!  Author:
!
!    Original FORTRAN77 version by R ONeill.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    John Nelder, Roger Mead,
!    A simplex method for function minimization,
!    Computer Journal,
!    Volume 7, 1965, pages 308-313.
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, external FN, the name of the function which evaluates
!    the function to be minimized.
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!    0 < N is required.
!
!    Input/output, real ( kind = 8 ) START(N).  On input, a starting point
!    for the iteration.  On output, this data may have been overwritten.
!
!    Output, real ( kind = 8 ) XMIN(N), the coordinates of the point which
!    is estimated to minimize the function.
!
!    Output, real ( kind = 8 ) YNEWLO, the minimum value of the function.
!
!    Input, real ( kind = 8 ) REQMIN, the terminating limit for the variance
!    of the function values.  0 < REQMIN is required.
!
!    Input, real ( kind = 8 ) STEP(N), determines the size and shape of the
!    initial simplex.  The relative magnitudes of its elements should reflect
!    the units of the variables.
!
!    Input, integer ( kind = 4 ) KONVGE, the convergence check is carried out
!    every KONVGE iterations. 0 < KONVGE is required.
!
!    Input, integer ( kind = 4 ) KCOUNT, the maximum number of function
!    evaluations.
!
!    Output, integer ( kind = 4 ) ICOUNT, the number of function evaluations
!    used.
!
!    Output, integer ( kind = 4 ) NUMRES, the number of restarts.
!
!    Output, integer ( kind = 4 ) IFAULT, error indicator.
!    0, no errors detected.
!    1, REQMIN, N, or KONVGE has an illegal value.
!    2, iteration terminated because KCOUNT was exceeded without convergence.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ), parameter :: ccoeff = 0.5D+00
  real ( kind = 8 ) del
  real ( kind = 8 ), parameter :: ecoeff = 2.0D+00
  real ( kind = 8 ), parameter :: eps = 0.001D+00
  real ( kind = 8 ), external :: fn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jcount
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) l
  integer ( kind = 4 ) numres
  real ( kind = 8 ) p(n,n+1)
  real ( kind = 8 ) p2star(n)
  real ( kind = 8 ) pbar(n)
  real ( kind = 8 ) pstar(n)
  real ( kind = 8 ), parameter :: rcoeff = 1.0D+00
  real ( kind = 8 ) reqmin
  real ( kind = 8 ) rq
  real ( kind = 8 ) start(n)
  real ( kind = 8 ) step(n)
  real ( kind = 8 ) x
  real ( kind = 8 ) xmin(n)
  real ( kind = 8 ) y(n+1)
  real ( kind = 8 ) y2star
  real ( kind = 8 ) ylo
  real ( kind = 8 ) ynewlo
  real ( kind = 8 ) ystar
  real ( kind = 8 ) z
!
!  Check the input parameters.
!
  if ( reqmin <= 0.0D+00 ) then
    ifault = 1
    return
  end if

  if ( n < 1 ) then
    ifault = 1
    return
  end if

  if ( konvge < 1 ) then
    ifault = 1
    return
  end if
!
!  Initialization.
!
  icount = 0
  numres = 0
  jcount = konvge
  del = 1.0D+00
  rq = reqmin * real ( n, kind = 8 )
!
!  Initial or restarted loop.
!
  do

    p(1:n,n+1) = start(1:n)
    y(n+1) = fn ( start )
    icount = icount + 1
!
!  Define the initial simplex.
!
    do j = 1, n
      x = start(j)
      start(j) = start(j) + step(j) * del
      p(1:n,j) = start(1:n)
      y(j) = fn ( start )
      icount = icount + 1
      start(j) = x
    end do
!
!  Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
!  the vertex of the simplex to be replaced.
!
    ilo = minloc ( y(1:n+1), 1 )
    ylo = y(ilo)
!
!  Inner loop.
!
    do while ( icount < kcount )
!
!  YNEWLO is, of course, the HIGHEST value???
!
      ihi = maxloc ( y(1:n+1), 1 )
      ynewlo = y(ihi)
!
!  Calculate PBAR, the centroid of the simplex vertices
!  excepting the vertex with Y value YNEWLO.
!
      do i = 1, n
        pbar(i) = ( sum ( p(i,1:n+1) ) - p(i,ihi) ) / real ( n, kind = 8 )
      end do
!
!  Reflection through the centroid.
!
      pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
      ystar = fn ( pstar )
      icount = icount + 1
!
!  Successful reflection, so extension.
!
      if ( ystar < ylo ) then

        p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
        y2star = fn ( p2star )
        icount = icount + 1
!
!  Retain extension or contraction.
!
        if ( ystar < y2star ) then
          p(1:n,ihi) = pstar(1:n)
          y(ihi) = ystar
        else
          p(1:n,ihi) = p2star(1:n)
          y(ihi) = y2star
        end if
!
!  No extension.
!
      else

        l = 0
        do i = 1, n + 1
          if ( ystar < y(i) ) then
            l = l + 1
          end if
        end do

        if ( 1 < l ) then

          p(1:n,ihi) = pstar(1:n)
          y(ihi) = ystar
!
!  Contraction on the Y(IHI) side of the centroid.
!
        else if ( l == 0 ) then

          p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
          y2star = fn ( p2star )
          icount = icount + 1
!
!  Contract the whole simplex.
!
          if ( y(ihi) < y2star ) then

            do j = 1, n + 1
              p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5D+00
              xmin(1:n) = p(1:n,j)
              y(j) = fn ( xmin )
              icount = icount + 1
            end do

            ilo = minloc ( y(1:n+1), 1 )
            ylo = y(ilo)

            cycle
!
!  Retain contraction.
!
          else
            p(1:n,ihi) = p2star(1:n)
            y(ihi) = y2star
          end if
!
!  Contraction on the reflection side of the centroid.
!
        else if ( l == 1 ) then

          p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
          y2star = fn ( p2star )
          icount = icount + 1
!
!  Retain reflection?
!
          if ( y2star <= ystar ) then
            p(1:n,ihi) = p2star(1:n)
            y(ihi) = y2star
          else
            p(1:n,ihi) = pstar(1:n)
            y(ihi) = ystar
          end if

        end if

      end if
!
!  Check if YLO improved.
!
      if ( y(ihi) < ylo ) then
        ylo = y(ihi)
        ilo = ihi
      end if

      jcount = jcount - 1

      if ( 0 < jcount ) then
        cycle
      end if
!
!  Check to see if minimum reached.
!
      if ( icount <= kcount ) then

        jcount = konvge

        x = sum ( y(1:n+1) ) / real ( n + 1, kind = 8 )
        z = sum ( ( y(1:n+1) - x )**2 )

        if ( z <= rq ) then
          exit
        end if

      end if

    end do
!
!  Factorial tests to check that YNEWLO is a local minimum.
!
    xmin(1:n) = p


핑백


최근 포토로그