← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pyem/_c_densities.py

  • 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/_c_densities.py
1
 
#! /usr/bin/python
2
 
#
3
 
# Copyrighted David Cournapeau
4
 
# Last Change: Thu Nov 09 05:00 PM 2006 J
5
 
 
6
 
# This module uses a C implementation through ctypes, for diagonal cases
7
 
# TODO:
8
 
#   - portable way to find/open the shared library
9
 
#   - full cov matrice
10
 
 
11
 
import numpy as N
12
 
import numpy.linalg as lin
13
 
from numpy.random import randn
14
 
from scipy.stats import chi2
15
 
import densities as D
16
 
 
17
 
import ctypes
18
 
from ctypes import cdll, c_uint, c_int, c_double, POINTER
19
 
from numpy.ctypeslib import ndpointer, load_library
20
 
 
21
 
ctypes_major    = int(ctypes.__version__.split('.')[0])
22
 
if ctypes_major < 1:
23
 
    msg =  "version of ctypes is %s, expected at least %s" \
24
 
            % (ctypes.__version__, '1.0.0')
25
 
    raise ImportError(msg)
26
 
 
27
 
# Requirements for diag gden
28
 
_gden   = load_library('c_gden.so', __file__)
29
 
arg1    = ndpointer(dtype=N.float64)
30
 
arg2    = c_uint
31
 
arg3    = c_uint
32
 
arg4    = ndpointer(dtype=N.float64)
33
 
arg5    = ndpointer(dtype=N.float64)
34
 
arg6    = ndpointer(dtype=N.float64)
35
 
_gden.gden_diag.argtypes    = [arg1, arg2, arg3, arg4, arg5, arg6]
36
 
_gden.gden_diag.restype     = c_int
37
 
 
38
 
# Error classes
39
 
class DenError(Exception):
40
 
    """Base class for exceptions in this module.
41
 
    
42
 
    Attributes:
43
 
        expression -- input expression in which the error occurred
44
 
        message -- explanation of the error"""
45
 
    def __init__(self, message):
46
 
        self.message    = message
47
 
    
48
 
    def __str__(self):
49
 
        return self.message
50
 
 
51
 
# The following function do all the fancy stuff to check that parameters
52
 
# are Ok, and call the right implementation if args are OK.
53
 
def gauss_den(x, mu, va, log = False):
54
 
    """ Compute multivariate Gaussian density at points x for 
55
 
    mean mu and variance va.
56
 
    
57
 
    Vector are row vectors, except va which can be a matrix
58
 
    (row vector variance for diagonal variance)
59
 
    
60
 
    If log is True, than the log density is returned 
61
 
    (useful for underflow ?)"""
62
 
    mu  = N.atleast_2d(mu)
63
 
    va  = N.atleast_2d(va)
64
 
    x   = N.atleast_2d(x)
65
 
    
66
 
    #=======================#
67
 
    # Checking parameters   #
68
 
    #=======================#
69
 
    if len(N.shape(mu)) != 2:
70
 
        raise DenError("mu is not rank 2")
71
 
        
72
 
    if len(N.shape(va)) != 2:
73
 
        raise DenError("va is not rank 2")
74
 
        
75
 
    if len(N.shape(x)) != 2:
76
 
        raise DenError("x is not rank 2")
77
 
        
78
 
    (n, d)      = x.shape
79
 
    (dm0, dm1)  = mu.shape
80
 
    (dv0, dv1)  = va.shape
81
 
    
82
 
    # Check x and mu same dimension
83
 
    if dm0 != 1:
84
 
        msg = "mean must be a row vector!"
85
 
        raise DenError(msg)
86
 
    if dm1 != d:
87
 
        msg = "x and mu not same dim"
88
 
        raise DenError(msg)
89
 
    # Check va and mu same size
90
 
    if dv1 != d:
91
 
        msg = "mu and va not same dim"
92
 
        raise DenError(msg)
93
 
    if dv0 != 1 and dv0 != d:
94
 
        msg = "va not square"
95
 
        raise DenError(msg)
96
 
 
97
 
    #===============#
98
 
    # Computation   #
99
 
    #===============#
100
 
    if d == 1:
101
 
        # scalar case
102
 
        return _scalar_gauss_den(x[:, 0], mu[0, 0], va[0, 0], log)
103
 
    elif dv0 == 1:
104
 
        # Diagonal matrix case
105
 
        return _diag_gauss_den(x, mu, va, log)
106
 
    elif dv1 == dv0:
107
 
        # full case
108
 
        return  _full_gauss_den(x, mu, va, log)
109
 
    else:
110
 
        raise DenError("variance mode not recognized, this is a bug")
111
 
 
112
 
# Those 3 functions do almost all the actual computation
113
 
