← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pyem/src/c_gmm.pyx

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pyem/src/c_gmm.pyx
1
 
# Last Change: Thu Jul 13 04:00 PM 2006 J
2
 
 
3
 
cimport c_numpy
4
 
cimport c_python
5
 
import numpy
6
 
 
7
 
c_numpy.import_array()
8
 
 
9
 
# pyrex version of _vq. Much faster in high dimension/high number of K 
10
 
# (ie K > 3-4)
11
 
def _vq(data, init):
12
 
    (n, d)  = data.shape
13
 
    label   = numpy.zeros(n, int)
14
 
    _imp_vq(data, init, label)
15
 
 
16
 
    return label
17
 
 
18
 
def _imp_vq(c_numpy.ndarray data, c_numpy.ndarray init, c_numpy.ndarray label):
19
 
    cdef int n
20
 
    cdef int d
21
 
    cdef int nc
22
 
    cdef int i
23
 
    cdef int j
24
 
    cdef int k
25
 
    cdef int *labeld
26
 
    cdef double *da, *code
27
 
    cdef double dist
28
 
    cdef double acc
29
 
 
30
 
    n   = data.dimensions[0]
31
 
    d   = data.dimensions[1]
32
 
    nc  = init.dimensions[0]
33
 
 
34
 
    if not data.dtype == numpy.dtype(numpy.float64):
35
 
        print '_vq not (yet) implemented for dtype %s'%dtype.name
36
 
        return
37
 
    da  = (<double*>data.data)
38
 
 
39
 
    if not init.dtype == numpy.dtype(numpy.float64):
40
 
        print '_vq not (yet) implemented for dtype %s'%dtype.name
41
 
        return
42
 
    code    = (<double*>init.data)
43
 
 
44
 
    if not label.dtype == numpy.dtype(numpy.int32):
45
 
        print '_vq not (yet) implemented for dtype %s'%dtype.name
46
 
        return
47
 
    labeld  = (<int*>label.data)
48
 
 
49
 
    for i from 0<=i<n:
50
 
        acc = 0
51
 
        lab = 0
52
 
        for j from 0<=j<d:
53
 
            acc = acc + (da[i * d + j] - code[j]) * \
54
 
                (da[i * d + j] - code[j])
55
 
        dist    = acc
56
 
        for k from 1<=k<nc:
57
 
            acc     = 0
58
 
            for j from 0<=j<d:
59
 
                acc = acc + (da[i * d + j] - code[k * d + j]) * \
60
 
                    (da[i * d + j] - code[k * d + j])
61
 
            if acc < dist:
62
 
                dist    = acc
63
 
                lab     = k
64
 
        labeld[i]   = lab
65
 
 
66
 
    return lab