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

Re: [Numpy-discussion] optimizing power() for complex and real cases

[<co...@ph...>]

Tim Hochberg <tim...@co...> writes:

> Here's a little progress report: I now have A**2 running as fast as
> square(A). This is done by special casing stuff in array_power so that
> A**2 acutally calls square(A) instead of going through power(A,2).
> Things still need a bunch of cleaning up (in fact right now A**1
> segfaults, but I know why and it should be an easy fix). However, I
> think I've discovered why you couldn't get your special cased power to
> run as fast for A**2 as square(A) or A*A. It appears that the overhead
> of creating a new array object from the integer 2 is the bottleneck. I
> was running into the same mysterious overhead, even when dispatching
> from array_power, until I special cased on PyInt to avoid the array
> creation in that case.

Hmm, if that's true about the overhead, that'll hit all computations
of the type op(x, <some scalar>). Something to look at. That ufunc code
for casting the arguments is pretty big and hairy, so I'm not going to
look at right now :-)

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

Re: [Numpy-discussion] cross product of two 3xn arrays

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

Charles R Harris wrote:

>On 2/17/06, Bill Baxter <wb...@gm...> wrote:
>  
>
>>For folks using quats to represent rotations (which is all I use them for,
>>anyway), if you're batch transforming a bunch of vectors by one quaternion,
>>it's a lot more efficient to convert the quat to a 3x3 matrix first and
>>transform using matrix multiply (9 mults per transform that way vs 21 or so
>>depending on the implementation of q*v*q^-1).  Given that, I can't see many
>>situations when I'd need a super speedy C version of quaternion multiply.
>>
>>--Bill
>>    
>>
>
>On a related note, indexing in numarray is some 3x faster than in
>NumPy and I'm wondering what needs to be done to speed that up.
>  
>
Please explain with a benchmark.  This is not true for all indexing 
operations.  But, it is possible that certain use cases are faster.    
We can't do anything without knowing what you are talking about exactly.


-Travis

Re: [Numpy-discussion] Dot products and casting

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

Ed Schofield wrote:

>Hi all,
>
>  
>
>The problem seems to be with the PyArray_ObjectType() calls in
>dotblas_matrixproduct(), which are returning typenum=PyArray_FLOAT, but
>this isn't sufficiently large for a safe cast from the int32 and int64
>arrays.  It seems like PyArray_ObjectType() should be returning
>PyArray_DOUBLE here instead.
>  
>
This sounds like an accurate diagnosis.  I'll have to look at the 
type-evaluation code a bit more to see why a suitable type is not being 
found --- unless someone else gets there first.  I won't have time for 
awhile today.

-Travis

Re: [Numpy-discussion] optimizing power() for complex and real cases

[Tim H. <tim...@co...>]

Here's a little progress report: I now have A**2 running as fast as 
square(A). This is done by special casing stuff in array_power so that 
A**2 acutally calls square(A) instead of going through power(A,2). 
Things still need a bunch of cleaning up (in fact right now A**1 
segfaults, but I know why and it should be an easy fix). However, I 
think I've discovered why you couldn't get your special cased power to 
run as fast for A**2 as square(A) or A*A. It appears that the overhead 
of creating a new array object from the integer 2 is the bottleneck. I 
was running into the same mysterious overhead, even when dispatching 
from array_power, until I special cased on PyInt to avoid the array 
creation in that case.

-tim

[Numpy-discussion] PyTables with support for NumPy 0.9.5

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

Hi,

I've just uploaded a new version of PyTables with (almost) complete
support for the recent NumPy 0.9.5. All the range of homogeneous and
heterogeneous (including those with nested fields) arrays using any
combination of data-types should be supported. The only exception is
the lack of support of unicode types (I have to figure out yet which
HDF5 datatype would be best to mapping them; suggestions are
welcome!).

You can fetch the tarball from:

http://pytables.carabos.com/download/preliminary/pytables-1.3beta2.tar.gz

Test it as much as you can, and if you find any strange quirk, do not
hesitate to report it back.

Regards,

=2D-=20
>0,0<   Francesc Altet =A0 =A0 http://www.carabos.com/
V   V   C=E1rabos Coop. V. =A0=A0Enjoy Data
 "-"

Re: [Numpy-discussion] ANN: Python Enthought Edition Version 0.9.2 Released

