# O-Matrix Benchmark (5 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.
#
# To start the test: menu File -> Run Program... select "c:/<dir>/OMatrix.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"
  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(1200, 1200)/10)
      b = a'
      dim b(600, 2400)
      a = b';
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(1, 1) = timing
  print "Creation, transp., deformation of a 1200x1200 matrix (sec): ", timing
  
  # (2)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = abs(snormal(1250, 1250)/2);
    tic 
      b = a^1000; 
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(2, 1) = timing
  print "1250x1250 normal distributed random matrix ^1000____ (sec): ", timing

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

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

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

  times = sort(times)
  print "                      ----------------------------------------------"
  print "                      Trimmed mean (2 extremes eliminated): ", colmean(double(times(2::4,1)))
  print " "

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

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

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

  # (3)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(750, 750)
    tic
      b = det(a)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(3, 2) = timing
  print "Determinant of a 750x750 random matrix______________ (sec): ", timing
  
  # (4)
  cumulate = 0; b = 0
  for i = 1 to runs begin
    a = snormal(1000, 1000)
    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 1000x1000 matrix________ (sec): ", timing

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

  times = sort(times)
  print "                      ----------------------------------------------"
  print "                      Trimmed mean (2 extremes eliminated): ", colmean(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(225000, 1))
    tic	  
      b = (phi^a - (-phi)^(-a)) / sqrt(5.)
    timing = toc
    cumulate = cumulate + timing
  end
  timing = cumulate/runs
  times(1, 3) = timing
  print "225,000 Fibonacci numbers calculation (vector calc)_ (sec): ", timing

  # (2)
  cumulate = 0; a = 1500; 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 1500x1500 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, 35000))
    b = ceil(1000 * rand(1, 35000))
    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 35,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(22, 22))
    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
  times(5, 3) = timing
  print "Escoufier's method on a 22x22 matrix (mixed)________ (sec): ", timing
  times = sort(times)
  print "                      ----------------------------------------------"
  print "                      Trimmed mean (2 extremes eliminated): ", colmean(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): ", sum(colmean(double(times(2::4, 1::3))))/3
  print "                      --- End of test ---" 
end
clear