module MINLP

  use Constants
  use BranchAndBound_Options

  logical, parameter :: RANDOM_DIRECTION = .true.

contains

  ! ************************************************************************
  ! ************************************************************************

  subroutine BranchAndBound(n,incumbent,lower,upper,m,lambda,equatn,linear,   &
       rho,constf,epsfeas,epsopk,maxit,optSolverIter,objectiveValue,g,gpsupn, &
       cnorm,cnormu,geninfo,inform,numBinary,numInteger,iteration,            &
       solutionFound,nodeSelection,branchingVariableStrategy,                 &
       multistartStrategy,maxSolverCalls,numNodesGenerated,numProblemsSolved, &
       iterationSolutionFound,variableTypes,minlpinfo,options)

    use Auxiliar
    use LinkedList
    use OptSolver
    use Problems
    use ProblemType
    use PseudoCosts
    use Random
    use Structures
    use VariableSelection

    implicit none

    ! SCALAR ARGUMENTS
    logical,         intent(in)     :: constf
    integer,         intent(in)     :: m,n,maxit,numBinary,numInteger
    integer,         intent(in)     :: maxSolverCalls
    integer,         intent(in)     :: nodeSelection
    integer,         intent(in)     :: branchingVariableStrategy,multistartStrategy
    integer,         intent(out)    :: iteration,geninfo,optSolverIter,minlpinfo
    integer,         intent(in out) :: inform
    real(kind=8),    intent(out)    :: objectiveValue,gpsupn
    real(kind=8),    intent(in)     :: epsfeas,epsopk
    real(kind=8),    intent(in out) :: cnorm,cnormu
    logical,         intent(out)    :: solutionFound
    integer,         intent(out)    :: iterationSolutionFound
    integer,         intent(out)    :: numNodesGenerated,numProblemsSolved

    ! ARRAY ARGUMENTS
    logical,      intent(in)     :: equatn(m),linear(m)
    integer,      intent(in)     :: variableTypes(n)
    real(kind=8), intent(in)     :: lambda(m),rho(m)
    real(kind=8), intent(in out) :: lower(n),upper(n),incumbent(n)
    real(kind=8), intent(out)    :: g(n)

    type(bb_options), intent(in) :: options

    ! Parameters of the subroutine:
    ! =============================
    !
    ! On Entry:
    ! =========
    !
    ! INCUMBENT real(kind=8) incumbent(n): Initial estimation of the solution.
    ! -----------------------------------
    !
    ! N integer: Number of variables.
    ! ---------
    !
    ! LOWER real(kind=8) lower(n): Lower bounds of the variables.
    ! ---------------------------
    !
    ! UPPER real(kind=8) upper(n): Upper bounds of the variables.
    ! ---------------------------
    !
    ! M integer: Number of constraints (equalities plus inequalities)
    ! ---------  without considering the bound constraints.
    !
    ! LAMBDA real(kind=8) lambda(m): 
    ! -----------------------------
    !
    ! Estimation of the Lagrange multipliers.
    !
    ! EQUATN logical equatn(m):
    ! ------------------------
    !
    ! Logical array that, for each constraint, indicates whether the 
    ! constraint is an equality constraint (TRUE) or an inequality 
    ! constraint (FALSE).
    !
    ! LINEAR logical linear(m)
    ! ------------------------
    !
    ! Logical array that, for each constraint, indicates whether the 
    ! constraint is a linear constraint (TRUE) or a nonlinear constraint 
    ! (FALSE).
    !
    ! NUMBINARY integer: Number of binary variables.
    ! -----------------
    !
    ! NUMINTEGER integer: Number of integer variables.
    ! ------------------
    !
    ! NODESELECTION integer:
    ! ---------------------
    !
    ! The rule for selecting the next node.
    !
    ! BRANCHINGVARIABLESTRATEGY integer: 
    ! ---------------------------------
    !
    ! The strategy to be used to select the branching variable. The
    ! possible strategies are:
    !
    ! 1) Choose the non-integer variable whose index is smaller
    ! (LOWEST_INDEX_FIRST).
    !
    ! 2) Choose the non-integer variable which is most fractional,
    ! that is, the variable whose value is farther than an integer
    ! (MOST_FRACTIONAL).
    !
    ! MULTISTARTSTRATEGY integer:
    ! --------------------------
    !
    ! The strategy to be used joint with multistart. It can be:
    !
    ! 1) Among all feasible solutions found for a subproblem, choose
    ! the best one (BEST_FEASIBLE).
    !
    ! 2) Solve a subproblem at most 'maxSolverCalls' times and stop
    ! with the first feasible solution (FIRST_FEASIBLE).
    !
    ! If no feasible solution was found for a subproblem, this
    ! subproblem will be declared infeasible.
    !
    ! MAXSOLVERCALLS integer: Maximum number of attempts to solve the
    ! ----------------------  same subproblem.
    !
    ! On Return:
    ! ==========
    !
    ! INCUMBENT real(kind=8) incumbent(n): The solution found for the problem.
    ! -----------------------------------
    !
    ! OBJECTIVEVALUE real(kind=8): Value of the solution found.
    ! ---------------------------
    !
    ! ITERATION integer: Number of iterations of the algorithm.
    ! -----------------
    !
    ! SOLUTIONFOUND logical: Indicates whether a solution was
    ! ---------------------  found (TRUE) or not (FALSE).
    !
    ! ITERATIONSOLUTIONFOUND integer: Iteration at which the solution was
    ! ------------------------------  found.
    !
    ! NUMNODESGENERATED integer: Number of generated nodes.
    ! -------------------------
    !
    ! NUMPROBLEMSSOLVED integer: Number of (sub)problems solved.
    ! -------------------------

    ! LOCAL SCALARS
    integer :: branchingIndex,m_
    logical :: satisfy_domain,feasibleSolutionFound
    real(kind=8) :: f,best_f,round_f,cutoff,boundTolerance
    integer :: listType

    real(kind=8), save :: seed = 101063.0d0
    real(kind=8), save :: seed_prob = 7919.0d0

    integer :: j,k,nextID,jj
    integer :: i
    integer :: multiStartIteration
    integer :: variableType
    integer :: maxSolverCallsInc,updatesBeforeTrust

    integer :: iter_,iter_2

    real(kind=8) :: f_,f_2,lowerBound,frac,S0

    ! Counters
    integer :: maxSize,pruningByBoundCount,pruningByStoredBound,&
               pruningByInfeas,pruningByOpt

    integer, save :: totalNumNodesGenerated = 0
    integer, save :: totalNumProblemsSolved = 0

    ! LOCAL ARRAYS
    real(kind=8) :: l(n),u(n),x(n),best_x(n),round_x(n),l_(n),u_(n)

    logical :: switched,circlePackingInitialPoint,solveRoundedSolution

    real :: time,startTime

    ! OTHER LOCAL VARIABLES
    type(data)    :: data_
    type(domain)  :: mainDomain,subdomain
    type(list)    :: openProblemsList
    type(problem) :: p,p1,p2

    type(set), allocatable :: discrete_sets(:)

    ! **************
    ! INITIALIZATION
    ! **************

    boundTolerance = getBoundTolerance(options)
    solveRoundedSolution = getSolveRoundedSolution(options)

    circlePackingInitialPoint = .false.

    call cpu_time(startTime)

    cutoff = getCutOff(options)

    select case(nodeSelection)
    case(DEPTH_FIRST)
       listType = STACK
    case(SWITCH_DF_BBR)
       listType = STACK
    case(BREADTH_FIRST)
       listType = QUEUE
    case(BEST_BOUND_RULE)
       listType = PRIORITY
    end select

    switched = .false.
    lowerBound = LARGEST_VALUE
    optSolverIter = 0
    iteration = 0
    minlpinfo = 0
    nextID    = 2
    solutionFound = .false.
    iterationSolutionFound = -1

    maxSolverCallsInc = getNumNodeResolutionsIncrease(options)

    updatesBeforeTrust = getNumberOfUpdatesBeforeTrust(options)

    ! Counters.
    maxSize = 1
    pruningByBoundCount = 0
    pruningByStoredBound = 0
    pruningByInfeas = 0
    pruningByOpt    = 0

    numNodesGenerated = 0
    numProblemsSolved = 0

    ! Initialize the variables pseudo-costs.
    if(branchingVariableStrategy .eq. PSEUDO_COSTS) then
       call initializePseudoCosts(n)
    end if

    call initialize(n,numBinary,numInteger,mainDomain,discrete_sets,&
         variableTypes)

    ! Create the main problem.
    ! An empty array is given by (/ (0, i=1,0) /).
    p = createProblem(n,incumbent,createDomain(0,(/ (0, jj=1,0) /),&
         (/ (BINARY, jj=1,0) /),(/ (0.0d0, jj=1,0) /) ),1,0,cutoff,&
         -1,0.0d0,.false.,0)

    ! Create the list of open problems and add the main problem to it.
    data_%problem = p
    openProblemsList = createList(listType)
    call add(openProblemsList,data_)

    ! Verify whether the initial point is a feasible solution.
    if(satisfyDomain(n,incumbent,mainDomain,discrete_sets) .and. &
         satisfyConstraints(n,incumbent,lower,upper)) then

       ! The initial point is an integer feasible solution.

       ! Evaluate the objective function.
       call sevalal(n,incumbent,m,lambda,rho,equatn,linear,f,inform)

       if(inform .lt. 0) then

          call deallocateDomain(mainDomain)
          call destroyList(openProblemsList)
          deallocate(discrete_sets)
          if(branchingVariableStrategy .eq. PSEUDO_COSTS) then
             call deallocatePseudoCosts()
          end if

          return

       end if

       if(isNumber(f)) then

          ! Update the incumbent solution.
          objectiveValue = f
          solutionFound  = .true.
          iterationSolutionFound = 0

       end if

    else

       ! Initialize the objective value of the best solution with
       ! "infinity".
       objectiveValue = LARGEST_VALUE

    end if

    ! *********
    ! MAIN LOOP
    ! *********
    !
    ! The problems are removed from the open problems list until it
    ! becomes empty.
    do while(.not. empty(openProblemsList))

       call flush()

       if(nodeSelection .eq. SWITCH_DF_BBR .and. & 
          iteration .eq. getNumIterationsForBestBoundRule(options)) then

          call changeListType(openProblemsList,PRIORITY)

       end if

       numNodesGenerated = numNodesGenerated + 1
       iteration = iteration + 1

       ! Get a subprolem.
       data_ = remove(openProblemsList)
       p = data_%problem

       ! Add the bounds on the variables.
       l = lower
       u = upper
       call addBounds(p%domain,l,u,n)

       ! Get the initial point for this subproblem.
       x = p%initialPoint

       ! PRUNING BY (STORED) BOUND
       ! -------------------------
       !
       ! Check whether the stored lower bound for this subproblem is
       ! greater than or equal to the value of the incumbent
       ! solution. If it occurs, the subproblem can be discarded,
       ! since it cannot provide a better solution.

       if(solutionFound .and. p%geninfoFather .eq. 0 .and. &
          boundTolerance * p%lowerBound .ge. objectiveValue) then

          pruningByStoredBound = pruningByStoredBound + 1
          call deallocateProblem(p)

          ! Continue with the next subproblem from the list.
          cycle

       end if

       ! *******************
       ! MULTISTART STRATEGY
       ! *******************
       !
       ! The relaxation of this subproblem is solved 'maxSolverCalls'
       ! times following the multistart strategy. Initial points are
       ! randomly generated according to the following rule:
       !
       ! We have two procedures to apply joint with multistart:
       !
       ! 1) The relaxation of the subproblem is solved at most
       ! 'maxSolverCalls' times and we accept the first feasible
       ! relaxed solution found. This strategy is activated by
       ! choosing FIRST_FEASIBLE as subproblem resolution strategy.
       !
       ! 2) The relaxation of the subproblem is solved exactly
       ! 'maxSolverCalls' times and we get the best feasible relaxed
       ! solution found among all solutions found. This strategy is
       ! activated by choosing BEST_FEASIBLE as subproblem resolution
       ! strategy.

       ! Initially, the value of the best relaxed solution found is
       ! set to "infinity".
       best_f = LARGEST_VALUE

       ! This variable indicates whether a solution for the relaxed
       ! problem was found.
       feasibleSolutionFound = .false.

       ! Initialize the seed to generate random initial points for
       ! this subproblem.
       if(iteration .gt. 1) then

          call getSeed(seed,l,u,n)

       else

          seed = real(84543787,kind=8)

       end if

       do multiStartIteration = 1,max(1,maxSolverCalls + maxSolverCallsInc*p%depth)

          ! Get a random initial point for this subproblem. If i = 1,
          ! we consider the initial point previous defined (stored).
          if(multiStartIteration .gt. 1) then

             call randomInitialPoint(n,x,seed)

          end if

          if(circlePackingInitialPoint) then

             call defCirclePackingProblem()

             call optimizationSolver(n,x,l,u,m,lambda,equatn,&
                  linear,rho,constf,epsfeas,epsopk,maxit,iter_2,f_2,g,&
                  gpsupn,cnorm,cnormu,geninfo,inform)

             optSolverIter = optSolverIter + iter_2
             call defRectanglePackingProblem()

          end if

          call solveProblem(p,n,x,l,u,m,lambda,equatn,linear,rho,constf,&
               epsfeas,epsopk,maxit,iter_,f,g,gpsupn,cnorm,cnormu,geninfo,&
               inform,numProblemsSolved)

          optSolverIter = optSolverIter + iter_

          ! Update lower bound on the objective function.
          if(f < lowerBound) then
             lowerBound = f
          end if

          ! Updating the best solution.
          if(satisfyConstraints(n,x,l,u)) then

             ! A feasible solution for the relaxed problem was found.
             feasibleSolutionFound = .true.

             if(f < best_f) then

                ! Update the best relaxed solution found.
                best_f = f
                best_x = x

             end if

             ! Verify whether the relaxed solution found is feasible
             ! for the main problem, that is, if it also satisfies the
             ! integrality constraints.

             if(satisfyDomain(n,x,mainDomain,discrete_sets) .and. &
                  isNumber(f)) then

                if(.not. solutionFound .or. f < objectiveValue) then

                   if(PRINT_) then
                      call cpu_time(time)
                      write(*,*) '   Solution updated:   ', f, ' time = ', time
                   end if

                   ! Update the incumbent solution.
                   objectiveValue = f
                   incumbent      = x
                   solutionFound  = .true.
                   iterationSolutionFound = iteration

                   if(objectiveValue .le. cutoff) then

                      exit

                   end if
                   
                end if
                
             else

                ! Round solution x.
                call round(n,x,round_x,mainDomain)

                if(solveRoundedSolution) then

                   ! FIX THE INTEGER VARIABLES AND SOLVE THE PROBLEM
                   ! ONLY FOR THE CONTINUOUS VARIABLES.
                   do i = 1,n

                      if(variableTypes(i) .eq. BINARY .or. &
                           variableTypes(i) .eq. INTEGER) then
          
                         l_(i) = round_x(i)
                         u_(i) = l_(i)

                      else

                         l_(i) = l(i)
                         u_(i) = u(i)

                      end if
       
                   end do

                   call optimizationSolver(n,round_x,l_,u_,m,lambda,equatn,&
                        linear,rho,constf,epsfeas,1.0D-8,maxit,iter_,f_,g,&
                        gpsupn,cnorm,cnormu,geninfo,inform)

                   numProblemsSolved = numProblemsSolved + 1

                end if

                ! Check whether the rounded solution 'round_x' is a
                ! feasible (integer) solution for the main problem.
                if(satisfyDomain(n,round_x,mainDomain,discrete_sets) .and. &
                     satisfyConstraints(n,round_x,l,u)) then

                   ! Evaluate the objective function at round_x.
                   call sevalal(n,round_x,m,lambda,rho,equatn,linear,&
                        round_f,inform)

                   if(inform .lt. 0) then
                      call deallocateDomain(mainDomain)
                      call destroyList(openProblemsList)
                      deallocate(discrete_sets)

                      if(branchingVariableStrategy .eq. PSEUDO_COSTS) then
                         call deallocatePseudoCosts()
                      end if

                      return
                   end if

                   ! Rounded solution is feasible for the main problem.
                   feasibleSolutionFound = .true.

                   ! Check whether the rounded solution is better than
                   ! the best relaxed solution found for this
                   ! subproblem.
                   if(round_f < best_f) then

                      ! Update the best solution found.
                      best_f = round_f
                      best_x = round_x

                   end if

                   ! Check whether the rounded solution is better than
                   ! the incumbent solution.
                   if(inform .eq. 0 .and. isNumber(f) .and. &
                     (.not. solutionFound .or. round_f < objectiveValue)) then

                      ! The rounded solution is feasible for the main
                      ! problem and it is better than the incumbent
                      ! solution.

                      if(PRINT_) then
                         call cpu_time(time)
                         write(*,*) '   Solution updated (rounded solution):',&
                              round_f, ' time = ', time
                      end if

                      ! Update the incumbent solution.                
                      objectiveValue = round_f
                      incumbent      = round_x
                      solutionFound  = .true.
                      iterationSolutionFound = iteration

                      if(objectiveValue .le. cutoff) then                         
                         exit
                      end if

                   end if

                end if ! End of testing rounded solution.

             end if

          end if

          ! Accept the first relaxed feasible solution if this
          ! strategy was chosen.
          if(feasibleSolutionFound .and. &
               multistartStrategy .eq. FIRST_FEASIBLE) then

             ! Exit from the multistart loop.
             exit

          end if

       end do

       ! best_x is the best feasible solution found fir the relaxed
       ! subproblem.
       f = best_f
       x = best_x

       if(iteration .eq. 1) then
          S0 = 0.0d0
          do i = 1,n
             if(variableTypes(i) .eq. BINARY .or. &
                  variableTypes(i) .eq. INTEGER) then
          
                frac = x(i) - real(floor(x(i)),kind=8)
                if(frac .le. 0.5d0) then
                   S0 = S0 + frac
                else
                   S0 = S0 + 1.0d0 - frac
                end if
             end if
          end do
       end if

! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!       EXIT
!       exit
! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

       ! **************************
       ! END OF MULTISTART STRATEGY
       ! **************************

       p%solved = .true.

       ! PRUNING BY INFEASIBILITY
       ! ------------------------
       !
       ! If no solution was found for the relaxed problem, this
       ! subproblem is declared to be infeasible.
       if(.not. feasibleSolutionFound) then

          pruningByInfeas = pruningByInfeas + 1

          call deallocateProblem(p)

          ! Continue with the next subproblem from the list.
          cycle

       end if

       ! UPDATE PSEUDO-COSTS
       ! -------------------
       if(branchingVariableStrategy .eq. PSEUDO_COSTS) then

          if(p%indexBranchingVariable .ge. 1 .and. &
               p%indexBranchingVariable .le. n) then

             call updatePseudoCost(p%indexBranchingVariable,&
                  p%branchingVariableValue,p%downBranch,&
                  p%lowerBound,f)

          end if

       end if

       if(objectiveValue .le. cutoff) then

          pruningByOpt = pruningByOpt + 1
          exit

       end if

       ! Check whether the solution found for the relaxed subproblem
       ! is a feasible solution for the main problem.
       satisfy_domain = satisfyDomain(n,x,mainDomain,discrete_sets)

       ! PRUNING BY OPTIMALITY
       ! ---------------------
       !
       ! If the feasible solution of the relaxed subproblem also
       ! satisfies the integrality constraints then it is a feasible
       ! solution for the main problem.
       if(isNumber(f) .and. satisfy_domain) then

          pruningByOpt = pruningByOpt + 1

          if(.not. solutionFound .or. f < objectiveValue) then

             if(PRINT_) then
                write(*,*) '   Iteration: ',iteration
                write(*,*) '   Solution updated:   ',f
             end if

             ! Update the incumbent solution.
             objectiveValue = f
             incumbent = x
             solutionFound = .true.
             iterationSolutionFound = iteration

             if(objectiveValue .le. cutoff) then

                exit

             end if
                   
          end if

          call deallocateProblem(p)

          ! Continue with the next subproblem from the list.
          cycle

       end if

       ! PRUNING BY BOUND
       ! ----------------
       !
       ! If the value of the solution of the relaxed subproblem is not
       ! less than the value of the incumbent solution, this node can
       ! be pruned.
       if(solutionFound .and. boundTolerance * f .ge. objectiveValue) then

          pruningByBoundCount = pruningByBoundCount + 1
          call deallocateProblem(p)

          ! Continue with the next subproblem from the list.
          cycle

       end if

       ! *******************
       ! BRANCHING PROCEDURE
       ! *******************

       ! Choose a branching variable following the strategy specified by
       ! 'branchingVariableStrategy' variable.
       call chooseBranchingVariable(x,mainDomain,discrete_sets,&
            variableType,branchingVariableStrategy,branchingIndex,options)

       m_ = size(p%domain%type) + 1

       ! BINARY VARIABLE
       ! ---------------
       if(variableType .eq. BINARY) then

          call branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID)

          if(RANDOM_DIRECTION .and. drand(seed_prob) .lt. 0.5d0) then

             data_%problem = p2
             call add(openProblemsList,data_)
             data_%problem = p1
             call add(openProblemsList,data_)

          else

             data_%problem = p1
             call add(openProblemsList,data_)
             data_%problem = p2
             call add(openProblemsList,data_)

          end if

       ! INTEGER VARIABLE
       ! ----------------
       else if(variableType .eq. INTEGER) then

          call branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID)

          if(RANDOM_DIRECTION .and. drand(seed_prob) .lt. 0.5d0) then
             data_%problem = p2
             call add(openProblemsList,data_)
             data_%problem = p1
             call add(openProblemsList,data_)
          else
             data_%problem = p1
             call add(openProblemsList,data_)
             data_%problem = p2
             call add(openProblemsList,data_)
          end if


       ! DISCRETE VARIABLE
       ! -----------------
       else if(variableType .eq. DISCRETE) then

          do j = 1,size(discrete_sets)

             if(discrete_sets(j)%index .eq. branchingIndex) then

                ! For each possible value for the branching variable, a
                ! subproblem is created and a constraint is added in
                ! order to fix the value of the variable.
                do k = 1,size(discrete_sets(j)%elements)

                   subdomain = createDomain(m_,&
                        (/ p%domain%indices,branchingIndex /),&
                        (/ p%domain%type,EQUAL /),&
                        (/ p%domain%value,discrete_sets(j)%elements(k) /))

                   p1 = createProblem(n,x,subdomain,nextID,p%id,f,&
                        branchingIndex,x(branchingIndex),.false.,p%depth + 1)
                   p1%geninfoFather = p%geninfo
                   data_%problem = p1
                   call add(openProblemsList,data_)

                   nextID = nextID + 1

                end do

                exit

             end if

          end do

       end if

       call deallocateProblem(p)

       maxSize = max(maxSize, getSize(openProblemsList))

    end do ! END OF MAIN LOOP

    call deallocateDomain(mainDomain)
    deallocate(discrete_sets)

    if(branchingVariableStrategy .eq. PSEUDO_COSTS) then
       call deallocatePseudoCosts()
    end if

    ! FIX THE INTEGER VARIABLES AND SOLVE THE PROBLEM ONLY FOR THE
    ! CONTINUOUS VARIABLES.
    do i = 1,n

       if(variableTypes(i) .eq. BINARY .or. variableTypes(i) .eq. INTEGER) then
          
          l(i) = incumbent(i)
          u(i) = l(i)

       else

          l(i) = lower(i)
          u(i) = upper(i)

       end if

    end do

    call optimizationSolver(n,incumbent,l,u,m,lambda,equatn,linear,rho,constf,&
         epsfeas,1.0D-8,maxit,iter_,objectiveValue,g,gpsupn,cnorm,cnormu,&
         geninfo,inform)

    totalNumNodesGenerated = totalNumNodesGenerated + numNodesGenerated
    totalNumProblemsSolved = totalNumProblemsSolved + numProblemsSolved

    call cpu_time(time)

    if(PRINT_) then
       write(*,*)
       write(*,*) 'Time:               ', time
       write(*,*) 'Solution found:     ', objectiveValue
       write(*,*) 'Nodes generated:    ', totalNumNodesGenerated
       write(*,*) 'Subproblems solved: ', totalNumProblemsSolved
       write(*,*)
    end if

