// Ox Benchmark (3 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.             *
//***************************************************************************
// and from the corresponding Ox Benchmark V.2 from S. Steinhaus
// 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.
//
// On the dos prompt, type:
//     cd c:\<dir>
//     oxl Ox
// to start the test

#include <oxstd.h>

// This is for test III.3
gcd2(a, b)
{
  if (b <= 1.0E-4)
    return a;
  else {
    b[vecindex(b, 0)] = a[vecindex(b, 0)];
    return gcd2(b, fmod(a, b));
  }
}

main()
{
  decl runs = 3;			// Number of times the tests are executed
  decl times, cumul, timeini;
  decl a, b, c, phi, i, j, k;
  decl x, x2, p, vt, vr, vrt, rvt;
  decl RV, R, Rxx, Ryy, Rxy, Ryx, Rvmax;
	
  times = zeros(5, 3);

  print("   Ox Benchmark\n");
  print("   ================\n");
  print("Number of times each test is run__________________________: ", runs, "\n\n");

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


  // (1)
  cumul = 0; a = 0; b = 0;
  for (i = 0; i < runs; i++) {
    timeini = timer();
      a = fabs(rann(1200, 1200)/10);
      b = a';
      a = reshape(b, 600, 2400);
      b = a';
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[0][0] = cumul;
  print("Creation, transp., deformation of a 1200x1200 matrix (sec): " , cumul, "\n");
  delete a; delete b;

  // (2)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = fabs(rann(1250, 1250)/2);
    timeini = timer(); 
      b = (a).^1000; 
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[1][0] = cumul;
  print("1250x1250 normal distributed random matrix ^1000____ (sec): ", cumul, "\n");
  delete a; delete b;

  // (3)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(1100000, 1);
    timeini = timer();
      b = sortc(a);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[2][0] = cumul;
  print("Sorting of 1,100,000 random values__________________ (sec): ", cumul, "\n");
  delete a; delete b; 

  //(4)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(550, 550);
    timeini = timer();
      b = a'a;
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[3][0] = cumul;
  print("550x550 cross-product matrix (b = a' * a)___________ (sec): ", cumul, "\n");
  delete a; delete b;

  // (5)
  cumul = 0; c = 0;
  for (i = 0; i < runs; i++) {
    a = rann(700, 700);
    b = range(1, 700);
    timeini = timer();
      olsr(b, a, &c);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[4][0] = cumul;
  print("Linear regression over a 700x700 matrix (c = a \ b') (sec): ", cumul, "\n");
  delete a; delete b; delete c;

  times = sortc(times);
  print("                      --------------------------------------------\n");
  print("                      Trimmed mean (2 extremes eliminated): ", meanc(times[1:3][0]), "\n\n");
  
  print("   II. Matrix functions\n");
  print("   --------------------\n");


  // (1)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(1, 900000);
    timeini = timer();
      b = fft(a);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[0][1] = cumul;
  print("FFT over 900,000 random values______________________ (sec): ", cumul, "\n");
  delete a; delete b;

  // (2)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(220, 220);
    timeini = timer();
      eigen(a, &b);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[1][1] = cumul;
  print("Eigenvalues of a 220x220 random matrix______________ (sec): ", cumul, "\n");
  delete a; delete b;

  // (3)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(750, 750)/10;
    timeini = timer();
      b = determinant(a);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[2][1] = cumul;
  print("Determinant of a 750x750 random matrix______________ (sec): ", cumul, "\n");
  delete a; delete b;

  // (4)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(1000, 1000);
    a = a'a;
    timeini = timer();
      b = choleski(a);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[3][1] = cumul;
  print("Cholesky decomposition of a 1000x1000 matrix________ (sec): ", cumul, "\n");
  delete a; delete b;

  // (5)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    a = rann(500, 500);
    timeini = timer();
      b = invert(a);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[4][1] = cumul;
  print("Inverse of a 500x500 random matrix__________________ (sec): ", cumul, "\n");
  delete a; delete b; 

  times = sortc(times);
  print("                      --------------------------------------------\n");
  print("                      Trimmed mean (2 extremes eliminated): ", meanc(times[1:3][1]), "\n\n");
  
  print("   III. Programmation\n");
  print("   ------------------\n");

  // (1)
  cumul = 0; a = 0; b = 0; phi = 1.6180339887498949;
  for (i = 0; i < runs; i++) {
    a = floor(1000 * ranu(225000, 1));
    timeini = timer();	  
      b = (phi^a - (-phi)^(-a)) / sqrt(5);
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[0][2] = cumul;
  print("225,000 Fibonacci numbers calculation (vector calc)_ (sec): ", cumul, "\n");
  delete a; delete b; delete phi;

  // (2) 
  cumul = 0; a = 1500; b = 0;
  for (i = 0; i < runs; i++) {
    timeini = timer();
      b = ones(a, a)./((range(1, a))' * ones(1, a) + ones(a, 1) * (range(0, a-1)));    
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[1][2] = cumul;
  print("Creation of a 1500x1500 Hilbert matrix (matrix calc) (sec): ", cumul, "\n");
  delete a; delete b;

  // (3)
  cumul = 0; c = 0;
  for (i = 0; i < runs; i++) {
    a = ceil(1000 * ranu(35000, 1));
    b = ceil(1000 * ranu(35000, 1));
    timeini = timer();	  
      c = gcd2(a, b);                            // gcd2 is a recursive function
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[2][2] = cumul;
  print("Grand common divisors of 35,000 pairs (recursion)___ (sec): ", cumul, "\n");
  delete a; delete b; delete c;

  // (4)
  cumul = 0; b = 0;
  for (i = 0; i < runs; i++) {
    b = zeros(220, 220);
    timeini = timer();
      for (j = 0; j < 220; j++) {
        for (k = 0; k < 220; k++) {
          b[k][j] = fabs(j - k) + 1;
        }
      }
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[3][2] = cumul;
  print("Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ", cumul, "\n");
  delete b; delete j; delete k;

  // (5)
  cumul = 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;
  for (i = 0; i < runs; i++) {
    x = fabs(rann(22, 22));
    timeini = timer();
      // Calculation of Escoufier's equivalent vectors
      p = sizec(x) - 1;
      vt = range(0, p);                          // Variables to test
      vr = ones(p+1, 1);                         // Result: ordered variables
      RV = ones(p+1, 1);                         // Result: correlations
      for (j = 0; j < p+1; j++) {                // loop on the variable number
        Rvmax = 0;
        for (k = 0; k < p-j+1; k++) {            // loop on the variables
          x2 = x~x[][vt[k]];                    // New table to test
          R = correlation(x2);                   // Correlations table
          Ryy = R[0:p][0:p];
          Rxx = R[p+1:p+j+1][p+1:p+j+1];
          Rxy = R[p+1:p+j+1][0:p];
          Ryx = Rxy';
          rvt = trace(Ryx*Rxy)/((trace(Ryy^2)*trace(Rxx^2))^0.5); // RV calculation
          if (rvt > Rvmax) {
            Rvmax = rvt;                         // test of RV
            vrt = vt[k];                         // temporary held variable
          }
        }
        vr[j][0] = vrt + 1;                      // +1 because we index variables from 1!
        RV[j][0] = Rvmax;                        // Result: correlation
        x = x~x[][vrt];                          // add the held variable to x
        vt = deletec(vt, vrt);                   // reidentify variables to test
      }
    cumul += (timer() - timeini) / 100;
  }
  cumul /= runs;
  times[4][2] = cumul;
  print("Escoufier''s method on a 22x22 matrix (mixed)________ (sec): ", cumul, "\n");
  delete x; delete p; delete vt; delete vr; delete vrt; delete rvt; delete RV; delete j; delete k;
  delete x2; delete R; delete Rxx; delete Ryy; delete Rxy; delete Ryx; delete Rvmax; 

  times = sortc(times);
  print("                      --------------------------------------------\n");
  print("                      Trimmed mean (2 extremes eliminated): ", meanc(times[1:3][2]), "\n\n");
  
  print("\n");
  print("Total time for all 15 tests_________________________ (sec): ", sumr(sumc(times)), "\n");
  print("Overall mean (sum of I, II and III trimmed means/3)_ (sec): ", meanr(meanc(times[1:3][0:2])), "\n");
  delete cumul; delete timeini; delete times; delete runs; delete i;
  print("                      --- End of test ---\n");  
}