Simplex Subroutine

private subroutine Simplex(varModel)

Uses

This subroutine implements the simplex method of function minimization of Nelder and Mead (1965). This ingenious method is efficient for non-linear mini- mization in a multiparameter space but is still easy to program and only requires evaluations of the objective function. The program stops the iterations whenever one of the two following criteria is reached: 1. Convergence is reached 2. The maximum number of iterations is reached.

Output: hfina: objective function at the minimum parameters of the estimated variogram

References:

Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Computer Journal 7, 308-313.

Arguments

Type IntentOptional Attributes Name
integer(kind=short), intent(in) :: varModel

type of variogram: 1 = spherical, 2 = exponential, 3 = gaussian, 4 = power

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: bmax
real(kind=float), public :: c(6)
real(kind=float), public, parameter :: cont = 0.5
real(kind=float), public, parameter :: eps = 0.01
real(kind=float), public, parameter :: expa = 2.0
integer(kind=short), public :: i
integer(kind=short), public :: imax
integer(kind=short), public :: imin
integer(kind=short), public :: j
integer, public, parameter :: maxiter = 200
real(kind=float), public :: p1(6)
real(kind=float), public :: p2(6)
real(kind=float), public :: pf(6)
real(kind=float), public :: pu(7,6)
real(kind=float), public, parameter :: refl = 0.5
real(kind=float), public :: val
real(kind=float), public :: y(7)
real(kind=float), public :: y1
real(kind=float), public :: y2
real(kind=float), public :: ymax
real(kind=float), public :: ymin

Source Code

SUBROUTINE Simplex &
!
(varModel)


USE Units, ONLY: &
    !imported parameters:
    pi

IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short), INTENT (IN) :: varModel !!type of variogram: 1 = spherical, 2 = exponential, 3 = gaussian, 4 = power


!local declarations:
INTEGER, PARAMETER :: maxiter = 200 !maximum number of iteration
REAL (KIND = float), PARAMETER :: eps = 0.01 !convergence criterion
REAL (KIND = float), PARAMETER :: refl = 0.5 !reflection factor
REAL (KIND = float), PARAMETER :: expa = 2.0 !expansion factor
REAL (KIND = float), PARAMETER :: cont = 0.5 !contraction factor
REAL (KIND = float) :: pu (7,6)
REAL (KIND = float) :: p1 (6), p2 (6), pf (6)
REAL (KIND = float) :: c (6)
REAL (KIND = float) :: y (7)
REAL (KIND = float) :: val, y1, y2
REAL (KIND = float) :: ymax, ymin, bmax
INTEGER (KIND = short) :: imax, imin
INTEGER (KIND = short) :: i, j

!--------------------------------end of declarations---------------------------

