[Scipy-svn] r2281 - in trunk/Lib/sandbox/pyem: . pyem pyem/profile_data pyem/src pyem/tests

scipy-svn at scipy.org scipy-svn at scipy.org
Thu Oct 12 08:48:25 CDT 2006


Author: cdavid
Date: 2006-10-12 08:48:08 -0500 (Thu, 12 Oct 2006)
New Revision: 2281

Added:
   trunk/Lib/sandbox/pyem/count.sh
   trunk/Lib/sandbox/pyem/pyem/profile_data/
   trunk/Lib/sandbox/pyem/pyem/profile_data/profile_densities.py
   trunk/Lib/sandbox/pyem/pyem/profile_data/profile_gmm.py
   trunk/Lib/sandbox/pyem/pyem/tests/test_kmean.py
Removed:
   trunk/Lib/sandbox/pyem/pyem/profile_densities.py
   trunk/Lib/sandbox/pyem/pyem/profile_gmm.py
Modified:
   trunk/Lib/sandbox/pyem/.bzrignore
   trunk/Lib/sandbox/pyem/Changelog
   trunk/Lib/sandbox/pyem/MANIFEST.in
   trunk/Lib/sandbox/pyem/pyem/__init__.py
   trunk/Lib/sandbox/pyem/pyem/_c_densities.py
   trunk/Lib/sandbox/pyem/pyem/densities.py
   trunk/Lib/sandbox/pyem/pyem/kmean.py
   trunk/Lib/sandbox/pyem/pyem/src/Makefile
   trunk/Lib/sandbox/pyem/pyem/src/c_gden.c
   trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py
   trunk/Lib/sandbox/pyem/setup.py
Log:
[pyem @ david at ar.media.kyoto-u.ac.jp-20061003092950-0daedfcee3ce611e]
MSG
David Cournapeau <david at ar.media.kyoto-u.ac.jp> | 2006-10-03 18:29:50 +0900 (Tue, 03 Oct 2006)

Modified: trunk/Lib/sandbox/pyem/.bzrignore
===================================================================
--- trunk/Lib/sandbox/pyem/.bzrignore	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/.bzrignore	2006-10-12 13:48:08 UTC (rev 2281)
@@ -23,3 +23,6 @@
 ../MSG
 MSG
 exinfo.py
+blop
+*.prog
+*.prof

Modified: trunk/Lib/sandbox/pyem/Changelog
===================================================================
--- trunk/Lib/sandbox/pyem/Changelog	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/Changelog	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,3 +1,9 @@
+pyem (0.5.2) Tue, 03 Oct 2006 18:28:13 +0900
+
+	* Update tests so that they work within the numpy test framework
+
+-- David Cournapeau <david at ar.media.kyoto-u.ac.jp> 
+
 pyem (0.5.2) Tue, 29 Aug 2006 14:53:49 +0900
 
 	* Add a class method to GM to create a model directly from

Modified: trunk/Lib/sandbox/pyem/MANIFEST.in
===================================================================
--- trunk/Lib/sandbox/pyem/MANIFEST.in	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/MANIFEST.in	2006-10-12 13:48:08 UTC (rev 2281)
@@ -3,4 +3,4 @@
 include Changelog
 include TODO
 include LICENSE.txt
-exclude pyem/_c_densities.py
+#exclude pyem/_c_densities.py

Added: trunk/Lib/sandbox/pyem/count.sh
===================================================================
--- trunk/Lib/sandbox/pyem/count.sh	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/count.sh	2006-10-12 13:48:08 UTC (rev 2281)
@@ -0,0 +1,23 @@
+#! /bin/sh
+# Last Change: Mon Sep 11 08:00 PM 2006 J
+
+n=0
+np=0
+nc=0
+
+files=`find . -name '*.py'`
+
+for i in $files; do
+    tp=`wc "$i" | tr -s " " | cut -f 2 -d" "`
+    let np="$np + $tp"
+done
+
+files=`find . -name '*.[ch]'`
+
+for i in $files; do
+    tp=`wc "$i" | tr -s " " | cut -f 2 -d" "`
+    let nc="$nc + $tp"
+done
+
+echo $nc
+echo $np

Modified: trunk/Lib/sandbox/pyem/pyem/__init__.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/__init__.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/__init__.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,7 +1,12 @@
 #! /usr/bin/env python
-# Last Change: Thu Aug 24 07:00 PM 2006 J
+# Last Change: Thu Sep 28 01:00 PM 2006 J
 