! Print formats.
!100 format('Iteration     Objective')

!200 format (I9,1X,1PE20.10)

  end subroutine BranchAndBound

  ! ************************************************************************
  ! ************************************************************************

  subroutine branch(p,p1,p2,n,x,branchingIndex,variableType,f,nextID)

    use Problems

    implicit none

    ! SCALAR ARGUMENTS
    integer,      intent(in)     :: n,branchingIndex,variableType
    integer,      intent(in out) :: nextID
    real(kind=8), intent(in)     :: x(n),f

    ! OTHER ARGUMENTS
    type(problem), intent(in)  :: p
    type(problem), intent(out) :: p1,p2

    ! LOCAL SCALARS
    integer      :: m_
    real(kind=8) :: boundValue

    ! OTHER SCALARS
    type(domain) :: subdomain

    ! p1 must be the down branch and p2 must be the up branch.

    ! Number of constraints on the subproblems.
    m_ = size(p%domain%type) + 1

    ! BINARY VARIABLE
    ! ---------------
    if(variableType .eq. BINARY) then

       ! First problem: we set x_i = 0.
       subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),&
            (/ p%domain%type,EQUAL /),(/ p%domain%value,0.0d0 /))

       p1 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,&
            x(branchingIndex),.true.,p%depth + 1)

       p1%geninfoFather = p%geninfo

       nextID = nextID + 1

       ! Second problem: we set x_i = 1.
       subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),&
            (/ p%domain%type,EQUAL /),(/ p%domain%value,1.0d0 /))

       p2 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,&
            x(branchingIndex),.false.,p%depth + 1)

       p2%geninfoFather = p%geninfo

       nextID = nextID + 1

    ! INTEGER VARIABLE
    ! ----------------
    else if(variableType .eq. INTEGER) then

       boundValue = real(floor(x(branchingIndex)),kind=8)

       ! First problem. We add the following constraint:
       ! x_i <= floor(x_i).
       subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),&
            (/ p%domain%type,LESS_EQUAL /),(/ p%domain%value,boundValue /))

       p1 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,&
            x(branchingIndex),.true.,p%depth + 1)

       p1%geninfoFather = p%geninfo

       nextID = nextID + 1

       ! Second problem. We add the following constraint:
       ! x_i >= floor(x_i) + 1.
       subdomain = createDomain(m_,(/ p%domain%indices,branchingIndex /),&
            (/ p%domain%type,GREATER_EQUAL /),&
            (/ p%domain%value,boundValue + 1.0d0 /))

       p2 = createProblem(n,x,subdomain,nextID,p%id,f,branchingIndex,&
            x(branchingIndex),.false.,p%depth + 1)

       p2%geninfoFather = p%geninfo

       nextID = nextID + 1

    ! DISCRETE VARIABLE
    ! -----------------
