numpy-discussion — Discussion list for all users of Numerical Python

[Numpy-discussion] Re: Re-implementation of Python Numerical arrays (Numeric) available for download

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

Perry Greenfield wrote:

> > One major comment that isn't directly addressed on the web page is the
> > ease of writing new functions, I suppose Ufuncs, although I don't
> > usually care if they work on anything other than Arrays. I hope the new
> > system will make it easier to write new ones. 
<snip>

> Absolutely. We will provide examples of how to write new ufuncs. It should
> be very simple in one sense (requiring few lines of code) if our code
> generator machinery is used (but context is important here so this
> is why examples or a template is extremely important). But it isn't
> particularly hard to do without the code generator. And such ufuncs
> will handle *all* the generality of arrays including slices, non-aligned
> arrays, byteswapped arrays, and type conversion. I'd like to provide
> examples of writing ufuncs within a few weeks (along with examples
> of other kinds of functions using the C-API as well).

This sounds great! The code generting machinery sound very promising,
and examples are, of course, key. I found digging through the NumPy
source to figure out how to do things very treacherous. Making writing
Ufuncs easy will enocourage a lot more C Ufuncs to be written which
should help perfomance.

> > Also, I can't help wondering if this could leverage more existing code.
> > The blitz++ package being used by Eric Jones in the SciPy.compiler
> > project looks very promising. It's probably too late, but I'm wondering
> > what the reasons are for re-inventing such a general purpose wheel.
> >
> I'm not sure which "wheel" you are talking about :-)

The wheel I'm talking about are multi-dimensional array objects...

> We certainly
> aren't trying to replilcate what Eric Jones has done with the
> SciPy.compiler approach (which is very interesting in its own right).

I know, I just think using an existing set of C++ classes for multiple
typed multidimansional arrays would make sense, although I imagine it is
too late now!

> If the issue is why we are redoing Numeric:

Actually, I think I had a pretty good idea why you were working on this.

> 1) it has to be rewritten to be acceptable to Guido before it can be
>    part of the Standard Library.
> 2) to add new types (e.g. unsigned) and representations (e.g., non-aligned,
>    byteswapped, odd strides, etc). Using memory mapped data requires some
>    of these.
> 3) to make it more memory efficient with large arrays.
> 4) to make it more generally extensible

I'm particualry excited about 1) and 4)

> > As a whole I have found that I would like the transition from Python to
> > Compiled laguages to be smoother. The standard answer to Python
> > perfomance is to profile, and then re-write the computationally intesive
> > pertions in C. This would be a whole lot easier if Python used datatypes
> > that are easy to use from C/C++ as well as Python. I hope NumPy2 can
> > move in this direction.
> >
> What do you see as missing in numarray in that sense? Aside from UInt32
> I'm not aware of any missing type that is available on all platforms.
> There is the issue of Float128 and such. Adding these is not hard.
> The real issue is how to deal with the platforms that don't support them.

I used Poor wording. When I wrote "datatypes", I meant data types in a
much higher order sense. Perhaps structures or classes would be a better
term. What I mean is that is should be easy to use an manipulate the
same multidimensional arrays from both Python and C/C++. In the current
Numeric, most folks generate a contiguous array, and then just use the
array->data pointer to get what is essentially a C array. That's fine if
you are using it in a traditional C way, with fixed dimension, one
datatype, etc. What I'm imagining is having an object in C or C++ that
could be easily used as a multidimentional array. I'm thinking C++ would
probably neccesary, and probably templates as well, which is why blitz++
looked promising. Of course, blitz++ only compiles with a few up-to-date
compilers, so you'd never get it into the standard library that way!

This could also lead the way to being able to compile NumPy code....<end
fantasy>

> I think it is pretty easy to install since it use distutils.

I agree, but from the newsgroup, it is clear that a lot of folks are
very reluctant to use something that is not part of the standard
library.

> > >    We estimate
> > >    that numarray is probably another order of magnitude worse,
> > >    i.e., that 20K element arrays are at half the asymptotic
> > >    speed. How much should this be improved?
> >
> > A lot. I use arrays smaller than that most of the time!
> >
> What is good enough? As fast as current Numeric?

As fast as current Numeric would be "good enough" for me. It would be a
shame to go backwards in performance!

> (IDL does much
> better than that for example).

My personal benchmark is MATLAB, which I imagine is similar to IDL in
performance.

> 10 element arrays will never be
> close to C speed in any array based language embedded in an
> interpreted environment.

Well, sure, I'm not expecting that

> 100, maybe, but will be very hard.
> 1000 should be possible with some work.

I suppose MATLAB has it easier, as all arrays are doubles, and, (untill
recently anyway), all variable where arrays, and all arrays were 2-d.
NumPy is a lot more flexible that that. Is is the type and size checking
that takes the time?
 
> Another approach is to try to cast many of the functions as being
> able to broadcast over repeated small arrays. After all, if one
> is only doing a computation on one small array, it seems unlikely
> that the overhead of Python will be objectionable. Only if you
> have many such arrays to repeat calculations on, should it be
> a problem (or am I wrong about that).

