C     =================================================================
C     File: tunn-main.f
C     =================================================================

      subroutine tun_lissajous(n,x,m,lambda,rho,l,u,retcode)

C     This subroutine try to find some y such that f(y) < f(x) by
C     using the Lissajous's curve defined by:
C
C     lis(t) = (cos( p1 * r + t1 ), ...  cos( pn * r + tn ) )
C
C     This algorithm was published on the paper:
C     "OTIMIZACAO GLOBAL USANDO TRAJETORIAS DENSAS E APLICACOES"
C     Mario Salvatierra Junior, UNICAMP
C

C     PARAMETERS
      integer ntrials
      double precision INFINITO
      parameter ( ntrials   = 100000 )
      parameter ( INFINITO  = 10000000000000000.0d0 )

C     SCALAR ARGUMENTS
      integer n,m,retcode

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

C     LOCAL SCALARS
      integer i,j,k,err,itrial
      double precision ind,num,den,signal,flis,r,delta,lim,fc,bl

C     LOCAL ARRAYS
      double precision t(n),p(n),lis(n)

C     Set the local scalars

      delta   = -1.0d0
      lim     = -0.0001d0
      retcode = 0
      bl      = INFINITO

C     SET INITIAL POINT

      do i = 1,n
          lis(i) = x(i)
      end do

C     SET f(x)

      flis    = INFINITO
C      call evalal(n,x,m,fc,err)
      call evalf(n,x,fc,err)

C     Set the curve's initial position

      do i = 1,n

          t( i ) = acos( ( 2.0d0 * x( i ) - l( i ) - u( i ) ) /
     +                 ( u( i ) - l( i ) ) )

      end do

C     Set the p's value in radians
      open(01,file="primes.txt")

      do i = 1,n

          read(01,*) k
          p( i ) = k
          p( i ) = sqrt( p( i ) )

      end do

      close(01)

C      open(01,file="l_1_2.txt" )
C      open(02,file="l_1_3.txt" )
C      open(03,file="l_1_4.txt" )
C      open(04,file="l_2_3.txt" )
C      open(05,file="l_2_4.txt" )
C      open(06,file="l_3_4.txt" )

C     TRY THE BEST SOLUTION MOVING ONE POINT EACH TIME

C      do j = 1,n

C         SET THE LOCAL SCALARS METHOD

          ind     = 1.0d0
          den     = ind
          num     = 1.0d0
          r       = 0.0d0
          signal  = 1.0d0
          itrial  = 0

C         ITERATE THE LISSAJOUS'S TUNNELING METHOD

          do while ( itrial .le. ntrials ) 
C     +               ( flis .gt. fc ) )

              itrial = itrial + 1

C             Set next r, where r is on the form:
C             r = 1, -1, 2/1, 1/2, -2/1, -1/2, ...
              r = signal * ( num / den )

              den = den - 1.0d0
              num = num + 1.0d0

              if ( den .le. 0.0d0 ) then

                  if ( signal .lt. 0.0d0 ) then

                      signal = 1.0d0
                      ind = ind * 2.0d0

                  else

                      signal = -1.0d0

                  end if

                  den = ind
                  num = 1.0d0

              end if

C             MOVE ALONG THE CURVE

              do j = 1,n

C                 Set a point lis in the Lissajous's curve

                  lis( j ) = 0.5d0 * ( l( j ) + u( j )
     +                     + ( u( j ) - l( j ) ) *
     +                     ( cos( p( j ) * r + t( j ) ) ) )

              end do

C             Set f on the lis, point:

C              write(01,0010) itrial,lis(1),lis(2)
C              write(02,0010) itrial,lis(1),lis(3)
C              write(03,0010) itrial,lis(1),lis(4)
C              write(04,0010) itrial,lis(2),lis(3)
C              write(05,0010) itrial,lis(2),lis(4)
C              write(06,0010) itrial,lis(3),lis(4)
C              call evalf(n,lis,flis,err)
C              call evalal(n,lis,m,flis,err)
              call evalf(n,lis,flis,err)

