← Back to branch summary

~jtaylor/ubuntu/precise/python-numpy/multiarch-fix-818867

« back to all changes in this revision

Viewing changes to numpy/random/mtrand/mtrand.pyx

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik, Riku Voipio, Tiziano Zito, Carlos Galisteo, Ondrej Certik
  • Date: 2008-07-08 15:08:16 UTC
  • mfrom: (0.1.21 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080708150816-ekf992jcp2k1eua3
Tags: 1:1.1.0-3
http://bugs.debian.org/489568
http://bugs.debian.org/489726
[ Riku Voipio ]
* debian/control: atlas is not available on armel, and after a quick look
  neither on alpha. I'd also suggest dropping
  libatlas-sse-dev|libatlas-sse2-dev|libatlas-3dnow-dev alternative combo
  away, these are potentially dangerous on buildd's. Ondrej: dropped.
  (Closes: #489568)

[ Tiziano Zito ]
* patch: build _dotblas.c when ATLAS is not installed, build-conflict with
  atlas, build-depend on blas+lapack only, as it used to be (Closes: #489726)

[ Carlos Galisteo ]
* debian/control
  - Added Homepage field.

[ Ondrej Certik ]
* Checked the package on i386 and amd64, both with and without atlas, all
  tests run and the numpy package is faster if atlas is around. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
numpy/random/mtrand/mtrand.pyx
9
9
# distribute, sublicense, and/or sell copies of the Software, and to
10
10
# permit persons to whom the Software is furnished to do so, subject to
11
11
# the following conditions:
12
 
# 
 
12
#
13
13
# The above copyright notice and this permission notice shall be included
14
14
# in all copies or substantial portions of the Software.
15
 
# 
 
15
#
16
16
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
17
17
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18
18
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
36
36
    ctypedef struct rk_state:
37
37
        unsigned long key[624]
38
38
        int pos
 
39
        int has_gauss
 
40
        double gauss
39
41
 
40
42
    ctypedef enum rk_error:
41
43
        RK_NOERR = 0
43
45
        RK_ERR_MAX = 2
44
46
 
45
47
    char *rk_strerror[2]
46
 
    
 
48
 
47
49
    # 0xFFFFFFFFUL
48
50
    unsigned long RK_MAX
49
51
 
61
63
    double rk_gauss(rk_state *state)
62
64
 
63
65
cdef extern from "distributions.h":
64
 
    
 
66
 
65
67
    double rk_normal(rk_state *state, double loc, double scale)
66
68
    double rk_standard_exponential(rk_state *state)
67
69
    double rk_exponential(rk_state *state, double scale)
86
88
    double rk_rayleigh(rk_state *state, double mode)
87
89
    double rk_wald(rk_state *state, double mean, double scale)
88
90
    double rk_triangular(rk_state *state, double left, double mode, double right)
89
 
    
 
91
 
90
92
    long rk_binomial(rk_state *state, long n, double p)
91
93
    long rk_binomial_btpe(rk_state *state, long n, double p)
92
94
    long rk_binomial_inversion(rk_state *state, long n, double p)
111
113
 
112
114
 
113
115
cdef extern from "initarray.h":
114
 
   void init_by_array(rk_state *self, unsigned long *init_key, 
 
116
   void init_by_array(rk_state *self, unsigned long *init_key,
115
117
                      unsigned long key_length)
116
118
 
117
119
# Initialize numpy
135
137
            array_data[i] = func(state)
136
138
        return array
137
139
 
 
140
 
 
141
cdef object cont1_array_sc(rk_state *state, rk_cont1 func, object size, double a):
 
142
    cdef double *array_data
 
143
    cdef ndarray array "arrayObject"
 
144
    cdef long length
 
145
    cdef long i
 
146
 
 
147
    if size is None:
 
148
        return func(state, a)
 
149
    else:
 
150
        array = <ndarray>_sp.empty(size, _sp.float64)
 
151
        length = PyArray_SIZE(array)
 
152
        array_data = <double *>array.data
 
153
        for i from 0 <= i < length:
 
154
            array_data[i] = func(state, a)
 
155
        return array
 
156
 
138
157
cdef object cont1_array(rk_state *state, rk_cont1 func, object size, ndarray oa):
139
158
    cdef double *array_data
140
159
    cdef double *oa_data
143
162
    cdef npy_intp i
144
163
    cdef flatiter itera
145
164
    cdef broadcast multi
146
 
    cdef int scalar
147
 
 
148
 
    scalar = 0
149
 
    if oa.nd == 0:
150
 
        oa_data = <double *>oa.data
151
 
        scalar = 1
152
165
 
153
166
    if size is None:
154
 
        if scalar:
155
 
            return func(state, oa_data[0])
156
 
        else:
157
 
            array = <ndarray>PyArray_SimpleNew(oa.nd, oa.dimensions, NPY_DOUBLE)
158
 
            length = PyArray_SIZE(array)
159
 
            array_data = <double *>array.data
160
 
            itera = <flatiter>PyArray_IterNew(<object>oa)
161
 
            for i from 0 <= i < length:
162
 
                array_data[i] = func(state, (<double *>(itera.dataptr))[0])
163
 
                PyArray_ITER_NEXT(itera)
 
167
        array = <ndarray>PyArray_SimpleNew(oa.nd, oa.dimensions, NPY_DOUBLE)
 
168
        length = PyArray_SIZE(array)
 
169
        array_data = <double *>array.data
 
170
        itera = <flatiter>PyArray_IterNew(<object>oa)
 
171
        for i from 0 <= i < length:
 
172
            array_data[i] = func(state, (<double *>(itera.dataptr))[0])
 
173
            PyArray_ITER_NEXT(itera)
164
174
    else:
165
175
        array = <ndarray>_sp.empty(size, _sp.float64)
166
176
        array_data = <double *>array.data
167
 
        if scalar:
168
 
            length = PyArray_SIZE(array)
169
 
            for i from 0 <= i < length:
170
 
                array_data[i] = func(state, oa_data[0])
171
 
        else:
172
 
            multi = <broadcast>PyArray_MultiIterNew(2, <void *>array,
173
 
                                                    <void *>oa)
174
 
            if (multi.size != PyArray_SIZE(array)): 
175
 
                raise ValueError("size is not compatible with inputs")
176
 
            for i from 0 <= i < multi.size:
177
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
178
 
                array_data[i] = func(state, oa_data[0])
179
 
                PyArray_MultiIter_NEXTi(multi, 1)
 
177
        multi = <broadcast>PyArray_MultiIterNew(2, <void *>array,
 
178
                                                <void *>oa)
 
179
        if (multi.size != PyArray_SIZE(array)):
 
180
            raise ValueError("size is not compatible with inputs")
 
181
        for i from 0 <= i < multi.size:
 
182
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
183
            array_data[i] = func(state, oa_data[0])
 
184
            PyArray_MultiIter_NEXTi(multi, 1)
180
185
    return array
181
186
 
 
187
cdef object cont2_array_sc(rk_state *state, rk_cont2 func, object size, double a,
 
188
                           double b):
 
189
    cdef double *array_data
 
190
    cdef ndarray array "arrayObject"
 
191
    cdef long length
 
192
    cdef long i
 
193
 
 
194
    if size is None:
 
195
        return func(state, a, b)
 
196
    else:
 
197
        array = <ndarray>_sp.empty(size, _sp.float64)
 
198
        length = PyArray_SIZE(array)
 
199
        array_data = <double *>array.data
 
200
        for i from 0 <= i < length:
 
201
            array_data[i] = func(state, a, b)
 
202
        return array
 
203
 
 
204
 
182
205
cdef object cont2_array(rk_state *state, rk_cont2 func, object size,
183
206
                        ndarray oa, ndarray ob):
184
207
    cdef double *array_data
188
211
    cdef npy_intp length
189
212
    cdef npy_intp i
190
213
    cdef broadcast multi
191
 
    cdef int scalar
192
214
 
193
 
    scalar = 0
194
 
    if oa.nd == 0 and ob.nd == 0:
195
 
        oa_data = <double *>oa.data
196
 
        ob_data = <double *>ob.data
197
 
        scalar = 1
198
 
        
199
215
    if size is None:
200
 
        if scalar:
201
 
            return func(state, oa_data[0], ob_data[0])
202
 
        else:
203
 
            multi = <broadcast> PyArray_MultiIterNew(2, <void *>oa, <void *>ob)
204
 
            array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE)
205
 
            array_data = <double *>array.data
206
 
            for i from 0 <= i < multi.size:
207
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 0)
208
 
                ob_data = <double *>PyArray_MultiIter_DATA(multi, 1)
209
 
                array_data[i] = func(state, oa_data[0], ob_data[0])
210
 
                PyArray_MultiIter_NEXT(multi)
 
216
        multi = <broadcast> PyArray_MultiIterNew(2, <void *>oa, <void *>ob)
 
217
        array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE)
 
218
        array_data = <double *>array.data
 
219
        for i from 0 <= i < multi.size:
 
220
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 0)
 
221
            ob_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
222
            array_data[i] = func(state, oa_data[0], ob_data[0])
 
223
            PyArray_MultiIter_NEXT(multi)
211
224
    else:
212
225
        array = <ndarray>_sp.empty(size, _sp.float64)
213
226
        array_data = <double *>array.data
214
 
        if scalar:
215
 
            length = PyArray_SIZE(array)
216
 
            for i from 0 <= i < length:
217
 
                array_data[i] = func(state, oa_data[0], ob_data[0])
218
 
        else:
219
 
            multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>oa, <void *>ob)
220
 
            if (multi.size != PyArray_SIZE(array)):
221
 
                raise ValueError("size is not compatible with inputs")                 
222
 
            for i from 0 <= i < multi.size:
223
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
224
 
                ob_data = <double *>PyArray_MultiIter_DATA(multi, 2)
225
 
                array_data[i] = func(state, oa_data[0], ob_data[0])