+from info import __doc__
+
+from numpy.testing import NumpyTest
+test = NumpyTest().test
+
 version = '0.5.2'
 
 from gauss_mix import GmParamError, GM
-from gmm_em import GmmParamError, GMM
+from gmm_em import GmmParamError, GMM, EM

Modified: trunk/Lib/sandbox/pyem/pyem/_c_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/_c_densities.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/_c_densities.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,7 +1,7 @@
 #! /usr/bin/python
 #
 # Copyrighted David Cournapeau
-# Last Change: Mon Aug 28 06:00 PM 2006 J
+# Last Change: Tue Oct 03 05:00 PM 2006 J
 
 # This module uses a C implementation through ctypes, for diagonal cases
 # TODO:
@@ -14,6 +14,27 @@
 from scipy.stats import chi2
 import densities as D
 
+import ctypes
+from ctypes import cdll, c_uint, c_int, c_double, POINTER
+from numpy.ctypeslib import ndpointer, load_library
+
+ctypes_major    = int(ctypes.__version__.split('.')[0])
+if ctypes_major < 1:
+    msg =  "version of ctypes is %s, expected at least %s" \
+            % (ctypes.__version__, '1.0.0')
+    raise Exception(msg)
+
+# Requirements for diag gden
+_gden   = load_library('c_gden.so', __file__)
+arg1    = ndpointer(dtype='<f8')
+arg2    = c_uint
+arg3    = c_uint
+arg4    = ndpointer(dtype='<f8')
+arg5    = ndpointer(dtype='<f8')
+arg6    = ndpointer(dtype='<f8')
+_gden.gden_diag.argtypes    = [arg1, arg2, arg3, arg4, arg5, arg6]
+_gden.gden_diag.restype     = c_int
+
 # Error classes
 class DenError(Exception):
     """Base class for exceptions in this module.
@@ -120,24 +141,34 @@
     d   = mu.size
     n   = x.shape[0]
     if not log:
-        from ctypes import cdll, c_uint, c_int, c_double, POINTER
-        _gden   = cdll.LoadLibrary('src/libgden.so')
-        _gden.gden_diag.restype     = c_int
-        _gden.gden_diag.argtypes    = [POINTER(c_double), c_uint, c_uint,
-                POINTER(c_double), POINTER(c_double), POINTER(c_double)]
+        y       = N.zeros(n)
+        vat     = va.copy()
+        # _gden.gden_diag(N.require(x, requirements = 'C'), n, d, 
+        #         N.require(mu, requirements = 'C'),
+        #         N.require(inva, requirements = 'C'),
+        #         N.require(y, requirements = 'C'))
+        x       = N.require(x, requirements = 'C')
+        mu      = N.require(mu, requirements = 'C')
+        vat     = N.require(vat, requirements = 'C')
+        y       = N.require(y, requirements = 'C')
+        _gden.gden_diag(x, n, d, mu, vat, y)
+        return y
+        # _gden.gden_diag.restype     = c_int
+        # _gden.gden_diag.argtypes    = [POINTER(c_double), c_uint, c_uint,
+        #         POINTER(c_double), POINTER(c_double), POINTER(c_double)]
 
-        y   = N.zeros(n)
-        inva= 1/va
-        _gden.gden_diag(x.ctypes.data_as(POINTER(c_double)),
-            n, d,
-            mu.ctypes.data_as(POINTER(c_double)),
-            inva.ctypes.data_as(POINTER(c_double)),
-            y.ctypes.data_as(POINTER(c_double)))
+        # y   = N.zeros(n)
+        # inva= 1/va
+        # _gden.gden_diag(x.ctypes.data_as(POINTER(c_double)),
+        #     n, d,
+        #     mu.ctypes.data_as(POINTER(c_double)),
+        #     inva.ctypes.data_as(POINTER(c_double)),
+        #     y.ctypes.data_as(POINTER(c_double)))
     else:
         y   = _scalar_gauss_den(x[:,0], mu[0,0], va[0,0], log)
         for i in range(1, d):
             y    +=  _scalar_gauss_den(x[:,i], mu[0,i], va[0,i], log)
-    return y
+        return y
 
 def _full_gauss_den(x, mu, va, log):
     """ This function is the actual implementation

Modified: trunk/Lib/sandbox/pyem/pyem/densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/densities.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/densities.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,7 +1,7 @@
 #! /usr/bin/python
 #
 # Copyrighted David Cournapeau
