[Numpy-discussion] API tricks any my new "clip" function

[Chris B. <chr...@ho...>]

Hi all,

I wrote last week looking for hints for writing a new clip function,
that would be faster than the existing one. I got some suggestions, and
have written the function. It works, but I'd love to get a little code
review from the more experienced folks here. I have a few specific
questions:

A) Are there any tricks to making a function work with multiple types? I
currently have it working with only double arrays, but it would be great
for it ot be universal (and it could then be added to the main NumPy
distribution, I suppose) I seem tohave to typecast the data to do the
comparison I want, so I can't figure out how to make it general. For
example, I have this line:

if (*(double *)(Array->data + i*elementsize) > *(double *)(Max->data +
i*elementsize))
{
          *(double *)(Array->data + i*elementsize) = *(double
*)(Max->data + i*elementsize) ;
}

How could I do that without the explicit typecasts? I'm pretty sure I
could make the rest of the routine general.

B) I was trying to figure out how to loop through all the elements of a
non-contiguous array. I got two hints:
"Rob W. W. Hooft" wrote:

> Never done this in C for numpy yet, but I normally program this
> along the following lines:
> 
> n=[2,3,2,3]
> x=[0,0,0,0]
> while 1:
>     print x
>     i=len(n)-1
>     while i>=0 and x[i]+1==n[i]:
>        x[i]=0
>        i=i-1
>     if i<0:
>        break
>     x[i]=x[i]+1
> 
> It is always going to be slower than the straight loop for contiguous
> arrays. 

This is what I did, and it is a LOT slower (about a factor of ten). See
the enclosed C code. There has got to be a way to optimize this!

>So if you want to do it properly, you'd have to check whether
> the array is contiguous, and if it is, use the fast way.

Since I had already written code for the contiguous case, it was easy to
special case it.

Pete Shinners wrote:
> it's easiest to use recursion than nested loops. then you can
> support any number of dimensions. that should point you in the
> right direction :]

When I read this, I thought: "well, duh!", but then I couldn't figure
out how to do it. Would there be any advantage over the above loops
anyway?

Is there either a more elgant or faster way to loop through all the
elements of these three arrays?

C) Where is the DL_EXPORT macro defined? I havn't been able to find it
with grep yet. Does it work on all NumPy supported platforms?

D) Is there an exising function or Macro for checking if two arrays are
the same shape? I wrote one, but I imagine it exists.

F) I needed to work with input that could be any Python number, or an
Array of the right size. What I did was call: PyArray_FromObject on the
min and max inputs, so that I would get a PyArray, then I could check if
it was a rank-o array, or the right size. I havn't been able to decide
if this is a nifty trcik, or a kludge. Any other suggestions?

E) Have I got reference counting right? I tried to Py_DECREF all my temp
arrays, but I'm not sure I did it right.

Sorry to enclose the code if you have no interest in it, but it's not
really that big.

-Chris

-- 
Christopher Barker,
Ph.D.                                                           
Chr...@ho...                 ---           ---           ---
http://members.home.net/barkerlohmann ---@@       -----@@       -----@@
                                   ------@@@     ------@@@     ------@@@
Oil Spill Modeling                ------   @    ------   @   ------   @
Water Resources Engineering       -------      ---------     --------    
Coastal and Fluvial Hydrodynamics --------------------------------------
------------------------------------------------------------------------