.packageName <- "subselect"
anneal<-function(mat, kmin, kmax=kmin, nsol=1, niter=1000,
exclude=NULL, include=NULL, improvement=TRUE, 
setseed = FALSE, cooling=0.05, temp=1,
coolfreq=1, criterion="RM", pcindices="first_k", initialsol=NULL){

###############################
# general validation of input #
###############################

        validation(mat, kmin, kmax, exclude, include, criterion, pcindices)

##########################################################################
# more specific validation of input for the anneal and improve functions #
##########################################################################

        implog<-FALSE
        validannimp(kmin, kmax, nsol, exclude, nexclude, include, ninclude, initialsol,implog)


#############################################################
# validation specific for the input to the anneal function  #
#############################################################

        if (!(as.integer(coolfreq) == coolfreq) | (coolfreq < 1)) stop("\n The cooling frequency must be a non-negative integer")
        if (cooling<=0 || cooling >=1) stop("\n values of cooling must be between 0 and 1")



##########################################
# initializations for Fortran subroutine #
##########################################

        if (setseed == TRUE) {set.seed(2,kind="default")} 
        valores<-rep(0.0,length(kmin:kmax)*nsol)    
        vars<-rep(0,nsol*length(kmin:kmax)*kmax)
        bestval<-rep(0.0,length(kmin:kmax))
        bestvar<-rep(0,kmax*length(kmin:kmax))

##################################
# call to the Fortran subroutine #
##################################

        Fortout<-.Fortran("anneal",as.integer(criterio),as.integer(p),
          as.double(as.vector(mat)),
          as.integer(kmin),as.integer(kmax),as.double(valores),
          as.integer(vars),as.double(bestval),as.integer(bestvar),
          as.integer(nexclude),as.integer(exc),as.integer(ninclude),
          as.integer(inc),as.integer(nsol),as.integer(niter),
          as.logical(improvement),
          as.double(cooling),as.double(temp),
          as.integer(coolfreq),as.integer(length(pcindices)),
          as.integer(pcindices),as.logical(esp),as.logical(silog),
          as.integer(as.vector(initialsol)),PACKAGE="subselect")

########################
# preparing the output #
########################

        valores<-matrix(nrow=nsol,ncol=length(kmin:kmax),Fortout[[6]])
        dimnames(valores)<-list(paste("Solution",1:nsol,sep=" "),paste("card.",kmin:kmax,sep=""))
        variaveis<-array(Fortout[[7]],c(nsol,kmax,length(kmin:kmax)))
        dimnames(variaveis)<-list(paste("Solution",1:nsol,sep=" "),paste("Var",1:kmax,sep="."),paste("Card",kmin:kmax,sep="."))
        bestval<-Fortout[[8]]
        names(bestval)<-paste("Card",kmin:kmax,sep=".")
        bestvar<-t(matrix(nrow=kmax,ncol=length(kmin:kmax),Fortout[[9]]))
        dimnames(bestvar)<-list(paste("Card",kmin:kmax,sep="."),paste("Var",1:kmax,sep="."))
        output<-list(variaveis,valores,bestval,bestvar,match.call())
        names(output)<-c("subsets","values","bestvalues","bestsets","call")
        output}


