R/essR.R

In MEIGOR: MEIGO - MEtaheuristics for bIoinformatics Global Optimization

essR <-
  function(problem,opts=list(maxeval=NULL,maxtime=NULL),...){
    # Function   : essR 
    # R code of the eSS optimization code from: Process Engineering Group IIM-CSIC	
    # Created on : 09/06/2011
    # Last Update: 25/07/2012
    # Email      : josea.egea@gmail.com
    # 
    # (c) CSIC, Spanish Council for Scientific Research, Spain
    # 
    # Global optimization algorithm for MINLPs based on Scatter Search
    # 
    #   essR attempts to solve problems of the form:
    #       min F(x)  subject to:  ceq(x) = 0 (equality constraints)
    #        x                     c_L <= c(x) <= c_U (inequality constraints)
    #                              x_L <= x <= x_U (bounds on the decision variables)
    # 
    # Constraint functions, if applicable, must be declared in the same script as the
    # objective function as a second output argument:
    # e.g., myfunction <- function(x){
    #				calculate fx - scalar containing the objective function value
    #				calculate gx - vector (or empty) containing the constraints values
    #				return(list(fx,gx))
    #				}
    # 
    # Please have a look at the manual before using essR
    # 
    #       USAGE: Results <- essR(problem,opts,p1,p2,....,pn);
    # 
    # INPUT PARAMETERS:
    # ****************
    #   problem - List containing problem settings
    #       problem$f   = Name of the file containing the objective
    #                     function (String)
    #       problem$x_L = Lower bounds of decision variables (vector)
    #       problem$x_U = Upper bounds of decision variables (vector)
    #       problem$x_0 = Initial point(s) (optional; vector or matrix)
    #       problem.$f_0 = Function values of initial point(s) (optional) These values MUST correspond to feasible points.
    #
    #       NOTE:The dimension of f_0 and x_0 may be different. For example, if
    #       we want to introduce 5 initial points but we only know the values
    #       for 3 of them, x_0 would have 5 rows whereas f_0 would have only 3
    #       elements. In this example, it is mandatory that the first 3 rows of x_0 correspond
    #       to the values of f_0
    # 
    #       Additionally, fill the following fields if your problem has
    #       non-linear constraints
    #       problem$neq     = Number of equality constraints (Integer; do not define it
    #                         if there are no equality constraints)
    #       problem$c_L     = Lower bounds of nonlinear inequality constraints
    #                        (vector)
    #       problem$c_U     = Upper bounds of nonlinear inequality constraints
    #                         (vector)
    #       problem$int_var = Number of integer variables (Integer)
    #       problem$bin_var = Number of binary variables  (Integer)
    #       problem$vtr     = Objective function value to be reached (optional)
    #
    # 
    # NOTE: The order of decision variables is x <- c(cont, int, bin)
    # 
    #   opts - List containing options (if set as opts <- numeric(0) default options
    #          will be loaded). See the script of default options to know their values
    # 
    #       User options
    #           opts$maxeval    = Maximum number of function evaluations
    #                             (Default 1000)
    #           opts$maxtime    = Maximum CPU time in seconds (Default 60)
    #           opts$iterprint  = Print each iteration on screen: 0-Deactivated
    #                             1-Activated (Default 1)
    #           opts$weight     = Weight that multiplies the penalty term added
    #                             to the objective function in constrained
    #                             problems (Default 1e6)
    #           opts$log_var    = Indexes of the variables which will be used
    #                             to generate diverse solutions in different
    #                             orders of magnitude (vector)
    #           opts$tolc       = Maximum absolute violation of the constraints
    #                             (Default 1e-5)
    #           opts$prob_bound = Probability (0-1) of biasing the search towards
    #                             the bounds (Default 0.5)
    #           opts$save_results = Saves final results in eSSR_report.RData. (Binary; Default = 1)
    #           opts$inter_save = Saves results of intermediate iterations in eSSR_report.RData
    #                             Useful for very long runs. (Binary; Default = 0)
    # 
    #       Global options
    #           opts$dim_refset         = Number of elements in Refset
    #                                     (Integer; automatically calculated)
    #           opts$ndiverse           = Number of solutions generated by the
    #                                     diversificator (Default 10*nvar)
    #           opts$combination        = Type of combination of Refset
    #                                     elements (Default 1)
    #                                     1: hyper-rectangles
    #                                     2: linear combinations
    # 
    #       Local options
    #           opts$local_solver               = Choose local solver
    #                                             0: Local search deactivated (Default),
    #                                             "NM", "BFGS", "CG", "LBFGSB",
    #                                             "SA","SOLNP", "DHC", "NLS2SOL"
    #           opts$local_tol                  = Level of tolerance in local
    #                                             search
    #           opts$local_iterprint            = Print each iteration of local
    #                                             solver on screen (Binary;
    #                                             default = 0)
    #           opts$local_n1                   = Number of iterations
    #                                             before applying
    #                                             local search for the 1st time
    #                                             (Default 1)
    #           opts$local_n2                   = Minimum number of iterations
    #                                             in the global
    #                                             phase between 2 local calls
    #                                             (Default 10)
    #           opts$local_balance              = Balances between quality (=0)
    #                                             and diversity (=1) for
    #                                             choosing initial points for
    #                                             the local search (default
    #                                             0.5)
    #           opts$local_finish               = Applies local search to the
    #                                             best solution found once the
    #                                             optimization if finished
    #                                             (same values as
    #                                             opts.local.solver)
    #           opts$local_bestx                = When activated (i.e. =1) only
    #                                             applies local search to the
    #                                             best solution found to
    #                                             date,ignoring filters
    #                                             (Default=0)
    # 
    # 
    #   p1,p2... :  optional input parameters to be passed to the objective
    #   function
    # 
    # 
    # OUTPUT PARAMETERS:
    # *****************
    # A file called "eSSR_report.Rdata" is generated containing.
    # 
    #   Results - Data frame containing results
    #   problem - Data frame containing problem settings
    #   opts    - Data frame containing all options
    #
    # 
    # Fields in Results
    #       Results$fbest                   = Best objective function value
    #                                         found after the optimization
    #       Results$xbest                   = Vector providing the best
    #                                         function value
    #       Results$cpu_time                = Time in seconds consumed in the
    #                                         optimization
    #       Results$f                       = Vector containing the best
    #                                         objective function value after each
    #                                         iteration
    #       Results$x                       = Matrix containing the best vector
    #                                         after each iteration
    #       Results$time                    = Vector containing the cpu time
    #                                         consumed after each iteration
    #       Results$neval                   = Vector containing the number of
    #                                         function evaluations after each
    #                                         iteration
    #       Results$numeval                 = Number of function evaluations
    #       Results$local_solutions         = Local solutions found by the
    #                                         local solver (in rows)
    #       Results$local_solutions_values  = Function values of the local
    #                                         solutions
    #       Results$end_crit                = Criterion to finish the
    #                                         optimization
    #                                         1: Maximal number of function
    #                                         evaluations achieved
    #                                         2: Maximum allowed CPU Time
    #                                         achieved
    #                                         3: Value to reach achieved
    # 
    #
    # REFERENCES
    # 
    # If you use essR and publish the results, please cite the following papers:
    #
    # Egea, J.A., Marti, R., Banga, J.R. (2010) An evolutionary method for
    # complex-process optimization. Computers & Operations Research 37(2):
    # 315-324.
    # 
    # Egea, J.A., Balsa-Canto, E., Garcia, M.S.G., Banga, J.R. (2009) Dynamic
    # optimization of nonlinear processes with an enhanced scatter search
    # method. Industrial & Engineering Chemistry Research 49(9): 4388-4401. 
			
	# Function   : eSS Release 2014A
	print('eSS R2014A - Enhanced Scatter Search');
	
    stopOptimization<-0;
    
    # Initialize time
    cpu_time=proc.time()[3];
    #Take the 3rd component of cpu_time for the total elapsed time
    
    #Initialize the number of function evaluations and a stop variable
    nfuneval <- 0;
    fin <- 0;

    if (!is.numeric(opts$maxeval) && !is.numeric(opts$maxtime)){
      
      cat("WARNING:Either opts.maxeval or opts.maxtime must be defined as a stop criterion \n")
      cat("Define at least one of these two options and rerun \n")
      Results <- numeric(0);
      return(Results)
      stop()
    } else{
      if (!is.numeric(opts$maxeval)){
        maxeval <- 1e12;
      } else{
        maxeval <- opts$maxeval;
      }
      if (!is.numeric(opts$maxtime)){
        maxtime <- 1e12;
      } else{
        maxtime <- opts$maxtime;
      }
      if(!(maxtime>0)){
        Results <- numeric(0);
        return(Results);
      }
      if(!(maxeval>0)){
        Results <- numeric(0);
        return(Results);
      }
    }
    
    cat("\n ------------------------------------------------------------------------------ \n")
    cat(" essR - Enhanced Scatter Search in R \n");
    cat("<c> IIM-CSIC, Vigo, Spain -  email: gingproc@iim.csic.es \n");
    cat("------------------------------------------------------------------------------ \n\n");
    
    #Load default values for the options
    
    default <- ssm_defaults();
    
    nargin <- length(as.list(match.call())) -1;
    
    if (nargin<2){opts <- numeric(0)};
    
    #Set all optiions
    opts<-ssm_optset(default,opts);
    
    x_U<-problem$x_U;
    x_L<-problem$x_L;
    
    #Check if the bounds have the same dimension
    if (length(x_U)!= length(x_L)){
      cat("Upper and lower bounds have different dimensions!!! \n")
      cat("EXITING")
      Results<-numeric(0);
      return(Results)
      stop()
    } else{
      #Number of decision variables
      nvar<-length(x_L);
    }
    
    prob_names<-names(problem);
    
    temp<-match("x_0",prob_names);
    if (is.na(temp)){x_0<-numeric(0)} else{x_0<-problem$x_0}
    
    temp<-match("f_0",prob_names);
    if (is.na(temp)){
      f_0<-numeric(0);
      fbest<-Inf;
      xbest<-rep(Inf,nvar)
    } else{
      f_0<-problem$f_0;
      fbest<-min(f_0);
      iii<-which.min(f_0);
      xbest<-x_0[iii,];
    }
    
    
    temp<-match("neq",prob_names);
    if (is.na(temp)){neq<-0} else{neq<-problem$neq}
    temp<-match("c_U",prob_names);
    if (is.na(temp)){c_L<-numeric(0); c_U<-numeric(0)} else{c_U<-problem$c_U; c_L<-problem$c_L}
    temp<-match("int_var",prob_names);
    if (is.na(temp)){int_var<-0} else{int_var<-problem$int_var}
    temp<-match("bin_var",prob_names);
    if (is.na(temp)){bin_var<-0} else{bin_var<-problem$bin_var}
    temp<-match("vtr",prob_names);
    if (is.na(temp)){vtr<-numeric(0)} else{vtr<-problem$vtr}
    
    #User options
    iterprint=opts$iterprint;
    prob_bound=opts$prob_bound;
    save_results=opts$save_results;
    inter_save=opts$inter_save;
    plot_results=opts$plot;
    weight=opts$weight;
    tolc=opts$tolc;
    log_var=opts$log_var;
    
    #Global options
    dim_refset=opts$dim_refset;
    if (is.character(dim_refset)){
      z<-c(-10*nvar/1,-1,1);
      nnn<-polyroot(z);
      iii<-which(Re(nnn)>0);
      dim_refset<-ceiling(Re(nnn[iii]));
      
      
      
      #%% is the modulus of a division e.g. 7%%4 provides 3 (like mod(7,4) in Matlab)
      if (dim_refset%%2){
        dim_refset<-dim_refset+1
        if (iterprint){
          cat("Refset size automatically calculated:",dim_refset, "\n");
        }
      }
    }
    
    
    
    ndiverse<-opts$ndiverse;
    if (is.character(ndiverse)){
      ndiverse=10*nvar;
      if (iterprint){
        cat("Number of diverse solutions automatically calculated:",ndiverse,"\n");
      }
    }
    if (ndiverse<dim_refset){ndiverse<-dim_refset}
    
    combin=opts$combination;
    
    #Local options
    local_solver=opts$local_solver;
    if(is.character(local_solver)){local_solver<-toupper(local_solver)}
    local_tol=opts$local_tol;
    local_iterprint=opts$local_iterprint;
    local_n1<-opts$local_n1;
    local_n2<-opts$local_n2;
    if (is.character(local_n1)){local_n1<-1}
    if (is.character(local_n2)){local_n2<-10}
    local_balance=opts$local_balance;
    
    if (!length(opts$local_finish)){
      local_finish<-local_solver;
    } else{
      local_finish=opts$local_finish;
    }
    
    
    
    local_bestx=opts$local_bestx;
    
    #If there are equality constraints
    if (neq){
      #Set the new bounds for all the constraints
      c_L<-c(rep(0,neq),c_L);
      c_U<-c(rep(0,neq),c_U)
    }
    
    
    nconst<-length(c_U);
    
    fobj<-problem$f;
    
    #Check if the objective function has 2 output arguments in constrained problems
    if (nconst){
      n_out_f<-do.call(fobj,list(x_L,...));
      if (length(n_out_f)<2){
        cat("For constrained problems the objective function must have at least 2 output arguments \n");
        cat("EXITING \n");
        Results<-numeric(0);
        return(Results)
        stop()
      }
    }		 
    
    #Change to log
    xl_log<-x_L;
    xu_log<-x_U;
    
    if (length(log_var)){
      aaa<-which(!xl_log[log_var]);
      if (length(aaa)){
        xl_log[log_var[aaa]]<-1e-10;
      }
      
      xl_log[log_var]<-log(xl_log[log_var]);
      xu_log[log_var]<-log(xu_log[log_var]);
    }
    
    #Check if the objective function has 3 output arguments when NL2SOL is invoked
    if(match(local_solver,"NL2SOL", nomatch=0) | match(local_finish,"NL2SOL", nomatch=0)){
      n_out_f<-do.call(fobj,list(x_L,...));
      nfuneval<-nfuneval+1;
      if (length(n_out_f)<3){
        cat("NL2SOL requires 3 output arguments: f, g, and the vector of residuals \n");
        cat("EXITING \n");
        Results<-numeric(0);
        return(Results)
        stop()
      }
    }
    
    
    # TRANSFORM THIS INTO R CODE WHEN THE LOCAL SEARCHES ARE IMPLEMENTED!!!!!!
    #	%Build auxiliary files in case local search is activated
    #	if local_solver
    #	ssm_aux_local(problem.f,x_L,x_U,c_L,c_U,neq,local_solver,nvar,varargin{:});
    #	end
    
    #Initialize variables
    stage_1<-1;
    n_minimo<-0;
    n_critico<-local_n1;
    
    
    local_solutions<-numeric(0);
    local_solutions_values<-numeric(0);
    initial_points<-numeric(0);
    
    Results<-list(f=numeric(0), x=numeric(0), time=numeric(0), neval=numeric(0));
    
    #Possible combinations among the dim_refset elements in Refset, taken in pairs
    ncomb<-t(combn(1:dim_refset,2));
    MaxSubSet<-(dim_refset^2-dim_refset)/2;
    MaxSubSet2<-2*MaxSubSet;
    
    
    #We generate 5 solutions
    solutions<-matrix(0,ndiverse+5,nvar);
    
    solutions[1:5,]<-(matrix(runif(nvar*5),5,nvar)+1:5-1)/5;
    
    solutions[6:(ndiverse+5),]<-matrix(runif(ndiverse*nvar),ndiverse,nvar);
    
    solutions<-solutions*kronecker(matrix(1,ndiverse+5,1),t((xu_log-xl_log)))+kronecker(matrix(1,ndiverse+5,1),t(xl_log));
    
    
    for (i in 1:(ndiverse+5)){
      solutions[i,log_var]=exp(solutions[i,log_var]);
    }
    
    l_f_0<-length(f_0);
    if (is.vector(x_0) && length(x_0)){l_x_0<-1} else if(is.matrix(x_0)){l_x_0<-nrow(x_0)} else{l_x_0<-0};
    if (!length(l_x_0)){l_x_0<-0}
    
    #Put x_0 WITHOUT their corresponding f_0 values at the beginning of solutions
    solutions<-rbind(x_0, solutions)
    if (l_f_0){solutions<-solutions[-seq(1,l_f_0),]}
    temp<-l_x_0-l_f_0;
    
    sol<-matrix(0,ndiverse+5+temp,nvar);
    sol_val<-rep(0,ndiverse+5+temp);
    sol_val_pen<-rep(0,ndiverse+5+temp);
    sol_pen<-rep(0,ndiverse+5+temp);
    sol_nlc<-matrix(0,ndiverse+5+temp,nconst);
    
    
    
    counter<-1;
    
  
    for (i in 1:(ndiverse+5+temp)){
      res<-ssm_evalfc(solutions[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
      val<-res[[1]];
      val_penalty<-res[[2]];
      pena<-res[[3]];
      nlc<-res[[4]];
      includ<-res[[5]];
      x<-res[[6]];
      
      nfuneval<-nfuneval+1;
      if (includ){
        sol[counter,]<-x;
        sol_val[counter]<-val;
        sol_val_pen[counter]<-val_penalty;
        sol_pen[counter]<-pena;
        sol_nlc[counter,]<-nlc;
        counter<-counter+1;
      }
    }
    
    
    nsol<-nrow(sol);
    
    if ((counter-1)<nsol){
      sol<-sol[-seq(counter-1,nsol),];
      sol_val<-sol_val[-seq(counter-1,nsol)];
      sol_val_pen<-sol_val_pen[-seq(counter-1,nsol)];
      sol_pen<-sol_pen[-seq(counter-1,nsol)];
      sol_nlc<-sol_nlc[-seq(counter-1,nsol),];
    }
    
    #Add points with f_0 values. We assume feasiblity
    if (l_f_0){
      sol<-rbind(x_0[1:l_f_0,], sol);
      sol_val<-c(f_0, sol_val);
      #We assume feasible points in f_0
      sol_val_pen<-c(f_0, sol_val_pen);
      sol_pen<-c(rep(0,l_f_0),sol_pen);
      #We do not know these values, but we do not care. We assume feasibility
      sol_nlc<-rbind(NaN*matrix(0,l_f_0,nconst),sol_nlc);     
    }
    
    #Initialize Refset
    Refset<-matrix(0,dim_refset,nvar);
    Refset_values<-rep(0,dim_refset);
    Refset_values_penalty<-rep(0,dim_refset);
    Refset_nlc<-numeric(0);
    penalty<-numeric(0);
    
    temp<-sort(sol_val_pen,index.return=TRUE);
    iii<-temp[[2]];
    first_members<-ceiling(dim_refset/2);
    last_members<-dim_refset-first_members;
    
    #Create the initial sorted Refset
    Refset<-sol[iii[1:first_members],];
    Refset_values<-sol_val[iii[1:first_members]];
    Refset_values_penalty<-sol_val_pen[iii[1:first_members]];
    penalty<-sol_pen[iii[1:first_members]];
    Refset_nlc<-sol_nlc[iii[1:first_members],];
    
    #The rest of Refset members are chosen randomly
    iii<-iii[-(1:first_members)];
    temp<-1:length(iii);
    ooo<-sample(temp);
    
    #The comment below comes from matlab code
    #ooo=randperm(ndiverse+2+size(x_0,1)-first_members);
    Refset<-rbind(Refset, sol[iii[ooo[1:last_members]],]);
    Refset_values<-c(Refset_values, sol_val[iii[ooo[1:last_members]]]);
    Refset_values_penalty<-c(Refset_values_penalty, sol_val_pen[iii[ooo[1:last_members]]]);
    penalty<-c(penalty, sol_pen[iii[ooo[1:last_members]]]);
    Refset_nlc<-rbind(Refset_nlc, sol_nlc[iii[ooo[1:last_members]],]);
    
    #Check best value
    ggg<-which(penalty<=tolc);
    if (length(ggg)){
      fbest<-min(Refset_values_penalty[ggg]);
      iii<-which.min(Refset_values_penalty[ggg]);
      xbest<-Refset[ggg[iii],];
    }
    
    
    iter<-0;
    if (iterprint){
      cat("Initial Pop: NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
    }
    
    Results$f=c(Results$f,fbest);
    Results$x=rbind(Results$x,xbest);
    Results$time=c(Results$time, proc.time()[3]-cpu_time);
    Results$neval=c(Results$neval, nfuneval);
    
    
    #Sort the Refset
    temp<-sort(Refset_values_penalty,index.return=TRUE);
    Refset_values_penalty<-temp[[1]];
    I<-temp[[2]];
    Refset_values<-Refset_values[I];
    penalty<-penalty[I];
    Refset<-Refset[I,];
    Refset_nlc<-Refset_nlc[I,];
    
    if (combin==1){nrand=nvar} else if (combin==2) {nrand=1}
    
    use_bestx<-0;
    
    index1<-ncomb[,1];
    index2=ncomb[,2];
    index<-c(index1,index2);
    diff_index<-index2-index1;
    
    #Minimum ppp
    lb_p<-1;
    #Defines maximum ppp. max(ppp)= lb_p+ st_p
    st_p<-0.75;
    #Example: lb_p=1.25 and st_p=0.5 varies ppp from 1.25 to 1.75
    ppp<-st_p*(diff_index-1)/(dim_refset-2)+lb_p;
    hyper_x_L<-kronecker(matrix(1,MaxSubSet,1),t(x_L));
    hyper_x_U<-kronecker(matrix(1,MaxSubSet,1),t(x_U));
    refset_change<-rep(0,dim_refset);
    while (!fin){
      child<-matrix(0,MaxSubSet2,nvar);
      child_values<-rep(0,MaxSubSet2);
      child_values_penalty<-rep(0,MaxSubSet2);
      child_penalty<-rep(0,MaxSubSet2);
      child_nlc<-matrix(0,MaxSubSet2,nconst);
      
      child_parent_index<-rep(0,MaxSubSet2);
      
      members_update<-rep(0,dim_refset);
      
      counter<-1;
      continuar<-0;
      
      while (!continuar){
        #Sort Refset
        #This is mandatory because the variable index is also sorted
        counter<-1;
        temp<-sort(Refset_values_penalty,index.return=TRUE);
        Refset_values_penalty<-temp[[1]];
        I<-temp[[2]];
        Refset_values<-Refset_values[I];
        penalty<-penalty[I];
        Refset<-Refset[I,];
        Refset_nlc<-Refset_nlc[I,];
        
        refset_change<-refset_change[I];
        
        parents_index1<-Refset[index1,];
        parents_index2<-Refset[index2,];
        parents_index1_values_penalty<-Refset_values_penalty[index1];
        parents_index2_values_penalty<-Refset_values_penalty[index2];
        
        denom22<-pmax(abs(parents_index1),abs(parents_index2));			
        denom22[which(!denom22)]<-1;
        
        AAA<-abs((parents_index1-parents_index2)/denom22);
        
        #Compute the maximum value of AAA by rows
        temp<-apply(AAA,1,max);
        #if any of these maximum values are lower than a threshold, then keep the indexes
        BBB<-which(temp<1e-3);
        
        
        if (length(BBB)){
          #Check the index of worst element among the twins of index2
          index_refset_out<-max(index2[BBB]);
          includ=0;
          while (!includ){
            #Generate a new random solution
            new_refset_member<-runif(nvar)*(xu_log-xl_log)+xl_log;
            new_refset_member[log_var]=exp(new_refset_member[log_var]);
            res<-ssm_evalfc(new_refset_member,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
            val<-res[[1]];
            val_penalty<-res[[2]];
            pena<-res[[3]];
            nlc<-res[[4]];
            includ<-res[[5]];
            x<-res[[6]];
            nfuneval<-nfuneval+1;
            if (includ){
              Refset[index_refset_out,]<-x;
              Refset_values[index_refset_out]<-val;
              Refset_values_penalty[index_refset_out]<-val_penalty;
              Refset_nlc[index_refset_out,]<-nlc;
              penalty[index_refset_out]<-pena;
              members_update[index_refset_out]<-0;
            }
          }
        } else{
          continuar=1;
        }
        
        candidate<-Refset;
        candidate_values<-Refset_values;
        candidate_values_penalty<-Refset_values_penalty;
        candidate_nlc<-Refset_nlc;
        candidate_penalty<-penalty;
        candidate_update<-rep(0,dim_refset);
        
        members_update<-rep(1,dim_refset);
        
        factor<-kronecker(matrix(1,1,nvar),ppp)*(parents_index2-parents_index1)/1.5;
        v1<-parents_index1-factor;
        v2<-parents_index2-factor;
        v3<-2*parents_index2-parents_index1-factor;
        
        aaa<-which(v1<hyper_x_L);
        if (length(aaa)){
          rand_aaa<-runif(length(aaa));
          AAA<-which(rand_aaa>prob_bound);
          aaa<-aaa[AAA];
          v1[aaa]=hyper_x_L[aaa];
        }
        
        bbb<-which(v1>hyper_x_U);
        if (length(bbb)){
          rand_bbb<-runif(length(bbb));
          BBB<-which(rand_bbb>prob_bound);
          bbb<-bbb[BBB];
          v1[bbb]=hyper_x_U[bbb];
        }
        
        ccc<-which(v3<hyper_x_L);
        if (length(ccc)){
          rand_ccc<-runif(length(ccc));
          CCC<-which(rand_ccc>prob_bound);
          ccc<-ccc[CCC];
          v3[ccc]=hyper_x_L[ccc];
        }
        
        ddd<-which(v3>hyper_x_U);
        if (length(ddd)){
          rand_ddd<-runif(length(ddd));
          DDD<-which(rand_ddd>prob_bound);
          ddd<-ddd[DDD];
          v3[ddd]=hyper_x_U[ddd];
        }
        
        if (nrand>1){
          new_comb1<-v1+(v2-v1)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
          new_comb2=v2+(v3-v2)*matrix(runif(MaxSubSet*nrand),MaxSubSet,nrand);
        } else{
          new_comb1=v1+(v2-v1)*runif(1);
          new_comb2=v2+(v3-v2)*runif(1);
        }
        
        new_comb<-rbind(new_comb1,new_comb2);
        
        for (i in 1:MaxSubSet2){
          #Evaluate
          res<-ssm_evalfc(new_comb[i,],x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
          val<-res[[1]];
          val_penalty<-res[[2]];
          pena<-res[[3]];
          nlc<-res[[4]];
          includ<-res[[5]];
          x<-res[[6]];
          
          nfuneval<-nfuneval+1;
          if (includ){
            child[counter,]<-x;
            child_values[counter]<-val;
            child_values_penalty[counter]<-val_penalty;
            child_penalty[counter]<-pena;
            child_nlc[counter,]<-nlc;
            child_parent_index[counter]<-index[i];
            
            counter<-counter+1;
            
            if (val_penalty<candidate_values_penalty[index[i]]){
              #Update candidate
              #This can be accelerated saving the index and replacing at the end with the best one
              candidate[index[i],]<-x;
              candidate_values_penalty[index[i]]<-val_penalty;
              candidate_values[index[i]]<-val;
              candidate_penalty[index[i]]<-pena;
              candidate_nlc[index[i],]<-nlc;
              candidate_update[index[i]]<-1;
              members_update[index[i]]<-0;
            }
          }
        }
        #Check the stop criterion
        if (nfuneval>= maxeval){
          fin<-1;
        } else if(proc.time()[3]-cpu_time>=maxtime){
          fin<-2;
        } else if(length(vtr)){
          if (fbest<=vtr){
            fin<-3;
          }
        }
        
        child[-(1:(counter-1)),]<-numeric(0);
        child_values[-(1:(counter-1))]<-numeric(0);
        child_values_penalty[-(1:(counter-1))]<-numeric(0);
        child_penalty[-(1:(counter-1))]<-numeric(0);
        child_parent_index[-(1:(counter-1))]<-numeric(0);
        
        
        members_to_update<-which(candidate_update==1);
        if (length(members_to_update)){
          for (i in 1:length(members_to_update)){
            
            res <- ssm_beyond(Refset[members_to_update[i],],candidate[members_to_update[i],],candidate_values_penalty[members_to_update[i]],
                              fobj,nrand,tolc,weight,x_L,x_U,c_L,c_U,nconst,int_var,bin_var,nfuneval,prob_bound,nvar,...)
            vector<-res[[1]];
            vector_value<-res[[2]];
            vector_penalty<-res[[3]];
            vector_value_penalty<-res[[4]];
            vector_nlc<-res[[5]];
            new_child<-res[[6]];
            new_child_value<-res[[7]];
            new_child_penalty<-res[[8]];
            new_child_value_penalty<-res[[9]];
            new_child_nlc<-res[[10]];
            nfuneval<-res[[11]];
            
            if (length(vector)){
              Refset[members_to_update[i],]<-vector;
              Refset_values[members_to_update[i]]<-vector_value;
              Refset_values_penalty[members_to_update[i]]<-vector_value_penalty;
              penalty[members_to_update[i]]<-vector_penalty;
              Refset_nlc[members_to_update[i],]<-vector_nlc;
            } else{
              Refset[members_to_update[i],]<-candidate[members_to_update[i],];
              Refset_values[members_to_update[i]]<-candidate_values[members_to_update[i]];
              Refset_values_penalty[members_to_update[i]]<-candidate_values_penalty[members_to_update[i]];
              penalty[members_to_update[i]]<-candidate_penalty[members_to_update[i]];
              Refset_nlc[members_to_update[i],]<-candidate_nlc[members_to_update[i],];
            }
          }
        }
        
        fbest_lastiter=fbest;
        
        ggg<-which(penalty<=tolc);
        if (length(ggg)){
          fbest_new<-min(Refset_values_penalty[ggg]);
          hhh<-which.min(Refset_values_penalty[ggg]);
          if (fbest_new<fbest){
            xbest<-Refset[ggg[hhh],];
            fbest<-fbest_new;
            use_bestx<-1;
            if (inter_save){
              Results$fbest<-fbest;
              Results$xbest<-xbest;
              Results$numeval<-nfuneval;
              Results$end_crit<-fin;
              Results$cpu_time<-cpu_time;
              Results$Refset$x<-Refset;                         
              Results$Refset$f<-Refset_values;                  
              Results$Refset$fpen<-Refset_values_penalty;       
              Results$Refset$const<-Refset_nlc;                
              Results$Refset$penalty<-penalty;
              if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
            }
          }
        }
        
        iter<-iter+1;
        n_minimo<-n_minimo+1;
        
        if (iterprint){
          cat("Iteration:",iter, "NFunEvals:", nfuneval, "Bestf:",fbest,"CPUTime:",proc.time()[3]-cpu_time,"Mean:", mean(Refset_values_penalty),"\n");
        }
        
        ddd<-which(members_update>0);
        eee<-which(!members_update);
        refset_change[ddd]<-refset_change[ddd]+1;
        refset_change[eee]<-0;
        fff<-which(refset_change>=20);
        if (length(fff)){
          to_replace<-length(fff);
          replaced<-0;
          while (replaced<to_replace){
            newsol<-runif(nvar)*(xu_log-xl_log)+xl_log;
            newsol[log_var]=exp(newsol[log_var]);
            res<-ssm_evalfc(newsol,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
            val<-res[[1]];
            val_penalty<-res[[2]];
            pena<-res[[3]];
            nlc<-res[[4]];
            includ<-res[[5]];
            x<-res[[6]];
            
            nfuneval<-nfuneval+1;
            if (includ){
              Refset[fff[replaced+1],]<-x;
              Refset_values[fff[replaced+1]]<-val;
              Refset_values_penalty[fff[replaced+1]]<-val_penalty;
              penalty[fff[replaced+1]]<-pena;
              Refset_nlc[fff[replaced+1],]<-nlc;
              refset_change[fff[replaced+1]]<-0;
              replaced<-replaced+1;
            }
          }
        }
        if (fbest<fbest_lastiter){
          Results$f<-c(Results$f,fbest);
          Results$x<-rbind(Results$x,xbest);
          Results$time<-c(Results$time, proc.time()[3]-cpu_time);
          Results$neval<-c(Results$neval,nfuneval);
          #if plot_results==1
          #           stairs(Results.time,Results.f)
          #           xlabel('CPU Time (s)');
          #            ylabel('Objective Function Value (-)');
          #            title(strcat('Best f:  ',num2str(fbest)))
          #            drawnow
          #			end
        }
        
        
        #******************************************************************************************
        #************************* LOCAL FILTERS ROUTINE ******************************************
        call_filters<-0;
        if(is.character(local_solver) && n_minimo>=n_critico){
          call_filters<-1;
          if (local_bestx && !use_bestx){
            call_filters<-0;
          }
        }
        
        if (call_filters){
          if (stage_1==1){
            #In this first stage we apply the local search over the best solution found so far
            minimo<-min(Refset_values_penalty);
            I<-which.min(Refset_values_penalty);
            x0<-Refset[I,];
            f0<-minimo;
            #Switch to stage 2
            n_critico<-local_n2;
            if (local_bestx){
              #not to go directly to stage 2 the next time we pass over here unless we have really improved fbest
              use_bestx<-0;    
            }
            #use_bestx is a flag indicating that the best solution has been improved
          }
          else if (stage_2==1){
            if (local_bestx && use_bestx){
              x0<-xbest;
              f0<-fbest;
              use_bestx<-0;
            } else{
              #Calculate distance between our children and local_solutions+initial points
              local_init<-rbind(local_solutions, initial_points);
              xuxl<-x_U-x_L;
              
              #child<-rbind(child,new_child);
              #						child_values<-c(child_values,new_child_value);
              #						child_values_penalty<-c(child_values_penalty,new_child_value_penalty);
              #						child_penalty<-c(child_penalty,new_child_penalty);
              #                       child_nlc<-rbind(child_nlc,new_child_nlc);
              
              child_norm<-child/kronecker(matrix(1,nrow(child),1),t(xuxl));        
              
              local_init_norm<-local_init/kronecker(matrix(1,nrow(local_init),1),t(xuxl));
              
              distance<-eucl_dist(child_norm,local_init_norm);
              
              #Now we need to sort the children according the their distances to local minima and intial points
              
              
              #Here we compute the minimum euclidean distance between child i and all the vector in initial_points+local_solutions
              temp<-apply(distance,1,min);
              #Now we sort those minimum distances and favour the highest ones to favour diverse points
              res<-sort(temp,decreasing=TRUE,index.return=TRUE);
              ooo<-res[[2]];
              
              res<-sort(child_values_penalty,index.return=TRUE);
              uuu<-res[[2]];
              
              
              n_child<-length(child_values_penalty)
              child_points<-rep(0,n_child);
              www<-local_balance;
              for (i in 1:n_child){
                child_points[ooo[i]]<-child_points[ooo[i]]+www*i;
                child_points[uuu[i]]<-child_points[uuu[i]]+(1-www)*i;
              }
              
              I<-which.min(child_points);
              x0<-child[I,];
              f0<-child_values_penalty[I];
            }
          }
          
          if (iterprint){
            cat("Call local solver:", toupper(local_solver),"\n")
            if (child_penalty[I]>tolc){
              cat("Initial point function value:",child_values[I], "(INFEASIBLE)","\n");
            } else{
              cat("Initial point function value:",f0,"\n");
            }
            tic<-proc.time()[3];
          }
          
          #AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
          res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_solver,local_iterprint,local_tol,weight,nconst,tolc,...);
          x<-res[[1]];
          fval<-res[[2]];
          numeval<-res[[3]];
          
          initial_points<-rbind(initial_points, x0);
          nfuneval<-nfuneval+numeval;
          
          #Evaluate the local solution
          res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
          val<-res[[1]];
          val_penalty<-res[[2]];
          pena<-res[[3]];
          nlc<-res[[4]];
          includ<-res[[5]];
          x<-res[[6]];
          nfuneval<-nfuneval+1;
          if (includ && pena<=tolc){
            if (iterprint){
              cat("Local solution function value:",val,"\n");
              
            }
            if (stage_1==1){
              stage_1=0;
              stage_2=1;
            }
              if (val_penalty<fbest){
              fbest<-val_penalty;
              xbest<-x;
              if (fbest<fbest_lastiter){
                Results$f<-c(Results$f,fbest);
                Results$x<-rbind(Results$x,xbest);
                Results$time<-c(Results$time,proc.time()[3]-cpu_time);
                Results$neval<-c(Results$neval,nfuneval);
                #if plot_results==1
                #							stairs(Results.time,Results.f)
                #							xlabel('CPU Time (s)');
                #						ylabel('Objective Function Value (-)');
                #						title(strcat('Best f:  ',num2str(fbest)))
                #						drawnow
                #						end
              }
            }
          } else{
            cat("Local solution function value:","UNFEASIBLE","\n");
            
          }
          cat("Number of function evaluations in the local search:",numeval,"\n");
          cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");
          
          #Check the stopping criterion
          if (nfuneval>= maxeval){
            fin<-1;
          } else if(proc.time()[3]-cpu_time>=maxtime){
            fin<-2;
          } else if(length(vtr)){
            if (fbest<=vtr){
              fin<-3;
            }
          }
          if (includ){
            if (length(local_solutions)){
              #Add the found local solution to the list (if it is a new one)
              adicionar_local<-1;
              res<-ssm_isdif2(x,local_solutions,1e-2,1);
              f<-res[[1]];
              ind<-res[[2]];
              #If it is too close to a previously found local solution, we do not add it
              if (f){
                adicionar_local<-0
              } 
            } else {
              #Add it if no local solutions have been found yet
              adicionar_local<-1;
              
            }
            
            if (adicionar_local){
              
              #choose either this (replace the worst element in Refset)
              
              #         [aaa,jjj]=max(Refset_values_penalty);
              #
              #         if val_penalty<Refset_values_penalty(jjj)
              #             Refset(jjj,:)=x;
              #             Refset_values(jjj)=val;
              #             Refset_values_penalty(jjj)=val_penalty;
              #							Refset_nlc(jjj,:)=nlc;
              #				%             penalty(jjj)=pena;
              #				%         end
              #				
              #				%or this (Replace the initial point)
              #				
              #				%         if val_penalty<f0
              #				%             child(I,:)=x;
              #				%             child_values(I)=val;
              #				%             child_values_penalty(I)=val_penalty;
              #				%             child_nlc(I,:)=nlc;
              #				%         end
              #				
              
              #or this (Replace the parent of the initial point)
              if (val_penalty<Refset_values_penalty[child_parent_index[I]]){
                Refset[child_parent_index[I],]<-x;
                Refset_values[child_parent_index[I]]<-val;
                Refset_values_penalty[child_parent_index[I]]<-val_penalty;
                Refset_nlc[child_parent_index[I],]<-nlc;
                penalty[child_parent_index[I]]<-pena;
                refset_change[child_parent_index[I]]<-0;
              }
              local_solutions<-rbind(local_solutions,x);
              local_solutions_values<-c(local_solutions_values,fval);
            }
          }
          
          n_minimo<-0
        }
        #************************* END OF LOCAL FILTERS ROUTINE ***********************************
        #******************************************************************************************
        
        
        
        fbest_lastiter<-fbest;
        
        #Check the stopping criterion
        if (nfuneval>= maxeval){
          fin<-1;
        } else if(proc.time()[3]-cpu_time>=maxtime){
          fin<-2;
        } else if(length(vtr)){
          if (fbest<=vtr){
            fin<-3;
          }
        }
        
        
        if (fin){
          cat("\n");
          if (is.character(local_finish) && fin<3){
            local_tol=3;
            
            if (is.numeric(xbest)){
              x0=xbest;
              f0=fbest;
              }else{
                x0=Refset[1,];
                f0=Refset_values_penalty[1];
            }
            if (iterprint){
              cat("Final local refinement: ", toupper(local_finish),"\n")
              cat("Initial point function value: ", f0,"\n");
              tic<-proc.time()[3];
            }
            
            res<-ssm_localsolver(x0,x_L,x_U,c_L,c_U,neq,int_var,bin_var,fobj,local_finish,local_iterprint,local_tol,weight,nconst,tolc,...);
            x<-res[[1]];
            fval<-res[[2]];
            numeval<-res[[3]];
            
            
            #Evaluate the local solution
            res<-ssm_evalfc(x,x_L,x_U,fobj,nconst,c_L,c_U,tolc,weight,int_var,bin_var,nvar,...);
            val<-res[[1]];
            val_penalty<-res[[2]];
            pena<-res[[3]];
            nlc<-res[[4]];
            includ<-res[[5]];
            x<-res[[6]];
            nfuneval<-nfuneval+1;
            if (includ && pena<=tolc){
              if (iterprint){
                cat("Local solution function value:",val,"\n");
                
              }
              if (val_penalty<fbest){
                fbest<-val_penalty;
                xbest<-x;
              }
            } else{
              cat("Local solution function value:","UNFEASIBLE","\n");
              
            }
            cat("Number of function evaluations in the local search:",numeval,"\n");
            cat("CPU Time of the local search:",proc.time()[3]-tic,"seconds \n\n");            
          }

          cpu_time<-proc.time()[3]-cpu_time;
          Results$f<-c(Results$f,fbest);
          Results$x<-rbind(Results$x, xbest);
          Results$neval<-c(Results$neval,nfuneval);
          Results$time<-c(Results$time,cpu_time);
          
          if (iterprint){
            msg<-switch(fin, "1" = "Maximum number of function evaluations achieved",
                        "2" = "Maximum computation time achieved", "3"="Desired function value achieved")
            cat(msg,"\n");
            cat("Best solution value", fbest, "\n");
            cat("Decision vector",xbest, "\n");
            cat("CPU time", cpu_time,"\n");
            cat("Number of function evaluations", nfuneval,"\n");
          }
          #	fclose all;
          Results$fbest<-fbest;
          Results$xbest<-xbest;
          Results$numeval<-nfuneval;
          Results$end_crit<-fin;
          Results$cpu_time<-cpu_time;
          Results$Refset$x<-Refset;
          Results$Refset$f<-Refset_values;
          Results$Refset$fpen<-Refset_values_penalty;
          Results$Refset$const<-Refset_nlc;
          Results$Refset$penalty<-penalty;
          if (save_results){save(opts, problem, Results, file = "eSSR_report.RData")}
          
          #if plot_results==1 |plot_results==2
          #					stairs(Results.time,Results.f)
          #					xlabel('CPU Time (s)');
          #					ylabel('Objective Function Value (-)');
          #					title(strcat('Best f:  ',num2str(fbest)))
          #					end
          #					save ess_report problem opts Results
          #					%Delete aux files
          #					if local_solver & isempty(strmatch(local_solver,local_finish));
          #					ssm_delete_files(local_solver,c_U);
          #					end
          #					%     keyboard
          return(Results)
          #					end
        }
      }
    }
  }