226
 
                PyArray_MultiIter_NEXTi(multi, 1)
227
 
                PyArray_MultiIter_NEXTi(multi, 2)
 
227
        multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>oa, <void *>ob)
 
228
        if (multi.size != PyArray_SIZE(array)):
 
229
            raise ValueError("size is not compatible with inputs")
 
230
        for i from 0 <= i < multi.size:
 
231
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
232
            ob_data = <double *>PyArray_MultiIter_DATA(multi, 2)
 
233
            array_data[i] = func(state, oa_data[0], ob_data[0])
 
234
            PyArray_MultiIter_NEXTi(multi, 1)
 
235
            PyArray_MultiIter_NEXTi(multi, 2)
228
236
    return array
229
237
 
230
 
cdef object cont3_array(rk_state *state, rk_cont3 func, object size, ndarray oa, 
 
238
cdef object cont3_array_sc(rk_state *state, rk_cont3 func, object size, double a,
 
239
                           double b, double c):
 
240
 
 
241
    cdef double *array_data
 
242
    cdef ndarray array "arrayObject"
 
243
    cdef long length
 
244
    cdef long i
 
245
 
 
246
    if size is None:
 
247
        return func(state, a, b, c)
 
248
    else:
 
249
        array = <ndarray>_sp.empty(size, _sp.float64)
 
250
        length = PyArray_SIZE(array)
 
251
        array_data = <double *>array.data
 
252
        for i from 0 <= i < length:
 
253
            array_data[i] = func(state, a, b, c)
 
254
        return array
 
255
 
 
256
cdef object cont3_array(rk_state *state, rk_cont3 func, object size, ndarray oa,
231
257
    ndarray ob, ndarray oc):
232
258
 
233
259
    cdef double *array_data
238
264
    cdef npy_intp length
239
265
    cdef npy_intp i
240
266
    cdef broadcast multi
241
 
    cdef int scalar
242
267
 
243
 
    scalar = 0
244
 
    if (oa.nd ==0 and ob.nd==0 and oc.nd == 0):
245
 
        oa_data = <double *>oa.data
246
 
        ob_data = <double *>ob.data
247
 
        oc_data = <double *>oc.data
248
 
        scalar = 1
249
 
        
250
268
    if size is None:
251
 
        if scalar:
252
 
            return func(state, oa_data[0], ob_data[0], oc_data[0])
253
 
        else:
254
 
            multi = <broadcast> PyArray_MultiIterNew(3, <void *>oa, <void *>ob, <void *>oc)
255
 
            array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE)
256
 
            array_data = <double *>array.data
257
 
            for i from 0 <= i < multi.size:
258
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 0)
259
 
                ob_data = <double *>PyArray_MultiIter_DATA(multi, 1)
260
 
                oc_data = <double *>PyArray_MultiIter_DATA(multi, 2)
261
 
                array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0])
262
 
                PyArray_MultiIter_NEXT(multi)
 
269
        multi = <broadcast> PyArray_MultiIterNew(3, <void *>oa, <void *>ob, <void *>oc)
 
270
        array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE)
 
271
        array_data = <double *>array.data
 
272
        for i from 0 <= i < multi.size:
 
273
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 0)
 
274
            ob_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
275
            oc_data = <double *>PyArray_MultiIter_DATA(multi, 2)
 
276
            array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0])
 
277
            PyArray_MultiIter_NEXT(multi)
263
278
    else:
264
279
        array = <ndarray>_sp.empty(size, _sp.float64)
265
 
        array_data = <double *>array.data        
266
 
        if scalar:
267
 
            length = PyArray_SIZE(array)
268
 
            for i from 0 <= i < length:
269
 
                array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0])
270
 
        else:
271
 
            multi = <broadcast>PyArray_MultiIterNew(4, <void*>array, <void *>oa,
272
 
                                                    <void *>ob, <void *>oc)
273
 
            if (multi.size != PyArray_SIZE(array)):
274
 
                raise ValueError("size is not compatible with inputs")                 
275
 
            for i from 0 <= i < multi.size:
276
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
277
 
                ob_data = <double *>PyArray_MultiIter_DATA(multi, 2)
278
 
                oc_data = <double *>PyArray_MultiIter_DATA(multi, 3)
279
 
                array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0])
280
 
                PyArray_MultiIter_NEXT(multi)                
 
280
        array_data = <double *>array.data
 
281
        multi = <broadcast>PyArray_MultiIterNew(4, <void*>array, <void *>oa,
 
282
                                                <void *>ob, <void *>oc)
 
283
        if (multi.size != PyArray_SIZE(array)):
 
284
            raise ValueError("size is not compatible with inputs")
 
285
        for i from 0 <= i < multi.size:
 
286
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
287
            ob_data = <double *>PyArray_MultiIter_DATA(multi, 2)
 
288
            oc_data = <double *>PyArray_MultiIter_DATA(multi, 3)
 
289
            array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0])
 
290
            PyArray_MultiIter_NEXT(multi)
281
291
    return array
282
292
 
283
293
cdef object disc0_array(rk_state *state, rk_disc0 func, object size):
296
306
            array_data[i] = func(state)
297
307
        return array
298
308
 
 
309
cdef object discnp_array_sc(rk_state *state, rk_discnp func, object size, long n, double p):
 
310
    cdef long *array_data
 
311
    cdef ndarray array "arrayObject"
 
312
    cdef long length
 
313
    cdef long i
 
314
 
 
315
    if size is None:
 
316
        return func(state, n, p)
 
317
    else:
 
318
        array = <ndarray>_sp.empty(size, int)
 
319
        length = PyArray_SIZE(array)
 
320
        array_data = <long *>array.data
 
321
        for i from 0 <= i < length:
 
322
            array_data[i] = func(state, n, p)
 
323
        return array
 
324
 
299
325
cdef object discnp_array(rk_state *state, rk_discnp func, object size, ndarray on, ndarray op):
300
326
    cdef long *array_data
301
327
    cdef ndarray array "arrayObject"
303
329
    cdef npy_intp i
304
330
    cdef double *op_data
305
331
    cdef long *on_data
306
 
    cdef int scalar
307
332
    cdef broadcast multi
308
333
 
309
 
    scalar = 0
310
 
    if (on.nd == 0 and op.nd == 0):
311
 
        on_data = <long *>on.data
312
 
        op_data = <double *>op.data
313
 
        scalar = 1
314
 
        
315
334
    if size is None:
316
 
        if (scalar):
317
 
            return func(state, on_data[0], op_data[0])
318
 
        else:
319
 
            multi = <broadcast> PyArray_MultiIterNew(2, <void *>on, <void *>op)
320
 
            array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG)
321
 
            array_data = <long *>array.data
322
 
            for i from 0 <= i < multi.size:
323
 
                on_data = <long *>PyArray_MultiIter_DATA(multi, 0)
324
 
                op_data = <double *>PyArray_MultiIter_DATA(multi, 1)
325
 
                array_data[i] = func(state, on_data[0], op_data[0])
326
 
                PyArray_MultiIter_NEXT(multi)
 
335
        multi = <broadcast> PyArray_MultiIterNew(2, <void *>on, <void *>op)
 
336
        array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG)
 
337
        array_data = <long *>array.data
 
338
        for i from 0 <= i < multi.size:
 
339
            on_data = <long *>PyArray_MultiIter_DATA(multi, 0)
 
340
            op_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
341
            array_data[i] = func(state, on_data[0], op_data[0])
 
342
            PyArray_MultiIter_NEXT(multi)
327
343
    else:
328
 
        array = <ndarray>_sp.empty(size, int)        
329
 
        if (scalar):
330
 
            length = PyArray_SIZE(array)
331
 
            array_data = <long *>array.data
332
 
            for i from 0 <= i < length:
333
 
                array_data[i] = func(state, on_data[0], op_data[0])
334
 
        else:
335
 
            multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>on, <void *>op)
336
 
            if (multi.size != PyArray_SIZE(array)):
337
 
                raise ValueError("size is not compatible with inputs")                 
338
 
            for i from 0 <= i < multi.size:
339
 
                on_data = <long *>PyArray_MultiIter_DATA(multi, 1)
340
 
                op_data = <double *>PyArray_MultiIter_DATA(multi, 2)
341
 
                array_data[i] = func(state, on_data[0], op_data[0])
342
 
                PyArray_MultiIter_NEXTi(multi, 1)
343
 
                PyArray_MultiIter_NEXTi(multi, 2)
344
 
            
 
344
        array = <ndarray>_sp.empty(size, int)
 
345
        array_data = <long *>array.data
 
346
        multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>on, <void *>op)
 
347
        if (multi.size != PyArray_SIZE(array)):
 
348
            raise ValueError("size is not compatible with inputs")
 
349
        for i from 0 <= i < multi.size:
 
350
            on_data = <long *>PyArray_MultiIter_DATA(multi, 1)
 
351
            op_data = <double *>PyArray_MultiIter_DATA(multi, 2)
 
352
            array_data[i] = func(state, on_data[0], op_data[0])
 
353
            PyArray_MultiIter_NEXTi(multi, 1)
 
354
            PyArray_MultiIter_NEXTi(multi, 2)
 
355
 
345
356
    return array
346
357
 
347
 