gcd.coef<-function(mat, indices, pcindices = NULL)
{
#   
#       calcula o GCD entre um subconjunto ("indices") de variaveis
#       e um subconjunto ("pcindices") das CPs de todas as variaveis, 
#       cuja matriz de covariancias e "mat".

#  error checking

         if  (sum(!(as.integer(indices) == indices)) > 0) stop("\n The variable indices must be integers")         
         if  (!is.null(pcindices) & (sum(!(as.integer(pcindices) == pcindices)) > 0)) stop("\n The PC indices must be integers")
     if (!is.matrix(mat)) {
         stop("Data is missing or is not given in matrix form")}
     if (dim(mat)[1] != dim(mat)[2]) {
         mat<-cov(mat)
         warning("Data must be given as a covariance or correlation matrix. \n It has been assumed that you wanted the covariance matrix of the \n data matrix which was supplied.")
     }        

# body of function

# initializations

      if (!is.null(pcindices)) {
       if (!is.vector(pcindices)) stop("If Principal Components are user-specified, only one set of PCs is allowed for each function call")
      }
  
     dvsmat <- svd(mat)
      tr<-function(mat){sum(diag(mat))}
      gcd.1d<-function(mat,indices, pcindices){
             if (is.null(pcindices)) {pcindices <- 1:sum(!indices==0)}
             indices<-indices[!indices == 0]
             svdapprox <- function(mat, indices) {
             t(dvsmat$v[, indices] %*% (t(dvsmat$u[, indices]) * dvsmat$d[indices]))
                           }
             sum(diag(solve(mat[indices, indices]) %*% svdapprox(mat, 
             pcindices)[indices, indices]))/sqrt(length(indices) * 
             length(pcindices))
        }
      dimension<-length(dim(indices))

# output for each dimension of input array

      if (dimension > 1){
         gcd.2d<-function(mat,subsets,pcindices){
             apply(subsets,1,function(indices){gcd.1d(mat,indices,pcindices)})
            }  
           if (dimension > 2) {               
            gcd.3d<-function(mat,array3d,pcindices){
             apply(array3d,3,function(subsets){gcd.2d(mat,subsets,pcindices)})
            }
            output<-gcd.3d(mat,indices,pcindices)
           }
           if (dimension == 2) {output<-gcd.2d(mat,indices,pcindices)}
      }

      if (dimension < 2) {output<-gcd.1d(mat,indices,pcindices)}
      output
}
genetic<-function(mat, kmin, kmax=kmin, popsize=100, nger=100,
mutate=FALSE, mutprob=0.01, maxclone=5, exclude=NULL, include=NULL,
improvement=TRUE, setseed= FALSE,  criterion="RM", pcindices="first_k",
initialpop=NULL){

###############################
# general validation of input #
###############################

        validation(mat, kmin, kmax, exclude, include, criterion, pcindices)

##############################################################
# more specific validation of input for the genetic function #
##############################################################

        validgenetic(kmin, kmax, popsize, mutprob, exclude, nexclude, include, ninclude, initialpop)


##############################################
# initializations for the Fortran subroutine #
##############################################

        if (setseed == TRUE) set.seed(2,kind="default")
        valores<-rep(0.0,length(kmin:kmax)*popsize)
        vars<-rep(0,popsize*length(kmin:kmax)*kmax)
        bestval<-rep(0.0,length(kmin:kmax))
        bestvar<-rep(0,kmax*length(kmin:kmax))
        kabort<-kmax+1

##################################
# call to the Fortran subroutine #
##################################

        Fortout<-.Fortran("genetic",as.integer(criterio),as.integer(p),
          as.double(as.vector(mat)),
          as.integer(kmin),as.integer(kmax),as.double(valores),
          as.integer(vars),as.double(bestval),as.integer(bestvar),
          as.integer(nexclude),as.integer(exc),as.integer(ninclude),
          as.integer(inc),as.integer(popsize),as.integer(nger),
          as.integer(maxclone),as.logical(mutate),
          as.double(mutprob),as.logical(improvement),
          as.integer(length(pcindices)),as.integer(pcindices),as.logical(esp),
          as.integer(kabort),as.logical(pilog),
          as.integer(as.vector(initialpop)),
          PACKAGE="subselect")

########################
# preparing the output #
########################

        kabort<-Fortout[[23]]
        valores<-matrix(nrow=popsize,ncol=length(kmin:kmax),Fortout[[6]])
        dimnames(valores)<-list(paste("Solution",1:popsize,sep=" "),paste("card.",kmin:kmax,sep=""))
        variaveis<-array(Fortout[[7]],c(popsize,kmax,length(kmin:kmax)))
        dimnames(variaveis)<-list(paste("Solution",1:popsize,sep=" "),paste("Var",1:kmax,sep="."),paste("Card",kmin:kmax,sep="."))
        bestval<-Fortout[[8]]
        names(bestval)<-paste("Card",kmin:kmax,sep=".")
        bestvar<-t(matrix(nrow=kmax,ncol=length(kmin:kmax),Fortout[[9]]))
        dimnames(bestvar)<-list(paste("Card",kmin:kmax,sep="."),paste("Var",1:kmax,sep="."))
        output<-list(variaveis[,1:(kabort-1),1:(kabort-kmin)],valores[,1:(kabort-kmin)],bestval[1:(kabort-kmin)],bestvar[1:(kabort-kmin),1:(kabort-1)],match.call())
        names(output)<-c("subsets","values","bestvalues","bestsets","call")
        if (kabort > kmin) output}


