<repo_name>Murilodsv/dssat_sscan_calib<file_sep>/bin/check_par_nonlogic.R
#--- Check parameters phyiological consistency 

check_par_nonlogic = function(calib){
  
  #--- check parameter physilogical logic
  out = F
  
  if(calib$pvalue[calib$Par_orig=="CHUPIBASE"] < calib$pvalue[calib$Par_orig=="TTPLNTEM"]){
    message("")
    message("Warning: CHUPIBASE lower than TTPLNTEM (Objective increased)")
    out = T
  }
  
  if(calib$pvalue[calib$Par_orig=="CHUPIBASE"] < calib$pvalue[calib$Par_orig=="TTRATNEM"]){
    message("")
    message("Warning: CHUPIBASE lower than TTRATNEM (Objective increased)")
    out = T
  }
  
  if(calib$pvalue[calib$Par_orig=="PI2"]   < calib$pvalue[calib$Par_orig=="PI1"]){
    message("")
    message("Warning: PI2 lower than PI1 (Objective increased)")
    out = T
  }
  
  if(calib$pvalue[calib$Par_orig=="PI2"]   < calib$pvalue[calib$Par_orig=="PI1"]){
    message("")
    message("Warning: PI2 lower than PI1 (Objective increased)")
    out = T
  }
  
  return(out)
  
}<file_sep>/bin/f_dssat_sccan_calib.R
dssat_sccan_calib = function(svalue){
  
  #------------------------------------------#
  #--- DSSAT-CANEGRO Calibration function ---#
  #------------------------------------------#
  
  #--- Goals:
  #------- 1) Update parameters file; 
  #------- 2) Run the model; 
  #------- 3) Compute performance based on observed data
  
  #--- reference: Functional, structural and agrohydrological sugarcane crop modelling: towards a simulation platform for Brazilian farming systems (PhD Dissertation, 2018)
  #--- contact: <NAME> (<EMAIL>)
  
  #--- Track time
  start_time  = Sys.time()
  
  p             = 4 # parameter precision (#decimals)

  #--- Convert scaled to parameter values (feature scales method)
  calib$pvalue = calib$Calib_values
  calib$pvalue[calib$Calibrate] = svalue * (calib$Calib_range_max[calib$Calibrate] - calib$Calib_range_min[calib$Calibrate]) + calib$Calib_range_min[calib$Calibrate]
  
  #--- cultivar ID to follow-up optimization
  cat("\014")
  message(paste("Calibrating cultivar: ",calib_id,sep=""))
  message("")
  
  #--- write new parameters on screen
  message("Using the following parameters:")
  for(i in unique(calib$Par_orig[calib$Calibrate])){ 
    
    plab = i
    sizep = max(nchar(as.character(calib$Par_orig)))
    while(nchar(plab) < sizep){
      plab = paste(plab," ",sep = "")
    }
    
    message(paste(plab," = ",sprintf(paste("%.",p,"f",sep=""),calib$pvalue[calib$Par_orig==i]),sep=""))
  }
  
  #--- check for physiological meaning
  penalize = check_par_nonlogic(calib)
  
  #--- Check if estimated parameters are within physiological-boundaries 
  if(uselimits){
    
    #--- check constraints
    for(i in 1:length(svalue)) {
      if(svalue[i] > 1 | svalue[i] < 0){
        
        #--- msg
        message(paste("Warning: Parameter ",calib$Par_orig[calib$Calibrate][i], " is out of min and max range (Objective Penalized)",sep = ""))
        
        #--- Increase objective function 
        penalize = T
      }
    }
  }
  
  
  #--- Read model parameters MASTER files
  par_cul_file  = readLines(paste(wd,"/templates/SCCAN047_M.CUL",sep = ""))
  par_eco_file  = readLines(paste(wd,"/templates/SCCAN047_M.ECO",sep = ""))
  rep_line      = 4 # line of file where the replacement will be done
  
  #--- re-build cultivar file (.CUL)
  l_ftype = c(".CUL",".ECO")
  
  for(ftype in l_ftype){
    
    #---Read model parameters MASTER file
    p_file = readLines(paste(wd,"/templates/SCCAN047_M",ftype,sep = ""))
    
    #--- list of parameters
    l_rp = unique(calib$Rep_ID[calib$File == ftype])
    
    for(rp in l_rp) {
      
      replace = sprintf(paste("%.",p,"f",sep=""),calib$pvalue[calib$Rep_ID==rp & calib$File==ftype])
      
      #--- check whether new parameter match size
      m_size = calib$Par_size[calib$Rep_ID==rp & calib$File==ftype] - nchar(replace)
      
      if(m_size > 0){
        #--- add spaces to match size
        while(m_size > 0){
          replace = paste(" ",replace,sep="")
          m_size = calib$Par_size[calib$Rep_ID==rp & calib$File==ftype] - nchar(replace)
        }
        
      }else if(m_size < 0){
        #--- remove precision to match size
        message(paste("Warning: Precision of parameter ",calib$Par_orig[calib$Rep_ID==rp & calib$File==ftype]," reduced to match file.CUL fmt",sep=""))
        
        red_p = p
        while(m_size < 0){
          red_p = red_p - 1
          
          if(red_p < 0){stop(paste("Parameter ",calib$Par_orig[calib$Rep_ID==rp & calib$File==ftype], " is too high for file.CUL fmt (try to reduce maximun range)",sep=""))}
          
          replace = sprintf(gsub(p,red_p,paste("%.",p,"f",sep="")),calib$pvalue[calib$Rep_ID==rp & calib$File==ftype])
          m_size = calib$Par_size[calib$Rep_ID==rp & calib$File==ftype] - nchar(replace)
        }
      }
      
      
      p_file[rep_line] = gsub(rp,
                              replace
                              ,p_file[rep_line])
      
    }
    
    #--- write parameter file 
    write(p_file,file =paste("C:/DSSAT",ds_v,"/Genotype/",parnm,ftype,sep = ""))
  }
  
  #--- set wd to run
  setwd(paste("C:/DSSAT",ds_v,"/",crop,"/",sep = ""))
  
  #--- write paramters used on the screen
  message("")
  message("Running DSSAT-Canegro...")
  
  #-----------------#
  #--- Run DSSAT ---#
  #-----------------#
  
  #--- Call DSSAT047.exe and run X files list within DSSBatch.v47
  system(paste("C:/DSSAT",ds_v,"/DSCSM0",ds_v,".EXE SCCAN0",ds_v," B ",paste("DSSBatch.v",ds_v,sep=""),sep=""),show.output.on.console = pdssat)
  
  #--- Read simulated data
  plant_lines = readLines("PlantGro.OUT")
  
  #--- Note: writing file is required to speed up! (for some reason is faster than reading directly from plant_lines variable)
  write.table(plant_lines[substr(plant_lines,2,3)=="19" | substr(plant_lines,2,3)=="20"],
              file = paste("PlantGro_",calib_id,".OUT",sep=""),
              row.names = F, col.names = F,quote = F)
  plant = read.table(file = paste("PlantGro_",calib_id,".OUT",sep=""))                   #Read numeric lines as data.frame
  
  #--- Columns name accordingly to DSSAT output name
  colnames(plant) = pgro_names
  
  #--- observed data
  obs = obs_df$obs
  
  #-- simulated data
  sim = plant[plant$dap %in% obs_df$dap,sscan_out]
  
  
  
  if(plotperf){
    
    png(paste(wd,"/results/perf_",calib_id,".png",sep=""),
        units="in", 
        width=24, 
        height=12, 
        pointsize=24, 
        res=300)
    
    par(mfrow=c(1,2), mar = c(4.5, 4.5, 0.5, 0.5), oma = c(0, 0, 0, 0))
    objective = mperf(sim,obs,sscan_out,plotperf,outidx)[,2]
    
    #--- write performance
    model_perf = mperf(sim,obs,sscan_out,F)
    
    perf_df = data.frame(cv = calib_id,model_perf)
    write.csv(perf_df, file = paste(wd,"/results/perf_",calib_id,".csv",sep=""),row.names = F)
    
    plot(plant[,sscan_out]~plant[,"dap"],
         type = "l",
         xlab = "DAP",
         ylab = sscan_out)
    points(obs_df$obs~obs_df$dap)
    
    legend("topleft",
           inset   = 0.02,
           legend  = c("Simulated","Observed"),
           col     = c("black","black"),
           lt      = c(1,0),
           pch     = c(NA,1),
           bg      = "grey",
           cex     = 1.0,
           box.lty = 1)
    
    dev.off()
    
  }else{
    #--- mperf compute several indexes, RMSE is outputed
    objective = mperf(sim,obs,sscan_out,plotperf,outidx)[,2]

  }
  
  #--- print this run objective
  message("")
  message(paste(outidx," is: ",objective,sep=""))
  
  if(plotdev){
    
    it_before = read.csv(file = paste(wd,"/results/optim_dev.csv",sep=""))
    
    obj_df = data.frame(n = (max(it_before$n)+1),obj = objective,inbounds = T)
    
    for(i in 1:length(svalue)) {
      df = data.frame(v = svalue[i] * (calib$Calib_range_max[calib$Calibrate][i] - calib$Calib_range_min[calib$Calibrate][i]) + calib$Calib_range_min[calib$Calibrate][i])
      colnames(df) = calib$Par_orig[calib$Calibrate][i]
      obj_df = cbind(obj_df,df)
      
      #--- Use constrained limits from par_set max and min?
      if(uselimits){
        
        if(penalize){
          #--- Flag values that are out of user-boundaries
          obj_df$inbounds = F
        }
      }
    }
    
    if(length(it_before$n) == 1){
      it_before = obj_df
      it_after = rbind(it_before,obj_df)
    }else{
      
      it_after = rbind(it_before,obj_df)
      if(length(dev.list())>0){dev.off()}
    }
    
    #--- plot on screen
    plot(it_after$obj~it_after$n, type = "l",ylab = outidx,xlab = "Number of iterations", ylim = c(0,max(it_after$obj)))
    lines(c(min(it_after$obj),min(it_after$obj))~c(-1000,max(it_after$n)*1000), lty = 3,col = "red")
    
    write.csv(it_after,file = paste(wd,"/results/optim_dev.csv",sep=""),row.names = F)
    
  }
  
  #--- Penalize objetive
  if(uselimits){
    if(penalize){
      if(plotdev){
        objective = max(it_before$obj)
      }else{
        penalty = 1000
        objective = objective * penalty
      }
    }
  }
  
  end_time    = Sys.time()
  message(paste("Elapsed time: ",round(end_time - start_time,digits = 2),"sec"))
  message("--------------------------------------------")
  
  
  objective
  
}


