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

  ! LOCAL ARRAYS
  real(kind=8) :: p(2)
  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

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

  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))  = (/ put(1,80), put(1,44), put(1,78), put(1,124), put(1,45), put(1,79), put(1,80) /)
  pdata%ypol( 1,1:pdata%nvpols(1))  = (/ put(2,80), put(2,44), put(2,78), put(2,124), put(2,45), put(2,79), put(2,80) /)

  pdata%xpol( 2,1:pdata%nvpols(2))  = (/ put(1,80), put(1,79), put(1,46), put(1,43), put(1,80) /)
  pdata%ypol( 2,1:pdata%nvpols(2))  = (/ put(2,80), put(2,79), put(2,46), put(2,43), put(2,80) /)

  pdata%xpol( 3,1:pdata%nvpols(3))  = (/ put(1,83), put(1,35), put(1,64), put(1,83) /)
  pdata%ypol( 3,1:pdata%nvpols(3))  = (/ put(2,83), put(2,35), put(2,64), put(2,83) /)

  pdata%xpol( 4,1:pdata%nvpols(4))  = (/ put(1,46), put(1,83), put(1,64), put(1,102), put(1,81), put(1,82), put(1,43), put(1,46) /)
  pdata%ypol( 4,1:pdata%nvpols(4))  = (/ put(2,46), put(2,83), put(2,64), put(2,102), put(2,81), put(2,82), put(2,43), put(2,46) /)

  pdata%xpol( 5,1:pdata%nvpols(5))  = (/ put(1,64), put(1,63), put(1,101), put(1,102), put(1,64) /)
  pdata%ypol( 5,1:pdata%nvpols(5))  = (/ put(2,64), put(2,63), put(2,101), put(2,102), put(2,64) /)

  pdata%xpol( 6,1:pdata%nvpols(6))  = (/ put(1,102), put(1,40), put(1,81), put(1,102) /)
  pdata%ypol( 6,1:pdata%nvpols(6))  = (/ put(2,102), put(2,40), put(2,81), put(2,102) /)

  pdata%xpol( 7,1:pdata%nvpols(7))  = (/ put(1,40), put(1,127), put(1,81), put(1,40) /)
  pdata%ypol( 7,1:pdata%nvpols(7))  = (/ put(2,40), put(2,127), put(2,81), put(2,40) /)

  pdata%xpol( 8,1:pdata%nvpols(8))  = (/ put(1,40), put(1,41), put(1,42), put(1,127), put(1,40) /)
  pdata%ypol( 8,1:pdata%nvpols(8))  = (/ put(2,40), put(2,41), put(2,42), put(2,127), put(2,40) /)

  pdata%xpol( 9,1:pdata%nvpols(9))  = (/ put(1,63), put(1,97), put(1,36), put(1,100), put(1,101), put(1,63) /)
  pdata%ypol( 9,1:pdata%nvpols(9))  = (/ put(2,63), put(2,97), put(2,36), put(2,100), put(2,101), put(2,63) /)

  pdata%xpol(10,1:pdata%nvpols(10)) = (/ put(1,100), put(1,37), put(1,38), put(1,39), put(1,101), put(1,100) /)
  pdata%ypol(10,1:pdata%nvpols(10)) = (/ put(2,100), put(2,37), put(2,38), put(2,39), put(2,101), put(2,100) /)

  pdata%xpol(11,1:pdata%nvpols(11)) = (/ put(1,36), put(1,99), put(1,100), put(1,36) /)
  pdata%ypol(11,1:pdata%nvpols(11)) = (/ put(2,36), put(2,99), put(2,100), put(2,36) /)

  pdata%xpol(12,1:pdata%nvpols(12)) = (/ put(1,35), put(1,120), put(1,56), put(1,121), put(1,33), put(1,105), put(1,106), &
                                         put(1,96), put(1,97), put(1,63), put(1,64), put(1,35) /)
  pdata%ypol(12,1:pdata%nvpols(12)) = (/ put(2,35), put(2,120), put(2,56), put(2,121), put(2,33), put(2,105), put(2,106), &
                                         put(2,96), put(2,97), put(2,63), put(2,64), put(2,35) /)

  pdata%xpol(13,1:pdata%nvpols(13)) = (/ put(1,106), put(1,60), put(1,96), put(1,106) /)
  pdata%ypol(13,1:pdata%nvpols(13)) = (/ put(2,106), put(2,60), put(2,96), put(2,106) /)

  pdata%xpol(14,1:pdata%nvpols(14)) = (/ put(1,62), put(1,98), put(1,61), put(1,60), put(1,62) /)
  pdata%ypol(14,1:pdata%nvpols(14)) = (/ put(2,62), put(2,98), put(2,61), put(2,60), put(2,62) /)

  pdata%xpol(15,1:pdata%nvpols(15)) = (/ put(1,106), put(1,94), put(1,62), put(1,60), put(1,106) /)
  pdata%ypol(15,1:pdata%nvpols(15)) = (/ put(2,106), put(2,94), put(2,62), put(2,60), put(2,106) /)

  pdata%xpol(16,1:pdata%nvpols(16)) = (/ put(1,106), put(1,34), put(1,94), put(1,106) /)
  pdata%ypol(16,1:pdata%nvpols(16)) = (/ put(2,106), put(2,34), put(2,94), put(2,106) /)

  pdata%xpol(17,1:pdata%nvpols(17)) = (/ put(1,106), put(1,117), put(1,107), put(1,85), put(1,34), put(1,106) /)
  pdata%ypol(17,1:pdata%nvpols(17)) = (/ put(2,106), put(2,117), put(2,107), put(2,85), put(2,34), put(2,106) /)

  pdata%xpol(18,1:pdata%nvpols(18)) = (/ put(1,107), put(1,86), put(1,113), put(1,114), put(1,89), put(1,85), put(1,107) /)
  pdata%ypol(18,1:pdata%nvpols(18)) = (/ put(2,107), put(2,86), put(2,113), put(2,114), put(2,89), put(2,85), put(2,107) /)

  pdata%xpol(19,1:pdata%nvpols(19)) = (/ put(1,114), put(1,31), put(1,32), put(1,87), put(1,89), put(1,114) /)
  pdata%ypol(19,1:pdata%nvpols(19)) = (/ put(2,114), put(2,31), put(2,32), put(2,87), put(2,89), put(2,114) /)

  pdata%xpol(20,1:pdata%nvpols(20)) = (/ put(1,87), put(1,115), put(1,88), put(1,103), put(1,89), put(1,87) /)
  pdata%ypol(20,1:pdata%nvpols(20)) = (/ put(2,87), put(2,115), put(2,88), put(2,103), put(2,89), put(2,87) /)

  pdata%xpol(21,1:pdata%nvpols(21)) = (/ put(1,33), put(1,118), put(1,119), put(1,105), put(1,33) /)
  pdata%ypol(21,1:pdata%nvpols(21)) = (/ put(2,33), put(2,118), put(2,119), put(2,105), put(2,33) /)

  pdata%xpol(22,1:pdata%nvpols(22)) = (/ put(1,118), put(1,104), put(1,116), put(1,119), put(1,118) /)
  pdata%ypol(22,1:pdata%nvpols(22)) = (/ put(2,118), put(2,104), put(2,116), put(2,119), put(2,118) /)

  pdata%xpol(23,1:pdata%nvpols(23)) = (/ put(1,31), put(1,92), put(1,91), put(1,90), put(1,32), put(1,31) /)
  pdata%ypol(23,1:pdata%nvpols(23)) = (/ put(2,31), put(2,92), put(2,91), put(2,90), put(2,32), put(2,31) /)

  pdata%xpol(24,1:pdata%nvpols(24)) = (/ put(1,92), put(1,93), put(1,30), put(1,29), put(1,91), put(1,92) /)
  pdata%ypol(24,1:pdata%nvpols(24)) = (/ put(2,92), put(2,93), put(2,30), put(2,29), put(2,91), put(2,92) /)

  pdata%xpol(25,1:pdata%nvpols(25)) = (/ put(1,49), put(1,50), put(1,51), put(1,52), put(1,49) /)
  pdata%ypol(25,1:pdata%nvpols(25)) = (/ put(2,49), put(2,50), put(2,51), put(2,52), put(2,49) /)

  pdata%xpol(26,1:pdata%nvpols(26)) = (/ put(1,28), put(1,27), put(1,128), put(1,59), put(1,68), put(1,76), put(1,48), &
                                         put(1,129), put(1,28) /)
  pdata%ypol(26,1:pdata%nvpols(26)) = (/ put(2,28), put(2,27), put(2,128), put(2,59), put(2,68), put(2,76), put(2,48), &
                                         put(2,129), put(2,28) /)

  pdata%xpol(27,1:pdata%nvpols(27)) = (/ put(1,73), put(1,77), put(1,71), put(1,72), put(1,73) /)
  pdata%ypol(27,1:pdata%nvpols(27)) = (/ put(2,73), put(2,77), put(2,71), put(2,72), put(2,73) /)

  pdata%xpol(28,1:pdata%nvpols(28)) = (/ put(1,129), put(1,47), put(1,125), put(1,29), put(1,30), put(1,28), put(1,129) /)
  pdata%ypol(28,1:pdata%nvpols(28)) = (/ put(2,129), put(2,47), put(2,125), put(2,29), put(2,30), put(2,28), put(2,129) /)

  pdata%xpol(29,1:pdata%nvpols(29)) = (/ put(1,125), put(1,17), put(1,19), put(1,123), put(1,24), put(1,29), put(1,125) /)
  pdata%ypol(29,1:pdata%nvpols(29)) = (/ put(2,125), put(2,17), put(2,19), put(2,123), put(2,24), put(2,29), put(2,125) /)

  pdata%xpol(30,1:pdata%nvpols(30)) = (/ put(1,19), put(1,122), put(1,123), put(1,19) /)
  pdata%ypol(30,1:pdata%nvpols(30)) = (/ put(2,19), put(2,122), put(2,123), put(2,19) /)

  pdata%xpol(31,1:pdata%nvpols(31)) = (/ put(1,59), put(1,130), put(1,57), put(1,58), put(1,2), put(1,69), put(1,59) /)
  pdata%ypol(31,1:pdata%nvpols(31)) = (/ put(2,59), put(2,130), put(2,57), put(2,58), put(2,2), put(2,69), put(2,59) /)

  pdata%xpol(32,1:pdata%nvpols(32)) = (/ put(1,69), put(1,67), put(1,59), put(1,69) /)
  pdata%ypol(32,1:pdata%nvpols(32)) = (/ put(2,69), put(2,67), put(2,59), put(2,69) /)

  pdata%xpol(33,1:pdata%nvpols(33)) = (/ put(1,67), put(1,66), put(1,4), put(1,59), put(1,67) /)
  pdata%ypol(33,1:pdata%nvpols(33)) = (/ put(2,67), put(2,66), put(2,4), put(2,59), put(2,67) /)

  pdata%xpol(34,1:pdata%nvpols(34)) = (/ put(1,4), put(1,3), put(1,5), put(1,68), put(1,59), put(1,4) /)
  pdata%ypol(34,1:pdata%nvpols(34)) = (/ put(2,4), put(2,3), put(2,5), put(2,68), put(2,59), put(2,4) /)
  
  pdata%poledges( 1,1:pdata%nvpols(1)-1)  = (/ 0, 0, 0, 0, 0, 1 /)
  pdata%poledges( 2,1:pdata%nvpols(2)-1)  = (/ 1, 0, 1, 0 /)
  pdata%poledges( 3,1:pdata%nvpols(3)-1)  = (/ 0, 1, 1 /)
  pdata%poledges( 4,1:pdata%nvpols(4)-1)  = (/ 0, 1, 1, 1, 0, 0, 1 /)
  pdata%poledges( 5,1:pdata%nvpols(5)-1)  = (/ 1, 1, 0, 1 /)
  pdata%poledges( 6,1:pdata%nvpols(6)-1)  = (/ 0, 1, 1 /)
  pdata%poledges( 7,1:pdata%nvpols(7)-1)  = (/ 1, 0, 1 /)
  pdata%poledges( 8,1:pdata%nvpols(8)-1)  = (/ 0, 0, 0, 1 /)
  pdata%poledges( 9,1:pdata%nvpols(9)-1)  = (/ 1, 0, 1, 1, 1 /)
  pdata%poledges(10,1:pdata%nvpols(10)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(11,1:pdata%nvpols(11)-1) = (/ 0, 0, 1 /)
  pdata%poledges(12,1:pdata%nvpols(12)-1) = (/ 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1 /)
  pdata%poledges(13,1:pdata%nvpols(13)-1) = (/ 1, 0, 1 /)
  pdata%poledges(14,1:pdata%nvpols(14)-1) = (/ 0, 0, 0, 1 /)
  pdata%poledges(15,1:pdata%nvpols(15)-1) = (/ 1, 0, 1, 1 /)
  pdata%poledges(16,1:pdata%nvpols(16)-1) = (/ 1, 0, 1 /)
  pdata%poledges(17,1:pdata%nvpols(17)-1) = (/ 0, 0, 1, 0, 1 /)
  pdata%poledges(18,1:pdata%nvpols(18)-1) = (/ 0, 0, 0, 1, 0, 1 /)
  pdata%poledges(19,1:pdata%nvpols(19)-1) = (/ 0, 1, 0, 1, 1 /)
  pdata%poledges(20,1:pdata%nvpols(20)-1) = (/ 0, 0, 0, 0, 1 /)
  pdata%poledges(21,1:pdata%nvpols(21)-1) = (/ 0, 1, 0, 1 /)
  pdata%poledges(22,1:pdata%nvpols(22)-1) = (/ 0, 0, 0, 1 /)
  pdata%poledges(23,1:pdata%nvpols(23)-1) = (/ 0, 1, 0, 0, 1 /)
  pdata%poledges(24,1:pdata%nvpols(24)-1) = (/ 0, 0, 1, 0, 1 /)
  pdata%poledges(25,1:pdata%nvpols(25)-1) = (/ 0, 0, 0, 0 /)
  pdata%poledges(26,1:pdata%nvpols(26)-1) = (/ 0, 0, 0, 1, 0, 0, 0, 1 /)
  pdata%poledges(27,1:pdata%nvpols(27)-1) = (/ 0, 0, 0, 0 /)
  pdata%poledges(28,1:pdata%nvpols(28)-1) = (/ 0, 0, 1, 1, 0, 1 /)
  pdata%poledges(29,1:pdata%nvpols(29)-1) = (/ 0, 0, 1, 0, 0, 1 /)
  pdata%poledges(30,1:pdata%nvpols(30)-1) = (/ 0, 0, 1 /)
  pdata%poledges(31,1:pdata%nvpols(31)-1) = (/ 0, 0, 0, 0, 0, 1 /)
  pdata%poledges(32,1:pdata%nvpols(32)-1) = (/ 0, 1, 1 /)
  pdata%poledges(33,1:pdata%nvpols(33)-1) = (/ 0, 0, 1, 1 /)
  pdata%poledges(34,1:pdata%nvpols(34)-1) = (/ 0, 0, 0, 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

  ! Code below is to show that the area of "Tierra del Fuego" (polygon
  ! A_{27}) is very small; and, if it is much smaller than the desired
  ! area, then every solution in which a cell covers "Tierra del
  ! Fuego" only is such that J is at least ( vol("Tierra del Fuego") -
  ! ideal cell's vol )^2 / nsites.
!!$  write(*,*) 'With the computed scaling factor, scale = ',scale
!!$  volt = 0.0d0
!!$  do j = 1,pdata%npols
!!$     volpj = 0.0d0
!!$     do k = 1,pdata%nvpols(j) - 1
!!$        volpj = volpj + pdata%xpol(j,k) * pdata%ypol(j,k+1) - pdata%xpol(j,k+1) * pdata%ypol(j,k)
!!$     end do
!!$     write(*,*) 'Volume of polygon ',j,' is equal to: ',0.5d0 * volpj
!!$     volt = volt + volpj
!!$  end do
!!$  write(*,*) 'And the ideal volume of a cell is ',volt / nsites
     
  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 = 0.35d0 / scale

  ! 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,"(A16,I0,A3)") 'fig-america-ini-',nsites,'.mp'
        call drawsol(n,x0,filename,c_loc(pdata))
  
        write(filename,"(A16,I0,A3)") 'fig-america-end-',nsites,'.mp'
        call drawsol(n,x,filename,c_loc(pdata))

        call compvols(n,x,vol,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, &
             pdata%volA / nsites, &
             count( vol(1:nsites) .lt. sqrt( epsilon( 1.0d0 ) ) ), &
             maxval( abs( vol(1:nsites) * nsites / pdata%volA - 1.0d0 ), vol(1:nsites) .ge. sqrt( epsilon( 1.0d0 ) ) )
        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,1X,1P,E12.5,1X,I6,1X,1P,E12.5)
  
end program optvorma
