!
! ----------------------------------
! Sample program of PSO optimization
! ----------------------------------
!
!   Using simplified wrapper routine (with standard control parameters)
!
!----------------------------------------------------------------------
! Main Program
!----------------------------------------------------------------------
program test1p
  use opti_pso
  integer, parameter :: nmax = 100
  integer  :: n, i
  real     :: xinit(1:nmax), xmin(1:nmax), xmax(1:nmax)
  real     :: x(1:nmax), r
  integer  :: np, it, ir                 ! PSO control parameters
  external :: evtest                     ! interface declaration is better

  !====== 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)
  np = 7                                 ! population size
  it = 400                               ! maximum iteration count
  ir = 50                                ! interval of message printing

  !====== Optimize ======
  call do_pso(n, xmin, xmax, xinit, np, it, ir, pso_maximize, evtest, x, r)

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

end program test1p

!----------------------------------------------------------------------
! Evaluation Function
!----------------------------------------------------------------------
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 :: t1, t2
  t1 = x(1)**4 - 24. * x(1)**3 + 193. * x(1)**2 - 570. * x(1) + 400.
  t2 = x(2)**4 - 21. * x(2)**3 + 151. * x(2)**2 - 411. * x(2) + 280.
  ev  = -3. / 2. * (t1 + t2)
  r = 0
end function evtest