improve<-function(mat, kmin, kmax=kmin, nsol=1, exclude=NULL,
include=NULL, setseed = FALSE, criterion="RM", pcindices="first_k",
initialsol=NULL){


###############################
# general validation of input #
###############################

        validation(mat, kmin, kmax, exclude, include, criterion, pcindices)

##########################################################################
# more specific validation of input for the anneal and improve functions #
##########################################################################

        implog<-TRUE
        validannimp(kmin, kmax, nsol, exclude, nexclude, include, ninclude, initialsol, implog)



###########################################
# initializations for Fortran subroutine  #
###########################################

        if (setseed == TRUE) set.seed(2,kind="default")
        valores<-rep(0.0,(kmax-kmin+1)*nsol)
        vars<-rep(0,nsol*length(kmin:kmax)*kmax)
        bestval<-rep(0.0,length(kmin:kmax))
        bestvar<-rep(0,kmax*length(kmin:kmax))

###############################
# call to Fortran subroutine  #
###############################

        Fortout<-.Fortran("improve",as.integer(criterio),as.integer(p),
          as.double(as.vector(mat)),
          as.integer(kmin),as.integer(kmax),as.double(valores),
          as.integer(vars),as.double(bestval),as.integer(bestvar),
          as.integer(nexclude),as.integer(exc),as.integer(ninclude),
          as.integer(inc),as.integer(nsol),
          as.integer(length(pcindices)),as.integer(pcindices),
          as.logical(esp),as.logical(silog),
          as.integer(as.vector(initialsol)),PACKAGE="subselect")

########################
# preparing the output #
########################

        valores<-matrix(ncol=length(kmin:kmax),nrow=nsol,Fortout[[6]])
        dimnames(valores)<-list(paste("Solution",1:nsol,sep=" "),paste("card.",kmin:kmax,sep=""))
        variaveis<-array(Fortout[[7]],c(nsol,kmax,length(kmin:kmax)))
        dimnames(variaveis)<-list(paste("Solution",1:nsol,sep=" "),paste("Var",1:kmax,sep="."),paste("Card",kmin:kmax,sep="."))
        bestval<-Fortout[[8]]
        names(bestval)<-paste("Card",kmin:kmax,sep=".")
        bestvar<-t(matrix(nrow=kmax,ncol=length(kmin:kmax),Fortout[[9]]))
        dimnames(bestvar)<-list(paste("Card",kmin:kmax,sep="."),paste("Var",1:kmax,sep="."))
        output<-list(variaveis,valores,bestval,bestvar,match.call())
        names(output)<-c("subsets","values","bestvalues","bestsets","call")
        output}



