'''
Benchmark for numpy, numarray and Numeric blas operations.
'''

import time
from numpy import *
import numarray, numarray.random_array
import Numeric, RandomArray

order = ['x.T*y', 'x*y.T', 'A*x', 'A*B', 'A.T*x', 'half', '2in2']

def benchOld(N, ra, d, rep):
    A = ra.random((d, d))
    B = ra.random((d, d))
    x = ra.random(d)
    y = ra.random(d)

    timings = []

    # Inner product
    t = time.clock()
    for i in range(d*rep): res = N.dot(x, y)
    timings.append(time.clock() - t)

    # Outer product
    y = N.reshape(y, (1, d)).copy()
    x = N.reshape(x, (d, 1)).copy()
    t = time.clock()
    for i in range(rep): res = N.dot(x, y)
    timings.append(time.clock() - t)
    
    # Matrix times vector
    x = N.reshape(x, (d,)).copy()
    t = time.clock()
    for i in range(rep): res = N.dot(A, x)
    timings.append(time.clock() - t)

    # Matrix times matrix
    t = time.clock()
    for i in range(rep): res = N.dot(A, B)
    timings.append(time.clock() - t)
    
    # With transpose 
    t = time.clock()
    for i in range(rep): res = N.dot(N.transpose(A), x)
    timings.append(time.clock() - t)

    # Select
    start = d/2
    t = time.clock()
    for i in range(rep): res = N.dot(A[start:,start:], x[start:])
    timings.append(time.clock() - t)

    # Select 2 in 2
    ind = N.arange(0, d, 2)
    t = time.clock()
    for i in range(rep): res = N.dot(A[::2,::2], x[::2])
    timings.append(time.clock() - t)

    return timings
    
def bench(d, rep):
    A = rand(d, d)
    B = rand(d, d)
    x = rand(d)
    y = rand(d)

    MA = mat(A)
    MB = mat(B)
    mx = mat(x).T.copy()
    my = mat(y).T.copy()

    arrayTime = []
    matrixTime = []

    # Inner product
    t = time.clock()
    for i in range(d*rep): res = dot(x, y)
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(d*rep): res = mx.T*my
    matrixTime.append(time.clock() - t)

    # Outer product
    y = y.reshape(1, d).copy()
    x = x.reshape(d, 1).copy()
    t = time.clock()
    for i in range(rep): res = dot(x, y)
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = mx*my.T
    matrixTime.append(time.clock() - t)
    
    # Matrix times vector
    t = time.clock()
    for i in range(rep): res = dot(A, x)
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = MA*mx
    matrixTime.append(time.clock() - t)

    # Matrix times matrix
    t = time.clock()
    for i in range(rep): res = dot(A, B)
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = MA*MB
    matrixTime.append(time.clock() - t)
    
    # With transpose 
    t = time.clock()
    for i in range(rep):
        M = transpose(A).copy()
        res = dot(M, x)
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = MA.T*mx
    matrixTime.append(time.clock() - t)

    # Select
    start = d/2
    t = time.clock()
    for i in range(rep): res = dot(A[start:,start:], x[start:])
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = MA[start:,start:]*mx[start:]
    matrixTime.append(time.clock() - t)

    # Select 2 in 2
    t = time.clock()
    for i in range(rep): res = dot(A[::2,::2], x[::2])
    arrayTime.append(time.clock() - t)
    t = time.clock()
    for i in range(rep): res = MA[::2,::2]*mx[::2]
    matrixTime.append(time.clock() - t)

    return arrayTime, matrixTime

if __name__ == '__main__':
    dims = [5, 50, 500]
    reps = [100000, 100000, 1000]
    print 'Tests  %+7s %+7s %+7s %+7s %+7s %+7s %+7s' % tuple(order)
    for i in range(len(dims)):
        print 'Dimension:', dims[i]
        arrayTime, matrixTime = bench(dims[i], reps[i])
        print 'Array  %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f' % tuple(arrayTime)
        print 'Matrix %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f' % tuple(matrixTime)
        numarrayTime = benchOld(numarray, numarray.random_array,
                                dims[i], reps[i])
        print 'NumArr %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f' % tuple(numarrayTime)
        numericTime = benchOld(Numeric, RandomArray,
                               dims[i], reps[i])
        print 'Numeri %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f' % tuple(numericTime)