You are probably right about that.

> If these repeated calculations
> can be "assembled"  into a higher dimensionality array (which
> I understand isn't always possible) and operated on in that sense,
> the efficiency issue can be dealt with.

I do that when possible, but it's not always possible.

> But I guess this can only
> be seen with specific existing examples and programs. I would
> be interested in seeing the kinds of applications you have now
> to gauge what the most effective solution would be.

One of the things I do a lot with are coordinates of points and
polygons. Sets if points I can handle easily as an NX2 array, but
polygons don't work so well, as each polgon has a different number of
points, so I use a list of arrays, which I have to loop over. Each
polygon can have from about 10 to thousands of points (mostly 10-20,
however). One way I have dealt with this is to store a polygon set as a
large array of all the points, and another array with the indexes of the
start and end of each polygon. That way I can transform the coordinates
of all the polygons in one operation. It works OK, but sometimes it is
more useful to have them in a sequence. 

> As mentioned,
> we tend to deal with large data sets and so I don't think we have
> a lot of such examples ourselves.

I know large datasets were one of your driving factors, but I really
don't want to make performance on smaller datasets secondary.

I hope I'll get a chance to play with it soon....

-Chris


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

[Numpy-discussion] Re: Re-implementation of Python Numerical arrays (Numeric) available for download

[Perry G. <pe...@st...>]

>   6) Should array properties be accessible as public attributes
>     instead of through accessor methods?
> 
>     We don't currently allow public array attributes to make
>     the Python code simpler and faster (otherwise we will
>     be forced to use __setattr__ and such). This results in
>     incompatibilty with previous code that uses such attributes.
> 
> 
> I prefer the use of public attributes over accessor methods.
> 
> 
> -- 
> Paul Barrett, PhD      Space Telescope Science Institute

The issue of efficiency may not be a problem with Python 2.2
or later since it provides new mechanisms that avoid the need
to use __setattr__ to solve this problem. (e.g. __slots__,
property, __get__, and __set__). So it becomes more
of an issue of which style people prefer rather than simplicity
and speed of the code.

Perry

Re: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Paul B. <Ba...@st...>]

Perry Greenfield wrote:

> 
> An early beta version is available on sourceforge as the
> package Numarray (http://sourceforge.net/projects/numpy/)
> 
> Information on the goals, changes in user interface, open issues,
> and design can be found at http://aten.stsci.edu/numarray


  6) Should array properties be accessible as public attributes
    instead of through accessor methods?

    We don't currently allow public array attributes to make
    the Python code simpler and faster (otherwise we will
    be forced to use __setattr__ and such). This results in
    incompatibilty with previous code that uses such attributes.


I prefer the use of public attributes over accessor methods.


-- 
Paul Barrett, PhD      Space Telescope Science Institute
Phone: 410-338-4475    ESS/Science Software Group
FAX:   410-338-4767    Baltimore, MD 21218

Re: [Numpy-discussion] Heigenvalues is broken

[Scott R. <ra...@ph...>]

On November 19, 2001 02:36 pm, Jeff Whitaker wrote:
>
> This is definitely a hardware/compiler dependant feature.  I get the
> "right" answer on Solaris (with the forte compiler) but the same "wrong"
> answer as Alessandro on MacOS X/gcc.  I've tried fiddling with compiler
> options on my OS X box, to no avail.

But seemingly it is even stranger than this.  Here are my results from Debian 
unstable using Lapack 3.0 on an Athlon system:

