!
! -------------------------------------
! Sample program of SCE-UA optimization
! -------------------------------------
!
!   Advanced example
!
!     1. Separate library independent, algorithm independent routines
!         library independent: sharable with other optimization libraries
!         algorithm independent: sharable with SCE-UA module
!     2. Use low-level primitive routine
!     3. Use user-defined callback function for
!         1) history saving
!         2) convergence checking as exit criterion
!
!----------------------------------------------------------------------
! Module for Common Data (library dependent, algorithm independent)
!----------------------------------------------------------------------
module common_data
  integer, parameter, public :: itmax = 100000
  integer, public :: it0(itmax)     ! history of iteration count
  integer, public :: iev0(itmax)    ! history of evaluation count
  real,    public :: ev0(itmax)     ! history of evaluation value
  integer, public :: nt = 0         ! number of history data
end module common_data

!----------------------------------------------------------------------
! Main Program (library independent, algorithm independent)
!----------------------------------------------------------------------
program test3s
  use common_data
  integer, parameter :: nmax = 100
  integer  :: n, i
  real     :: xinit(1:nmax), xmin(1:nmax), xmax(1:nmax)
  real     :: x(1:nmax), r

  !====== Initialize parameters ======
  n = 2                                  ! number of parameter
  xmin(1) = 0.0                          ! lower bound of x(1)
  xmax(1) = 12.0                         ! upper bound of x(1)
  xmin(2) = 0.0                          ! lower bound of x(2)
  xmax(2) = 12.0                         ! upper bound of x(2)
  xinit(1) = 2.43                        ! initial value of x(1)
  xinit(2) = 7.57                        ! initial value of x(2)

  !====== Optimize ======
  call optimize(n, xmin, xmax, xinit, x, r)

  !====== Print result ======
  print *
  do i=1, n
    print *, 'P', i, x(i)
  enddo
  print *, 'E', r
  print *

  !====== Print history ======
  do i=1, nt
    print *, '---', it0(i), ev0(i), iev0(i)
  end do
  print *

end program test3s

!----------------------------------------------------------------------
! Objective Function (library independent, algorithm independent)
!
!     If you want to share some data with other routines,
!       make a module which contains public variables like "common_data",
!       referable from other routines.
!
!     However, when you enable OpenMP, this function is parallelized.
!       Never update shared data without locking from this function.
!
!     When OpenMP is enabled, you can write to only local (automatic)
!       variables allocated in this function like "t1" and "t2".
!
!----------------------------------------------------------------------
function f(x, y) result(z)
  real, intent(in) :: x, y
  real             :: z
  real :: t1, t2
  t1 = x**4 - 24. * x**3 + 193. * x**2 - 570. * x + 400.
  t2 = y**4 - 21. * y**3 + 151. * y**2 - 411. * y + 280.
  z  = -3. / 2. * (t1 + t2)
end function f

!----------------------------------------------------------------------
! Evaluation Function (library dependent, algorithm independent)
!----------------------------------------------------------------------
function evtest(x, ev) result(r)
  real, intent(in)  :: x(:)      ! parameter list
  real, intent(out) :: ev        ! evaluation value
  integer           :: r         ! number of constrained condition
  real :: f
  ev = f(x(1), x(2))
  r = 0
end function evtest

!----------------------------------------------------------------------
! Callback Function (library dependent, algorithm independent)
!
!     This function is not parallelized during OpenMP enabled,
!       you can update shared data safely.
!
!     If you don't want to use user-defined callback function,
!       please use pre-defined function "sce_cbdummy" in module "opti_sceua".
!
!----------------------------------------------------------------------
function mycallback(it, iev, ev, x) result(r)
  use common_data
  integer, intent(in) :: it      ! current iteration counter
  integer, intent(in) :: iev     ! current evaluation counter
  real,    intent(in) :: ev      ! current evaluation value
  real,    intent(in) :: x(:)    ! current parameter value
  integer             :: r       ! continuation flag (>=0:cont., <0:exit)
  real, parameter :: eps = 1.e-5

  !====== Update history ======
  nt = nt + 1
  it0(nt) = it
  iev0(nt) = iev
  ev0(nt) = ev

  !====== Check Convergence ======
  if(nt > 2 .and. abs(ev - ev0(max(nt - 2, 1))) < eps) then
    r = -1
  else
    r = 1
  end if

  !====== Print intermediate status ======
  print *, it, ev, '(', x(1), ',', x(2), ')'

end function mycallback

!----------------------------------------------------------------------
! Optimization Routine (library dependent, algorithm dependent)
!----------------------------------------------------------------------
subroutine optimize(n, xmin, xmax, xinit, x, r)
  use opti_sceua
  integer, intent(in)  :: n
  real,    intent(in)  :: xmin(1:n), xmax(1:n), xinit(1:n)
  real,    intent(out) :: x(1:n), r
  integer  :: nc, nm, ne, alpha, beta, it, ie, ir, ic, flag
  external :: evtest, mycallback     ! interface declaration is better

  !====== Setup control parameters =====
  nc = 3
  nm = 2 * n + 1           ! member in a complex (default:2*n+1)
  ne = n + 1               ! evolving member in a complex (default:n+1)
  alpha = 1                ! number of alpha iteration (default:1)
  beta = 2 * n + 1         ! number of beta iteration (default:2*n+1)
  it = 200                 ! maximum iteration count (<=0 means unlimitted)
  ie = 0                   ! maximum evaluation count (<=0 means unlimitted)
  ir = 0                   ! interval of message printing (<=0 means no message)
  ic = 10                  ! interval of callback (<=0 means no callback)
  flag = 0                 ! control flag
  flag = flag + f_quiet        ! restrain messages
  flag = flag + f_initp        ! use initial parameters
  flag = flag + f_widemutate   ! wide range mutation (recommended)
! flag = flag + f_qmcinit      ! initialize by Quasi-Monte Carlo

  !====== Optimize ======
  call sceua(n, xmin, xmax, xinit, nc, nm, ne, alpha, beta, &
&            it, ie, ir, ic, flag, sce_maximize, evtest, mycallback, x, r)

end subroutine optimize