dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday 07 December 2009 09:32:57 Anders Logg wrote:
> On Mon, Dec 07, 2009 at 05:26:17PM +0000, Harish Narayanan wrote:
> > Johan Hake wrote:
> > > On Monday 07 December 2009 04:11:28 Garth N. Wells wrote:
> > >> Johan Hake wrote:
> > >> > On Sunday 06 December 2009 07:23:27 Garth N. Wells wrote:
> > >> >> Anders Logg wrote:
> > >> >>> On Sun, Dec 06, 2009 at 12:46:53PM +0000, Garth N. Wells wrote:
> > >> >>>> Anders Logg wrote:
> > >> >>>>> Shouldn't this be templated? This could be used also for passing
> > >
> > > uint
> > >
> > >> >>>>> arrays to and from Python.
> > >> >>>>
> > >> >>>> It probably should be, but it may make things more complicated on
> > >> >>>> the
> > >> >>>>
> > >> >>>> SWIG side. Perhaps Johan H. could comment - not on whether it's
> > >> >>>>
> > >> >>>> possible (I'm sure for Johan that it is), but on whether or not
> > >> >>>> it's
> > >> >>>>
> > >> >>>> easy.
> > >> >>>>
> > >> >>>>
> > >> >>>>
> > >> >>>> Templating it would be better on the C++ side, not just for
> > >> >>>> double,
> > >> >>>>
> > >> >>>> int, etc, but to handle const pointers correctly.
> > >> >>>
> > >> >>> ok.
> > >> >>>
> > >> >>>>> I imagine we can typedef
> > >> >>>>>
> > >> >>>>>
> > >> >>>>>
> > >> >>>>> DoubleArray
> > >> >>>>>
> > >> >>>>> UintArray
> > >> >>>>
> > >> >>>> This is a bit ugly. Might as well just use Array<double>,
> > >
> > > Array<uint>,
> > >
> > >> >>>> etc.
> > >> >>>
> > >> >>> I agree it's ugly. I didn't suggest it for use in C++ but that it
> > >> >>>
> > >> >>> might simplify for the wrapper generators but when I think of it
> > >> >>>
> > >> >>> again, it probably makes no difference for SWIG.
> > >> >>
> > >> >> As a start, we can declare the templates with these names in SWIG,
> > >> >> and
> > >> >>
> > >> >> SWIG will wrap the member functions and it should work out of the
> > >> >> box
> > >> >>
> > >> >> (with a few tweaks for [] and =).
> > >> >
> > >> > We can make the Array class templated. We can also specialize SWIG
> > >> >
> > >> > interfaces for the a particular Array<Foo> class dependent on what
> > >> > Foo
> > >> >
> > >> > is. If it is a primitive (double, uint) we can expand the class with
> > >> > a
> > >> >
> > >> > set of constructors, utility methods, and typemaps which should make
> > >> > it
> > >> >
> > >> > easy to interact with NumPy. However we need to cover all possible
> > >> > cases
> > >> >
> > >> > of such extensions wrt to memory handling and difference in intent
> > >> > use.
> > >>
> > >> Can you give us a short tutorial on NumPy and memory management,
> > >> namely:
> > >>
> > >>
> > >>
> > >> 1. If a Numpy object is created from an point to a double, will NumPy
> > >>
> > >> ever destroy the data?
> > >
> > > The long story:
> > >
> > > Here are lines 487-505 in ndarrayobject.h from NumPy in karmic. The
> > > lines describes the quite straightforward Capi of a PyArrayObject:
> > >
> > > typedef struct PyArrayObject {
> > >
> > > PyObject_HEAD
> > >
> > > char *data; /* pointer to raw data buffer */
> > >
> > > int nd; /* number of dimensions, also called ndim */
> > >
> > > npy_intp *dimensions; /* size in each dimension */
> > >
> > > npy_intp *strides; /* bytes to jump to get to the
> > >
> > > next element in each dimension */
> > >
> > > PyObject *base; /* This object should be decref'd
> > >
> > > upon deletion of array */
> > >
> > > /* For views it points to the original array */
> > >
> > > /* For creation from buffer object it points
> > >
> > > to an object that shold be decref'd on
> > >
> > > deletion */
> > >
> > > /* For UPDATEIFCOPY flag this is an array
> > >
> > > to-be-updated upon deletion of this one */
> > >
> > > PyArray_Descr *descr; /* Pointer to type structure */
> > >
> > > int flags; /* Flags describing array -- see below*/
> > >
> > > PyObject *weakreflist; /* For weakreferences */
> > >
> > > } PyArrayObject;
> > >
> > > We are mostly interested in data, nd, dimensions, descr and flags. The
> > > first three I think are self explanatory, descr describes the data type
> > > (strides tell us the size in bytes between the different data segments
> > > for each dimension. Interesting enough but I have never used it because
> > > macros handles low level iterations). flags holds a number of important
> > > flags for the state of the data. These are explained in the same file:
> > > lines 530-584. The one we are most interested in is:
> > >
> > > /* If set, the array owns the data: it will be free'd when the array
> > >
> > > is deleted. */
> > >
> > > #define NPY_OWNDATA 0x0004
> > >
> > > You can also print the most important flags in Python:
> > >
> > > z = Zeros(10)
> > >
> > > z.flags
> > >
> > > <http://docs.scipy.org/doc/numpy/reference/c-api.array.html#creating-ar
> > >rays>
> > >
> > > describes a bunch of macros for creating a NumPy array. If one use:
> > >
> > > PyArray_SimpleNewFromData
> > >
> > > which we use in for example mesh_pre.i for coordinates, the flag will
> > > be set to not owe the data. We can can use the more elaborated macro:
> > >
> > > PyArray_NewFromDescr
> > >
> > > to both pass a data structure and set the flags, to for example to owe
> > > the data. I suppose one can operate directly on an NumPy object to set
> > > the flags, but I would consider this dangerous.
> > >
> > >> 2. If the pointer to the NumyPy data is used to create say an Array
> > >>
> > >> object, can the Numpy object be stopped from destroying it, or do we
> > >>
> > >> need to make a copy?
> > >
> > > See above.
> > >
> > >> > In particular will we need to decide what we want the argument:
> > >> >
> > >> > Array<Foo> & to do. If we know a method will not resize the
> > >> > argument, we
> > >> >
> > >> > can just map the content to the NumPy pointer.
> > >>
> > >> I think that we should allow resizing.
> > >
> > > Me too, but I suppose we will not "allow" it for all methods. For
> > > example we do not want to let values change in the eval method? If we
> > > have such convention for some methods we could create a typemap that
> > > just pass the pointer from the NumPy array to the Array<Foo>. If we are
> > > not sure whether the data will be resized we cannot pass the pointer
> > > from the NumPy array.
> > >
> > >> > If we allow resize we can either turn the
> > >> >
> > >> > argument into a pure out argument, i.e., we construct the Array<Foo>
> > >> > in
> > >> >
> > >> > the interface, pass it to the C++ method, and when it is returned,
> > >
> > > create
> > >
> > >> > a NumPy array which takes ownership of the data and pass this to
> > >> > Python.
> > >>
> > >> This sounds easy.
> > >
> > > Easy enough.
> > >
> > >> Or we can let
> > >>
> > >> > the user pass an instance of ArrayFoo to the method, and use
> > >> >
> > >> > ArrayFoo.array() to get a NumPy view of the data.
> > >>
> > >> This would involve extending Array with the member function 'array' in
> > >>
> > >> the SWIG code?
> > >
> > > yes.
> > >
> > >> Would the function array create a NumPy view of the data, or would the
> > >>
> > >> NumPy object already exist (i.e., be created in the interface when
> > >>
> > >> moving the object from C++ to Python)?
> > >
> > > Depending on how we do it in the array method it can be a view.
> > >
> > >> This option sounds more complicated, but could be nice if it works
> > >> such
> > >>
> > >> that if a NumPy object is not required, it is never constructed.
> > >
> > > Exactly. I think we should add it just as an improvement of the Python
> > > interface of Array. See comment on copy versus view below.
> > >
> > >> > Note that we can have both of these typemaps (for the same Foo) at
> > >> > the
> > >> >
> > >> > same time in PyDOLFIN, by either using specific argument names for
> > >> > the
> > >> >
> > >> > two different cases:
> > >> >
> > >> >
> > >> >
> > >> > Array<Foo> & values
> > >> >
> > >> > Array<Foo> & values_resize
> > >> >
> > >> >
> > >> >
> > >> > or we need to construct module specific typemaps, which might get
> > >> > out of
> > >> >
> > >> > hands, when things get complicated.
> > >> >
> > >> >
> > >> >
> > >> > Garth do you suggest we also use Array<const Foo*> to pass what we
> > >> > now
> > >> >
> > >> > use std::vector for (DirichleBC, GenericFunction, aso)? Would this
> > >> > work
> > >> >
> > >> > from the C++ side? We loose the push_back functionality from
> > >
> > > std::vector.
> > >
> > >> > If you think we go for Array<const Foo*> it is no big deal creating
> > >> >
> > >> > specific typemaps for this.
> > >>
> > >> I'm not sure.
> > >
> > > Ok.
> > >
> > >> > As soon as I am able I can add NumPy typemaps for the most
> > >> >
> > >> > straightforward cases.
> > >>
> > >> Great. We can tweak the Array class as needed. I expect that we'll
> > >> need
> > >>
> > >> some functions to deal with memory management in the wrapper code.
> > >
> > > Maybe, or we could just control the deeper memory handling in the
> > > typemaps, preventing a user to segfault his program just because he
> > > missused the Python interface of ArrayFoo.
> > >
> > > Using a view of data is of course dangerous as it is. Try:
> > >
> > > mesh = UnitSquare(1,1)
> > >
> > > c = mesh.coordinates()
> > >
> > > del mesh
> > >
> > > Now will c's data variable will point to a deleted memory. If we want
> > > to avoid this we need to copy the data.
> >
> > EEEEEeek! KMail HTML mail.
> >
> > (With high-quality markup!)
> 
> +1

To all mutt users: 

  Move on!

Johan