C     =================================================================
C     File: grasp.f
C     =================================================================

      program grasp

C     This program try to find some y such that f(y) <= f(x) by
C     using the GRASP algorithm, described in:
C
C     C.N. MENESES, P.M. PARDALOS, AND M.G.C. RESENDE. GRASP FOR
C     NONLINEAR OPTIMIZATION. AT&T Labs Research Technical
C     Report TD-6DUTRG.

C     PARAMETERS
      integer MaxNumIterNoImprov, MaxIter, MaxNoFunctionImprov,nmax,mmax
      double precision INFINITO,alpha,ZERO
      parameter ( MaxNumIterNoImprov  = 20         )
      parameter ( MaxIter             = 2000       )
      parameter ( INFINITO = 10000000000000000.0d0 )
      parameter ( ZERO     = 0.000000000001d0      )
      parameter ( alpha               = 0.3d0      )
      parameter ( MaxNoFunctionImprov = 100        )
      parameter ( nmax    = 1000 )
      parameter ( mmax    =    1 )

C     LOCAL SCALARS
      real etime,ntime,btime
      integer n,m,i,iter,NumIterNoImprov,err,itbest,NoFunctionImprov,
     + nsamples,bnsamples,totevalf,bevalf,fcnmb
      double precision seed,drand,h,BestSolValue,SolValue,vdummy

C     LOCAL ARRAYS
      real dum(2)
      double precision BestSolution(nmax),xlocalopt(nmax),y(nmax),
     + l(nmax),u(nmax),
     + lambda(mmax),rho(mmax),equatn(mmax),linear(mmax)

C     SET LOCAL SCALARS

      h            = 1.0d0
      SolValue     = INFINITO
      BestSolValue = INFINITO

C     SET LOCAL COUNTERS

      totevalf    = 0
      nsamples = 0

C     SETUP PROBLEM

      call inip(n,y,l,u,m,lambda,rho,equatn,linear)

C     SET INITIAL POINT

      do i = 1,n
          y(i) = l(i) + ( u(i) - l(i) ) * 0.5d0
      end do
      nsamples = nsamples + 1

C     START METHOD

C      do while ( ( iter .lt. MaxIter ) .and. 
C     +           ( f .le. BestSolValue ) )

      iter         = 1
      NoFunctionImprov = 0
      
      open(11,file='fbest.txt')

      do while ( iter .le. MaxIter )

          call ConstructGreedyRandomizedSolution(y,n,h,l,u,seed,alpha,
     +         fcnmb)
          totevalf = totevalf + fcnmb
C          call tun_lissajous(n,y,m,lambda,rho,l,u,SolValue,err)

          if ( n .le. 4 )  then 


              call LocalSearch(y,h,n,l,u,xlocalopt,fcnmb,err)
C              call evalal(n,xlocalopt,m,SolValue,err)
              call evalf(n,xlocalopt,SolValue,err)
	      totevalf = totevalf + fcnmb + 1

C              call tun_lissajous(n,y,m,lambda,rho,l,u,err)
C              call evalal(n,y,m,SolValue,err)

C              if ( SolValue .lt. BestSolValue ) then
	
C                  do i = 1,n
C                      y(i) = xlocalopt(i)
C                  end do

C              end if

          else

C             call EvaluateFunction(x,SolValue)
C              call evalal(n,y,m,SolValue,err)
              call evalf(n,y,SolValue,err)
	      totevalf = totevalf + 1

          end if
 
C         CHECK IF ITS BETTER THAN f(x)

          ntime = etime(dum) 

	  if( BestSolValue .gt. ( SolValue + 0.0001 * abs( SolValue ) 
     +        + 1.0E-6 ) ) then	  

C             UPDATES FUNCTION IMPROVEMENT ONLY FOR RELEVANT VALUES

              NoFunctionImprov = 0
              NumIterNoImprov  = 0

