MODULE de_opt
  USE de_kind, ONLY : prec  ! Laskentatarkkuus
  IMPLICIT NONE
  PRIVATE
  INTEGER, PARAMETER :: step = 10
  PUBLIC :: initialize, optimize, prec
CONTAINS
  SUBROUTINE initialize(f, pop, cost, low, up)
    IMPLICIT NONE
    INTERFACE
      FUNCTION f(x) RESULT(res)
        USE de_kind, ONLY : prec
        REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: x
        REAL(KIND=prec) :: res
      END FUNCTION f
    END INTERFACE
    REAL(KIND=prec), DIMENSION(:,:), INTENT(OUT) :: pop
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: cost
    REAL(KIND=prec), INTENT(IN) :: low, up
    INTEGER :: pop_size, i

    pop_size = SIZE(pop,2)
    IF (SIZE(cost) /= pop_size) &
        STOP 'Kustannusvektori väärän kokoinen!'
    CALL RANDOM_NUMBER(pop)
    pop = low + (up - low)*pop
    DO i = 1, pop_size
      cost(i) = f(pop(:,i))
    END DO
  END SUBROUTINE initialize

  SUBROUTINE optimize(f, pop, cost, max_iter, coeff, crossover)
    IMPLICIT NONE
    INTERFACE
      FUNCTION f(x) RESULT(res)
        USE de_kind, ONLY : prec
        REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: x
        REAL(KIND=prec) :: res
      END FUNCTION f
    END INTERFACE
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: pop
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: cost
    INTEGER, INTENT(IN) :: max_iter
    REAL(KIND=prec), INTENT(IN) :: coeff, crossover
    INTEGER :: n, pop_size, allocstat, iter, i, idx1, idx2, idx3
    REAL(KIND=prec), DIMENSION(:,:), ALLOCATABLE :: new_pop
    REAL(KIND=prec), DIMENSION(SIZE(pop,1)) :: trial
    REAL(KIND=prec) :: score
    LOGICAL, DIMENSION(SIZE(pop,1)) :: mask

    n = SIZE(pop,1)
    pop_size = SIZE(pop,2)
    WRITE (*,'(A,I2,/,A,I4,/,A,F5.3,/,A,F5.3)') &
        'Dimensio: ', n, &
        'Pop. koko: ', pop_size, &
        'Kerroin: ', coeff, &
        'Risteytystodenn: ', crossover
    IF (SIZE(cost) /= pop_size) &
        STOP 'Vektori cost väärää kokoa!'
    ALLOCATE(new_pop(n,pop_size), STAT=allocstat)
    IF (allocstat /= 0) STOP 'Virhe muistivarauksessa!'

    WRITE (*,'(A,I4,G14.6)') &
        'Iteraatio, minimi: ', 0, MINVAL(cost)

    DO iter = 1, max_iter
      DO i = 1, pop_size
        CALL triple(i, pop_size, idx1, idx2, idx3)
        mask = rnd(n) < crossover
        mask(idx(n)) = .TRUE.
        WHERE (mask)
          trial = pop(:,idx3) + coeff*(pop(:,idx1)-pop(:,idx2))
        ELSEWHERE
          trial = pop(:,i)
        END WHERE
        score = f(trial)
        IF (score < cost(i)) THEN
          new_pop(:,i) = trial
          cost(i) = score
        ELSE
          new_pop(:,i) = pop(:,i)
        END IF
      END DO
      pop = new_pop
      IF (MOD(iter,step) == 0) THEN
        WRITE (*,'(A,I4,G14.6)') 'Iteraatio, minimi: ', &
            iter, MINVAL(cost)
      END IF
    END DO
  END SUBROUTINE optimize

  SUBROUTINE triple(i, n, i1, i2, i3)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, n
    INTEGER, INTENT(OUT) :: i1, i2, i3
    DO
      i1 = idx(n)
      IF (i1 /= i) EXIT
    END DO
    DO
      i2 = idx(n)
      IF (i2 /= i .AND. i2 /= i1) EXIT
    END DO
    DO
      i3 = idx(n)
      IF (ALL(i3 /= (/i,i1,i2/))) EXIT
    END DO
  END SUBROUTINE triple

  FUNCTION rnd(n)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n
    REAL(KIND=prec), DIMENSION(n) :: rnd
    CALL RANDOM_NUMBER(rnd)
  END FUNCTION rnd

  INTEGER FUNCTION idx(n)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n
    REAL(KIND=prec) :: x
    CALL RANDOM_NUMBER(x)
    idx = n*x + 1
  END FUNCTION idx
END MODULE de_opt