<file_sep>/ReadMe.md
---------------------------------
# DSSAT/CANEGRO Calibration
<NAME> (Jun-2018)
----------------------------------

# Main Goal:
Calibrate .CUL and .ECO crop parameters of DSSAT/CANEGRO (PBM)

# Methods:
It uses the general purpose optimization function "optim()" embedded in R environment. User can set up the optimization method, objective function, observed data, and parameters to be calibrated. Charts and tables are provided as outputs with calibrated parameters in folder 'results'.

![alt text](https://github.com/Murilodsv/dssat_sscan_calib/blob/master/framework.png)

# Warning:
This subroutine will replace your .CUL and .ECO, so please make sure you backed up your original coefficients before using it.  

# Example Run:
An example is already set in the table dssat_canegro_calib_par.csv and file dssat_sccan_calib.R
To run it make sure to copy the files SCGO0001.SCX and GOGO.WTH into the DSSAT47/Sugarcane folder.

Then, open dssat_sccan_calib.R and follow the comented instructions.

# How to Use:
1) Set the dssat_canegro_calib_par.csv file with initial parameters and boundaries:
Specify in the column 'Calibrate' which parameter will be calibrated by writing T or F
In Columns 'Init_values' 'Calib_range_min'	'Calib_range_max' specify the initial parameters values, the minimum and maximum range.

