## Rlab Benchmark (2 April 2001)
## Author : Philippe Grosjean
## eMail  : phgrosjean@sciviews.org
## Web    : http://www.sciviews.org
## License: GPL 2 or above at your convenience (see: http://www.gnu.org)
##
## Several tests are adapted from:
##***************************************************************************
##* Matlab Benchmark program version 2.0                                    *
##* Author : Stefan Steinhaus                                               *
##* EMAIL  : stst@informatik.uni-frankfurt.de                               * 
##* This program is public domain. Feel free to copy it freely.             *
##***************************************************************************
## Escoufier's equivalents vectors (III.5) is adapted from Planque & Fromentin, 1996 
## Ref: Escoufier Y., 1970. Echantillonnage dans une population de variables
## aleatoires réelles. Publ. Inst. Statis. Univ. Paris 19 Fasc 4, 1-47.
##
## Type: load("c:/<dir>/Rlab.r"); to start the test

format(32);
srand();
runs = 3;			## Number of times the tests are executed
times = zeros(5, 3);

printf("   Rlab Benchmark\n");
printf("   ==============\n");
printf("Number of times each test is run__________________________: %f\n", runs, "\n");
printf(" \n");

printf("   I. Matrix calculation\n");
printf("   ---------------------\n");

## (1)
cumulate = 0; a = 0; b = 0;
for (i in 1:runs) {
  tic();
    rand("normal", 0, 0.1);
    a = abs(rand(1200, 1200));
    b = a';
    a = reshape(b, 600, 2400);
    b = a';
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[1; 1] = timethis;
printf("Creation, transp., deformation of a 1200x1200 matrix (sec): %f \n", timethis, "\n");
clear(a, b, timethis);

## (2)
cumulate = 0; b = 0;
rand("uniform", 0, 1);
for (i in 1:runs) {
  a = abs(rand(1250, 1250));
  tic(); 
    b = a.^1000; 
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[2; 1] = timethis;
printf("1250x1250 normal distributed random matrix ^1000____ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (3)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(1100000, 1);
  tic();
    b = sort(a).val;
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[3; 1] = timethis;
printf("Sorting of 1,100,000 random values__________________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (4)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(550, 550);
  tic();
    b = a'*a;
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[4; 1] = timethis;
printf("550x550 cross-product matrix (b = a' * a)___________ (sec): %f\n", timethis, "\n");
clear (a, b, timethis);

## (5)
cumulate = 0; c = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(700, 700);
  b = 1:700;
  tic();
    c = a\b';
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[5; 1] = timethis;
printf("Linear regression over a 700x700 matrix (c = a \ b') (sec): %f\n", timethis, "\n");
clear(a, b, c, timethis);

times = sort(times).val;
printf("                      ----------------------------------------------\n");
printf("                      Trimmed mean (2 extremes eliminated): %f\n", sum(times[2:4;1])/3);
printf(" \n");

printf("   II. Matrix functions\n");
printf("   --------------------\n");

## (1)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(900000, 1);
  tic();
    b = fft(a);
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[1; 2] = timethis;
printf("FFT over 900,000 random values______________________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (2)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(220, 220);
  tic();
    b = eig(a).val;
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[2; 2] = timethis;
printf("Eigenvalues of a 220x220 random matrix______________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (3)
cumulate = 0; b = 0;
rand("uniform", 0, 0.1);
for (i in 1:runs) {
  a = rand(750, 750);
  tic();
    b = det(a);
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[3; 2] = timethis;
printf("Determinant of a 750x750 random matrix______________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (4)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(1000, 1000);
  a = a'*a;
  tic();
    b = chol(a);
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[4; 2] = timethis;
printf("Cholesky decomposition of a 1000x1000 matrix________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (5)
cumulate = 0; b = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  a = rand(500, 500);
  tic();
    b = inv(a);
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[5; 2] = timethis;
printf("Inverse of a 500x500 random matrix__________________ (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

times = sort(times).val;
printf("                      ----------------------------------------------\n");
printf("                      Trimmed mean (2 extremes eliminated): %f\n", sum(times[2:4;2])/3);
printf(" \n");

printf("   III. Programmation\n");
printf("   ------------------\n");

## (1)
cumulate = 0; a = 0; b = 0; phi = 1.6180339887498949;
rand("uniform", 0, 1000);
for (i in 1:runs) {
  a = floor(rand(225000, 1));
  tic();	  
    b = (phi.^a - (-phi).^(-a)) / sqrt(5);
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[1; 3] = timethis;
printf("225,000 Fibonacci numbers calculation (vector calc)_ (sec): %f\n", timethis, "\n");
clear(a, b, phi, timethis);

## (2)
cumulate = 0; a = 1500; b = 0;
for (i in 1:runs) {
  tic();
    b = ones(a, a)./((1:a)' * ones(1, a) + ones(a, 1) * (0:(a-1)));    
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[2; 3] = timethis;
printf("Creation of a 1500x1500 Hilbert matrix (matrix calc) (sec): %f\n", timethis, "\n");
clear(a, b, timethis);

## (3)
gcd2 = function(a, b)
{
  if (any(b > 1.0E-4)) {
    for (i in 1:b.n) {
	  if (b[i] == 0) {b[i] = a[i];}
    }
    return $self(b, mod(a, b));
  else
    return a;
  }
};

cumulate = 0; c = 0;
rand("uniform", 0, 1000);
for (i in 1:runs) {
  a = ceil(rand(35000, 1));
  b = ceil(rand(35000, 1));
  tic();	  
    c = gcd2(a, b);                            ## gcd2 is a recursive function
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[3; 3] = timethis;
printf("Grand common divisors of 35,000 pairs (recursion)___ (sec): %f\n", timethis, "\n");
clear (a, b, c, timethis);

## (4)
cumulate = 0; b = 0;
for (i in 1:runs) {
  b = zeros(220, 220);
  tic();
    for (j in 1:220) {
      for (k in 1:220) {
        b[k; j] = abs(j - k) + 1;
      }
    }
  timethis = toc();
  cumulate = cumulate + timethis;
}
timethis = cumulate/runs;
times[4; 3] = timethis;
printf("Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): %f\n", timethis, "\n");
clear(b, j, k, timethis);

## (5)
cumulate = 0; p = 0; vt = 0; vr = 0; vrt = 0; rvt = 0; RV = 0; j = 0; k = 0;
x2 = 0; R = 0; Rxx = 0; Ryy = 0; Rxy = 0; Ryx = 0; Rvmax = 0; f = 0;
rand("normal", 0, 1);
for (i in 1:runs) {
  x = abs(rand(22, 22));
  tic();
    ## Calculation of Escoufier's equivalent vectors
    p = x.nc;
    vt = [1:p];                                ## Variables to test
    vr = [];                                   ## Result: ordered variables
    RV = [1:p];                                ## Result: correlations
    for (j in 1:p) {                           ## loop on the variable number
      Rvmax = 0;
      for (k in 1:(p-j+1)) {                   ## loop on the variables
        x2 = [x, x[ ; vr], x[ ; vt[k]]];       ## New table to test
        ## R = corrcoef(x2);                   ## Correlations table
        ## Not in Rlab, so the 5 following lines do it
        q = x2.nr;
        x2 = x2 - ones(q, 1) * sum(x2) / q;
		c = conj(x2' * x2 / (q - 1));
		d=diag(c);
		R = c ./ sqrt(d * d');
        Ryy = R[1:p; 1:p];
        Rxx = R[p+1:p+j; p+1:p+j];
        Rxy = R[p+1:p+j; 1:p];
        Ryx = Rxy';
        rvt = trace(Ryx*Rxy)/((trace(Ryy*Ryy)*trace(Rxx*Rxx))^0.5); ## RV calculation
        if (rvt > Rvmax) {
          Rvmax = rvt;                         ## test of RV
          vrt = vt[k];                         ## temporary held variable
        }
      }
      vr[j] = vrt;                             ## Result: variable
      RV[j] = Rvmax;                           ## Result: correlation
      f = find(vt!=vr[j]);                     ## identify the held variable
      vt = vt[f];                              ## reidentify variables to test
    }
  timethis = toc();
  cumulate = cumulate + timethis;
}
times[5; 3] = timethis;
printf("Escoufier's method on a 22x22 matrix (mixed)________ (sec): %f\n", timethis, "\n");
clear(x, p, vt, vr, vrt, rvt, RV, j, k, timethis);
clear(x2, R, Rxx, Ryy, Rxy, Ryx, Rvmax, f); 

times = sort(times).val;
printf("                      ----------------------------------------------\n");
printf("                      Trimmed mean (2 extremes eliminated): %f\n", sum(times[2:4;3])/3);
printf(" \n");

printf(" \n");
printf("Total time for all 15 tests_________________________ (sec): %f\n", sum(sum(times)), "\n");
printf("Overall mean (sum of I, II and III trimmed means/3)_ (sec): %f\n", sum(sum(times[2:4; 1:3]))/9);
clear(cumulate, times, runs, i);
printf("                      --- End of test ---\n");