!    else if(variableType .eq. DISCRETE) then
       ! TODO

    end if

  end subroutine branch

  ! ************************************************************************
  ! ************************************************************************

  subroutine solveProblem(p,n,x,l,u,m,lambda,equatn,linear,rho,constf,&
       epsfeas,epsopk,maxit,iter,f,g,gpsupn,cnorm,cnormu,geninfo,inform,&
       numProblemsSolved)

    use OptSolver
    use Problems

    implicit none

    ! SCALAR ARGUMENTS
    logical,      intent(in)     :: constf
    integer,      intent(in)     :: m,n,maxit
    integer,      intent(out)    :: geninfo,iter
    integer,      intent(in out) :: inform
    real(kind=8), intent(in)     :: epsfeas,epsopk
    real(kind=8), intent(out)    :: f,gpsupn
    real(kind=8), intent(in out) :: cnorm,cnormu
    integer,      intent(in out) :: numProblemsSolved

    ! ARRAY ARGUMENTS
    logical,      intent(in)     :: equatn(m),linear(m)
    real(kind=8), intent(in)     :: lambda(m),rho(m)
    real(kind=8), intent(in)     :: l(n),u(n)
    real(kind=8), intent(out)    :: g(n)
    real(kind=8), intent(in out) :: x(n)

    ! OTHER ARGUMENTS
    type(problem), intent(in out) :: p

    ! Call the solver for the relaxation of this subproblem.
    if(.not. p%solved) then

       call optimizationSolver(n,x,l,u,m,lambda,equatn,linear,rho,constf,&
            epsfeas,epsopk,maxit,iter,f,g,gpsupn,cnorm,cnormu,geninfo,inform)

       numProblemsSolved = numProblemsSolved + 1

    else

       ! This subproblem has already been solved.
       x = p%initialPoint
       f = p%solutionValue
       geninfo = p%geninfo

    end if

  end subroutine solveProblem

  ! ************************************************************************
  ! ************************************************************************

  subroutine round(n,x,round_x,domain_)

    use Problems

    implicit none

    ! SCALAR ARGUMENT
    integer, intent(in) :: n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in)  :: x(n)
    real(kind=8), intent(out) :: round_x(n)

    type(domain), intent(in) :: domain_

    ! LOCAL SCALARS
    integer      :: i, index
    real(kind=8) :: floor_x, ceil_x

    round_x = x

    do i = 1, size(domain_%type)

       index = domain_%indices(i)

       if(domain_%type(i) .eq. BINARY .or. domain_%type(i) .eq. INTEGER) then

          floor_x = real(floor(x(index)), kind=8)
          ceil_x  = real(ceiling(x(index)), kind=8)

          ! Round to the closest value.
          if(x(index) - floor_x .lt. ceil_x - x(index)) then

             round_x(index) = floor_x

          else

             round_x(index) = ceil_x

          end if

       end if

    end do
       
  end subroutine round

  ! ************************************************************************
  ! ************************************************************************

  subroutine printSolution(objectiveValue,x,n,iterations,time,solutionFound,&
       numNodesGenerated,numProblemsSolved)

    implicit none

    ! SCALAR ARGUMENTS
    integer,      intent(in) :: iterations,n
    real(kind=8), intent(in) :: objectiveValue
    real,         intent(in) :: time
    logical,      intent(in) :: solutionFound
    integer,      intent(in) :: numNodesGenerated,numProblemsSolved

    ! ARRAY ARGUMENT
    real(kind=8) :: x(n)

    ! LOCAL SCALAR
    integer :: i

    if(PRINT_) then

       if(solutionFound) then

          print *,'\nBest solution found:\nx ='
          do i = 1,n
             write(*,300) x(i)
          end do
          write(*,400) objectiveValue,iterations,time

       else 

          print *,'\nSolution not found.'
          write(*,500) iterations,time

       end if

       write(*,600) numNodesGenerated,numProblemsSolved

    end if

