module optimi
! funktion minimin haku Nelderin-Meadin polytooppi- eli simpleksimenetelmällä
! algoritmi kirjasta Numerical recipes
! kutsu opt(f, p, ndim)
! f(p, ndim) on optimoitava kohdefunktio
! ndim on avaruuden dimensio eli muuttujien määrä
! p(ndim+1, ndim)  sisältää simpleksin ndim+1 kärkipisteen
!   koordinaatit; alkuarvot on syytä valita niin, että
!   minimi on pisteiden rajoittaman simpleksin sisällä
!   aliohjelman päättyessä simpleksi on kuroutunut hyvin
!   pieneksi alueeksi optimin ympärillä

implicit none
private
real, parameter ::    FTOL = 0.001, & ! lopetustarkkuus
                      ALPHA = 1.0, &  ! peilauskerroin
                      BETA = 0.5, &   ! pienennyskerroin
                      GAMMA = 2.0     ! venytyskerroin
integer, parameter :: NMAX = 1000     ! askelten maksimimaara
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         ! avaruuden dimensio
  real, dimension(ndim+1, ndim) :: p
  real ::  ytry, ysave, sum, rtol, psum(ndim), &
       y(ndim+1)

  integer, parameter ::  mpts = 3     ! simpleksin karkien maara
  integer :: i,j

  ! funktion arvot simpleksin karjissa
  do i=1,mpts
    y(i) = funk(p(i,:), ndim)
  end do

  nfunk=0

  ! koordinaattien summat; karjen i vastakkaisen 
  ! sivun painopiste on (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
    ! haetaan paras (ilo), huonoin (ihi) ja
    ! toiseksi huonoin (inhi) karkipiste
    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

    ! jos huonoin ja paras arvo lahes samoja, lopetetaan
    rtol=abs(y(ihi)-y(ilo)) 

    ! alkuperäisessä ohjelmassa lasketaan suhteellinen virhe
    ! jakamalla rtol luvulla (abs(y(ihi))+abs(y(ilo)))
    ! tämä ei toimi, jos y=0

    if (rtol < FTOL) then
           write(6,'("rtol=",f10.5)') rtol
           write(6,'("nfunk=",i5)') nfunk
           exit main
    end if
    ! jos iteraatiokierroksia liikaa, lopetetaan
    if (nfunk >= NMAX) then
           write(6,'("nfunk=",i5)') nfunk
           exit main
    end if

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

    ! jos tulos parani, siirretaan uutta karkea viela kauemmas
    if (ytry <= y(ilo)) then
       ytry=amotry(GAMMA,  funk, p, psum, y, ndim)
    else if (ytry >= y(inhi)) then
      ! uusi piste huonompi kuin toiseksi huonoin;
      ! pienennetaan simpleksia
      ysave=y(ihi)
      ytry=amotry(BETA,  funk, p, psum, y, ndim)
      if (ytry >= ysave) then
        ! edelleen liian iso; pienennetaan simpleksia 
        ! parhaan pisteen suhteen
        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
        ! paivitetaan koordinaattien summat
        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)
  ! simpleksin muunnosaskel  
  ! huonoin karki peilataan vastakkaisen sivun suhteen ja
  ! etaisyytta muutetaan kertoimella 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         ! avaruuden dimensio
    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