cdef object discnmN_array(rk_state *state, rk_discnmN func, object size, 
 
358
cdef object discnmN_array_sc(rk_state *state, rk_discnmN func, object size,
 
359
    long n, long m, long N):
 
360
    cdef long *array_data
 
361
    cdef ndarray array "arrayObject"
 
362
    cdef long length
 
363
    cdef long i
 
364
 
 
365
    if size is None:
 
366
        return func(state, n, m, N)
 
367
    else:
 
368
        array = <ndarray>_sp.empty(size, int)
 
369
        length = PyArray_SIZE(array)
 
370
        array_data = <long *>array.data
 
371
        for i from 0 <= i < length:
 
372
            array_data[i] = func(state, n, m, N)
 
373
        return array
 
374
 
 
375
cdef object discnmN_array(rk_state *state, rk_discnmN func, object size,
348
376
    ndarray on, ndarray om, ndarray oN):
349
377
    cdef long *array_data
350
378
    cdef long *on_data
354
382
    cdef npy_intp length
355
383
    cdef npy_intp i
356
384
    cdef broadcast multi
357
 
    cdef int scalar
358
 
 
359
 
    scalar = 0
360
 
    if (on.nd == 0 and om.nd == 0 and oN.nd==0):
361
 
        scalar = 1
362
 
        on_data = <long *>on.data
363
 
        om_data = <long *>om.data
364
 
        oN_data = <long *>oN.data        
365
385
 
366
386
    if size is None:
367
 
        if (scalar):
368
 
            return func(state, on_data[0], om_data[0], oN_data[0])
369
 
        else:
370
 
            multi = <broadcast> PyArray_MultiIterNew(3, <void *>on, <void *>om, <void *>oN)
371
 
            array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG)
372
 
            array_data = <long *>array.data
373
 
            for i from 0 <= i < multi.size:
374
 
                on_data = <long *>PyArray_MultiIter_DATA(multi, 0)
375
 
                om_data = <long *>PyArray_MultiIter_DATA(multi, 1)
376
 
                oN_data = <long *>PyArray_MultiIter_DATA(multi, 2)                
377
 
                array_data[i] = func(state, on_data[0], om_data[0], oN_data[0])
378
 
                PyArray_MultiIter_NEXT(multi)
 
387
        multi = <broadcast> PyArray_MultiIterNew(3, <void *>on, <void *>om, <void *>oN)
 
388
        array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG)
 
389
        array_data = <long *>array.data
 
390
        for i from 0 <= i < multi.size:
 
391
            on_data = <long *>PyArray_MultiIter_DATA(multi, 0)
 
392
            om_data = <long *>PyArray_MultiIter_DATA(multi, 1)
 
393
            oN_data = <long *>PyArray_MultiIter_DATA(multi, 2)
 
394
            array_data[i] = func(state, on_data[0], om_data[0], oN_data[0])
 
395
            PyArray_MultiIter_NEXT(multi)
379
396
    else:
380
397
        array = <ndarray>_sp.empty(size, int)
381
398
        array_data = <long *>array.data
382
 
        if (scalar):
383
 
            length = PyArray_SIZE(array)
384
 
            for i from 0 <= i < length:
385
 
                array_data[i] = func(state, on_data[0], om_data[0], oN_data[0])
386
 
        else:
387
 
            multi = <broadcast>PyArray_MultiIterNew(4, <void*>array, <void *>on, <void *>om,
388
 
                                                    <void *>oN)
389
 
            if (multi.size != PyArray_SIZE(array)):
390
 
                raise ValueError("size is not compatible with inputs") 
391
 
            for i from 0 <= i < multi.size:
392
 
                on_data = <long *>PyArray_MultiIter_DATA(multi, 1)
393
 
                om_data = <long *>PyArray_MultiIter_DATA(multi, 2)
394
 
                oN_data = <long *>PyArray_MultiIter_DATA(multi, 3)
395
 
                array_data[i] = func(state, on_data[0], om_data[0], oN_data[0])
396
 
                PyArray_MultiIter_NEXT(multi)
 
399
        multi = <broadcast>PyArray_MultiIterNew(4, <void*>array, <void *>on, <void *>om,
 
400
                                                <void *>oN)
 
401
        if (multi.size != PyArray_SIZE(array)):
 
402
            raise ValueError("size is not compatible with inputs")
 
403
        for i from 0 <= i < multi.size:
 
404
            on_data = <long *>PyArray_MultiIter_DATA(multi, 1)
 
405
            om_data = <long *>PyArray_MultiIter_DATA(multi, 2)
 
406
            oN_data = <long *>PyArray_MultiIter_DATA(multi, 3)
 
407
            array_data[i] = func(state, on_data[0], om_data[0], oN_data[0])
 
408
            PyArray_MultiIter_NEXT(multi)
397
409
 
398
410
    return array
399
411
 
 
412
cdef object discd_array_sc(rk_state *state, rk_discd func, object size, double a):
 
413
    cdef long *array_data
 
414
    cdef ndarray array "arrayObject"
 
415
    cdef long length
 
416
    cdef long i
 
417
 
 
418
    if size is None:
 
419
        return func(state, a)
 
420
    else:
 
421
        array = <ndarray>_sp.empty(size, int)
 
422
        length = PyArray_SIZE(array)
 
423
        array_data = <long *>array.data
 
424
        for i from 0 <= i < length:
 
425
            array_data[i] = func(state, a)
 
426
        return array
 
427
 
400
428
cdef object discd_array(rk_state *state, rk_discd func, object size, ndarray oa):
401
429
    cdef long *array_data
402
430
    cdef double *oa_data
405
433
    cdef npy_intp i
406
434
    cdef broadcast multi
407
435
    cdef flatiter itera
408
 
    cdef int scalar 
409
 
 
410
 
    scalar = 0
411
 
    if (oa.nd == 0):
412
 
        oa_data = <double *>oa.data
413
 
        scalar =1
414
436
 
415
437
    if size is None:
416
 
        if (scalar):
417
 
            return func(state, oa_data[0])
418
 
        else:
419
 
            array = <ndarray>PyArray_SimpleNew(oa.nd, oa.dimensions, NPY_LONG)
420
 
            length = PyArray_SIZE(array)
421
 
            array_data = <long *>array.data
422
 
            itera = <flatiter>PyArray_IterNew(<object>oa)
423
 
            for i from 0 <= i < length:
424
 
                array_data[i] = func(state, (<double *>(itera.dataptr))[0])
425
 
                PyArray_ITER_NEXT(itera)
 
438
        array = <ndarray>PyArray_SimpleNew(oa.nd, oa.dimensions, NPY_LONG)
 
439
        length = PyArray_SIZE(array)
 
440
        array_data = <long *>array.data
 
441
        itera = <flatiter>PyArray_IterNew(<object>oa)
 
442
        for i from 0 <= i < length:
 
443
            array_data[i] = func(state, (<double *>(itera.dataptr))[0])
 
444
            PyArray_ITER_NEXT(itera)
426
445
    else:
427
446
        array = <ndarray>_sp.empty(size, int)
428
447
        array_data = <long *>array.data
429
 
        if (scalar):
430
 
            length = PyArray_SIZE(array)
431
 
            for i from 0 <= i < length:
432
 
                array_data[i] = func(state, oa_data[0])
433
 
        else:
434
 
            multi = <broadcast>PyArray_MultiIterNew(2, <void *>array, <void *>oa)
435
 
            if (multi.size != PyArray_SIZE(array)): 
436
 
                raise ValueError("size is not compatible with inputs")                
437
 
            for i from 0 <= i < multi.size:
438
 
                oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
439
 
                array_data[i] = func(state, oa_data[0])
440
 
                PyArray_MultiIter_NEXTi(multi, 1)
 
448
        multi = <broadcast>PyArray_MultiIterNew(2, <void *>array, <void *>oa)
 
449
        if (multi.size != PyArray_SIZE(array)):
 
450
            raise ValueError("size is not compatible with inputs")
 
451
        for i from 0 <= i < multi.size:
 
452
            oa_data = <double *>PyArray_MultiIter_DATA(multi, 1)
 
453
            array_data[i] = func(state, oa_data[0])
 
454
            PyArray_MultiIter_NEXTi(multi, 1)
441
455
    return array
442
456
 
443
457
cdef double kahan_sum(double *darr, long n):
498
512
            errcode = rk_randomseed(self.internal_state)
499
513
        elif type(seed) is int:
500
514
            rk_seed(seed, self.internal_state)
 
515
        elif isinstance(seed, _sp.integer):
 
516
            iseed = int(seed)
 
517
            rk_seed(iseed, self.internal_state)
501
518
        else:
502
519
            obj = <ndarray>PyArray_ContiguousFromObject(seed, NPY_LONG, 1, 1)
