1
/* Copyright (c) 1996, 1997, The Regents of the University of California.
2
* All rights reserved. See Legal.htm for full text and disclaimer. */
4
#include "numpy/arrayobject.h"
5
/*#include "hlevel.h"*/
9
#define MAX_INTERP_DIMS 6
11
static PyObject *ErrorObject;
13
/* Define 2 macros for error handling:
15
If BOOLEAN is FALSE, assume the error object has
16
been set and return NULL
18
Py_Assert(BOOLEAN,ERROBJ,MESS)
19
If BOOLEAN is FALSE set the error object to
20
ERROBJ, and the message to MESS
24
static char * errstr = NULL ;
26
#define Py_Try(BOOLEAN) {if (!(BOOLEAN)) return NULL;}
27
#define Py_Assert(BOOLEAN,ERROBJ,MESS) {if (!(BOOLEAN)) { \
28
PyErr_SetString((ERROBJ), (MESS)); \
31
#define A_DATA(a) (((PyArrayObject *)a)->data)
32
#define A_SIZE(a) PyArray_Size((PyObject *) a)
33
#define A_TYPE(a) (int)(((PyArrayObject *)a)->descr->type_num)
34
#define isARRAY(a) ((a) && PyArray_Check((PyArrayObject *)a))
35
#define A_NDIM(a) (((PyArrayObject *)a)->nd)
36
#define A_DIM(a,i) (((PyArrayObject *)a)->dimensions[i])
37
#define GET_ARR(ap,op,type,dim) \
38
Py_Try(ap=(PyArrayObject *)PyArray_ContiguousFromObject(op,type,dim,dim))
39
#define GET_ARR2(ap,op,type,min,max) \
40
Py_Try(ap=(PyArrayObject *)PyArray_ContiguousFromObject(op,type,min,max))
41
#define ERRSS(s) ((PyObject *)(PyErr_SetString(ErrorObject,s),(void *)0))
42
#define SETERR(s) if(!PyErr_Occurred()) ERRSS(errstr ? errstr : s)
43
#define DECREF_AND_ZERO(p) do{Py_XDECREF(p);p=0;}while(0)
46
/* ----------------------------------------------------- */
48
static char arr_histogram__doc__[] =
52
static int mxx ( int * i , int len)
54
/* find the index of the maximum element of an integer array */
55
int mx = 0, max = i [0] ;
57
for ( j = 1 ; j < len; j ++ )
64
static int mnx ( int * i , int len)
66
/* find the index of the minimum element of an integer array */
67
int mn = 0, min = i [0] ;
69
for ( j = 1 ; j < len; j ++ )
77
arr_histogram(PyObject *self, PyObject *args)
79
/* histogram accepts one or two arguments. The first is an array
80
* of non-negative integers and the second, if present, is an
81
* array of weights, which must be promotable to double.
82
* Call these arguments list and weight. Both must be one-
83
* dimensional. len (weight) >= max (list) + 1.
84
* If weight is not present:
85
* histogram (list) [i] is the number of occurrences of i in list.
86
* If weight is present:
87
* histogram (list, weight) [i] is the sum of all weight [j]
88
* where list [j] == i. */
89
/* self is not used */
90
PyObject * list = NULL, * weight = NULL ;
91
PyArrayObject *lst, *wts , *ans;
92
int * numbers, *ians, len , mxi, mni, i, ans_size;
93
double * weights , * dans ;
95
Py_Try(PyArg_ParseTuple(args, "O|O", &list, &weight));
96
GET_ARR(lst,list,PyArray_INT,1);
98
numbers = (int *) A_DATA (lst) ;
99
mxi = mxx (numbers, len) ;
100
mni = mnx (numbers, len) ;
101
if (numbers [mni] < 0)
102
{SETERR ("First argument of histogram must be nonnegative.");
105
ans_size = numbers [mxi] + 1 ;
109
(PyArrayObject *) PyArray_FromDims (1, &ans_size, PyArray_INT));
110
ians = (int *) A_DATA (ans) ;
111
for (i = 0 ; i < len ; i++)
112
ians [numbers [i]] += 1 ;
117
GET_ARR(wts,weight,PyArray_DOUBLE, 1);
118
weights = (double *) A_DATA (wts) ;
119
if (A_SIZE (wts) != len)
120
{SETERR ("histogram: length of weights does not match that of list.");
125
(PyArrayObject *) PyArray_FromDims (1, &ans_size, PyArray_DOUBLE));
126
dans = (double *) A_DATA (ans);
127
for (i = 0 ; i < len ; i++) {
128
dans [numbers [i]] += weights [i];
134
return PyArray_Return (ans);
137
static char arr_array_set__doc__[] =
142
arr_array_set(PyObject *self, PyObject *args)
144
/* array_set accepts three arguments. The first is an array of
145
* numerics (Python characters, integers, or floats), and the
146
* third is of the same type. The second is an array of integers
147
* which are valid subscripts into the first. The third array
148
* must be at least long enough to supply all the elements
149
* called for by the subscript array. (It can also be a scalar,
150
* in which case its value will be broadcast.) The result is that
151
* elements of the third array are assigned in order to elements
152
* of the first whose subscripts are elements of the second.
153
* arr_array_set (vals1, indices, vals2)
154
* is equivalent to the Yorick assignment vals1 (indices) = vals2.
155
* I have generalized this so that the source and target arrays
156
* may be two dimensional; the second dimensions must match.
157
* Then the array of subscripts is assumed to apply to the first
158
* subscript only of the target. The target had better be contiguous. */
159
/* self is not used */
160
PyObject * tararg, * subsarg, *srcarg;
161
PyArrayObject * tararr, * subsarr, * srcarr = NULL;
162
double * dtar, * dsrc, ds=0.0;
163
float * ftar, * fsrc, fs=0.0;
164
char * ctar, * csrc, cs='\0';
165
unsigned char * utar, * usrc, us=0;
166
int * itar, * isrc, * isubs;
169
int i, j, len, mn, mx;
170
int scalar_source = 0;
171
char scalar_type = 'x';
172
int nd, d1; /* number of dimensions and value of second dim. */
174
Py_Try(PyArg_ParseTuple(args, "OOO", &tararg, &subsarg, &srcarg));
176
nd = A_NDIM (tararg) ;
177
if (PyFloat_Check (srcarg)) {
180
ds = PyFloat_AS_DOUBLE ( (PyFloatObject *) srcarg) ;
182
else if (PyInt_Check (srcarg)) {
185
is = PyInt_AS_LONG ( (PyIntObject *) srcarg) ;
187
else if (PyString_Check (srcarg)) {
190
cs = PyString_AS_STRING ( (PyStringObject *) srcarg) [0] ;
193
d1 = A_DIM (tararg, 1) ;
194
if (A_NDIM (srcarg) != 2 || A_DIM (srcarg,1) != d1) {
195
SETERR ("array_set: dimension mismatch between source and target.");
200
SETERR ("array_set: target must have one or two dimensions.");
203
GET_ARR(subsarr,subsarg,PyArray_INT,1);
204
isubs = (int *)A_DATA(subsarr);
205
len = A_SIZE(subsarr);
206
mn = mnx (isubs, len);
208
{SETERR ("array_set: negative subscript specified.");
211
mx = mxx (isubs, len);
212
switch (A_TYPE(tararg)) {
214
GET_ARR(tararr,tararg,PyArray_UBYTE,nd);
215
if (d1 * isubs [mx] > A_SIZE(tararr))
216
{SETERR ("array_set: a subscript is out of range.");
220
if (! scalar_source) {
221
GET_ARR(srcarr,srcarg,PyArray_UBYTE,1);
222
if (A_SIZE(srcarr) < d1 * len)
224
("array_set: source is too short for number of subscripts.");
229
utar = (unsigned char *)A_DATA(tararr);
230
usrc = (unsigned char *)A_DATA(srcarr);
231
for (i = 0; i < len; i++ )
232
for (j = 0; j < d1; j++)
233
utar [d1 * isubs [i] + j] = usrc [i * d1 + j];
236
switch (scalar_type) {
238
us = (unsigned char) cs ;
241
us = (unsigned char) is ;
244
us = (unsigned char) ds ;
247
utar = (unsigned char *)A_DATA(tararr);
248
for (i = 0; i < len; i++ )
249
for (j = 0; j < d1; j++)
250
utar [d1 * isubs [i] + j] = us;
254
GET_ARR(tararr,tararg,PyArray_CHAR,nd);
255
if (d1 * isubs [mx] > A_SIZE(tararr))
256
{SETERR ("array_set: a subscript is out of range.");
260
if (! scalar_source) {
261
GET_ARR(srcarr,srcarg,PyArray_CHAR,nd);
262
if (A_SIZE(srcarr) < d1 * len)
264
("array_set: source is too short for number of subscripts.");
269
ctar = (char *)A_DATA(tararr);
270
csrc = (char *)A_DATA(srcarr);
271
for (i = 0; i < len; i++ )
272
for (j = 0; j < d1; j++)
273
ctar [isubs [i] * d1 + j] = csrc [i * d1 + j];
276
switch (scalar_type) {
280
cs = (unsigned char) is ;
283
cs = (unsigned char) ds ;
286
ctar = (char *)A_DATA(tararr);
287
for (i = 0; i < len; i++ )
288
for (j = 0; j < d1; j++)
289
ctar [d1 * isubs [i] + j] = cs;
293
GET_ARR(tararr,tararg,PyArray_INT,nd);
294
if (isubs [mx] * d1 > A_SIZE(tararr))
295
{SETERR ("array_set: a subscript is out of range.");
299
if (! scalar_source) {
300
GET_ARR(srcarr,srcarg,PyArray_INT,nd);
301
if (A_SIZE(srcarr) < len * d1)
303
("array_set: source is too short for number of subscripts.");
308
itar = (int *)A_DATA(tararr);
309
isrc = (int *)A_DATA(srcarr);
310
for (i = 0; i < len; i++ )
311
for (j = 0; j < d1; j++)
312
itar [isubs [i] * d1 + j] = isrc [i * d1 + j];
315
switch (scalar_type) {
325
itar = (int *)A_DATA(tararr);
326
for (i = 0; i < len; i++ )
327
for (j = 0; j < d1; j++)
328
itar [d1 * isubs [i] + j] = is;
332
GET_ARR(tararr,tararg,PyArray_LONG,nd);
333
if (isubs [mx] * d1 > A_SIZE(tararr))
334
{SETERR ("array_set: a subscript is out of range.");
338
if (! scalar_source) {
339
GET_ARR(srcarr,srcarg,PyArray_LONG,nd);
340
if (A_SIZE(srcarr) < len * d1)
342
("array_set: source is too short for number of subscripts.");
347
ltar = (long *)A_DATA(tararr);
348
lsrc = (long *)A_DATA(srcarr);
349
for (i = 0; i < len; i++ )
350
for (j = 0; j < d1; j++)
351
ltar [isubs [i] * d1 + j] = lsrc [i * d1 + j];
354
switch (scalar_type) {
364
ltar = (long *)A_DATA(tararr);
365
for (i = 0; i < len; i++ )
366
for (j = 0; j < d1; j++)
367
ltar [d1 * isubs [i] + j] = is;
371
GET_ARR(tararr,tararg,PyArray_FLOAT,nd);
372
if (isubs [mx] * d1 > A_SIZE(tararr))
373
{SETERR ("array_set: a subscript is out of range.");
377
if (! scalar_source) {
378
GET_ARR(srcarr,srcarg,PyArray_FLOAT,nd);
379
if (A_SIZE(srcarr) < len * d1)
381
("array_set: source is too short for number of subscripts.");
386
ftar = (float *)A_DATA(tararr);
387
fsrc = (float *)A_DATA(srcarr);
388
for (i = 0; i < len; i++ )
389
for (j = 0; j < d1; j++)
390
ftar [isubs [i] * d1 + j] = fsrc [i * d1 + j];
393
switch (scalar_type) {
404
ftar = (float *)A_DATA(tararr);
405
for (i = 0; i < len; i++ )
406
for (j = 0; j < d1; j++)
407
ftar [d1 * isubs [i] + j] = fs;
411
GET_ARR(tararr,tararg,PyArray_DOUBLE,nd);
412
if (isubs [mx] * d1 > A_SIZE(tararr))
413
{SETERR ("array_set: a subscript is out of range.");
417
if (! scalar_source) {
418
GET_ARR(srcarr,srcarg,PyArray_DOUBLE,nd);
419
if (A_SIZE(srcarr) < len * d1)
421
("array_set: source is too short for number of subscripts.");
426
dtar = (double *)A_DATA(tararr);
427
dsrc = (double *)A_DATA(srcarr);
428
for (i = 0; i < len; i++ )
429
for (j = 0; j < d1; j++)
430
dtar [isubs [i] * d1 + j] = dsrc [i * d1 + j];
433
switch (scalar_type) {
443
dtar = (double *)A_DATA(tararr);
444
for (i = 0; i < len; i++ )
445
for (j = 0; j < d1; j++)
446
dtar [d1 * isubs [i] + j] = ds;
450
SETERR("array_set: Not implemented for this type.");
462
static void adjust (double * k, int * list, int i, int n)
464
/* adjust the binary tree k with root list [i] to satisfy the heap
465
* property. The left and right subtrees of list [i], with roots
466
* list [2 * i + 1] and list [2 * i + 2], already satisfy the heap
467
* property. No node has index greater than n.
468
* Horowitz & Sahni, p. 358. */
469
double kt; /* will contain root value which may need to be moved */
470
int kj , /* root value is at kj position in list */
472
lowj ; /* parent of current j node */
480
if (j < n - 1 && k [list [j]] < k [list [j + 1]])
481
/* make list [j] point to the larger child */
483
if ( kt >= k [list [j]] )
488
list [lowj] = list [j] ;
495
static char arr_index_sort__doc__[] =
500
arr_index_sort(PyObject *self, PyObject *args)
502
/* index_sort accepts one array of some numerical type and returns
503
* an integer array of the same length whose entries are the
504
* subscripts of the elements of the original array arranged
505
* in increasing order. I chose to use heap sort because its
506
* worst behavior is n*log(n), unlike quicksort, whose worst
507
* behavior is n**2. */
508
/* self is not used */
510
PyArrayObject * alist, * ilist;
512
int len, i, * isubs, itmp;
514
Py_Try(PyArg_ParseTuple(args, "O", &list));
515
GET_ARR(alist,list,PyArray_DOUBLE,1);
517
Py_Try(ilist = (PyArrayObject *) PyArray_FromDims (1, &len, PyArray_INT));
518
isubs = (int *) A_DATA (ilist);
519
for ( i = 0 ; i < len ; i ++ )
522
data = (double *) A_DATA(alist) ;
523
/* now do heap sort on subscripts */
524
for (i = len / 2; i >= 0; i--) {
525
adjust (data, isubs, i, len) ;
527
for (i = len - 1; i >= 0; i-- )
530
isubs [i] = isubs [0] ;
532
adjust (data, isubs, 0, i) ;
536
return (PyObject *) ilist ;
540
binary_search(double dval, double dlist [], int len)
542
/* binary_search accepts three arguments: a numeric value and
543
* a numeric array and its length. It assumes that the array is sorted in
544
* increasing order. It returns the index of the array's
545
* largest element which is <= the value. It will return -1 if
546
* the value is less than the least element of the array. */
547
/* self is not used */
548
int bottom , top , middle, result;
550
if (dval < dlist [0])
555
while (bottom < top) {
556
middle = (top + bottom) / 2 ;
557
if (dlist [middle] < dval)
558
bottom = middle + 1 ;
559
else if (dlist [middle] > dval)
564
if (dlist [bottom] > dval)
565
result = bottom - 1 ;
574
binary_searchf(float dval, float dlist [], int len)
576
/* binary_search accepts three arguments: a numeric value and
577
* a numeric array and its length. It assumes that the array is sorted in
578
* increasing order. It returns the index of the array's
579
* largest element which is <= the value. It will return -1 if
580
* the value is less than the least element of the array. */
581
/* self is not used */
582
int bottom , top , middle, result;
584
if (dval < dlist [0])
589
while (bottom < top) {
590
middle = (top + bottom) / 2 ;
591
if (dlist [middle] < dval)
592
bottom = middle + 1 ;
593
else if (dlist [middle] > dval)
598
if (dlist [bottom] > dval)
599
result = bottom - 1 ;
606
/* return float, rather than double */
609
arr_interpf(PyObject *self, PyObject *args)
611
/* interp (y, x, z) treats (x, y) as a piecewise linear function
612
* whose value is y [0] for x < x [0] and y [len (y) -1] for x >
613
* x [len (y) -1]. An array of floats the same length as z is
614
* returned, whose values are ordinates for the corresponding z
615
* abscissae interpolated into the piecewise linear function. */
616
/* self is not used */
617
PyObject * oy, * ox, * oz ;
618
PyArrayObject * ay, * ax, * az , * _interp;
619
float * dy, * dx, * dz , * dres, * slopes;
620
int leny, lenz, i, left ;
622
PyObject *tpo = Py_None; /* unused argument, we've already parsed it*/
624
Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz, &tpo));
625
GET_ARR(ay,oy,PyArray_FLOAT,1);
626
GET_ARR(ax,ox,PyArray_FLOAT,1);
627
if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) {
628
SETERR ("interp: x and y are not the same length.");
632
GET_ARR2(az,oz,PyArray_FLOAT,1,MAX_INTERP_DIMS);
634
dy = (float *) A_DATA (ay);
635
dx = (float *) A_DATA (ax);
636
dz = (float *) A_DATA (az);
637
/* create output array with same size as 'Z' input array */
638
Py_Try (_interp = (PyArrayObject *) PyArray_FromDims
639
(A_NDIM(az), az->dimensions, PyArray_FLOAT));
640
dres = (float *) A_DATA (_interp) ;
641
slopes = (float *) malloc ( (leny - 1) * sizeof (float)) ;
642
for (i = 0 ; i < leny - 1; i++) {
643
slopes [i] = (dy [i + 1] - dy [i]) / (dx [i + 1] - dx [i]) ;
645
for ( i = 0 ; i < lenz ; i ++ )
647
left = binary_searchf (dz [i], dx, leny) ;
650
else if (left >= leny - 1)
651
dres [i] = dy [leny - 1] ;
653
dres [i] = slopes [left] * (dz [i] - dx [left]) + dy [left];
660
return PyArray_Return (_interp);
663
static char arr_interp__doc__[] =
664
"interp(y, x, z [,resulttypecode]) = y(z) interpolated by treating y(x) as piecewise fcn."
668
arr_interp(PyObject *self, PyObject *args)
670
/* interp (y, x, z) treats (x, y) as a piecewise linear function
671
* whose value is y [0] for x < x [0] and y [len (y) -1] for x >
672
* x [len (y) -1]. An array of floats the same length as z is
673
* returned, whose values are ordinates for the corresponding z
674
* abscissae interpolated into the piecewise linear function. */
675
/* self is not used */
676
PyObject * oy, * ox, * oz ;
677
PyArrayObject * ay, * ax, * az , * _interp;
678
double * dy, * dx, * dz , * dres, * slopes;
679
int leny, lenz, i, left ;
680
PyObject *tpo = Py_None;
681
char type_char = 'd';
682
char *type = &type_char;
684
Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz,&tpo));
685
if (tpo != Py_None) {
686
type = PyString_AsString(tpo);
693
return arr_interpf(self, args);
694
} else if (*type != 'd') {
695
SETERR ("interp: unimplemented typecode.");
698
GET_ARR(ay,oy,PyArray_DOUBLE,1);
699
GET_ARR(ax,ox,PyArray_DOUBLE,1);
700
if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) {
701
SETERR ("interp: x and y are not the same length.");
705
GET_ARR2(az,oz,PyArray_DOUBLE,1,MAX_INTERP_DIMS);
707
dy = (double *) A_DATA (ay);
708
dx = (double *) A_DATA (ax);
709
dz = (double *) A_DATA (az);
710
/* create output array with same size as 'Z' input array */
711
Py_Try (_interp = (PyArrayObject *) PyArray_FromDims
712
(A_NDIM(az), az->dimensions, PyArray_DOUBLE));
713
dres = (double *) A_DATA (_interp) ;
714
slopes = (double *) malloc ( (leny - 1) * sizeof (double)) ;
715
for (i = 0 ; i < leny - 1; i++) {
716
slopes [i] = (dy [i + 1] - dy [i]) / (dx [i + 1] - dx [i]) ;
718
for ( i = 0 ; i < lenz ; i ++ )
720
left = binary_search (dz [i], dx, leny) ;
723
else if (left >= leny - 1)
724
dres [i] = dy [leny - 1] ;
726
dres [i] = slopes [left] * (dz [i] - dx [left]) + dy [left];
733
return PyArray_Return (_interp);
736
static int incr_slot_ (float x, double *bins, int lbins)
739
for ( i = 0 ; i < lbins ; i ++ )
745
static int decr_slot_ (double x, double * bins, int lbins)
748
for ( i = lbins - 1 ; i >= 0; i -- )
754
static int monotonic_ (double * a, int lena)
759
("digitize: If a vector, second argument must have at least 2 elements.");
761
if (a [0] <= a [1]) /* possibly monotonic increasing */
763
for (i = 1 ; i < lena - 1; i ++)
764
if (a [i] > a [i + 1]) return 0 ;
767
else /* possibly monotonic decreasing */
769
for (i = 1 ; i < lena - 1; i ++)
770
if (a [i] < a [i + 1]) return 0 ;
775
static char arr_digitize__doc__[] =
780
arr_zmin_zmax(PyObject *self, PyObject *args)
782
/* zmin_zmax (z, ireg) returns a 2-tuple which consists
783
of the minimum and maximum values of z on the portion of the
784
mesh where ireg is not zero. z is a 2d array of Float and ireg
785
is an array of the same shape of Integer. By convention the first
786
row and column of ireg are zero, and the remaining entries are
787
used to determine which region each cell belongs to. A zero
788
entry says that this cell is excluded from the mesh. */
789
PyObject * zobj, * iregobj;
790
PyArrayObject * zarr, * iregarr;
791
double * z, zmin=0.0, zmax=0.0;
793
int have_min_max = 0;
796
Py_Try(PyArg_ParseTuple(args, "OO", &zobj, &iregobj));
797
GET_ARR (zarr, zobj, PyArray_DOUBLE, 2);
798
if (! (iregarr = (PyArrayObject *)PyArray_ContiguousFromObject(iregobj,
799
PyArray_INT, 2, 2))) {
803
n = A_DIM (iregarr, 0);
804
m = A_DIM (iregarr, 1);
805
if (n != A_DIM (zarr, 0) || m != A_DIM (zarr, 1)) {
806
SETERR ("zmin_zmax: z and ireg do not have the same shape.");
811
ireg = (int *) A_DATA (iregarr);
812
z = (double *) A_DATA (zarr);
813
k = 0; /* k is i * m + j */
814
for ( i = 0; i < n; i++) {
815
for (j = 0; j < m; j++) {
816
if ( (ireg [k] != 0) ||
817
(i != n - 1 && j != m - 1 &&
818
(ireg [k + m] != 0 || ireg [k + 1] != 0 || ireg [k + m + 1] != 0))) {
819
if (! have_min_max ) {
824
if (z [k] < zmin) zmin = z [k];
825
else if (z [k] > zmax) zmax = z [k];
834
SETERR ("zmin_zmax: unable to calculate zmin and zmax!");
837
return Py_BuildValue ("dd", zmin, zmax);
840
static char arr_zmin_zmax__doc__[] =
845
arr_digitize(PyObject *self, PyObject *args)
847
/* digitize (x, bins) returns an array of python integers the same
848
length of x (if x is a one-dimensional array), or just an integer
849
(if x is a scalar). The values i returned are such that
850
bins [i - 1] <= x < bins [i] if bins is monotonically increasing,
851
or bins [i - 1] > x >= bins [i] if bins is monotonically decreasing.
852
Beyond the bounds of bins, returns either i = 0 or i = len (bins)
854
/* self is not used */
855
PyObject * ox, * obins ;
856
PyArrayObject *ax=NULL, *abins=NULL, *aret ;
857
double x=0.0, bins=0.0; /* if either or both is a scalar */
858
double *dx=NULL, *dbins=NULL; /* if either or both is a vector */
859
int lbins=0, lx ; /* lengths, if vectors */
862
int x_is_scalar, bins_is_scalar ;
864
Py_Try(PyArg_ParseTuple(args, "OO", & ox, & obins));
865
x_is_scalar = ! isARRAY (ox) ;
866
bins_is_scalar = ! isARRAY (obins) ;
869
GET_ARR(ax,ox,PyArray_DOUBLE,1);
870
if (A_NDIM (ax) > 1) {
871
SETERR ("digitize: first argument has too many dimensions.") ;
875
dx = (double *) A_DATA (ax) ;
879
if (PyInt_Check (ox))
880
x = (double) PyInt_AsLong (ox) ;
881
else if (PyFloat_Check (ox))
882
x = PyFloat_AS_DOUBLE ((PyFloatObject *)ox) ;
884
SETERR ("digitize: bad type for first argument.") ;
887
if ( !bins_is_scalar )
889
GET_ARR(abins,obins,PyArray_DOUBLE,1);
890
if (A_NDIM (abins) > 1) {
891
SETERR ("digitize: second argument has too many dimensions.") ;
895
lbins = A_SIZE (abins) ;
896
dbins = (double *) A_DATA (abins) ;
900
if (PyInt_Check (obins))
901
bins = (double) PyInt_AsLong (obins) ;
902
else if (PyFloat_Check (obins))
903
bins = PyFloat_AS_DOUBLE ((PyFloatObject *)obins) ;
905
SETERR ("digitize: bad type for second argument.") ;
912
return PyInt_FromLong (0) ;
914
return PyInt_FromLong (1) ;
917
aret = (PyArrayObject *) PyArray_FromDims (1, &lx, PyArray_LONG) ;
918
iret = (long *) A_DATA (aret) ;
919
for ( i = 0 ; i < lx ; i ++ )
921
iret [i] = (long) 1 ;
925
m = monotonic_ (dbins, lbins) ;
929
return PyInt_FromLong (decr_slot_ ((float)x, dbins, lbins)) ;
930
aret = (PyArrayObject *) PyArray_FromDims (1, &lx, PyArray_LONG) ;
931
iret = (long *) A_DATA (aret) ;
932
for ( i = 0 ; i < lx ; i ++ )
933
iret [i] = (long) decr_slot_ (dx [i], dbins, lbins) ;
938
return PyInt_FromLong (incr_slot_ ((float)x, dbins, lbins)) ;
939
aret = (PyArrayObject *) PyArray_FromDims (1, &lx, PyArray_LONG) ;
940
iret = (long *) A_DATA (aret) ;
941
for ( i = 0 ; i < lx ; i ++ )
942
iret [i] = (long) incr_slot_ ((float)dx [i], dbins, lbins) ;
946
SETERR ("digitize: Second argument must be monotonic.") ;
955
return PyArray_Return (aret) ;
958
static char arr_reverse__doc__[] =
963
arr_reverse(PyObject *self, PyObject *args)
965
/* reverse (x, n) returns a PyFloat Matrix the same size and shape as
966
x, but with the elements along the nth dimension reversed.
967
x must be two-dimensional. */
968
/* self is not used */
971
PyArrayObject * ax, * ares ;
973
int d0, d1, dims [2] ;
976
Py_Try(PyArg_ParseTuple(args, "Oi", &ox, &n));
977
if (n != 0 && n != 1) {
978
SETERR("reverse: Second argument must be 0 or 1.");
980
GET_ARR(ax,ox,PyArray_DOUBLE,2);
981
dx = (double *) A_DATA (ax);
982
d0 = dims [0] = A_DIM (ax, 0);
983
d1 = dims [1] = A_DIM (ax, 1);
984
Py_Try(ares = (PyArrayObject *) PyArray_FromDims (2, dims, PyArray_DOUBLE));
985
dres = (double *) A_DATA (ares);
986
if (n == 0) /* reverse the columns */
987
for (i = 0; i < d1 ; i ++ )
988
{for (jl = i, jh = (d0 - 1) * d1 + i; jl < jh; jl += d1, jh -= d1)
990
dres [jl] = dx [jh] ;
991
dres [jh] = dx [jl] ;
993
if (jl == jh) dres [jl] = dx [jl] ;
995
else /* reverse the rows */
996
for (i = 0; i < d0 ; i ++ )
997
{for (jl = i * d1, jh = (i + 1) * d1 - 1; jl < jh; jl +=1, jh -= 1)
999
dres [jl] = dx [jh] ;
1000
dres [jh] = dx [jl] ;
1002
if (jl == jh) dres [jl] = dx [jl] ;
1006
return PyArray_Return (ares) ;
1009
static char arr_span__doc__[] =
1014
arr_span(PyObject *self, PyObject *args)
1016
/* span (lo, hi, num, d2 = 0) returns an array of num equally
1017
spaced PyFloats starting with lo and ending with hi. if d2 is
1018
not zero, it will return a two-dimensional array, each of the
1019
d2 rows of which is the array of equally spaced numbers. */
1020
/* self is not used */
1024
double lo, hi, * drow, * dres;
1026
PyArrayObject * arow, * ares ;
1028
Py_Try(PyArg_ParseTuple(args, "ddi|i", &lo, &hi, &num, &d2));
1031
Py_Try(arow = (PyArrayObject *) PyArray_FromDims (1, &num, PyArray_DOUBLE));
1032
drow = (double *) A_DATA (arow) ;
1033
for (i = 0; i < num; i++)
1034
drow [i] = lo + ( (double) i) * (hi - lo) / ( (double) (num - 1)) ;
1036
return PyArray_Return (arow) ;
1040
(ares = (PyArrayObject *) PyArray_FromDims (2, dims, PyArray_DOUBLE));
1041
dres = (double *) A_DATA (ares) ;
1042
for ( id2 = 0 ; id2 < num * d2 ; id2 += num )
1043
for (j = 0; j < num; j ++ )
1044
dres [id2 + j] = drow [j] ;
1048
return PyArray_Return (ares) ;
1051
static char arr_nz__doc__ [] =
1056
arr_nz (PyObject *self, PyObject *args)
1058
/* nz_ (x): x is an array of unsigned bytes. If x
1059
ends with a bunch of zeros, this returns with the index of
1060
the first zero element after the last nonzero element.
1061
It returns the length of the array if its last element
1062
is nonzero. This is essentially the "effective length"
1064
/* self is not used */
1066
unsigned char * cdat;
1068
PyArrayObject * adat;
1070
Py_Try(PyArg_ParseTuple(args, "O", &odat)) ;
1071
GET_ARR(adat,odat,PyArray_UBYTE,1) ;
1072
cdat = (unsigned char *) A_DATA (adat) ;
1073
len = A_SIZE (adat) ;
1074
for (i = len; i > 0; i --)
1075
if (cdat [i - 1] != (unsigned char) 0) break ;
1077
return PyInt_FromLong ( (long) i) ;
1080
static char arr_find_mask__doc__ [] =
1084
static PyObject * arr_find_mask (PyObject * self, PyObject * args)
1086
/* find_mask (fs, node_edges): This function is used to calculate
1087
a mask whose corresponding entry is 1 precisely if an edge
1088
of a cell is cut by an isosurface, i. e., if the function
1089
fs is one on one of the two vertices of an edge and zero
1090
on the other (fs = 1 represents where some function on
1091
the mesh was found to be negative by the calling routine).
1092
fs is ntotal by nv, where nv is the number of vertices
1093
of a cell (4 for a tetrahedren, 5 for a pyramid, 6 for a prism).
1094
node_edges is a nv by ne array, where ne is the number of
1095
edges on a cell (6 for a tet, 8 for a pyramid, 9 for a prism).
1096
The entries in each row are 1 precisely if the corresponding edge
1097
is incident on the vertex. The exclusive or of the rows
1098
which correspond to nonzero entries in fs contains 1 in
1099
entries corresponding to edges where fs has opposite values
1102
PyObject * fso, * node_edgeso ;
1103
PyArrayObject * fsa, * node_edgesa, * maska ;
1104
int * fs, * node_edges, * mask ;
1105
int i, j, k, l, ll, ifs, imask, ntotal, ne, nv, ans_size ;
1107
Py_Try (PyArg_ParseTuple ( args, "OO", & fso, & node_edgeso ) ) ;
1108
GET_ARR (fsa, fso, PyArray_INT, 2) ;
1109
GET_ARR (node_edgesa, node_edgeso, PyArray_INT, 2) ;
1110
fs = (int *) A_DATA (fsa) ;
1111
node_edges = (int *) A_DATA (node_edgesa) ;
1112
ntotal = A_DIM (fsa, 0) ;
1113
nv = A_DIM (fsa, 1) ;
1114
if ( nv != A_DIM (node_edgesa, 0) ) {
1115
SETERR ("2nd dimension of 1st arg and 1st dimension of 2nd not equal.");
1117
Py_DECREF (node_edgesa) ;
1120
ne = A_DIM (node_edgesa, 1) ;
1121
ans_size = ntotal * ne ;
1122
Py_Try (maska = (PyArrayObject *) PyArray_FromDims
1123
(1, & ans_size, PyArray_INT)) ;
1124
mask = (int *) A_DATA (maska) ;
1126
for (i = 0, ifs = 0, imask = 0 ; i < ntotal ;
1127
i ++, imask += ne, ifs += nv) {
1128
for (j = ifs, k = 0; k < nv; j ++, k ++) {
1130
for ( l = imask , ll = 0; ll < ne ; l ++ , ll ++) {
1131
mask [l] ^= node_edges [j % nv * ne + ll] ;
1137
return PyArray_Return (maska) ;
1141
/* Data for construct3 and walk3 */
1142
int start_face4 [] = {0, 0, 1, 0, 2, 1} ;
1143
int start_face5 [] = {0, 0, 1, 2, 0, 1, 2, 3} ;
1144
int start_face6 [] = {1, 1, 0, 0, 2, 2, 0, 0, 1} ;
1145
int start_face8 [] = {0, 1, 0, 1, 0, 1, 0, 1, 2, 3, 2, 3} ;
1146
static int * start_face [4] = {start_face4, start_face5, start_face6,
1149
int ef0 [] = {0, 1} ;
1150
int ef1 [] = {0, 2} ;
1151
int ef2 [] = {0, 3} ;
1152
int ef3 [] = {0, 4} ;
1153
int ef4 [] = {0, 5} ;
1154
int ef5 [] = {1, 2} ;
1155
int ef6 [] = {1, 3} ;
1156
int ef7 [] = {1, 4} ;
1157
int ef8 [] = {1, 5} ;
1158
int ef9 [] = {2, 3} ;
1159
int ef10 [] = {2, 4} ;
1160
int ef11 [] = {2, 5} ;
1161
int ef12 [] = {3, 4} ;
1162
int ef13 [] = {3, 5} ;
1163
int * edge_faces4 [] = {ef0, ef1, ef5, ef2, ef9, ef6} ;
1164
int * edge_faces5 [] = {ef2, ef0, ef5, ef9, ef3, ef7, ef10, ef12} ;
1165
int * edge_faces6 [] = {ef6, ef7, ef2, ef3, ef9, ef10, ef0, ef1, ef5} ;
1166
int * edge_faces8 [] = {ef1, ef5, ef2, ef6, ef3, ef7, ef4, ef8, ef10,
1168
static int ** edge_faces [] = {edge_faces4, edge_faces5, edge_faces6,
1171
int fe40 [] = {0, 1, 3} ;
1172
int fe41 [] = {0, 5, 2} ;
1173
int fe42 [] = {1, 2, 4} ;
1174
int fe43 [] = {3, 4, 5} ;
1175
int * face_edges4 [] = {fe40, fe41, fe42, fe43} ;
1176
int fe50 [] = {0, 1, 4} ;
1177
int fe51 [] = {1, 2, 5} ;
1178
int fe52 [] = {2, 3, 6} ;
1179
int fe53 [] = {0, 7, 3} ;
1180
int fe54 [] = {4, 5, 6, 7} ;
1181
int * face_edges5 [] = {fe50, fe51, fe52, fe53, fe54} ;
1182
int fe60 [] = {2, 7, 3, 6} ;
1183
int fe61 [] = {0, 6, 1, 8} ;
1184
int fe62 [] = {4, 8, 5, 7} ;
1185
int fe63 [] = {0, 4, 2} ;
1186
int fe64 [] = {1, 3, 5} ;
1187
int * face_edges6 [] = {fe60, fe61, fe62, fe63, fe64} ;
1188
int fe80 [] = {0, 6, 2, 4} ;
1189
int fe81 [] = {1, 5, 3, 7} ;
1190
int fe82 [] = {0, 8, 1, 10} ;
1191
int fe83 [] = {2, 11, 3, 9} ;
1192
int fe84 [] = {4, 9, 5, 8} ;
1193
int fe85 [] = {6, 10, 7, 11} ;
1194
int * face_edges8 [] = {fe80, fe81, fe82, fe83, fe84, fe85} ;
1195
static int ** face_edges [] = {face_edges4, face_edges5, face_edges6,
1199
int lens4 [] = {3, 3, 3, 3} ;
1200
int lens5 [] = {3, 3, 3, 3, 4} ;
1201
int lens6 [] = {4, 4, 4, 3, 3} ;
1202
int lens8 [] = {4, 4, 4, 4, 4, 4} ;
1203
static int * lens [] = {lens4, lens5, lens6, lens8} ;
1205
static int no_edges [4] = {6, 8, 9, 12} ;
1206
/* static int no_verts [4] = {4, 5, 6, 8} ; */
1207
static int powers [4] = {14, 30, 62, 254} ;
1211
static void walk3 ( int * permute, int * mask, int itype, int pt )
1225
for (i = 0; i < 12; i++) splits [i] = 0 ;
1226
/* list = mask.nonzero () */
1227
for (i = 0, nlist = 0 ; i < no_edges [itype] ; i ++)
1229
list [nlist++] = i ;
1234
face = start_face [itype] [edge] ;
1235
for (i = 0 ; i < nlist - 1 ; i ++ )
1237
/* record this edge */
1238
permute [edge * pt] = i ;
1239
splits [edge] = split ;
1241
/* advance to next edge */
1242
ttry = face_edges [itype] [face] ;
1243
lttry = lens [itype] [face] ;
1245
for ( j = 1 ; j < lttry ; j ++ )
1246
if (abs (ttry [now] - edge) > abs (ttry [j] - edge))
1249
/* test opposite edge first (make X, never // or \\) */
1250
edge = ttry [(now + 1) % lttry] ;
1253
/* otherwise one or the other opposite edges must be set */
1254
edge = ttry [now % lttry] ;
1257
edge = ttry [ (now + 2) % lttry] ;
1261
for (edge = 0 ; edge < no_edges [itype] ; edge++)
1263
if ( mask [edge] != 0 )
1271
ttry = edge_faces [itype] [edge] ;
1272
face = (face == ttry [0]) ? ttry [1] : ttry [0] ;
1274
permute [edge * pt] = nlist - 1 ;
1275
splits [edge] = split ;
1278
for ( i = 0 , j = 0 ; i < no_edges [itype] ; i ++ , j += pt)
1280
permute [j] += no_edges [itype] * splits [i] ;
1285
static char arr_construct3__doc__ [] =
1289
arr_construct3 (PyObject * self, PyObject * args)
1290
{ /* construct3 (mask, itype) computes how the cut
1291
edges of a particular type of cell must be ordered so
1292
that the polygon of intersection can be drawn correctly.
1293
itype = 0 for tetrahedra; 1 for pyramids; 2 for prisms;
1294
3 for hexahedra. Suppose nv is the number of vertices
1295
of the cell type, and ne is the number of edges. Mask
1296
has been ravelled so that it was flat, originally it
1297
had 2**nv-2 rows, each with ne entries. Each row is
1298
ne long, and has an entry of 1 corresponding to each
1299
edge that is cut when the set of vertices corresponding
1300
to the row index has negative values. (The binary number
1301
for the row index + 1 has a one in position i if vertex
1302
i has a negative value.) The return array permute is
1303
ne by 2**nv-2, and the rows of permute tell how
1304
the edges should be ordered to draw the polygon properly. */
1307
PyArrayObject * permutea, * maska ;
1308
int itype, ne, pt, nm ;
1309
int * permute, * mask ;
1310
int permute_dims [2] ;
1313
/* dbg = fopen ("dbg","w"); */
1314
Py_Try (PyArg_ParseTuple ( args, "Oi", & masko, & itype ) ) ;
1315
GET_ARR (maska, masko, PyArray_INT, 1) ;
1316
mask = (int *) A_DATA (maska) ;
1317
permute_dims [0] = ne = no_edges [itype] ;
1318
permute_dims [1] = pt = powers [itype] ;
1319
nm = A_DIM (maska, 0) ;
1320
if ( ne * pt != nm ) {
1321
SETERR ("permute and mask must have same number of elements.") ;
1326
(PyArrayObject *) PyArray_FromDims (2, permute_dims, PyArray_INT));
1327
permute = (int *) A_DATA (permutea) ;
1328
for ( i = 0 ; i < pt ; i ++, permute ++, mask += ne )
1330
walk3 (permute, mask, itype, pt) ;
1333
/* fclose (dbg) ; */
1334
return PyArray_Return (permutea) ;
1337
static char arr_to_corners__doc__ [] =
1341
arr_to_corners (PyObject * self, PyObject * args)
1343
/* This routine takes an array of floats describing cell-centered
1344
values and expands it to node-centered values. */
1345
PyObject * oarr, * onv;
1347
PyArrayObject * aarr, *anv, *ares;
1348
int * nv , i, j, snv, jtot;
1349
double * arr, * res;
1351
Py_Try (PyArg_ParseTuple (args, "OOi", & oarr, & onv, & sum_nv));
1352
GET_ARR (aarr, oarr, PyArray_DOUBLE, 1) ;
1354
GET_ARR (anv, onv, PyArray_INT, 1) ;
1357
ERRSS ("The second argument must be an Int array") ;
1361
snv = A_SIZE (anv) ;
1362
if (snv != A_SIZE (aarr)) {
1363
ERRSS ("The first and second arguments must be the same size.") ;
1368
if ( ! (ares = (PyArrayObject *)
1369
PyArray_FromDims (1, & sum_nv, PyArray_DOUBLE))) {
1370
ERRSS ("Unable to create result array.") ;
1375
res = (double *) A_DATA (ares) ;
1376
arr = (double *) A_DATA (aarr) ;
1377
nv = (int *) A_DATA (anv) ;
1379
for ( i = 0; i < snv; i++) {
1380
for (j = 0; j < nv [i]; j++) {
1381
res [j + jtot] = arr [i];
1383
jtot = jtot + nv [i];
1389
return PyArray_Return (ares) ;
1392
/* List of methods defined in the module */
1394
static struct PyMethodDef arr_methods[] = {
1395
{"histogram", arr_histogram, 1, arr_histogram__doc__},
1396
{"array_set", arr_array_set, 1, arr_array_set__doc__},
1397
{"index_sort", arr_index_sort, 1, arr_index_sort__doc__},
1398
{"interp", arr_interp, 1, arr_interp__doc__},
1399
{"digitize", arr_digitize, 1, arr_digitize__doc__},
1400
{"zmin_zmax", arr_zmin_zmax, 1, arr_zmin_zmax__doc__},
1401
{"reverse", arr_reverse, 1, arr_reverse__doc__},
1402
{"span", arr_span, 1, arr_span__doc__},
1403
{"nz", arr_nz, 1, arr_nz__doc__},
1404
{"find_mask",arr_find_mask, 1, arr_find_mask__doc__},
1405
{"construct3", arr_construct3, 1, arr_construct3__doc__},
1406
{"to_corners", arr_to_corners, 1, arr_to_corners__doc__},
1408
{NULL, NULL} /* sentinel */
1412
/* Initialization function for the module (*must* be called initarrayfns) */
1414
static char arrayfns_module_documentation[] =
1423
/* Create the module and add the functions */
1424
m = Py_InitModule4("gistfuncs", arr_methods,
1425
arrayfns_module_documentation,
1427
PYTHON_API_VERSION);
1429
/* Add some symbolic constants to the module */
1430
d = PyModule_GetDict(m);
1431
ErrorObject = PyErr_NewException("gistfuncs.error", NULL, NULL);
1432
PyDict_SetItemString(d, "error", ErrorObject);
1434
/* XXXX Add constants here */
added
removed
scipy/sandbox/xplt/pygist/gistfuncsmodule.c