[Travis N. V. <tr...@en...>]

Arnd Baecker wrote:
> On Thu, 16 Feb 2006, Travis N. Vaught wrote:
>
>   
>> Enthought is pleased to announce the release of Python Enthought Edition
>> Version 0.9.2 (http://code.enthought.com/enthon/) -- a python
>> distribution for Windows.  This is a kitchen-sink-included Python
>> distribution including the following packages/tools out of the box:
>>
>> Numeric 24.2
>> SciPy 0.3.3
>> IPython 0.6.15
>> Enthought Tool Suite 1.0.2
>> wxPython 2.6.1.0
>> PIL 1.1.4
>> mingw 20030504-1
>> f2py 2.45.241_1926
>> MayaVi 1.5
>> Scientific Python 2.4.5
>> VTK 4.4
>> and many more...
>>     
>
> Brilliant - many thanks for the effort!
>
> I was just about to ask for the plans about numpy/scipy,
> but the changelog at
> http://code.enthought.com/release/changelog-enthon0.9.2.shtml
> shows quite a bit of activity in this direction!
>
> Do you have an estimate about when a numpy/scipy version
> of the Enthought Edition might happen?
>
> Many thanks,
>
> Arnd
>   
It's a bit difficult to say with much accuracy, so I'll be transparent 
but imprecise.  Our release of Enthon versions typically tracks the 
state of the platform we are using for the custom software development 
we do to pay the bills.  Thus, our current project code typically has to 
be ported to build and run on a cobbled-together build of the newer 
versions before we do a release.  I realize this is a drag on the 
release schedule for Enthon, but it's how we allocate resources to the 
builds.

Enough excuses, though--we are working on the migration of our project 
code now (Pearu Peterson) and I expect in weeks (rather than months) 
we'll have an Enthon release candidate with Python 2.4.2, and the latest 
SciPy and NumPy on Windows.

Robert Kern is already working on a project that is based on this tool 
chain, so the wedge is in place.

Thanks for the interest! (and sorry for the cross-post)

Travis

Re: [Numpy-discussion] cross product of two 3xn arrays

[Charles R H. <cha...@gm...>]

On 2/17/06, Bill Baxter <wb...@gm...> wrote:
> For folks using quats to represent rotations (which is all I use them for=
,
> anyway), if you're batch transforming a bunch of vectors by one quaternio=
n,
> it's a lot more efficient to convert the quat to a 3x3 matrix first and
> transform using matrix multiply (9 mults per transform that way vs 21 or =
so
> depending on the implementation of q*v*q^-1).  Given that, I can't see ma=
ny
> situations when I'd need a super speedy C version of quaternion multiply.
>
> --Bill

True. On the other hand, I have files containing 20,000 quaternions,
each of which needs to be converted to a rotation matrix and applied
to some 400 vectors, so a c version would have its place. I have a
python quaternion class that I use for such things but profiling shows
that it is one of the prime time bandits, so I am tempted to use c for
that class anyway. Note that the current NumPy cross product is
implemented in Python.

On a related note, indexing in numarray is some 3x faster than in
NumPy and I'm wondering what needs to be done to speed that up.

Chuck

Re: [Numpy-discussion] ANN: Python Enthought Edition Version 0.9.2 Released

[Arnd B. <arn...@we...>]

On Thu, 16 Feb 2006, Travis N. Vaught wrote:

> Enthought is pleased to announce the release of Python Enthought Edition
> Version 0.9.2 (http://code.enthought.com/enthon/) -- a python
> distribution for Windows.  This is a kitchen-sink-included Python
> distribution including the following packages/tools out of the box:
>
> Numeric 24.2
> SciPy 0.3.3
> IPython 0.6.15
> Enthought Tool Suite 1.0.2
> wxPython 2.6.1.0
> PIL 1.1.4
> mingw 20030504-1
> f2py 2.45.241_1926
> MayaVi 1.5
> Scientific Python 2.4.5
> VTK 4.4
> and many more...

Brilliant - many thanks for the effort!

I was just about to ask for the plans about numpy/scipy,
but the changelog at
http://code.enthought.com/release/changelog-enthon0.9.2.shtml
shows quite a bit of activity in this direction!

Do you have an estimate about when a numpy/scipy version
of the Enthought Edition might happen?

Many thanks,

Arnd

Re: [Numpy-discussion] storage for records

[Stefan v. d. W. <st...@su...>]

I am probably trying to do something silly, but still:

In [1]: import numpy as N

In [2]: N.__version__
Out[2]: '0.9.6.2127'

In [3]: P = N.array(N.zeros((2,2)), N.dtype((('f4',3), {'names': ['x','y','z'], 'formats': ['f4','f4','f4']})))
*** glibc detected *** malloc(): memory corruption: 0x0830bb48 ***
Aborted

Regards
Stéfan

On Thu, Feb 16, 2006 at 02:06:15PM -0700, Travis Oliphant wrote:
> Stefan van der Walt wrote:
> 
> >Is there any way to control the underlying storage for a record?
> >
> >I am trying to use Travis' earlier example of an image with named fields:
> >
> >dt = N.dtype('<f12', [('r','<f4'),('g','<f4'),('b','<f4')])
> >img = N.array(N.empty((rows,columns)), dtype=dt)

[Numpy-discussion] Dot products and casting

[<sch...@ft...>]

Hi all,

I think there's a bug in dot() that prevents it from operating on two
arrays, neither of which can be safely cast to the other.  Here's an
example:

>>> from numpy import *
>>> a = arange(10, dtype=float32)
>>> b = arange(10, dtype=float64)
>>> c = arange(10, dtype=int64)
>>> d = arange(10, dtype=int32)
>>> e = arange(10, dtype=int16)

# Dot products between b and either c or d work fine:
>>> dot(b,c)
285.0
>>> dot(b,d)
285.0

# Dot products with e also work fine:
>>> dot(a,e)
285.0
>>> dot(b,e)
285.0

# But dot products between a and either c or d don't work:
>>> dot(a,c)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: array cannot be safely cast to required type
>>> dot(a,d)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: array cannot be safely cast to required type


The problem seems to be with the PyArray_ObjectType() calls in
dotblas_matrixproduct(), which are returning typenum=PyArray_FLOAT, but
this isn't sufficiently large for a safe cast from the int32 and int64
arrays.  It seems like PyArray_ObjectType() should be returning
PyArray_DOUBLE here instead.

Here's another example:

>>> f = arange(10, dtype=complex64)
>>> dot(b, f)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: array cannot be safely cast to required type

So it seems like the problem isn't isolated to float32 arrays, but occurs
elsewhere when we need to find a minimum data type of two arrays when
*both* need to be upcasted.


-- Ed

[Numpy-discussion] maxima of masked arrays

[George N. <ag...@no...>]

I am trying to get the n-1 dimensional array of maxima of an array  
taken over a given axis.

with ordinary arrays this works fine.

E.g.
In [49]: a = arange(1,13).reshape(3,4)

In [50]: a
Out[50]:
array([[ 1,  2,  3,  4],
        [ 5,  6,  7,  8],
        [ 9, 10, 11, 12]])

The maximum over all elements is:
In [51]: print a.max()
12

& the array of maxima over the 0-axis is OK too:
In [52]: print a.max(0)
[ 9 10 11 12]

But with a masked array there are problems.

In [54]: amask = ma.masked_where(a < 5,a)

In [55]: amask
Out[55]:
array(data =
[[999999 999999 999999 999999]
[     5      6      7      8]
[     9     10     11     12]],
       mask =
[[True True True True]
[False False False False]
[False False False False]],
       fill_value=999999)


The maximum over all elements is fine:
In [56]: amask.max()
Out[56]: 12

but trying to get an array of maxima over the 0-axis fails:

n [57]: amask.max(0)
Out[57]:
array(data =
[[999999 999999 999999 999999]
[     5      6      7      8]
[     9     10     11     12]],
       mask =
[[True True True True]
.......


Are there any workarounds for this?

-George.

[Numpy-discussion] Re: calculating on matrix indices

[Brian B. <bb...@br...>]

Robert Kern wrote:
> The traceback tells you exactly what's wrong:
> 
> In [7]: x[idx] = exp(-t[idx]/tau)
> ---------------------------------------------------------------------------
> exceptions.TypeError                                 Traceback (most recent call
> last)
> 

yes, I saw that, but all of the types (i.e. type(x)) came out to be the same, so I 
figured the problem was with the indexing, and that was causing a typecast problem. 
I didn't know about dtype

> In [13]: x = zeros(len(t), float)

well that is confusing!  zeros(5,'f') is single precision, zeros(5,'d') is double, 
and zeros(5,float) is double!  that's where I got whacked, because I remembered that 
"float" was "double" in python...but I guess, not always.


thanks for your help!



				bb
-- 
-----------------

             bb...@br...
             http://web.bryant.edu/~bblais

Re: [Numpy-discussion] Problems compiling on Cygwin

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

Jose Gomez-Dans wrote:

>Hi!
>Yesterday I posted on the scipy mailing list that I could not compile NumPy on
>Cygwin. I would like to provide some more information on what the problems are,
>as I would really like to be able to use it on Cygwin.
>  
>
I looked in to how people at cygwin ports got the IEEE math stuff done.  
They borrowed it from BSD basically.  So, I've taken their patch and 
placed it in the main tree.

Jose, could you check out the latest SVN version of numpy and try to 
build and install it on cygwin to see if I made the right changes?

-Travis

Re: [Numpy-discussion] Problems compiling on Cygwin

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

Jose Gomez-Dans wrote:

>Hi!
>Yesterday I posted on the scipy mailing list that I could not compile NumPy on
>Cygwin. I would like to provide some more information on what the problems are,
>as I would really like to be able to use it on Cygwin.
>  
>
Thanks Jose.

It looks like we are not doing the right thing in the platform-specific 
section of code here.

But the right thing can potentially be done. 

Look here:

ftp://sunsite.dk/projects/cygwinports/release/python/numpy/

It looks like somebody figured out how to make it work with cygwin (one 
option, of course is to just disable the IEEE error-setting modes for 
cygwin).


-Travis

[Numpy-discussion] Re: ANN: Release of NumPy 0.9.5

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

Thomas Gellekum wrote:
> "Travis E. Oliphant" <oli...@ie...> writes:
> 
> 
>>  -  Improvements to numpy.distutils
> 
> 
> Stupid questions: is it really necessary to keep your own copy of
> distutils and even install it? What's wrong with the version in the
> Python distribution? How can I disable the colorized output to get
> something more readable than yellow on white (well, seashell1)?
> 

Yes --- distutils does not provide enough functionality.

Besides, it's not a *copy* of distutils.  It's enhancements to 
distutils.  It builds on top of standard distutils.

I don't know the answer to the colorized output question.  Please post 
to num...@li... to get more help and/or make 
suggestions.

-Travis

[Numpy-discussion] Problems compiling on Cygwin

[Jose Gomez-D. <jos...@gm...>]

Hi!
Yesterday I posted on the scipy mailing list that I could not compile NumPy on
Cygwin. I would like to provide some more information on what the problems are,
as I would really like to be able to use it on Cygwin.

I got the 0.9.5 tarball, and uncompress it, and type python setup.py build. The
process starts, there is an indication that it finds BLAS and LAPACK (cygwin
versions). It stops when linking the umath.dll, complaining about missing
references. Here's an extract:
"gcc options: '-fno-strict-aliasing -DNDEBUG -g -O3 -Wall -Wstrict-prototypes'
compile options: '-Ibuild/src/numpy/core/src -Inumpy/core/include
-Ibuild/src/numpy/core -Inumpy/core/src -Inumpy/core/include
-I/usr/include/python2.4 -c'
gcc -shared -Wl,--enable-auto-image-base
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o
-L/usr/lib/python2.4/config -lpython2.4 -o
build/lib.cygwin-1.5.19-i686-2.4/numpy/core/umath.dll
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x2f5e):
referencia a `_feraiseexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x2fe3):
referencia a `_feraiseexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x3081):
referencia a `_feraiseexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x3139):
referencia a `_feraiseexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x320c):
referencia a `_feraiseexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x129ef):
referencia a `_fetestexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x12a1d):
referencia a `_feclearexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x12b1b):
referencia a `_fetestexcept' sin definir
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o:
umathmodule.c:(.text+0x12b27):
referencia a `_feclearexcept' sin definir
collect2: ld devolvi'o el estado de salida 1
error: Command "gcc -shared -Wl,--enable-auto-image-base
build/temp.cygwin-1.5.19-i686-2.4/build/src/numpy/core/src/umathmodule.o
-L/usr/lib/python2.4/config -lpython2.4 -o
build/lib.cygwin-1.5.19-i686-2.4/numpy/core/umath.dll" failed with exit 
status 1"

(yes, I have the Spanish locale set :D). The functions needed are responsible
for setting exceptions, and presumably, only need a simple addition to the
library linking path. Is this correct? does anyone know how to deal with this?

Many thanks!
Jose

[Numpy-discussion] Re: [matplotlib-devel] Numpy 0.9.5 + pylab 0.86.2 = python.exe crash

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

Bill Baxter wrote:

> After upgrading to numpy 0.9.5 (numpy-0.9.5.win32-py2.4.exe) running 
> the following script causes a crash of python.exe:
>
> -------------test.py-----------------
> import pylab
> ----------------------------------------
>
>
> In my .matplotlibrc file I have the following:
> -------------
> backend : WXAgg
> numerix : numpy
> -------------
>
>
> Reinstalling numpy 0.9.4 fixed the problem.
>
> Matplotlib version is 0.86.2-win32-py2.4
> I also tried reinstalling matplotlib, but that didn't help.
>
You have to re-compile the matplotlib extension.   There are warnings 
present now so that hopefully in the future such needs will be 
comunicated better.

-Travis

Re: [Numpy-discussion] Numpy 0.9.5 + pylab 0.86.2 = python.exe crash

[Bill B. <wb...@gm...>]

Ew.  Ok.  No thanks.  :-)  I'll just stick with numpy 0.9.4 for now.

I appreciate the speedy response.

--bb

On 2/17/06, David M. Cooke <co...@ph...> wrote:
>
> Bill Baxter <wb...@gm...> writes:
>
> > After upgrading to numpy 0.9.5 (numpy-0.9.5.win32-py2.4.exe) running th=
e
> > following script causes a crash of python.exe:
> >
> > -------------test.py-----------------
> > import pylab
> > ----------------------------------------
> >
> >
> > In my .matplotlibrc file I have the following:
> > -------------
> > backend : WXAgg
> > numerix : numpy
> > -------------
> >
> >
> > Reinstalling numpy 0.9.4 fixed the problem.
> >
> > Matplotlib version is 0.86.2-win32-py2.4
> > I also tried reinstalling matplotlib, but that didn't help.
>
> You'll have to recompile matplotlib against the newer numpy.
>
> --
> |>|\/|<
>
> /------------------------------------------------------------------------=
--\
> |David M. Cooke
> http://arbutus.physics.mcmaster.ca/dmc/
> |co...@ph...
>



--
William V. Baxter III
OLM Digital
Kono Dens Building Rm 302
1-8-8 Wakabayashi Setagaya-ku
Tokyo, Japan  154-0023
+81 (3) 3422-3380

Re: [Numpy-discussion] Numpy 0.9.5 + pylab 0.86.2 = python.exe crash

[<co...@ph...>]

Bill Baxter <wb...@gm...> writes:

> After upgrading to numpy 0.9.5 (numpy-0.9.5.win32-py2.4.exe) running the
> following script causes a crash of python.exe:
>
> -------------test.py-----------------
> import pylab
> ----------------------------------------
>
>
> In my .matplotlibrc file I have the following:
> -------------
> backend : WXAgg
> numerix : numpy
> -------------
>
>
> Reinstalling numpy 0.9.4 fixed the problem.
>
> Matplotlib version is 0.86.2-win32-py2.4
> I also tried reinstalling matplotlib, but that didn't help.

You'll have to recompile matplotlib against the newer numpy.

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

[Numpy-discussion] Numpy 0.9.5 + pylab 0.86.2 = python.exe crash

[Bill B. <wb...@gm...>]

After upgrading to numpy 0.9.5 (numpy-0.9.5.win32-py2.4.exe) running the
following script causes a crash of python.exe:

-------------test.py-----------------
import pylab
----------------------------------------


In my .matplotlibrc file I have the following:
-------------
backend : WXAgg
numerix : numpy
-------------


Reinstalling numpy 0.9.4 fixed the problem.

Matplotlib version is 0.86.2-win32-py2.4
I also tried reinstalling matplotlib, but that didn't help.

--Bill Baxter

Re: [Numpy-discussion] cross product of two 3xn arrays

[Bill B. <wb...@gm...>]

For folks using quats to represent rotations (which is all I use them for,
anyway), if you're batch transforming a bunch of vectors by one quaternion,
it's a lot more efficient to convert the quat to a 3x3 matrix first and
transform using matrix multiply (9 mults per transform that way vs 21 or so
depending on the implementation of q*v*q^-1).  Given that, I can't see many
situations when I'd need a super speedy C version of quaternion multiply.

--Bill

Re: [Numpy-discussion] optimizing power() for complex and real cases

[<co...@ph...>]

Tim Hochberg <tim...@co...> writes:
> Gary Ruben wrote:
> So, optimize:
>
> A**-1, A**0, A**1, A**2, etc.
>
>>
>> and
>>
>> A**array([1,2,3])
>>
>> but not
>>
>> A**array[1.0, 2.0, 3.0], A**2.0, A**(2.0+0j)
>
>
> That makes sense. It's safer and easier to explain: "numpy optimizes
> raising matrices (and possibly scalars) to integer powers)". The only
> sticking point that I see is if David is still interested in
> optimizing A**0.5, that's not going to mesh with this. On the other
> hand, perhaps he can be persuaded that sqrt(A) is just as good. After
> all, it's only one more character long ;)

sigh, ok :-)

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

Re: [Numpy-discussion] optimizing power() for complex and real cases

[<co...@ph...>]

Tim Hochberg <tim...@co...> writes:
> David M. Cooke wrote:
>
>>[1] which brings up another point. Would using the 3-multiplication
>>version for complex multiplication be good? There might be some
>>effects with cancellation errors due to the extra subtractions...
>>
>>
> I'm inclined to leave this be for now. Both because I'm unsure of the
> rounding issues and because I'm not sure it would actually be faster.
> It has one less multiplication, but several more additions, so it
> would depend on the relative speed add/sub with multiplication and how
> things end up getting scheduled in the FP pipeline. At some point it's
> probably worth trying; if it turns out to be signifigantly faster we
> can think about rounding then. If it's not faster then no need to
> think.

I did some thinking, and looked up how to analyse it. 3M goes like
this:

xy = (a+bi)(c+di) = (ac-bd) + ((a+c)(b+d)-ac-bd)i

Consider x = y = t + i/t, for which x**2 = (t**2-1/t**2) + 2i, then

xy=x^2 = t*t - 1/t*1/t + ((t+1/t)(t+1/t) - t**2 - 1/t**2)i

Consider when t is large enough that (t+1/t)**2 = t**2 in floating
point; then Im fl(xy) will be -1/t**2, instead of 2.

So...let's leave it as is.

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

Re: [Numpy-discussion] Creating arrays with fromfile

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

Charles R Harris wrote:

>Hi Travis,
>
>I notice that the fromfile function in NumPy no longer accepts the
>shape keyword that the numarray version has. The functionalitiy can be
>duplicated by reshaping the array after creating it, but I think the
>shape keyword is a bit more convenient for that. Thoughts?
>
>  
>
It was just that much more effort to implement correctly in C and since 
it can be easily done using

fromfile(....).reshape(dim1,dim2,dim3,...)

I didn't think it critical. 

Perhaps numarray compatibility functions should be placed in a numcompat 
module.

-Travis

Re: [Numpy-discussion] cross product of two 3xn arrays

[Charles R H. <cha...@gm...>]

Bill,

On 2/16/06, Bill Baxter <wb...@gm...> wrote:
> Quaternions using which convention?  [s,x,y,z] or [x,y,z,w]?
> The docstring should make it very clear.  Perhaps support a flag for
> choosing which, unless there's some python-wide standard for quaternions
> that I'm not aware of.
>
> --Bill
>
>
> On 2/17/06, Charles R Harris <cha...@gm...> wrote:
> > Would anyone be interested in a quaternion version of this for nx4
> > arrays with nx3 as a special case where the scalar part =3D=3D 0? Looki=
ng
> > at the the cross product implementation, it shouldn't be to hard to
> > duplicate this for quaternions. What should such a product be called?
> > Something like qprod?
> >
> > Chuck
> >
> >

I like to put the scalar last, but I am open to putting it first if
anyone has strong feelings about it. As far as I know, there is no
scipy convention on this. Hmm, maybe a flag would be useful just
because folks are likely to have files sitting around full of
quaternions using both conventions. Maybe one more scalar type to add
to the NumPy mix? I must admit that dtype=3DQuaternion512 seems a bit
much. Anyway, I am open to suggestions.