503
520
            init_by_array(self.internal_state, <unsigned long *>(obj.data),
506
523
    def get_state(self):
507
524
        """Return a tuple representing the internal state of the generator.
508
525
 
509
 
        get_state() -> ('MT19937', int key[624], int pos)
 
526
        get_state() -> ('MT19937', int key[624], int pos, int has_gauss, float cached_gaussian)
510
527
        """
511
528
        cdef ndarray state "arrayObject_state"
512
 
        state = <ndarray>_sp.empty(624, int)
513
 
        memcpy(<void*>(state.data), self.internal_state.key, 624*sizeof(long))
514
 
        return ('MT19937', state, self.internal_state.pos)
515
 
        
 
529
        state = <ndarray>_sp.empty(624, _sp.uint)
 
530
        memcpy(<void*>(state.data), <void*>(self.internal_state.key), 624*sizeof(long))
 
531
        state = <ndarray>_sp.asarray(state, _sp.uint32)
 
532
        return ('MT19937', state, self.internal_state.pos,
 
533
            self.internal_state.has_gauss, self.internal_state.gauss)
 
534
 
516
535
    def set_state(self, state):
517
536
        """Set the state from a tuple.
518
 
        
 
537
 
 
538
        state = ('MT19937', int key[624], int pos, int has_gauss, float cached_gaussian)
 
539
 
 
540
        For backwards compatibility, the following form is also accepted
 
541
        although it is missing some information about the cached Gaussian value.
 
542
 
519
543
        state = ('MT19937', int key[624], int pos)
520
 
        
 
544
 
521
545
        set_state(state)
522
546
        """
523
547
        cdef ndarray obj "arrayObject_obj"
525
549
        algorithm_name = state[0]
526
550
        if algorithm_name != 'MT19937':
527
551
            raise ValueError("algorithm must be 'MT19937'")
528
 
        key, pos = state[1:]
529
 
        obj = <ndarray>PyArray_ContiguousFromObject(key, NPY_LONG, 1, 1)
 
552
        key, pos = state[1:3]
 
553
        if len(state) == 3:
 
554
            has_gauss = 0
 
555
            cached_gaussian = 0.0
 
556
        else:
 
557
            has_gauss, cached_gaussian = state[3:5]
 
558
        try:
 
559
            obj = <ndarray>PyArray_ContiguousFromObject(key, NPY_ULONG, 1, 1)
 
560
        except TypeError:
 
561
            # compatibility -- could be an older pickle
 
562
            obj = <ndarray>PyArray_ContiguousFromObject(key, NPY_LONG, 1, 1)
530
563
        if obj.dimensions[0] != 624:
531
564
            raise ValueError("state must be 624 longs")
532
 
        memcpy(self.internal_state.key, <void*>(obj.data), 624*sizeof(long))
 
565
        memcpy(<void*>(self.internal_state.key), <void*>(obj.data), 624*sizeof(long))
533
566
        self.internal_state.pos = pos
534
 
    
 
567
        self.internal_state.has_gauss = has_gauss
 
568
        self.internal_state.gauss = cached_gaussian
 
569
 
535
570
    # Pickling support:
536
571
    def __getstate__(self):
537
572
        return self.get_state()
580
615
        diff = hi - lo - 1
581
616
        if diff < 0:
582
617
            raise ValueError("low >= high")
583
 
    
 
618
 
584
619
        if size is None:
585
 
            return rk_interval(diff, self.internal_state) + lo
 
620
            return <long>rk_interval(diff, self.internal_state) + lo
586
621
        else:
587
622
            array = <ndarray>_sp.empty(size, int)
588
623
            length = PyArray_SIZE(array)
597
632
        bytes(length) -> str
598
633
        """
599
634
        cdef void *bytes
600
 
        bytes = PyMem_Malloc(length)
 
635
        bytestring = PyString_FromStringAndSize(NULL, length)
 
636
        bytes = PyString_AS_STRING(bytestring)
601
637
        rk_fill(bytes, length, self.internal_state)
602
 
        bytestring = PyString_FromStringAndSize(<char*>bytes, length)
603
 
        PyMem_Free(bytes)
604
638
        return bytestring
605
639
 
606
640
    def uniform(self, low=0.0, high=1.0, size=None):
608
642
 
609
643
        uniform(low=0.0, high=1.0, size=None) -> random values
610
644
        """
611
 
        cdef ndarray olow
612
 
        cdef ndarray ohigh
613
 
        cdef ndarray odiff
 
645
        cdef ndarray olow, ohigh, odiff
 
646
        cdef double flow, fhigh
614
647
        cdef object temp
 
648
 
 
649
        flow = PyFloat_AsDouble(low)
 
650
        fhigh = PyFloat_AsDouble(high)
 
651
        if not PyErr_Occurred():
 
652
            return cont2_array_sc(self.internal_state, rk_uniform, size, flow, fhigh-flow)
 
653
        PyErr_Clear()
615
654
        olow = <ndarray>PyArray_FROM_OTF(low, NPY_DOUBLE, NPY_ALIGNED)
616
655
        ohigh = <ndarray>PyArray_FROM_OTF(high, NPY_DOUBLE, NPY_ALIGNED)
617
656
        temp = _sp.subtract(ohigh, olow)
621
660
        return cont2_array(self.internal_state, rk_uniform, size, olow, odiff)
622
661
 
623
662
    def rand(self, *args):
624
 
        """Return an array of the given dimensions which is initialized to 
 
663
        """Return an array of the given dimensions which is initialized to
625
664
        random numbers from a uniform distribution in the range [0,1).
626
665
 
627
666
        rand(d0, d1, ..., dn) -> random values
629
668
        Note:  This is a convenience function. If you want an
630
669
                    interface that takes a tuple as the first argument
631
670
                    use numpy.random.random_sample(shape_tuple).
632
 
        
 
671
 
633
672
        """
634
673
        if len(args) == 0:
635
674
            return self.random_sample()
637
676
            return self.random_sample(size=args)
638
677
 
639
678
    def randn(self, *args):
640
 
        """Returns zero-mean, unit-variance Gaussian random numbers in an 
 
679
        """Returns zero-mean, unit-variance Gaussian random numbers in an
641
680
        array of shape (d0, d1, ..., dn).
642
681
 
643
682
        randn(d0, d1, ..., dn) -> random values
676
715
 
677
716
        normal(loc=0.0, scale=1.0, size=None) -> random values
678
717
        """
679
 
        cdef ndarray oloc
680
 
        cdef ndarray oscale
 
718
        cdef ndarray oloc, oscale
 
719
        cdef double floc, fscale
 
720
 
 
721
        floc = PyFloat_AsDouble(loc)
 
722
        fscale = PyFloat_AsDouble(scale)
 
723
        if not PyErr_Occurred():
 
724
            if fscale <= 0:
 
725
                raise ValueError("scale <= 0")
 
726
            return cont2_array_sc(self.internal_state, rk_normal, size, floc, fscale)
 
727
 
 
728
        PyErr_Clear()
 
729
 
681
730
        oloc = <ndarray>PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ALIGNED)
682
731
        oscale = <ndarray>PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
683
732
        if _sp.any(_sp.less_equal(oscale, 0)):
689
738
 
690
739
        beta(a, b, size=None) -> random values
691
740
        """
692
 
        cdef ndarray oa
693
 
        cdef ndarray ob
 
741
        cdef ndarray oa, ob
 
742
        cdef double fa, fb
 
743
 
 
744
        fa = PyFloat_AsDouble(a)
 
745
        fb = PyFloat_AsDouble(b)
 
746
        if not PyErr_Occurred():
 
747
            if fa <= 0:
 
748
                raise ValueError("a <= 0")
 
749
            if fb <= 0:
 
750
                raise ValueError("b <= 0")
 
751
            return cont2_array_sc(self.internal_state, rk_beta, size, fa, fb)
 
752
 
 
753
        PyErr_Clear()
 
754
 
694
755
        oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ALIGNED)
695
756
        ob = <ndarray>PyArray_FROM_OTF(b, NPY_DOUBLE, NPY_ALIGNED)
696
 
 
697
757
        if _sp.any(_sp.less_equal(oa, 0)):
698
758
            raise ValueError("a <= 0")
699
 
        if _sp.any(_sp.less_equal(ob, 0)):        
 
759
        if _sp.any(_sp.less_equal(ob, 0)):
700
760
            raise ValueError("b <= 0")
701
761
        return cont2_array(self.internal_state, rk_beta, size, oa, ob)
702
762
 
706
766
        exponential(scale=1.0, size=None) -> random values
707
767
        """
708
768
        cdef ndarray oscale
 
769
        cdef double fscale
 
770
 
 
771
        fscale = PyFloat_AsDouble(scale)
 
772
        if not PyErr_Occurred():
 
773
            if fscale <= 0:
 
774
                raise ValueError("scale <= 0")
 
775
            return cont1_array_sc(self.internal_state, rk_exponential, size, fscale)
 
776
 
 
777
        PyErr_Clear()
 
778
 
709
779
        oscale = <ndarray> PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
710
780
        if _sp.any(_sp.less_equal(oscale, 0.0)):
711
781
            raise ValueError("scale <= 0")
724
794
        standard_gamma(shape, size=None) -> random values
725
795
        """
726
796
        cdef ndarray oshape
 
797
        cdef double fshape
 
798
 
 
799
        fshape = PyFloat_AsDouble(shape)
 
800
        if not PyErr_Occurred():
 
801
            if fshape <= 0:
 
802
                raise ValueError("shape <= 0")
 
803
            return cont1_array_sc(self.internal_state, rk_standard_gamma, size, fshape)
 
804
 
 
805
        PyErr_Clear()
727
806
        oshape = <ndarray> PyArray_FROM_OTF(shape, NPY_DOUBLE, NPY_ALIGNED)
728
807
        if _sp.any(_sp.less_equal(oshape, 0.0)):
729
808
            raise ValueError("shape <= 0")
734
813
 
735
814
        gamma(shape, scale=1.0, size=None) -> random values
736
815
        """
737
 
        cdef ndarray oshape
738
 
        cdef ndarray oscale
 
816
        cdef ndarray oshape, oscale
 
817
        cdef double fshape, fscale
 
818
 
 
819
        fshape = PyFloat_AsDouble(shape)
 
820
        fscale = PyFloat_AsDouble(scale)
 
821
        if not PyErr_Occurred():
 
822
            if fshape <= 0:
 
823
                raise ValueError("shape <= 0")
 
824
            if fscale <= 0:
 
825
                raise ValueError("scale <= 0")
 
826
            return cont2_array_sc(self.internal_state, rk_gamma, size, fshape, fscale)
 
827
 
 
828
        PyErr_Clear()
739
829
        oshape = <ndarray>PyArray_FROM_OTF(shape, NPY_DOUBLE, NPY_ALIGNED)
740
830
        oscale = <ndarray>PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
741
831
        if _sp.any(_sp.less_equal(oshape, 0.0)):
742
832
            raise ValueError("shape <= 0")
743
 
        if _sp.any(_sp.less_equal(oscale, 0.0)):        
 
833
        if _sp.any(_sp.less_equal(oscale, 0.0)):
744
834
            raise ValueError("scale <= 0")
745
835
        return cont2_array(self.internal_state, rk_gamma, size, oshape, oscale)
746
836
 
749
839
 
750
840
        f(dfnum, dfden, size=None) -> random values
