A. M. Archibald wrote:
Why not write the algorithm in C?
I did just that a while back, for Numeric. I've enclosed the code for
reference.
Unfortunately, I never did figure out an efficient way to write this
sort of thing for all types, so it only does doubles. Also, it does a
bunch of special casing for discontiguous vs. contiguous arrays, and
clipping to an array vs a scaler for the min and max arguments.
Do what you will with the code...
look for NumericExtras_fastclip in with the rest of the stuff in there
-- much of which is now obsolete -- I haven't touched this is a few years.
-Chris
NOTE: is there a clean and efficient way to write custom functions of
this sort for all types and both contiguous and discontiguous arrays? I
guess I'm imaging some smart iterators and a micro version of C++
templates -- or maybe use C++ templates themselves.
--
Christopher Barker, Ph.D.
Oceanographer
Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception
[EMAIL PROTECTED]
#include "Python.h"
#include "Numeric/arrayobject.h"
// These are little macros I use to access array values for specific rank arrays
#define ARRAYVAL0(aType,a) ( *(aType *)(a->data))
#define ARRAYVAL1(aType,a,i) ( *(aType *)(a->data + (i)*a->strides[0]))
#define ARRAYVAL2(aType,a,i,j) ( *(aType *)(a->data + (i)*a->strides[0] +
(j)*a->strides[1]))
#define ARRAYVAL3(aType,a,i,j,k) ( *(aType *)(a->data + (i)*a->strides[0] +
(j)*a->strides[1] + (k)*a->strides[2]))
// These are some macros borrowed from arrayobject.c
//#define SIZE(mp) (_PyArray_multiply_list((mp)->dimensions, (mp)->nd))
int AreSameShape(PyArrayObject *Array1, PyArrayObject *Array2){
int i;
if (Array1->nd != Array2->nd){
return 0;
}
for (i = 0; i < Array1->nd; i++){
if (Array1->dimensions[i] != Array2->dimensions[i]){
return 0;
}
}
return 1;
}
int ComputeOffset(PyArrayObject *Array, int *indexes){
int i, offset = 0;
for (i = 0; i < Array->nd; i++)
{
offset += indexes[i]*Array->strides[i];
}
return offset;
}
void ReplaceValue(PyArrayObject *Array, PyArrayObject *Min, PyArrayObject *Max,
int *indexes){
// This function does the actual replacement of the values.
// it allows Min and Max to be one element, or the same shape as Array
// the offsets are computed separately, so as long as they are the same shape
// it should work even if some are contiguous, and some not, for example
double min, max;
int offset;
offset = ComputeOffset(Array,indexes);
if (PyArray_Size((PyObject *)Min) == 1){
min = ARRAYVAL0(double,Min);
}
else {
min = *(double *)(Min->data + ComputeOffset(Min,indexes));
}
if (PyArray_Size((PyObject *)Max) == 1){
max = ARRAYVAL0(double,Max);
}
else {
max = *(double *)(Max->data + ComputeOffset(Max,indexes));
}
if (*(double *)(Array->data + offset) > max){
*(double *)(Array->data + offset) = max ;
}
else if (*(double *)(Array->data + offset) < min){
*(double *)(Array->data + offset) = min;
}
return;
}
static char doc_fastclip[] =
"fastclip(a,min,max): a is clipped so that there are no values \n"
"less than min, or greater than max.\n\n"
"It only works on PyArrays of type Float.\n"
"min and max can be either scalar or the same size as a.\n"
"If a, min, and max are all contiguous (or scalar) then it is a LOT faster\n\n"
"Note that if min > max, the results are undetermined";
static PyObject * NumericExtras_fastclip(PyObject *self, PyObject *args)
{
int NumArray, elementsize, i;
int *indexes;
PyObject *InMin;
PyObject *InMax;
PyArrayObject *Array;
PyArrayObject *Min;
PyArrayObject *Max;
if (!PyArg_ParseTuple(args, "O!OO",
&PyArray_Type, &Array,
&InMin,
&InMax))
{
return NULL;
}
// check types of input
if (Array->descr->type_num != PyArray_DOUBLE){
PyErr_SetString (PyExc_ValueError,
"a must be a NumPy array of type Float");
return NULL;
}
// convert min and max to double arrays:
// if they don't convert, the input values are not valid
// also check if they are the right size
Min = (PyArrayObject *) PyArray_FromObject(InMin, PyArray_DOUBLE, 0, 0);
if (Min == NULL){
PyErr_SetString (PyExc_ValueError,
"min must be an object that can be converted to an array
of Floats");
return NULL;
}
if (!((PyArray_Size((PyObject *)Min) == 1) || AreSameShape(Min, Array))){
PyErr_SetString (PyExc_ValueError,
"min must be either a scalar or the same size as a");
Py_DECREF(Min);
return NULL;
}
Max = (PyArrayObject *) PyArray_FromObject(InMax, PyArray_DOUBLE, 0, 0);
if (Max == NULL){
PyErr_SetString (PyExc_ValueError,
"max must be an object that can be converted to an array
of Floats");
Py_DECREF(Min);
return NULL;
}
if (!((PyArray_Size((PyObject *)Max) == 1) || AreSameShape(Max, Array))){
PyErr_SetString (PyExc_ValueError,
"max must be either a scalar or the same size as a");
Py_DECREF(Min);
Py_DECREF(Max);
return NULL;
}
// After all that input checking, switch between the contiguous and
non-contiguous cases.
if (PyArray_ISCONTIGUOUS(Array) && PyArray_ISCONTIGUOUS(Min) &&
PyArray_ISCONTIGUOUS(Max)){ //we have the contiguous case
// This seems pretty kludgy to have the four cases, but it was the first
// thing that came to mind that I was sure would work, and would accomidate
// either array or scalar arguments for min and max.
// it's also faster than checking if they are scalar inside the loop
NumArray = PyArray_Size((PyObject *)Array);
elementsize = sizeof(double);
if (PyArray_Size((PyObject *)Min) == 1 && PyArray_Size((PyObject *)Max) ==
1){ // both limits are scalar
for (i = 0; i < NumArray; i++) { // loop over each element
if (*(double *)(Array->data + i*elementsize) > ARRAYVAL0(double,Max)){
*(double *)(Array->data + i*elementsize) = ARRAYVAL0(double,Max);
}
else if (*(double *)(Array->data + i*elementsize) <
ARRAYVAL0(double,Min)){
*(double *)(Array->data + i*elementsize) = ARRAYVAL0(double,Min);
}
}
}
else if (PyArray_Size((PyObject *)Min) == 1 && PyArray_Size((PyObject
*)Max) == NumArray){ // Min is scalar
for (i = 0; i < NumArray; i++) { // loop over each element
if (*(double *)(Array->data + i*elementsize) > *(double *)(Max->data +
i*elementsize)){
*(double *)(Array->data + i*elementsize) = *(double *)(Max->data +
i*elementsize) ;
}
else if (*(double *)(Array->data + i*elementsize) <
ARRAYVAL0(double,Min)){
*(double *)(Array->data + i*elementsize) = ARRAYVAL0(double,Min);
}
}
}
else if (PyArray_Size((PyObject *)Max) == 1 && PyArray_Size((PyObject
*)Min) == NumArray){ // Max is scalar
for (i = 0; i < NumArray; i++) { // loop over each element
if (*(double *)(Array->data + i*elementsize) > ARRAYVAL0(double,Max)){
*(double *)(Array->data + i*elementsize) = ARRAYVAL0(double,Max);
}
else if (*(double *)(Array->data + i*elementsize) < *(double
*)(Min->data + i*elementsize)){
*(double *)(Array->data + i*elementsize) = *(double *)(Min->data +
i*elementsize) ;
}
}
}
else if (PyArray_Size((PyObject *)Min) == NumArray &&
PyArray_Size((PyObject *)Max) == NumArray){ // Neither is scalar
for (i = 0; i < NumArray; i++) { // loop over each element
if (*(double *)(Array->data + i*elementsize) > *(double *)(Max->data +
i*elementsize)){
*(double *)(Array->data + i*elementsize) = *(double *)(Max->data +
i*elementsize) ;
}
else if (*(double *)(Array->data + i*elementsize) < *(double
*)(Min->data + i*elementsize)){
*(double *)(Array->data + i*elementsize) = *(double *)(Min->data +
i*elementsize) ;
}
}
}
else { // The size of Min and/or Max don't match Array: we should never get
here, do to previous checks.
PyErr_SetString (PyExc_ValueError,
"min and max must be either scalar or the same shape as
a");
Py_DECREF(Min);
Py_DECREF(Max);
return NULL;
}
}
else { // this now the non-contiguous case: it is much slower, but presumably
could be optimised : suggestions gratefully excepted!
indexes = calloc(Array->nd,sizeof(int));
while (1){
ReplaceValue(Array, Min, Max, indexes);
i = (Array->nd) - 1;
while ((i >= 0) && (indexes[i]+1 == Array->dimensions[i])){
indexes[i] = 0;
i -= 1;
}
if (i<0) {
break;
}
indexes[i] = indexes[i]+1;
}
free(indexes);
}
Py_DECREF(Min);
Py_DECREF(Max);
return Py_BuildValue("");
}
static void byte_swap_data(void *p, int n, int size) {
/* p is a pointer to the data block
n is the size of the data block
size is the size of the elements in the data block
This code borrowed from arrayobject.c in Numeric
It is used by byteswap
*/
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_byteswap[] = "byteswap(m). Swaps the bytes in the contents of
array m.\n"
"Useful for reading data written by a machine with
a different byte order.\n"
"Operates on the array in place. Only works on
contiguous arrays.\n";
static PyObject * NumericExtras_byteswap(PyObject *self, PyObject *args) {
PyArrayObject *Array;
if (!PyArg_ParseTuple(args, "O!",
&PyArray_Type, &Array)) {
return NULL;
}
if (!PyArray_ISCONTIGUOUS(Array)){
PyErr_SetString (PyExc_ValueError,
"m must be contiguous");
return NULL;
}
if (Array->descr->type_num < PyArray_CFLOAT) {
byte_swap_data(Array->data, PyArray_SIZE(Array), Array->descr->elsize);
} else {
byte_swap_data(Array->data, PyArray_SIZE(Array)*2,
Array->descr->elsize/2);
}
return Py_BuildValue("");
}
static char doc_changetype[] =
"changetype(m,typecode). Changes the type, at the C level, of the block of\n"
"bytes in the data array of m. It does this in place. This is useful for\n"
"manipulating binary data that may have been read in from a file, without\n"
"making multiple copies of the data.\n\n"
"The array returned is always of rank 1"
"changetype(m,typecode) should have the same result as:\n\n"
"m = fromstring(m.tostring(),typecode)\n\n"
"it only works on contiguous arrays\n";
static PyObject * NumericExtras_changetype(PyObject *self, PyObject *args) {
PyArrayObject *Array;
//int *new_strides;
//int *new_dimensions;
PyArray_Descr *NewDescription;
int TotalBytes;
int NewElementSize;
char typecode;
if (!PyArg_ParseTuple(args, "O!c",
&PyArray_Type, &Array,
&typecode)) {
return NULL;
}
if (!PyArray_ISCONTIGUOUS(Array)){
PyErr_SetString (PyExc_ValueError,
"m must be contiguous");
return NULL;
}
TotalBytes = PyArray_NBYTES(Array);
// This will crash if an invalid typecode is passed in.
NewDescription = PyArray_DescrFromType(typecode);
NewElementSize = NewDescription->elsize;
if (TotalBytes % NewElementSize > 0){
PyErr_SetString (PyExc_ValueError,
"The input array is the wrong size for the requested
type");
return NULL;
}
Array->descr = NewDescription;
// does the old descr need to be freed?? how??
Array->nd = 1;
Array->strides[0] = NewElementSize;
Array->dimensions[0] = TotalBytes / NewElementSize;
return Py_BuildValue("");
}
static PyMethodDef NumericExtrasMethods[] = {
{"fastclip", NumericExtras_fastclip, METH_VARARGS, doc_fastclip},
{"byteswap", NumericExtras_byteswap, METH_VARARGS, doc_byteswap},
{"changetype", NumericExtras_changetype, METH_VARARGS, doc_changetype},
{NULL, NULL} /* Sentinel */
};
void initNumericExtras(void){
(void) Py_InitModule("NumericExtras", NumericExtrasMethods);
import_array()
}
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