Python 2.1.1 (#1, Nov 11 2001, 18:19:24)
[GCC 2.95.4 20011006 (Debian prerelease)] on linux2
Type "copyright", "credits" or "license" for more information.
>>> from LinearAlgebra import *
>>> a=array([[1,0],[0,1]])
>>> b=array([[0,1],[-1,0]])
>>> M=a+b*complex(0,1.0)
>>> Heigenvalues(M)
array([ 0.,  2.])

Scott


> On Sun, 18 Nov 2001, Travis Oliphant wrote:
> > On Sunday 18 November 2001 09:40 am, Alessandro Mirone wrote:
> > > Is it a problem of lapack3.0 of of
> > > LinearAlgebra.py?
> > > ..................... ==> (Eigenvalues should be (0,2))
> > >
> > > >>> a=array([[1,0],[0,1]])
> > > >>> b=array([[0,1],[-1,0]])
> > > >>> M=a+b*complex(0,1.0)
> > > >>> Heigenvalues(M)
> >
> > I suspect it is your lapack.  On an Athlon running Mandrake Linux with
> > the lapack-3.0-9 package, I get.
> >
> > >>> a=array([[1,0],[0,1]])
> > >>> b=array([[0,1],[-1,0]])
> > >>> M=a+b*complex(0,1.0)
> > >>> Heigenvalues(M)
> >
> > array([ 0.,  2.])
>
> This is definitely a hardware/compiler dependant feature.  I get the
> "right" answer on Solaris (with the forte compiler) but the same "wrong"
> answer as Alessandro on MacOS X/gcc.  I've tried fiddling with compiler
> options on my OS X box, to no avail.
>
> -Jeff
>
>  --
> Jeffrey S. Whitaker         Phone  : (303)497-6313
> Meteorologist               FAX    : (303)497-6449
> NOAA/OAR/CDC  R/CDC1        Email  : js...@cd...
> 325 Broadway                Web    : www.cdc.noaa.gov/~jsw
> Boulder, CO, USA 80303-3328 Office : Skaggs Research Cntr 1D-124
>
>
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

-- 
Scott M. Ransom              Address:  McGill Univ. Physics Dept.
Phone:  (514) 398-6492                 3600 University St., Rm 338
email:  ra...@ph...       Montreal, QC  Canada H3A 2T8 
GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989

Re: [Numpy-discussion] Heigenvalues is broken

[Jeff W. <js...@cd...>]

On Sun, 18 Nov 2001, Travis Oliphant wrote:

> On Sunday 18 November 2001 09:40 am, Alessandro Mirone wrote:
> > Is it a problem of lapack3.0 of of
> > LinearAlgebra.py?
> > ..................... ==> (Eigenvalues should be (0,2))
> >
> > >>> a=array([[1,0],[0,1]])
> > >>> b=array([[0,1],[-1,0]])
> > >>> M=a+b*complex(0,1.0)
> > >>> Heigenvalues(M)
>
> I suspect it is your lapack.  On an Athlon running Mandrake Linux with the
> lapack-3.0-9 package, I get.
>
> >>> a=array([[1,0],[0,1]])
> >>> b=array([[0,1],[-1,0]])
> >>> M=a+b*complex(0,1.0)
> >>> Heigenvalues(M)
> array([ 0.,  2.])

This is definitely a hardware/compiler dependant feature.  I get the
"right" answer on Solaris (with the forte compiler) but the same "wrong"
answer as Alessandro on MacOS X/gcc.  I've tried fiddling with compiler
options on my OS X box, to no avail.

-Jeff

 --
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/CDC  R/CDC1        Email  : js...@cd...
325 Broadway                Web    : www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80303-3328 Office : Skaggs Research Cntr 1D-124

[Numpy-discussion] Re: Numpy-discussion digest, Vol 1 #345 - 4 msgs

[Joe H. <jh...@oo...>]

Just to fill in the blanks, here's what IDL does:

IDL> a = [[1,2,3,4], $
IDL>      [5,6,7,8], $
IDL>      [9,10,11,12], $
IDL>      [13,14,15,16]]
IDL> print,a
       1       2       3       4
       5       6       7       8
       9      10      11      12
      13      14      15      16
IDL> print, a[[1,3],[0,3]] 
       2      16

--jh--

[Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Todd A. P. Ph.D. <tp...@ac...>]

Thanks for all of your work.  Things seem to be shaping up nicely.  I
just wanted to second some of the concerns below:

> Complex Types:
> ==============
> 
> 1)  I don't like the idea of complex types being a separate subclass of 
> ndarray.  This makes them "different."   Unless this "difference" can be 
> completely hidden (which I doubt), I would prefer complex types to be on the 
> same level as other numeric types.
> 
> 2)  Also,  in your C-API, you have a different pointer to the imaginary data. 
>   I much prefer the way it is done currently to have complex numbers 
> represented as an 8-byte, or 16-byte chunk of contiguous memory. 
> 

The second comment above is really critical for accessing utility
available in a very large number of numerical libraries.  In my view
this would "break" the utility of numpy severely -- recopying arrays
both on the way out and the way in would be extremely cumbersome.

-Todd Alan Pitts

Re: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Peter V. <Pet...@em...>]

On Saturday 17 November 2001 00:11 am, you wrote:
> note that my main use of numpy is as a pixel buffer for images. some of
> the changes like avoiding type promotion sounds really good to me :]

I have exactly the same application so I agree with this.

> 7) necessary to add other types?
> yes. i really want unsigned int16 and unsigned int32. all my operations
> are on pixel data, and things can just get messy when i need to treat
> packed color values as signed integers.

Yes please! One of the things that irritates me most on the original Numeric 
is that some types are lacking. I think the whole range of data types should 
be supported, even if some may be seldom used by most people.

> one other thing i'd like there to be a little focus on is adding my own
> new ufunc operators. for image manipulation i'd like new ufunc operators
> that clamp the results to legal values. i'd be happy to do this myself,
> but i don't believe it's possible with the current Numeric.

I write functions in C that directly access the numeric data. I don't use the 
ufunc api. One reason that I do that is that I want my libary of routines to 
be useful independent of  Numeric, so I only have a tiny glue between my C 
routines and Numeric. I hope that it will be still possible to do this in the 
new version.

> the last thing i really really want is for this to be rolled into the
> standard python distribution. that is perhaps the most important aspect
> for me. i do not like requiring the extra dependency for generic numeric
> arrays. :]

I second that!

