# O-Matrix benchmark test version 5.0
# Stefan Steinhaus stefan@steinhaus-net.de
#

clear

include .\stats\vcov.oms

function [r,v,e,w] = pca(z,k,tp) begin
	[m,C0,R0,Cn,Rn,G0,Q0] = vcov(z,k)
	if tp == 0 then C = C0
	else if tp == 1 then C = R0
	else if tp == 2 then C = G0
	else if tp == 3 then C = Q0
	else begin
		error("[pca] error: tp flag should be between 0 and 3")
		halt
	end
	e = novalue
	r = eigsym(C,e)
	r = reverse(r)
	r = abs(r)
	d = prod(r)
	v = cumsum(r)/colsum(r)
	e = reverse(e')'
	w = z*e
end

# 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)

begin
	runs=3
	result=0
	print
	print
	print "                  !!! O-Matrix - Benchmark program !!!"
	print "                  ===================================="
	print
	print
	print "Function / timing                                      result (sec.)"
	print "-----------------------------------------------------------------------"

	# Misc. operation

	for counter=1 to runs begin
	  t0 = time
	  file="D:\\Mathematik\\Ncrunch5\\Currency2.txt"
	  header=read(file, "char",1)
	  datmat=read(file, "double",3159,38)
	  close(file)
	  Limits_Area={1,261,522,784,1045,1305,1565,1827,2088,2349,2610,2870,3131}
	  Diff_Area={260,261,262,261,260,260,262,261,261,261,260,260,27}
	  for i=1 to 2000 begin
	    for j=1 to 13 begin
		Year_Data=datmat.blk(Limits_Area(j),5,Diff_Area(j),34)
		for k=1 to 34 begin
		  Min_Year=mins(Year_Data.col(k))
		  Max_Year=max(Year_Data.col(k))
		  Mean_Year=colmean(Year_Data.col(k))
		  GainLoss_Year=100-(Year_Data(1,k)+0.00000001)/(Year_Data(Diff_Area(j),k)+0.00000001)*100
		end
	    end
	  end
	  t1 = time - t0
	  result=result+t1
	end
	result=result/runs
	print "IO test & descriptive statistics___________________: ",result

	result=0
	N = 15000
	for counter=1 to runs begin
	  a=1
	  t0 = time
	  for i = 1 to N begin
	    for j = 1 to N begin
		a=a+i+j
	    end
	  end
	  t1 = time - t0
	  result=result+t1
	end
	result=result/runs
	print "Loop test__________________________________________: ",result

	N = 2000
	result=0
	for counter=1 to runs begin
        for counter2=1 to 200 begin
          x=unirnd(double(1.0),double(1.01),N,N)
	    t0 = time
	    x = x^1000.
	    result = result + time-t0
        end
	end
	result=result/runs
	print "2000x2000 random matrix^1000_______________________: ",result

	N = 1000000
	result=0
	for counter=1 to runs begin
        for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,1)
	    t0 = time
	    X = sort(y)
	    result = result + time-t0
        end
	end
	result=result/runs
	print "1000000 values sorted ascending____________________: ", result

	# Analysis
 
	N = 1048576
	result=0
	for counter=1 to runs begin
        for counter2=1 to 100 begin
          y=complex(unirnd(double(0.0),double(1.0),N,1))
	    t0 = time
	    X = fft(y)
	    result = result + time-t0
        end
	end
	result=result/runs
	print "FFT over 1048576 (2^20) values_____________________: ", result

	# Algebra

	N = 1500
	result=0
	for counter=1 to runs begin
        for counter2=1 to 100 begin
	    x=unirnd(double(0.0),double(1.0),N,N)
	    t0 = time
	    y=det(x)
	    result = result + time-t0
        end
	end
	result=result/runs
	print "Determinant of a 1500x1500 random matrix___________: ",result

	N = 1500
	result=0
	for counter=1 to runs begin
        for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,N)
	    t0 = time
	    X = inv(y)	
	    result = result + time-t0
        end
	end
	result=result/runs
	print "Inverse of 1500x1500 random matrix_________________: ", result

	N = 1200
	result=0
	for counter=1 to runs begin
	  y=unirnd(double(0.0),double(1.0),N,N)
	  t0 = time
	  X = eigen(y)
	  t1 = time - t0
	  result=result+t1
	end
	result=result/runs
	print "Eigenvalues of a 1200x1200 random matrix___________: ",result

	N = 1500
	result=0
	for counter=1 to runs begin
        for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,N)
	    y = y'*y
          t0 = time
	    X = cholesky(y)
	    result = result + time-t0
	  end
	end
	result=result/runs
	print "Cholesky decompositon of 1500x1500 matrix__________: ",result

	N = 1500
	result=0
	for counter=1 to runs begin
	  for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,N)
          t0 = time
	    X = y'*y
	    result = result + time-t0
	  end
	end
	result=result/runs
	print "1500x1500 cross-product matrix_____________________: ",result

	# Number theory

	N = 10000000
	result=0
	phi = 1.6180339887498949
	for counter=1 to runs begin
	  a=unirnd(double(100.0),double(1000.0),N,1)
	  t0=time
	  b=(phi^a-(-phi)^(-a))/sqrt(5.)
	  result = result + time-t0
	end
	result=result/runs
	print "Calculation of 10000000 fibonacci numbers__________: ",result

	# Stochastic-statistic

	N = 10000
	result=0
	for counter=1 to runs begin
		a=unirnd(double(0.0),double(1.0),N,N/10)
		t0=time
		[r,v,e,w] = pca(a,0,0)
		result = result + time-t0
	end
	result=result/runs
	print "Principal component analysis over 10000x1000 matrix: ",result

	# O-Matrix implements ln(gamma(x)) instead of gamma(x),
	# since the latter will overflow many computers' floating-point
	# representation at modest values of x.   
	N = 1500
	result=0
	for counter=1 to runs begin
	  for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,N)
          t0=time
	    X = exp(gammln(y))
	    result = result + time-t0
	  end
	end
	result=result/runs
	print "Gamma function over a 1500x1500 matrix_____________: ",result

	N = 1500
	result=0
	for counter=1 to runs begin
	  for counter2=1 to 100 begin
	    y=unirnd(double(0.0),double(1.0),N,N)
          t0=time
	    X = erf(y)
	    result = result + time-t0
	  end
	end
	result=result/runs
	print "Gaussian error function over a 1500x1500 matrix____: ",result
	
	N = 1000
	result=0
	for counter=1 to runs begin
	  for counter2=1 to 100 begin
	    a = unirnd(double(0.0),double(1.0),N,N)
	    b = 1::N
          t0=time
	    X = a\b;
	    result = result + time-t0
	  end
	end
	result=result/runs
	print "Linear regression over 1000x1000 matrix____________: ",result
end