numpy-discussion

[Numpy-discussion] Please help - pointer to slice

[Sebastian H. <ha...@ms...>]

Hi,
I am using SWIG to wrap some C++ function operating on numarray arrays.
Everything works quite well -- up until yesterday...

I have some 3d data and was thinking if the actual operation I want to do is
done "slice by slice", I could just loop over Z in python and have my C
function
defined like this:
void DoSomething(float * array2d, int nx, int ny)

The python loop looks like this:
arr = numarray.array(shape=(100,100,100), type=numarray.Float32)
for z in range(nz):
     DoSomething( arr[z] )

SOMEHOW the pointer that get transfered to the C++ side is always the same
as
for the  full (3d) array.

if I try to debug this by doing:
for z in range(nz):
    print repr(arr[z]._data)

that also tells be that python / numarray thinks all slices are located at
the same address, that is:
every slice looks just like the full array.

Please help,  I don't know what to do...

Sebastian Haase

Re: [Numpy-discussion] Please help - pointer to slice

[Francesc A. <fa...@op...>]

A Divendres 07 Mar=E7 2003 01:42, Sebastian Haase va escriure:
> if I try to debug this by doing:
> for z in range(nz):
>     print repr(arr[z]._data)
>
> that also tells be that python / numarray thinks all slices are located=
 at
> the same address, that is:
> every slice looks just like the full array.

You can access the slices by taking into account the byte offset
attribute. Look at the next example:

In [20]: a=3Dnumarray.arange(100, shape=3D(10,10))

In [21]: a._data
Out[21]: <memory at 083e58e8 with size:400 held by object at 083e58d0>

In [22]: a[1]._data
Out[22]: <memory at 083e58e8 with size:400 held by object at 083e58d0>

as you already know, both memory buffers point to the same address. I gue=
ss
this is implemented in that way so as to not copy data unnecessarily.

now, look at:

In [23]: a._byteoffset
Out[23]: 0

In [24]: a[1]._byteoffset
Out[24]: 40

So, you can know where the data actually starts by looking at the
_byteoffset property. In fact, you should always look at it in order to g=
et
proper results!.