Cheers, Peter
-- 
Dr. Peter J. Verveer
Bastiaens Group
Cell Biology and Cell Biophysics Programme
EMBL
Meyerhofstrasse 1
D-69117 Heidelberg
Germany
Tel. : +49 6221 387245
Fax  : +49 6221 387242
Email: Pet...@em...

Re: [Numpy-discussion] Heigenvalues is broken

[Nils W. <nw...@me...>]

Alessandro Mirone schrieb:
> 
> Is it a problem of lapack3.0 of of
> LinearAlgebra.py?
> ..................... ==> (Eigenvalues should be (0,2))
> 
> >>> a=array([[1,0],[0,1]])
> >>> b=array([[0,1],[-1,0]])
> >>> M=a+b*complex(0,1.0)
> >>> Heigenvalues(M)
> array([-2.30277564,  1.30277564])
> >>> print M
> [[ 1.+0.j  0.+1.j]
>  [ 0.-1.j  1.+0.j]]
> >>>
> 
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

On an Athlon running SuSE Linux 7.3 with the  
lapack-3.0-0 package, I get.

[-2.30277564  1.30277564]

Nils

Re: [Numpy-discussion] Heigenvalues is broken

[Travis O. <oli...@ie...>]

On Sunday 18 November 2001 09:40 am, Alessandro Mirone wrote:
> Is it a problem of lapack3.0 of of
> LinearAlgebra.py?
> ..................... ==> (Eigenvalues should be (0,2))
>
> >>> a=array([[1,0],[0,1]])
> >>> b=array([[0,1],[-1,0]])
> >>> M=a+b*complex(0,1.0)
> >>> Heigenvalues(M)

I suspect it is your lapack.  On an Athlon running Mandrake Linux with the  
lapack-3.0-9 package, I get.

>>> a=array([[1,0],[0,1]])
>>> b=array([[0,1],[-1,0]])
>>> M=a+b*complex(0,1.0)
>>> Heigenvalues(M)
array([ 0.,  2.])


-Travis

[Numpy-discussion] Heigenvalues is broken

[Alessandro M. <ale...@wa...>]

Is it a problem of lapack3.0 of of
LinearAlgebra.py?
..................... ==> (Eigenvalues should be (0,2)) 


>>> a=array([[1,0],[0,1]])
>>> b=array([[0,1],[-1,0]])
>>> M=a+b*complex(0,1.0)
>>> Heigenvalues(M)
array([-2.30277564,  1.30277564])
>>> print M
[[ 1.+0.j  0.+1.j]
 [ 0.-1.j  1.+0.j]]
>>>

Re: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Perry G. <gre...@ho...>]

From: Pete Shinners <pe...@sh...>

> 7) necessary to add other types?
> yes. i really want unsigned int16 and unsigned int32. all my operations 
> are on pixel data, and things can just get messy when i need to treat 
> packed color values as signed integers.
> 
Unsigned int16 is already supported. UInt32 could be done, but
raises some interesting issues with regard to combining with
Int32. I don't believe the current implementation prevents you
from carrying around unsigned data in Int32 arrays. If you
are using them as packed color values, do you ever do any
arithmetic operations on them other than to pack and unpack 
them?

> one other thing i'd like there to be a little focus on is adding my own 
> new ufunc operators. for image manipulation i'd like new ufunc operators 
> that clamp the results to legal values. i'd be happy to do this myself, 
> but i don't believe it's possible with the current Numeric.
> 
It will be possible for users to add their own ufuncs. We will 
eventually document how to do so (and it should be fairly simple
to do once we give a few example templates).

Perry
> 

RE: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Perry G. <gre...@ho...>]

> > I think that we also don't like that, and after doing the original,
> > somewhat incomplete, implementation using the subarray approach,
> > I began to feel that implementing it in C (albiet using a different
> > approach for the code generation) was probably easier and more
> > elegant than what was done here. So you are very likely to see
> > it integrated as a regular numeric type, with a more C-based
> > implementation.
>
> Sounds good.   Is development going to take place on the CVS
> tree.  If so, I
> could help out by comitting changes directly.
>
> >
> > > 2)  Also,  in your C-API, you have a different pointer to the
> > > imaginary data.
> > >   I much prefer the way it is done currently to have complex numbers
> > > represented as an 8-byte, or 16-byte chunk of contiguous memory.
> >
> > Any reason not to allow both? (The pointer to the real can be
> interpreted
> > as either a pointer to 8-byte or 16-byte quantities). It is true
> > that figuring out the imaginary pointer from the real is trivial
> > so I suppose it really isn't necessary.
>
> I guess the way you've structured the ndarray, it is possible.  I figured
> some operations might be faster, but perhaps not if you have two pointers
> running at the same time, anyway.
>
Well, the C implementation I was thinking of would only use
one pointer. The API could supply both if some algorithms would
find it useful to just access the imaginary data alone. But as
mentioned, I don't think it is important to include, so we
could easily get rid of it (and probably should)

