[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 09:48:25 EDT 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