751
841
        """
752
 
        cdef ndarray odfnum
753
 
        cdef ndarray odfden
 
842
        cdef ndarray odfnum, odfden
 
843
        cdef double fdfnum, fdfden
 
844
 
 
845
        fdfnum = PyFloat_AsDouble(dfnum)
 
846
        fdfden = PyFloat_AsDouble(dfden)
 
847
        if not PyErr_Occurred():
 
848
            if fdfnum <= 0:
 
849
                raise ValueError("shape <= 0")
 
850
            if fdfden <= 0:
 
851
                raise ValueError("scale <= 0")
 
852
            return cont2_array_sc(self.internal_state, rk_f, size, fdfnum, fdfden)
 
853
 
 
854
        PyErr_Clear()
 
855
 
754
856
        odfnum = <ndarray>PyArray_FROM_OTF(dfnum, NPY_DOUBLE, NPY_ALIGNED)
755
857
        odfden = <ndarray>PyArray_FROM_OTF(dfden, NPY_DOUBLE, NPY_ALIGNED)
756
858
        if _sp.any(_sp.less_equal(odfnum, 0.0)):
757
 
            raise ValueError("dfnum <= 0")            
758
 
        if _sp.any(_sp.less_equal(odfden, 0.0)):        
 
859
            raise ValueError("dfnum <= 0")
 
860
        if _sp.any(_sp.less_equal(odfden, 0.0)):
759
861
            raise ValueError("dfden <= 0")
760
862
        return cont2_array(self.internal_state, rk_f, size, odfnum, odfden)
761
863
 
764
866
 
765
867
        noncentral_f(dfnum, dfden, nonc, size=None) -> random values
766
868
        """
767
 
        cdef ndarray odfnum
768
 
        cdef ndarray odfden
769
 
        cdef ndarray ononc
 
869
        cdef ndarray odfnum, odfden, ononc
 
870
        cdef double fdfnum, fdfden, fnonc
 
871
 
 
872
        fdfnum = PyFloat_AsDouble(dfnum)
 
873
        fdfden = PyFloat_AsDouble(dfden)
 
874
        fnonc = PyFloat_AsDouble(nonc)
 
875
        if not PyErr_Occurred():
 
876
            if fdfnum <= 1:
 
877
                raise ValueError("dfnum <= 1")
 
878
            if fdfden <= 0:
 
879
                raise ValueError("dfden <= 0")
 
880
            if fnonc < 0:
 
881
                raise ValueError("nonc < 0")
 
882
            return cont3_array_sc(self.internal_state, rk_noncentral_f, size,
 
883
                                  fdfnum, fdfden, fnonc)
 
884
 
 
885
        PyErr_Clear()
 
886
 
770
887
        odfnum = <ndarray>PyArray_FROM_OTF(dfnum, NPY_DOUBLE, NPY_ALIGNED)
771
888
        odfden = <ndarray>PyArray_FROM_OTF(dfden, NPY_DOUBLE, NPY_ALIGNED)
772
 
        ononc = <ndarray>PyArray_FROM_OTF(nonc, NPY_DOUBLE, NPY_ALIGNED)        
773
 
        
 
889
        ononc = <ndarray>PyArray_FROM_OTF(nonc, NPY_DOUBLE, NPY_ALIGNED)
 
890
 
774
891
        if _sp.any(_sp.less_equal(odfnum, 1.0)):
775
 
            raise ValueError("dfnum <= 1")            
776
 
        if _sp.any(_sp.less_equal(odfden, 0.0)):        
 
892
            raise ValueError("dfnum <= 1")
 
893
        if _sp.any(_sp.less_equal(odfden, 0.0)):
777
894
            raise ValueError("dfden <= 0")
778
895
        if _sp.any(_sp.less(ononc, 0.0)):
779
896
            raise ValueError("nonc < 0")
786
903
        chisquare(df, size=None) -> random values
787
904
        """
788
905
        cdef ndarray odf
 
906
        cdef double fdf
 
907
 
 
908
        fdf = PyFloat_AsDouble(df)
 
909
        if not PyErr_Occurred():
 
910
            if fdf <= 0:
 
911
                raise ValueError("df <= 0")
 
912
            return cont1_array_sc(self.internal_state, rk_chisquare, size, fdf)
 
913
 
 
914
        PyErr_Clear()
 
915
 
789
916
        odf = <ndarray>PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ALIGNED)
790
 
        if _sp.any(_sp.less_equal(odf, 0.0)):        
 
917
        if _sp.any(_sp.less_equal(odf, 0.0)):
791
918
            raise ValueError("df <= 0")
792
919
        return cont1_array(self.internal_state, rk_chisquare, size, odf)
793
920
 
796
923
 
797
924
        noncentral_chisquare(df, nonc, size=None) -> random values
798
925
        """
799
 
        cdef ndarray odf
800
 
        cdef ndarray ononc
 
926
        cdef ndarray odf, ononc
 
927
        cdef double fdf, fnonc
 
928
        fdf = PyFloat_AsDouble(df)
 
929
        fnonc = PyFloat_AsDouble(nonc)
 
930
        if not PyErr_Occurred():
 
931
            if fdf <= 1:
 
932
                raise ValueError("df <= 0")
 
933
            if fnonc <= 0:
 
934
                raise ValueError("nonc <= 0")
 
935
            return cont2_array_sc(self.internal_state, rk_noncentral_chisquare,
 
936
                                  size, fdf, fnonc)
 
937
 
 
938
        PyErr_Clear()
 
939
 
801
940
        odf = <ndarray>PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ALIGNED)
802
941
        ononc = <ndarray>PyArray_FROM_OTF(nonc, NPY_DOUBLE, NPY_ALIGNED)
803
942
        if _sp.any(_sp.less_equal(odf, 0.0)):
804
943
            raise ValueError("df <= 1")
805
 
        if _sp.any(_sp.less_equal(ononc, 0.0)):        
 
944
        if _sp.any(_sp.less_equal(ononc, 0.0)):
806
945
            raise ValueError("nonc < 0")
807
946
        return cont2_array(self.internal_state, rk_noncentral_chisquare, size,
808
947
            odf, ononc)
809
 
    
 
948
 
810
949
    def standard_cauchy(self, size=None):
811
950
        """Standard Cauchy with mode=0.
812
951
 
820
959
        standard_t(df, size=None)
821
960
        """
822
961
        cdef ndarray odf
 
962
        cdef double fdf
 
963
 
 
964
        fdf = PyFloat_AsDouble(df)
 
965
        if not PyErr_Occurred():
 
966
            if fdf <= 0:
 
967
                raise ValueError("df <= 0")
 
968
            return cont1_array_sc(self.internal_state, rk_standard_t, size, fdf)
 
969
 
 
970
        PyErr_Clear()
 
971
 
823
972
        odf = <ndarray> PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ALIGNED)
824
 
        if _sp.any(_sp.less_equal(odf, 0.0)):        
 
973
        if _sp.any(_sp.less_equal(odf, 0.0)):
825
974
            raise ValueError("df <= 0")
826
975
        return cont1_array(self.internal_state, rk_standard_t, size, odf)
827
976
 
831
980
 
832
981
        vonmises(mu, kappa, size=None)
833
982
        """
834
 
        cdef ndarray omu
835
 
        cdef ndarray okappa
 
983
        cdef ndarray omu, okappa
 
984
        cdef double fmu, fkappa
 
985
 
 
986
        fmu = PyFloat_AsDouble(mu)
 
987
        fkappa = PyFloat_AsDouble(kappa)
 
988
        if not PyErr_Occurred():
 
989
            if fkappa < 0:
 
990
                raise ValueError("kappa < 0")
 
991
            return cont2_array_sc(self.internal_state, rk_vonmises, size, fmu, fkappa)
 
992
 
 
993
        PyErr_Clear()
 
994
 
836
995
        omu = <ndarray> PyArray_FROM_OTF(mu, NPY_DOUBLE, NPY_ALIGNED)
837
996
        okappa = <ndarray> PyArray_FROM_OTF(kappa, NPY_DOUBLE, NPY_ALIGNED)
838
 
        if _sp.any(_sp.less_equal(okappa, 0.0)):        
 
997
        if _sp.any(_sp.less(okappa, 0.0)):
839
998
            raise ValueError("kappa < 0")
840
999
        return cont2_array(self.internal_state, rk_vonmises, size, omu, okappa)
841
1000
 
845
1004
        pareto(a, size=None)
846
1005
        """
847
1006
        cdef ndarray oa
 
1007
        cdef double fa
 
1008
 
 
1009
        fa = PyFloat_AsDouble(a)
 
1010
        if not PyErr_Occurred():
 
1011
            if fa <= 0:
 
1012
                raise ValueError("a <= 0")
 
1013
            return cont1_array_sc(self.internal_state, rk_pareto, size, fa)
 
1014
 
 
1015
        PyErr_Clear()
 
1016
 
848
1017
        oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ALIGNED)
849
 
        if _sp.any(_sp.less_equal(oa, 0.0)):        
 
1018
        if _sp.any(_sp.less_equal(oa, 0.0)):
850
1019
            raise ValueError("a <= 0")
851
1020
        return cont1_array(self.internal_state, rk_pareto, size, oa)
852
1021
 
853
 
    def weibull(self, double a, size=None):
 
1022
    def weibull(self, a, size=None):
854
1023
        """Weibull distribution.
855
1024
 
856
1025
        weibull(a, size=None)
857
1026
        """
858
1027
        cdef ndarray oa
 
1028
        cdef double fa
 
1029
 
 
1030
        fa = PyFloat_AsDouble(a)
 
1031
        if not PyErr_Occurred():
 
1032
            if fa <= 0:
 
1033
                raise ValueError("a <= 0")
 
1034
            return cont1_array_sc(self.internal_state, rk_weibull, size, fa)
 
1035
 
 
1036
        PyErr_Clear()
 
1037
 
859
1038
        oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ALIGNED)
860
 
        if _sp.any(_sp.less_equal(oa, 0.0)):                
 
1039
        if _sp.any(_sp.less_equal(oa, 0.0)):
861
1040
            raise ValueError("a <= 0")