C             UPDATES THE BEST SOLUTION

              do i = 1,n
                  BestSolution(i) = y(i)
              end do

              itbest    = iter
              btime     = ntime
              bevalf    = totevalf
	      bnsamples = nsamples

          write(*,*) "--------------- ",iter," -------------------"
          write(*,*) " "
          write(*,*) "h,inoi,nfci | F,bt | x0,totf > ", h,
     +               NumIterNoImprov,NoFunctionImprov,' | ',
     +               SolValue,btime,' | ',bnsamples,bevalf

              BestSolValue = SolValue

          else

              NumIterNoImprov  = NumIterNoImprov + 1
              NoFunctionImprov = NoFunctionImprov + 1

          end if

C         SET A NEW SAMPLE POINT

          if ( NoFunctionImprov .ge. MaxNoFunctionImprov ) then

              seed = iter+1
              vdummy = drand(seed)

              do i = 1,n
                  y(i) = l(i) + ( u(i) - l(i) ) * drand(seed)
              end do
	      nsamples = nsamples + 1

C             RETURNS TO ORIGINAL GRID

              h = 1.0d0
              NumIterNoImprov = 0
	      NoFunctionImprov = 0

          else

C             MAKES THE GRID MORE DENSE

              if( NumIterNoImprov .ge. MaxNumIterNoImprov ) then

                  h = h / 2.0d0
                  NumIterNoImprov = 0

              end if

          end if
	  
	  write(11,0011) iter,SolValue,BestSolValue,ntime

C         ITERATE

          iter = iter + 1

      end do
      
      close(11)

      write(*,*) "............... FINAL ..................."
      write(*,*) " "
      write(*,*) 'f(*) : ',BestSolValue,' | it(*) : ',itbest 
      write(*,*) 'bt : ',btime,'( ',ntime,' ) '
      write(*,*) 'evalf : ',bevalf,' | tot(x0) : ',bnsamples
      do i = 1,n
          write(*,*) 'x (',i,') = ',BestSolution(i)
      end do
      write(*,*) "........................................"
      do i = 1,n
          write(*,*) 'x0 (',i,') = ',l(i) + ( u(i) - l(i) ) * 0.5d0
      end do

C     NON-EXECUTABLE STATEMENTS
 0010 format('f = ',D25.16)
 0011 format(I7,1X,D25.16,1X,D25.16,1X,F8.2)

      end;

C     ******************************************************************
C     ******************************************************************

      subroutine ConstructGreedyRandomizedSolution(x,n,h,l,u,seed,
     +                                             alpha,totevalf)
C     This subroutine try to find some y such that f(y) <= f(x) by
C     using the GRASP algorithm, described in:
C
C     C.N. MENESES, P.M. PARDALOS, AND M.G.C. RESENDE. GRASP FOR
C     NONLINEAR OPTIMIZATION. AT&T Labs Research Technical
C     Report TD-6DUTRG.