! set initial values
SELECT CASE (ipara)
    CASE (1)
        pu (1,1) = parRangeMin (1)
        pu (2,1) = parRangeMax (1)
    CASE (2) 
        pu (1,1) = parRangeMin (1)
        pu (1,2) = parRangeMin (2)
        pu (2,1) = parRangeMin (1)
        pu (2,2) = parRangeMax (2)
        pu (3,1) = parRangeMax (1)
        pu (3,2) = parRangeMax (2)
    CASE (3)
        pu (1,1) = parRangeMin (1)
        pu (1,2) = parRangeMin (2)
        pu (1,3) = parRangeMin (3)
        pu (2,1) = parRangeMin (1)
        pu (2,2) = parRangeMin (2)
        pu (2,3) = parRangeMax (3)
        pu (3,1) = parRangeMin (1)
        pu (3,2) = parRangeMax (2)
        pu (3,3) = parRangeMin (3)
        pu (4,1) = parRangeMax (1)
        pu (4,2) = parRangeMax (2)
        pu (4,3) = parRangeMax (3)
        IF ( varModel == 4) THEN
            pu (4,1) = parRangeMax (1) / 2.0
            pu (4,2) = parRangeMax (2) / 2.0
            pu (4,3) = parRangeMin (3)
        END IF   
    CASE (4)
         pu (1,1) = parRangeMin (1)
         pu (1,2) = parRangeMin (2)
         pu (1,3) = parRangeMin (3)
         pu (1,4) = parRangeMin (4)
         pu (2,1) = parRangeMin (1)
         pu (2,2) = parRangeMin (2)
         pu (2,3) = parRangeMin (3)
         pu (2,4) = parRangeMax (4)
         pu (3,1) = parRangeMin (1)
         pu (3,2) = parRangeMin (2)
         pu (3,3) = parRangeMax (3)
         pu (3,4) = parRangeMin (4)
         pu (4,1) = parRangeMin (1)
         pu (4,2) = parRangeMax (2)
         pu (4,3) = parRangeMin (3)
         pu (4,4) = parRangeMin (4)
         pu (5,1) = parRangeMax (1)
         pu (5,2) = parRangeMax (2)
         pu (5,3) = parRangeMax (3)
         pu (5,4) = parRangeMax (4)
      CASE DEFAULT 
         pu (1,1) = parRangeMin (1)
         pu (1,2) = parRangeMin (2)
         pu (1,3) = parRangeMin (3)
         pu (1,4) = parRangeMin (4)
         pu (1,5) = parRangeMin (5)
         pu (2,1) = parRangeMin (1)
         pu (2,2) = parRangeMin (2)
         pu (2,3) = parRangeMin (3)
         pu (2,4) = parRangeMin (4)
         pu (2,5) = parRangeMax (5)
         pu (3,1) = parRangeMin (1)
         pu (3,2) = parRangeMin (2)
         pu (3,3) = parRangeMin (3)
         pu (3,4) = parRangeMax (4)
         pu (3,5) = parRangeMin (5)
         pu (4,1) = parRangeMin (1)
         pu (4,2) = parRangeMin (2)
         pu (4,3) = parRangeMax (3)
         pu (4,4) = parRangeMin (4)
         pu (4,5) = parRangeMin (5)
         pu (5,1) = parRangeMin (1)
         pu (5,2) = parRangeMax (2)
         pu (5,3) = parRangeMin (3)
         pu (5,4) = parRangeMin (4)
         pu (5,5) = parRangeMin (5)
         pu (6,1) = parRangeMax (1)
         pu (6,2) = parRangeMax (2)
         pu (6,3) = parRangeMax (3)
         pu (6,4) = parRangeMax (4)
         pu (6,5) = parRangeMin (5)
 END SELECT
DO  i = 1, ipara + 1
  DO j = 1, ipara
     pf (j) = pu (i,j)
  END DO
  y (i) = Valf (pf, varmodel)
END DO  

! simplex algorithm
200 ymax = -1.0E+15
ymin = 1.0E+15
bmax = ymax

DO i = 1, ipara + 1
    
    IF ( y(i) > ymax ) THEN
        ymax = y(i)
        imax = i
    END IF
    
    IF ( y(i) < ymin ) THEN
        ymin = y(i)
        imin = i
    END IF
END DO

DO  i = 1, ipara + 1
    IF (i /= imax) THEN
       IF (y(i) > bmax) bmax = y(i)
    END IF
END DO


!     200 ITERATIONS MAXIMUM

      IF ( ieva > maxiter ) GOTO 300

!    CONTROID

DO  I=1,IPARA
    VAL=0.0
    DO  J=1,IPARA+1
        IF (J.NE.IMAX) VAL=VAL+PU(J,I)
    END DO
    C(I)=VAL/IPARA
END DO

!     P*
  
DO  I=1,IPARA
     P1(I)=(1.0+REFL)*C(I)-REFL*PU(IMAX,I)