leaps <- function(mat,kmin=1,kmax=ncol(mat)-1,nsol=1,exclude=NULL,include=NULL,
                  criterion="RM",pcindices=NULL,timelimit=15)
{

#############################################################
# validation specific for the input to the leaps function  #
#############################################################

 if ((criterion == "3") || (criterion == "GCD") || (criterion == "gcd") || (criterion == "Gcd"  || (criterion == 3))) {  
        if (!is.null(pcindices)) {if (sum(pcindices == "first_k") > 0) stop("\n The 'first_k' option is not available in the 'leaps' function: \n PC indices must be explicitely set. \n")}}
	if (timelimit <= 0) {stop("\n The time limit argument must be a positive real number")} 


###############################
# general validation of input #
###############################

        validation(mat, kmin, kmax, exclude, include, criterion, pcindices)



##########################################
# Initializations for the C++ subroutine #
##########################################

         klength    <- kmax-kmin+1

         subsets    <- integer(nsol*kmax*klength);   dim(subsets)    <- c(nsol,kmax,klength)
         values     <- double(nsol*klength);         dim(values)     <- c(nsol,klength)
         bestvalues <- double(klength);             
         bestsets   <- integer(klength*kmax);        dim(bestsets)   <- c(klength,kmax) 

         dimnames(subsets)   <-  list(paste("Solution",1:nsol,sep=" "),paste("Var",1:kmax,sep="."),paste("Card",kmin:kmax,sep="."))
         dimnames(values)    <-  list(paste("Solution",1:nsol,sep=" "),paste("card.",kmin:kmax,sep=""))
         names(bestvalues)   <-  paste("Card",kmin:kmax,sep=".")
         dimnames(bestsets)  <-  list(paste("Card",kmin:kmax,sep="."),paste("Var",1:kmax,sep="."))


############################
# Call to the C subroutine #
############################

 	 Cout <- .C("leaps",
            as.double (mat),
            as.integer(kmin),
            as.integer(kmax),
            as.integer(nsol),
            as.integer(exclude),
            as.integer(include),
            as.integer(nexclude),
            as.integer(ninclude),
            as.character(criterion),
            as.integer(pcindices),
            as.integer(length(pcindices)),
            as.integer(p),
	    as.double(timelimit),
	    found = logical(1),	    
            subsets,
            values,
            bestvalues,
            bestsets,
            PACKAGE="subselect"
        ) 

#######################################
# Preparing and returning the output  #
#######################################
        
	 if (Cout$found == FALSE) {
	    stop("\n Leaps was not able to complete the search within the specified time limit.\n Either increase this limit or try one of the available meta-heuristics \n") }
         output <- c(Cout[15:18],match.call())
         names(output) <- c("subsets","values","bestvalues","bestsets","call")
         output

}


