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
	
  ! LOCAL SCALARS
  integer :: allocerr,fcnt,i,inform,iprint,iter,itrial,j,maxfc,maxit,maxnv,n, &
       nsites,ntrials,spginfo, nsegs
  real(kind=8) :: area,epsopt,f,fbest,finish,ftarget,gpsupn,marg,scale,seed,start
  type(pdata_type), target :: pdata

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

  ! FUNCTIONS
  real(kind=8) :: drand

  ! Problem data
  
  ! Set sites

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

  nsegs = 4
  iter  = 1
  
  vxb(1:nsegs) = (/0.0d0, 1.0d0, 1.0d0, 0.0d0 /)
  vxb(5:64) = 0.0d0
  
  vyb(1:nsegs) = (/ 0.0d0, 0.0d0, 1.0d0, 1.0d0 /)
  vyb(5:64) = 0.0d0
  
  call includes_triang_segs(vxb,vyb,nsegs,iter)

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

  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
  
  pdata%xpol( 1,1:pdata%nvpols(1))  = (/ vxb(1), vxb(2), vxb(3), vxb(63), vxb(64), vxb(1) /)
  pdata%ypol( 1,1:pdata%nvpols(1))  = (/ vyb(1), vyb(2), vyb(3), vyb(63), vyb(64), vyb(1) /)
  
  pdata%xpol( 2,1:pdata%nvpols(2))  = (/ vxb(3), vxb(4), vxb(5), vxb(6), vxb(7), vxb(3) /)
  pdata%ypol( 2,1:pdata%nvpols(2))  = (/ vyb(3), vyb(4), vyb(5), vyb(6), vyb(7), vyb(3) /)
  
  pdata%xpol( 3,1:pdata%nvpols(3))  = (/ vxb(59), vxb(60), vxb(61), vxb(62), vxb(63), vxb(59) /)
  pdata%ypol( 3,1:pdata%nvpols(3))  = (/ vyb(59), vyb(60), vyb(61), vyb(62), vyb(63), vyb(59) /)
  
  pdata%xpol( 4,1:pdata%nvpols(4))  = (/ vxb(63), vxb(3), vxb(7), vxb(59), vxb(63) /)
  pdata%ypol( 4,1:pdata%nvpols(4))  = (/ vyb(63), vyb(3), vyb(7), vyb(59), vyb(63) /)
  
  pdata%xpol( 5,1:pdata%nvpols(5))  = (/ vxb(59), vxb(7), vxb(8), vxb(9), vxb(57), vxb(58), vxb(59) /)
  pdata%ypol( 5,1:pdata%nvpols(5))  = (/ vyb(59), vyb(7), vyb(8), vyb(9), vyb(57), vyb(58), vyb(59) /)
  
  pdata%xpol( 6,1:pdata%nvpols(6))  = (/ vxb(11), vxb(12), vxb(13), vxb(14), vxb(15), vxb(11) /)
  pdata%ypol( 6,1:pdata%nvpols(6))  = (/ vyb(11), vyb(12), vyb(13), vyb(14), vyb(15), vyb(11) /)
  
  pdata%xpol( 7,1:pdata%nvpols(7))  = (/ vxb(15), vxb(16), vxb(17), vxb(18), vxb(19), vxb(15) /)
  pdata%ypol( 7,1:pdata%nvpols(7))  = (/ vyb(15), vyb(16), vyb(17), vyb(18), vyb(19), vyb(15) /)
  
  pdata%xpol( 8,1:pdata%nvpols(8))  = (/ vxb(19), vxb(20), vxb(21), vxb(22), vxb(23), vxb(19) /)
  pdata%ypol( 8,1:pdata%nvpols(8))  = (/ vyb(19), vyb(20), vyb(21), vyb(22), vyb(23), vyb(19) /)
  
  pdata%xpol( 9,1:pdata%nvpols(9))  = (/ vxb(19), vxb(23), vxb(11), vxb(15), vxb(19) /)
  pdata%ypol( 9,1:pdata%nvpols(9))  = (/ vyb(19), vyb(23), vyb(11), vyb(15), vyb(19) /)
  
  pdata%xpol(10,1:pdata%nvpols(10)) = (/ vxb(9), vxb(10), vxb(11), vxb(23), vxb(24), vxb(25), vxb(9) /)
  pdata%ypol(10,1:pdata%nvpols(10)) = (/ vyb(9), vyb(10), vyb(11), vyb(23), vyb(24), vyb(25), vyb(9) /)
  
  pdata%xpol(11,1:pdata%nvpols(11)) = (/ vxb(27), vxb(28), vxb(29), vxb(30), vxb(31), vxb(27) /)
  pdata%ypol(11,1:pdata%nvpols(11)) = (/ vyb(27), vyb(28), vyb(29), vyb(30), vyb(31), vyb(27) /)
  
  pdata%xpol(12,1:pdata%nvpols(12)) = (/ vxb(31), vxb(32), vxb(33), vxb(34), vxb(35), vxb(31) /)
  pdata%ypol(12,1:pdata%nvpols(12)) = (/ vyb(31), vyb(32), vyb(33), vyb(34), vyb(35), vyb(31) /)
  
  pdata%xpol(13,1:pdata%nvpols(13)) = (/ vxb(35), vxb(36), vxb(37), vxb(38), vxb(39), vxb(35) /)
  pdata%ypol(13,1:pdata%nvpols(13)) = (/ vyb(35), vyb(36), vyb(37), vyb(38), vyb(39), vyb(35) /)
  
  pdata%xpol(14,1:pdata%nvpols(14)) = (/ vxb(35), vxb(39), vxb(27), vxb(31), vxb(35) /)
  pdata%ypol(14,1:pdata%nvpols(14)) = (/ vyb(35), vyb(39), vyb(27), vyb(31), vyb(35) /)
  
  pdata%xpol(15,1:pdata%nvpols(15)) = (/ vxb(27), vxb(39), vxb(40), vxb(41), vxb(25), vxb(26), vxb(27) /)
  pdata%ypol(15,1:pdata%nvpols(15)) = (/ vyb(27), vyb(39), vyb(40), vyb(41), vyb(25), vyb(26), vyb(27) /)
  
  pdata%xpol(16,1:pdata%nvpols(16)) = (/ vxb(43), vxb(44), vxb(45), vxb(46), vxb(47), vxb(43) /)
  pdata%ypol(16,1:pdata%nvpols(16)) = (/ vyb(43), vyb(44), vyb(45), vyb(46), vyb(47), vyb(43) /)
  
  pdata%xpol(17,1:pdata%nvpols(17)) = (/ vxb(47), vxb(48), vxb(49), vxb(50), vxb(51), vxb(47) /)
  pdata%ypol(17,1:pdata%nvpols(17)) = (/ vyb(47), vyb(48), vyb(49), vyb(50), vyb(51), vyb(47) /)
  
  pdata%xpol(18,1:pdata%nvpols(18)) = (/ vxb(51), vxb(52), vxb(53), vxb(54), vxb(55), vxb(51) /)
  pdata%ypol(18,1:pdata%nvpols(18)) = (/ vyb(51), vyb(52), vyb(53), vyb(54), vyb(55), vyb(51) /)
  
  pdata%xpol(19,1:pdata%nvpols(19)) = (/ vxb(51), vxb(55), vxb(43), vxb(47), vxb(51) /)
  pdata%ypol(19,1:pdata%nvpols(19)) = (/ vyb(51), vyb(55), vyb(43), vyb(47), vyb(51) /)
  
  pdata%xpol(20,1:pdata%nvpols(20)) = (/ vxb(43), vxb(55), vxb(56), vxb(57), vxb(41), vxb(42), vxb(43) /)
  pdata%ypol(20,1:pdata%nvpols(20)) = (/ vyb(43), vyb(55), vyb(56), vyb(57), vyb(41), vyb(42), vyb(43) /)
  
  pdata%xpol(21,1:pdata%nvpols(21)) = (/ vxb(41), vxb(57), vxb(9), vxb(25), vxb(41) /)
  pdata%ypol(21,1:pdata%nvpols(21)) = (/ vyb(41), vyb(57), vyb(9), vyb(25), vyb(41) /)
  
  pdata%poledges( 1,1:pdata%nvpols(1)-1)  = (/ 0, 0, 1, 0, 0 /)
  pdata%poledges( 2,1:pdata%nvpols(2)-1)  = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges( 3,1:pdata%nvpols(3)-1)  = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges( 4,1:pdata%nvpols(4)-1)  = (/ 1, 1, 1, 1 /)
  pdata%poledges( 5,1:pdata%nvpols(5)-1)  = (/ 1, 0, 0, 1, 0, 0 /)
  pdata%poledges( 6,1:pdata%nvpols(6)-1)  = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges( 7,1:pdata%nvpols(7)-1)  = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges( 8,1:pdata%nvpols(8)-1)  = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges( 9,1:pdata%nvpols(9)-1)  = (/ 1, 1, 1, 1 /)
  pdata%poledges(10,1:pdata%nvpols(10)-1) = (/ 0, 0, 1, 0, 0, 1 /)
  pdata%poledges(11,1:pdata%nvpols(11)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(12,1:pdata%nvpols(12)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(13,1:pdata%nvpols(13)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(14,1:pdata%nvpols(14)-1) = (/ 1, 1, 1, 1 /)
  pdata%poledges(15,1:pdata%nvpols(15)-1) = (/ 1, 0, 0, 1, 0, 0 /)
  pdata%poledges(16,1:pdata%nvpols(16)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(17,1:pdata%nvpols(17)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(18,1:pdata%nvpols(18)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(19,1:pdata%nvpols(19)-1) = (/ 1, 1, 1, 1 /)
  pdata%poledges(20,1:pdata%nvpols(20)-1) = (/ 1, 0, 0, 1, 0, 0 /)
  pdata%poledges(21,1:pdata%nvpols(21)-1) = (/ 1, 1, 1, 1 /)

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

  pdata%volA = area
  
  do j = 1,pdata%npols
     pdata%xpol(j,1:pdata%nvpols(j)) = scale * pdata%xpol(j,1:pdata%nvpols(j))
     pdata%ypol(j,1:pdata%nvpols(j)) = scale * pdata%ypol(j,1:pdata%nvpols(j))
  end do
  
  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)) )
  do j = 2,pdata%npols
     pdata%xmin = min( pdata%xmin, minval( pdata%xpol(j,1:pdata%nvpols(j)) ) )
     pdata%xmax = max( pdata%xmax, maxval( pdata%xpol(j,1:pdata%nvpols(j)) ) )
     pdata%ymin = min( pdata%ymin, minval( pdata%ypol(j,1:pdata%nvpols(j)) ) )
     pdata%ymax = max( pdata%ymax, maxval( pdata%ypol(j,1:pdata%nvpols(j)) ) )
  end do
  
  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 = 5.0d0 / scale

  ! Goal function data

  pdata%Jweight(0:5) = 0.0d0
  pdata%Jweight(0)  = 10.0d0
  pdata%Jweight(1)  =  1.0d0
  
  pdata%J0delta = 0.1d0
  
  ! Optimize
  
  write(*,*) 'Enter ntrials>0: '
  read(*,*) ntrials
  
  n = 2 * nsites
  
  allocate(x(n),x0(n),g(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

     x0(1:n) = x(1:n)

     epsopt  = 1.0d-08
     ftarget = 1.0d-08
     maxit = 10000
     maxfc = 10000
     iprint = 0
  
     call cpu_time(start)
  
     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-cesaro-ini-',nsites,'.mp'
        call drawsol(n,x0,filename,c_loc(pdata))
  
        write(filename,"(A15,I0,A3)") 'fig-cesaro-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,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)
  
contains

  ! ******************************************************************
  ! ******************************************************************
  recursive subroutine includes_triang_segs(vxb,vyb,nsegs,iter)

    implicit none
        
    ! PARAMETERS
    real(kind=8), parameter :: PI = dacos( - 1.0d0 )
        
    ! SCALAR ARGUMENTS
    integer, intent(in) :: iter,nsegs

    ! ARRAY ARGUMENTS
    real(kind=8), intent(inout) :: vxb(64),vyb(64)

    ! LOCAL SCALARS
    integer :: i,ind
    real(kind=8) :: rotv(2)

    ! LOCAL ARRAYS
    real(kind=8) :: mthirdx(2),mthirdy(2),auxvxb(64),auxvyb(64)
        
    auxvxb(1:64) = vxb(1:64)
    auxvyb(1:64) = vyb(1:64)
        
    ind = 1
    do i = 1,nsegs
       if ( i .lt. nsegs ) then
          mthirdx(1) = auxvxb(i) + (1.0d0/3.0d0)*(auxvxb(i+1) - auxvxb(i))
          mthirdx(2) = auxvxb(i) + (2.0d0/3.0d0)*(auxvxb(i+1) - auxvxb(i))
          mthirdy(1) = auxvyb(i) + (1.0d0/3.0d0)*(auxvyb(i+1) - auxvyb(i))
          mthirdy(2) = auxvyb(i) + (2.0d0/3.0d0)*(auxvyb(i+1) - auxvyb(i))
                
          call vector_rotate_base_2d( (/ mthirdx(2), mthirdy(2) /), (/ mthirdx(1), mthirdy(1) /), PI / 3.0d0, rotv )
                
       else
          mthirdx(1) = auxvxb(i) + (1.0d0/3.0d0)*(auxvxb(1) - auxvxb(i))
          mthirdx(2) = auxvxb(i) + (2.0d0/3.0d0)*(auxvxb(1) - auxvxb(i))
          mthirdy(1) = auxvyb(i) + (1.0d0/3.0d0)*(auxvyb(1) - auxvyb(i))
          mthirdy(2) = auxvyb(i) + (2.0d0/3.0d0)*(auxvyb(1) - auxvyb(i))
                
          call vector_rotate_base_2d( (/ mthirdx(2), mthirdy(2) /), (/ mthirdx(1), mthirdy(1) /), PI / 3.0d0, rotv )
       end if
            
       vxb(ind) = auxvxb(i)
       vyb(ind) = auxvyb(i)
       ind = ind + 1

       vxb(ind) = mthirdx(1)
       vyb(ind) = mthirdy(1)
       ind = ind + 1

       vxb(ind) = rotv(1)
       vyb(ind) = rotv(2)
       ind = ind + 1

       vxb(ind) = mthirdx(2)
       vyb(ind) = mthirdy(2)
       ind = ind + 1
    end do
        
    if ( iter .lt. 2 ) then
       call includes_triang_segs(vxb,vyb,nsegs*4,iter+1)
    end if
        
  end subroutine includes_triang_segs

end program optvorma