862
1041
        return cont1_array(self.internal_state, rk_weibull, size, oa)
863
1042
 
864
 
    def power(self, double a, size=None):
 
1043
    def power(self, a, size=None):
865
1044
        """Power distribution.
866
1045
 
867
1046
        power(a, size=None)
868
1047
        """
869
1048
        cdef ndarray oa
 
1049
        cdef double fa
 
1050
 
 
1051
        fa = PyFloat_AsDouble(a)
 
1052
        if not PyErr_Occurred():
 
1053
            if fa <= 0:
 
1054
                raise ValueError("a <= 0")
 
1055
            return cont1_array_sc(self.internal_state, rk_power, size, fa)
 
1056
 
 
1057
        PyErr_Clear()
 
1058
 
870
1059
        oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ALIGNED)
871
 
        if _sp.any(_sp.less_equal(oa, 0.0)):       
 
1060
        if _sp.any(_sp.less_equal(oa, 0.0)):
872
1061
            raise ValueError("a <= 0")
873
1062
        return cont1_array(self.internal_state, rk_power, size, oa)
874
1063
 
875
1064
    def laplace(self, loc=0.0, scale=1.0, size=None):
876
1065
        """Laplace distribution.
877
 
        
 
1066
 
878
1067
        laplace(loc=0.0, scale=1.0, size=None)
879
1068
        """
880
 
        cdef ndarray oloc
881
 
        cdef ndarray oscale
 
1069
        cdef ndarray oloc, oscale
 
1070
        cdef double floc, fscale
 
1071
 
 
1072
        floc = PyFloat_AsDouble(loc)
 
1073
        fscale = PyFloat_AsDouble(scale)
 
1074
        if not PyErr_Occurred():
 
1075
            if fscale <= 0:
 
1076
                raise ValueError("scale <= 0")
 
1077
            return cont2_array_sc(self.internal_state, rk_laplace, size, floc, fscale)
 
1078
 
 
1079
        PyErr_Clear()
882
1080
        oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ALIGNED)
883
1081
        oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
884
 
        if _sp.any(_sp.less_equal(oscale, 0.0)):        
 
1082
        if _sp.any(_sp.less_equal(oscale, 0.0)):
885
1083
            raise ValueError("scale <= 0")
886
1084
        return cont2_array(self.internal_state, rk_laplace, size, oloc, oscale)
887
 
    
 
1085
 
888
1086
    def gumbel(self, loc=0.0, scale=1.0, size=None):
889
1087
        """Gumbel distribution.
890
 
        
 
1088
 
891
1089
        gumbel(loc=0.0, scale=1.0, size=None)
892
1090
        """
893
 
        cdef ndarray oloc
894
 
        cdef ndarray oscale
 
1091
        cdef ndarray oloc, oscale
 
1092
        cdef double floc, fscale
 
1093
 
 
1094
        floc = PyFloat_AsDouble(loc)
 
1095
        fscale = PyFloat_AsDouble(scale)
 
1096
        if not PyErr_Occurred():
 
1097
            if fscale <= 0:
 
1098
                raise ValueError("scale <= 0")
 
1099
            return cont2_array_sc(self.internal_state, rk_gumbel, size, floc, fscale)
 
1100
 
 
1101
        PyErr_Clear()
895
1102
        oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ALIGNED)
896
1103
        oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
897
 
        if _sp.any(_sp.less_equal(oscale, 0.0)):        
 
1104
        if _sp.any(_sp.less_equal(oscale, 0.0)):
898
1105
            raise ValueError("scale <= 0")
899
1106
        return cont2_array(self.internal_state, rk_gumbel, size, oloc, oscale)
900
 
    
 
1107
 
901
1108
    def logistic(self, loc=0.0, scale=1.0, size=None):
902
1109
        """Logistic distribution.
903
 
        
 
1110
 
904
1111
        logistic(loc=0.0, scale=1.0, size=None)
905
1112
        """
906
 
        cdef ndarray oloc
907
 
        cdef ndarray oscale
 
1113
        cdef ndarray oloc, oscale
 
1114
        cdef double floc, fscale
 
1115
 
 
1116
        floc = PyFloat_AsDouble(loc)
 
1117
        fscale = PyFloat_AsDouble(scale)
 
1118
        if not PyErr_Occurred():
 
1119
            if fscale <= 0:
 
1120
                raise ValueError("scale <= 0")
 
1121
            return cont2_array_sc(self.internal_state, rk_logistic, size, floc, fscale)
 
1122
 
 
1123
        PyErr_Clear()
908
1124
        oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ALIGNED)
909
1125
        oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
910
 
        if _sp.any(_sp.less_equal(oscale, 0.0)):        
 
1126
        if _sp.any(_sp.less_equal(oscale, 0.0)):
911
1127
            raise ValueError("scale <= 0")
912
1128
        return cont2_array(self.internal_state, rk_logistic, size, oloc, oscale)
913
1129
 
914
1130
    def lognormal(self, mean=0.0, sigma=1.0, size=None):
915
1131
        """Log-normal distribution.
916
 
        
917
 
        Note that the mean parameter is not the mean of this distribution, but 
 
1132
 
 
1133
        Note that the mean parameter is not the mean of this distribution, but
918
1134
        the underlying normal distribution.
919
 
        
 
1135
 
920
1136
            lognormal(mean, sigma) <=> exp(normal(mean, sigma))
921
 
        
 
1137
 
922
1138
        lognormal(mean=0.0, sigma=1.0, size=None)
923
1139
        """
924
 
        cdef ndarray omean
925
 
        cdef ndarray osigma
 
1140
        cdef ndarray omean, osigma
 
1141
        cdef double fmean, fsigma
 
1142
 
 
1143
        fmean = PyFloat_AsDouble(mean)
 
1144
        fsigma = PyFloat_AsDouble(sigma)
 
1145
 
 
1146
        if not PyErr_Occurred():
 
1147
            if fsigma <= 0:
 
1148
                raise ValueError("sigma <= 0")
 
1149
            return cont2_array_sc(self.internal_state, rk_lognormal, size, fmean, fsigma)
 
1150
 
 
1151
        PyErr_Clear()
 
1152
 
926
1153
        omean = PyArray_FROM_OTF(mean, NPY_DOUBLE, NPY_ALIGNED)
927
1154
        osigma = PyArray_FROM_OTF(sigma, NPY_DOUBLE, NPY_ALIGNED)
928
1155
        if _sp.any(_sp.less_equal(osigma, 0.0)):
929
1156
            raise ValueError("sigma <= 0.0")
930
1157
        return cont2_array(self.internal_state, rk_lognormal, size, omean, osigma)
931
 
    
 
1158
 
932
1159
    def rayleigh(self, scale=1.0, size=None):
933
1160
        """Rayleigh distribution.
934
 
        
 
1161
 
935
1162
        rayleigh(scale=1.0, size=None)
936
1163
        """
937
1164
        cdef ndarray oscale
 
1165
        cdef double fscale
 
1166
 
 
1167
        fscale = PyFloat_AsDouble(scale)
 
1168
 
 
1169
        if not PyErr_Occurred():
 
1170
            if fscale <= 0:
 
1171
                raise ValueError("scale <= 0")
 
1172
            return cont1_array_sc(self.internal_state, rk_rayleigh, size, fscale)
 
1173
 
 
1174
        PyErr_Clear()
 
1175
 
938
1176
        oscale = <ndarray>PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
939
 
        if _sp.any(_sp.less_equal(oscale, 0.0)):        
 
1177
        if _sp.any(_sp.less_equal(oscale, 0.0)):
940
1178
            raise ValueError("scale <= 0.0")
941
1179
        return cont1_array(self.internal_state, rk_rayleigh, size, oscale)
942
 
            
 
1180
 
943
1181
    def wald(self, mean, scale, size=None):
944
1182
        """Wald (inverse Gaussian) distribution.
945
 
        
 
1183
 
946
1184
        wald(mean, scale, size=None)
947
1185
        """
948
 
        cdef ndarray omean
949
 
        cdef ndarray oscale
 
1186
        cdef ndarray omean, oscale
 
1187
        cdef double fmean, fscale
 
1188
 
 
1189
        fmean = PyFloat_AsDouble(mean)
 
1190
        fscale = PyFloat_AsDouble(scale)
 
1191
        if not PyErr_Occurred():
 
1192
            if fmean <= 0:
 
1193
                raise ValueError("mean <= 0")
 
1194
            if fscale <= 0:
 
1195
                raise ValueError("scale <= 0")
 
1196
            return cont2_array_sc(self.internal_state, rk_wald, size, fmean, fscale)
 
1197
 
 
1198
        PyErr_Clear()
950
1199
        omean = PyArray_FROM_OTF(mean, NPY_DOUBLE, NPY_ALIGNED)
951
1200
        oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ALIGNED)
952
1201
        if _sp.any(_sp.less_equal(omean,0.0)):
955
1204
            raise ValueError("scale <= 0.0")
956
1205
        return cont2_array(self.internal_state, rk_wald, size, omean, oscale)
957
1206
 
 
1207
 
 
1208
 
958
1209
    def triangular(self, left, mode, right, size=None):
959
 
        """Triangular distribution starting at left, peaking at mode, and 
 
1210
        """Triangular distribution starting at left, peaking at mode, and
960
1211
        ending at right (left <= mode <= right).
961
 
        
 
1212
 
962
1213
        triangular(left, mode, right, size=None)