300 format('    ', 1PE30.20)

400 format(/, 'Objective function value = ', 0PF30.20, &
           /, 'Number of iterations     = ', I14,      &
           /, 'Runtime                  = ', 0PF30.20)

500 format(/, 'Number of iterations     = ', I14,      &
           /, 'Runtime                  = ', 0PF30.20)

600 format(/, 'Nodes generated          = ', I14,      &
           /, 'Problems solved          = ', I14)

  end subroutine printSolution

  ! ************************************************************************
  ! ************************************************************************

  subroutine initialize(n,numBinary,numInteger,domain_,discrete_sets,&
       variableTypes)
  
    use Problems
    use Structures

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n,numBinary,numInteger

    ! ARRAY ARGUMENTS
    integer,   intent(in)               :: variableTypes(n)
    type(set), allocatable, intent(out) :: discrete_sets(:)

    ! OTHER ARGUMENTS
    type(domain), intent(out) :: domain_

    ! LOCAL SCALARS
    integer :: i,count

    ! Allocate memory for the domain structure.
    allocate(domain_%indices(numBinary + numInteger))
    allocate(domain_%type   (numBinary + numInteger))
    allocate(domain_%value  (numBinary + numInteger))

    allocate(discrete_sets(0))

    count = 1

    do i = 1,n

       if(variableTypes(i) .eq. BINARY) then

          domain_%type(count)    = BINARY
          domain_%indices(count) = i
          count = count + 1

       else if(variableTypes(i) .eq. INTEGER) then

          domain_%type(count)    = INTEGER
          domain_%indices(count) = i
          count = count + 1

       end if

    end do

  end subroutine initialize

  ! ************************************************************************
  ! ************************************************************************

  logical function satisfyConstraints(n, x, lower, upper)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n), lower(n), upper(n)

    satisfyConstraints = satisfyBounds(n, x, lower, upper)

  end function satisfyConstraints

  ! ************************************************************************
  ! ************************************************************************

  logical function satisfyBounds(n, x, lower, upper)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n), lower(n), upper(n)

    ! LOCAL SCALARS
    integer      :: i

    do i = 1, n

       if(x(i) .lt. lower(i) - EPSILON_FEAS .or. &
          x(i) .gt. upper(i) + EPSILON_FEAS) then

          satisfyBounds = .false.
          return

       end if

    end do

    satisfyBounds = .true.

  end function satisfyBounds

  ! ************************************************************************
  ! ************************************************************************

  logical function satisfyDomain(n, x, domain_, discrete_sets)

    use Problems
    use Structures
    use Auxiliar

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    type(set),    intent(in) :: discrete_sets(:)

    ! OTHER ARGUMENTS
    type(domain), intent(in) :: domain_

    integer :: i, j, k, index
    logical :: satisfyDiscrete

    satisfyDomain = .true.
    
    do i = 1, size(domain_%type)

       index = domain_%indices(i)

       ! INTEGER variable.
       if(domain_%type(i) .eq. INTEGER) then

          if(.not. isInteger(x(index))) then

             satisfyDomain = .false.
             return

          end if

       ! BINARY variable.
       else if(domain_%type(i) .eq. BINARY) then

          if(.not. isBinary(x(index))) then

             satisfyDomain = .false.
             return

          end if

       ! DISCRETE variable.
       else if(domain_%type(i) .eq. DISCRETE) then

          satisfyDiscrete = .false.

          ! What is the set of variable x_index?
          do j = 1, size(discrete_sets)

             if(discrete_sets(j)%index .eq. index) then

                do k = 1, size(discrete_sets(j)%elements)

                   if(abs(x(index) - discrete_sets(j)%elements(k)) &
                        .le. EPSILON_INT) then

                      satisfyDiscrete = .true.
                      exit

                   end if

                end do

                if(.not. satisfyDiscrete) then

                   satisfyDomain = .false.
                   return

                end if

             end if

          end do

       end if

    end do

  end function satisfyDomain

  ! ************************************************************************
  ! ************************************************************************

  subroutine addBounds(domain_, lower, upper, n)

    use Problems

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in out) :: lower(n), upper(n)

    ! OTHER ARGUMENTS
    type(domain), intent(in) :: domain_

    ! LOCAL SCALARS
    integer :: i, index

    do i = 1, size(domain_%type)

       index = domain_%indices(i)

       if(domain_%type(i) .eq. EQUAL) then

          lower(index) = domain_%value(i)
          upper(index) = lower(index)

       else if(domain_%type(i) .eq. LESS_EQUAL) then

          upper(index) = domain_%value(i)

       else if(domain_%type(i) .eq. GREATER_EQUAL) then

          lower(index) = domain_%value(i)

       end if

    end do

  end subroutine addBounds

  ! ************************************************************************
  ! ************************************************************************

  subroutine getSeed(seed, lower, upper, n)

    use Problems

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    double precision, intent(out) :: seed

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lower(n), upper(n)

    ! LOCAL SCALARS
    integer :: i

    seed = 0.0d0

    do i = 1, n

       seed = seed + real(i,kind=8) * (abs(lower(i)) + abs(upper(i)))

    end do

    seed = real(floor(seed),kind=8)

    if(seed .lt. 0.0d0 .or. seed .gt. 2.0d0**31.0d0 - 1.0d0) then

       seed = 84543787.0d0

    end if

  end subroutine getSeed

  ! ************************************************************************
  ! ************************************************************************

  logical function isNumber(number)

    implicit none

    ! This function verifies whether the specified number is valid,
    ! that is, if the number is not +infinity, -infinity and NaN.

    ! SCALAR ARGUMENT
    real(kind=8), intent(in) :: number

    ! Parameters of the subroutine:
    ! =============================
    !
    ! On Entry:
    ! =========
    !
    ! NUMBER real(kind=8): The number to be verified.
    ! -------------------
    !
    ! Return:
    ! =======
    !
    ! Return false if the number is either +infinity, -infinity or NaN.
    
    if(.not. (abs(number) .le. LARGEST_VALUE)) then

       isNumber = .false.

    else 

       isNumber = .true.

    end if    

  end function isNumber

! ************************************************************************
! ************************************************************************

end module MINLP