-# Last Change: Wed Aug 30 12:00 PM 2006 J
+# Last Change: Fri Sep 29 06:00 PM 2006 J
 
 import numpy as N
 import numpy.linalg as lin

Modified: trunk/Lib/sandbox/pyem/pyem/kmean.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/kmean.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/kmean.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,6 +1,9 @@
 # /usr/bin/python
-# Last Change: Tue Aug 29 12:00 PM 2006 J
+# Last Change: Thu Sep 28 01:00 PM 2006 J
 
+#TODO:
+#   - a demo for kmeans
+
 import numpy as N
 
 def _py_vq(data, code):
@@ -69,39 +72,5 @@
 
     return code, label
 
-# Test functions usable for now
-def test_kmean():
-    X   = N.array([[3.0, 3], [4, 3], [4, 2],
-            [9, 2], [5, 1], [6, 2], [9, 4], 
-            [5, 2], [5, 4], [7, 4], [6, 5]])
-
-    initc   = N.concatenate(([[X[0]], [X[1]], [X[2]]])) 
-
-    codet1  = N.array([[3.0000, 3.0000],
-            [6.2000, 4.0000], 
-            [5.8000, 1.8000]])
-            
-    codet2  = N.array([[11.0/3, 8.0/3], 
-            [6.7500, 4.2500],
-            [6.2500, 1.7500]])
-
-    code    = initc.copy()
-
-    code1   = kmean(X, code, 1)[0]
-    code2   = kmean(X, code, 2)[0]
-
-    import numpy.testing as testing
-    try:
-        testing.assert_array_almost_equal(code1, codet1)
-        print "kmean test 1 succeded"
-    except AssertionError:
-        print "kmean test 1 failed"
-
-    try:
-        testing.assert_array_almost_equal(code2, codet2)
-        print "kmean test 2 succeded"
-    except AssertionError:
-        print "kmean test 2 failed"
-
 if __name__ == "__main__":
-    test_kmean()
+    pass

Added: trunk/Lib/sandbox/pyem/pyem/profile_data/profile_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/profile_data/profile_densities.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/profile_data/profile_densities.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -0,0 +1,54 @@
+import numpy as N
+from numpy.random import randn
+from pyem import densities as D
+from pyem import _c_densities as DC
+import tables
+
+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())
+
+    X   = randn(n, d)
+    for i in range(niter):
+        Y   = func(X, mu, va)
+
+def benchpy():
+    bench(D.gauss_den)
+
+def benchc():
+    bench(DC.gauss_den)
+
+def benchpyfull():
+    bench(D.gauss_den, 'full')
+
+def benchcfull():
+    bench(DC.gauss_den, 'full')
+
+if __name__ == "__main__":
+    import hotshot, hotshot.stats
+    profile_file    = 'gdenpy.prof'
+    prof    = hotshot.Profile(profile_file, lineevents=1)
+    prof.runcall(benchpy)
+    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(benchc)
+    p = hotshot.stats.load(profile_file)
+    print p.sort_stats('cumulative').print_stats(20)
+    prof.close()

Added: trunk/Lib/sandbox/pyem/pyem/profile_data/profile_gmm.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/profile_data/profile_gmm.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/profile_data/profile_gmm.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -0,0 +1,57 @@
+import numpy as N
+from pyem import GM, GMM
+import copy
+
+def bench1(mode = 'diag'):
+    #===========================================
+    # GMM of 20 comp, 20 dimension, 1e4 frames
+    #===========================================
+    d       = 15
+    k       = 30
+    nframes = 1e5
+    niter   = 10
+    mode    = 'diag'
+
+    #+++++++++++++++++++++++++++++++++++++++++++
+    # Create an artificial GMM model, samples it
+    #+++++++++++++++++++++++++++++++++++++++++++
+    print "Generating the mixture"
+    # Generate a model with k components, d dimensions
+    w, mu, va   = GM.gen_param(d, k, mode, spread = 3)
+    # gm          = GM(d, k, mode)
+    # gm.set_param(w, mu, va)
+    gm          = GM.fromvalues(w, mu, va)
+
+    # Sample nframes frames  from the model
+    data    = gm.sample(nframes)
+
+    #++++++++++++++++++++++++
+    # Learn the model with EM
+    #++++++++++++++++++++++++
+
+    # Init the model
+    print "Init a model for learning, with kmean for initialization"
+    lgm = GM(d, k, mode)
+    gmm = GMM(lgm, 'kmean')
+    
+    gmm.init(data)
+    # Keep the initialized model for drawing
+    gm0 = copy.copy(lgm)
+
+    # The actual EM, with likelihood computation
+    like    = N.zeros(niter)
+
+    print "computing..."
+    for i in range(niter):
+        print "iteration %d" % i
+        g, tgd  = gmm.sufficient_statistics(data)
+        like[i] = N.sum(N.log(N.sum(tgd, 1)))
+        gmm.update_em(data, g)
+
+if __name__ == "__main__":
+    import profile
+    profile.run('bench1()', 'gmmprof')
+    import pstats
+    p = pstats.Stats('gmmprof')
+    print p.sort_stats('cumulative').print_stats(20)
+