> >
> > > Index Arrays:
> > > ===========
> > >
> > > 1)  For what it's worth, my initial reaction to your indexing
> scheme is
> > > negative.  I would prefer that if
> > >
> > > a = [[1,2,3,4],
> > >       [5,6,7,8],
> > >       [9,10,11,12],
> > >       [13,14,15,16]]
> > >
> > > then
> > >
> > > a[[1,3],[0,3]] returns the sub-matrix:
> > >
> > > [[   4,  6],
> > >  [ 12, 14]
> > >
> > > i.e. the cross-product of [1,3] x [0,3]   This is the way MATLAB
> > > works.  I'm
> > > not sure what IDL does.
> >
> > I'm afraid I don't understand the example. Could you elaborate
> > a bit more how this is supposed to work? (Or is it possible
> > there is an error? I would understand it if the result were
> > [[5, 8],[13,16]] corresponding to the index pairs
> > [[(1,0),(1,3)],[(3,0),(3,3)]])
> >
>
> The idea is to consider indexing with arrays of integers to be a
> generalization of slice index notation.   Simply interpret the
> slice as an
> array of integers that would be formed by using the range operator.
>
> For example, I would like to see
>
> a[1:5,1:3] be the same thing  as  a[[1,2,3,4],[1,2]]
>
> a[1:5,1:3] selects the 2-d subarray consisting of rows 1 to 4 and
> columns 1
> to 2 (inclusive starting with the first row being row 0).  In
> other words,
> the indices used to select the elements of a are ordered-pairs
> taken from the
> cross-product of the index set:
>
> [1,2,3,4] x [1,2] = [(1,1), (1,2), (2,1), (2,2), (3,1), (3,2),
> (4,1), (4,2)]
> and these selected elements are structured as a 2-d array of shape (4,2)
>
> Does this make more sense?  Indexing would be a natural extension of this
> behavior but allowing sets that can't be necessarily formed from
> the range
> function.
>
I understand this (but is the example in the first message
consistent with this?). This is certainly a reasonable
interpetation. But if this is the way multiple index arrays
are interpreted, how does one easily specify scattered points
in a multidimensional array? The only other alternative I can
think of is to use some of the dimensions of a multidimensional
index array as indicies for each of the dimensions. For example,
if one wanted to index random points in a 2d array, then
supplying an nx2 array would provide a list of n such points.
But I see this as a more limiting way to do this (and there
are often benefits to being able to keep the indices for
different dimensions in separate arrays.

But I think doing what you would like to do is straightforward
even with the existing implementation. For example, if x is a
2d array we could easily develop a function such that:

x[outer_index_product([1,3,4],[1,5])]
# with a better function name!

The function outer_index_product would return a tuple of two
index arrays each with a shape of 3x2. These arrays
would not take up more space than the original
arrays even though they appear to have a much
larger size (the one dimension is replicated by
use of a 0 stride size so the data buffer is
the same as the original). Would this be acceptable?

In the end, all these indexing behaviors can be provided
by different functions. So it isn't really a question of
which one to have and which not to have. The question is
what is supported by the indexing notation? For us, the
behavior we have implemented is far more useful for our
applications than the one you propose. But perhaps we are
in the minority, so I'd be very interested in hearing which
indexing interpretation is most useful to the general
community.

> > Why not:
> >
> > ravel(a)[[9,10,11]] ?
>
> sure, that would work, especially if ravel doesn't make a copy of
> the data
> (which I presume it does not).
>
Correct.

Perry

FW: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Perry G. <gre...@ho...>]

-----Original Message-----
>
> What I've seen looks great.   You've all done some good work here.
>
Thanks, you were origin of some of the ideas used.

> Of course, I do have some feedback.  I haven't looked at
> everything, these
> points have just caught my eye.
>
> Complex Types:
> ==============
>
> 1)  I don't like the idea of complex types being a separate subclass of
> ndarray.  This makes them "different."   Unless this "difference" can be
> completely hidden (which I doubt), I would prefer complex types
> to be on the
> same level as other numeric types.
>
I think that we also don't like that, and after doing the original,
somewhat incomplete, implementation using the subclassed approach,
I began to feel that implementing it in C (albeit using a different
approach for the code generation) was probably easier and more
elegant than what was done here. So you are very likely to see
it integrated as a regular numeric type, with a more C-based
implementation.