END DO
      Y1=VALF(P1, varmodel)
      IF (Y1.GT.YMIN.AND.Y1.LT.BMAX) THEN
         DO  I=1,IPARA
            PU(IMAX,I)=P1(I)
         END DO
         Y(IMAX)=Y1
         IF (FCONV(Y,IPARA).LT.EPS) GOTO 300
         GOTO 200
      END IF
      IF (Y1.LT.YMIN) THEN
         DO 11 I=1,IPARA
            P2(I)=(1.0+EXPA)*P1(I)-EXPA*C(I)
 11      CONTINUE
         Y2=VALF(P2, varmodel)
         IF (Y2.LT.YMIN) THEN
            DO 12 I=1,IPARA
               PU(IMAX,I)=P2(I)
 12         CONTINUE
            Y(IMAX)=Y2
            IF (FCONV(Y,IPARA).LT.EPS) GOTO 300
            GOTO 200
         ELSE
            DO 13 I=1,IPARA
               PU(IMAX,I)=P1(I)
 13         CONTINUE
            Y(IMAX)=Y1
            IF (FCONV(Y,IPARA).LT.EPS) GOTO 300
            GOTO 200
         END IF
      END IF

      IF (Y1.LT.YMAX) THEN
         DO 15 I=1,IPARA
            PU(IMAX,I)=P1(I)
 15      CONTINUE
         Y(IMAX)=Y1
      END IF

      DO 16 I=1,IPARA
         P2(I)=CONT*PU(IMAX,I)+(1.0-CONT)*C(I)
 16   CONTINUE
      Y2=VALF(P2, varmodel)
      IF (Y2.GT.YMAX) THEN
         DO 17 I=1,IPARA+1
            DO 18 J=1,IPARA
               PU(I,J)=(PU(I,J)+PU(IMIN,J))/2.0
 18         CONTINUE
            DO 19 J=1,IPARA
               PF(J)=PU(I,J)
 19         CONTINUE
            Y(I)=VALF(PF, varmodel)
 17      CONTINUE
         IF (FCONV(Y,IPARA).LT.EPS) GOTO 300
         GOTO 200
      ELSE
         DO 1 I=1,IPARA
            PU(IMAX,I)=P2(I)
 1       CONTINUE
         Y(IMAX)=Y2
         IF (FCONV(Y,IPARA).LT.EPS) GOTO 300
         GOTO 200
      END IF

300    IF (IPARA.EQ.1) THEN
          range=PU(IMIN,1)
          anisotropyAngle=0.0
          anisotropyRatio=1.0
          partialSill=VAR
          nugget=0.0
       ELSE IF (IPARA.EQ.2) THEN
          IF (NPEP.EQ.1) THEN
            range=PU(IMIN,1)
            nugget=PU(IMIN,2)
            partialSill=VAR-nugget
            anisotropyAngle=0.0
            anisotropyRatio=1.0
          ELSE
            range=PU(IMIN,1)
            partialSill=PU(IMIN,2)
            nugget=0.0
            anisotropyAngle=0.0
            anisotropyRatio=1.0
            VAR=partialSill
          END IF
       ELSE IF (IPARA.EQ.3) THEN
          IF (NPEP.EQ.1) THEN
            range=PU(IMIN,1)
            partialSill=PU(IMIN,2)
            nugget=PU(IMIN,3)
            VAR=partialSill+nugget
            anisotropyAngle=0.0
            anisotropyRatio=1.0
          ELSE
            range=PU(IMIN,1)
            anisotropyAngle=PU(IMIN,2)
            anisotropyRatio=PU(IMIN,3)
            nugget=0.0
            partialSill=VAR
          END IF
       ELSE IF (IPARA.EQ.4) THEN
          IF (NPEP.EQ.1) THEN
            range=PU(IMIN,1)
            anisotropyAngle=PU(IMIN,2)
            anisotropyRatio=PU(IMIN,3)
            nugget=PU(IMIN,4)
            partialSill=VAR-nugget
          ELSE
            range=PU(IMIN,1)
            anisotropyAngle=PU(IMIN,2)
            anisotropyRatio=PU(IMIN,3)
            partialSill=PU(IMIN,4)
            nugget=0.0
            VAR=partialSill
          END IF
       ELSE
          range=PU(IMIN,1)
          anisotropyAngle=PU(IMIN,2)
          anisotropyRatio=PU(IMIN,3)
          partialSill=PU(IMIN,4)
          nugget=PU(IMIN,5)
          VAR=partialSill+nugget
       END IF
        HFINA=YMIN


RETURN
END SUBROUTINE Simplex