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

[Numpy-discussion] Release tomorrow

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

I just made some changes to the system_info file to parse site.cfg from 
the three standard locations (in order)

1) System-wide (the location of numpy/distutils/system_info.py)
2) Users HOME directory
3) Current working directory

I'm assuming the config parser will update the appropriate sections as 
later files are read.

I'd like to make a release tomorrow but would like a code-review on my 
changes in this section before hand.  If somebody who uses site.cfg 
could try out the SVN version and see if it works as expected.  That 
would be great.

-Travis

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

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

Ian Harrison wrote:

>Hello,
>
>I have two groups of 3x1 arrays that are arranged into two larger 3xn
>arrays. Each of the 3x1 sub-arrays represents a vector in 3D space. In
>Matlab, I'd use the function cross() to calculate the cross product of
>the corresponding 'vectors' from each array. In other words:
>  
>


Help on function cross in module numpy.core.numeric:

cross(a, b, axisa=-1, axisb=-1, axisc=-1)
    Return the cross product of two (arrays of) vectors.

    The cross product is performed over the last axis of a and b by default,
    and can handle axes with dimensions 2 and 3. For a dimension of 2,
    the z-component of the equivalent three-dimensional cross product is
    returned.

It's the axisa, axisb, and axisc  that you are interested in.

The default is to assume you have Nx3 arrays and return an Nx3 array.  
But you can change the axis used to find vectors.

cross(A,B,axisa=0,axisb=0,axisc=0)

will do what you want.  I suppose, a single axis= argument might be 
useful as well for the common situation of having all the other axis 
arguments be the same.

-Travis










>if ai and bj are 3x1 column vectors:
>
>A = [ a1  a2  a3 ]
>
>B = [ b1  b2  b3 ]
>
>C = A x B = [ (a1 x b1)  (a2 x b2)  (a3 x b3) ]
>
>Could someone suggest a clean way to do this? I suppose I could write
>a for loop to cycle through each pair of vectors and send them to
>numpy's cross(), but since I'm new to python/scipy/numpy, I'm guessing
>that there's probably a better method that I'm overlooking.
>
>Thanks,
>Ian
>
>
>-------------------------------------------------------
>This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
>for problems?  Stop!  Download the new AJAX search engine that makes
>searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
>http://sel.as-us.falkag.net/sel?cmd=k&kid3432&bid#0486&dat1642
>_______________________________________________
>Numpy-discussion mailing list
>Num...@li...
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>  
>

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

[Gary R. <gr...@bi...>]

This *almost* does what you want, I think. I can't see a neat way to 
give column vectors in the solution:

In [21]: a=array([[[1],[2],[3]],[[4],[5],[6]],[[7],[8],[9]]])

In [22]: b=array([[[1],[2],[4]],[[4],[5],[7]],[[7],[8],[10]]])

In [24]: cross(a.transpose(),b.transpose())
Out[24]:
array([[[ 2, -1,  0],
         [ 5, -4,  0],
         [ 8, -7,  0]]])

Gary R.

Ian Harrison wrote:
> Hello,
> 
> I have two groups of 3x1 arrays that are arranged into two larger 3xn
> arrays. Each of the 3x1 sub-arrays represents a vector in 3D space. In
> Matlab, I'd use the function cross() to calculate the cross product of
> the corresponding 'vectors' from each array. In other words:
> 
> if ai and bj are 3x1 column vectors:
> 
> A = [ a1  a2  a3 ]
> 
> B = [ b1  b2  b3 ]
> 
> C = A x B = [ (a1 x b1)  (a2 x b2)  (a3 x b3) ]
> 
> Could someone suggest a clean way to do this? I suppose I could write
> a for loop to cycle through each pair of vectors and send them to
> numpy's cross(), but since I'm new to python/scipy/numpy, I'm guessing
> that there's probably a better method that I'm overlooking.
> 
> Thanks,
> Ian

[Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

[<sk...@po...>]

    Travis> Yes, that is a bad place.  We need some suggestions as to w=
here
    Travis> site.cfg should be read from.

First place to look should be `pwd`.

    Travis> I think you can set the environment variable ATLAS to 'None=
' and
    Travis> it will ignore ATLAS...

Thank you, thank you, thank you.  I now have numpy built...  I'll tackl=
e the
rest of scipy ma=F1ana.

Skip

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

[Gary R. <gr...@bi...>]

Tim Hochberg wrote:

>> optimize:
>>
>> A**2, A**(2+0j), etc
> 
> 
> I'm guessing here that you did not mean to include (2+0j) on both lists 
> and that, in fact you wanted not to optimize on complex exponents:

Oops. The complex index in the optimise list has integer parts. I 
assumed Python distinguished between complex numbers with integer and 
real parts, but it doesn't, so you're correct about which cases I'd vote 
to optimise.

> 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 ;)
> 
> -tim

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

[Ian H. <har...@gm...>]

Hello,

I have two groups of 3x1 arrays that are arranged into two larger 3xn
arrays. Each of the 3x1 sub-arrays represents a vector in 3D space. In
Matlab, I'd use the function cross() to calculate the cross product of
the corresponding 'vectors' from each array. In other words:

if ai and bj are 3x1 column vectors:

A =3D [ a1  a2  a3 ]

B =3D [ b1  b2  b3 ]

C =3D A x B =3D [ (a1 x b1)  (a2 x b2)  (a3 x b3) ]

Could someone suggest a clean way to do this? I suppose I could write
a for loop to cycle through each pair of vectors and send them to
numpy's cross(), but since I'm new to python/scipy/numpy, I'm guessing
that there's probably a better method that I'm overlooking.

Thanks,
Ian

[Numpy-discussion] Re: New release of NumPy coming

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

Francesc Altet wrote:
> A Dimecres 15 Febrer 2006 11:18, Pearu Peterson va escriure:
> 
>>On Wed, 15 Feb 2006, Francesc Altet wrote:
>>
>>>I've been trying to see how to correctly load the unicode tests, but
>>>failed miserably. Perhaps Pearu can tell us about the correct way to
>>>do that.
>>
>>I have fixed it in svn. When importing modules from tests/ directory, one
>>must surround the corresponding import statements with set_local_path()
>>and restore_path() calls.
> 
> 
> Ah, ok. Is there any place where this is explained or we have to use
> the source to figure out these sort of things?
>

Good thing to put in the doc subdirectory....  I was obviously not 
certain about the use of set_local_path() myself until Pearu corrected me.

-Travis

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

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

Gary Ruben wrote:

> Tim Hochberg wrote:
> <big snip>
>
>> As I've been thinking about this some more, I think the correct thing 
>> to do is not to mess with  the power ufuncs at all. Rather in 
>> x.__pow__ (since I don't know that there's anywhere else to do it), 
>> after the above checks check the types of the array and in the cases 
>> where the first argument is a float or complex and the second 
>> argument is some sort of integer array. This would be dispatched to 
>> some other helper function instead of the normal pow_ufunc. In other 
>> words, optimize:
>>
>> A**2, A**2.0, A**(2.0+0j), etc
>>
>> and
>>
>> A**array([1,2,3])
>>
>> but not
>>
>> A**array[1.0, 2.0, 3.0]
>>
>> I think that this takes care of the optimization slowing down power 
>> for general floats and optimizes the only array-array case that 
>> really matter.
>
>
> I think this might still be a tiny bit dangerous despite not observing 
> monotonicity problems and would be a bit more conservative and change 
> it to:
>
> optimize:
>
> A**2, A**(2+0j), etc


I'm guessing here that you did not mean to include (2+0j) on both lists 
and that, in fact you wanted not to optimize on complex exponents:
.
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 ;)

-tim

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

[Gary R. <gr...@bi...>]

Tim Hochberg wrote:
<big snip>
> As I've been thinking about this some more, I think the correct thing to 
> do is not to mess with  the power ufuncs at all. Rather in x.__pow__ 
> (since I don't know that there's anywhere else to do it), after the 
> above checks check the types of the array and in the cases where the 
> first argument is a float or complex and the second argument is some 
> sort of integer array. This would be dispatched to some other helper 
> function instead of the normal pow_ufunc. In other words, optimize:
> 
> A**2, A**2.0, A**(2.0+0j), etc
> 
> and
> 
> A**array([1,2,3])
> 
> but not
> 
> A**array[1.0, 2.0, 3.0]
> 
> I think that this takes care of the optimization slowing down power for 
> general floats and optimizes the only array-array case that really matter.