963
1214
        """
964
 
        cdef ndarray oleft
965
 
        cdef ndarray omode
966
 
        cdef ndarray oright
 
1215
        cdef ndarray oleft, omode, oright
 
1216
        cdef double fleft, fmode, fright
 
1217
 
 
1218
        fleft = PyFloat_AsDouble(left)
 
1219
        fright = PyFloat_AsDouble(right)
 
1220
        fmode = PyFloat_AsDouble(mode)
 
1221
        if not PyErr_Occurred():
 
1222
            if fleft > fmode:
 
1223
                raise ValueError("left > mode")
 
1224
            if fmode > fright:
 
1225
                raise ValueError("mode > right")
 
1226
            if fleft == fright:
 
1227
                raise ValueError("left == right")
 
1228
            return cont3_array_sc(self.internal_state, rk_triangular, size, fleft,
 
1229
                                  fmode, fright)
 
1230
 
 
1231
        PyErr_Clear()
967
1232
        oleft = <ndarray>PyArray_FROM_OTF(left, NPY_DOUBLE, NPY_ALIGNED)
968
1233
        omode = <ndarray>PyArray_FROM_OTF(mode, NPY_DOUBLE, NPY_ALIGNED)
969
 
        oright = <ndarray>PyArray_FROM_OTF(right, NPY_DOUBLE, NPY_ALIGNED)        
970
 
        
 
1234
        oright = <ndarray>PyArray_FROM_OTF(right, NPY_DOUBLE, NPY_ALIGNED)
 
1235
 
971
1236
        if _sp.any(_sp.greater(oleft, omode)):
972
1237
            raise ValueError("left > mode")
973
1238
        if _sp.any(_sp.greater(omode, oright)):
974
1239
            raise ValueError("mode > right")
975
1240
        if _sp.any(_sp.equal(oleft, oright)):
976
1241
            raise ValueError("left == right")
977
 
        return cont3_array(self.internal_state, rk_triangular, size, oleft, 
 
1242
        return cont3_array(self.internal_state, rk_triangular, size, oleft,
978
1243
            omode, oright)
979
1244
 
980
1245
    # Complicated, discrete distributions:
983
1248
 
984
1249
        binomial(n, p, size=None) -> random values
985
1250
        """
986
 
        cdef ndarray on
987
 
        cdef ndarray op
 
1251
        cdef ndarray on, op
 
1252
        cdef long ln
 
1253
        cdef double fp
 
1254
 
 
1255
        fp = PyFloat_AsDouble(p)
 
1256
        ln = PyInt_AsLong(n)
 
1257
        if not PyErr_Occurred():
 
1258
            if ln <= 0:
 
1259
                raise ValueError("n <= 0")
 
1260
            if fp < 0:
 
1261
                raise ValueError("p < 0")
 
1262
            elif fp > 1:
 
1263
                raise ValueError("p > 1")
 
1264
            return discnp_array_sc(self.internal_state, rk_binomial, size, ln, fp)
 
1265
 
 
1266
        PyErr_Clear()
 
1267
 
988
1268
        on = <ndarray>PyArray_FROM_OTF(n, NPY_LONG, NPY_ALIGNED)
989
 
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)        
 
1269
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)
990
1270
        if _sp.any(_sp.less_equal(n, 0)):
991
1271
            raise ValueError("n <= 0")
992
1272
        if _sp.any(_sp.less(p, 0)):
1002
1282
        """
1003
1283
        cdef ndarray on
1004
1284
        cdef ndarray op
 
1285
        cdef long ln
 
1286
        cdef double fp
 
1287
 
 
1288
        fp = PyFloat_AsDouble(p)
 
1289
        ln = PyInt_AsLong(n)
 
1290
        if not PyErr_Occurred():
 
1291
            if ln <= 0:
 
1292
                raise ValueError("n <= 0")
 
1293
            if fp < 0:
 
1294
                raise ValueError("p < 0")
 
1295
            elif fp > 1:
 
1296
                raise ValueError("p > 1")
 
1297
            return discnp_array_sc(self.internal_state, rk_negative_binomial,
 
1298
                                   size, ln, fp)
 
1299
 
 
1300
        PyErr_Clear()
 
1301
 
1005
1302
        on = <ndarray>PyArray_FROM_OTF(n, NPY_LONG, NPY_ALIGNED)
1006
 
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)        
 
1303
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)
1007
1304
        if _sp.any(_sp.less_equal(n, 0)):
1008
1305
            raise ValueError("n <= 0")
1009
1306
        if _sp.any(_sp.less(p, 0)):
1019
1316
        poisson(lam=1.0, size=None) -> random values
1020
1317
        """
1021
1318
        cdef ndarray olam
 
1319
        cdef double flam
 
1320
        flam = PyFloat_AsDouble(lam)
 
1321
        if not PyErr_Occurred():
 
1322
            if lam < 0:
 
1323
                raise ValueError("lam < 0")
 
1324
            return discd_array_sc(self.internal_state, rk_poisson, size, flam)
 
1325
 
 
1326
        PyErr_Clear()
 
1327
 
1022
1328
        olam = <ndarray>PyArray_FROM_OTF(lam, NPY_DOUBLE, NPY_ALIGNED)
1023
1329
        if _sp.any(_sp.less(olam, 0)):
1024
1330
            raise ValueError("lam < 0")
1026
1332
 
1027
1333
    def zipf(self, a, size=None):
1028
1334
        """Zipf distribution.
1029
 
        
 
1335
 
1030
1336
        zipf(a, size=None)
1031
1337
        """
1032
1338
        cdef ndarray oa
 
1339
        cdef double fa
 
1340
 
 
1341
        fa = PyFloat_AsDouble(a)
 
1342
        if not PyErr_Occurred():
 
1343
            if fa <= 1.0:
 
1344
                raise ValueError("a <= 1.0")
 
1345
            return discd_array_sc(self.internal_state, rk_zipf, size, fa)
 
1346
 
 
1347
        PyErr_Clear()
 
1348
 
1033
1349
        oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ALIGNED)
1034
1350
        if _sp.any(_sp.less_equal(oa, 1.0)):
1035
1351
            raise ValueError("a <= 1.0")
1036
1352
        return discd_array(self.internal_state, rk_zipf, size, oa)
1037
 
    
 
1353
 
1038
1354
    def geometric(self, p, size=None):
1039
1355
        """Geometric distribution with p being the probability of "success" on
1040
1356
        an individual trial.
1041
 
        
 
1357
 
1042
1358
        geometric(p, size=None)
1043
1359
        """
1044
1360
        cdef ndarray op
 
1361
        cdef double fp
 
1362
 
 
1363
        fp = PyFloat_AsDouble(p)
 
1364
        if not PyErr_Occurred():
 
1365
            if fp < 0.0:
 
1366
                raise ValueError("p < 0.0")
 
1367
            if fp > 1.0:
 
1368
                raise ValueError("p > 1.0")
 
1369
            return discd_array_sc(self.internal_state, rk_geometric, size, fp)
 
1370
 
 
1371
        PyErr_Clear()
 
1372
 
 
1373
 
1045
1374
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)
1046
 
        if _sp.any(_sp.less(op, 0.0)):        
 
1375
        if _sp.any(_sp.less(op, 0.0)):
1047
1376
            raise ValueError("p < 0.0")
1048
1377
        if _sp.any(_sp.greater(op, 1.0)):
1049
1378
            raise ValueError("p > 1.0")
1050
1379
        return discd_array(self.internal_state, rk_geometric, size, op)
1051
 
    
 
1380
 
1052
1381
    def hypergeometric(self, ngood, nbad, nsample, size=None):
1053
1382
        """Hypergeometric distribution.
1054
 
        
1055
 
        Consider an urn with ngood "good" balls and nbad "bad" balls. If one 
1056
 
        were to draw nsample balls from the urn without replacement, then 
1057
 
        the hypergeometric distribution describes the distribution of "good" 
 
1383
 
 
1384
        Consider an urn with ngood "good" balls and nbad "bad" balls. If one
 
1385
        were to draw nsample balls from the urn without replacement, then
 
1386
        the hypergeometric distribution describes the distribution of "good"
1058
1387
        balls in the sample.
1059
 
        
1060
 
        hypergeometric(ngood, nbad, nsample, size=None)        
 
1388
 
 
1389
        hypergeometric(ngood, nbad, nsample, size=None)
1061
1390
        """
1062
 
        cdef ndarray ongood
1063
 
        cdef ndarray onbad
1064
 
        cdef ndarray onsample
 
1391
        cdef ndarray ongood, onbad, onsample
 
1392
        cdef long lngood, lnbad, lnsample
 
1393
 
 
1394
        lngood = PyInt_AsLong(ngood)
 
1395
        lnbad = PyInt_AsLong(nbad)
 
1396
        lnsample = PyInt_AsLong(nsample)
 
1397
        if not PyErr_Occurred():
 
1398
            if ngood < 1:
 
1399
                raise ValueError("ngood < 1")
 
1400
            if nbad < 1:
 
1401
                raise ValueError("nbad < 1")
 
1402
            if nsample < 1:
 
1403
                raise ValueError("nsample < 1")
 
1404
            if ngood + nbad < nsample:
 
1405
                raise ValueError("ngood + nbad < nsample")
 
1406
            return discnmN_array_sc(self.internal_state, rk_hypergeometric, size,
 
1407
                                    lngood, lnbad, lnsample)
 
1408
 
 
1409
 
 
1410
        PyErr_Clear()
1065
1411
 
1066
1412
        ongood = <ndarray>PyArray_FROM_OTF(ngood, NPY_LONG, NPY_ALIGNED)
1067
1413
        onbad = <ndarray>PyArray_FROM_OTF(nbad, NPY_LONG, NPY_ALIGNED)
1070
1416
            raise ValueError("ngood < 1")
1071
1417
        if _sp.any(_sp.less(onbad, 1)):
1072
1418
            raise ValueError("nbad < 1")
 
1419
        if _sp.any(_sp.less(onsample, 1)):
 
1420
            raise ValueError("nsample < 1")
1073
1421
        if _sp.any(_sp.less(_sp.add(ongood, onbad),onsample)):
1074
1422
            raise ValueError("ngood + nbad < nsample")
1075
 
        if _sp.any(_sp.less(onsample, 1)):
1076
 
            raise ValueError("nsample < 1")