2) Open dssat_sccan_calib.R and follow the comented instructions

# Sample Results:

Optimization progress of RMSE

![opmin_progress](https://github.com/Murilodsv/dssat_sscan_calib/blob/master/results/optimization_CTC2.png)

Best fit results
![bf_res](https://github.com/Murilodsv/dssat_sscan_calib/blob/master/results/perf_CTC2.png)
<file_sep>/dssat_sccan_calib.R
#---------------------------------#
#--- DSSAT CANEGRO Calibration ---#
#---------------------------------#

#--- Goal: Calibrate crop parameters of DSSAT/CANEGRO 
#---       (this is not GLUE or GenCalc)

#--- History
#--- Jun-2018: Created (<NAME>)
#--- Jan-2020: Updated for new CANEGRO version (<NAME>)

#----------------------------#
#--- Running this example ---#
#----------------------------#

#--- Before start copy the below files to DSSATv47/Sugarcane folder:
#--- 'repository'/db/SCGO0001.SCX     (X file for this example)
#--- 'repository'/db/GOGO.WTH         (Weather file for this example)

#--- This example will optimize the tillering coefficients of DSSAT/CANEGRO for
#--- different sugarcane varieties using the variety Nco376 as base
#--- Results will be placed in the 'results' of this folder

#--------------------#
#--- Script Setup ---#
#--------------------#
wd        = "C:/Murilo/dssat_sscan_calib" # Working directory (repository)            
ds_v      = 47                            # DSSAT version
crop      = "Sugarcane"                   # Crop name
xfile     = "SCGO0001.SCX"                # X file
parnm     = "SCCAN047"                    # Model name SCCAN047 = CANEGRO
pdssat    = F                             # Print DSSAT echo while running
savepng   = T                             # save optimization dev to png file?

#--------------------------#
#--- Optimization Setup ---#
#--------------------------#
op_reltol     = 0.01          # Relative convergence tolerance
method.opt    = "Nelder-Mead" # More options can be set for other methods (see ?optim())

#--- Number of optimization repetitions 
#--- A new set of initial parameters is randomly pick for each repetition  
#--- The higher this number, the higher is the probability to achieve global minimum objective
nopt          =  5            # We set 5 for the sake of this example, but a higher must be used for serious calibration

#--- Observed data used for calibration
#--- Note: here we are using one variable but you can add as many as needed
#---  As long as the model outputs that...
used.data     = "Tillering_n_m-2"

#--- model output to be compared with observation
sscan_out     = "t.ad" # is possible to add as many outputs as needed with minor change in the code

#--- Statistical index used as objective function by the optimization (for different indexes see mperf())
outidx        = "rmse" # Using Root Mean Squared Error (RMSE). See mperf() function for more options

#--- Load Functions (~/bin/) 
invisible(sapply(list.files(path = paste0(wd,"/bin/"),full.names = T),
                               function(x) source(x)))

#--- Create results folder
dir.create(paste0(wd,"/results"), showWarnings = FALSE)

#--- Reading plant observations and PlantGro Header
obs_raw   = read.csv(paste0(wd,"/db/gogo_field_data.csv"))
pgro_head = read.csv(paste0(wd,"/db/PlantGro_Head.csv"))

#--- PlantGro header
pgro_names  = pgro_head$R_head

#--- prepare batch call
bfile = readLines(paste(wd,"/templates/DSSBatch_Master.v47",sep=""))
bfile[4] = gsub("<calib_xfile>",xfile,bfile[4])

#--- write batch in Crop folder
write(bfile,file = paste("C:/DSSAT",ds_v,"/",crop,"/","DSSBatch.v",ds_v,sep = ""))

#--- Read parameters set up
par_set   = read.csv(paste(wd,"/dssat_canegro_calib_par.csv",sep=""))
calib     = par_set

l_cv = unique(obs_raw$Cultivar)

for(cv in l_cv){

  #--- logical settings
  plotperf = F # plot sim x obs?
  plotdev  = T # follow-up optimization?
  uselimits= T # use min and max boudaries to drive optimization? (in set_par max and min) 
  
  #--- current cultivar
  calib_id  = cv

  message(paste("Start of optimization for ",cv,sep=""))
  
  #--- Feature scaling (0-1)
  calib$svalue = (calib$Init_values - calib$Calib_range_min) /  (calib$Calib_range_max - calib$Calib_range_min)
  
  #--- Parameters to be calibrated
  svalue = calib$svalue[calib$Calibrate]
  
  #--- number of parameters
  npar = length(svalue)
  
  #--- wipe optimization file
  obj_df = data.frame(n = 0, obj = 0,inbounds = T)
  
  for(i in 1:length(svalue)) {
    df = data.frame(v = calib$Init_values[calib$Calibrate][i])
    colnames(df) = calib$Par_orig[calib$Calibrate][i]
    obj_df = cbind(obj_df,df)
  }
  
  #--- write a temporary file with optimization progression
  write.csv(obj_df,file = paste(wd,"/results/optim_dev.csv",sep=""),row.names = F)
  
  #--- read observed data for cultivar (cv)
  obs_df      = data.frame(cv  = cv,
                           dap = obs_raw$DAP[obs_raw$Cultivar==calib_id & obs_raw$Type==used.data],
                           obs = obs_raw$Value[obs_raw$Cultivar==calib_id & obs_raw$Type==used.data])
  
  
for(i in 1:nopt){
  
  #--------------------#
  #--- Optimization ---#
  #--------------------#
  
  #--- Optimize 
  optim(svalue,dssat_sccan_calib,control = list(reltol = op_reltol),
        method = method.opt)
  
  #--- restart iniial conditions (try to fall on global minimum)
  new_val = rnorm(1000,0.5,0.25)
  new_val = new_val[new_val>=0 & new_val<=1]
  
  #--- pick randomized points within the distribution
  new_val = new_val[abs(rnorm(length(svalue),500,100))]
  
  #--- check if the number of paramters match with requirements
  if(length(new_val) < npar){
    ntry = 1
    while(length(new_val) == npar){
      new_val = new_val[abs(rnorm(length(svalue),500,100))]
      ntry = ntry + 1
      if(ntry > 1000){stop("Failed to find new inital conditions, please review the initial parameters setup.")}
    }
  }
  
  #--- Reinit
  if(i != nopt){svalue = new_val}
  
}
  
  message(paste("End of optimization for ",cv,sep=""))
  
  #--- write optimization parameters
  opt_par = read.csv(paste(wd,"/results/optim_dev.csv",sep=""))
  write.csv(opt_par,file = paste(wd,"/results/optim_dev_",cv,".csv",sep=""),row.names = F)
  
  #--- save to png file (there is room for using ggplot here)
  if(savepng){
    png(paste(wd,"/results/optimization_",calib_id,".png",sep=""),
        units="in", 
        width=24, 
        height=12, 
        pointsize=24, 
        res=300)
    plot(opt_par$obj~opt_par$n, type = "l",ylab = outidx,xlab = "Number of iterations", ylim = c(0,max(opt_par$obj)))
    lines(c(min(opt_par$obj),min(opt_par$obj))~c(-1000,max(opt_par$n)*1000), lty = 3,col = "red")
    
    dev.off()
    
  }
  
  
  #--- Best set of parameters
  svalue_df =  opt_par[opt_par$obj==min(opt_par$obj[opt_par$inbounds]),as.character(calib$Par_orig[calib$Calibrate])]
  
  if(length(svalue_df[,as.character(calib$Par_orig[calib$Calibrate][1])]) > 1){
    #--- pick up the median minimun values
    svalue_df = sapply(as.character(calib$Par_orig[calib$Calibrate]),function(x) median(svalue_df[,x]))
  }
  
  svalue = svalue_df
  
  #--- scale paramters
  svalue = (svalue - calib$Calib_range_min[calib$Calibrate]) / (calib$Calib_range_max[calib$Calibrate] - calib$Calib_range_min[calib$Calibrate])
  
  #--- check best parameters performance
  plotperf = T # plot sim x obs?
  plotdev  = F # follow-up optimization?
  dssat_sccan_calib(t(svalue))

}

<file_sep>/bin/mperf.R
# Model Performance -------------------------------------------------------

mperf     = function(sim,obs,vnam,dchart,outidx){
  
  #--------------------------------------------------#
  #------------- Performance function ---------------#
  #--- Compute statistical indexes of performance ---#
  #--------------------------------------------------#
  #
  # Decription:
  # sim     - Simulated values          [Real]
  # obs     - Observed values           [Real]
  # vnam    - Name of variable          [String]
  # dchart  - Display Chart?            [T or F]
  # outidx  - Output peformance indexes [List]
  # 
  # <NAME>
  # source: https://github.com/Murilodsv/R-scripts/blob/master/mperf.r
  #
  # Literature: 
  # <NAME>., <NAME>., <NAME>., & <NAME>. (2006). 
  # Working with dynamic crop models: evaluation, analysis, 
  # parameterization, and applications. Elsevier.
  #--------------------------------------------------#
  
  if(missing(sim)){stop("Missing sim argument")}
  if(missing(obs)){stop("Missing obs argument")}
  if(missing(dchart)){dchart = T}
  if(missing(outidx)){outidx = c("bias","mse","rmse","mae","rrmse","rmae","ef","r","r2","d","cc","a","b","mi.sim","sd.sim","cv.sim","mi.obs","sd.obs","cv.obs","n")}
  if(missing(vnam)){
    warning("Missing vnam argument: vnam set to none.")
    vnam = ""
  }
  
  #--- all outputs
  if(outidx[1] == "all"){outidx = c("bias","mse","rmse","mae","rrmse","rmae","ef","r","r2","d","cc","a","b","mi.sim","sd.sim","cv.sim","mi.obs","sd.obs","cv.obs","n")}
  
  #--- Check Input data
  sim = as.numeric(sim)
  obs = as.numeric(obs)
  
  if(length(sim) != length(obs)){stop("Vector length of Simulated and Observed do not match.")}
  if(length(sim[is.na(sim)]) > 0 |  length(obs[is.na(obs)]) > 0){
    
    #--- remove NA values
    df_so = na.omit(data.frame(sim,obs))
    sim = df_so$sim
    obs = df_so$obs
    warning("NA values ignored in Simulated or Observed data")
  }
  
  #--- Statistical indexes
  n     = length(obs)
  mi.sim= mean(sim)
  mi.obs= mean(obs)
  fit   = lm(sim~obs)
  bias  = (1/n) * sum(sim-obs)
  mse   = (1/n) * sum((sim-obs)^2)
  rmse  = sqrt(mse)
  mae   = (1/n) * sum(abs(sim-obs))
  rrmse = rmse / mean(obs)
  rmae  = (1/length(obs[obs>0])) * sum(abs(sim[obs>0]-obs[obs>0])/abs(obs[obs>0]))
  ef    = 1 - (sum((sim-obs)^2) / sum((obs-mean(obs))^2))
  r     = sum((obs-mean(obs))*(sim-mean(sim)))/sqrt(sum((obs-mean(obs))^2)*sum((sim-mean(sim))^2))
  r2    = r^2
  d     = 1 - (sum((sim-obs)^2) / sum((abs(sim-mean(obs))+abs(obs-mean(obs)))^2))
  if(length(unique(sim)) > 1){
    a     = fit$coefficients[1]
    b     = fit$coefficients[2]
    if(any(is.na(c(a,b)))){
      outidx = outidx[outidx != "a" & outidx != "b"]
    }
  }else{
    outidx = outidx[outidx != "a" & outidx != "b"]
    a = NA
    b = NA
  }
  
  sigma.obs.sim = (1 / n) * sum((obs - mi.obs) * (sim - mi.sim))
  sigma.obs.2   = (1 / n) * sum((obs - mi.obs) ^ 2)
  sigma.sim.2   = (1 / n) * sum((sim - mi.sim) ^ 2)
  cc    = (2 * sigma.obs.sim) / (sigma.obs.2 + sigma.sim.2 + (mi.sim - mi.obs)^2)
  
  sd.obs = sd(obs)
  sd.sim = sd(sim)
  cv.obs = sd.obs / mi.obs
  cv.sim = sd.sim / mi.sim
  
  if(dchart){
    #--- Chart Sim ~ Obs
    varlab = vnam 
    
    mindt = min(obs,sim)
    maxdt = max(obs,sim)
    #--- Ploting limits 
    pllim = c(mindt-0.1*(maxdt-mindt),maxdt+0.1*(maxdt-mindt))
    xx = seq(min(obs),max(obs),length = (max(obs)-min(obs))*1000)
    z = summary(fit)
    
    plot(sim~obs,
         ylab = paste("Sim - ",varlab,sep = ""),
         xlab = paste("Obs - ",varlab,sep = ""),
         ylim = pllim,
         xlim = pllim)
    
    lines(xx, predict(fit, data.frame(obs=xx)),
          col = "black",
          lty = 1,
          lwd = 1.5)
    
    l11 = seq(pllim[1]-0.5*(maxdt-mindt), pllim[2] + 0.5 * (maxdt-mindt),length = 1000)
    
    lines(l11*1~l11,
          col = "red",
          lty = 2,
          lwd = 1.5)
  }
  
  perf = data.frame(model = vnam,
                    bias,
                    mse,
                    rmse,
                    mae,
                    rrmse,
                    rmae,
                    ef,
                    r,
                    r2,
                    d,
                    cc,
                    a,
                    b,
                    mi.sim,
                    sd.sim,
                    cv.sim,
                    mi.obs,
                    sd.obs,
                    cv.obs,
                    n)
  
  rownames(perf) = c()
  
  return(perf[,c("model",outidx)])
  
}