I think this might still be a tiny bit dangerous despite not observing 
monotonicity problems and would be a bit more conservative and change it to:

optimize:

A**2, A**(2+0j), etc

and

A**array([1,2,3])

but not

A**array[1.0, 2.0, 3.0], A**2.0, A**(2.0+0j)

--
Gary R.

Re: [Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

[David M. C. <co...@ph...>]

On Wed, Feb 15, 2006 at 03:13:56PM -0700, Travis Oliphant wrote:
> Robert Kern wrote:
> I've started this process already.  I think a useful search order is
> 
> 1) next to current setup.py  --- os.getcwd() is probably better than 
> what I did (backing up the frame until you can't go back anymore and 
> getting the __file__ from that frame).  Incidentally, it looks like a 
> site.cfg present there is already copied to numpy/distutils on install 
> --- it's looks like its just not used for the numpy build itself.
> 
> 2) in the compilers "HOME" directory --- not sure how to implement that.

Have a look at distutils.dist for the Distribution.find_config_files
method. Also, the parse_config_files method reads the config options in
a way that keeps which filenames they come from.

> 3) in the system-wide directory (what is currently done --- except when 
> you are installing numpy that means it has to be in 
> numpy/distutils/site.cfg).

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

Re: [Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

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

David M. Cooke wrote:

>On Wed, Feb 15, 2006 at 03:45:33PM -0600, Robert Kern wrote:
>  
>
>>Travis Oliphant wrote:
>>
>>    
>>
>>>So, in particular, does this mean that it is read from (relative to the 
>>>location of the main setup.py file)
>>>
>>>numpy/distutils/site.cfg ??
>>>      
>>>
>>Yes. And in the case of scipy, it needs to be in the *installed* numpy/distutils
>>directory.
>>
>>    
>>
>>>Yes, that is a bad place.   We need some suggestions as to where 
>>>site.cfg should be read from.
>>>      
>>>
>>Next to the setup.py that *invokes* numpy.distutils. AFAICT using os.getcwd() in
>>system_info.py will give you this directory even if you start running the script
>>from a different directory (e.g. "python ~/svn/scipy/setup.py install"). I can
>>check this in if we agree that this is what we want.
>>    
>>
>
>As a note: Python's distutils looks for distutils.cfg in it's installed
>location (/usr/lib/python2.4/distutils or whatever), then in
>~/.pydistutils.cfg (or $HOME/pydistutils.cfg on non-Posix systems like
>Windows), then for setup.cfg in the current directory. Keys in later files
>override ones in earlier files.
>
>  
>
I think this is a good plan.  However, what I've started doesn't 
implement the over-riding process properly.  Anybody want to take a stab 
at that?   It would be nice if it could get into this next release.

-Travis

Re: [Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

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

Robert Kern wrote:

>Travis Oliphant wrote:
>
>  
>
>>So, in particular, does this mean that it is read from (relative to the 
>>location of the main setup.py file)
>>
>>numpy/distutils/site.cfg ??
>>    
>>
>
>Yes. And in the case of scipy, it needs to be in the *installed* numpy/distutils
>directory.
>
>  
>
>>Yes, that is a bad place.   We need some suggestions as to where 
>>site.cfg should be read from.
>>    
>>
>
>Next to the setup.py that *invokes* numpy.distutils. AFAICT using os.getcwd() in
>system_info.py will give you this directory even if you start running the script
>from a different directory (e.g. "python ~/svn/scipy/setup.py install"). I can
>check this in if we agree that this is what we want.
>
>  
>
I've started this process already.  I think a useful search order is

1) next to current setup.py  --- os.getcwd() is probably better than 
what I did (backing up the frame until you can't go back anymore and 
getting the __file__ from that frame).  Incidentally, it looks like a 
site.cfg present there is already copied to numpy/distutils on install 
--- it's looks like its just not used for the numpy build itself.

2) in the compilers "HOME" directory --- not sure how to implement that.

3) in the system-wide directory (what is currently done --- except when 
you are installing numpy that means it has to be in 
numpy/distutils/site.cfg).

I created a get_site_cfg() function in system_info where this searching 
can be done.   Feel free to change it as appropriate. 

-Travis