1077
1423
        return discnmN_array(self.internal_state, rk_hypergeometric, size,
1078
1424
            ongood, onbad, onsample)
1079
1425
 
1080
1426
    def logseries(self, p, size=None):
1081
1427
        """Logarithmic series distribution.
1082
 
        
 
1428
 
1083
1429
        logseries(p, size=None)
1084
1430
        """
1085
1431
        cdef ndarray op
 
1432
        cdef double fp
 
1433
 
 
1434
        fp = PyFloat_AsDouble(p)
 
1435
        if not PyErr_Occurred():
 
1436
            if fp < 0.0:
 
1437
                raise ValueError("p < 0.0")
 
1438
            if fp > 1.0:
 
1439
                raise ValueError("p > 1.0")
 
1440
            return discd_array_sc(self.internal_state, rk_logseries, size, fp)
 
1441
 
 
1442
        PyErr_Clear()
 
1443
 
1086
1444
        op = <ndarray>PyArray_FROM_OTF(p, NPY_DOUBLE, NPY_ALIGNED)
1087
 
        if _sp.any(_sp.less(op, 0.0)):        
 
1445
        if _sp.any(_sp.less(op, 0.0)):
1088
1446
            raise ValueError("p < 0.0")
1089
1447
        if _sp.any(_sp.greater(op, 1.0)):
1090
1448
            raise ValueError("p > 1.0")
1129
1487
        # Create a matrix of independent standard normally distributed random
1130
1488
        # numbers. The matrix has rows with the same length as mean and as
1131
1489
        # many rows are necessary to form a matrix of shape final_shape.
1132
 
        x = standard_normal(_sp.multiply.reduce(final_shape))
 
1490
        x = self.standard_normal(_sp.multiply.reduce(final_shape))
1133
1491
        x.shape = (_sp.multiply.reduce(final_shape[0:len(final_shape)-1]),
1134
1492
                   mean.shape[0])
1135
1493
        # Transform matrix of standard normals into matrix where each row
1138
1496
        # Then the matrix products of the rows of x and A has the desired
1139
1497
        # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value
1140
1498
        # decomposition of cov is such an A.
1141
 
        
 
1499
 
1142
1500
        from numpy.dual import svd
1143
1501
        # XXX: we really should be doing this by Cholesky decomposition
1144
1502
        (u,s,v) = svd(cov)
1151
1509
 
1152
1510
    def multinomial(self, long n, object pvals, size=None):
1153
1511
        """Multinomial distribution.
1154
 
        
 
1512
 
1155
1513
        multinomial(n, pvals, size=None) -> random values
1156
1514
 
1157
1515
        pvals is a sequence of probabilities that should sum to 1 (however, the
1169
1527
        parr = <ndarray>PyArray_ContiguousFromObject(pvals, NPY_DOUBLE, 1, 1)
1170
1528
        pix = <double*>parr.data
1171
1529
 
1172
 
        if kahan_sum(pix, d-1) > 1.0:
1173
 
            raise ValueError("sum(pvals) > 1.0")
 
1530
        if kahan_sum(pix, d-1) > (1.0 + 1e-12):
 
1531
            raise ValueError("sum(pvals[:-1]) > 1.0")
1174
1532
 
1175
1533
        if size is None:
1176
1534
            shape = (d,)
1199
1557
 
1200
1558
        return multin
1201
1559
 
 
1560
    def dirichlet(self, object alpha, size=None):
 
1561
        """dirichlet(alpha, size=None)
 
1562
 
 
1563
        Draw `size` samples of dimension k from a Dirichlet distribution. A
 
1564
        Dirichlet-distributed random variable can be seen as a multivariate
 
1565
        generalization of a Beta distribution. Dirichlet pdf is the conjugate
 
1566
        prior of a multinomial in Bayesian inference.
 
1567
 
 
1568
        Parameters
 
1569
        ----------
 
1570
        alpha : array
 
1571
            Parameter of the distribution (k dimension for sample of
 
1572
            dimension k).
 
1573
        size : array
 
1574
            Number of samples to draw.
 
1575
 
 
1576
        Notes
 
1577
        -----
 
1578
        .. math:: X \\approx \\prod_{i=1}^{k}{x^{\\alpha_i-1}_i}
 
1579
 
 
1580
        Uses the following property for computation: for each dimension,
 
1581
        draw a random sample y_i from a standard gamma generator of shape
 
1582
        `alpha_i`, then
 
1583
        :math:`X = \\frac{1}{\\sum_{i=1}^k{y_i}} (y_1, \\ldot, y_n)` is
 
1584
        Dirichlet distributed.
 
1585
 
 
1586
        References
 
1587
        ----------
 
1588
        .. [1] David McKay, "Information Theory, Inference and Learning
 
1589
               Algorithms," chapter 23,
 
1590
               http://www.inference.phy.cam.ac.uk/mackay/
 
1591
 
 
1592
        """
 
1593
 
 
1594
        #=================
 
1595
        # Pure python algo
 
1596
        #=================
 
1597
        #alpha   = N.atleast_1d(alpha)
 
1598
        #k       = alpha.size
 
1599
 
 
1600
        #if n == 1:
 
1601
        #    val = N.zeros(k)
 
1602
        #    for i in range(k):
 
1603
        #        val[i]   = sgamma(alpha[i], n)
 
1604
        #    val /= N.sum(val)
 
1605
        #else:
 
1606
        #    val = N.zeros((k, n))
 
1607
        #    for i in range(k):
 
1608
        #        val[i]   = sgamma(alpha[i], n)
 
1609
        #    val /= N.sum(val, axis = 0)
 
1610
        #    val = val.T
 
1611
 
 
1612
        #return val
 
1613
 
 
1614
        cdef long       k
 
1615
        cdef long       totsize
 
1616
        cdef ndarray    alpha_arr, val_arr
 
1617
        cdef double     *alpha_data, *val_data
 
1618
        cdef long       i, j
 
1619
        cdef double     acc, invacc
 
1620
 
 
1621
        k           = len(alpha)
 
1622
        alpha_arr   = <ndarray>PyArray_ContiguousFromObject(alpha, NPY_DOUBLE, 1, 1)
 
1623
        alpha_data  = <double*>alpha_arr.data
 
1624
 
 
1625
        if size is None:
 
1626
            shape = (k,)
 
1627
        elif type(size) is int:
 
1628
            shape = (size, k)
 
1629
        else:
 
1630
            shape = size + (k,)
 
1631
 
 
1632
        diric   = _sp.zeros(shape, _sp.float64)
 
1633
        val_arr = <ndarray>diric
 
1634
        val_data= <double*>val_arr.data
 
1635
 
 
1636
        i = 0
 
1637
        totsize = PyArray_SIZE(val_arr)
 
1638
        while i < totsize:
 
1639
            acc = 0.0
 
1640
            for j from 0 <= j < k:
 
1641
                val_data[i+j]   = rk_standard_gamma(self.internal_state, alpha_data[j])
 
1642
                acc             = acc + val_data[i+j]
 
1643
            invacc  = 1/acc
 
1644
            for j from 0 <= j < k:
 
1645
                val_data[i+j]   = val_data[i+j] * invacc
 
1646
            i = i + k
 
1647
 
 
1648
        return diric
 
1649
 
1202
1650
    # Shuffling and permutations:
1203
1651
    def shuffle(self, object x):
1204
1652
        """Modify a sequence in-place by shuffling its contents.
1205
 
        
 
1653
 
1206
1654
        shuffle(x)
1207
1655
        """
1208
1656
        cdef long i, j
 
1657
        cdef int copy
1209
1658
 
1210
 
        # adaptation of random.shuffle()
1211
1659
        i = len(x) - 1
1212
 
        while i > 0:
1213
 
            j = rk_interval(i, self.internal_state)
1214
 
            x[i], x[j] = x[j], x[i]
1215
 
            i = i - 1
1216
 
                
 
1660
        try:
 
1661
            j = len(x[0])
 
1662
        except:
 
1663
            j = 0
 
1664
 
 
1665
        if (j == 0):
 
1666
            # adaptation of random.shuffle()
 
1667
            while i > 0:
 
1668
                j = rk_interval(i, self.internal_state)
 
1669
                x[i], x[j] = x[j], x[i]
 
1670
                i = i - 1
 
1671
        else:
 
1672
            # make copies
 
1673
            copy = hasattr(x[0], 'copy')
 
1674
            if copy:
 
1675
                while(i > 0):
 
1676
                    j = rk_interval(i, self.internal_state)
 
1677
                    x[i], x[j] = x[j].copy(), x[i].copy()
 
1678
                    i = i - 1
 
1679
            else:
 
1680
                while(i > 0):
 
1681
                    j = rk_interval(i, self.internal_state)
 
1682
                    x[i], x[j] = x[j][:], x[i][:]
 
1683
                    i = i - 1
 
1684
 
1217
1685
    def permutation(self, object x):
1218
 
        """Given an integer, return a shuffled sequence of integers >= 0 and 
 
1686
        """Given an integer, return a shuffled sequence of integers >= 0 and
1219
1687
        < x; given a sequence, return a shuffled array copy.
1220
1688
 
1221
1689
        permutation(x)
1222
1690
        """
1223
 
        if type(x) is int:
 
1691
        if isinstance(x, (int, _sp.integer)):
1224
1692
            arr = _sp.arange(x)
1225
1693
        else:
1226
1694
            arr = _sp.array(x)
1227
1695
        self.shuffle(arr)
1228
 
        return arr    
 
1696
        return arr
1229
1697
 
1230
1698
_rand = RandomState()
1231
1699
seed = _rand.seed
1273
1741
 
1274
1742
multivariate_normal = _rand.multivariate_normal
1275
1743
multinomial = _rand.multinomial
 
1744
dirichlet = _rand.dirichlet
1276
1745
 
1277
1746
shuffle = _rand.shuffle
1278
1747
permutation = _rand.permutation