select.cp.aic.bic.fun_function(x,y,keepit=NULL)
{
  ## this program works for a search of models with main effects and no interactions
  ## command is:  select.cp.aic.bic.fun(cbind(x1,x2,x3,x4,x5),y)
  ## the output is ordered according to BIC.
  n_length(y)
  fullmod <- lm(y~x)
  s2full <- (summary(fullmod)$sigma)^2
  sstot_crossprod(y-mean(y))
  nullmod <- lm(y~1)
  null.p <- 1
  null.cp <- 1
  null.s2 <- summary(nullmod)$sigma^2
  null.bic <- n*log(null.s2)+null.p*log(n)
  null.aic <- n*log(null.s2)+2
  null.r2 <- summary(nullmod)$r.squared
  null.r2adj <- 1-null.s2/(sstot/(n-1))
  r2info_leaps(x,y,method="r2",keep=keepit)
  mstot <- sstot/n
  p_r2info$size
  r2 <- r2info$r2/100
  r2adj <- (r2-(1-r2)*(p-1)/(n-p))
  ssres_(1-r2info$r2/100)*sstot
  sigma2hat_ssres/(n-p)
  Cp <- n-p*sigma2hat/s2full+(2*p-n)
  AIC_n*log(sigma2hat)+2*p
  
  BIC_n*log(sigma2hat)+p*log(n)
  ord <- order(BIC)
  output_data.frame(model=r2info$label[ord],p=r2info$size[ord],Cp=round(Cp,3)[ord],AIC=round(AIC,3)[ord],BIC=round(BIC,3)[ord],R2=round(r2info$r2/100,4)[ord],R2adj=round(r2adj,4)[ord])
  print("Subset selection criteria results: [[1]] Model Names, [[2]] Criteria")
  nullmod <- data.frame("NULL",round(null.p,0),round(null.cp,3),round(null.aic,3),round(null.bic,3),round(null.r2,4),round(null.r2adj,4))
  names(nullmod) <- names(output)
printout_rbind(output,nullmod)
return(criteria=printout[,2:7],model.names=printout[,1])
}