def _scalar_gauss_den(x, mu, va, log):
114
 
    """ This function is the actual implementation
115
 
    of gaussian pdf in scalar case. It assumes all args
116
 
    are conformant, so it should not be used directly
117
 
    
118
 
    ** Expect centered data (ie with mean removed) **
119
 
 
120
 
    Call gauss_den instead"""
121
 
    d       = mu.size
122
 
    inva    = 1/va
123
 
    fac     = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
124
 
    y       = ((x-mu) ** 2) * -0.5 * inva
125
 
    if not log:
126
 
        y   = fac * N.exp(y)
127
 
    else:
128
 
        y   = y + log(fac)
129
 
 
130
 
    return y
131
 
    
132
 
def _diag_gauss_den(x, mu, va, log):
133
 
    """ This function is the actual implementation
134
 
    of gaussian pdf in scalar case. It assumes all args
135
 
    are conformant, so it should not be used directly
136
 
    
137
 
    ** Expect centered data (ie with mean removed) **
138
 
 
139
 
    Call gauss_den instead"""
140
 
    # Diagonal matrix case
141
 
    d   = mu.size
142
 
    n   = x.shape[0]
143
 
    if not log:
144
 
        y       = N.zeros(n)
145
 
        vat     = va.copy()
146
 
        # _gden.gden_diag(N.require(x, requirements = 'C'), n, d, 
147
 
        #         N.require(mu, requirements = 'C'),
148
 
        #         N.require(inva, requirements = 'C'),
149
 
        #         N.require(y, requirements = 'C'))
150
 
        x       = N.require(x, requirements = 'C')
151
 
        mu      = N.require(mu, requirements = 'C')
152
 
        vat     = N.require(vat, requirements = 'C')
153
 
        y       = N.require(y, requirements = 'C')
154
 
        _gden.gden_diag(x, n, d, mu, vat, y)
155
 
        return y
156
 
        # _gden.gden_diag.restype     = c_int
157
 
        # _gden.gden_diag.argtypes    = [POINTER(c_double), c_uint, c_uint,
158
 
        #         POINTER(c_double), POINTER(c_double), POINTER(c_double)]
159
 
 
160
 
        # y   = N.zeros(n)
161
 
        # inva= 1/va
162
 
        # _gden.gden_diag(x.ctypes.data_as(POINTER(c_double)),
163
 
        #     n, d,
164
 
        #     mu.ctypes.data_as(POINTER(c_double)),
165
 
        #     inva.ctypes.data_as(POINTER(c_double)),
166
 
        #     y.ctypes.data_as(POINTER(c_double)))
167
 
    else:
168
 
        y   = _scalar_gauss_den(x[:,0], mu[0,0], va[0,0], log)
169
 
        for i in range(1, d):
170
 
            y    +=  _scalar_gauss_den(x[:,i], mu[0,i], va[0,i], log)
171
 
        return y
172
 
 
173
 
def _full_gauss_den(x, mu, va, log):
174
 
    """ This function is the actual implementation
175
 
    of gaussian pdf in full matrix case. 
176
 
    
177
 
    It assumes all args are conformant, so it should 
178
 
    not be used directly Call gauss_den instead
179
 
    
180
 
    ** Expect centered data (ie with mean removed) **
181
 
 
182
 
    Does not check if va is definite positive (on inversible 
183
 
    for that matter), so the inverse computation and/or determinant
184
 
    would throw an exception."""
185
 
    d       = mu.size
186
 
    inva    = lin.inv(va)
187
 
    fac     = 1 / N.sqrt( (2*N.pi) ** d * N.fabs(lin.det(va)))
188
 
 
189
 
    # we are using a trick with sum to "emulate" 
190
 
    # the matrix multiplication inva * x without any explicit loop
191
 
    y   = N.dot((x-mu), inva)
192
 
    y   = -0.5 * N.sum(y * (x-mu), 1)
193
 
 
194
 
    if not log:
195
 
        y   = fac * N.exp(y)
196
 
    else:
197
 
        y   = y + N.log(fac)
198
 
 
199
 
    return y
200
 
 
201
 
if __name__ == "__main__":
202
 
    #=========================================
203
 
    # Test accuracy between pure and C python
204
 
    #=========================================
205
 
    mu  = N.array([2.0, 3])
206
 
    va  = N.array([5.0, 3])
207
 
 
208
 
    # Generate a multivariate gaussian of mean mu and covariance va
209
 
    nframes = 1e4
210
 
    X       = randn(nframes, 2)
211
 
    Yc      = N.dot(N.diag(N.sqrt(va)), X.transpose())
212
 
    Yc      = Yc.transpose() + mu
213
 
 
214
 
    Y   = D.gauss_den(Yc, mu, va)
215
 
    Yt  = gauss_den(Yc, mu, va)
216
 
 
217
 
    print "Diff is " + str(N.sqrt(N.sum((Y-Yt) ** 2))/nframes/2)