rm.coef<-function(mat, indices)
{

#   Computes the matrix correlation between data matrices and their 
#   regression on a subset of their variables. Expected input is a
#   variance-covariance (or correlation) matrix. 

#  error checking


  if  (sum(!(as.integer(indices) == indices)) > 0) stop("\n The variable indices must be integers")
  if (!is.matrix(mat)) {
         stop("Data is missing or is not given in matrix form")}
     if (dim(mat)[1] != dim(mat)[2]) {
         mat<-cov(mat)
         warning("Data must be given as a covariance or correlation matrix. \n It has been assumed that you wanted the covariance matrix of the \n data matrix which was supplied.")
       }
      tr<-function(mat){sum(diag(mat))}
      rm.1d<-function(mat,indices){
        sqrt(tr((mat %*% mat)[indices,indices] %*% solve(mat[indices,indices]))/tr(mat))
        }
      dimension<-length(dim(indices))
      if (dimension > 1){
         rm.2d<-function(mat,subsets){
             apply(subsets,1,function(indices){rm.1d(mat,indices)})
            }  
             if (dimension > 2) {
               rm.3d<-function(mat,array3d){
                   apply(array3d,3,function(subsets){rm.2d(mat,subsets)})
                 }
               output<-rm.3d(mat,indices)
              }
              if (dimension == 2) {output<-rm.2d(mat,indices)}
      }

      if (dimension < 2) {output<-rm.1d(mat,indices)}
      output
}
rv.coef<-function(mat, indices)
{
#    Computes Escoufier's RV-coefficient for the configuration of
#    points defined by n observations of a set of p variables, and by
#    the regression of all variables on a subset of k variables
#    (given by \code{indices}). 

#  error checking

     if  (sum(!(as.integer(indices) == indices)) > 0) stop("\n The variable indices must be integers")
     if (!is.matrix(mat)) {
         stop("Data is missing or is not given in matrix form")}
     if (dim(mat)[1] != dim(mat)[2]) {
         mat<-cov(mat)
         warning("Data must be given as a covariance or correlation matrix. \n It has been assumed that you wanted the covariance matrix of the \n data matrix which was supplied.")
       }
      tr<-function(mat){sum(diag(mat))}
      rv.1d<-function(mat,indices){
             mat2 <- (mat %*% mat)[indices, indices]
             invmatk <- solve(mat[indices, indices])
             sqrt(tr(mat2 %*% invmatk %*% mat2 %*% invmatk)/tr(mat %*% mat))
        }
      dimension<-length(dim(indices))
      if (dimension > 1){
         rv.2d<-function(mat,subsets){
             apply(subsets,1,function(indices){rv.1d(mat,indices)})
            }  
             if (dimension > 2) {
               rv.3d<-function(mat,array3d){
                   apply(array3d,3,function(subsets){rv.2d(mat,subsets)})
                 }
               output<-rv.3d(mat,indices)
              }
              if (dimension == 2) {output<-rv.2d(mat,indices)}
      }

      if (dimension < 2) {output<-rv.1d(mat,indices)}
      output
}
validannimp<-function(kmin, kmax, nsol, exclude, nexclude, include, ninclude, initialsol, implog){


##################################################################
# validation of input that is specific to the anneal and improve #
# functions (which admit similar possible initial solutions)     #
##################################################################



##################################################################
# initializations when initial solutions are specified; checking #
# nature of initialsol (and how to interpret it) and for         #
# conflicts with the exclude and include requirements (that must #
# be respected by the initial solutions)                         #
##################################################################

        if  (!is.null(initialsol) & (sum(!(as.integer(initialsol) == initialsol)) > 0)) stop("\n The initial solutions must be specified as integers indicating variable numbers")
        if (is.null(initialsol)) {silog <- FALSE}
        else {silog <- TRUE  # initial solutions have been specified by user          
#########################################################
# checking for the presence of variables that are to be #
# forcefully excluded                                   #
#########################################################

            if ((nexclude != 0) & (sum(exclude == rep(initialsol,rep(length(exclude),length(as.vector(initialsol))))) !=0)) stop("\n the specified initial solutions contain variables that are to be excluded")

#########################################################
# how to deal with various formats of input (of initial #
# solutions) for a single cardinality                   #
#########################################################

          dimsol<-dim(initialsol)
              
          if (length(dimsol) > 3) {
            stop("\n Can't handle arrays of more than 3 dimensions")}
          if (kmin == kmax) {

#############################################################
# inital solution is a vector: must be repeated if nsol > 1 #
#############################################################


               if (is.vector(initialsol)) {
                   if (length(initialsol) != kmax) 
                     {stop("\n The specified initial and final solutions have different cardinalities")}
                   else
                     {if (nsol > 1) {
                     initialsol<-matrix(nrow=nsol,ncol=kmax,rep(initialsol,rep(nsol,kmax)))
                     if (implog==TRUE) {warning("\n a single initial subset necessarily produces nsol identical \n final solutions in the improve algorithm")}}
                    }
               }


##############################
# initial solution is array? #
##############################

           if (is.array(initialsol)) {

###########################
# initialsol is 3-d array #
###########################

           if (length(dimsol) == 3) {
            if (dim(initialsol)[[3]] > 1) stop("\n The input array of initial solutions must have dimensions as \n (1 or nsol) x k x 1 when a single cardinality is requested")
            else initialsol<-matrix(nrow=dim(initialsol)[[1]],ncol=dim(initialsol)[[2]],initialsol)
            }
           else

##################################################
# initialsol is matrix: must be of form nsol x k #
##################################################

                 {if (dim(initialsol)[[2]] != kmax) stop("\n Input matrix of initial solutions must have as many columns as variables in the requested subset")
                  else 
                      {if  (dim(initialsol)[[1]] != nsol) {
                         if (dim(initialsol)[[1]] == 1) {
                          initialsol<-matrix(nrow=nsol,ncol=dim(initialsol)[[2]],rep(initialsol,rep(nsol,kmax))) 
                          if (implog==TRUE) {warning("\n a single initial subset necessarily produces nsol identical \n final solutions in the improve algorithm")}
                         }
                           else {stop("\n The number of initial solutions can only be 1 or nsol (number of final solutions requested)")}}}}
              }}



##########################################################
# how to deal with various formats of input (initialsol) #
# if more than one cardinality is requested              #
##########################################################

             else  # (if kmax > kmin), i.e., more than one cardinality requested

###############################################
# initial solution cannot be a single vector  #
###############################################

              {if (is.vector(initialsol)) stop("\n There must be initial solutions for all cardinalities requested")

               if (is.array(initialsol)) {

###############################
# initial solution 3-d array? #
###############################


                  if (length(dim(initialsol)) == 3) {
                   if ((dim(initialsol)[[3]] != length(kmin:kmax)) | (dim(initialsol)[[2]] != kmax) | ((dim(initialsol)[[1]] != nsol) & (dim(initialsol)[[1]] > 1))) stop("\n The input array of initial solutions must have dimensions as \n nsol x kmax x no. of different cardinalities requested")
                   else 
                    { if (dim(initialsol)[[1]] != nsol) {
                       initialsol<-array(dim=c(nsol,kmax,length(kmin:kmax)),rep(initialsol,rep(nsol,kmax*length(kmin:kmax))))
                      if (implog==TRUE) {warning("\n a single initial subset necessarily produces nsol identical \n final solutions in the improve algorithm")}
                      } 
                  }}
                else 


############################
# if initalsol is a matrix #
############################


                    {
                    if ((dim(initialsol)[[2]] != kmax) | (dim(initialsol)[[1]] != length(kmin:kmax))) stop("\n A matrix of initial solutions for more than one cardinality must be of type no. of cardinalities x kmax")
                    if (nsol > 1) {
                       initialsol<-array(dim=c(nsol,kmax,length(kmin:kmax)),rep(t(initialsol),rep(nsol,kmax*length(kmin:kmax)))) 
                    if (implog==TRUE) {warning("\n a single initial subset necessarily produces nsol identical \n final solutions in the improve algorithm")}
                    }
}}}
                 if ((ninclude != 0) & (sum(include == rep(initialsol,rep(length(include),length(as.vector(initialsol))))) != nsol*length(kmin:kmax)*length(include))) stop("\n Not all the specified initial solutions contain the variables that are to be included") 
}


##################################
#  assigning any changed values  #
##################################

         assign("silog",silog,pos=parent.frame())          
         assign("initialsol",initialsol,pos=parent.frame())

}
validation<-function(mat, kmin, kmax, exclude, include, criterion, pcindices){

##########################################################
#  general validation of input for all search functions  #
##########################################################


####################################################################
# checking for an input matrix that must be square, of full rank,  #
#    symmetric, positive definite                                  #
####################################################################

         if (!is.matrix(mat)) 
          stop("Data is missing or is not given in matrix form\n") 
         p<-dim(mat)[2] 
         if (qr(mat)$rank != p) stop("\n The covariance/correlation matrix supplied is not of full rank\n")
         if (dim(mat)[1] != p) {
         mat <- var(mat)
         warning("Data must be given as a covariance or correlation matrix.\n It has been assumed that you wanted the covariance matrix of the \n data matrix supplied\n")}
         if ( max( abs(mat-t(mat)) ) > 1E-6 ) 
	    stop("\n The covariance/correlation matrix supplied is not symmetric\n")
  	 if ( eigen(mat,only.values=TRUE)$values[p] < 0 )
            stop("\n The covariance/correlation matrix supplied is not positive definite\n")

#################################################
# checking acceptability of criterion requested #
#################################################
         
         labelsrm<-c("RM","Rm","rm","1",1)
         labelsrv<-c("RV","Rv","rv","2",2)
         labelsgcd<-c("GCD","Gcd","gcd","3",3)
         if (sum(criterion == c(labelsrm,labelsrv,labelsgcd)) == 0) stop("criterion requested is not catered for, or has been misspecified\n")

         if (sum(criterion == labelsrm) > 0) criterio<-1
         if (sum(criterion == labelsrv) > 0) criterio<-2
         if (sum(criterion == labelsgcd) > 0) criterio<-3


         
######################################################################
#  checking consistency of the requested values of kmin, kmax        #
#  (extreme sizes of variable subsets). The validation that kmax<=p  #
#   is made later, considering the no. of excluded variables         #
######################################################################         

         if (!is.numeric(kmin) || !is.numeric(kmax)) stop("\n Arguments kmin and kmax must be numeric.\n Perhaps unnamed arguments in the wrong order?")
         if (kmax < kmin) {
          aux<-kmin
          kmin<-kmax
          kmax<-aux
          warning("the argument kmin should precede the argument kmax.\n Since the value of kmin exceeded that of kmax, they have been swapped \n")}
         if (kmin >= p) {
             kmin<-p-1
             warning("\n The value of kmin requested is equal to or exceeds the number \n of variables. It has been set at p-1. \n")
            }
         if (kmax >= p) {
                         kmax<-p-1
                         warning("\n The value of kmax requested is equal to or exceeds the number \n of variables. It has been set at p-1. \n")
            }

#############################################################
# checking for consistency of requests to exclude/include   #
#  certain variables.                                       #
#############################################################

         if (sum(duplicated(c(exclude,include))) > 0)
          stop("\n You have requested that the same variable be both included in, and excluded \n from, the subset.\n")
         nexclude<-length(exclude)
         if (nexclude !=0) {
         if (kmax >= p-nexclude) {
                                  kmax<-p-nexclude-1
                                  warning("\n Cardinalities requested are too large for the requested number of excluded \n variables, and kmax has been set at p-nexclude-1 \n")}
         exclude<-sort(exclude)}  
         exc<-c(0,exclude)
         ninclude<-length(include)
         if (ninclude !=0){
         if (kmin <= ninclude) {kmin<-ninclude+1
                                warning("\n Cardinalities requested are too small for the requested number of included \n variables, and kmin has been set to ninclude + 1 \n")}
         include<-sort(include)}
         inc<-c(0,include)
         if (kmax<kmin) stop("\n After trying to adapt to the requests for exclusion and inclusion of variables, \n kmax is now smaller than kmin. \n There must be a mistake\n")

######################
# Checking pcindices #
######################


         
# Value always assigned to esp and to non-numeric pcindices to enable passing to parent.frame and to Fortran or C++ routines, for criteria other than GCD

         esp<-FALSE
         if (criterio == 3) { 
                            if (is.null(pcindices)) stop("\n For criterion GCD, argument pcindices must be explicitely set in the leaps \n function, must be non-NULL in other search functions. \n")
                            if (is.numeric(pcindices))  {esp<-TRUE
                                    if  (sum(!(as.integer(pcindices) == pcindices)) > 0) stop("\n The PC indices must be integers.\n")
                                    if (max(pcindices)  >  p) stop("\n PCs of rank larger than the data set were requested. \n")
                                                        }
                            else {if (pcindices != "first_k")
                                    {stop("\n unrecognized value for 'pcindices' argument \n")}}}
         if ((pcindices == "first_k") || is.null(pcindices)) {pcindices<-1:kmax}

#################################
# assigning any changed values  #
#################################
         
         assign("mat",mat,pos=parent.frame())         
         assign("kmax",kmax,pos=parent.frame())         
         assign("kmin",kmin,pos=parent.frame())         
         assign("exclude",exclude,pos=parent.frame())         
         assign("include",include,pos=parent.frame())         
         assign("nexclude",nexclude,pos=parent.frame())         
         assign("ninclude",ninclude,pos=parent.frame())         
         assign("criterio",criterio,pos=parent.frame())    
         assign("pcindices",pcindices,pos=parent.frame())         
         assign("p",p,pos=parent.frame())                  
         assign("exc",exc,pos=parent.frame())         
         assign("inc",inc,pos=parent.frame())         
         assign("esp",esp,pos=parent.frame())
       }