> 2)  Also,  in your C-API, you have a different pointer to the
> imaginary data.
>   I much prefer the way it is done currently to have complex numbers
> represented as an 8-byte, or 16-byte chunk of contiguous memory.
>
Any reason not to allow both? (The pointer to the real can be interpreted
as either a pointer to 8-byte or 16-byte quantities). It is true
that figuring out the imaginary pointer from the real is trivial
so I suppose it really isn't necessary.
>
> Index Arrays:
> ===========
>
> 1)  For what it's worth, my initial reaction to your indexing scheme is
> negative.  I would prefer that if
>
> a = [[1,2,3,4],
>       [5,6,7,8],
>       [9,10,11,12],
>       [13,14,15,16]]
>
> then
>
> a[[1,3],[0,3]] returns the sub-matrix:
>
> [[   4,  6],
>  [ 12, 14]
>
> i.e. the cross-product of [1,3] x [0,3]   This is the way MATLAB
> works.  I'm
> not sure what IDL does.
>
I'm afraid I don't understand the example. Could you elaborate
a bit more how this is supposed to work? (Or is it possible
there is an error? I would understand it if the result were
[[5, 8],[13,16]] corresponding to the index pairs
[[(1,0),(1,3)],[(3,0),(3,3)]])

> If I understand your scheme, right now, then I would have to
> append an extra
> dimension to my indexing arrays to get this behavior, right?
>
> 2) I would like to be able to index the array in a flattenned
> sense as well
> (is that possible?)  in other words, it would be nice if
> a[flat(9,10,11)] or
> something got me the elements 9,10,11 in a one-dimensional
> interpretation of
> the array.
>
Why not:

ravel(a)[[9,10,11]] ?

> 3) Why can't you combine slice notation and indexing?  Just interpret the
> slice as index array that would be created from using tha range
> operator on
> the same start, stop, and step objects.  Is this the plan?
>
I think that allowing slicing could be possible. But things were
getting pretty complex as they were, and we wanted to see if
there was agreement on how it was being done so far. It could
be extended to handle slices, if there was a well defined
interpretation. (I think there may be at least two possible
interpretations considered). As for the above, sure, but of
course the slice would have to be shape consistent with
the other index arrays (under the current scheme).

> That's all for now.  I don't mean to be critical, I'm really
> impressed with
> what works so far.  These are just some concerns I have right now.
>
> -Travis Oliphant
>
Thanks Travis, we're looking for constructive feedback, positive
or negative.

Perry

[Numpy-discussion] Numeric Python EM Project may need mirror

[Rob <eu...@ho...>]

Hi all,

I just got an email from @home yesterday, saying that all customers
should back up their web pages, email, etc etc.  I know they are in
bankruptcy, but this email sounded ominous.  I'm wondering if there is
some kindly soul who would want to mirror this site.   I'd really love
to have this site on Starship Python, but haven't had any responses to
emails to them.  

I'm continuously working on more code for the site so I'd hate to see it
go down, even if temporarily.

Sincerely,  Rob.
-- 
The Numeric Python EM Project

www.members.home.net/europax

Re: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Travis O. <oli...@ie...>]

>
> While we think this version is not yet mature enough for
> most to use in everyday projects, we are interested in
> feedback on the user interface and the open issues (see
> the documents on the web page shown below). We also welcome
> those who would like to contribute to this effort by helping
> with the development or adding libraries.
>

What I've seen looks great.   You've all done some good work here.  

Of course, I do have some feedback.  I haven't looked at everything, these 
points have just caught my eye.

Complex Types:
==============

1)  I don't like the idea of complex types being a separate subclass of 
ndarray.  This makes them "different."   Unless this "difference" can be 
completely hidden (which I doubt), I would prefer complex types to be on the 
same level as other numeric types.

2)  Also,  in your C-API, you have a different pointer to the imaginary data. 
  I much prefer the way it is done currently to have complex numbers 
represented as an 8-byte, or 16-byte chunk of contiguous memory. 


Index Arrays:
===========

1)  For what it's worth, my initial reaction to your indexing scheme is 
negative.  I would prefer that if

a = [[1,2,3,4],
      [5,6,7,8],
      [9,10,11,12],
      [13,14,15,16]]

then 

a[[1,3],[0,3]] returns the sub-matrix:

[[   4,  6],
 [ 12, 14]

i.e. the cross-product of [1,3] x [0,3]   This is the way MATLAB works.  I'm 
not sure what IDL does.

If I understand your scheme, right now, then I would have to append an extra 
dimension to my indexing arrays to get this behavior, right?

2) I would like to be able to index the array in a flattenned sense as well 
(is that possible?)  in other words, it would be nice if a[flat(9,10,11)] or 
something got me the elements 9,10,11 in a one-dimensional interpretation of 
the array.

3) Why can't you combine slice notation and indexing?  Just interpret the 
slice as index array that would be created from using tha range operator on 
the same start, stop, and step objects.  Is this the plan?

That's all for now.  I don't mean to be critical, I'm really impressed with 
what works so far.  These are just some concerns I have right now.

-Travis Oliphant

Re: [Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Pete S. <pe...@sh...>]

Perry Greenfield wrote:

> An early beta version is available on sourceforge as the
> package Numarray (http://sourceforge.net/projects/numpy/)
> 
> Information on the goals, changes in user interface, open issues,
> and design can be found at http://aten.stsci.edu/numarray

you ask a few questions on the information website, here are some of my 
answers for things i "care" about.

note that my main use of numpy is as a pixel buffer for images. some of 
the changes like avoiding type promotion sounds really good to me :]

5) should the implementation be bulletproof for private vars?
i don't think you should worry about this. as long as the interface is 
well defined, i wouldn't worry about protecting users from themselves. i 
this it will be the rare numarray user will be in a situation where they 
need to modify the internal C data.

