# ---------------------------------------------------------------------
#
#   Use Latin Hypercube sampling to populate the GAMS parameter
#
#   p_doe(draws,factors)
#
#   and output it to GDX
#
#   where draws   are the design draws
#         factors are the orthogonal factor
#
#   # of draws and set factor are read from GDX
#
# ---------------------------------------------------------------------
sink(stdout(), type="message")
#options(echo=TRUE)


#install.packages("d:\r\R-2.15.1\library\mc2d_0.1-13.zip",repos=NULL);
#install.packages("d:\r\R-2.15.1\library\mvtnorm_0.9-9992.zip",repos=NULL);
#install.packages("d:\r\R-2.15.1\library\lhs_0.10.zip",repos=NULL);
#install.packages("d:\\temp\\gclus_1.3.1.zip",repos=NULL);
# install.packages("t:\\britz\\gdxrrw_0.0-2.zip",repos=NULL);
# install.packages("D:\\temp\\gdxrrw_1.0.2.zip", repos = NULL, type="source");
 library(lhs);
 library(gdxrrw);
 igdx("N:/soft/gams24.7new/24.7");
 library(mc2d);
 library(mvtnorm);
 library(Matrix);
 library(gclus);
 library(reshape2);
#useCorr    <- Sys.getenv("useCorr")
#useColors  <- Sys.getenv("useColors")
#inputFile  <- Sys.getenv("inputFile");
#plotFile   <- Sys.getenv("plotFile")
#outputFile <- Sys.getenv("outputFile");
#maxRunTime <- as.numeric(Sys.getenv("maxRunTime"))


rm(list=ls(all=TRUE));
incFile <- commandArgs(trailingOnly = TRUE);

#incFile <- "d:\\temp\\toLHS.r";

if ( is.na(incFile[1]) )
  incFile[1] <- "d:\\temp\\tolhs.r"

source(incFile[1], echo=TRUE);


 if ( exists("inputFile") ){


 if (!file.exists(inputFile)) stop(paste("inputFile ",inputFile," file doesn't exist! Check the code!"));


#
# --- determine GDX file with input data
#

#
# --- number of draws as requested by GAMS
#
      n = rgdx.scalar(inputFile,"p_n");

#
# --- name of the set which comprises the scenario name
#
      scen_name = rgdx.set(inputFile,"scen_name",compress="true");
      scenName  <- levels(scen_name[scen_name$i == 'name',2])
#
# --- name of the set which comprises the factors
#
      factor_name = rgdx.set(inputFile,"factor_name",compress="true");
      set_to_load <- levels(factor_name[factor_name$i == 'name',2])
#
# --- load that set
#
      factors = rgdx.set(inputFile,set_to_load,compress="true");

      k = nrow(factors);
 } else {
  k = length(factors);
 }

 print(paste("factor ",k));
#
# --- the number of factors is equal to the number of set elements
#
 LHSType <- "improvedLHS";



 if ( LHSType == "optimumLHS" ){
    out <- optimumLHS(n,k,2,0.01)
    reportDraws = 1;
 } else {
    out <- improvedLHS(n,k)
    reportDraws = 1;
 }

 print(useCorr);
 if ( useCorr == "true" )
 {

    begTime <- as.numeric(Sys.time(),units="seconds");
    runTime <- 0;
#
#   --- load correlation matrix from GAMS
#
    t <- rgdx.param(inputFile,"p_cor",names=c("f1","f2"),compress="true");
    t
    t<-acast(t, f1~f2, value.var="value")
    t<-as.matrix(t);

#
#   --- find nearest positive definite matrix
#

    t1<-nearPD(t);

    t <- as.matrix(t1$mat);

    bestFit <- 10;
    iDraw <- 0;


    while( runTime < maxRunTime ){

       iDraw <- iDraw + 1;

       if ( LHSType == "optimumLHS" ) {
          out1 <- optimumLHS(n,k,2,0.01);

          print("shit");

       } else {
          out1 <- improvedLHS(n,k);
       }

#
# --- use cornode to apply Iman & conover 1982 to impose correlation
#     t on the LHS matrix out
#
       out1  <- cornode(out1,target=t)
       c <- cor(out1);

       fit = 0;
       for ( i in 1:k )
            for ( j in 1:k )
                if ( fit < bestFit )

       if ( fit < bestFit )
       {
          out     <-out1;
          bestFit <-fit;

          meanDev = sqrt(fit/k)*100;
       };

       if ( iDraw %% reportDraws == 0){
          curTime <- as.numeric(Sys.time(),units="seconds");
          runTime <- curTime - begTime;
          print(paste(" draw :",iDraw," runTime ",round(runTime)," of ",maxRunTime,"seconds, mean sqrt of squared diff between given corr and best draw: ",round(meanDev,2),"%"));
       }

# --- meanDev threshold: if mean % dev. of randomized correlation matrix is lower than "1%" the draw is taken for the sample and repetition of sample draws is stopped

       if ( meanDev < 1 ){
            print(paste(" draw :",iDraw," runTime ",round(runTime)," of ",maxRunTime,"seconds, mean sqrt of squared diff between given corr and best draw: ",round(meanDev,2),"%"));
            runTime <- maxRunTime;
       }

    }
 }

 c <- cor(out);
 c;

 test <- summary(out);

 p <- c(0.05,0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95);
 t(apply(out, 2, quantile, probs = p))

#
# --- convert matrix to data frame
#
 out <- data.frame(out)
#
# --- use the set elements as names for the data frame
#
 if ( exists("inputFile") ){
   colnames(out) <- factors$i
 } else {
   colnames(out) <- factors
 }

#
# --- scatter plot parts