* Optimization Code for NOLH: Mixed integer program for constructing nearly orthogonal Latin Hypercube designs
* Copyright (c) 2008 Alejandro S. Hernandez  

* This optimization code within a heuristic is distributed in the hope that it will be useful 
* for developing new nearly orthogonal Latin hypercube designs, but WITHOUT ANY WARRANTY; 
* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  

* This code may be redistributed and/or modified under the terms of the GNU Lesser General Public License; 
* see the GNU Lesser General Public License for more details.  
* The author retains all intellectual property rights in relation to this code.

* Hernandez, A.S., T.W. Lucas, and M. Carlyle. 2011. Enabling nearly orthogonal Latin hypercube construction for any non-saturated run-variable combination. Under review.


$TITLE ITERATIVE LATIN HYPERCUBE GENERATION 20JUL12

$INLINECOM { }
 OPTIONS
   SOLPRINT =     OFF,
   LIMCOL   =       8,
   LIMROW   =      8,
   RESLIM   =     600, {max seconds}
   ITERLIM  =    99999, {max pivots}
   OPTCR    =    0.001,  {relative integrality tolerance}
   LP       =   CPLEX,
   RMIP     =   CPLEX,
   MIP      =   CPLEX
 ;
*{OSL, CPLEX, XA, ... }
 SETS
   k       factors /k1*k64/
   n       runs    /n1*n65/
 ;

 alias
  (i,n)
  (kp,k)
 ;

 TABLE
   columns(n,k)

*Identify file that contains the initial design matrix.
$ONDELIM

$INCLUDE c:\##.csv

$OFFDELIM
 ;

 SCALARS
   xbar
   sumx
   offset
   currcol
   numer
   denom
   counter
   maxcount
   maxrms
   limitrho
   found
   rms
 ;
*Compute constants.
 xbar = (1+card(n))/2;
 sumx = (card(n)*(card(n)+1))/2;
 offset = -2*sumx*xbar + card(n)*xbar*xbar;
 denom = card(n)**3 - card(n);
 maxcount = 10*(card(n)*card(k));

*Create initial correlation table.
 PARAMETER
 rhotab(k,kp);
 loop(k,
   loop(kp $(ord(kp) > ord(k)),
   numer = sum(n, sqr(columns(n,k) - columns(n,kp)));
   rhotab(k,kp) = 1 - 6*(numer/denom);
   rhotab(kp,k) = rhotab(k,kp);
  );
 );

*Select the next column to optimize, using mean square.
 maxrms = 0;
 loop(k ,
  rms = sqrt((sum(kp, sqr(rhotab(k,kp)))) / (card(k) - 1) );
  if (rms >= maxrms,
   maxrms = rms;
   currcol = ord(k);
   );
  );

 VARIABLE
   Z
   RHO
 ;

 BINARY VARIABLES
   X(n,i)   run n is set to level i
 ;

 EQUATIONS
   OBJECTIVE
   ROWS(n)
   NUMBERS(n)
   COV_BOUND_PLUS(k)
   COV_BOUND_MINUS(k)
 ;

 OBJECTIVE..
   Z =e= RHO
 ;

 ROWS(n)..
   SUM(i,X(n,i)) =e= 1
 ;

 NUMBERS(n)..
   SUM(i,X(i,n)) =e= 1
 ;

 COV_BOUND_PLUS(k)$(ord(k) <> currcol)..
   RHO =g= SUM(n, columns(n,k)*SUM(i,ORD(i)*X(n,i))) + offset
 ;

 COV_BOUND_MINUS(k)$(ord(k) <> currcol)..
   RHO =g= -SUM(n, columns(n,k)*SUM(i,ORD(i)*X(n,i))) - offset
 ;

 MODEL ONECOL /
   OBJECTIVE
   ROWS
   NUMBERS
   COV_BOUND_PLUS
   COV_BOUND_MINUS
 /;
*ONECOL.optfile = 1;

*Provide initial values for the variable X.
 LOOP(k$(ORD(k)=currcol),
   LOOP(n,
     LOOP(i,
       IF(ORD(i)=columns(n,k),
         X.l(n,i)=1;
       ELSE
         X.l(n,i)=0;
       );
     );
   );
 );

 SOLVE ONECOL USING MIP MINIMIZING Z;
*Update the columns(n,k) table with the optimized values of X.
 loop(k$(ord(k) = currcol),
  loop(n,
   loop(i,
    if(X.l(n,i) = 1,
     columns(n,k) = ord(i);
    );
   );
  );
 );
*Update correlation table based on optimized column.
  loop(k $(ord(k) = currcol),
   loop(kp $(ord(kp) <> ord(k)),
   numer = sum(n, sqr(columns(n,k) - columns(n,kp)));
   rhotab(k,kp) = 1 - 6*(numer/denom);
   rhotab(kp,k) = rhotab(k,kp);
   );
  );
*Select the next column to optimize, using mean square.
 maxrms = 0;
 loop(k ,
  rms = sqrt((sum(kp, sqr( rhotab(k,kp) ) ) )/(card(k) - 1));
  if (rms >= maxrms,
   maxrms = rms;
   currcol = ord(k);
  );
 );
*Iterate process until columns(n,k) is best design matrix.
counter = 0;
limitrho = 0.01;
DISPLAY "The current correlation root mean square is:", maxrms;
DISPLAY "Next column to check is,", currcol;
 while((maxrms > limitrho) and (counter < maxcount),
         SOLVE ONECOL USING MIP MINIMIZING Z;
*Update the columns(n,k) table with the optimized values of X.
         loop(k$(ord(k) = currcol),
                 loop(n,
                         loop(i,
                                 if(X.l(n,i) = 1,
                                 columns(n,k) = ord(i);
                                 );
                         );
                 );
         );
*Update correlation table based on optimized column.
        loop(k ,
                 loop(kp $(ord(kp) <> ord(k)),
                 numer = sum(n, sqr(columns(n,k) - columns(n,kp)));
                 rhotab(k,kp) = 1 - 6*(numer/denom);
                 rhotab(kp,k) = rhotab(k,kp);
                 );
         );
*Select the next column to optimize, using mean square.
         maxrms = 0;
         loop(k ,
         rms = sqrt((sum(kp, sqr(rhotab(k,kp))))/(card(k) - 1));
                 if (rms >= maxrms,
                 maxrms = rms;
                 currcol = ord(k);
                 );
         );
 counter = counter + 1;
 DISPLAY "The current correlation mean square is:", maxrms;
 DISPLAY "Next column to check is,", currcol;
 );
 DISPLAY "Best design matrix." , columns;
 DISPLAY "Final correlation matrix for design." , rhotab;