Deleted: trunk/Lib/sandbox/pyem/pyem/profile_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/profile_densities.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/profile_densities.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,63 +0,0 @@
-import numpy as N
-from numpy.random import randn
-import densities as D
-import _c_densities as DC
-import tables
-
-def bench(func, mode = 'diag'):
-    #===========================================
-    # Diag Gaussian of dimension 20
-    #===========================================
-    d       = 30
-    n       = 1e5
-    niter   = 100
-
-    # 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())
-
-    X   = randn(n, d)
-    print "Compute %d times densities, %d dimension, %d frames" % (niter, d, n)
-    for i in range(niter):
-        Y   = func(X, mu, va)
-    
-    # Check values
-    h5file  = tables.openFile('diag.dat', "r")
-    X   = h5file.root.input.read()
-    mu  = h5file.root.mu.read()
-    va  = h5file.root.va.read()
-    Yt  = h5file.root.output.read()
-    Y   = func(X, mu, va)
-
-    try:
-        N.testing.assert_array_almost_equal(Y, Yt) 
-    except AssertionError:
-        print N.sum(Y)
-        print N.sqrt(N.sum((Y-Yt) **2)) / n
-        raise "Not accurate !!!!"
-
-def benchpy():
-    bench(D.gauss_den)
-
-def benchc():
-    bench(DC.gauss_den)
-
-def benchpyfull():
-    bench(D.gauss_den, 'full')
-
-def benchcfull():
-    bench(DC.gauss_den, 'full')
-
-if __name__ == "__main__":
-    import profile
-    import pstats
-    profile.run('benchpy()', 'gdenprof')
-    p = pstats.Stats('gdenprof')
-    print p.sort_stats('cumulative').print_stats(20)
-    profile.run('benchc()', 'gdenprof')
-    p = pstats.Stats('gdenprof')
-    print p.sort_stats('cumulative').print_stats(20)

Deleted: trunk/Lib/sandbox/pyem/pyem/profile_gmm.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/profile_gmm.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/profile_gmm.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,56 +0,0 @@
-import numpy as N
-from gmm_em import GM, GMM
-import copy
-
-def bench1(mode = 'diag'):
-    #===========================================
-    # GMM of 20 comp, 20 dimension, 1e4 frames
-    #===========================================
-    d       = 20
-    k       = 20
-    nframes = 1e4
-    niter   = 10
-    mode    = 'full'
-
-    #+++++++++++++++++++++++++++++++++++++++++++
-    # Create an artificial GMM model, samples it
-    #+++++++++++++++++++++++++++++++++++++++++++
-    print "Generating the mixture"
-    # Generate a model with k components, d dimensions
-    w, mu, va   = GM.gen_param(d, k, mode, spread = 3)
-    gm          = GM(d, k, mode)
-    gm.set_param(w, mu, va)
-
-    # Sample nframes frames  from the model
-    data    = gm.sample(nframes)
-
-    #++++++++++++++++++++++++
-    # Learn the model with EM
-    #++++++++++++++++++++++++
-
-    # Init the model
-    print "Init a model for learning, with kmean for initialization"
-    lgm = GM(d, k, mode)
-    gmm = GMM(lgm, 'kmean')
-    
-    gmm.init(data)
-    # Keep the initialized model for drawing
-    gm0 = copy.copy(lgm)
-
-    # The actual EM, with likelihood computation
-    like    = N.zeros(niter)
-
-    print "computing..."
-    for i in range(niter):
-        print "iteration %d" % i
-        g, tgd  = gmm.sufficient_statistics(data)
-        like[i] = N.sum(N.log(N.sum(tgd, 1)))
-        gmm.update_em(data, g)
-
-if __name__ == "__main__":
-    import profile
-    profile.run('bench1()', 'gmmprof')
-    import pstats
-    p = pstats.Stats('gmmprof')
-    print p.sort_stats('cumulative').print_stats(20)
-