C             DEBUG
              if (bl .gt. flis ) then

                  bl = flis
                  write(*,*) 'BLIS -> [',itrial,']',bl,'/',fc

              end if

             if ( fc .gt. flis ) then

                  fc = flis
                  do j =1,n
                      x(j) = lis(j)
                  end do

              end if

          end do
C          close(01)
C          close(02)
C          close(03)
C          close(04)
C          close(05)
C          close(06)

C         For the next iteration choose the best point found

C          lis( j ) = x( j )

C      end do 
       read(*,*) itrial

C     NON-EXECUTABLE STATEMENTS
 0010 format(i7,D20.12,1X,D20.12)

      end

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

      subroutine tun_grasp(n,x,m,l,u,f,itrial,g,gpsupn,inform)

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
      integer MaxNumIterNoImprov, MaxIter, MaxNoFunctionImprov
      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        )

C     SCALAR ARGUMENTS
      integer n,m,itrial
      double precision f

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

C     SCALAR ARGUMENTS
      integer inform
      double precision gpsupn

C     ARRAY ARGUMENTS
      double precision g(n)

C     LOCAL SCALARS
      real etime,ntime,btime
      integer i,iter,NumIterNoImprov,err,itbest,NoFunctionImprov
      double precision seed,h,BestSolValue,SolValue

C     LOCAL ARRAYS
      real dum(2)
      double precision BestSolution(n),y(n),xlocalopt(n)

C     SET LOCAL SCALARS

      iter         = itrial
      h            = 1.0d0
      SolValue     = INFINITO
      BestSolValue = f
      seed         = 0.0d0
      NumIterNoImprov = 0

C     SET INITIAL POINT

      do i = 1,n

          y(i) = x(i)
          BestSolution = x(i)

      end do

C     START METHOD

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

      NoFunctionImprov = 0

      do while ( ( iter .lt. MaxIter ) .and.
     +           ( NoFunctionImprov .le. MaxNoFunctionImprov ) )
C         Set new iteration

          iter = iter + 1

          call ConstructGreedyRandomizedSolution(y,n,h,l,u,seed,alpha,m)
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,
     +                         m,err)
C              call evalal(n,xlocalopt,m,SolValue,err)
               call evalf(n,xlocalopt,SolValue,err)

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)

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

          ntime = etime(dum) 

          if( SolValue .lt. BestSolValue ) then

C             UPDATES FUNCTION IMPROVEMENT ONLY FOR RELEVANT VALUES

              if( abs( SolValue - BestSolValue ) .gt. 1.0E-5 ) then
                  NoFunctionImprov = 0
              end if

C             UPDATES THE BEST SOLUTION

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

              BestSolValue = SolValue
              itbest       = iter
              NumIterNoImprov = 0

              btime = ntime

          else

              NumIterNoImprov  = NumIterNoImprov + 1
              NoFunctionImprov = NoFunctionImprov + 1

          end if

          write(*,*) "=============== ",iter," ==================="
          write(*,*) " "
          write(*,*) "h,i,i2 | f,t | F,bt > ", h,NumIterNoImprov,
     +               NoFunctionImprov,' | ',
     +               SolValue,ntime,' | ',BestSolValue,btime

C         MAKES THE GRID MORE DENSE

          if( NumIterNoImprov .ge. MaxNumIterNoImprov ) then

              h = h / 2.0d0
              NumIterNoImprov = 0

          end if

      end do

      itrial = iter

      if ( f .gt. BestSolValue) then

          f     = BestSolValue

          do i = 1,n

              x(i) = BestSolution(i)

          end do

      end if

      write(*,*) "=============== FINAL ==================="
      write(*,*) " "
      write(*,*) "i,f*,bt,tt ",itbest,BestSolValue,btime,ntime
      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

      read(*,*)

      end;

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

      subroutine ConstructGreedyRandomizedSolution(x,n,h,l,u,seed,
     +                                             alpha,m)
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,m
      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

      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)),m)

              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,nran

C     LOCAL SCALARS
      integer delta
      double precision drand

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

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

      return

      end;

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

      subroutine Minimize(x,h,indx,n,l,u,xibest,fbest,m)
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,m,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

      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)

          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,
     +                       m,err)

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

C     SCALAR ARGUMENTS
      integer n,m,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)

      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)

          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)

                  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