[Scipy-svn] r3070 - in trunk/Lib/sandbox/pyem: . profile_data

[scipy-svn at scipy.org]

Author: cdavid
Date: 2007-06-03 23:32:56 -0500 (Sun, 03 Jun 2007)
New Revision: 3070

Added:
   trunk/Lib/sandbox/pyem/data/
Modified:
   trunk/Lib/sandbox/pyem/densities2.py
   trunk/Lib/sandbox/pyem/profile_data/profile_densities.py
Log:
More benchmarking for basic operations in row vs col

Modified: trunk/Lib/sandbox/pyem/densities2.py
===================================================================
--- trunk/Lib/sandbox/pyem/densities2.py	2007-06-02 22:12:32 UTC (rev 3069)
+++ trunk/Lib/sandbox/pyem/densities2.py	2007-06-04 04:32:56 UTC (rev 3070)
@@ -1,7 +1,7 @@
 #! /usr/bin/python
 #
 # Copyrighted David Cournapeau
-# Last Change: Wed Dec 06 09:00 PM 2006 J
+# Last Change: Sat Jun 02 07:00 PM 2007 J
 
 # New version, with default numpy ordering.
 

Modified: trunk/Lib/sandbox/pyem/profile_data/profile_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/profile_data/profile_densities.py	2007-06-02 22:12:32 UTC (rev 3069)
+++ trunk/Lib/sandbox/pyem/profile_data/profile_densities.py	2007-06-04 04:32:56 UTC (rev 3070)
@@ -1,42 +1,78 @@
 import numpy as N
 from numpy.random import randn
-from scipy.sandbox.pyem import densities as D
-from scipy.sandbox.pyem import _c_densities as DC
-#import tables
 
+from numpy.ctypeslib import load_library, ndpointer
+from ctypes import cdll, c_uint, c_int, c_double, POINTER
+
+lib = load_library("blop.so", "file")
+
+arg1    = ndpointer(dtype=N.float64)
+arg2    = c_uint
+arg3    = c_uint
+arg4    = ndpointer(dtype=N.float64)
+arg5    = ndpointer(dtype=N.float64)
+
+lib.compute.argtypes    = [arg1, arg2, arg3, arg4, arg5]
+lib.compute.restype     = c_int
+# Compare computing per component likelihood for frame per row vs frame per column
+def component_likelihood(x, mu, va, log = False):
+    """expect one frame to be one row (rank 2). mu and var are rank 1 array."""
+    d = mu.size
+
+    return N.exp(N.sum((x - mu) ** 2, 1))
+
+def component_likelihood2(x, mu, va, log = False):
+    """expect one frame to be one column (rank 2). mu and var are rank 1 array."""
+    d = mu.size
+
+    y = (x[0] - mu[0]) ** 2
+    for i in range(1, d):
+        y += (x[i] - mu[i]) ** 2
+
+    return N.exp(y)
+
+def component_likelihood3(x, mu, va, log = False):
+    """expect one frame to be one row (rank 2). mu and var are rank 1 array."""
+    d = mu.size
+
+    y = N.empty(x.shape[0], x.dtype)
+    return lib.compute(x, x.shape[0], d, mu, y)
+
 def bench(func, mode = 'diag'):
-    #===========================================
-    # Diag Gaussian of dimension 20
-    #===========================================
     d       = 30
     n       = 1e5
     niter   = 10
 
     print "Compute %d times densities, %d dimension, %d frames" % (niter, d, n)
-    # Generate a model with k components, d dimensions
-    mu  = randn(1, d)
-    if mode == 'diag':
-        va  = abs(randn(1, d))
-    elif mode == 'full':
-        va  = randn(d, d)
-        va  = N.dot(va, va.transpose())
-
+    mu  = randn(d)
+    va  = abs(randn(d))
+    
     X   = randn(n, d)
     for i in range(niter):
         Y   = func(X, mu, va)
 
+def bench2(func, mode = 'diag'):
+    d       = 30
+    n       = 1e5
+    niter   = 10
+
+    print "Compute %d times densities, %d dimension, %d frames" % (niter, d, n)
+    mu  = randn(d)
+    va  = abs(randn(d))
+    
+    X   = randn(d, n)
+    for i in range(niter):
+        Y   = func(X, mu, va)
+
 def benchpy():
-    bench(D.gauss_den)
+    bench(component_likelihood)
 
-def benchc():
-    bench(DC.gauss_den)
+def benchpy3():
+    bench(component_likelihood3)
 
-def benchpyfull():
-    bench(D.gauss_den, 'full')
+def benchpy2():
+    bench2(component_likelihood2)
 
-def benchcfull():
-    bench(DC.gauss_den, 'full')
-
 if __name__ == "__main__":
     import hotshot, hotshot.stats
     profile_file    = 'gdenpy.prof'
@@ -48,7 +84,14 @@
 
     profile_file    = 'gdenc.prof'
     prof    = hotshot.Profile(profile_file, lineevents=1)
-    prof.runcall(benchc)
+    prof.runcall(benchpy2)
     p = hotshot.stats.load(profile_file)
     print p.sort_stats('cumulative').print_stats(20)
     prof.close()
+
+    profile_file    = 'gdenc.prof'
+    prof    = hotshot.Profile(profile_file, lineevents=1)
+    prof.runcall(benchpy3)
+    p = hotshot.stats.load(profile_file)
+    print p.sort_stats('cumulative').print_stats(20)
+    prof.close()