Modified: trunk/Lib/sandbox/pyem/pyem/src/Makefile
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/src/Makefile	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/src/Makefile	2006-10-12 13:48:08 UTC (rev 2281)
@@ -5,12 +5,12 @@
 
 PYTHONINC	= -I/usr/include/python2.4
 NUMPYINC	= -I/usr/lib/python2.4/site-packages/numpy/core/include
-OPTIMS		= -O3 -funroll-all-loops -ffast-math -march=pentium4
+OPTIMS		= -O3 -funroll-all-loops -march=pentium4 -msse2
 WARN		= -W -Wall
 
 CFLAGS	= $(PYTHONINC) $(NUMPYINC) $(OPTIMS) $(WARN)
 
-libgden.so: c_gden.o
+c_gden.so: c_gden.o
 	$(LD) -shared -o $@ $< -Wl,-soname,$@
 
 c_gden.o: c_gden.c

Modified: trunk/Lib/sandbox/pyem/pyem/src/c_gden.c
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/src/c_gden.c	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/src/c_gden.c	2006-10-12 13:48:08 UTC (rev 2281)
@@ -2,8 +2,11 @@
 #include <stdio.h>
 #include <math.h>
 
+/*
+ * use va for caching, so its content is modified
+ */
 int gden_diag(const double *in, const size_t n, const size_t d, 
-        const double* mu, const double* inva, double* out)
+        const double* mu, double* va, double* out)
 {
     size_t  j, i;
     double  fac = 1.0;
@@ -13,8 +16,9 @@
      * Cache some precomputing
      */
     for(j = 0; j < d; ++j) {
-        fac     *= sqrt(inva[j] / (2 * M_PI) );
-        //printf("inva[%d] is %f, fac %f\n", j, inva[j], fac);
+        va[j]   = 0.5/va[j];
+        fac     *= sqrt(va[j] / (M_PI) );
+        va[j]   *= -1.0;
     }
 
     /*
@@ -23,11 +27,11 @@
     for(i = 0; i < n; ++i) {
         acc = 0;
         for(j = 0; j < d; ++j) {
-            tmp = (in[d * i + j] - mu[j]);
-            acc += tmp * tmp * -0.5 * inva[j];
+            tmp = (in[i *d + j] - mu[j]);
+            acc += tmp * tmp * va[j];
         }
         out[i]  = exp(acc) * fac;
     }
-
+ 
     return 0;
 }

Modified: trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,145 +1,181 @@
-def generate_test_data(n, d, mode = 'diag', file='test.dat'):
-    """Generate a set of data of dimension d, with n frames,
-    that is input data, mean, var and output of gden, so that
-    other implementations can be tested against"""
-    mu  = randn(1, d)
-    if mode == 'diag':
-        va  = abs(randn(1, d))
-    elif mode == 'full':
-        va  = randn(d, d)
-        va  = dot(va, va.transpose())
+#! /usr/bin/env python
+# Last Change: Fri Sep 29 06:00 PM 2006 J
 
-    input   = randn(n, d)
-    output  = gauss_den(input, mu, va)
+import sys
+from numpy.testing import *
 
-    import tables
-    h5file  = tables.openFile(file, "w")
+import numpy as N
 
-    h5file.createArray(h5file.root, 'input', input)
-    h5file.createArray(h5file.root, 'mu', mu)
-    h5file.createArray(h5file.root, 'va', va)
-    h5file.createArray(h5file.root, 'output', output)
+set_package_path()
+from pyem.densities import gauss_den
+from pyem._c_densities import gauss_den as c_gauss_den
 
-    h5file.close()
+restore_path()
 
-def test_gauss_den():
-    """"""
-    # import tables
-    # import numpy as N
-    # 
-    # filename    = 'dendata.h5'
+#Optional:
+set_local_path()
+# import modules that are located in the same directory as this file.
+restore_path()
 
-    # # # Dimension 1
-    # # d   = 1
-    # # mu  = 1.0
-    # # va  = 2.0
+class test_densities(NumpyTestCase):
+    def _generate_test_data_1d(self):
+        self.va     = 2.0
+        self.mu     = 1.0
+        self.X      = N.linspace(-2, 2, 10)[:, N.newaxis]
 
-    # # X   = randn(1e3, 1)
+        self.Yt     = N.array([0.02973257230591, 0.05512079811082, 0.09257745306945, 
+                0.14086453882683,
+                0.19418015562214, 0.24250166773127, 0.27436665745048, 0.28122547107069,
+                0.26114678964743, 0.21969564473386])
 
-    # # Y   = gauss_den(X, mu, va)
+    def _generate_test_data_2d_diag(self):
+        #============================
+        # Small test in 2d (diagonal)
+        #============================
+        self.mu  = N.atleast_2d([-1.0, 2.0])
+        self.va  = N.atleast_2d([2.0, 3.0])
+        
+        self.X  = N.zeros((10, 2))
+        self.X[:,0] = N.linspace(-2, 2, 10)
+        self.X[:,1] = N.linspace(-1, 3, 10)
 
-    # # h5file      = tables.openFile(filename, "w")
+        self.Yt  = N.array([0.01129091565384, 0.02025416837152, 0.03081845516786, 
+                0.03977576221540, 0.04354490552910, 0.04043592581117, 
+                0.03184994053539, 0.02127948225225, 0.01205937178755, 
+                0.00579694938623 ])
 
-    # # h5file.createArray(h5file.root, 'X', X)
-    # # h5file.createArray(h5file.root, 'mu', mu)
-    # # h5file.createArray(h5file.root, 'va', va)
-    # # h5file.createArray(h5file.root, 'Y', Y)
 
-    # # h5file.close()
+    def _generate_test_data_2d_full(self):
+        #============================
+        # Small test in 2d (full mat)
+        #============================
+        self.mu = N.array([[0.2, -1.0]])
+        self.va = N.array([[1.2, 0.1], [0.1, 0.5]])
+        X1      = N.linspace(-2, 2, 10)[:, N.newaxis]
+        X2      = N.linspace(-3, 3, 10)[:, N.newaxis]
+        self.X  = N.concatenate(([X1, X2]), 1)
+        
+        self.Yt = N.array([0.00096157109751, 0.01368908714856,
+            0.07380823191162, 0.15072050533842, 
+            0.11656739937861, 0.03414436965525,
+            0.00378789836599, 0.00015915297541, 
+            0.00000253261067, 0.00000001526368])
 
-    # # # Dimension 2, diag
-    # # d   = 2
-    # # mu  = N.array([1.0, -2.0])
-    # # va  = N.array([1.0, 2.0])
+    def _check_py(self, level, decimal = 12):
+        Y   = gauss_den(self.X, self.mu, self.va)
+        assert_array_almost_equal(Y, self.Yt, decimal)
 
-    # # X   = randn(1e3, 2)
+    def _check_c(self, level, decimal = 12):
+        Y   = c_gauss_den(self.X, self.mu, self.va)
+        assert_array_almost_equal(Y, self.Yt, decimal)
 
-    # # Y   = gauss_den(X, mu, va)
+    def check_py_1d(self, level = 1):
+        self._generate_test_data_1d()
+        self._check_py(level)
 
-    # # h5file      = tables.openFile(filename, "w")
+    def check_py_2d_diag(self, level = 1):
+        self._generate_test_data_2d_diag()
+        self._check_py(level)
 
-    # # h5file.createArray(h5file.root, 'X', X)
-    # # h5file.createArray(h5file.root, 'mu', mu)
-    # # h5file.createArray(h5file.root, 'va', va)
-    # # h5file.createArray(h5file.root, 'Y', Y)
+    def check_py_2d_full(self, level = 1):
+        self._generate_test_data_2d_full()
+        self._check_py(level)
 
-    # # Dimension 2, full
-    # d   = 2
-    # mu  = N.array([[0.2, -1.0]])
-    # va  = N.array([[1.2, 0.1], [0.1, 0.5]])
+    def check_c_1d(self, level = 1):
+        self._generate_test_data_1d()
+        self._check_c(level)
 
-    # X   = randn(1e3, 2)
+    def check_c_2d_diag(self, level = 1):
+        self._generate_test_data_2d_diag()
+        self._check_c(level)
 
-    # Y   = gauss_den(X, mu, va)
+    def check_c_2d_full(self, level = 1):
+        self._generate_test_data_2d_full()
+        self._check_c(level)
 
-    # h5file      = tables.openFile(filename, "w")
+if __name__ == "__main__":
+    NumpyTest().run()
 
-    # h5file.createArray(h5file.root, 'X', X)
-    # h5file.createArray(h5file.root, 'mu', mu)
-    # h5file.createArray(h5file.root, 'va', va)
-    # h5file.createArray(h5file.root, 'Y', Y)
-
-    # h5file.close()
-
-    import numpy.testing as testing
-    #=================
-    # Small test in 1d
-    #=================
-    va  = 2.0
-    mu  = 1.0
-    X   = N.linspace(-2, 2, 10)[:, N.NewAxis]
-
-    Yt  = N.array([0.02973257230591, 0.05512079811082, 0.09257745306945, 
-            0.14086453882683,
-            0.19418015562214, 0.24250166773127, 0.27436665745048, 0.28122547107069,
-            0.26114678964743, 0.21969564473386])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "1d test succeded"
-    except AssertionError:
-        print "test fails in 1d"
-
-    #============================
-    # Small test in 2d (diagonal)
-    #============================
-    mu  = N.atleast_2d([-1.0, 2.0])
-    va  = N.atleast_2d([2.0, 3.0])
-    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
-    X2  = N.linspace(-1, 3, 10)[:, N.NewAxis]
-    X   = N.concatenate(([X1, X2]), 1)
-    
-    Yt  = N.array([0.01129091565384, 0.02025416837152, 0.03081845516786, 
-            0.03977576221540, 0.04354490552910, 0.04043592581117, 
-            0.03184994053539, 0.02127948225225, 0.01205937178755, 
-            0.00579694938623 ])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "2d diag test succeded"
-    except AssertionError:
-        print "test fails in 2d diag"
-
-    #============================
-    # Small test in 2d (full mat)
-    #============================
-    mu  = N.array([[0.2, -1.0]])
-    va  = N.array([[1.2, 0.1], [0.1, 0.5]])
-    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
-    X2  = N.linspace(-3, 3, 10)[:, N.NewAxis]
-    X   = N.concatenate(([X1, X2]), 1)
-    
-    Yt  = N.array([0.00096157109751, 0.01368908714856,
-        0.07380823191162, 0.15072050533842, 
-        0.11656739937861, 0.03414436965525,
-        0.00378789836599, 0.00015915297541, 
-        0.00000253261067, 0.00000001526368])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "2d full test succeded"
-    except AssertionError:
-        print "test fails in 2d full"
+# def generate_test_data(n, d, mode = 'diag', file='test.dat'):
+#     """Generate a set of data of dimension d, with n frames,
+#     that is input data, mean, var and output of gden, so that
+#     other implementations can be tested against"""
+#     mu  = randn(1, d)
+#     if mode == 'diag':
+#         va  = abs(randn(1, d))
+#     elif mode == 'full':
+#         va  = randn(d, d)
+#         va  = dot(va, va.transpose())
+# 
+#     input   = randn(n, d)
+#     output  = gauss_den(input, mu, va)
+# 
+#     import tables
+#     h5file  = tables.openFile(file, "w")
+# 
+#     h5file.createArray(h5file.root, 'input', input)
+#     h5file.createArray(h5file.root, 'mu', mu)
+#     h5file.createArray(h5file.root, 'va', va)
+#     h5file.createArray(h5file.root, 'output', output)
+# 
+#     h5file.close()
+# 
+# def test_gauss_den():
+#     """"""
+#     # import tables
+#     # import numpy as N
+#     # 
+#     # filename    = 'dendata.h5'
+# 
+#     # # # Dimension 1
+#     # # d   = 1
+#     # # mu  = 1.0
+#     # # va  = 2.0
+# 
+#     # # X   = randn(1e3, 1)
+# 
+#     # # Y   = gauss_den(X, mu, va)
+# 
+#     # # h5file      = tables.openFile(filename, "w")
+# 
+#     # # h5file.createArray(h5file.root, 'X', X)
+#     # # h5file.createArray(h5file.root, 'mu', mu)
+#     # # h5file.createArray(h5file.root, 'va', va)
+#     # # h5file.createArray(h5file.root, 'Y', Y)
+# 
+#     # # h5file.close()
+# 
+#     # # # Dimension 2, diag
+#     # # d   = 2
+#     # # mu  = N.array([1.0, -2.0])
+#     # # va  = N.array([1.0, 2.0])
+# 
+#     # # X   = randn(1e3, 2)
+# 
+#     # # Y   = gauss_den(X, mu, va)
+# 
+#     # # h5file      = tables.openFile(filename, "w")
+# 
+#     # # h5file.createArray(h5file.root, 'X', X)
+#     # # h5file.createArray(h5file.root, 'mu', mu)
+#     # # h5file.createArray(h5file.root, 'va', va)
+#     # # h5file.createArray(h5file.root, 'Y', Y)
+# 
+#     # # Dimension 2, full
+#     # d   = 2
+#     # mu  = N.array([[0.2, -1.0]])
+#     # va  = N.array([[1.2, 0.1], [0.1, 0.5]])
+# 
+#     # X   = randn(1e3, 2)
+# 
+#     # Y   = gauss_den(X, mu, va)
+# 
+#     # h5file      = tables.openFile(filename, "w")
+# 
+#     # h5file.createArray(h5file.root, 'X', X)
+#     # h5file.createArray(h5file.root, 'mu', mu)
+#     # h5file.createArray(h5file.root, 'va', va)
+#     # h5file.createArray(h5file.root, 'Y', Y)
+# 
+#     # h5file.close()
+# 

Added: trunk/Lib/sandbox/pyem/pyem/tests/test_kmean.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/tests/test_kmean.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/pyem/tests/test_kmean.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -0,0 +1,46 @@
+#! /usr/bin/env python
+# Last Change: Thu Sep 28 01:00 PM 2006 J
+
+import sys
+from numpy.testing import *
+
+import numpy as N
+
+set_package_path()
+from pyem.kmean import kmean
+restore_path()
+
+#Optional:
+set_local_path()
+# import modules that are located in the same directory as this file.
+restore_path()
+
+# Global data
+X   = N.array([[3.0, 3], [4, 3], [4, 2],
+        [9, 2], [5, 1], [6, 2], [9, 4], 
+        [5, 2], [5, 4], [7, 4], [6, 5]])
+
+codet1  = N.array([[3.0000, 3.0000],
+        [6.2000, 4.0000], 
+        [5.8000, 1.8000]])
+        
+codet2  = N.array([[11.0/3, 8.0/3], 
+        [6.7500, 4.2500],
+        [6.2500, 1.7500]])
+
+class test_kmean(NumpyTestCase):
+    def check_iter1(self, level=1):
+        initc   = N.concatenate(([[X[0]], [X[1]], [X[2]]])) 
+        code    = initc.copy()
+        code1   = kmean(X, code, 1)[0]
+
+        assert_array_almost_equal(code1, codet1)
+    def check_iter2(self, level=1):
+        initc   = N.concatenate(([[X[0]], [X[1]], [X[2]]])) 
+        code    = initc.copy()
+        code2   = kmean(X, code, 2)[0]
+
+        assert_array_almost_equal(code2, codet2)
+
+if __name__ == "__main__":
+    NumpyTest().run()

Modified: trunk/Lib/sandbox/pyem/setup.py
===================================================================
--- trunk/Lib/sandbox/pyem/setup.py	2006-10-12 13:47:57 UTC (rev 2280)
+++ trunk/Lib/sandbox/pyem/setup.py	2006-10-12 13:48:08 UTC (rev 2281)
@@ -1,5 +1,5 @@
 #! /usr/bin/env python
-# Last Change: Mon Aug 21 12:00 PM 2006 J
+# Last Change: Fri Sep 29 07:00 PM 2006 J
 # TODO:
 #   - check how to handle cmd line build options with distutils and use
 #   it in the building process
@@ -12,6 +12,7 @@
 # distutils does not update MANIFEST correctly, removes it
 import os
 if os.path.exists('MANIFEST'): os.remove('MANIFEST')
+from os.path import join
 
 import re
 
@@ -59,7 +60,8 @@
 class SetupOption:
     def __init__(self):
         self.kmean      = 'py'
-        self.ext_modules= []
+        self.ext_modules= [Extension(join('pyem', 'c_gden'),
+                              sources=[join('pyem', 'src', 'c_gden.c')]) ]
         self.cmdclass   = {}
         self.subsdic     = {'%KMEANIMPORT%': []}
 
@@ -95,7 +97,7 @@
             #self.subsdic['%KMEANIMPORT%']   = pyrex_kmean
     def setup(self):
         self._config_kmean()
-        import time
+        #import time
         #do_subst_in_file('pyem/kmean.py.in', 'pyem/kmean.py', self.subsdic)
         setup(name      = DISTNAME,
             version     = VERSION,
@@ -103,7 +105,7 @@
             author      = AUTHOR,
             author_email= AUTHOR_EMAIL,
             url         = URL,
-            packages    = ['pyem'],
+            packages    = ['pyem', 'pyem.tests', 'pyem.profile_data'],
             ext_modules = self.ext_modules,
             cmdclass    = self.cmdclass)
 



More information about the Scipy-svn mailing list