module optimi

! Find the minimum of a function using the Nelder-Mead polytope
! or simplex method. Program adapted from Numerical recipes.
! call opt(f, p, ndim)
! f(p, ndim) is the target function to be optimized.
! ndim is the dimension of the space, i.e. number of variables.
! p(ndim+1, ndim)  contains the coordinates of the ndim+1 apices
! of the simplex; the initial values should be chosen in such a way
! that the minimum is inside the simplex limited by the apices.
! At return the simplex has shrunk to a very small region
! surrounding the optimum.

implicit none
private
real, parameter ::    FTOL = 0.001, & ! termination accuracy
                      ALPHA = 1.0, &  ! reflection coefficient
                      BETA = 0.5, &   ! reduction coefficient
                      GAMMA = 2.0     ! expansion coefficient
integer, parameter :: NMAX = 1000     ! max number of steps
integer :: ilo,ihi,inhi, nfunk

public :: opt

contains

subroutine opt(funk, p, ndim)
  interface
     real function funk(p, ndim)
       integer, intent(in) :: ndim
       real, intent(in) :: p(ndim)
     end function
  end interface

  integer, intent(in) :: ndim         ! dimension the space
  real, dimension(ndim+1, ndim) :: p
  real ::  ytry, ysave, sum, rtol, psum(ndim), &
       y(ndim+1)

  integer, parameter ::  mpts = 3     ! nr of apices of the simplex
  integer :: i,j

  ! function values at the apices of the simplex
  do i=1,mpts
    y(i) = funk(p(i,:), ndim)
  end do

  nfunk=0

  ! sums of coordinates; the centre of the mass of the side
  ! opposite to the apex i is (psum-p(i))/ndim
  do j=1,ndim
    sum=0.0
    do i=1,mpts
      sum = sum+p(i,j) 
    end do
    psum(j) = sum
  end do

  main: do
    ! find the best (ilo), worst (ihi) and
    ! second worst (inhi) apex
    ilo=1
    if (y(1) > y(2)) then
       inhi=2; ihi=1
    else
       inhi=1; ihi=2
    end if
    do i=1,mpts
      if (y(i) < y(ilo)) ilo=i
      if (y(i) > y(ihi)) then
         inhi=ihi; ihi=i
      else if (y(i) > y(inhi)) then
         if (i /= ihi) inhi=i
      end if
    end do

    ! if the worst and best almost equal, finish
    rtol=abs(y(ihi)-y(ilo)) 

    ! the original program calculates the relative error
    ! by dividing rtol lby (abs(y(ihi))+abs(y(ilo)))
    ! this does not work, if y=0

    if (rtol < FTOL) then
           write(6,'("rtol=",f10.5)') rtol
           write(6,'("nfunk=",i5)') nfunk
           exit main
    end if
    ! if too many iterations, terminate
    if (nfunk >= NMAX) then
           write(6,'("nfunk=",i5)') nfunk
           exit main
    end if

    ! reflection
    ytry=amotry(-ALPHA, funk, p, psum, y, ndim)

    ! if the result was improved move the new apex even further
    if (ytry <= y(ilo)) then
       ytry=amotry(GAMMA,  funk, p, psum, y, ndim)
    else if (ytry >= y(inhi)) then
      ! new apex worse than the second worse;
      ! shrink the simplex
      ysave=y(ihi)
      ytry=amotry(BETA,  funk, p, psum, y, ndim)
      if (ytry >= ysave) then
        ! still too big; shrink the simplex w.r.t the best apex
        do i=1,mpts
          if (i /= ilo) then
            do j=1,ndim
              psum(j)=0.5*(p(i,j)+p(ilo,j))
              p(i,j)=psum(j)
            end do
            y(i)=funk(psum, ndim)              
          end if
        end do
        nfunk = nfunk+ndim
        ! update sums of coordinates
        do j=1,ndim
          sum=0.0
          do i=1,mpts
            sum = sum+p(i,j) 
          end do
          psum(j) = sum
        end do
      end if
    end if 
  end do main

  end subroutine

  real function amotry(fac, funk, p, psum, y, ndim)
  ! modification step 
  ! the worst apex is reflected w.r.t the opposite side
  ! and the distance is changed by the factor fac
    real, intent(in) :: fac
    interface
       real function funk(p, ndim)
         integer, intent(in) :: ndim
         real, intent(in) :: p(ndim)
       end function
    end interface
    integer, intent(in) :: ndim         ! dimension of the space
    real, dimension(ndim+1, ndim) :: p
    real, dimension(ndim) ::  psum(ndim)
    real, dimension (ndim+1) :: y

    real fac2
    integer j
    real fac1, ytry, ptry(ndim)

    fac1=(1.0-fac)/ndim
    fac2=fac-(1-fac)/ndim
    do j=1,ndim
       ptry(j)=psum(j)*fac1+p(ihi,j)*fac2
    end do
    ytry=funk(ptry, ndim)
    nfunk = nfunk+1
    if (ytry < y(ihi)) then
      y(ihi)=ytry
      do j=1,ndim
        psum(j) = psum(j)+ptry(j)-p(ihi,j)
        p(ihi,j)=ptry(j)
      end do
    end if
    amotry=ytry
  end function

end module