Re: [Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

[David M. C. <co...@ph...>]

On Wed, Feb 15, 2006 at 03:45:33PM -0600, Robert Kern wrote:
> Travis Oliphant wrote:
> 
> > So, in particular, does this mean that it is read from (relative to the 
> > location of the main setup.py file)
> > 
> > numpy/distutils/site.cfg ??
> 
> Yes. And in the case of scipy, it needs to be in the *installed* numpy/distutils
> directory.
> 
> > Yes, that is a bad place.   We need some suggestions as to where 
> > site.cfg should be read from.
> 
> Next to the setup.py that *invokes* numpy.distutils. AFAICT using os.getcwd() in
> system_info.py will give you this directory even if you start running the script
> from a different directory (e.g. "python ~/svn/scipy/setup.py install"). I can
> check this in if we agree that this is what we want.

As a note: Python's distutils looks for distutils.cfg in it's installed
location (/usr/lib/python2.4/distutils or whatever), then in
~/.pydistutils.cfg (or $HOME/pydistutils.cfg on non-Posix systems like
Windows), then for setup.cfg in the current directory. Keys in later files
override ones in earlier files.

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

[Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

[Robert K. <rob...@gm...>]

Travis Oliphant wrote:

> So, in particular, does this mean that it is read from (relative to the 
> location of the main setup.py file)
> 
> numpy/distutils/site.cfg ??

Yes. And in the case of scipy, it needs to be in the *installed* numpy/distutils
directory.

> Yes, that is a bad place.   We need some suggestions as to where 
> site.cfg should be read from.

Next to the setup.py that *invokes* numpy.distutils. AFAICT using os.getcwd() in
system_info.py will give you this directory even if you start running the script
from a different directory (e.g. "python ~/svn/scipy/setup.py install"). I can
check this in if we agree that this is what we want.

-- 
Robert Kern
rob...@gm...

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter

[Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

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

sk...@po... wrote:

>I'm trying to build numpy an svn sandbox (just updated a couple minutes
>ago).  If I grub around in numpy/distutils/system_info.py it says something
>about creating a site.cfg file with (for example) information about locating
>atlas.  It says nothing about where this file belongs.  I took a stab and
>placed it in my numpy source tree, right next to setup.py, with these lines:
>
>Failing all this, is there some way to build numpy/scipy without atlas?  At
>this point I just want the damn thing to build.  I'll worry about
>performance later (if at all).
>  
>

Yes.

You need to set the appropriate environment variables to 'None'

In particular, on my system (which is multithreaded and has a BLAS 
picked up from /usr/lib and an unthreaded ATLAS that the system will find)

export PTATLAS='None'
export ATLAS='None'
export BLAS='None'

did the trick.

-Travis

[Numpy-discussion] Re: [SciPy-user] Where do I put site.cfg? What should it contain? Can I build w/o atlas?

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

Robert Kern wrote:

>sk...@po... wrote:
>  
>
>>I'm trying to build numpy an svn sandbox (just updated a couple minutes
>>ago).  If I grub around in numpy/distutils/system_info.py it says something
>>about creating a site.cfg file with (for example) information about locating
>>atlas.  It says nothing about where this file belongs. 
>>    
>>
>
>Sure it does. "The file 'site.cfg' in the same directory as this module is read
>for configuration options." I think it's a really bad place for it to be, but
>that is the state of affairs right now.
>
>  
>
So, in particular, does this mean that it is read from (relative to the 
location of the main setup.py file)

numpy/distutils/site.cfg ??

Yes, that is a bad place.   We need some suggestions as to where 
site.cfg should be read from.

I think you can set the environment variable ATLAS to 'None' and it will 
ignore ATLAS...  I believe this is true of any of the configuration 
options.

-Travis

Re: [Numpy-discussion] numpy/scipy transition

[Christopher B. <Chr...@no...>]

Arnd Baecker wrote:
> - we need scipy.xplt
>   (matplotlib is still no option at this point)

Why not? just curious.

-Chris

-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

[Numpy-discussion] Numeric.identify(4) failure

[Randall H. <vi...@ch...>]

Here on 64-bit Linux, I get strange Python errors with some Numeric functions
(see below).  Has this been fixed in more recent versions of Numeric?  Currently
running python-numeric-23.7-3 which is what comes with SuSE 9.3.

Thanks,

Randall

Python 2.4 (#1, Mar 22 2005, 18:42:42)
[GCC 3.3.5 20050117 (prerelease) (SUSE Linux)] on linux2
Type "help", "copyright", "credits" or "license" for more information.

 import Numeric
 Numeric.identity(4)

Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "/usr/lib64/python2.4/site-packages/Numeric/Numeric.py", line 604, in
identity
    return resize(array([1]+n*[0],typecode=typecode), (n,n))
  File "/usr/lib64/python2.4/site-packages/Numeric/Numeric.py", line 398, in resize
    return reshape(a, new_shape)
ValueError: total size of new array must be unchanged

Re: [Numpy-discussion] possible bug in dtype for records

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

Stefan van der Walt wrote:

>Using
>
>In [3]: numpy.__version__
>Out[3]: '0.9.5.2024'
>
>I see the following:
>
>In [4]: import numpy as N
>
>In [5]: ctype = N.dtype({'names': ('x', 'y', 'z'), 'formats' : [N.float32, N.float32, N.float32]})
>
>In [6]: ctype
>Out[6]: dtype([('x', '<f4'), ('y', '<f4'), ('z', '<f4')])
>
>In [7]: N.array([(1,2,3), (4,5,6)], dtype=ctype)
>Segmentation fault
>
>However, when I use a mutable list for defining dtype, i.e.
>
>'names': ['x', 'y', 'z'] instead of
>'names': ('x', 'y', 'z')
>
>it works fine.
>
>Is this expected behaviour?
>  
>

Got it.  Basically, the VOID_setitem code expected the special -1 entry 
in the fields dictionary to be a list, but it was never enforced.  So, 
when you entered a tuple like this problems arose.

I changed it in SVN so that the -1 entry is always a tuple (you can't 
change the fields names anyway unless you define a new data-type) so its 
more aptly described as a tuple.

Thanks again for finding this...

-Travis

Re: [Numpy-discussion] possible bug in dtype for records

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

Stefan van der Walt wrote:

>Using
>
>In [3]: numpy.__version__
>Out[3]: '0.9.5.2024'
>
>I see the following:
>
>In [4]: import numpy as N
>
>In [5]: ctype = N.dtype({'names': ('x', 'y', 'z'), 'formats' : [N.float32, N.float32, N.float32]})
>
>In [6]: ctype
>Out[6]: dtype([('x', '<f4'), ('y', '<f4'), ('z', '<f4')])
>
>In [7]: N.array([(1,2,3), (4,5,6)], dtype=ctype)
>Segmentation fault
>  
>
A segmentation fault is never expected behavior.  Thanks for pointing 
this out.  I'll see if I can figure out what is wrong.

-Travis

[Numpy-discussion] possible bug in dtype for records

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

Using

In [3]: numpy.__version__
Out[3]: '0.9.5.2024'

I see the following:

In [4]: import numpy as N

In [5]: ctype =3D N.dtype({'names': ('x', 'y', 'z'), 'formats' : [N.float=
32, N.float32, N.float32]})

In [6]: ctype
Out[6]: dtype([('x', '<f4'), ('y', '<f4'), ('z', '<f4')])

In [7]: N.array([(1,2,3), (4,5,6)], dtype=3Dctype)
Segmentation fault

However, when I use a mutable list for defining dtype, i.e.

'names': ['x', 'y', 'z'] instead of
'names': ('x', 'y', 'z')

it works fine.

Is this expected behaviour?

Regards
St=E9fan

Re: [Numpy-discussion] numpy/scipy transition

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

Arnd Baecker wrote:

>Reasons not to switch
>- there is no enthought edition yet (right?)
>- there are only packages for a few platforms/distribution
>- we need scipy.xplt
>  (matplotlib is still no option at this point)
>
>Discussion/Background:
>
>To us the two main show-stoppers are scipy.xplt and the question
>about an Enthought Edition  for Windows.
>For the Pool of PCs where the tutorial groups are to be held,
>it won't be a problem to install numpy/scipy in such
>a way that scipy.sandbox.xplt is visible as scipy.xplt
>(at least I hope).
>However, the students will either have windows (around 80%)
>or Linux at home. For windows users we have used
>the Enthought Edition (http://code.enthought.com/enthon/)
>and linux users were pointed to
>available packages for their machines or to install Numeric/scipy
>themselves.
>
>  
>
As long as there are binaries for all the packages.  Just having a list 
of Windows installers can also work.    Were you using all of what is in 
the Enthon edition?

>Concerning xplt another option might be to
>install scipy.sandbox.xplt in such a way
>that a `import xplt` would work. If that is possible we could
>try to supply `xplt` separately for some of the distributions,
>and maybe also for windows (which I don't use, so I have
>no idea how difficult that would be).
>  
>
I don't think that would be hard at all.  You can just run python 
setup.py bdist_inst from within the sandbox/xplt directory and get a 
windows installer.  

>If something like this was possible, the main question is
>whether a new enthon distribution with new numpy/scipy/ipython
>and all the other niceties of mayavi/VTK/wxPython/....
>will come out in the near future?
>  
>
I have no idea about that one.   But, it sounds like the guy (Joe) at 
Enthought who did most of the work on the Enthon distribution is no 
longer as available for them, so I'm not sure...

-Travis

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

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

David M. Cooke wrote:

>[changed subject to reflect this better]
>
>Tim Hochberg <tim...@co...> writes:
>  
>
>>David M. Cooke wrote:
>>    
>>
>>>Tim Hochberg <tim...@co...> writes:
>>>      
>>>
>>>>David M. Cooke wrote:
>>>> 2, Don't distinguish between scalars and arrays -- that just makes
>>>>things harder to explain.
>>>>        
>>>>
>>>Makes the optimizations better, though.
>>>
>>>
>>>      
>>>
>>Ah, Because you can hoist all the checks for what type of optimization
>>to do, if any, out of the core loop, right? That's a good point. Still
>>I'm not keen on a**b having different performance *and* different
>>results depending on whether b is a scalar or matrix. The first thing
>>to do is to measure how much overhead doing the optimization element
>>by element is going to add. Assuming that it's signifigant that leaves
>>us with the familiar dilema: fast, simple or general purpose; pick any
>>two.
>>
>>1. Do what I've proposed: optimize things at the c_pow level. This is
>>general purpose and relatively simple to implement (since we can steal
>>most of the code from complexobject.c). It may have a signifigant
>>speed penalty versus 2 though:
>>
>>2. Do what you've proposed: optimize things at the ufunc level. This
>>fast and relatively simple to implement. It's more limited in scope
>>and a bit harder to explain than 2.
>>
>>3. Do both. This is straightforward, but adds a bunch of extra code
>>paths with all the attendant required testing and possibility for
>>bugs. So, fast, general purpose, but not simple.
>>    
>>
>
>Start with #1, then try #2. The problem with #2 is that you still have
>to include #1: if you're doing x**y when y is an array, then you have
>do if (y==2) etc. checks in your inner loop anyways. In that case, you
>might as well do it in nc_pow. At that point, it may be better to move
>the #1 optimization to the level of x.__pow__ (see below).
>  
>
OK.

>>>The speed comes from the inlining of how to do the power _outside_ the
>>>inner loop. The reason x**2, etc. are slower currently is there is a
>>>function call in the inner loop. Your's and mine C library's pow()
>>>function mostly likely does something like I have above, for a single
>>>case: pow(x, 2.0) is calculated as x*x. However, each time through it
>>>has do decide _how_ to do it.
>>>
>>>
>>>      
>>>
>>Part of our difference in perspective comes from the fact that I've
>>just been staring at the guts of complex power. In this case you
>>always have function calls at present, even for s*s. (At least I'm
>>fairly certain that doesn't get inlined although I haven't checked).
>>Since much of the work I do is with complex matrices,  it's
>>appropriate that I focus on this.
>>    
>>
>
>Ah, ok, now things are clicking. Complex power is going to be harder,
>because making sure that going from x**2.001 to x**2 doesn't do some
>funny complex branch cut stuff (I work in reals all the time :-)
>  
>
We're always dealing with the principle branch though, so probably we 
can just ignore any branch cut issues. We'll see I suppose.

>For the real numbers, these type of optimizations *are* a big win, and
>don't have the same type of continuity problems. I'll put them into numpy soon.
>  
>
Complex power is something like thirty times slower than s*s, so there 
is some room for optimization there. Some peril too though, as you note.

>  
>
>>Have you measured the effect of a function call on the speed here, or
>>is that just an educated guess. If it's an educated guess, it's
>>probably worth determining how much of speed hit the function call
>>actually causes.  I was going to try to get a handle on this by
>>comparing multiplication of Complex numbers (which requires a function
>>call plus more math), with multiplication of Floats which does not.
>>Perversly, the Complex multiplication came out marginally faster,
>>which is hard to explain whichever way you look at it.
>>
>>    
>>
>>>>>timeit.Timer("a*b", "from numpy import arange; a =
>>>>>          
>>>>>
>>arange(10000)+0j; b = arange(10000)+0j").time
>>it(10000)
>>3.2974959107959876
>>    
>>
>>>>>timeit.Timer("a*b", "from numpy import arange; a = arange(10000);
>>>>>b
>>>>>          
>>>>>
>>= arange(10000)").timeit(100
>>00)
>>3.4541194481425919
>>    
>>
>
>You're not multiplying floats in the last one: you're multiplying
>integers. You either need to use a = arange(10000.0), or a =
>arange(10000.0, dtype=float) (to be more specific).
>  
>
Doh! Now that's embarassing. Well, when I actually measure float 
multiplication, it's between two and ten times as fast. For small arrays 
(N=1000) the difference is relatively small (3.5x), I assume because the 
setup overhead starts to dominate. For midsized array (N=10,000) the 
difference is larger (10x). For large arrays (N=100,000) the difference 
becomes small (2x). Presumably the memory is no longer fitting in the 
cache and I'm  having memory bandwidth issues.

>Your integer numbers are about 3x better than mine, though (difference
>in architecture, maybe? I'm on an Athlon64).
>  
>
I'm on a P4.

>  
>
>>>That's why I limited the optimization to scalar exponents: array
>>>exponents would mean it's about as slow as the pow() call, even if the
>>>checks were inlined into the loop. It would probably be even slower
>>>for the non-optimized case, as you'd check for the special exponents,
>>>then call pow() if it fails (which would likely recheck the exponents).
>>>
>>>
>>>      
>>>
>>Again, here I'm thinking of the complex case. In that case at least, I
>>don't think that the non-optimized case would take a noticeable speed
>>hit. I would put it into pow itself, which already special cases a==0
>>and b==0. For float pow it might, but that's already slow, so I doubt
>>that it would make much difference.
>>    
>>
>
>It does make a bit of difference with float pow: the general case
>slows down a bit.
>  
>
OK.  I was hoping that the difference would not be noticeable. I suspect 
that in the complex pow case, that will be the case since complex pow is 
so slow to begin with and since it already is doing some testing on the 
exponent.

>>>Maybe a simple way to add this is to rewrite x.__pow__() as something
>>>like the C equivalent of
>>>
>>>def __pow__(self, p):
>>>   if p is not a scalar:
>>>      return power(self, p)
>>>   elif p == 1:
>>>      return p
>>>   elif p == 2:
>>>      return square(self)
>>>   elif p == 3:
>>>      return cube(self)
>>>   elif p == 4:
>>>      return power_4(self)
>>>   elif p == 0:
>>>      return ones(self.shape, dtype=self.dtype)
>>>   elif p == -1:
>>>      return 1.0/self
>>>   elif p == 0.5:
>>>      return sqrt(self)
>>>
>>>and add ufuncs square, cube, power_4 (etc.).
>>>      
>>>
>>It sounds like we need to benchmark some stuff and see what we come up
>>with. One approach would be for each of us to implement this for one
>>time (say float) and see how the approaches compare speed wise. That's
>>not entirely fair as my approach will do much better at complex than
>>float I believe, but it's certainly easier.
>>    
>>
>
>The way the ufuncs are templated, we can split out the complex
>routines easily enough.
>
>Here's what I propose:
>
>- add a square() ufunc, where square(x) == x*x (but faster of course)
>- I'll fiddle with the floats
>- you fiddle with the complex numbers :-)
>
>I've created a new branch in svn, at
>http://svn.scipy.org/svn/numpy/branches/power_optimization
>to do this fiddling. The changes below I mention are all checked in as
>revision 2104 (http://projects.scipy.org/scipy/numpy/changeset/2104).
>
>I've added a square() ufunc to the power_optimization branch because
>I'd argue that it's probably *the* most common use of **. I've
>implemented it, and it's as fast as a*a for reals, and runs in 2/3 the
>time as a*a for complex (which makes sense: squaring a complex numbers
>has 3 real multiplications, while multiplying has 4 in the (simple)
>scheme [1]).
>
>At least with square(), there's no argument about continuity, as it
>only squares :-).
>  
>
Actually that's not entirely true,  This get's back to the odd 
inaccuracy that started this thread:

 >>> array([1234567j])**2
array([ -1.52415568e+12+0.00018665j])

If you special case this the extraneous imaginary value will vanish, but 
raising things to the 2.000001 power or 1.999999 power will still be off 
a similar amount. I played with this a bunch though and I couldn't come 
up with a plausible way for this to make things break. I suspect I could 
come up with an implausible one though.


[some time passes while I sleep and otherwise try to live a normal life....]

>The next step I'd suggest is special-casing x.__pow__, like I suggest
>above. We could just test for integer scalar exponents (0, 1, 2), and
>just special-case those (returning ones(), x.copy(), square(x)), and
>leave all the rest to power().
>  
>
As I've been thinking about this some more, I think the correct thing to 
do is not to mess with  the power ufuncs at all. Rather in x.__pow__ 
(since I don't know that there's anywhere else to do it), after the 
above checks check the types of the array and in the cases where the 
first argument is a float or complex and the second argument is some 
sort of integer array. This would be dispatched to some other helper 
function instead of the normal pow_ufunc. In other words, optimize:

A**2, A**2.0, A**(2.0+0j), etc

and

A**array([1,2,3])

but not

A**array[1.0, 2.0, 3.0]

I think that this takes care of the optimization slowing down power for 
general floats and optimizes the only array-array case that really matter.

>I've also checked in code to the power_optimization branch that
>special cases power(x, <scalar exponent>), or anytime the basic ufunc
>gets called with a stride of 0 for the exponent. It doesn't do complex
>x, so no problems on your side, but it's a good chunk faster for this
>case than what we've got now. One reason I'm also looking at adding
>square() is because my optimization of power() makes x**2 run (only)
>1.5 slower than x*x (and I can't for the life of me see where that 0.5
>is coming from! It should be 1.0 like square()!).
>  
>
I just checked out your branch and I'll fiddle with the complex stuff as 
I've got time. I've got relatives in town this week, so my extra cycles 
just dropped precipitously.

>[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.

-tim

Re: [Numpy-discussion] New release of NumPy coming

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

Travis Oliphant wrote:

>
> I'd like to make a new release of NumPy in the next day or two.   If 
> there are any outstanding issues, please let me know.
>
Just a datapoint: I just compiled a clean checkout here using VC7 on 
windows XP and it compiled and passed all tests.

-tim

Re: [Numpy-discussion] avoiding the matrix copy performance hit

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

Ed Schofield wrote:

>Travis Oliphant wrote:
>
>  
>
>>I think I originally tried to make mat *not* return a copy, but this
>>actually broke code in SciPy.  So, I left the default as it was as a
>>copy on input.  There is an *asmatrix* command that does not return a
>>copy...
>>    
>>
>
>All SciPy's unit tests actually pass with a default of copy=False in the
>matrix constructor.  So SciPy needn't be the blocker here.  I'd like to
>cast a vote for not copying by default, in the interests of efficiency
>and, as Bill Baxter argued, usefulness.
>  
>
I would like to cast a vote for keeping the behaviour the same. Note that:

mat([[1,2,3], [3,4,5]])

will always create a copy of its data by necesesity. Which means that 
changing the default copy=False means that some data will be copied, 
while others will not be, potentially leading to subtle bugs. I strongly 
disaprove of this sort of inconstent behaviour (don't get me started on 
reshape!).

In this situation, people should just use asmatrix.

-tim


>-- Ed
>
>
>-------------------------------------------------------
>This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
>for problems?  Stop!  Download the new AJAX search engine that makes
>searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
>http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642
>_______________________________________________
>Numpy-discussion mailing list
>Num...@li...
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>
>  
>