# O-Matrix Benchmark 2 (8 March 2003)
# version 2, scaled to get 1 +/- 0.1 sec with R 1.6.2
# using the standard ATLAS library (Rblas.dll)
# on a Pentium IV 1.6 Ghz with 1 Gb Ram on Win XP pro
# 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.
#
# To start the test: menu File -> Run Program... select "c:/<dir>/OMatrix2.oms"

clear
# By default O-Matrix multi-threads the various windows within the
# application.  This allows editing, graphics etc. to continue
# while a calculation is running.  This reduces performance
# about 10-15% so we are turning it off here
interrupt(0)
clc

begin
  runs = 3			# Number of times the tests are executed
  times = zeros(5, 3)

  print "   O-Matrix Benchmark 2"
  print "   ===================="
  print "Number of times each test is run__________________________: ", runs
  print " "


  print "   I. Matrix calculation"
  print "   ---------------------"


  # (1)
  cumulate = 0; a = 0; b = 0
  for i = 1 to runs begin
    tic
      a = abs(snormal(1500, 1500)/10)
      b = a'
      dim b(750, 3000)
      a = b';
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(1, 1) = timing
  print "Creation, transp., deformation of a 1500x1500 matrix (sec): ", timing
  
  # (2)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = abs(snormal(800, 800)/2);
    tic 
      b = a^1000; 
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(2, 1) = timing
  print "800x800 normal distributed random matrix ^1000______ (sec): ", timing

  # (3)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(2000000, 1)
    tic
      b = sort(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(3, 1) = timing
  print "Sorting of 2,000,000 random values__________________ (sec): ", timing

  # (4)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(700, 700)
    tic
      b = a'*a
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(4, 1) = timing
  print "700x700 cross-product matrix (b = a' * a)___________ (sec): ", timing

  # (5)
  cumulate = 0; c = 0
  for i = 1 to runs begin
    a = snormal(600, 600)
    b = 1::600
    tic
      c = a\b
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(5, 1) = timing
  print "Linear regression over a 600x600 matrix (c = a \ b') (sec): ", timing

  times(1::5,1) = sort(times(1::5,1))
  print "                ----------------------------------------------------"
  print "                Trimmed geom. mean (2 extremes eliminated): ", exp(colmean(log(double(times(2::4,1)))))
  print " "

  print "   II. Matrix functions"
  print "   --------------------"

  # (1)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(800000, 1)
    tic
      b = fft(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(1, 2) = timing
  print "FFT over 800,000 random values______________________ (sec): ", timing

  # (2)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(320, 320)
    tic
      b = eigen(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(2, 2) = timing
  print "Eigenvalues of a 320x320 random matrix______________ (sec): ", timing

  # (3)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(650, 650)
    tic
      b = det(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(3, 2) = timing
  print "Determinant of a 650x650 random matrix______________ (sec): ", timing
  
  # (4)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(900, 900)
    a = a'*a
    tic
      b = cholesky(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(4, 2) = timing
  print "Cholesky decomposition of a 900x900 matrix__________ (sec): ", timing

  # (5)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(400, 400)
    tic
      b = inv(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(5, 2) = timing
  print "Inverse of a 400x400 random matrix__________________ (sec): ", timing

  times(1::5,2) = sort(times(1::5, 2))
  print "                ----------------------------------------------------"
  print "                Trimmed geom. mean (2 extremes eliminated): ", exp(colmean(log(double(times(2::4,2)))))
  print " "

  print "   III. Programmation"
  print "   ------------------"

  # (1)
  cumulate = 0; a = 0; b = 0; phi = 1.6180339887498949
  for i = 1 to runs begin
    a = floor(1000 * rand(750000, 1))
    tic	  
      b = (phi^a - (-phi)^(-a)) / sqrt(5.)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(1, 3) = timing
  print "750,000 Fibonacci numbers calculation (vector calc)_ (sec): ", timing

  # (2)
  cumulate = 0; a = 2250; b = 0
  for i = 1 to runs begin
    tic
      b = ones(a, a)/double((1::a) * ones(1, a) + ones(a, 1) * (0::(a-1))')    
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(2, 3) = timing
  print "Creation of a 2250x2250 Hilbert matrix (matrix calc) (sec): ", timing

  # (3)
  cumulate = 0; c = 0
  local function gcd2(a, b) begin
    if sum(b > 1.0E-4) == 0 then return a
    b([find(b == 0)]) = a([find(b == 0)])
    return gcd2(b, rem(a, b))
  end
  for i = 1 to runs begin
    a = ceil(1000 * rand(1, 70000))
    b = ceil(1000 * rand(1, 70000))
    tic	  
      c = gcd2(a, b)                            # gcd2 is a recursive function
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(3, 3) = timing
  print "Grand common divisors of 70,000 pairs (recursion)___ (sec): ", timing
      
  # (4)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    b = zeros(220, 220)
    tic
      for j = 1 to 220 begin
        for k = 1 to 220 begin
          b(k,j) = abs(j - k) + 1
        end
      end
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(4, 3) = timing
  print "Creation of a 220x220 Toeplitz matrix (loops)_______ (sec): ", timing  

  # (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
  for i = 1 to runs begin
    x = abs(snormal(37, 37))
    tic
      # Calculation of Escoufier's equivalent vectors
      p = size(x, 2)
      vt = [1::p]                              # Variables to test
      vr = zeros(p, 1)                         # Result: ordered variables
      RV = [1::p]                              # Result: correlations
      for j = 1 to p begin                     # loop on the variable number
        Rvmax = 0
        for k = 1 to (p-j+1) begin             # loop on the variables
          if (j == 1) then x2 = [x, x(:, vt(k))]
          else x2 = [x, x(:, nonzeros(vr)), x(:, vt(k))] # New table to test
          R = colcor(x2)                       # Correlations table
          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 then begin
            Rvmax = rvt                        # test of RV
            vrt = vt(k)                        # temporary held variable
          end
        end
        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
      end
    timing = toc
    cumulate = cumulate + timing;
  end
  timing = cumulate/runs
  times(5, 3) = timing
  print "Escoufier's method on a 37x37 matrix (mixed)________ (sec): ", timing
  times(1::5,3) = sort(times(1::5,3))
  print "                ----------------------------------------------------"
  print "                Trimmed geom. mean (2 extremes eliminated): ", exp(colmean(log(double(times(2::4,3)))))
  print " "

  print " "
  print "Total time for all 15 tests_________________________ (sec): ", sum(sum(double(times)))
  print "Overall mean (sum of I, II and III trimmed means/3)_ (sec): ", exp(sum(colmean(log(double(times(2::4, 1::3)))))/3)
  print "                      --- End of test ---" 
end
clear