C     PARAMETERS
      double precision INFINITO
      parameter ( INFINITO = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,totevalf,fcnmb
      double precision seed,alpha,h

C     ARRAY ARGUMENTS
      double precision x(n),l(n),u(n)

C     LOCAL SCALARS
      integer i,j,k,RandomlySelectElement
      double precision vmax,vmin,xibest

C     LOCAL ARRAYS
      integer S(n+1),RCL(n+1)
      double precision z(n),g(n),y(n)

C     SET LOCAL SCALARS

      totevalf = 0
      fcnmb = 0

      do i = 1,n

          S(i)   = i
          RCL(i) = -1
          y(i)   = x(i)

      end do

      S(n+1)   = n
      RCL(n+1) = 0

C     TRY TO FIND SOME POINT z SUCH THAT f(x) < f(x)

      do while( S(n+1) .gt. 0 )

          vmin = INFINITO
          vmax = ( -1.0d0 ) * INFINITO

          do i = 1,S(n+1)

C             MINIMIZE COMPONENT y( S(i) ) WITH STEP h

C                  Minimize(x,h,indx,n,l,u,xibest,fbest,m)
              call Minimize(y,h,S(i),n,l,u,xibest,g(S(i)),fcnmb)
	      totevalf = totevalf + fcnmb

              z(S(i)) = xibest

              if( vmin .gt. g(S(i)) ) then

                  vmin = g(S(i))

              end if

              if( vmax .lt. g(S(i)) ) then

                  vmax = g(S(i))

              end if

          end do

C         SET SOLUTION CANDIDATES

          RCL(n+1) = 0

          do i=1,S(n+1)

              if( g( S(i) ) .le. ( vmin + alpha * ( vmax - vmin ) ) )
     +        then

                  RCL(n+1) = RCL(n+1) + 1
                  RCL( RCL(n+1) ) = i

              end if

          end do

C         RANDOMLY SELECT A RCL ELEMENT

          seed = ( S(n+1) + seed ) / 2.0d0
          j = RandomlySelectElement(seed,1,RCL(n+1))

C         UPDATES x,y

          x( RCL(j) ) = z( RCL(j) )
          y( RCL(j) ) = z( RCL(j) )

C         S = S \ {j}

          k           = S( RCL(j) )
          S( RCL(j) ) = S( S(n+1) )
          S( S(n+1) ) = k
          S(n+1)      = S(n+1) - 1

      end do

      end

C     ******************************************************************
C     ******************************************************************

      integer function RandomlySelectElement(seed,l,u)
C     This subroutine returns a random integer between l and u

C     Parameters of the subroutine:
C
C     On Entry:
C
C     seed     double precision,
C
C     l        integer,
C              intervals lower bound,
C
C     u        integer,
C              intervals upper bound,
C
C     SCALAR ARGUMENTS
      integer l,u
      double precision seed

C     LOCAL SCALARS
      integer delta
      double precision drand,nran

      delta = u - l + 1
      nran  = drand(seed)
      nran  = drand(seed)

      RandomlySelectElement = max(nint( nran * delta ),1)

      return

      end;

C     ******************************************************************
C     ******************************************************************

      subroutine Minimize(x,h,indx,n,l,u,xibest,fbest,totevalf)
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     x        double precision x(n),
C              initial point,
C
C     h        double precision,
C              grid density
C
C     indx     integer,
C              index to be moved,
C
C     n        integer,
C              number of variables,
C
C     l        double precision l(n),
C              lower bounds on x,
C
C     u        double precision u(n),
C              upper bounds on x,
C
C     On Return
C
C     xibest   double precision,
C              best position found
C
C     fbest    double precision
C              best f's value (with xibest)

C     PARAMETERS
      double precision INFINITO
      parameter ( INFINITO = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,totevalf,indx
      double precision h,xibest,fbest

C     ARRAY ARGUMENTS
      double precision x(n),l(n),u(n)

C     LOCAL SCALARS
      integer i,err
      double precision vmin,vle

C     LOCAL ARRAYS
      double precision t(n)

C     SET LOCAL SCALARS

      totevalf = 0

      do i = 1,n

          t(i) = x(i)

      end do

      t(indx) = l(indx)

      vmin = INFINITO
C      write(*,*) 'xa -----> ',t(1),t(2),t(3),t(4)

      do while ( t(indx) .le. u(indx) )

C          call EvaluateFunction(t,vle)
C          call evalal(n,t,m,vle,err)
          call evalf(n,t,vle,err)
	  totevalf = totevalf + 1

          if( vmin .gt. vle ) then

              vmin = vle
              xibest = t(indx)

          end if

C          write(*,*) 'h,indx,t,xibest,vle,vmin > ', h,indx,t(indx),
C     +     xibest,vle,vmin
          t(indx) = t(indx) + h

      end do

      fbest = vmin
C      write(*,*) 'xd -----> ',t(1),t(2),t(3),t(4)
C      write(*,*) 'FB -> ',fbest
C      read(*,*) i

      end

C     ******************************************************************
C     ******************************************************************

      subroutine LocalSearch(x,h,n,l,u,xlocalopt,
     +                       totevalf,err)

C     PARAMETERS
      double precision INFINITO
      parameter ( INFINITO = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,totevalf,err
      double precision h

C     ARRAY ARGUMENTS
      double precision x(n),xlocalopt(n),
     +                 l(n),u(n)

C     LOCAL SCALARS
      integer i,maxldir,ldir
      double precision xlocaloptValue,xprimeValue
      logical Feasible,Improved

C     ARRAY SCALARS
      double precision all_dir(324),dir(n),xprime(n)
      
      totevalf = 0

      if ( n .gt. 4 ) then

          err = 1

      else

C         SET LOCAL SCALARS
          Improved = .true.

          do i = 1,n

              xlocalopt(i) = x(i)

          end do
          call GenerateDir(n,all_dir,maxldir,err)
          ldir = 0

C          call evalal(n,xlocalopt,m,xlocaloptValue,err)
          call evalf(n,xlocalopt,xlocaloptValue,err)
	  totevalf = totevalf + 1

          do while ( Improved )

              Improved = .false.

              do while ( ( ldir .lt. maxldir ) .and.
     +                   ( .not. ( Improved ) ) )

C                 SET NEW DIRECTION
                  do i=1,n

                      dir(i) = all_dir( ldir + i )

                  end do

                  ldir = ldir + n

                  do i=1,n

                      xprime(i) = xlocalopt(i) + h * dir(i)

                  end do

C                  call evalal(n,xprime,m,xprimeValue,err)
                  call evalf(n,xprime,xprimeValue,err)
		  totevalf = totevalf + 1

                  if ( ( Feasible(xprime,n,l,u) ) .and.
     +                 ( xprimeValue .lt. xlocaloptValue ) ) then

                      do i=1,n

                          xlocalopt(i) = xprime(i)

                      end do

                      xlocaloptValue = xprimeValue
                      Improved = .true.
                      ldir      = 0

                  end if

              end do

          end do

      end if

      end


C     ******************************************************************
C     ******************************************************************

      subroutine GenerateDir(n,dir,ldir,err)
C     This subroutine generates a set of directions indicated by dir.
C     For n < 5 only.


C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              problem dimension
C
C     dir      integer,
C              intervals lower bound,
C
C     On Return:
C
C     dir      double precision dir(321),
C              directions
C
C     err      integer,
C              return status: -1 error, else ok
C
C     SCALAR ARGUMENTS
      integer n,ldir,err

C     ARRAY ARGUMENTS
      double precision dir(324)

C     LOCAL SCALARS
      integer i,j,ilocal,tset,addone


C     LOCAL ARRAYS
      double precision next(4)

C     SET dir index
      tset = ( 3**n )
      ldir     =  tset * n

C     SET LOCAL SCALARS
      next(1) = -1.0d0
      next(2) = -1.0d0
      next(3) = -1.0d0
      next(4) = -1.0d0

      err = -1

C     VERIFY IF n < 5

      if ( n .lt. 5 ) then

          do i = 1,ldir

              ilocal = mod(i,n)

C             Verify if the series finished
              if ( ilocal .eq. 0 ) then

                  dir(i) = next(n)

                  addone = 1.0d0

                  do j = 1,n

                      next(j) = next(j) + addone

                      if( next(j) .gt. 1.0d0 ) then

                          next(j) = -1.0d0

                      else	

                          addone = 0.0d0

                      end if

                  end do

              else

                  dir(i) = next(ilocal)

              end if

          end do

          err = 0

      end if

      end;

C     ******************************************************************
C     ******************************************************************

      logical function Feasible(x,n,l,u)

C     PARAMETERS

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n),l(n),u(n)

C     LOCAL SCALARS
      integer i
      logical feas

      feas = .true.

      i = 1

      do while ( ( i .le. n ) .and. ( feas ) )

          if ( ( x(i) .lt. l(i) ) .or. ( x(i) .gt. u(i) ) ) then

              feas = .false.

          end if

          i = i + 1

      end do

      Feasible = feas

      return

      end

C     ******************************************************************
C     ******************************************************************

      double precision function drand(ix)
      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return

      end