In C, you can access to this information by looking at the byteoffset fie=
ld
in the PyArrayObject struct (look at chapter 10 in the User's Manual).

Hope that helps,

--=20
Francesc Alted

Re: [Numpy-discussion] Please help - pointer to slice

[Sebastian H. <ha...@ms...>]

Thanks, SO MUCH ...
I almost went crazy yesterday - and was really desperate by the time I wrote
that email.
Somehow I newer needed that offset field until now.  So: when I do the
NA_InputArray call  I get a "proper" C-array
just that it does NOT necessarily start at  NAimg->data
but rather at
NAimg->data + NAimg->byteoffset

Again: thanks so much.
BTW: Is there general interest in my SWIG typemaps. (SWIG is maybe the
easiest way to wrap C/C++ functions (and classes)
into Python (and/or Perl, Java, Ruby,...) ?  I think especially for
numerical stuff, that "link" in of interest.

Sebastian


----- Original Message -----
From: "Francesc Alted" <fa...@op...>
To: "Sebastian Haase" <ha...@ms...>;
<Num...@li...>
Sent: Thursday, March 06, 2003 9:40 PM
Subject: Re: [Numpy-discussion] Please help - pointer to slice


A Divendres 07 Març 2003 01:42, Sebastian Haase va escriure:
> if I try to debug this by doing:
> for z in range(nz):
>     print repr(arr[z]._data)
>
> that also tells be that python / numarray thinks all slices are located at
> the same address, that is:
> every slice looks just like the full array.

You can access the slices by taking into account the byte offset
attribute. Look at the next example:

In [20]: a=numarray.arange(100, shape=(10,10))

In [21]: a._data
Out[21]: <memory at 083e58e8 with size:400 held by object at 083e58d0>

In [22]: a[1]._data
Out[22]: <memory at 083e58e8 with size:400 held by object at 083e58d0>

as you already know, both memory buffers point to the same address. I guess
this is implemented in that way so as to not copy data unnecessarily.

now, look at:

In [23]: a._byteoffset
Out[23]: 0

In [24]: a[1]._byteoffset
Out[24]: 40

So, you can know where the data actually starts by looking at the
_byteoffset property. In fact, you should always look at it in order to get
proper results!.

In C, you can access to this information by looking at the byteoffset field
in the PyArrayObject struct (look at chapter 10 in the User's Manual).

Hope that helps,

--
Francesc Alted

Re: [Numpy-discussion] Please help - pointer to slice

[Fernando P. <fp...@co...>]

Sebastian Haase wrote:
> Thanks, SO MUCH ...
> I almost went crazy yesterday - and was really desperate by the time I wrote
> that email.
> Somehow I newer needed that offset field until now.  So: when I do the
> NA_InputArray call  I get a "proper" C-array
> just that it does NOT necessarily start at  NAimg->data
> but rather at
> NAimg->data + NAimg->byteoffset

Well, in case you (or others) find it useful, I'm including here a little 
library I wrote for accessing general Numeric 2-d arrays (contiguous or not) 
in an easy manner.

Here's a snippet of a simple (included) example of a function to print an 
integer array:

static PyObject *idisp(PyObject *self, PyObject *args)
{
     PyArrayObject *array;
     int **arr;  // for the data area
     int i,j,cs;

     if (!PyArg_ParseTuple(args, "O!",&PyArray_Type,&array ) )
	return NULL;

     arr = imatrix_data(array,&cs);
     for (i=0;i<array->dimensions[0];++i) {
       for (j=0;j<cs*array->dimensions[1];j+=cs)
	printf("%5d ",arr[i][j]);
       printf("\n");
     }
     free(arr);

     Py_INCREF(Py_None);
     return Py_None;
}

You get the **arr pointer and you can then manipulate it as a[i][j] 
conveniently.  The supplied example file may be enough for many to write their 
Numpy C extensions without much trouble.

> Again: thanks so much.
> BTW: Is there general interest in my SWIG typemaps. (SWIG is maybe the
> easiest way to wrap C/C++ functions (and classes)
> into Python (and/or Perl, Java, Ruby,...) ?  I think especially for
> numerical stuff, that "link" in of interest.

I'd love to see them.  So far I've either used high-level stuff like 
weave.inline() or just written the extensions by hand (as in the code I'm 
supplying here).  I'd like to see this, especially if there is an easy way to 
handle contiguity issues with it.

Best,

f.

Re: [Numpy-discussion] Please help - pointer to slice

[Sebastian H. <ha...@ms...>]

<snip>

> Well, in case you (or others) find it useful, I'm including here a little
> library I wrote for accessing general Numeric 2-d arrays (contiguous or
not)
> in an easy manner.
>
> Here's a snippet of a simple (included) example of a function to print an
> integer array:
>
> static PyObject *idisp(PyObject *self, PyObject *args)
> {
>      PyArrayObject *array;
>      int **arr;  // for the data area
>      int i,j,cs;
>
>      if (!PyArg_ParseTuple(args, "O!",&PyArray_Type,&array ) )
> return NULL;
>
>      arr = imatrix_data(array,&cs);
>      for (i=0;i<array->dimensions[0];++i) {
>        for (j=0;j<cs*array->dimensions[1];j+=cs)
> printf("%5d ",arr[i][j]);
>        printf("\n");
>      }
>      free(arr);
>
>      Py_INCREF(Py_None);
>      return Py_None;
> }
>
> You get the **arr pointer and you can then manipulate it as a[i][j]
> conveniently.  The supplied example file may be enough for many to write
their
> Numpy C extensions without much trouble.
>
> > Again: thanks so much.
> > BTW: Is there general interest in my SWIG typemaps. (SWIG is maybe the
> > easiest way to wrap C/C++ functions (and classes)
> > into Python (and/or Perl, Java, Ruby,...) ?  I think especially for
> > numerical stuff, that "link" in of interest.
>
> I'd love to see them.  So far I've either used high-level stuff like
> weave.inline() or just written the extensions by hand (as in the code I'm
> supplying here).  I'd like to see this, especially if there is an easy way
to
> handle contiguity issues with it.

You are using Numeric  not numarray , right?

In any case :
My thinking goes as follows:
I was looking for a wrapping-solution that would be most transparent to the
people that would write the C/C++ code (or even Fortran, BTW)
Since I already had some experience using SWIG, that's what I wanted to use
to handle the "numarray binding".
SWIG has a quite strong (meaning: flexible, general) way of handling data
types. They call it "type maps":
Once you have that in place (later more) the whole "binding" looks like
this: (an example)

C side:
 double sebFunc(float *arr, int nx, int ny, int nz) { double a =0; for(int
i=0;i<nx*ny*nz;i++) a+= arr[i]; return a; }

Then there is the SWIG interface file (suggested file extension ".i"):
double sebFunc(float *array3d, int nx, int ny, int nz);

((This has (in firts approx.) the same syntax as a normal C-header file -
just that SWIG "bites" on the argument-variable names,
in this case (float *array3d, int nx, int ny, int nz)  and realizes (because
of my "type map") that this should be a python-(numeric) array

Then you call it from python just like that...
import myModule
print myModule.sebFunc(array)

So the magic is all hidden in the the typemap:
For that I have a separate SWIG-interface that gets #included into the above
mentioned one.
Here is  that interface file ( just the part for 3d float arrays)

%typecheck(SWIG_TYPECHECK_FLOAT) (float *array2d, int nx, int ny) {
  if(!PyArray_Check($input)) $1=0;
  else if(((PyArrayObject*)$input)->descr->type_num != tFloat32) $1=0;
  else if(((PyArrayObject*)$input)->nd != 2) $1=0;
  else $1=1;
}

%typemap(python,in)
  (float *array3d, int nx, int ny, int nz)
  (PyArrayObject *temp=NULL)
{
  debugPrintf("debug: float *array3d  -> NA_InputArray\n");
  PyArrayObject *NAimg  = NA_InputArray($input, tFloat32, NUM_C_ARRAY);

  if (!NAimg) {
 printf("**** no (float) numarray *****\n");
 return 0;
  }
  temp = NAimg;

  $1 = (float *) (NAimg->data + NAimg->byteoffset);
  switch(NAimg->nd) {
  case 1:
    $2 = NAimg->dimensions[0];
    $3=1;
    break;
  case 2:
    $2 = NAimg->dimensions[1];
    $3=NAimg->dimensions[0];
    break;
  default:
    debugPrintf(" **** numarray dim >2 (ns=%d)\n", NAimg->nd);
    _SWIG_exception(SWIG_RuntimeError, "numarray dim > 2");
    return 0;
}

%typemap(python,freearg)
  (float *array3d, int nx, int ny, int nz)
{
  Py_XDECREF(temp$argnum);
}


The first part (%typecheck) is needed so that SWIG call even handle
overloaded functioned correctly. (BTW, SWIG happily wraps my template
functions to; just with classes (that keep a reference to an array) I not
that sure yet [reference counting])

I think I don't have all necessary error handling parts implemented yet.

But this is already a VERY useful tool for me
- and I am in fact working right now on convincing some ( C and Fortran
only) people to try this  ;-)
 [ once they see how cool it is to call their (Fortran or C) directly from
Python, maybe they start to realize that there is "something new out there"
(emm, I meant Python ]


So what is the status on this list: Are people familiar with (know and/or
using) SWIG ?

Cheers,
Sebastian

Re: [Numpy-discussion] Please help - pointer to slice

[Sebastian H. <ha...@ms...>]

> > You are using Numeric  not numarray , right?
>
> Correct.  I need this stuff for production use, so I just can't quite move
to
> numarray yet.
>
> > In any case :
> > My thinking goes as follows:
> > I was looking for a wrapping-solution that would be most transparent to
the
> > people that would write the C/C++ code (or even Fortran, BTW)
> > Since I already had some experience using SWIG, that's what I wanted to
use
> > to handle the "numarray binding".
> > SWIG has a quite strong (meaning: flexible, general) way of handling
data
> > types. They call it "type maps":
> > Once you have that in place (later more) the whole "binding" looks like
> > this: (an example)
>
> [Snip]
>
> Does your stuff handle non-contiguity transparently?  I couldn't quite
find a
> solution for that, which means that a **pointer I get can't be blindly
passed
> to -say- a Numerical Recipes routine which expects to treat it like a
matrix.
>   If it's contiguous all is fine, but otherwise one needs to either make a
> copy or handle the 'skipping' by hand in the column-index loops (as shown
in
> my examples).
>
> I tried to come up with a solution in C, but managed to convince myself
that
> it just can't be done.  I was thinking of some typecasting black magic
where
> the pointers were recast for the pointer arithmetic to hide the column
> offsets, with the proper type still being returned at dereferencing time.
I
> finally decided this was just impossible in C because the compiler wants
to
> know the pointer sizes at compile time.  But maybe I'm wrong and just
didn't
> think about the problem hard enough.
>
> At any rate, I'd be interested in seeing your typemap stuff.  Especially
if it
> either works with Numeric or it can be modified for Numeric use with
> reasonably limited effort.
>
Fernando,
Sorry, I forgot to comment on the "non-contiguity " question...
First: My last post contained practically everything there is to typemapping
needed when you use SWIG. ("practically" means that I have many versions
of that for handling other types (like short, int, complex, ...) and fro
handling 2d
arrays.) Do YOU know SWIG ?

"handle non-contiguity transparently"
What are exactly the possible scenarios of "non-contiguous" ?
If, for example having a 3d array, only going along the first axis is non-c.
but along the other ones still is, I could imagine that it would not be to
problematic to "forward" that scheme into the C side, that is: in a
non-transparent
way.
The only alternative I see would be introducing some macros to translate
each memory
pointer. But that would also not be transparent ( since the C side would
have to use
these "wired" macros) and I is like to be slower, because of all the new
pointer arithmetic (inside those macros)
I would favor the the first scenario.

Regards,
Sebastian

Re: [Numpy-discussion] Please help - pointer to slice

[Todd M. <jm...@st...>]

Sebastian Haase wrote:

>Thanks, SO MUCH ...
>I almost went crazy yesterday - and was really desperate by the time=
 I wrote
>that email.
>Somehow I newer needed that offset field until now.  So: when I do t=
he
>NA_InputArray call  I get a "proper" C-array
>just that it does NOT necessarily start at  NAimg->data
>but rather at
>NAimg->data + NAimg->byteoffset
> =20
>
The numarray-0.4 manual (available at:=20
http://prdownloads.sourceforge.net/numpy/numarray-0.4.pdf?download)=
=20
documents how to write Python stubs using the "high level" API on Pag=
e=20
70.  The thing you appear to be missing is a call to NA_OFFSETDATA wh=
ich=20
gets a pointer to the data in an array by adding the buffer pointer a=
nd=20
byteoffset together.

Todd

>Again: thanks so much.
>BTW: Is there general interest in my SWIG typemaps. (SWIG is maybe t=
he
>easiest way to wrap C/C++ functions (and classes)
>into Python (and/or Perl, Java, Ruby,...) ?  I think especially for
>numerical stuff, that "link" in of interest.
>
>Sebastian
>
>
>----- Original Message -----
>From: "Francesc Alted" <fa...@op...>
>To: "Sebastian Haase" <ha...@ms...>;
><Num...@li...>
>Sent: Thursday, March 06, 2003 9:40 PM
>Subject: Re: [Numpy-discussion] Please help - pointer to slice
>
>
>A Divendres 07 Mar=E7 2003 01:42, Sebastian Haase va escriure:
> =20
>
>>if I try to debug this by doing:
>>for z in range(nz):
>>    print repr(arr[z]._data)
>>
>>that also tells be that python / numarray thinks all slices are loc=
ated at
>>the same address, that is:
>>every slice looks just like the full array.
>>   =20
>>
>
>You can access the slices by taking into account the byte offset
>attribute. Look at the next example:
>
>In [20]: a=3Dnumarray.arange(100, shape=3D(10,10))
>
>In [21]: a._data
>Out[21]: <memory at 083e58e8 with size:400 held by object at 083e58d=
0>
>
>In [22]: a[1]._data
>Out[22]: <memory at 083e58e8 with size:400 held by object at 083e58d=
0>
>
>as you already know, both memory buffers point to the same address. =
I guess
>this is implemented in that way so as to not copy data unnecessarily=
.
>
>now, look at:
>
>In [23]: a._byteoffset
>Out[23]: 0
>
>In [24]: a[1]._byteoffset
>Out[24]: 40
>
>So, you can know where the data actually starts by looking at the
>_byteoffset property. In fact, you should always look at it in order=
 to get
>proper results!.
>
>In C, you can access to this information by looking at the byteoffse=
t field
>in the PyArrayObject struct (look at chapter 10 in the User's Manual=
).
>
>Hope that helps,
>
>--
>Francesc Alted
>
>
>
>
>
>-------------------------------------------------------
>This SF.net email is sponsored by: Etnus, makers of TotalView, The d=
ebugger=20
>for complex code. Debugging C/C++ programs can leave you feeling los=
t and=20
>disoriented. TotalView can help you find your way. Available on majo=
r UNIX=20
>and Linux platforms. Try it free. www.etnus.com
>_______________________________________________
>Numpy-discussion mailing list
>Num...@li...
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> =20
>