validgenetic<-function(kmin, kmax, popsize, mutprob, exclude,
nexclude, include, ninclude, initialpop){



################################################################
# validation of input that is specific to the genetic function #
################################################################
# WARNING: Some initializations look similar to those of the   #
# anneal and improvement functions, but here input MUST be a   #
# 3-dimensional array.                                         #
################################################################
        
######################################
# validation of mutation probability #
######################################

        if ((mutprob < 0) | (mutprob > 1)) stop("\n The mutation probability parameter (mutprob) must be between 0 and 1 \n")


################################################################
# initializations when initial population has been specified;  #
# checking the nature of initialpop (and how to interpret it)  #
# and for conflicts with the exclude and include requirements  #
# (which initial the population must respect)                  #
################################################################

        if (is.null(initialpop)) {pilog <- FALSE}
        else {pilog <- TRUE  # initial population has been specified by user          
################################################
# initial solution is a vector: not acceptable #
################################################

            if (is.vector(initialpop)) {
                     {stop("\n The specified initial population must have different k-subsets \n")}
                    }
################################################
# checking for the presence of variables that  #
# are to be forcefully excluded                #
################################################

            if ((nexclude != 0) & (sum(exclude == rep(initialpop,rep(length(exclude),length(as.vector(initialpop))))) !=0)) stop("\n the specified initial population contains variables that are to be excluded \n")

#########################################################
# how to deal with various formats of input (of initial #
#  population) for a single cardinality                 #
#########################################################

            dimpop<-dim(initialpop)
            if (length(dimpop) > 3) stop("\n Can't handle arrays of more than 3 dimensions \n")
            if (kmin == kmax) {

##############################
# initial solution is array? #
##############################

              if (is.array(initialpop)) {

###########################
# initialpop is 3-d array #
###########################

              if (length(dimpop) == 3) {
                 if (dimpop[[3]] > 1) stop("\n The input array of initial population must have dimensions as \n popsize x k x 1 when a single cardinality is requested \n")
                 else initialpop<-matrix(nrow=dimpop[[1]],ncol=dimpop[[2]],initialpop)
              }
              else # initialpop is matrix: must be of form popsize x k
 
                 {if (dimpop[[2]] != kmax) stop("\n Input matrix of initial solutions must have as many columns as variables in the requested subset \n")
                  else 
                      {if  (dimpop[[1]] != popsize) {
                        stop("\n The number of initial solutions can only be 1 or popsize (number of final solutions requested) \n")}}}
              }}

##########################################################
# how to deal with various formats of input (initialpop) #
# if more than one cardinality is requested              #
##########################################################

             else  #(if kmax > kmin), i.e., more than one cardinality requested

#########################################
# initial population must be 3-d array  #
#########################################


              {if (length(dimpop) < 3) stop("\n There must be initial populations for all cardinalities requested \n")
               else 
               {if ((dimpop[[3]] != length(kmin:kmax)) | (dimpop[[2]] != kmax) | (dimpop[[1]] != popsize)) stop("\n The input array of initial solutions must have dimensions as \n popsize x kmax x no. of different cardinalities requested \n")
                  }
                 if ((ninclude != 0) & (sum(include == rep(initialpop,rep(length(include),length(as.vector(initialpop))))) != popsize*length(kmin:kmax)*length(include))) stop("\n Not all the specified initial solutions contain the variables that are to be included \n") 
}}

##################################
#  assigning any changed values  #
##################################

         assign("pilog",pilog,pos=parent.frame())          
         assign("initialpop",initialpop,pos=parent.frame())

}


.First.lib <- function(lib, pkg)
{
    library.dynam("subselect", pkg, lib)
}