7) necessary to add other types?
yes. i really want unsigned int16 and unsigned int32. all my operations 
are on pixel data, and things can just get messy when i need to treat 
packed color values as signed integers.

8) negative and out-of-range indices?
i'd prefer them to be kept as similar to python as can be. the current 
implementation in Numeric is nice for me.


one other thing i'd like there to be a little focus on is adding my own 
new ufunc operators. for image manipulation i'd like new ufunc operators 
that clamp the results to legal values. i'd be happy to do this myself, 
but i don't believe it's possible with the current Numeric.

the last thing i really really want is for this to be rolled into the 
standard python distribution. that is perhaps the most important aspect 
for me. i do not like requiring the extra dependency for generic numeric 
arrays. :]

[Numpy-discussion] Re-implementation of Python Numerical arrays (Numeric) available for download

[Perry G. <pe...@st...>]

We have been working on a reimplementation of Numeric, the
numeric array manipulation extension module for Python. 
The reimplementation is virtually a complete rewrite
and because it is not completely backwards compatible
with Numeric, we have dubbed it numarray to prevent
confusion.

While we think this version is not yet mature enough for
most to use in everyday projects, we are interested in
feedback on the user interface and the open issues (see 
the documents on the web page shown below). We also welcome
those who would like to contribute to this effort by helping
with the development or adding libraries.

An early beta version is available on sourceforge as the
package Numarray (http://sourceforge.net/projects/numpy/)

Information on the goals, changes in user interface, open issues,
and design can be found at http://aten.stsci.edu/numarray

Re: [Numpy-discussion] Numpy, BLAS, LAPACK, f2py

[Gerard V. <gve...@la...>]

Hi,

Try to link in the blas library (there is a dcopy_ in my blas library,
but better check the README first).

best regards -- Gerard

On Thursday 15 November 2001 11:15, Nils Wagner wrote:
> Hi,
>
> I have installed f2py on my system for wrapping existing FORTRAN 77
> codes to Python.
> Then I have gone through the following steps
>
> An example for using a TLS (total least squares routine)
> http://www.netlib.org/vanhuffel/
>
> 2) Get dtsl.f with dependencies
> 3) Run
>    f2py dtsl.f -m foo -h foo.pyf only: dtsl
>         \         \      \       \________ just wrap dtsl function
>          \         \      \______ create signature file
>           \         \____ python module name
>            \_____ Fortran 77 code
> 4) Edit foo.pyf to your specific needs (optional)
> 5) Run
>    f2py foo.pyf
>    \_____________ this will create Python C/API module foomodule.c
> 6) Run
>    make -f Makefile-foo
>    \_____________ this will build the module
> 7) In python:
>
> Python 2.1.1 (#1, Sep 24 2001, 05:28:47)
> [GCC 2.95.3 20010315 (SuSE)] on linux2
> Type "copyright", "credits" or "license" for more information.
>
> >>> import foo
>
> Traceback (most recent call last):
>   File "<stdin>", line 1, in ?
> ImportError: ./foomodule.so: undefined symbol: dcopy_
>
>
> Any suggestions to solve this problem ?
>
> Nils
>
> There are prebuilt libraries of LAPACK and BLAS in /usr/lib
>
> -rw-r--r--    1 root     root       657706 Sep 24 01:00 libblas.a
> lrwxrwxrwx    1 root     root           12 Okt 22 19:27 libblas.so ->
> libblas.so.2
> lrwxrwxrwx    1 root     root           16 Okt 22 19:27 libblas.so.2 ->
> libblas.so.2.2.0
> -rwxr-xr-x    1 root     root       559600 Sep 24 01:01 libblas.so.2.2.0
> -rw-r--r--    1 root     root      5763150 Sep 24 01:00 liblapack.a
> lrwxrwxrwx    1 root     root           14 Okt 22 19:27 liblapack.so ->
> liblapack.so.3
> lrwxrwxrwx    1 root     root           18 Okt 22 19:27 liblapack.so.3
> -> liblapack.so.3.0.0
> -rwxr-xr-x    1 root     root      4826626 Sep 24 01:01
> liblapack.so.3.0.0
>
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

[Numpy-discussion] Numpy, BLAS, LAPACK, f2py

[Nils W. <nw...@me...>]

Hi,

I have installed f2py on my system for wrapping existing FORTRAN 77
codes to Python.
Then I have gone through the following steps

An example for using a TLS (total least squares routine)
http://www.netlib.org/vanhuffel/

2) Get dtsl.f with dependencies
3) Run
   f2py dtsl.f -m foo -h foo.pyf only: dtsl
        \         \      \       \________ just wrap dtsl function      
         \         \      \______ create signature file
          \         \____ python module name
           \_____ Fortran 77 code
4) Edit foo.pyf to your specific needs (optional)
5) Run
   f2py foo.pyf
   \_____________ this will create Python C/API module foomodule.c
6) Run
   make -f Makefile-foo
   \_____________ this will build the module
7) In python:

