program optvorma

  use iso_c_binding, only: c_ptr,c_loc,c_f_pointer
  use optvor, only: pdata_type,evalfg,proj,approxg,drawsol,InA,compscale
  
  implicit none

  ! PARAMETERS
  integer, parameter :: nvert = 20
  real(kind=8), parameter :: PI = acos( - 1.0d0 ) 

  ! LOCAL SCALARS
  integer :: allocerr,fcnt,i,inform,iprint,iter,itrial,maxfc,maxit,maxnv,n, &
       nsites,ntrials,spginfo
  real(kind=8) :: ang,area,epsopt,f,fbest,finish,ftarget,gpsupn,marg,scale,seed,start
  type(pdata_type), target :: pdata

  ! LOCAL ARRAYS
  real(kind=8) :: p(2)
  real(kind=8), allocatable :: g(:),gtmp(:),x(:),x0(:),vol(:)
  character(len=80) :: filename

  ! FUNCTIONS
  real(kind=8) :: drand
  
  ! Set sites

  write(*,*) 'Enter nsites>0: '
  read(*,*) nsites

  ! Set disjoint convex polygons Aj whose union define A 
  
  pdata%npols = 1
  
  allocate(pdata%nvpols(pdata%npols),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  pdata%nvpols(1:pdata%npols) = (/ nvert + 1 /)
  
  maxnv = maxval( pdata%nvpols(1:pdata%npols) )
  
  allocate(pdata%poledges(pdata%npols,maxnv),pdata%xpol(pdata%npols,maxnv),pdata%ypol(pdata%npols,maxnv),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  do i = 1,nvert + 1
     ang = 2.0d0 * PI * ( i - 1 ) / nvert
     pdata%xpol(1,i) = cos( ang )
     pdata%ypol(1,i) = sin( ang )
  end do
  
  pdata%poledges(1,1:nvert) = 0

  ! Scale A so cells have ideal area near unity
  
  call compscale(nsites,scale,area,c_loc(pdata))

  pdata%volA = area 
  pdata%xpol(1,1:pdata%nvpols(1)) = scale * pdata%xpol(1,1:pdata%nvpols(1))
  pdata%ypol(1,1:pdata%nvpols(1)) = scale * pdata%ypol(1,1:pdata%nvpols(1))
  
  pdata%xmin = minval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%xmax = maxval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%ymin = minval( pdata%ypol(1,1:pdata%nvpols(1)) )
  pdata%ymax = maxval( pdata%ypol(1,1:pdata%nvpols(1)) )
  
  marg = 0.1d0 * max( 1.0d0, abs( pdata%xmin ), abs( pdata%xmax ), abs( pdata%ymin), abs( pdata%ymax ) )

  pdata%xmin = pdata%xmin - marg
  pdata%xmax = pdata%xmax + marg
  pdata%ymin = pdata%ymin - marg
  pdata%ymax = pdata%ymax + marg
  
  pdata%c(1:2) = 0.5d0 * (/ pdata%xmax + pdata%xmin, pdata%ymax + pdata%ymin /)
  pdata%rmax = norm2(pdata%c - (/ pdata%xmin, pdata%ymin /))
  pdata%drawscale = 2.5d0 / scale

  if ( nsites .ge. 10000 ) then
     pdata%drawscale = 5.0d0 / scale
  end if

  ! Goal function data

! pdata%Jweight(0:5) = 0.0d0
! pdata%Jweight(0)  = 10.0d0
! pdata%Jweight(1)  =  1.0d0
  pdata%Jweight(2)  =  1.0d0
! pdata%Jweight(3)  =  1.0d0
! pdata%Jweight(4)  =  1.0d0
! pdata%Jweight(5)  =  1.0d-04

! write(*,*) 'Enter pdata%Jweight(5): '
! read(*,*) pdata%Jweight(5)
  
! pdata%J0delta = 0.1d0
  
  pdata%J2c(1:2) = pdata%c(1:2)
  pdata%J2r      = scale

! write(*,*) 'Enter ctol for J3 (between 0 and 1): '
! read(*,*) pdata%J3ctol
  
! write(*,*) 'Enter ctol for J4 (between 0 and 1): '
! read(*,*) pdata%J4ctol
  
  ! Optimize
  
  write(*,*) 'Enter ntrials>0: '
  read(*,*) ntrials

  n = 2 * nsites
  
  allocate(x(n),x0(n),g(n),gtmp(n),vol(nsites),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if
        
  fbest = huge( 1.0d0 )
  
  do itrial = 1,ntrials
     seed = 123456.0d0 * itrial

     write(*,*) 'Processing trial itrial = ',itrial
     
     do i = 1,nsites
        p(1) = pdata%xmin + ( pdata%xmax - pdata%xmin ) * drand(seed)
        p(2) = pdata%ymin + ( pdata%ymax - pdata%ymin ) * drand(seed)

        do while ( .not. InA(p,c_loc(pdata)) )
           p(1) = pdata%xmin + ( pdata%xmax - pdata%xmin ) * drand(seed)
           p(2) = pdata%ymin + ( pdata%ymax - pdata%ymin ) * drand(seed)
        end do
       
        x(2*i-1) = p(1)
        x(2*i)   = p(2)
     end do

!!$     call evalfg(n,x,f,g,inform,c_loc(pdata))
!!$
!!$     write(*,*) 'f = ',f
!!$     write(*,*) 'g = ',g(1:n)
!!$
!!$     call approxg(n,x,gtmp,1.0d-06,inform,c_loc(pdata))
!!$     write(*,*) 'approxg = ',gtmp(1:n)
!!$
!!$     do i = 1,n
!!$        write(*,*) i,g(i),gtmp(i),abs(g(i)-gtmp(i))
!!$     end do
!!$     
!!$     write(*,*) 'Maximum absolute error = ',maxval( abs( g(1:n) - gtmp(1:n) ) )
!!$     stop
  
     x0(1:n) = x(1:n)

     epsopt  = 1.0d-08
     ftarget = 1.0d-08
     maxit =  50000
     maxfc = 500000
     iprint = 1
  
     call cpu_time(start)

!!$     call evalfg(n,x,f,g,inform,c_loc(pdata))
!!$
!!$     write(*,*) 'f = ',f
!!$     write(*,*) 'g = ',g(1:n)
!!$
!!$     call approxg(n,x,gtmp,1.0d-06,inform,c_loc(pdata))
!!$     write(*,*) 'approxg = ',gtmp(1:n)
!!$
!!$     do i = 1,n
!!$        write(*,*) i,g(i),gtmp(i),abs(g(i)-gtmp(i))
!!$     end do
!!$
!!$     write(*,*) 'Maximum absolute error = ',maxval( abs( g(1:n) - gtmp(1:n) ) )
!!$     stop
     
     call spg(n,x,epsopt,ftarget,maxit,maxfc,iprint,f,gpsupn,iter,fcnt,spginfo,inform, &
          evalfg,proj,c_loc(pdata))
  
     call cpu_time(finish)

     write(*,*) 'Statistics: f = ',f,' gpsupn = ',gpsupn,' it = ',iter,' fcnt = ',fcnt, &
          ' spginfo = ',spginfo,' time = ',finish-start
     
     if ( f .lt. fbest ) then
        fbest = f

        write(*,*) 'Solution was improved. New best f = ',fbest

        write(filename,"(A15,I0,A3)") 'fig-regpol-ini-',nsites,'.mp'
        call drawsol(n,x0,filename,c_loc(pdata))
  
        write(filename,"(A15,I0,A3)") 'fig-regpol-end-',nsites,'.mp'
        call drawsol(n,x,filename,c_loc(pdata))

        open(unit=10,file='tabline.txt')   
        write(10,1000) scale,pdata%volA,pdata%npols,nsites,itrial,f,gpsupn,iter,fcnt,spginfo,finish-start
        close(10)
        
        if ( fbest .le. ftarget ) then
           write(*,*) 'A global minimizer has been found!'
           exit
        end if
     end if
  end do
  
  deallocate(x,x0,g,gtmp,vol,stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Deallocation error.'
     stop
  end if
  
  stop

1000 format(1P,E12.5,1X,1P,E12.5,1X,I3,1X,I6,1X,I6,1X,1P,E12.5,1X,1P,E7.1,2(1X,I6),1X,I2,0P,F14.6)
  
end program optvorma