[Matrix-SIG] Re: NumPy changes
Paul F. Dubois
dubois1@llnl.gov
Tue, 17 Aug 1999 09:13:59 -0700
Lee Taylor sent me some changes we discussed so that PyArray_From...
routines are factored so that one has the ability to create the array
descriptor oneself. This proved useful to Lee in a database handler.
I updated the repo with Lee's changes after testing on Windows. However, I
did not include his change about DEFINE_EXTENSIONCLASS because they don't
work on Windows. I think the difference is that Lee loads his code with
Numerical statically, while David and I have been testing dynamically. The
static load gets multiple declarations of the same external.
Clearly something is wrong here. Maybe we need help. Konrad, this CAPI trick
was your idea, wasn't it? Do you see what to do? Maybe this is related to
the bug report about reloads causing trouble. (For the purposes of
understanding the problem, the beta release is fine; Lee's changes are
unrelated.)
----- Original Message -----
From: Lee Taylor <taylor@rhino.llnl.gov>
To: Paul Dubois <dubois1@llnl.gov>
Sent: Friday, August 13, 1999 4:15 PM
Subject: NumPy changes
>
> Paul,
> Attached are the changes I talked about to allow arrays
> of structures to be index via NumPy.
>
> I don't think that this will address the bit array problem.
> All of the indexing is done with byte addresses because that's
> the smallest item a pointer can point at. The bit array would
> require some shifting and masking.
>
> Also, the IBM compiler comments about line 1607 in arrayobject.c
> having a /* within a comment.
>
> Finally, these files have the DEFINE_EXTENSIONCLASS removed
> from arrayobject.h and _numpymodule.c as describe in my mail
> to the MATRIX-sig.
>
> Lee
>
>
----------------------------------------------------------------------------
----
> #include "Python.h"
>
> #define _ARRAY_MODULE
> #include "arrayobject.h"
> #include "ExtensionClassMacros.h"
> #define _UFUNC_MODULE
> #include "ufuncobject.h"
>
> /* Table of functions defined in the module */
>
> static PyMethodDef numpy_methods[] = {
> {NULL, NULL} /* sentinel */
> };
>
> /* Module initialization */
>
> void
> init_numpy()
> {
> PyObject *m, *d;
> static void *PyArray_API[PyArray_API_pointers];
> static void *Py_UFunc_API[PyUFunc_API_pointers];
>
> /* Create the module and add the functions */
> PyArray_Type.ob_type = &PyType_Type;
> PyUFunc_Type.ob_type = &PyType_Type;
>
> m = Py_InitModule("_numpy", numpy_methods);
> d = PyModule_GetDict(m);
>
> /* Initialize C API pointer arrays and store them in module */
> PyArray_API[PyArray_Type_NUM] = (void *)&PyArray_Type;
> PyArray_API[PyArray_SetNumericOps_NUM] = (void *)&PyArray_SetNumericOps;
> PyArray_API[PyArray_INCREF_NUM] = (void *)&PyArray_INCREF;
> PyArray_API[PyArray_XDECREF_NUM] = (void *)&PyArray_XDECREF;
> #if 0
> PyArray_API[PyArrayError_NUM] = (void *)&PyArrayError;
> #endif
> PyArray_API[PyArray_SetStringFunction_NUM] =
> (void *)&PyArray_SetStringFunction;
> PyArray_API[PyArray_DescrFromType_NUM] = (void *)&PyArray_DescrFromType;
> PyArray_API[PyArray_Cast_NUM] = (void *)&PyArray_Cast;
> PyArray_API[PyArray_CanCastSafely_NUM] = (void *)&PyArray_CanCastSafely;
> PyArray_API[PyArray_ObjectType_NUM] = (void *)&PyArray_ObjectType;
> PyArray_API[_PyArray_multiply_list_NUM] = (void
*)&_PyArray_multiply_list;
> PyArray_API[PyArray_Size_NUM] = (void *)&PyArray_Size;
> PyArray_API[PyArray_FromDims_NUM] = (void *)&PyArray_FromDims;
> PyArray_API[PyArray_FromDimsAndData_NUM] = (void
*)&PyArray_FromDimsAndData;
> PyArray_API[PyArray_FromDimsAndDataAndDescr_NUM] = (void
*)&PyArray_FromDimsAndDataAndDescr;
> PyArray_API[PyArray_ContiguousFromObject_NUM] =
> (void *)&PyArray_ContiguousFromObject;
> PyArray_API[PyArray_CopyFromObject_NUM] = (void
*)&PyArray_CopyFromObject;
> PyArray_API[PyArray_FromObject_NUM] = (void *)&PyArray_FromObject;
> PyArray_API[PyArray_Return_NUM] = (void *)&PyArray_Return;
> PyArray_API[PyArray_Reshape_NUM] = (void *)&PyArray_Reshape;
> PyArray_API[PyArray_Copy_NUM] = (void *)&PyArray_Copy;
> PyArray_API[PyArray_Take_NUM] = (void *)&PyArray_Take;
> PyArray_API[PyArray_As1D_NUM] = (void *)&PyArray_As1D;
> PyArray_API[PyArray_As2D_NUM] = (void *)&PyArray_As2D;
> PyArray_API[PyArray_Free_NUM] = (void *)&PyArray_Free;
> PyDict_SetItemString(d, "_ARRAY_API",
> PyCObject_FromVoidPtr((void *)PyArray_API, NULL));
>
> Py_UFunc_API[PyUFunc_Type_NUM] = (void *)&PyUFunc_Type;
> Py_UFunc_API[PyUFunc_FromFuncAndData_NUM] = (void
*)&PyUFunc_FromFuncAndData;
> Py_UFunc_API[PyUFunc_GenericFunction_NUM] = (void
*)&PyUFunc_GenericFunction;
> Py_UFunc_API[PyUFunc_f_f_As_d_d_NUM] = (void *)&PyUFunc_f_f_As_d_d;
> Py_UFunc_API[PyUFunc_d_d_NUM] = (void *)&PyUFunc_d_d;
> Py_UFunc_API[PyUFunc_F_F_As_D_D_NUM] = (void *)&PyUFunc_F_F_As_D_D;
> Py_UFunc_API[PyUFunc_D_D_NUM] = (void *)&PyUFunc_D_D;
> Py_UFunc_API[PyUFunc_O_O_NUM] = (void *)&PyUFunc_O_O;
> Py_UFunc_API[PyUFunc_ff_f_As_dd_d_NUM] = (void *)&PyUFunc_ff_f_As_dd_d;
> Py_UFunc_API[PyUFunc_dd_d_NUM] = (void *)&PyUFunc_dd_d;
> Py_UFunc_API[PyUFunc_FF_F_As_DD_D_NUM] = (void *)&PyUFunc_FF_F_As_DD_D;
> Py_UFunc_API[PyUFunc_DD_D_NUM] = (void *)&PyUFunc_DD_D;
> Py_UFunc_API[PyUFunc_OO_O_NUM] = (void *)&PyUFunc_OO_O;
> Py_UFunc_API[PyUFunc_O_O_method_NUM] = (void *)&PyUFunc_O_O_method;
> #if 0
> Py_UFunc_API[PyArray_Map_NUM] = (void *)&PyArray_Map;
> #endif
> PyDict_SetItemString(d, "_UFUNC_API",
> PyCObject_FromVoidPtr((void *)Py_UFunc_API, NULL));
>
> PyExtensionClass_Export(d, "multiarray", PyArray_Type);
>
> PyDict_SetItemString(d, "_UseExtensionClass", PyInt_FromLong(1));
> /* Check for errors */
> if (PyErr_Occurred())
> Py_FatalError("can't initialize module _numpy");
> }
>
----------------------------------------------------------------------------
----
> /*
> Python Array Object -- Provide multidimensional arrays as a basic
> object type in python.
>
> Copyright (c) 1995, 1996, 1997 Jim Hugunin, hugunin@mit.edu
> See file COPYING for details.
>
>
> These arrays are primarily designed for supporting multidimensional,
> homogeneous arrays of basic C numeric types. They also can support
> arrays of arbitrary Python Objects, if you are willing to sacrifice
> performance for heterogeneity.
> */
>
> /* $Id: arrayobject.c,v 0.16 1999/07/11 21:42:29 barnard Exp $ */
>
> #include "Python.h"
> /*Silly trick to make dll library work right under Win32*/
> #ifdef MS_WIN32
> #undef DL_IMPORT
> #define DL_IMPORT(RTYPE) __declspec(dllexport) RTYPE
> #endif
> #define _ARRAY_MODULE
> #define _UFUNC_MODULE
> #include "arrayobject.h"
> #include "ExtensionClassMacros.h"
> #include "ufuncobject.h"
>
> #define OBJECT(O) ((PyObject*)(O))
>
> /* There are several places in the code where an array of dimensions is */
> /* allocated statically. This is the size of that static allocation. I
*/
> /* can't come up with any reasonable excuse for a larger array than this.
*/
>
> #define MAX_DIMS 40
>
> /* Helper Functions */
> extern int _PyArray_multiply_list(int *l1, int n) {
> int s=1, i=0;
> while (i < n) s *= l1[i++];
> return s;
> }
> extern int _PyArray_compare_lists(int *l1, int *l2, int n) {
> int i;
> for(i=0;i<n;i++) {
> if (l1[i] != l2[i]) return 0;
> }
> return 1;
> }
>
> /* These can be cleaned up */
> #define SIZE(mp) (_PyArray_multiply_list((mp)->dimensions, (mp)->nd))
> #define NBYTES(mp) ((mp)->descr->elsize * SIZE(mp))
> /* Obviously this needs some work. */
> #define ISCONTIGUOUS(m) ((m)->flags & CONTIGUOUS)
>
> #define PyArray_CONTIGUOUS(m) (ISCONTIGUOUS(m) ? Py_INCREF(m), m : \
> (PyArrayObject *)(PyArray_ContiguousFromObject((PyObject *)(m),
(m)->descr->type_num, 0,0)))
>
> int do_sliced_copy(char *dest, int *dest_strides, int *dest_dimensions,
int dest_nd,
> char *src, int *src_strides, int *src_dimensions, int src_nd,
> int elsize, int copies) {
> int i, j;
>
> if (src_nd == 0 && dest_nd == 0) {
> for(j=0; j<copies; j++) {
> memcpy(dest, src, elsize);
> dest += elsize;
> }
> return 0;
> }
>
> if (dest_nd > src_nd) {
> for(i=0; i<*dest_dimensions; i++, dest += *dest_strides) {
> if (do_sliced_copy(dest, dest_strides+1, dest_dimensions+1, dest_nd-1,
> src, src_strides, src_dimensions, src_nd, elsize, copies) == -1)
return -1;
> }
> return 0;
> }
>
> if (dest_nd == 1) {
> if (*dest_dimensions != *src_dimensions) {
> PyErr_SetString(PyExc_ValueError, "matrices are not aligned for copy");
> return -1;
> }
> for(i=0; i<*dest_dimensions; i++, src += *src_strides) {
> for(j=0; j<copies; j++) {
> memcpy(dest, src, elsize);
> dest += *dest_strides;
> }
> }
> } else {
> for(i=0; i<*dest_dimensions; i++, dest += *dest_strides, src +=
*src_strides) {
> if (do_sliced_copy(dest, dest_strides+1, dest_dimensions+1, dest_nd-1,
> src, src_strides+1, src_dimensions+1, src_nd-1, elsize, copies) == -1)
return -1;
> }
> }
> return 0;
> }
>
> int optimize_slices(int **dest_strides, int **dest_dimensions, int
*dest_nd,
> int **src_strides, int **src_dimensions, int *src_nd,
> int *elsize, int *copies) {
> while (*src_nd > 0) {
> if (((*dest_strides)[*dest_nd-1] == *elsize) && ((*src_strides)[*src_nd-1]
== *elsize)) {
> if ((*dest_dimensions)[*dest_nd-1] != (*src_dimensions)[*src_nd-1]) {
> PyErr_SetString(PyExc_ValueError, "matrices are not aligned for copy");
> return -1;
> }
> *elsize *= (*dest_dimensions)[*dest_nd-1];
> *dest_nd-=1; *src_nd-=1;
> } else {
> break;
> }
> }
> if (*src_nd == 0) {
> while (*dest_nd > 0) {
> if (((*dest_strides)[*dest_nd-1] == *elsize)) {
> *copies *= (*dest_dimensions)[*dest_nd-1];
> *dest_nd-=1;
> } else {
> break;
> }
> }
> }
> return 0;
> }
>
> static char *contiguous_data(PyArrayObject *src) {
> int dest_strides[MAX_DIMS], *dest_strides_ptr;
> int *dest_dimensions=src->dimensions;
> int dest_nd=src->nd;
> int *src_strides = src->strides;
> int *src_dimensions=src->dimensions;
> int src_nd=src->nd;
> int elsize=src->descr->elsize;
> int copies=1;
> int ret, i;
> int stride=elsize;
> char *new_data;
>
> for(i=dest_nd-1; i>=0; i--) {
> dest_strides[i] = stride;
> stride *= dest_dimensions[i];
> }
>
> dest_strides_ptr = dest_strides;
>
> if (optimize_slices(&dest_strides_ptr, &dest_dimensions, &dest_nd,
&src_strides, &src_dimensions, &src_nd,
> &elsize, &copies) == -1) return NULL;
>
> new_data = (char *)malloc(stride);
>
> ret = do_sliced_copy(new_data, dest_strides_ptr, dest_dimensions, dest_nd,
> src->data, src_strides, src_dimensions, src_nd,
> elsize, copies);
>
> if (ret != -1) { return new_data; }
> else { free(new_data); return NULL; }
> }
>
>
> /* Used for arrays of python objects to increment the reference count of
*/
> /* every python object in the array. */
> extern int PyArray_INCREF(PyArrayObject *mp) {
> int i, n;
> PyObject **data;
>
> if (mp->descr->type_num != PyArray_OBJECT) return 0;
>
> if (ISCONTIGUOUS(mp)) {
> data = (PyObject **)mp->data;
> } else {
> if ((data = (PyObject **)contiguous_data(mp)) == NULL) return -1;
> }
>
> n = SIZE(mp);
> for(i=0; i<n; i++, data++) Py_XINCREF(*data);
>
> if (!ISCONTIGUOUS(mp)) free(data);
>
> return 0;
> }
>
> extern int PyArray_XDECREF(PyArrayObject *mp) {
> int i, n;
> PyObject **data;
>
> if (mp->descr->type_num != PyArray_OBJECT) return 0;
>
> if (ISCONTIGUOUS(mp)) {
> data = (PyObject **)mp->data;
> } else {
> if ((data = (PyObject **)contiguous_data(mp)) == NULL) return -1;
> }
>
> n = SIZE(mp);
> for(i=0; i<n; i++, data++) Py_XDECREF(*data);
>
> if (!ISCONTIGUOUS(mp)) free(data);
>
> return 0;
> }
>
> /* Including a C file is admittedly ugly. Look at the C file if you want
to */
> /* see something even uglier. This is computer generated code. */
> #include "arraytypes.c"
>
> char *index2ptr(PyArrayObject *mp, int i) {
> if (i==0 && (mp->nd == 0 || mp->dimensions[0] > 0)) return mp->data;
>
> if (mp->nd>0 && i>0 && i < mp->dimensions[0]) {
> return mp->data+i*mp->strides[0];
> }
> PyErr_SetString(PyExc_IndexError,"index out of bounds");
> return NULL;
> }
>
> extern int PyArray_Size(PyObject *op) {
> if (PyArray_Check(op)) {
> return SIZE((PyArrayObject *)op);
> } else {
> return 0;
> }
> }
>
> int PyArray_CopyArray(PyArrayObject *dest, PyArrayObject *src) {
> int *dest_strides=dest->strides;
> int *dest_dimensions=dest->dimensions;
> int dest_nd=dest->nd;
> int *src_strides = src->strides;
> int *src_dimensions=src->dimensions;
> int src_nd=src->nd;
> int elsize=src->descr->elsize;
> int copies=1;
> int ret;
>
> if (src->nd > dest->nd) {
> PyErr_SetString(PyExc_ValueError, "array too large for destination");
> return -1;
> }
> if (dest->descr->type_num != src->descr->type_num) {
> PyErr_SetString(PyExc_ValueError, "can only copy from a array of the same
type.");
> return -1;
> }
>
> if (optimize_slices(&dest_strides, &dest_dimensions, &dest_nd,
&src_strides, &src_dimensions, &src_nd,
> &elsize, &copies) == -1) return -1;
>
> ret = do_sliced_copy(dest->data, dest_strides, dest_dimensions, dest_nd,
> src->data, src_strides, src_dimensions, src_nd,
> elsize, copies);
>
> if (ret != -1) { ret = PyArray_INCREF(dest); }
> return ret;
> }
>
> int PyArray_CopyObject(PyArrayObject *dest, PyObject *src_object) {
> PyArrayObject *src;
> PyObject *tmp;
> int ret, n_new, n_old;
> char *new_string;
>
> /* Special function added here to try and make arrays of strings work out.
*/
> if ((dest->descr->type_num == PyArray_CHAR) && dest->nd > 0 &&
PyString_Check(src_object)) {
> n_new = dest->dimensions[dest->nd-1];
> n_old = PyString_Size(src_object);
> if (n_new > n_old) {
> new_string = (char *)malloc(n_new*sizeof(char));
> memcpy(new_string, PyString_AS_STRING((PyStringObject *)src_object),
n_old);
> memset(new_string+n_old, ' ', n_new-n_old);
> tmp = PyString_FromStringAndSize(new_string, n_new);
> free(new_string);
> src_object = tmp;
> }
> }
> src = (PyArrayObject *)PyArray_FromObject(src_object,
> dest->descr->type_num, 0,
> dest->nd);
> if (src == NULL) return -1;
>
> ret = PyArray_CopyArray(dest, src);
> Py_DECREF(src);
> return ret;
> }
>
> /* This is the basic array allocation function. */
> PyObject *PyArray_FromDimsAndDataAndDescr(int nd, int *d, PyArray_Descr
*descr,
> char *data) {
> PyArrayObject *self;
> int i,sd;
> int *dimensions, *strides;
> int flags=CONTIGUOUS | OWN_DIMENSIONS | OWN_STRIDES;
>
> dimensions = strides = NULL;
>
> if (nd < 0) {
> PyErr_SetString(PyExc_ValueError, "number of dimensions must be >= 0");
> return NULL;
> }
>
> if (nd > 0) {
> if ((dimensions = (int *)malloc(nd*sizeof(int))) == NULL) {
> PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
> goto fail;
> }
> if ((strides = (int *)malloc(nd*sizeof(int))) == NULL) {
> PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
> goto fail;
> }
> memcpy(dimensions, d, sizeof(int)*nd);
> }
>
> sd = descr->elsize;
> for(i=nd-1;i>=0;i--) {
> if (flags & OWN_STRIDES) strides[i] = sd;
> if (dimensions[i] < 0) {
> PyErr_SetString(PyExc_ValueError, "negative dimensions are not allowed");
> goto fail;
> }
> /*
> This may waste some space, but it seems to be
> (unsuprisingly) unhealthy to allow strides that are
> longer than sd.
> */
> sd *= dimensions[i] ? dimensions[i] : 1;
> /* sd *= dimensions[i]; */
> }
>
> /* Make sure we're alligned on ints. */
> sd += sizeof(int) - sd%sizeof(int);
>
> if (data == NULL) {
> if ((data = (char *)malloc(sd)) == NULL) {
> PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
> goto fail;
> }
> flags |= OWN_DATA;
> }
>
> /*
> if((self = PyObject_NEW(PyArrayObject, &PyArray_Type)) == NULL) goto fail;
> */
>
> if ((self = (PyArrayObject
*)PyObject_CallObject(OBJECT(&PyArray_Type),
> Py_BuildValue("()")))
> == NULL) goto fail;
> if (flags & OWN_DATA) memset(data, 0, sd);
>
> self->data=data;
> self->dimensions = dimensions;
> self->strides = strides;
> self->nd=nd;
> self->descr=descr;
> self->base = (PyObject *)NULL;
> self->flags = flags;
>
> return (PyObject*)self;
>
> fail:
>
> if (flags & OWN_DATA) free(data);
> if (dimensions != NULL) free(dimensions);
> if (strides != NULL) free(strides);
> return NULL;
> }
>
> PyObject *PyArray_FromDimsAndData(int nd, int *d, int type, char *data) {
> PyArray_Descr *descr;
> if ((descr = PyArray_DescrFromType(type)) == NULL) return NULL;
> return PyArray_FromDimsAndDataAndDescr(nd, d, descr, data);
> }
>
> PyObject *PyArray_FromDims(int nd, int *d, int type) {
> PyArray_Descr *descr;
> if ((descr = PyArray_DescrFromType(type)) == NULL) return NULL;
> return PyArray_FromDimsAndDataAndDescr(nd, d, descr, NULL);
> }
>
>
> extern PyObject *PyArray_Copy(PyArrayObject *m1) {
> PyArrayObject *ret =
> (PyArrayObject *)PyArray_FromDims(m1->nd, m1->dimensions,
m1->descr->type_num);
>
> if (PyArray_CopyArray(ret, m1) == -1) return NULL;
>
> return (PyObject *)ret;
> }
>
> static void array_dealloc(PyArrayObject *self) {
> if(self->base) Py_DECREF(self->base);
>
> if (self->flags & OWN_DATA) {
> PyArray_XDECREF(self);
> free(self->data);
> }
>
> if (self->flags & OWN_DIMENSIONS && self->dimensions != NULL)
free(self->dimensions);
> if (self->flags & OWN_STRIDES && self->strides != NULL)
free(self->strides);
>
> PyMem_DEL(self);
> }
>
> /* Code to handle accessing Array objects as sequence objects */
> static int array_length(PyArrayObject *self) {
> if (self->nd != 0) {
> return self->dimensions[0];
> } else {
> return 1; /* Because a[0] works on 0d arrays. */
> }
> }
>
> static PyObject *array_item(PyArrayObject *self, int i) {
> char *item;
>
> if ((item = index2ptr(self, i)) == NULL) return NULL;
>
> if(self->nd > 1) {
> PyArrayObject *r;
> if ((r = (PyArrayObject *)PyObject_CallObject(OBJECT(&PyArray_Type),
NULL)) == NULL) return NULL;
>
> r->nd = self->nd-1;
> r->dimensions = self->dimensions+1;
> r->strides = self->strides+1;
> r->descr = self->descr;
> r->data = item;
> r->base = (PyObject *)self;
> r->flags = self->flags & CONTIGUOUS;
> Py_INCREF(self);
> return (PyObject*)r;
> } else {
> return self->descr->getitem(item);
> }
> }
>
> extern PyObject * PyArray_Item(PyObject *op, int i) {
> if (PyArray_Check(op)) {
> return array_item((PyArrayObject *)op, i);
> } else {
> PyErr_SetString(PyExc_ValueError, "Not an array object");
> return NULL;
> }
> }
>
> extern PyObject *PyArray_Return(PyArrayObject *mp) {
> PyObject *op;
>
> if (PyErr_Occurred()) {
> if (mp != NULL) Py_DECREF(mp);
> return NULL;
> }
> if (mp->nd == 0) {
> op = array_item(mp, 0);
> Py_DECREF(mp);
> return op;
> } else {
> return (PyObject *)mp;
> }
> }
>
> static PyObject * array_slice(PyArrayObject *self, int ilow, int ihigh) {
> PyArrayObject *r;
> int l;
> char *data;
>
> if (self->nd == 0) {
> PyErr_SetString(PyExc_ValueError, "can't slice a scalar");
> return NULL;
> }
>
> l=self->dimensions[0];
> if (ihigh < 0) ihigh += l;
> if (ilow < 0) ilow += l;
> if (ilow < 0) ilow = 0;
> else if (ilow > l) ilow = l;
> if (ihigh < 0) ihigh = 0;
> else if (ihigh > l) ihigh = l;
> if (ihigh < ilow) ihigh = ilow;
>
> if (ihigh != ilow) {
> data = index2ptr(self, ilow);
> if (data == NULL) return NULL;
> } else {
> data = self->data;
> }
>
> self->dimensions[0] = ihigh-ilow;
> r = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(self->nd,
self->dimensions, self->descr, data);
> self->dimensions[0] = l;
>
> if (!ISCONTIGUOUS(self)) r->flags &= ~CONTIGUOUS;
> memcpy(r->strides, self->strides, sizeof(int)*self->nd);
> r->base = (PyObject *)self;
> Py_INCREF(self);
>
> return (PyObject *)r;
> }
>
> /* These will need some serious work when I feel like it. */
> static int array_ass_item(PyArrayObject *self, int i, PyObject *v) {
> PyObject *c=NULL;
> PyArrayObject *tmp;
> char *item;
> int ret;
>
> if (v == NULL) {
> PyErr_SetString(PyExc_ValueError, "Can't delete array elements.");
> return -1;
> }
>
> if (i < 0) i = i+self->dimensions[0];
>
> if (self->nd > 1) {
> if((tmp = (PyArrayObject *)array_item(self, i)) == NULL) return -1;
> ret = PyArray_CopyObject(tmp, v);
> Py_DECREF(tmp);
> return ret;
> }
>
> if ((item = index2ptr(self, i)) == NULL) return -1;
>
> if(self->descr->type_num != PyArray_OBJECT && PyString_Check(v) &&
PyObject_Length(v) == 1) {
> char *s;
> if ((s=PyString_AsString(v)) == NULL) return -1;
> if(self->descr->type == 'c') {
> ((char*)self->data)[i]=*s;
> return 0;
> }
> if(s) c=PyInt_FromLong((long)*s);
> if(c) v=c;
> }
>
> self->descr->setitem(v, item);
>
> if(c) Py_DECREF(c);
> if(PyErr_Occurred()) return -1;
> return 0;
> }
>
>
> static int array_ass_slice(PyArrayObject *self, int ilow, int ihigh,
PyObject *v) {
> int ret;
> PyArrayObject *tmp;
>
> if (v == NULL) {
> PyErr_SetString(PyExc_ValueError, "Can't delete array elements.");
> return -1;
> }
> if ((tmp = (PyArrayObject *)array_slice(self, ilow, ihigh)) == NULL)
return -1;
> ret = PyArray_CopyObject(tmp, v);
> Py_DECREF(tmp);
>
> return ret;
> }
>
> /* -------------------------------------------------------------- */
> static int slice_GetIndices(PySliceObject *r, int length,
> int *start, int *stop, int *step)
> {
> if (r->step == Py_None) {
> *step = 1;
> } else {
> if (!PyInt_Check(r->step)) return -1;
> *step = PyInt_AsLong(r->step);
> }
> if (r->start == Py_None) {
> *start = *step < 0 ? length-1 : 0;
> } else {
> if (!PyInt_Check(r->start)) return -1;
> *start = PyInt_AsLong(r->start);
> if (*start < 0) *start += length;
> }
> if (r->stop == Py_None) {
> *stop = *step < 0 ? -1 : length;
> } else {
> if (!PyInt_Check(r->stop)) return -1;
> *stop = PyInt_AsLong(r->stop);
> if (*stop < 0) *stop += length;
> }
>
> if (*start > (length-1)) *start = length;
> if (*start < 0) *start = 0;
> if (*stop < -1) *stop = -1;
> else if (*stop > length) *stop = length;
> return 0;
> }
>
>
> static int get_slice(PyObject *op, int max, int *np, int *sp) {
> int start, stop, step;
>
> if (PySlice_Check(op)) {
> if (slice_GetIndices((PySliceObject *)op, max,
> &start, &stop, &step) == -1) return -1;
>
> if (step != 0) {
> if (step < 0) *np = (stop-start+1+step)/step;
> else *np = (stop-start-1+step)/step;
> } else {
> if (stop == start) {
> *np = 0; step = 1;
> }
> else return -1;
> }
> if (*np < 0) *np = 0;
> *sp = step;
> return start;
> }
> return -1;
> }
>
> #define PseudoIndex -1
> #define RubberIndex -2
> #define SingleIndex -3
>
> static int parse_subindex(PyObject *op, int *step_size, int *n_steps, int
max) {
> int i, tmp;
>
> if (op == Py_None) {
> *n_steps = PseudoIndex;
> return 0;
> }
>
> if (op == Py_Ellipsis) {
> *n_steps = RubberIndex;
> return 0;
> }
>
> if (PySlice_Check(op)) {
> if ((i = get_slice(op, max, n_steps, step_size)) >= 0) {
> return i;
> } else {
> PyErr_SetString(PyExc_IndexError, "invalid slice");
> return -1;
> }
> }
>
> if (PyInt_Check(op)) {
> *n_steps=SingleIndex;
> *step_size=0;
> tmp = PyInt_AsLong(op);
> if (tmp < 0) tmp += max;
> if (tmp >= max || tmp < 0) {
> PyErr_SetString(PyExc_IndexError, "invalid index");
> return -1;
> }
> return tmp;
> }
> PyErr_SetString(PyExc_IndexError,
> "each subindex must be either a slice, an integer, Ellipsis, or NewAxis");
> return -1;
> }
>
> static int parse_index(PyArrayObject *self, PyObject *op,
> int *dimensions, int *strides, int *offset_ptr) {
> int i, j, n;
> int nd_old, nd_new, start, offset, n_add, n_pseudo;
> int step_size, n_steps;
> PyObject *op1;
> int is_slice;
>
> if (PySlice_Check(op) || op == Py_Ellipsis) {
> n = 1;
> op1 = op;
> Py_INCREF(op); /* this relies on the fact that n==1 for loop below */
> is_slice = 1;
> } else {
> if (!PySequence_Check(op)) {
> PyErr_SetString(PyExc_IndexError, "index must be either an int or a
sequence");
> return -1;
> }
> n = PySequence_Length(op);
> is_slice = 0;
> }
>
> nd_old = nd_new = 0;
>
> offset = 0;
> for(i=0; i<n; i++) {
> if (!is_slice) {
> if (!(op1=PySequence_GetItem(op, i))) {
> PyErr_SetString(PyExc_IndexError, "invalid index");
> return -1;
> }
> }
>
> start = parse_subindex(op1, &step_size, &n_steps,
> nd_old < self->nd ? self->dimensions[nd_old] : 0);
> Py_DECREF(op1);
> if (start == -1) break;
>
> if (n_steps == PseudoIndex) {
> dimensions[nd_new] = 1; strides[nd_new] = 0; nd_new++;
> } else {
> if (n_steps == RubberIndex) {
> for(j=i+1, n_pseudo=0; j<n; j++) {
> op1 = PySequence_GetItem(op, j);
> if (op1 == Py_None) n_pseudo++;
> Py_DECREF(op1);
> }
> n_add = self->nd-(n-i-n_pseudo-1+nd_old);
> if (n_add < 0) {
> PyErr_SetString(PyExc_IndexError, "too many indices");
> return -1;
> }
> for(j=0; j<n_add; j++) {
> dimensions[nd_new] = self->dimensions[nd_old];
> strides[nd_new] = self->strides[nd_old];
> nd_new++; nd_old++;
> }
> } else {
> if (nd_old >= self->nd) {
> PyErr_SetString(PyExc_IndexError, "too many indices");
> return -1;
> }
> offset += self->strides[nd_old]*start;
> nd_old++;
> if (n_steps != SingleIndex) {
> dimensions[nd_new] = n_steps;
> strides[nd_new] = step_size*self->strides[nd_old-1];
> nd_new++;
> }
> }
> }
> }
> if (i < n) return -1;
> n_add = self->nd-nd_old;
> for(j=0; j<n_add; j++) {
> dimensions[nd_new] = self->dimensions[nd_old];
> strides[nd_new] = self->strides[nd_old];
> nd_new++; nd_old++;
> }
> *offset_ptr = offset;
> return nd_new;
> }
>
> static PyObject *array_subscript(PyArrayObject *self, PyObject *op) {
> int dimensions[MAX_DIMS], strides[MAX_DIMS];
> int nd, offset, i, elsize;
> PyArrayObject *other;
>
> if (PyInt_Check(op)) {
> i = PyInt_AsLong(op);
> if (i < 0 && self->nd > 0) i = i+self->dimensions[0];
> return array_item(self, i);
> }
>
> if ((nd = parse_index(self, op, dimensions, strides, &offset)) == -1) {
> return NULL;
> }
>
> if ((other = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(nd,
> dimensions,
> self->descr,
> self->data+offset)) == NULL) {
> return NULL;
> }
> memcpy(other->strides, strides, sizeof(int)*other->nd);
> other->base = (PyObject *)self;
> Py_INCREF(self);
>
> elsize=other->descr->elsize;
> for (i=other->nd-1; i>=0; i--) {
> if (other->strides[i] == elsize) {
> elsize *= other->dimensions[i];
> } else {
> break;
> }
> }
> if (i >= 0) other->flags &= ~CONTIGUOUS;
>
> return (PyObject *)other;
> }
>
> static PyObject *array_subscript_nice(PyArrayObject *self, PyObject *op) {
> PyObject *ret;
>
> if ((ret = array_subscript(self, op)) == NULL) return NULL;
> if (PyArray_Check(ret)) return PyArray_Return((PyArrayObject *)ret);
> else return ret;
> }
>
> /* Another assignment hacked by using CopyObject. */
>
> static int array_ass_sub(PyArrayObject *self, PyObject *index, PyObject
*op) {
> int ret;
> PyArrayObject *tmp;
>
> if (op == NULL) {
> PyErr_SetString(PyExc_ValueError, "Can't delete array elements.");
> return -1;
> }
>
> if (PyInt_Check(index)) return array_ass_item(self, PyInt_AsLong(index),
op);
>
> if ((tmp = (PyArrayObject *)array_subscript(self, index)) == NULL)
return -1;
> ret = PyArray_CopyObject(tmp, op);
> Py_DECREF(tmp);
>
> return ret;
> }
>
> static PyMappingMethods array_as_mapping = {
> (inquiry)array_length, /*mp_length*/
> (binaryfunc)array_subscript_nice, /*mp_subscript*/
> (objobjargproc)array_ass_sub, /*mp_ass_subscript*/
> };
>
> /* -------------------------------------------------------------- */
>
> typedef struct {
> PyUFuncObject *add, *subtract, *multiply, *divide, *remainder, *power,
*negative, *absolute;
> PyUFuncObject *invert, *left_shift, *right_shift, *bitwise_and,
*bitwise_xor,
> *bitwise_or;
> } NumericOps;
>
> static NumericOps n_ops;
>
> #define GET(op) n_ops.op = (PyUFuncObject *)PyDict_GetItemString(dict,
#op)
>
> int PyArray_SetNumericOps(PyObject *dict) {
> GET(add);
> GET(subtract);
> GET(multiply);
> GET(divide);
> GET(remainder);
> GET(power);
> GET(negative);
> GET(absolute);
> GET(invert);
> GET(left_shift);
> GET(right_shift);
> GET(bitwise_and);
> GET(bitwise_or);
> GET(bitwise_xor);
> return 0;
> }
>
> static int array_coerce(PyArrayObject **pm, PyObject **pw) {
> PyObject *new_op;
> if ((new_op = PyArray_FromObject(*pw, PyArray_NOTYPE, 0, 0)) == NULL)
return -1;
> Py_INCREF(*pm);
> *pw = new_op;
> return 0;
>
> /**Eliminate coercion at this stage. Handled by umath now
> PyObject *new_op;
> char t1, t2;
>
> t1 = (*pm)->descr->type_num;
> t2 = PyArray_ObjectType(*pw, 0);
>
> if (PyArray_CanCastSafely(t2, t1))
> if ((new_op = PyArray_FromObject(*pw, t1, 0, 0)) == NULL) return -1;
> Py_INCREF(*pm);
> *pw = new_op;
> } else {
> if (PyArray_CanCastSafely(t1, t2)) {
> *pm = (PyArrayObject *)PyArray_Cast(*pm, t2);
> if ((new_op = PyArray_FromObject(*pw, t2, 0, 0)) == NULL) return -1;
> *pw = new_op;
> } else {
> PyErr_SetString(PyExc_TypeError, "cannot perform this operation on these
types");
> return -1;
> }
> }
> return 0;
> ***/
> }
>
> static PyObject *PyUFunc_BinaryFunction(PyUFuncObject *s, PyArrayObject
*mp1, PyObject *mp2) {
> PyObject *arglist;
> PyArrayObject *mps[3];
>
> arglist = Py_BuildValue("(OO)", mp1, mp2);
>
> mps[0] = mps[1] = mps[2] = NULL;
> if (PyUFunc_GenericFunction(s, arglist, mps) == -1) {
> Py_XDECREF(mps[0]); Py_XDECREF(mps[1]); Py_XDECREF(mps[2]);
> return NULL;
> }
>
> Py_DECREF(mps[0]); Py_DECREF(mps[1]);
> Py_DECREF(arglist);
> return PyArray_Return(mps[2]);
> }
>
> static PyObject *PyUFunc_UnaryFunction(PyUFuncObject *s, PyArrayObject
*mp1) {
> PyObject *arglist;
> PyArrayObject *mps[3];
>
> arglist = Py_BuildValue("(O)", mp1);
>
> mps[0] = mps[1] = NULL;
> if (PyUFunc_GenericFunction(s, arglist, mps) == -1) {
> Py_XDECREF(mps[0]); Py_XDECREF(mps[1]);
> return NULL;
> }
>
> Py_DECREF(mps[0]);
> Py_DECREF(arglist);
> return PyArray_Return(mps[1]);
> }
>
> static PyObject *array_add(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.add, m1, m2);
> }
> static PyObject *array_subtract(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.subtract, m1, m2);
> }
> static PyObject *array_multiply(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.multiply, m1, m2);
> }
> static PyObject *array_divide(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.divide, m1, m2);
> }
> static PyObject *array_remainder(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.remainder, m1, m2);
> }
> static PyObject *array_power(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.power, m1, m2);
> }
> static PyObject *array_negative(PyArrayObject *m1)
> return PyUFunc_UnaryFunction(n_ops.negative, m1);
> }
> static PyObject *array_absolute(PyArrayObject *m1)
> return PyUFunc_UnaryFunction(n_ops.absolute, m1);
> }
> static PyObject *array_invert(PyArrayObject *m1)
> return PyUFunc_UnaryFunction(n_ops.invert, m1);
> }
> static PyObject *array_left_shift(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.left_shift, m1, m2);
> }
> static PyObject *array_right_shift(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.right_shift, m1, m2);
> }
> static PyObject *array_bitwise_and(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.bitwise_and, m1, m2);
> }
> static PyObject *array_bitwise_or(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.bitwise_or, m1, m2);
> }
> static PyObject *array_bitwise_xor(PyArrayObject *m1, PyObject *m2) {
> return PyUFunc_BinaryFunction(n_ops.bitwise_xor, m1, m2);
> }
>
>
> static int array_nonzero(PyArrayObject *mp) {
> char *zero;
> PyArrayObject *self;
> char *data;
> int i, s, elsize;
>
> self = PyArray_CONTIGUOUS(mp);
> zero = self->descr->zero;
>
> s = SIZE(self);
> elsize = self->descr->elsize;
> data = self->data;
> for(i=0; i<s; i++, data+=elsize) {
> if (memcmp(zero, data, elsize) != 0) break;
> }
>
> Py_DECREF(self);
>
> return i!=s;
> }
>
> static PyObject *array_divmod(PyArrayObject *op1, PyObject *op2) {
> PyObject *divp, *modp;
>
> divp = array_divide(op1, op2);
> if (divp == NULL) return NULL;
> modp = array_remainder(op1, op2);
> if (modp == NULL) {
> Py_DECREF(divp);
> return NULL;
> }
> return Py_BuildValue("OO", divp, modp);
> }
>
>
> static PyNumberMethods array_as_number = {
> (binaryfunc)array_add, /*nb_add*/
> (binaryfunc)array_subtract, /*nb_subtract*/
> (binaryfunc)array_multiply, /*nb_multiply*/
> (binaryfunc)array_divide, /*nb_divide*/
> (binaryfunc)array_remainder, /*nb_remainder*/
> (binaryfunc)array_divmod, /*nb_divmod*/
> (ternaryfunc)array_power, /*nb_power*/
> (unaryfunc)array_negative,
> (unaryfunc)PyArray_Copy, /*nb_pos*/
> (unaryfunc)array_absolute, /* (unaryfunc)array_abs, */
> (inquiry)array_nonzero, /*nb_nonzero*/
> (unaryfunc)array_invert, /*nb_invert*/
> (binaryfunc)array_left_shift, /*nb_lshift*/
> (binaryfunc)array_right_shift, /*nb_rshift*/
> (binaryfunc)array_bitwise_and, /*nb_and*/
> (binaryfunc)array_bitwise_xor, /*nb_xor*/
> (binaryfunc)array_bitwise_or, /*nb_or*/
> (coercion)array_coerce, /*nb_coerce*/
> 0, /*nb_int*/
> 0, /*nb_long*/
> 0, /*nb_float*/
> 0, /*nb_oct*/
> 0, /*nb_hex*/
> };
>
> static PySequenceMethods array_as_sequence = {
> (inquiry)array_length, /*sq_length*/
> (binaryfunc)array_add, /*nb_add, concat is numeric add*/
> (intargfunc)NULL, /*nb_multiply, repeat is numeric multiply*/
> (intargfunc)array_item, /*sq_item*/
> (intintargfunc)array_slice, /*sq_slice*/
> (intobjargproc)array_ass_item, /*sq_ass_item*/
> (intintobjargproc)array_ass_slice, /*sq_ass_slice*/
> };
>
> /* -------------------------------------------------------------- */
>
> /*This ugly function quickly prints out an array. It's likely to
disappear */
> /*in the near future. */
> static int dump_data(char **string, int *n, int *max_n,
> char *data, int nd, int *dimensions,
> int *strides, PyArray_Descr *descr) {
> PyObject *op, *sp;
> char *ostring;
> int i, N;
>
> #define CHECK_MEMORY if (*n >= *max_n-16) { *max_n *= 2; *string = (char
*)realloc(*string, *max_n); }
>
> if (nd == 0) {
>
> if ((op = descr->getitem(data)) == NULL) return -1;
> sp = PyObject_Repr(op);
> if (sp == NULL) {Py_DECREF(op); return -1;}
> ostring = PyString_AsString(sp);
> N = PyString_Size(sp)*sizeof(char);
> *n += N;
> CHECK_MEMORY
> memcpy(*string+(*n-N), ostring, N);
> Py_DECREF(sp);
> Py_DECREF(op);
> return 0;
> } else {
> if (nd == 1 && descr->type_num == PyArray_CHAR) {
> N = (dimensions[0]+2)*sizeof(char);
> *n += N;
> CHECK_MEMORY
> (*string)[*n-N] = '"';
> memcpy(*string+(*n-N+1), data, N-2);
> (*string)[*n-1] = '"';
> return 0;
> } else {
> CHECK_MEMORY
> (*string)[*n] = '[';
> *n += 1;
> for(i=0; i<dimensions[0]; i++) {
> if (dump_data(string, n, max_n, data+(*strides)*i, nd-1, dimensions+1,
strides+1, descr) < 0) return -1;
> CHECK_MEMORY
> if (i<dimensions[0]-1) {
> (*string)[*n] = ',';
> (*string)[*n+1] = ' ';
> *n += 2;
> }
> }
> CHECK_MEMORY
> (*string)[*n] = ']'; *n += 1;
> return 0;
> }
> }
>
> #undef CHECK_MEMORY
> }
>
> static PyObject * array_repr_builtin(PyArrayObject *self) {
> PyObject *ret;
> char *string;
> int n, max_n;
>
> max_n = NBYTES(self)*4*sizeof(char) + 7;
>
> if ((string = (char *)malloc(max_n)) == NULL) {
> PyErr_SetString(PyExc_MemoryError, "out of memory");
> return NULL;
> }
>
> n = 6;
> sprintf(string, "array(");
>
> if (dump_data(&string, &n, &max_n, self->data, self->nd, self->dimensions,
> self->strides, self->descr) < 0) { free(string); return NULL; }
>
> sprintf(string+n, ", '%c')", self->descr->type);
>
> ret = PyString_FromStringAndSize(string, n+6);
> free(string);
> return ret;
> }
>
> static PyObject *PyArray_StrFunction=NULL;
> static PyObject *PyArray_ReprFunction=NULL;
>
> extern void PyArray_SetStringFunction(PyObject *op, int repr) {
> if (repr) {
> Py_XDECREF(PyArray_ReprFunction); /* Dispose of previous callback */
> Py_XINCREF(op); /* Add a reference to new callback */
> PyArray_ReprFunction = op; /* Remember new callback */
> } else {
> Py_XDECREF(PyArray_StrFunction); /* Dispose of previous callback */
> Py_XINCREF(op); /* Add a reference to new callback */
> PyArray_StrFunction = op; /* Remember new callback */
> }
> }
>
> static PyObject *array_repr(PyArrayObject *self) {
> PyObject *s, *arglist;
>
> if (PyArray_ReprFunction == NULL) {
> s = array_repr_builtin(self);
> } else {
> arglist = Py_BuildValue("(O)", self);
> s = PyEval_CallObject(PyArray_ReprFunction, arglist);
> Py_DECREF(arglist);
> }
> return s;
> }
>
> static PyObject *array_str(PyArrayObject *self) {
> PyObject *s, *arglist;
>
> if (PyArray_StrFunction == NULL) {
> s = array_repr(self);
> } else {
> arglist = Py_BuildValue("(O)", self);
> s = PyEval_CallObject(PyArray_StrFunction, arglist);
> Py_DECREF(arglist);
> }
> return s;
> }
>
> /* No longer implemented here. Use the default in object.c
> static int array_print(PyArrayObject *self, FILE *fp, int flags) {
> PyObject *s;
> int r=-1, n;
>
> s = array_str(self);
>
> if (s == NULL || !PyString_Check(s)) {
> Py_XDECREF(s);
> return -1;
> }
> n = PyString_Size(s);
> r = fwrite(PyString_AsString(s), sizeof(char), n, fp);
> Py_DECREF(s);
> return r==n ? 0 : -1;
> }
> */
>
> static void byte_swap_vector(void *p, int n, int size) {
> char *a, *b, c;
>
> switch(size) {
> case 2:
> for (a = (char*)p ; n > 0; n--, a += 1) {
> b = a + 1;
> c = *a; *a++ = *b; *b = c;
> }
> break;
> case 4:
> for (a = (char*)p ; n > 0; n--, a += 2) {
> b = a + 3;
> c = *a; *a++ = *b; *b-- = c;
> c = *a; *a++ = *b; *b = c;
> }
> break;
> case 8:
> for (a = (char*)p ; n > 0; n--, a += 4) {
> b = a + 7;
> c = *a; *a++ = *b; *b-- = c;
> c = *a; *a++ = *b; *b-- = c;
> c = *a; *a++ = *b; *b-- = c;
> c = *a; *a++ = *b; *b = c;
> }
> break;
> default:
> break;
> }
> }
>
> static char doc_byteswapped[] = "m.byteswapped(). Swaps the bytes in the
contents of array m. Useful for reading data written by a machine with a
different byte order. Returns a new array with the byteswapped data.";
>
> static PyObject * array_byteswap(PyArrayObject *self, PyObject *args) {
> PyArrayObject *ret;
>
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> if ((ret = (PyArrayObject *)PyArray_Copy(self)) == NULL) return NULL;
>
> if (self->descr->type_num < PyArray_CFLOAT) {
> byte_swap_vector(ret->data, SIZE(self), self->descr->elsize);
> } else {
> byte_swap_vector(ret->data, SIZE(self)*2, self->descr->elsize/2);
> }
>
> return (PyObject *)ret;
> }
>
> static char doc_tostring[] = "m.tostring(). Copy the data portion of the
array to a python string and return that string.";
>
> static PyObject *array_tostring(PyArrayObject *self, PyObject *args) {
> PyObject *so;
> PyArrayObject *mo;
>
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> if ((mo = PyArray_CONTIGUOUS(self)) ==NULL) return NULL;
>
> so = PyString_FromStringAndSize(mo->data, NBYTES(mo));
>
> Py_DECREF(mo);
>
> return so;
> }
>
>
> static PyObject *PyArray_ToList(PyObject *self) {
> PyObject *lp;
> PyObject *v;
> int sz, i;
>
> if (!PyArray_Check(self)) return self;
>
> if (((PyArrayObject *)self)->nd == 0) return array_item((PyArrayObject
*)self, 0);
>
> sz = ((PyArrayObject *)self)->dimensions[0];
> lp = PyList_New(sz);
>
> for (i=0; i<sz; i++) {
> v=array_item((PyArrayObject *)self, i);
> PyList_SetItem(lp, i, PyArray_ToList(v));
> if (((PyArrayObject *)self)->nd>1){
> Py_DECREF(v);
> }
> }
>
> return lp;
> }
>
>
> static char doc_tolist[] = "m.tolist(). Copy the data portion of the
array to a hierarchical python list and return that list.";
>
> static PyObject *array_tolist(PyArrayObject *self, PyObject *args) {
> if (!PyArg_ParseTuple(args, "")) return NULL;
> if (self->nd <= 0) {
> PyErr_SetString(PyExc_ValueError, "Can't convert a 0d array to a list");
> return NULL;
> }
>
> return PyArray_ToList((PyObject *)self);
> }
>
> static char doc_cast[] = "m.cast(t). Cast array m to type t. t can be
either a string representing a typecode, or a python type object of type
int, float, or complex.";
>
> PyObject *array_cast(PyArrayObject *self, PyObject *args) {
> PyObject *op;
> char typecode;
>
> if (!PyArg_ParseTuple(args, "O", &op)) return NULL;
>
> if (PyString_Check(op) && (PyString_Size(op) == 1)) {
> return PyArray_Cast(self, PyString_AS_STRING((PyStringObject *)op)[0]);
> }
>
> if (PyType_Check(op)) {
> typecode = 'O';
> if (((PyTypeObject *)op) == &PyInt_Type) {
> typecode = PyArray_LONG;
> }
> if (((PyTypeObject *)op) == &PyFloat_Type) {
> typecode = PyArray_DOUBLE;
> }
> if (((PyTypeObject *)op) == &PyComplex_Type) {
> typecode = PyArray_CDOUBLE;
> }
> return PyArray_Cast(self, typecode);
> }
> PyErr_SetString(PyExc_ValueError,
> "type must be either a 1-length string, or a python type object");
> return NULL;
> }
>
> static char doc_typecode[] = "m.typecode(). Return the single character
code indicating the type of the elements of m.";
>
> PyObject *array_typecode(PyArrayObject *self, PyObject *args) {
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> return PyString_FromStringAndSize(&(self->descr->type), 1);
> }
>
> static char doc_itemsize[] = "m.itemsize(). Return the size in bytes of a
single element of the array m.";
>
> PyObject *array_itemsize(PyArrayObject *self, PyObject *args) {
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> return PyInt_FromLong(self->descr->elsize);
> }
>
> static char doc_contiguous[] = "m.contiguous(). Return true if the memory
used by the array m is contiguous. A non-contiguous array can be converted
to a contiguous one by m.copy().";
>
> PyObject *array_contiguous(PyArrayObject *self, PyObject *args) {
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> return PyInt_FromLong(ISCONTIGUOUS(self));
> }
>
> static char doc_copy[] = "m.__copy__().";
>
> PyObject *array_copy(PyArrayObject *self, PyObject *args) {
> if (!PyArg_ParseTuple(args, "")) return NULL;
>
> return PyArray_Copy(self);
> }
>
> int array_compare(PyArrayObject *self, PyObject *other) {
> PyErr_SetString(PyExc_TypeError,
> "Comparison of multiarray objects is not implemented.");
> return -1;
> }
>
> static char doc_init[] = "Array(sequence, copy=1, typecode=None) will
return a new Array ExtensionClass instance formed from the given
(potentially nested) sequence and of type given by typecode. If no typecode
is given, then the type will be determined as the minimum type required to
hold the objects in sequence.If copy is zero and the the sequence in of type
array, then the sequence is not copied, i.e., the new array and the sequence
share the same data. Note that keywords cannot be used and that the order of
the last two arguments is reversed for that in the function array.";
>
> static PyObject *array_init(PyArrayObject *self, PyObject *args) {
> PyObject *op=Py_None, *tpo=Py_None;
> PyArrayObject *ret;
> int sd;
> PyArray_Descr *descr;
> char *data, *tp=NULL;
> int type, copy=1;
> int *dimensions, *strides;
> int flags=0;
> int nd = 0;
>
> if(!PyArg_ParseTuple(args, "|OiO", &op, ©, &tpo)) return NULL;
>
> /* Have data */
> if (op != Py_None) {
> if (tpo == Py_None) {
> type = PyArray_NOTYPE;
> } else {
> tp = PyString_AsString(tpo);
> if (tp[0] == 0) type = PyArray_NOTYPE;
> else type = tp[0];
> }
>
> if (copy) {
> if ((ret = (PyArrayObject *)PyArray_CopyFromObject(op, type, 0, 0)) ==
NULL)
> return NULL;
> } else {
> if ((ret = (PyArrayObject *)PyArray_FromObject(op, type, 0, 0)) == NULL)
> return NULL;
> }
> /* Copy ret to self. */
> self->data = ret->data;
> self->dimensions = ret->dimensions;
> self->strides = ret->strides;
> self->nd = ret->nd;
> self->descr= ret->descr;
>
> self->base = ret->base;
> self->flags = ret->flags;
>
> ret->flags = 0;
> Py_INCREF(Py_None);
> return Py_None;
> }
>
> /* Don't have any data */
> flags=CONTIGUOUS | OWN_DIMENSIONS | OWN_STRIDES;
> data = NULL;
> dimensions = strides = NULL;
> type = PyArray_DOUBLE;
>
> if ((descr = PyArray_DescrFromType(type)) == NULL) return NULL;
> sd = descr->elsize;
> /* Make sure we're alligned on ints. */
> sd += sizeof(int) - sd%sizeof(int);
>
> if ((data = (char *)malloc(sd)) == NULL) {
> PyErr_SetString(PyExc_MemoryError, "can't allocate memory for array");
> goto fail;
> }
> flags |= OWN_DATA;
> if (flags & OWN_DATA) memset(data, 0, sd);
>
> self->data=data;
> self->dimensions = dimensions;
> self->strides = strides;
> self->nd=nd;
> self->descr=descr;
> self->base = (PyObject *)NULL;
> self->flags = flags;
> Py_INCREF(Py_None);
> return Py_None;
>
> fail:
> if (flags & OWN_DATA) free(data);
> return NULL;
> }
>
>
>
> static struct PyMethodDef array_methods[] = {
> {"tostring", (PyCFunction)array_tostring, 1, doc_tostring},
> {"tolist", (PyCFunction)array_tolist, 1, doc_tolist},
> {"byteswapped", (PyCFunction)array_byteswap, 1, doc_byteswapped},
> {"astype", (PyCFunction)array_cast, 1, doc_cast},
> {"typecode", (PyCFunction)array_typecode, 1, doc_typecode},
> {"itemsize", (PyCFunction)array_itemsize, 1, doc_itemsize},
> {"iscontiguous", (PyCFunction)array_contiguous, 1, doc_contiguous},
> {"__copy__", (PyCFunction)array_copy, 1, doc_copy},
> {"__init__", (PyCFunction)array_init, METH_VARARGS, doc_init},
> {NULL, NULL} /* sentinel */
> };
>
> /* ---------- */
> static PyObject *array_getattr(PyArrayObject *self, char *name) {
> PyArrayObject *ret;
>
> if (strcmp(name, "shape") == 0) {
> PyObject *s, *o;
> int i;
>
> if ((s=PyTuple_New(self->nd)) == NULL) return NULL;
>
> for(i=self->nd; --i >= 0;) {
> if ((o=PyInt_FromLong(self->dimensions[i])) == NULL) return NULL;
> if (PyTuple_SetItem(s,i,o) == -1) return NULL;
> }
> return s;
> }
>
> if (strcmp(name, "real") == 0) {
> if (self->descr->type_num == PyArray_CFLOAT ||
> self->descr->type_num == PyArray_CDOUBLE) {
> ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
> self->dimensions,
> self->descr->type_num-2,
> self->data);
> if (ret == NULL) return NULL;
> memcpy(ret->strides, self->strides, sizeof(int)*ret->nd);
> ret->flags &= ~CONTIGUOUS;
> Py_INCREF(self);
> ret->base = (PyObject *)self;
> return (PyObject *)ret;
> } else {
> ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
> self->dimensions,
> self->descr->type_num,
> self->data);
> if (ret == NULL) return NULL;
> Py_INCREF(self);
> ret->base = (PyObject *)self;
> }
> }
>
> if ((strcmp(name, "imaginary") == 0) || (strcmp(name, "imag") == 0))
> if (self->descr->type_num == PyArray_CFLOAT ||
> self->descr->type_num == PyArray_CDOUBLE) {
> ret = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
> self->dimensions,
> self->descr->type_num-2,
> self->data+self->descr->elsize/2);
> if (ret == NULL) return NULL;
> memcpy(ret->strides, self->strides, sizeof(int)*ret->nd);
> ret->flags &= ~CONTIGUOUS;
> Py_INCREF(self);
> ret->base = (PyObject *)self;
> return (PyObject *)ret;
> } else {
> PyErr_SetString(PyExc_ValueError, "No imaginary part to real array");
> return NULL;
> }
> }
>
> if (strcmp(name, "flat") == 0) {
> int n;
> n = SIZE(self);
> if (!ISCONTIGUOUS(self)) {
> PyErr_SetString(PyExc_ValueError, "flattened indexing only available for
contiguous array");
> return NULL;
> }
> ret = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(1, &n,
> self->descr,
> self->data);
> if (ret == NULL) return NULL;
> Py_INCREF(self);
> ret->base = (PyObject *)self;
> return (PyObject *)ret;
> }
>
> return Py_FindMethod(array_methods, (PyObject *)self, name);
> }
>
> static int array_setattr(PyArrayObject *self, char *name, PyObject *op) {
> PyArrayObject *ap;
> int ret;
>
> if (strcmp(name, "shape") == 0) {
> /* This can be made more efficient by copying code from array_reshape if
needed */
> if ((ap = (PyArrayObject *)PyArray_Reshape(self, op)) == NULL) return -1;
> if(self->flags & OWN_DIMENSIONS) free(self->dimensions);
> self->dimensions = ap->dimensions;
> if(self->flags & OWN_STRIDES) free(self->strides);
> self->strides = ap->strides;
> self->nd = ap->nd;
> self->flags &= ~(OWN_DIMENSIONS | OWN_STRIDES);
> self->flags |= ap->flags & (OWN_DIMENSIONS | OWN_STRIDES);
> ap->flags &= ~(OWN_DIMENSIONS | OWN_STRIDES);
> Py_DECREF(ap);
> return 0;
> }
>
> if (strcmp(name, "real") == 0) {
> if (self->descr->type_num == PyArray_CFLOAT ||
> self->descr->type_num == PyArray_CDOUBLE) {
> ap = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
> self->dimensions,
> self->descr->type_num-2,
> self->data);
> if (ap == NULL) return -1;
> memcpy(ap->strides, self->strides, sizeof(int)*ap->nd);
> ap->flags &= ~CONTIGUOUS;
> ret = PyArray_CopyObject(ap, op);
> Py_DECREF(ap);
> return ret;
> } else {
> return PyArray_CopyObject(self, op);
> }
> }
>
> if ((strcmp(name, "imaginary") == 0) || (strcmp(name, "imag") == 0))
> if (self->descr->type_num == PyArray_CFLOAT ||
> self->descr->type_num == PyArray_CDOUBLE) {
> ap = (PyArrayObject *)PyArray_FromDimsAndData(self->nd,
> self->dimensions,
> self->descr->type_num-2,
> self->data+self->descr->elsize/2);
> if (ap == NULL) return -1;
> memcpy(ap->strides, self->strides, sizeof(int)*ap->nd);
> ap->flags &= ~CONTIGUOUS;
> ret = PyArray_CopyObject(ap, op);
> Py_DECREF(ap);
> return ret;
> } else {
> PyErr_SetString(PyExc_ValueError, "No imaginary part to real array");
> return -1;
> }
> }
>
> return PyDict_SetItemString(INSTANCE_DICT(self), name, op);
> }
>
> static char Arraytype__doc__[] =
> "A array object represents a multidimensional, homogeneous array of basic
values. It has the folowing data members, m.shape (the size of each
dimension in the array), m.itemsize (the size (in bytes) of each element of
the array), and m.typecode (a character representing the type of the
matrices elements). Matrices are sequence, mapping and numeric objects.
Sequence indexing is similar to lists, with single indices returning a
reference that points to the old matrices data, and slices returning by
copy. A array is also allowed to be indexed by a sequence of sequences or
ints. Each member of the sequence indexes the corresponding dimension of
the array. If it is a sequence, it returns the elments along that dimension
listed in the sequence, if a single number, it returns just that element
(reducing the dimensionality of the returned array by 1. Numeric operations
operate on matrices in an element-wise fashion.";
>
>
> /* PyTypeObject PyArray_Type = { */
> PyExtensionClass PyArray_Type = {
> PyObject_HEAD_INIT(0)
> 0, /*ob_size*/
> "array", /*tp_name*/
> sizeof(PyArrayObject), /*tp_basicsize*/
> 0, /*tp_itemsize*/
> /* methods */
> (destructor)array_dealloc, /*tp_dealloc*/
> (printfunc)NULL, /*tp_print*/
> (getattrfunc)array_getattr, /*tp_getattr*/
> (setattrfunc)array_setattr, /*tp_setattr*/
> (cmpfunc)array_compare, /*tp_compare*/
> (reprfunc)array_repr, /*tp_repr*/
> &array_as_number, /*tp_as_number*/
> &array_as_sequence, /*tp_as_sequence*/
> &array_as_mapping, /*tp_as_mapping*/
> (hashfunc)0, /*tp_hash*/
> (ternaryfunc)0, /*tp_call*/
> (reprfunc)array_str, /*tp_str*/
>
> /* Space for future expansion */
> 0L,0L,0L,0L,
> Arraytype__doc__, /* Documentation string */
> METHOD_CHAIN(array_methods) /* Required by ExtensionClass */
> };
>
>
> /* The rest of this code is to build the right kind of array from a python
*/
> /* object. */
>
> static int discover_depth(PyObject *s, int max, int stop_at_string) {
> int d=0;
> PyObject *e;
>
> if(max < 1) return -1;
>
> /*if(! PySequence_Check(s) || PyInstance_Check(s)) return 0;
> /*
> Relax this to allow class instances. Dangerous? */
> if(! PySequence_Check(s) || PySequence_Length(s) < 0) {
> PyErr_Clear(); return 0;
> }
> if(PyString_Check(s)) return stop_at_string ? 0:1;
> if (PySequence_Length(s) == 0) return 1;
>
> if ((e=PySequence_GetItem(s,0)) == NULL) return -1;
> if(e!=s)
> {
> d=discover_depth(e,max-1, stop_at_string);
> if(d >= 0) d++;
> }
> Py_DECREF(e);
> return d;
> }
>
> static int discover_dimensions(PyObject *s, int nd, int *d) {
> PyObject *e;
> int r, n, i, n_lower;
>
> n=PyObject_Length(s);
> *d = n;
> if(*d < 0) return -1;
> if(nd <= 1) return 0;
> n_lower = 0;
> for(i=0; i<n; i++) {
> if ((e=PySequence_GetItem(s,i)) == NULL) return -1;
> r=discover_dimensions(e,nd-1,d+1);
> Py_DECREF(e);
>
> if (r == -1) return -1;
> if (d[1] > n_lower) n_lower = d[1];
> }
> d[1] = n_lower;
>
> return 0;
> }
>
> #ifndef max
> #define max(a,b) (a)>(b)?(a):(b);
> #endif
>
> int PyArray_ObjectType(PyObject *op, int minimum_type) {
> int l;
> PyObject *ip;
>
> if (minimum_type == -1) return -1;
>
> if (PyArray_Check(op)) return max((int)((PyArrayObject
*)op)->descr->type_num, minimum_type);
>
> if (PyInstance_Check(op)) {
> if (PyObject_HasAttrString(op, "__array__")) {
> PyObject *ap, *arglist;
>
> arglist = Py_BuildValue("()");
> ap = PyObject_GetAttrString(op, "__array__");
> ip = PyEval_CallObject(ap, arglist);
> Py_DECREF(arglist);
>
> return max(minimum_type, (int)((PyArrayObject *)ip)->descr->type_num);
> } else {
> if (PySequence_Length(op)<0)
> PyErr_Clear();
> return (int)PyArray_OBJECT;
> }
> }
>
> if (PyString_Check(op)) {
> return max(minimum_type, (int)PyArray_CHAR);
> }
>
> if (PySequence_Check(op)) {
> l = PyObject_Length(op);
> if (l == 0 && minimum_type == 0) minimum_type = PyArray_LONG;
> while (--l >= 0) {
> ip = PySequence_GetItem(op, l);
> minimum_type = PyArray_ObjectType(ip, minimum_type);
> Py_DECREF(ip);
> }
> return minimum_type;
> }
>
> if (PyInt_Check(op)) {
> return max(minimum_type, (int)PyArray_LONG);
> } else {
> if (PyFloat_Check(op)) {
> return max(minimum_type, (int)PyArray_DOUBLE);
> } else {
> if (PyComplex_Check(op)) {
> return max(minimum_type, (int)PyArray_CDOUBLE);
> } else {
> return (int)PyArray_OBJECT;
> }
> }
> }
> }
>
> int Assign_Array(PyArrayObject *self, PyObject *v) {
> PyObject *e;
> int l, r;
>
> if (!PySequence_Check(v)) {
> PyErr_SetString(PyExc_ValueError,"assignment from non-sequence");
> return -1;
> }
>
> l=PyObject_Length(v);
> if(l < 0) return -1;
>
> while(--l >= 0)
> {
> e=PySequence_GetItem(v,l);
> if (e == NULL) return -1;
> r = PySequence_SetItem((PyObject*)self,l,e);
> Py_DECREF(e);
> if(r == -1) return -1;
> }
> return 0;
> }
>
> PyObject *
> Array_FromSequence(PyObject *s, int type, int min_depth, int max_depth)
> {
> PyArrayObject *r;
> int nd, *d;
>
> if (!PySequence_Check(s)) {
> PyErr_SetString(PyExc_ValueError,"expect source sequence");
> return NULL;
> }
> if (!((nd=discover_depth(s,99, type == PyArray_OBJECT)) > 0)) {
> PyErr_SetString(PyExc_ValueError, "invalid input sequence");
> return NULL;
> }
> if ((max_depth && nd > max_depth) || (min_depth && nd < min_depth)) {
> PyErr_SetString(PyExc_ValueError, "invalid number of dimensions");
> return NULL;
> }
> if ((d=(int *)malloc(nd*sizeof(int))) == NULL) {
> PyErr_SetString(PyExc_MemoryError,"out of memory");
> }
> if(discover_dimensions(s,nd,d) == -1)
> {
> free(d);
> return 0;
> }
>
> if (type == PyArray_CHAR && nd > 0 && d[nd-1] == 1) {
> nd = nd-1;
> }
>
> r=(PyArrayObject*)PyArray_FromDims(nd,d,type);
> free(d);
> if(! r) return NULL;
> if(Assign_Array(r,s) == -1)
> {
> Py_DECREF(r);
> return NULL;
> }
> return (PyObject*)r;
> }
>
> PyObject *PyArray_FromScalar(PyObject *op, int type) {
> PyArrayObject *ret;
>
> if ((ret = (PyArrayObject *)PyArray_FromDims(0, NULL, type)) == NULL)
> return NULL;
>
> ret->descr->setitem(op, ret->data);
>
> if (PyErr_Occurred()) {
> array_dealloc(ret);
> return NULL;
> } else {
> return (PyObject *)ret;
> }
> }
>
> PyObject *PyArray_Cast(PyArrayObject *mp, int type) {
> PyArrayObject *rp, *tmp;
>
> if (mp->descr->type_num == PyArray_OBJECT) {
> return PyArray_FromObject((PyObject *)mp, type, mp->nd, mp->nd);
> }
>
> if ((tmp = PyArray_CONTIGUOUS(mp)) == NULL) return NULL;
>
> rp = (PyArrayObject *)PyArray_FromDims(tmp->nd, tmp->dimensions, type);
> mp->descr->cast[(int)rp->descr->type_num](tmp->data, 1, rp->data, 1,
SIZE(mp));
>
> Py_DECREF(tmp);
> return (PyObject *)rp;
> }
>
> #define ByCopy 1
> #define Contiguous 2
>
> PyObject *array_fromobject(PyObject *op_in, int type, int min_depth, int
max_depth, int flags) {
> PyObject *r, *op;
>
> r = NULL;
>
> if (!PyArray_Check(op_in) && PyObject_HasAttrString(op_in, "__array__")) {
> PyObject *ap, *arglist;
>
> if (type == PyArray_NOTYPE) {
> arglist = Py_BuildValue("()");
> } else {
> arglist = Py_BuildValue("(c)", type);
> }
> ap = PyObject_GetAttrString(op_in, "__array__");
> op = PyEval_CallObject(ap, arglist);
> Py_DECREF(ap);
> Py_DECREF(arglist);
>
> if (op == NULL) return NULL;
> } else {
> op = op_in;
> Py_INCREF(op);
> }
>
> if (type == PyArray_NOTYPE) {
> type = PyArray_ObjectType(op, 0);
> }
>
> if (PyArray_Check(op) && (((PyArrayObject *)op)->descr->type_num !=
PyArray_OBJECT ||
> type == PyArray_OBJECT || type == 'O')) {
> if ((((PyArrayObject *)op)->descr->type_num == type ||
> ((PyArrayObject *)op)->descr->type == type)) {
> if ((flags & ByCopy) || (flags&Contiguous && !ISCONTIGUOUS((PyArrayObject
*)op))) {
> r = PyArray_Copy((PyArrayObject *)op);
> } else {
> Py_INCREF(op);
> r = op;
> }
> } else {
> if (type > PyArray_NTYPES) type = PyArray_DescrFromType(type)->type_num;
> if (PyArray_CanCastSafely(((PyArrayObject *)op)->descr->type_num, type))
> r = PyArray_Cast((PyArrayObject *)op, type);
> else {
> PyErr_SetString(PyExc_TypeError, "Array can not be safely cast to required
type");
> r = NULL;
> }
> }
> } else {
> r = Array_FromSequence(op, type, min_depth,max_depth);
> if (r == NULL && min_depth <= 0) {
> PyErr_Clear();
> r = PyArray_FromScalar(op, type);
> }
> }
> Py_DECREF(op);
> if (r == NULL) return NULL;
>
> if(!PyArray_Check(r)) {
> PyErr_SetString(PyExc_ValueError, "Internal error array_fromobject not
producing an array");
> return NULL;
> }
> if (min_depth != 0 && ((PyArrayObject *)r)->nd < min_depth) {
> Py_DECREF(r);
> PyErr_SetString(PyExc_ValueError,
> "Object of too small depth for desired array");
> return NULL;
> }
> if (max_depth != 0 && ((PyArrayObject *)r)->nd > max_depth) {
> Py_DECREF(r);
> PyErr_SetString(PyExc_ValueError,
> "Object too deep for desired array");
> return NULL;
> }
> return r;
> }
>
>
>
> PyObject *PyArray_FromObject(PyObject *op, int type, int min_depth, int
max_depth) {
> return array_fromobject(op, type, min_depth, max_depth, 0);
> }
> PyObject *PyArray_ContiguousFromObject(PyObject *op, int type, int
min_depth, int max_depth) {
> return array_fromobject(op, type, min_depth, max_depth, Contiguous);
> }
> PyObject *PyArray_CopyFromObject(PyObject *op, int type, int min_depth,
int max_depth) {
> return array_fromobject(op, type, min_depth, max_depth, ByCopy);
> }
>
> extern int PyArray_CanCastSafely(int fromtype, int totype) {
> if (fromtype == totype) return 1;
> if (totype == PyArray_OBJECT) return 1;
> if (fromtype == PyArray_OBJECT) return 0;
>
> switch(fromtype) {
> case PyArray_CHAR:
> return 0;
> case PyArray_UBYTE:
> return (totype >= PyArray_SHORT);
> case PyArray_SBYTE:
> case PyArray_SHORT:
> return (totype > fromtype);
> case PyArray_INT:
> case PyArray_LONG:
> /*Allow Longs to be cast to Ints, not safe on 64-bit machines!*/
> return (totype >= PyArray_INT) && (totype != PyArray_FLOAT);
> case PyArray_FLOAT:
> return (totype > PyArray_FLOAT);
> case PyArray_DOUBLE:
> return (totype == PyArray_CDOUBLE);
> case PyArray_CFLOAT:
> return (totype == PyArray_CDOUBLE);
> case PyArray_CDOUBLE:
> return 0;
> default:
> return 0;
> }
> }
>
> extern int PyArray_As1D(PyObject **op, char **ptr, int *d1, int typecode)
{
> PyArrayObject *ap;
>
> if ((ap = (PyArrayObject *)PyArray_ContiguousFromObject(*op, typecode, 1,
> 1)) == NULL)
> return -1;
>
> *op = (PyObject *)ap;
> *ptr = ap->data;
> *d1 = ap->dimensions[0];
> return 0;
> }
>
> extern int PyArray_As2D(PyObject **op, char ***ptr, int *d1, int *d2, int
typecode) {
> PyArrayObject *ap;
> int i, n, size;
> char **data;
>
> if ((ap = (PyArrayObject *)PyArray_ContiguousFromObject(*op, typecode, 2,
2)) == NULL)
> return -1;
>
> n = ap->dimensions[0];
> data = (char **)malloc(n*sizeof(char *));
> size = ap->descr->elsize * ap->dimensions[1];
> for(i=0; i<n; i++) {
> data[i] = ap->data + i*ap->strides[0];
> }
> *op = (PyObject *)ap;
> *ptr = data;
> *d1 = ap->dimensions[0];
> *d2 = ap->dimensions[1];
> return 0;
> }
>
>
> extern int PyArray_Free(PyObject *op, char *ptr) {
> PyArrayObject *ap = (PyArrayObject *)op;
> int i, n;
>
> if (ap->nd > 2) return -1;
> if (ap->nd == 3) {
> n = ap->dimensions[0];
> for (i=0; i<n; i++) {
> free(((char **)ptr)[i]);
> }
> }
> if (ap->nd >= 2) {
> free(ptr);
> }
> Py_DECREF(ap);
> return 0;
> }
>
> extern PyObject * PyArray_Reshape(PyArrayObject *self, PyObject *shape) {
> int i, n, s_original, i_unknown, s_known;
> int *dimensions;
> PyArrayObject *ret;
>
> if (!PyArray_ISCONTIGUOUS(self)) {
> PyErr_SetString(PyExc_ValueError, "reshape only works on contiguous
matrices");
> return NULL;
> }
>
> if (PyArray_As1D(&shape, (char **)&dimensions, &n, PyArray_INT) == -1)
return NULL;
>
> s_known = 1;
> i_unknown = -1;
> for(i=0; i<n; i++) {
> if (dimensions[i] < 0) {
> if (i_unknown == -1) {
> i_unknown = i;
> } else {
> PyErr_SetString(PyExc_ValueError, "can only specify one unknown
dimension");
> goto fail;
> }
> } else {
> s_known *= dimensions[i];
> }
> }
>
> s_original = PyArray_SIZE(self);
>
> if (i_unknown >= 0) {
> if ((s_known == 0) || (s_original % s_known != 0)) {
> PyErr_SetString(PyExc_ValueError, "total size of new array must be
unchanged");
> goto fail;
> }
> dimensions[i_unknown] = s_original/s_known;
> } else {
> if (s_original != s_known) {
> PyErr_SetString(PyExc_ValueError, "total size of new array must be
unchanged");
> goto fail;
> }
> }
> if ((ret = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(n, dimensions,
> self->descr,
> self->data)) == NULL)
> goto fail;
>
> Py_INCREF(self);
> ret->base = (PyObject *)self;
>
> PyArray_Free(shape, (char *)dimensions);
> return (PyObject *)ret;
>
> fail:
> PyArray_Free(shape, (char *)dimensions);
> return NULL;
> }
>
> extern PyObject *PyArray_Take(PyObject *self0, PyObject *indices0, int
axis) {
> PyArrayObject *self, *indices, *ret;
> int nd, shape[MAX_DIMS];
> int i, j, chunk, n, m, max_item, tmp;
> char *src, *dest;
>
> indices = ret = NULL;
> self = (PyArrayObject *)PyArray_ContiguousFromObject(self0,
PyArray_NOTYPE, 1, 0);
> if (self == NULL) return NULL;
>
> if (axis < 0) axis = axis + self->nd;
> if ((axis < 0) || (axis >= self->nd)) {
> PyErr_SetString(PyExc_ValueError, "Invalid axis for this array");
> goto fail;
> }
>
> indices = (PyArrayObject *)PyArray_ContiguousFromObject(indices0,
PyArray_LONG, 1, 0);
> if (indices == NULL) goto fail;
>
> n = m = chunk = 1;
> nd = self->nd + indices->nd - 1;
> for (i=0; i< nd; i++) {
> if (i < axis) {
> shape[i] = self->dimensions[i];
> n *= shape[i];
> } else {
> if (i < axis+indices->nd) {
> shape[i] = indices->dimensions[i-axis];
> m *= shape[i];
> } else {
> shape[i] = self->dimensions[i-indices->nd+1];
> chunk *= shape[i];
> }
> }
> }
> ret = (PyArrayObject *)PyArray_FromDims(nd, shape, self->descr->type_num);
>
> if (ret == NULL) goto fail;
>
> max_item = self->dimensions[axis];
> chunk = chunk * ret->descr->elsize;
> src = self->data;
> dest = ret->data;
>
> for(i=0; i<n; i++) {
> for(j=0; j<m; j++) {
> tmp = ((long *)(indices->data))[j];
> if (tmp < 0) tmp = tmp+max_item;
> if ((tmp < 0) || (tmp >= max_item)) {
> PyErr_SetString(PyExc_IndexError, "Index out of range for array");
> goto fail;
> }
> memcpy(dest, src+tmp*chunk, chunk);
> dest += chunk;
> }
> src += chunk*max_item;
> }
>
> PyArray_INCREF(ret);
>
> Py_XDECREF(indices);
> Py_XDECREF(self);
>
> return (PyObject *)ret;
>
>
> fail:
> Py_XDECREF(ret);
> Py_XDECREF(indices);
> Py_XDECREF(self);
> return NULL;
> }
>
----------------------------------------------------------------------------
----
> #ifndef Py_ARRAYOBJECT_H
> #define Py_ARRAYOBJECT_H
> #ifdef __cplusplus
> extern "C" {
> #endif
>
> #ifndef _ARRAY_MODULE
> #ifdef COMPILING_NUMPY
> #include "ExtensionClassMacros.h"
> #else
> #include "Numeric/ExtensionClassMacros.h"
> #endif
> #endif
>
> #ifdef COMPILING_NUMPY
> #include "ExtensionClass.h"
> #else
> #include "Numeric/ExtensionClass.h"
> #endif
>
>
> #define MAX_ELSIZE 16
>
> enum PyArray_TYPES { PyArray_CHAR, PyArray_UBYTE, PyArray_SBYTE,
> PyArray_SHORT, PyArray_INT, PyArray_LONG,
> PyArray_FLOAT, PyArray_DOUBLE,
> PyArray_CFLOAT, PyArray_CDOUBLE,
> PyArray_OBJECT,
> PyArray_NTYPES, PyArray_NOTYPE};
>
> typedef void (PyArray_VectorUnaryFunc) Py_FPROTO((char *, int, char *,
int, int));
>
> typedef PyObject * (PyArray_GetItemFunc) Py_FPROTO((char *));
> typedef int (PyArray_SetItemFunc) Py_FPROTO((PyObject *, char *));
>
> typedef struct {
> PyArray_VectorUnaryFunc *cast[PyArray_NTYPES]; /* Functions to cast to
*/
> /* all other types */
> PyArray_GetItemFunc *getitem;
> PyArray_SetItemFunc *setitem;
>
> int type_num, elsize;
> char *one, *zero;
> char type;
>
> } PyArray_Descr;
>
> #define CONTIGUOUS 1
> #define OWN_DIMENSIONS 2
> #define OWN_STRIDES 4
> #define OWN_DATA 8
>
> typedef struct {
> PyObject_HEAD
> char *data;
> int nd;
> int *dimensions, *strides;
> PyObject *base;
> PyArray_Descr *descr;
> int flags;
> #ifndef NUMPY_NOEXTRA
> PyObject* attributes; /* for user-defined information */
> #endif
> } PyArrayObject;
>
>
> /*
> * C API
> */
>
> /* Type definitions */
>
> #define PyArray_Type_NUM 0
>
> /* Function definitions */
>
> /* The following are not intended for use in user code, only needed by
umath. */
> /* If you write your own match library, you might want this function. */
> #define PyArray_SetNumericOps_RET int
> #define PyArray_SetNumericOps_PROTO Py_FPROTO((PyObject *))
> #define PyArray_SetNumericOps_NUM 1
>
> #define PyArray_INCREF_RET int
> #define PyArray_INCREF_PROTO Py_FPROTO((PyArrayObject *ap))
> #define PyArray_INCREF_NUM 2
>
> #define PyArray_XDECREF_RET int
> #define PyArray_XDECREF_PROTO Py_FPROTO((PyArrayObject *ap))
> #define PyArray_XDECREF_NUM 3
>
> /* Export the array error object. Is this a good idea? */
> #define PyArrayError_RET PyObject *
> #define PyArrayError_PROTO Py_FPROTO(())
> #define PyArrayError_NUM 4
>
> /* Set the array print function to be a python function */
> #define PyArray_SetStringFunction_RET void
> #define PyArray_SetStringFunction_PROTO Py_FPROTO((PyObject *op, int
repr))
> #define PyArray_SetStringFunction_NUM 5
>
> /* Get the PyArray_Descr structure for a typecode */
> #define PyArray_DescrFromType_RET PyArray_Descr *
> #define PyArray_DescrFromType_PROTO Py_FPROTO((int))
> #define PyArray_DescrFromType_NUM 6
>
> /* Cast an array to a different type */
> #define PyArray_Cast_RET PyObject *
> #define PyArray_Cast_PROTO Py_FPROTO((PyArrayObject *, int))
> #define PyArray_Cast_NUM 7
>
> /* Check the type coercion rules */
> #define PyArray_CanCastSafely_RET int
> #define PyArray_CanCastSafely_PROTO Py_FPROTO((int fromtype, int totype))
> #define PyArray_CanCastSafely_NUM 8
>
> /* Return the typecode to use for an object if it was an array */
> #define PyArray_ObjectType_RET int
> #define PyArray_ObjectType_PROTO Py_FPROTO((PyObject *, int))
> #define PyArray_ObjectType_NUM 9
>
> #define _PyArray_multiply_list_RET int
> #define _PyArray_multiply_list_PROTO Py_FPROTO((int *lp, int n))
> #define _PyArray_multiply_list_NUM 10
>
>
> /* The following defines the C API for the array object for most users */
>
> #define PyArray_SIZE(mp) (_PyArray_multiply_list((mp)->dimensions,
(mp)->nd))
> #define PyArray_NBYTES(mp) ((mp)->descr->elsize * PyArray_SIZE(mp))
> /* Obviously this needs some work. */
> #define PyArray_ISCONTIGUOUS(m) ((m)->flags & CONTIGUOUS)
>
>
> /* Return the size in number of items of an array */
> #define PyArray_Size_RET int
> #define PyArray_Size_PROTO Py_FPROTO((PyObject *))
> #define PyArray_Size_NUM 11
>
>
> /* Array creation functions */
> /* new_array = PyArray_FromDims(n_dimensions, dimensions[n_dimensions],
item_type); */
> #define PyArray_FromDims_RET PyObject *
> #define PyArray_FromDims_PROTO Py_PROTO((int, int *, int))
> #define PyArray_FromDims_NUM 12
>
> /* array_from_existing_data = PyArray_FromDimsAndData(n_dims,
dims[n_dims], item_type, old_data); */
> /* WARNING: using PyArray_FromDimsAndData is not reccommended! It should
only be used to refer to */
> /* global arrays that will never be freed (like FORTRAN common blocks). */
> #define PyArray_FromDimsAndData_RET PyObject *
> #define PyArray_FromDimsAndData_PROTO Py_PROTO((int, int *, int, char *))
> #define PyArray_FromDimsAndData_NUM 13
>
> /* Initialize from a python object. */
>
> /* PyArray_ContiguousFromObject(object, typecode, min_dimensions,
max_dimensions) */
> /* if max_dimensions = 0, then any number of dimensions are allowed. */
> /* If you want an exact number of dimensions, you should use
max_dimensions */
> /* = min_dimensions. */
>
> #define PyArray_ContiguousFromObject_RET PyObject *
> #define PyArray_ContiguousFromObject_PROTO Py_FPROTO((PyObject *, int,
int, int))
> #define PyArray_ContiguousFromObject_NUM 14
>
> /* Same as contiguous, except guarantees a copy of the original data */
> #define PyArray_CopyFromObject_RET PyObject *
> #define PyArray_CopyFromObject_PROTO Py_FPROTO((PyObject *, int, int,
int))
> #define PyArray_CopyFromObject_NUM 15
>
> /* Shouldn't be used unless you know what you're doing and are not scared
by discontiguous arrays */
> #define PyArray_FromObject_RET PyObject *
> #define PyArray_FromObject_PROTO Py_FPROTO((PyObject *, int, int, int))
> #define PyArray_FromObject_NUM 16
>
> /* Return either an array, or if passed a 0d array return the appropriate
python scalar */
> #define PyArray_Return_RET PyObject *
> #define PyArray_Return_PROTO Py_FPROTO((PyArrayObject *))
> #define PyArray_Return_NUM 17
>
> #define PyArray_Reshape_RET PyObject *
> #define PyArray_Reshape_PROTO Py_FPROTO((PyArrayObject *ap, PyObject
*shape))
> #define PyArray_Reshape_NUM 18
>
> #define PyArray_Copy_RET PyObject *
> #define PyArray_Copy_PROTO Py_FPROTO((PyArrayObject *ap))
> #define PyArray_Copy_NUM 19
>
> #define PyArray_Take_RET PyObject *
> #define PyArray_Take_PROTO Py_FPROTO((PyObject *ap, PyObject *items, int
axis))
> #define PyArray_Take_NUM 20
>
> /*Getting arrays in various useful forms. */
> #define PyArray_As1D_RET int
> #define PyArray_As1D_PROTO Py_FPROTO((PyObject **op, char **ptr, int *d1,
int typecode))
> #define PyArray_As1D_NUM 21
>
> #define PyArray_As2D_RET int
> #define PyArray_As2D_PROTO Py_FPROTO((PyObject **op, char ***ptr, int *d1,
int *d2, int typecode))
> #define PyArray_As2D_NUM 22
>
> #define PyArray_Free_RET int
> #define PyArray_Free_PROTO Py_FPROTO((PyObject *op, char *ptr))
> #define PyArray_Free_NUM 23
>
> /* array_from_existing_data = PyArray_FromDimsAndDataAndDescr(n_dims,
dims[n_dims], descr, old_data); */
> /* WARNING: using PyArray_FromDimsAndDataAndDescr is not reccommended! It
should only be used to refer to */
> /* global arrays that will never be freed (like FORTRAN common blocks). */
> #define PyArray_FromDimsAndDataAndDescr_RET PyObject *
> #define PyArray_FromDimsAndDataAndDescr_PROTO Py_PROTO((int, int *,
PyArray_Descr *, char *))
> #define PyArray_FromDimsAndDataAndDescr_NUM 24
>
> /* Total number of C API pointers */
> #define PyArray_API_pointers 25
>
>
> #ifdef _ARRAY_MODULE
> extern PyExtensionClass PyArray_Type;
>
> /* define PyArray_Check(op) ((op)->ob_type == &PyArray_Type) */
> #define PyArray_Check(op) (ExtensionClassSubclassInstance_Check((op),
&(PyArray_Type)))
>
>
> extern PyArray_SetNumericOps_RET PyArray_SetNumericOps \
> PyArray_SetNumericOps_PROTO;
> extern PyArray_INCREF_RET PyArray_INCREF PyArray_INCREF_PROTO;
> extern PyArray_XDECREF_RET PyArray_XDECREF PyArray_XDECREF_PROTO;
> extern PyArrayError_RET PyArrayError PyArrayError_PROTO;
> extern PyArray_SetStringFunction_RET PyArray_SetStringFunction \
> PyArray_SetStringFunction_PROTO;
> extern PyArray_DescrFromType_RET PyArray_DescrFromType \
> PyArray_DescrFromType_PROTO;
> extern PyArray_Cast_RET PyArray_Cast PyArray_Cast_PROTO;
> extern PyArray_CanCastSafely_RET PyArray_CanCastSafely \
> PyArray_CanCastSafely_PROTO;
> extern PyArray_ObjectType_RET PyArray_ObjectType PyArray_ObjectType_PROTO;
> extern _PyArray_multiply_list_RET _PyArray_multiply_list \
> _PyArray_multiply_list_PROTO;
> extern PyArray_Size_RET PyArray_Size PyArray_Size_PROTO;
> extern PyArray_FromDims_RET PyArray_FromDims PyArray_FromDims_PROTO;
> extern PyArray_FromDimsAndData_RET PyArray_FromDimsAndData \
> PyArray_FromDimsAndData_PROTO;
> extern PyArray_FromDimsAndDataAndDescr_RET PyArray_FromDimsAndDataAndDescr
\
> PyArray_FromDimsAndDataAndDescr_PROTO;
> extern PyArray_ContiguousFromObject_RET PyArray_ContiguousFromObject \
> PyArray_ContiguousFromObject_PROTO;
> extern PyArray_CopyFromObject_RET PyArray_CopyFromObject \
> PyArray_CopyFromObject_PROTO;
> extern PyArray_FromObject_RET PyArray_FromObject PyArray_FromObject_PROTO;
> extern PyArray_Return_RET PyArray_Return PyArray_Return_PROTO;
> extern PyArray_Reshape_RET PyArray_Reshape PyArray_Reshape_PROTO;
> extern PyArray_Copy_RET PyArray_Copy PyArray_Copy_PROTO;
> extern PyArray_Take_RET PyArray_Take PyArray_Take_PROTO;
> extern PyArray_As1D_RET PyArray_As1D PyArray_As1D_PROTO;
> extern PyArray_As2D_RET PyArray_As2D PyArray_As2D_PROTO;
> extern PyArray_Free_RET PyArray_Free PyArray_Free_PROTO;
>
>
> #else
>
> /* C API address pointer */
> #if defined(NO_IMPORT) || defined(NO_IMPORT_ARRAY)
> extern void **PyArray_API;
> #else
> void **PyArray_API;
> #endif
>
> #define PyArray_Check(op) \
> MakeSureExtensionClassImported && ExtensionClassSubclassInstance_Check(op,
(PyExtensionClass *)PyArray_API[PyArray_Type_NUM])
> #define PyArray_Type *(PyExtensionClass *)PyArray_API[PyArray_Type_NUM]
>
> #define PyArray_SetNumericOps \
> (*(PyArray_SetNumericOps_RET (*)PyArray_SetNumericOps_PROTO) \
> PyArray_API[PyArray_SetNumericOps_NUM])
> #define PyArray_INCREF \
> (*(PyArray_INCREF_RET (*)PyArray_INCREF_PROTO) \
> PyArray_API[PyArray_INCREF_NUM])
> #define PyArray_XDECREF \
> (*(PyArray_XDECREF_RET (*)PyArray_XDECREF_PROTO) \
> PyArray_API[PyArray_XDECREF_NUM])
> #define PyArrayError \
> (*(PyArrayError_RET (*)PyArrayError_PROTO) \
> PyArray_API[PyArrayError_NUM])
> #define PyArray_SetStringFunction \
> (*(PyArray_SetStringFunction_RET (*)PyArray_SetStringFunction_PROTO) \
> PyArray_API[PyArray_SetStringFunction_NUM])
> #define PyArray_DescrFromType \
> (*(PyArray_DescrFromType_RET (*)PyArray_DescrFromType_PROTO) \
> PyArray_API[PyArray_DescrFromType_NUM])
> #define PyArray_Cast \
> (*(PyArray_Cast_RET (*)PyArray_Cast_PROTO) \
> PyArray_API[PyArray_Cast_NUM])
> #define PyArray_CanCastSafely \
> (*(PyArray_CanCastSafely_RET (*)PyArray_CanCastSafely_PROTO) \
> PyArray_API[PyArray_CanCastSafely_NUM])
> #define PyArray_ObjectType \
> (*(PyArray_ObjectType_RET (*)PyArray_ObjectType_PROTO) \
> PyArray_API[PyArray_ObjectType_NUM])
> #define _PyArray_multiply_list \
> (*(_PyArray_multiply_list_RET (*)_PyArray_multiply_list_PROTO) \
> PyArray_API[_PyArray_multiply_list_NUM])
> #define PyArray_Size \
> (*(PyArray_Size_RET (*)PyArray_Size_PROTO) \
> PyArray_API[PyArray_Size_NUM])
> #define PyArray_FromDims \
> (*(PyArray_FromDims_RET (*)PyArray_FromDims_PROTO) \
> PyArray_API[PyArray_FromDims_NUM])
> #define PyArray_FromDimsAndData \
> (*(PyArray_FromDimsAndData_RET (*)PyArray_FromDimsAndData_PROTO) \
> PyArray_API[PyArray_FromDimsAndData_NUM])
> #define PyArray_FromDimsAndDataAndDescr \
> (*(PyArray_FromDimsAndDataAndDescr_RET
(*)PyArray_FromDimsAndDataAndDescr_PROTO) \
> PyArray_API[PyArray_FromDimsAndDataAndDescr_NUM])
> #define PyArray_ContiguousFromObject \
> (*(PyArray_ContiguousFromObject_RET
(*)PyArray_ContiguousFromObject_PROTO) \
> PyArray_API[PyArray_ContiguousFromObject_NUM])
> #define PyArray_CopyFromObject \
> (*(PyArray_CopyFromObject_RET (*)PyArray_CopyFromObject_PROTO) \
> PyArray_API[PyArray_CopyFromObject_NUM])
> #define PyArray_FromObject \
> (*(PyArray_FromObject_RET (*)PyArray_FromObject_PROTO) \
> PyArray_API[PyArray_FromObject_NUM])
> #define PyArray_Return \
> (*(PyArray_Return_RET (*)PyArray_Return_PROTO) \
> PyArray_API[PyArray_Return_NUM])
> #define PyArray_Reshape \
> (*(PyArray_Reshape_RET (*)PyArray_Reshape_PROTO) \
> PyArray_API[PyArray_Reshape_NUM])
> #define PyArray_Copy \
> (*(PyArray_Copy_RET (*)PyArray_Copy_PROTO) \
> PyArray_API[PyArray_Copy_NUM])
> #define PyArray_Take \
> (*(PyArray_Take_RET (*)PyArray_Take_PROTO) \
> PyArray_API[PyArray_Take_NUM])
> #define PyArray_As1D \
> (*(PyArray_As1D_RET (*)PyArray_As1D_PROTO) \
> PyArray_API[PyArray_As1D_NUM])
> #define PyArray_As2D \
> (*(PyArray_As2D_RET (*)PyArray_As2D_PROTO) \
> PyArray_API[PyArray_As2D_NUM])
> #define PyArray_Free \
> (*(PyArray_Free_RET (*)PyArray_Free_PROTO) \
> PyArray_API[PyArray_Free_NUM])
>
> #define import_array() \
> { \
> PyObject *numpy = PyImport_ImportModule("_numpy"); \
> if (numpy != NULL) { \
> PyObject *module_dict = PyModule_GetDict(numpy); \
> PyObject *c_api_object = PyDict_GetItemString(module_dict,
"_ARRAY_API"); \
> if (PyCObject_Check(c_api_object)) { \
> PyArray_API = (void **)PyCObject_AsVoidPtr(c_api_object); \
> } \
> } \
> }
>
> #endif
>
> #ifdef __cplusplus
> }
> #endif
> #endif /* !Py_ARRAYOBJECT_H */
>