Python 2.1.1 (#1, Sep 24 2001, 05:28:47)
[GCC 2.95.3 20010315 (SuSE)] on linux2
Type "copyright", "credits" or "license" for more information.
>>> import foo
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ImportError: ./foomodule.so: undefined symbol: dcopy_
>>>

Any suggestions to solve this problem ?

Nils

There are prebuilt libraries of LAPACK and BLAS in /usr/lib

-rw-r--r--    1 root     root       657706 Sep 24 01:00 libblas.a
lrwxrwxrwx    1 root     root           12 Okt 22 19:27 libblas.so ->
libblas.so.2
lrwxrwxrwx    1 root     root           16 Okt 22 19:27 libblas.so.2 ->
libblas.so.2.2.0
-rwxr-xr-x    1 root     root       559600 Sep 24 01:01 libblas.so.2.2.0
-rw-r--r--    1 root     root      5763150 Sep 24 01:00 liblapack.a
lrwxrwxrwx    1 root     root           14 Okt 22 19:27 liblapack.so ->
liblapack.so.3
lrwxrwxrwx    1 root     root           18 Okt 22 19:27 liblapack.so.3
-> liblapack.so.3.0.0
-rwxr-xr-x    1 root     root      4826626 Sep 24 01:01
liblapack.so.3.0.0

Re: [Numpy-discussion] Total least squares problem

[Nils W. <nw...@me...>]

Konrad Hinsen schrieb:
> 
> Nils Wagner <nw...@me...> writes:
> 
> > There is a difference between classical least squares (Numpy)
> > and TLS (total least squares).
> 
> Algorithmically speaking it is even a very different problem. I'd say
> the only reasonable (i.e. efficient) solution for NumPy is to
> implement the TLS algorithm in a C subroutine calling LAPACK routines
> for SVD etc.
> 
> Konrad.
> --
There are two Fortran implementations of the TLS algorithm already
available via
http://www.netlib.org/vanhuffel/      
.
Moreover there is a tool called f2py that generates Python C/API modules
for wrapping Fortran 77/90/95 codes to Python.

Unfortunately I am not very familar with this tool.
 
Therefore I need some advice for this.

Thanks in advance

 Nils
 
  
> -------------------------------------------------------------------------------
> Konrad Hinsen                            | E-Mail: hi...@cn...
> Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
> Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
> 45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
> France                                   | Nederlands/Francais
> -------------------------------------------------------------------------------

Re: [Numpy-discussion] Total least squares problem

[Konrad H. <hi...@cn...>]

Nils Wagner <nw...@me...> writes:

> There is a difference between classical least squares (Numpy)
> and TLS (total least squares).

Algorithmically speaking it is even a very different problem. I'd say
the only reasonable (i.e. efficient) solution for NumPy is to
implement the TLS algorithm in a C subroutine calling LAPACK routines
for SVD etc.

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hi...@cn...
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

Re: [Numpy-discussion] Total least squares problem

[Nils W. <nw...@me...>]

Travis Oliphant schrieb:
> 
> >
> > How do I solve a Total Least Squares problem in Numpy ?
> > A small example would be appreciated.
> >
> > The TLS problem assumes an overdetermined set of linear equations
> > AX = B, where both the data matrix A as well as the observation
> > matrix B are inaccurate:
> 
> X, resids, rank, s = LinearAlgebra.linear_least_squares(A,B)
> 
> -Travis

Travis,

There is a difference between classical least squares (Numpy)
and TLS (total least squares).
I am attaching a small example for illustration.

Nils

[Numpy-discussion] BLAS and innerproduct

[<R.M...@ex...>]

Hello,

So far as I can tell Numeric.dot(), which uses innerproduct() from
multiarraymodule.c doesn't call the BLAS, even if Numeric was compiled
against native BLAS.   This means (at least on my machine) that 

X = ones((150, 16384), 'd')
C = dot(X, tranpose(X))

is about 15 times as slow as the comparable operations in Matlab (v6),
which does, I think, use the native BLAS.

I guess that multiarray.c is not particularly optimised to use the
BLAS because of the difficulties of coping with all sorts of types
(float32, int64 etc), and with non-contiguous arrays.  The
innerproduct is so basic to most of the work I use Numeric for that a
speed up here would make a big difference.  I'm thinking of patching
multiarray.c to use the BLAS when it can, but before I start are there
good reasons for doing something different?

Any advice gratefully received!

Cheers,

Richard.


-- 
Department of Computer Science, Exeter University    Voice: +44 1392 264065
R.M...@ex...                         Secretary: +44 1392 264061
http://www.dcs.ex.ac.uk/people/reverson                Fax: +44 1392 264067

Re: [Numpy-discussion] Total least squares problem

[Travis O. <oli...@ee...>]

>
> How do I solve a Total Least Squares problem in Numpy ?
> A small example would be appreciated.
>
> The TLS problem assumes an overdetermined set of linear equations
> AX = B, where both the data matrix A as well as the observation
> matrix B are inaccurate:


X, resids, rank, s = LinearAlgebra.linear_least_squares(A,B)

-Travis