Re: [Numpy-discussion] Numpy 2.0 schedule

2011-01-26 Thread Ralf Gommers
On Wed, Jan 26, 2011 at 12:28 PM, Mark Wiebe mwwi...@gmail.com wrote:

 On Tue, Jan 25, 2011 at 5:18 PM, Charles R Harris 
 charlesr.har...@gmail.com wrote:



 On Tue, Jan 25, 2011 at 1:13 PM, Travis Oliphant 
 oliph...@enthought.comwrote:


  It may make sense for a NumPy 1.6 to come out in March / April in the
 interim.


 Pulling out the changes to attain backward compatibility isn't getting any
 easier. I'd rather shoot for 2.0 in June. What can the rest of us do to help
 move things along?



Focusing on 2.0 makes sense to me too. Besides that, March/April is bad
timing for me so someone else should volunteer to be the release manager if
we go for a 1.6.

I took a shot at fixing the ABI compatibility, and if PyArray_ArrFunc was
 the main issue, then that might be done.  An ABI compatible 1.6 with the
 datetime and half types should be doable, just some extensions might get
 confused if they encounter arrays made with the new data types.

 Even if you fixed the ABI incompatibility (I don't know enough about the
issue to confirm that), I'm not sure how much value there is in a release
with as main new feature two dtypes that are not going to work well with
scipy/other binaries compiled against 1.5.

Cheers,
Ralf
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] tril, triu, document/ implementation conflict

2011-01-26 Thread eat
Hi,

I just noticed a document/ implementation conflict with tril and triu.
According tril documentation it should return of same shape and data-type as

called. But this is not the case at least with dtype bool.

The input shape is referred as (M, N) in tril and triu, but as (N, M) in
tri.
Inconsistent?

Also I'm not very happy with the performance, at least dtype bool can be
accelerated as follows.

In []: M= ones((2000, 3000), dtype= bool)
In []: timeit triu(M)
10 loops, best of 3: 173 ms per loop
In []: timeit triu_(M)
10 loops, best of 3: 107 ms per loop

In []: M= asarray(M, dtype= int)
In []: timeit triu(M)
10 loops, best of 3: 160 ms per loop
In []: timeit triu_(M)
10 loops, best of 3: 163 ms per loop

In []: M= asarray(M, dtype= float)
In []: timeit triu(M)
10 loops, best of 3: 195 ms per loop
In []: timeit triu_(M)
10 loops, best of 3: 157 ms per loop

I have attached a crude 'fix' incase someone is interested.
Regards,
eat


twodim_base_fix.py
Description: Binary data
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] tril, triu, document/ implementation conflict

2011-01-26 Thread josef . pktd
On Wed, Jan 26, 2011 at 7:22 AM, eat e.antero.ta...@gmail.com wrote:
 Hi,

 I just noticed a document/ implementation conflict with tril and triu.
 According tril documentation it should return of same shape and data-type as
 called. But this is not the case at least with dtype bool.

 The input shape is referred as (M, N) in tril and triu, but as (N, M) in
 tri.
 Inconsistent?

 Also I'm not very happy with the performance, at least dtype bool can be
 accelerated as follows.

 In []: M= ones((2000, 3000), dtype= bool)
 In []: timeit triu(M)
 10 loops, best of 3: 173 ms per loop
 In []: timeit triu_(M)
 10 loops, best of 3: 107 ms per loop

 In []: M= asarray(M, dtype= int)
 In []: timeit triu(M)
 10 loops, best of 3: 160 ms per loop
 In []: timeit triu_(M)
 10 loops, best of 3: 163 ms per loop

 In []: M= asarray(M, dtype= float)
 In []: timeit triu(M)
 10 loops, best of 3: 195 ms per loop
 In []: timeit triu_(M)
 10 loops, best of 3: 157 ms per loop

 I have attached a crude 'fix' incase someone is interested.

You could open a ticket for this.

just one comment:
I don't think this is readable, especially if we only look at the
source of the function with np.source

out= mul(ge(so(ar(m.shape[0]), ar(m.shape[1])), -k), m)

from np.source(np.tri) with numpy 1.5.1
m = greater_equal(subtract.outer(arange(N), arange(M)),-k)

Josef


 Regards,
 eat
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] tril, triu, document/ implementation conflict

2011-01-26 Thread eat
Hi,

On Wed, Jan 26, 2011 at 2:35 PM, josef.p...@gmail.com wrote:

  On Wed, Jan 26, 2011 at 7:22 AM, eat e.antero.ta...@gmail.com wrote:
  Hi,
 
  I just noticed a document/ implementation conflict with tril and triu.
  According tril documentation it should return of same shape and data-type
 as
  called. But this is not the case at least with dtype bool.
 
  The input shape is referred as (M, N) in tril and triu, but as (N, M) in
  tri.
  Inconsistent?

Any comments about the names for rows and cols. I prefer (M, N).

  
  Also I'm not very happy with the performance, at least dtype bool can be
  accelerated as follows.
 
  In []: M= ones((2000, 3000), dtype= bool)
  In []: timeit triu(M)
  10 loops, best of 3: 173 ms per loop
  In []: timeit triu_(M)
  10 loops, best of 3: 107 ms per loop
 
  In []: M= asarray(M, dtype= int)
  In []: timeit triu(M)
  10 loops, best of 3: 160 ms per loop
  In []: timeit triu_(M)
  10 loops, best of 3: 163 ms per loop
 
  In []: M= asarray(M, dtype= float)
  In []: timeit triu(M)
  10 loops, best of 3: 195 ms per loop
  In []: timeit triu_(M)
  10 loops, best of 3: 157 ms per loop
 
  I have attached a crude 'fix' incase someone is interested.

 You could open a ticket for this.

 just one comment:
 I don't think this is readable, especially if we only look at the
 source of the function with np.source

 out= mul(ge(so(ar(m.shape[0]), ar(m.shape[1])), -k), m)

 from np.source(np.tri) with numpy 1.5.1
 m = greater_equal(subtract.outer(arange(N), arange(M)),-k)

I agree, thats why I called it crude. Before opening a ticket I'll try to
figure out if there exists somewhere in numpy .astype functionality, but not
copying if allready proper dtype.

Also I'm afraid that I can't produce sufficient testing.

Regards, eat


 Josef

 
  Regards,
  eat
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion
 
 
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy 2.0 schedule

2011-01-26 Thread David Cournapeau
On Wed, Jan 26, 2011 at 6:47 PM, Dag Sverre Seljebotn
da...@student.matnat.uio.no wrote:
 On 01/26/2011 02:05 AM, David wrote:
 On 01/26/2011 01:42 AM, Charles R Harris wrote:

 Hi All,

 Just thought it was time to start discussing a release schedule for
 numpy 2.0 so we have something to aim at. I'm thinking sometime in the
 period April-June might be appropriate. There is a lot coming with the
 next release: the Enthought's numpy refactoring, Mark's float16 and
 iterator work, and support for IronPython. How do things look to the
 folks involved in those projects?

 One thing which I was wondering about numpy 2.0: what's the story for
 the C-API compared to 1.x for extensions. Is it fundamentally different
 so that extensions will need to be rewritten ? I especially wonder about
 scipy and cython's codegen backend,


 For CPython, my understanding is that extensions that access struct
 fields directly without accessor macros need to be changed, but not much
 else. There's a backwards-compatability PyArray_* API for CPython.

 That doesn't work for .NET, but neither does anything else in C
 extensions. So in the SciPy port to .NET there's my efforts to replace
 f2py with fwrap/Cython, and many SciPy C extensions will be rewritten in
 Cython. These will use the Npy_* interface (or backwards-compatability
 PyArray_* wrappers in numpy.pxd, but these only work in Cython under
 .NET, not in C, due to typing issues (what is object and so on)).

Ok, good to know. A good test would be to continuously build numpy +
scipy on top of it ASAP. Do you think cython (or is it sage) could
donate some CPU resources on the cython CI server for numpy ? I could
spend some time to make that work,

cheers,

David

 Dag Sverre
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] How to tell if I succeeded to build numpy with amd, umfpack and lapack

2011-01-26 Thread Samuel John
Hi there!

I have successfully built numpy 1.5 on ubuntu lucid (32 for now).
I think I got ATLAS/lapack/BLAS support, and if I
  ldd linalg/lapack_lite.so
I see that my libptf77blas.so etc. are successfully linked. :-)

However, how to I find out, if (and where) libamd.a and libumfpack.a
have been found and (statically) linked.
As far as I understand, I they are not present, a fallback in pure
python is used, right?

Is there a recommended way, I can query against which libs numpy has
been built?
So I can be sure numpy uses my own compiled versions of libamd, lapack
and so forth.

And the fftw3 is no longer supported, I guess (even if it is still
mentioned in the site.cfg.example)

Bests,
 Samuel


-- 
Dipl.-Inform. Samuel John
- - - - - - - - - - - - - - - - - - - - - - - - -
PhD student, CoR-Lab(.de) and
Neuroinformatics Group, Faculty
of Technology, D33594 Bielefeld
in cooperation with the HONDA
Research Institute Europe GmbH


jabber: samuelj...@jabber.org
- - - - - - - - - - - - - - - - - - - - - - - - -
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy 2.0 schedule

2011-01-26 Thread Bruce Southey

On 01/25/2011 10:28 PM, Mark Wiebe wrote:
On Tue, Jan 25, 2011 at 5:18 PM, Charles R Harris 
charlesr.har...@gmail.com mailto:charlesr.har...@gmail.com wrote:




On Tue, Jan 25, 2011 at 1:13 PM, Travis Oliphant
oliph...@enthought.com mailto:oliph...@enthought.com wrote:


On Jan 25, 2011, at 10:42 AM, Charles R Harris wrote:

 Hi All,

 Just thought it was time to start discussing a release
schedule for numpy 2.0 so we have something to aim at. I'm
thinking sometime in the period April-June might be
appropriate. There is a lot coming with the next release: the
Enthought's numpy refactoring, Mark's float16 and iterator
work, and support for IronPython. How do things look to the
folks involved in those projects?


My suggestion is to do a 1.6 relatively soon, as the current trunk 
feels pretty stable to me, and it would be nice to release the 
features without having to go through the whole merging process.



I would target June / July at this point ;-)I know I
deserve a I told you so from Chuck --- I will take it.


How much remains to get done?


My changes probably make merging the refactor more challenging too.


There is a bit of work that Mark is doing that would be good
to include, also some modifications to the re-factoring that
will support better small array performance.


Not everything needs to go into first release as long as the
following releases are backward compatible. So the ABI needs it's
final form as soon as possible. Is it still in flux?


I would suggest it is - there are a number of things I think could be 
improved in it, and it would be nice to bake in the underlying support 
features to make lazy/deferred evaluation of array expressions possible.


It may make sense for a NumPy 1.6 to come out in March / April
in the interim.


Pulling out the changes to attain backward compatibility isn't
getting any easier. I'd rather shoot for 2.0 in June. What can the
rest of us do to help move things along?


I took a shot at fixing the ABI compatibility, and if PyArray_ArrFunc 
was the main issue, then that might be done.  An ABI compatible 1.6 
with the datetime and half types should be doable, just some 
extensions might get confused if they encounter arrays made with the 
new data types.


-Mark

I do understand that it may take time for the 'dust to settle' but there 
is the opportunity to implement aspects that may require 'significant' 
notification or least start the process for any appropriate changes.


So, would it be possible to start developing some strategic plan of the 
changes that will occur?


The type of things I think are in terms of:

1) Notifying/warning users of the API changes that will occur. I agree 
with Chuck that other 'eyes' need to see it.


2) Add any desired depreciation warnings but I do not know of any. 
Perhaps the files in numpy/oldnumeric and numpy/numarray - if these are 
still important then these should have a better home since both have not 
had a release since mid 2006.


3) Changes or reorganization of the namespace.
My personal one is my ticket 1051 (Renaming and removing NaN and related 
IEEE754 special cases):

http://projects.scipy.org/numpy/ticket/1051
Hopefully some of it will be applied.

4) Changes in functions.
Examples:
Ticket 1262 (genfromtxt: dtype should be None by default)
http://projects.scipy.org/numpy/ticket/1262

Tickets 465 and 518 related to the accumulator dtype argument issues 
because this topic keeps appearing on the list.

http://projects.scipy.org/numpy/ticket/518
http://projects.scipy.org/numpy/ticket/465
For example, perhaps changing the default arguments of mean in 
numpy/core/fromnumeric.py as that allows the old behavior to remain by 
changing the dtype argument:

Change:
def mean(a, axis=None, dtype=None, out=None):
To:
def mean(a, axis=None, dtype=float, out=None):

5) Adding any enhancement patches like median of Ticket 1213
http://projects.scipy.org/numpy/ticket/1213

Bruce
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] 3d plane to point cloud fitting using SVD

2011-01-26 Thread Huan Liu
Hi,

I just confirmed Stefan's answer on one of the examples in 
http://www.mathworks.co.jp/matlabcentral/newsreader/view_thread/262996

matlab:

  A = randn(100,2)*[2 0;3 0;-1 2]';
  A = A + randn(size(A))/3;
  [U,S,V] = svd(A);
  X = V(:,end)

python:

  from numpy import *
  A = random.randn(100,2)*mat([[2,3,-1],[0,0,2]])
  A = A + random.randn(100,3)/3.0
  u,s,vh = linalg.linalg.svd(A)
  v = vh.conj().transpose()
  print v[:,-1] 

It works! Thanks Peter for bringing this up and Stefan for answering!

Huan




___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] Numpy 2.0 schedule

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 2:23 AM, Ralf Gommers
ralf.gomm...@googlemail.comwrote:

 On Wed, Jan 26, 2011 at 12:28 PM, Mark Wiebe mwwi...@gmail.com wrote:

 On Tue, Jan 25, 2011 at 5:18 PM, Charles R Harris 
 charlesr.har...@gmail.com wrote:

 On Tue, Jan 25, 2011 at 1:13 PM, Travis Oliphant oliph...@enthought.com
  wrote:


  It may make sense for a NumPy 1.6 to come out in March / April in the
 interim.


 Pulling out the changes to attain backward compatibility isn't getting
 any easier. I'd rather shoot for 2.0 in June. What can the rest of us do to
 help move things along?



 Focusing on 2.0 makes sense to me too. Besides that, March/April is bad
 timing for me so someone else should volunteer to be the release manager if
 we go for a 1.6.


I think sooner than March/April might be a possibility.  I've gotten the ABI
working so this succeeds on my machine:

* Build SciPy against NumPy 1.5.1
* Build NumPy trunk
* Run NumPy trunk with the 1.5.1-built SciPy - all tests pass except for one
(PIL image resize, which tests all float types and half lacks the precisions
necessary)

I took a shot at fixing the ABI compatibility, and if PyArray_ArrFunc was
 the main issue, then that might be done.  An ABI compatible 1.6 with the
 datetime and half types should be doable, just some extensions might get
 confused if they encounter arrays made with the new data types.

 Even if you fixed the ABI incompatibility (I don't know enough about the
 issue to confirm that), I'm not sure how much value there is in a release
 with as main new feature two dtypes that are not going to work well with
 scipy/other binaries compiled against 1.5.


I've recently gotten the faster ufunc NEP implementation finished except for
generalized ufuncs, and most things work the same or faster with it. Below
are some timings of 1.5.1 vs the new_iterator branch.  In particular, the
overhead on small arrays hasn't gotten worse, but the output memory layout
speeds up some operations by a lot.

To exercise the iterator a bit, and try to come up with a better approach
than the generalized ufuncs, I came up with a new function, 'einsum' for the
Einstein summation convention.  I'll send another email about it, but it for
instance solves the problem discussed here:

http://mail.scipy.org/pipermail/numpy-discussion/2006-May/020506.html

as c = np.einsum('rij,rjk-rik', a, b)

-Mark

The timings:

In [1]: import numpy as np
In [2]: np.version.version
Out[2]: '1.5.1'
In [3]: a = np.arange(9.).reshape(3,3); b = a.copy()
In [4]: timeit a + b
10 loops, best of 3: 3.48 us per loop
In [5]: timeit 2 * a
10 loops, best of 3: 6.07 us per loop
In [6]: timeit np.sum(a)
10 loops, best of 3: 7.19 us per loop
In [7]: a = np.arange(100).reshape(100,100,100); b = a.copy()
In [8]: timeit a + b
100 loops, best of 3: 17.1 ms per loop
In [9]: a = np.arange(1920*1080*3).reshape(1080,1920,3).swapaxes(0,1)
In [10]: timeit a * a
1 loops, best of 3: 794 ms per loop

In [1]: import numpy as np
In [2]: np.version.version
Out[2]: '2.0.0.dev-c97e9d5'
In [3]: a = np.arange(9.).reshape(3,3); b = a.copy()
In [4]: timeit a + b
10 loops, best of 3: 3.24 us per loop
In [5]: timeit 2 * a
10 loops, best of 3: 6.12 us per loop
In [6]: timeit np.sum(a)
10 loops, best of 3: 6.6 us per loop
In [7]: a = np.arange(100).reshape(100,100,100); b = a.copy()
In [8]: timeit a + b
100 loops, best of 3: 17 ms per loop
In [9]: a = np.arange(1920*1080*3).reshape(1080,1920,3).swapaxes(0,1)
In [10]: timeit a * a
10 loops, best of 3: 116 ms per loop
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


[Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
I wrote a new function, einsum, which implements Einstein summation
notation, and I'd like comments/thoughts from people who might be interested
in this kind of thing.

In testing it, it is also faster than many of NumPy's built-in functions,
except for dot and inner.  At the bottom of this email you can find the
documentation blurb I wrote for it, and here are some timings:

In [1]: import numpy as np
In [2]: a = np.arange(25).reshape(5,5)

In [3]: timeit np.einsum('ii', a)
10 loops, best of 3: 3.45 us per loop
In [4]: timeit np.trace(a)
10 loops, best of 3: 9.8 us per loop

In [5]: timeit np.einsum('ii-i', a)
100 loops, best of 3: 1.19 us per loop
In [6]: timeit np.diag(a)
10 loops, best of 3: 7 us per loop

In [7]: b = np.arange(30).reshape(5,6)

In [8]: timeit np.einsum('ij,jk', a, b)
1 loops, best of 3: 11.4 us per loop
In [9]: timeit np.dot(a, b)
10 loops, best of 3: 2.8 us per loop

In [10]: a = np.arange(1.)

In [11]: timeit np.einsum('i-', a)
1 loops, best of 3: 22.1 us per loop
In [12]: timeit np.sum(a)
1 loops, best of 3: 25.5 us per loop

-Mark

The documentation:

einsum(subscripts, *operands, out=None, dtype=None, order='K',
casting='safe')

Evaluates the Einstein summation convention on the operands.

Using the Einstein summation convention, many common multi-dimensional
array operations can be represented in a simple fashion.  This function
provides a way compute such summations.

The best way to understand this function is to try the examples below,
which show how many common NumPy functions can be implemented as
calls to einsum.

The subscripts string is a comma-separated list of subscript labels,
where each label refers to a dimension of the corresponding operand.
Repeated subscripts labels in one operand take the diagonal.  For
example,
``np.einsum('ii', a)`` is equivalent to ``np.trace(a)``.

Whenever a label is repeated, it is summed, so ``np.einsum('i,i', a,
b)``
is equivalent to ``np.inner(a,b)``.  If a label appears only once,
it is not summed, so ``np.einsum('i', a)`` produces a view of ``a``
with no changes.

The order of labels in the output is by default alphabetical.  This
means that ``np.einsum('ij', a)`` doesn't affect a 2D array, while
``np.einsum('ji', a)`` takes its transpose.

The output can be controlled by specifying output subscript labels
as well.  This specifies the label order, and allows summing to be
disallowed or forced when desired.  The call ``np.einsum('i-', a)``
is equivalent to ``np.sum(a, axis=-1)``, and
``np.einsum('ii-i', a)`` is equivalent to ``np.diag(a)``.

It is also possible to control how broadcasting occurs using
an ellipsis.  To take the trace along the first and last axes,
you can do ``np.einsum('i...i', a)``, or to do a matrix-matrix
product with the left-most indices instead of rightmost, you can do
``np.einsum('ij...,jk...-ik...', a, b)``.

When there is only one operand, no axes are summed, and no output
parameter is provided, a view into the operand is returned instead
of a new array.  Thus, taking the diagonal as ``np.einsum('ii-i', a)``
produces a view.

Parameters
--
subscripts : string
Specifies the subscripts for summation.
operands : list of array_like
These are the arrays for the operation.
out : None or array
If provided, the calculation is done into this array.
dtype : None or data type
If provided, forces the calculation to use the data type specified.
Note that you may have to also give a more liberal ``casting``
parameter to allow the conversions.
order : 'C', 'F', 'A', or 'K'
Controls the memory layout of the output. 'C' means it should
be Fortran contiguous. 'F' means it should be Fortran contiguous,
'A' means it should be 'F' if the inputs are all 'F', 'C' otherwise.
'K' means it should be as close to the layout as the inputs as
is possible, including arbitrarily permuted axes.
casting : 'no', 'equiv', 'safe', 'same_kind', 'unsafe'
Controls what kind of data casting may occur.  Setting this to
'unsafe' is not recommended, as it can adversely affect
accumulations.
'no' means the data types should not be cast at all. 'equiv' means
only byte-order changes are allowed. 'safe' means only casts
which can preserve values are allowed. 'same_kind' means only
safe casts or casts within a kind, like float64 to float32, are
allowed.  'unsafe' means any data conversions may be done.

Returns
---
output : ndarray
The calculation based on the Einstein summation convention.

See Also

dot, inner, outer, tensordot


Examples


 a = np.arange(25).reshape(5,5)
 b = np.arange(5)
 c = np.arange(6).reshape(2,3)

 

Re: [Numpy-discussion] einsum

2011-01-26 Thread Joshua Holbrook
On Wed, Jan 26, 2011 at 11:27 AM, Mark Wiebe mwwi...@gmail.com wrote:
 I wrote a new function, einsum, which implements Einstein summation
 notation, and I'd like comments/thoughts from people who might be interested
 in this kind of thing.

This sounds really cool! I've definitely considered doing something
like this previously, but never really got around to seriously
figuring out any sensible API.

Do you have the source up somewhere? I'd love to try it out myself.

--Josh
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 1:36 PM, Joshua Holbrook josh.holbr...@gmail.comwrote:

 On Wed, Jan 26, 2011 at 11:27 AM, Mark Wiebe mwwi...@gmail.com wrote:
  I wrote a new function, einsum, which implements Einstein summation
  notation, and I'd like comments/thoughts from people who might be
 interested
  in this kind of thing.

 This sounds really cool! I've definitely considered doing something
 like this previously, but never really got around to seriously
 figuring out any sensible API.

 Do you have the source up somewhere? I'd love to try it out myself.


You can check out the new_iterator branch from here:

https://github.com/m-paradox/numpy

$ git clone https://github.com/m-paradox/numpy.git
Cloning into numpy...

-Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Joshua Holbrook
On Wed, Jan 26, 2011 at 12:48 PM, Mark Wiebe mwwi...@gmail.com wrote:
 On Wed, Jan 26, 2011 at 1:36 PM, Joshua Holbrook josh.holbr...@gmail.com
 wrote:

 On Wed, Jan 26, 2011 at 11:27 AM, Mark Wiebe mwwi...@gmail.com wrote:
  I wrote a new function, einsum, which implements Einstein summation
  notation, and I'd like comments/thoughts from people who might be
  interested
  in this kind of thing.

 This sounds really cool! I've definitely considered doing something
 like this previously, but never really got around to seriously
 figuring out any sensible API.

 Do you have the source up somewhere? I'd love to try it out myself.

 You can check out the new_iterator branch from here:
 https://github.com/m-paradox/numpy
 $ git clone https://github.com/m-paradox/numpy.git
 Cloning into numpy...
 -Mark


Thanks for the link!

How closely coupled is this new code with numpy's internals? That is,
could you factor it out into its own package? If so, then people could
have immediate use out of it without having to integrate it into numpy
proper.

--Josh
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 2:01 PM, Joshua Holbrook josh.holbr...@gmail.comwrote:

 snip
 How closely coupled is this new code with numpy's internals? That is,
 could you factor it out into its own package? If so, then people could
 have immediate use out of it without having to integrate it into numpy
 proper.


The code depends heavily on the iterator I wrote, and I think the idea
itself depends on having a good dynamic multi-dimensional array library.
 When the numpy-refactor branch is complete, this would be part of
libndarray, and could be used directly from C without depending on Python.

-Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Robert Kern
On Wed, Jan 26, 2011 at 16:43, Mark Wiebe mwwi...@gmail.com wrote:
 On Wed, Jan 26, 2011 at 2:01 PM, Joshua Holbrook josh.holbr...@gmail.com
 wrote:

 snip
 How closely coupled is this new code with numpy's internals? That is,
 could you factor it out into its own package? If so, then people could
 have immediate use out of it without having to integrate it into numpy
 proper.

 The code depends heavily on the iterator I wrote, and I think the idea
 itself depends on having a good dynamic multi-dimensional array library.
  When the numpy-refactor branch is complete, this would be part of
 libndarray, and could be used directly from C without depending on Python.

It think his real question is whether einsum() and the iterator stuff
can live in a separate module that *uses* a released version of numpy
rather than a development branch.

-- 
Robert Kern

I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth.
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Joshua Holbrook

 It think his real question is whether einsum() and the iterator stuff
 can live in a separate module that *uses* a released version of numpy
 rather than a development branch.

 --
 Robert Kern


Indeed, I would like to be able to install and use einsum() without
having to install another version of numpy. Even if it depends on
features of a new numpy, it'd be nice to have it be a separate module.

--Josh
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Hanno Klemm

Mark,

interesting idea. Given the fact that in 2-d euclidean metric, the  
Einstein summation conventions are only a way to write out  
conventional matrix multiplications, do you consider at some point to  
include a non-euclidean metric in this thing? (As you have in special  
relativity, for example)

Something along the lines:

eta = np.diag(-1,1,1,1)

a = np.array(1,2,3,4)
b = np.array(1,1,1,1)

such that

einsum('i,i', a,b, metric=eta) = -1 + 2 + 3 + 4

I don't know how useful it would be, just a thought,
Hanno


Am 26.01.2011 um 21:27 schrieb Mark Wiebe:

 I wrote a new function, einsum, which implements Einstein summation  
 notation, and I'd like comments/thoughts from people who might be  
 interested in this kind of thing.

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Gael Varoquaux
On Thu, Jan 27, 2011 at 12:18:30AM +0100, Hanno Klemm wrote:
 interesting idea. Given the fact that in 2-d euclidean metric, the  
 Einstein summation conventions are only a way to write out  
 conventional matrix multiplications, do you consider at some point to  
 include a non-euclidean metric in this thing? (As you have in special  
 relativity, for example)

In my experience, Einstein summation conventions are quite
incomprehensible for people who haven't studies relativity (they aren't
used much outside some narrow fields of physics). If you start adding
metrics, you'll make it even harder for people to follow.

My 2 cents,

Gaël
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 3:05 PM, Joshua Holbrook josh.holbr...@gmail.comwrote:

 
  It think his real question is whether einsum() and the iterator stuff
  can live in a separate module that *uses* a released version of numpy
  rather than a development branch.
 
  --
  Robert Kern
 

 Indeed, I would like to be able to install and use einsum() without
 having to install another version of numpy. Even if it depends on
 features of a new numpy, it'd be nice to have it be a separate module.

 --Josh


Ah, sorry for misunderstanding.  That would actually be very difficult, as
the iterator required a fair bit of fixes and adjustments to the core.  The
new_iterator branch should be 1.5 ABI compatible, if that helps.

-Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 3:18 PM, Hanno Klemm kl...@phys.ethz.ch wrote:


 Mark,

 interesting idea. Given the fact that in 2-d euclidean metric, the
 Einstein summation conventions are only a way to write out
 conventional matrix multiplications, do you consider at some point to
 include a non-euclidean metric in this thing? (As you have in special
 relativity, for example)

 Something along the lines:

 eta = np.diag(-1,1,1,1)

 a = np.array(1,2,3,4)
 b = np.array(1,1,1,1)

 such that

 einsum('i,i', a,b, metric=eta) = -1 + 2 + 3 + 4


This particular example is already doable as follows:

 eta = np.diag([-1,1,1,1])
 eta
array([[-1,  0,  0,  0],
   [ 0,  1,  0,  0],
   [ 0,  0,  1,  0],
   [ 0,  0,  0,  1]])
 a = np.array([1,2,3,4])
 b = np.array([1,1,1,1])
 np.einsum('i,j,ij', a, b, eta)
8

I think that's right, did I understand you correctly?

Cheers,
Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Hanno Klemm


Am 27.01.2011 um 00:29 schrieb Mark Wiebe:

On Wed, Jan 26, 2011 at 3:18 PM, Hanno Klemm kl...@phys.ethz.ch  
wrote:


Mark,

interesting idea. Given the fact that in 2-d euclidean metric, the
Einstein summation conventions are only a way to write out
conventional matrix multiplications, do you consider at some point to
include a non-euclidean metric in this thing? (As you have in special
relativity, for example)

Something along the lines:

eta = np.diag(-1,1,1,1)

a = np.array(1,2,3,4)
b = np.array(1,1,1,1)

such that

einsum('i,i', a,b, metric=eta) = -1 + 2 + 3 + 4

This particular example is already doable as follows:

 eta = np.diag([-1,1,1,1])
 eta
array([[-1,  0,  0,  0],
   [ 0,  1,  0,  0],
   [ 0,  0,  1,  0],
   [ 0,  0,  0,  1]])
 a = np.array([1,2,3,4])
 b = np.array([1,1,1,1])
 np.einsum('i,j,ij', a, b, eta)
8

I think that's right, did I understand you correctly?

Cheers,
Mark


Yes, that's what I had in mind. Thanks.


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Benjamin Root
On Wednesday, January 26, 2011, Gael Varoquaux
gael.varoqu...@normalesup.org wrote:
 On Thu, Jan 27, 2011 at 12:18:30AM +0100, Hanno Klemm wrote:
 interesting idea. Given the fact that in 2-d euclidean metric, the
 Einstein summation conventions are only a way to write out
 conventional matrix multiplications, do you consider at some point to
 include a non-euclidean metric in this thing? (As you have in special
 relativity, for example)

 In my experience, Einstein summation conventions are quite
 incomprehensible for people who haven't studies relativity (they aren't
 used much outside some narrow fields of physics). If you start adding
 metrics, you'll make it even harder for people to follow.

 My 2 cents,

 Gaël


Just to dispel the notion that Einstein notation is only used in the
study of relativity, I can personally attest that Einstein notation is
used in the field of fluid dynamics and some aspects of meteorology.
This is really a neat idea and I support the idea of packaging it as a
separate module.

Ben Root
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Joshua Holbrook
 Ah, sorry for misunderstanding.  That would actually be very difficult,
 as the iterator required a fair bit of fixes and adjustments to the core.
 The new_iterator branch should be 1.5 ABI compatible, if that helps.

I see. Perhaps the fixes and adjustments can/should be included with
numpy standard, even if the Einstein notation package is made a
separate module.

 Just to dispel the notion that Einstein notation is only used in the
 study of relativity, I can personally attest that Einstein notation is
 used in the field of fluid dynamics and some aspects of meteorology.

Einstein notation is also used in solid mechanics.

--Josh
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread T J
On Wed, Jan 26, 2011 at 5:02 PM, Joshua Holbrook
josh.holbr...@gmail.com wrote:
 Ah, sorry for misunderstanding.  That would actually be very difficult,
 as the iterator required a fair bit of fixes and adjustments to the core.
 The new_iterator branch should be 1.5 ABI compatible, if that helps.

 I see. Perhaps the fixes and adjustments can/should be included with
 numpy standard, even if the Einstein notation package is made a
 separate module.

snip
 Indeed, I would like to be able to install and use einsum() without
 having to install another version of numpy. Even if it depends on
 features of a new numpy, it'd be nice to have it be a separate module.

I don't really understand the desire to have this single function
exist in a separate package.  If it requires the new version of NumPy,
then you'll have to install/upgrade either way...and if it comes as
part of that new NumPy, then you are already set.  Doesn't a separate
package complicate things unnecessarily?  It make sense to me if
einsum consisted of many functions (such as Bottleneck).
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread josef . pktd
On Wed, Jan 26, 2011 at 7:35 PM, Benjamin Root ben.r...@ou.edu wrote:
 On Wednesday, January 26, 2011, Gael Varoquaux
 gael.varoqu...@normalesup.org wrote:
 On Thu, Jan 27, 2011 at 12:18:30AM +0100, Hanno Klemm wrote:
 interesting idea. Given the fact that in 2-d euclidean metric, the
 Einstein summation conventions are only a way to write out
 conventional matrix multiplications, do you consider at some point to
 include a non-euclidean metric in this thing? (As you have in special
 relativity, for example)

 In my experience, Einstein summation conventions are quite
 incomprehensible for people who haven't studies relativity (they aren't
 used much outside some narrow fields of physics). If you start adding
 metrics, you'll make it even harder for people to follow.

 My 2 cents,

 Gaël


 Just to dispel the notion that Einstein notation is only used in the
 study of relativity, I can personally attest that Einstein notation is
 used in the field of fluid dynamics and some aspects of meteorology.
 This is really a neat idea and I support the idea of packaging it as a
 separate module.

So, if I read the examples correctly we finally get dot along an axis

np.einsum('ijk,ji-', a, b)
np.einsum('ijk,jik-k', a, b)

or something like this.

the notation might require getting used to but it doesn't look worse
than figuring out what tensordot does.

The only disadvantage I see, is that choosing the axes to operate on
in a program or function requires string manipulation.

Josef



 Ben Root
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Jonathan Rocher
Nice function, and wonderful that it speeds some tasks up.

some feedback: the following notation is a little counter intuitive to me:
 np.einsum('i...-', a)
array([50, 55, 60, 65, 70])
 np.sum(a, axis=0)
array([50, 55, 60, 65, 70])
 Since there is nothing after the -, I expected a scalar not a vector. I
might suggest 'i...-...'

Just noticed also a typo in the doc:
 order : 'C', 'F', 'A', or 'K'
Controls the memory layout of the output. 'C' means it should
be Fortran contiguous. 'F' means it should be Fortran contiguous,
should be changed to
order : 'C', 'F', 'A', or 'K'
Controls the memory layout of the output. 'C' means it should
be C contiguous. 'F' means it should be Fortran contiguous,


Hope this helps,
Jonathan

On Wed, Jan 26, 2011 at 2:27 PM, Mark Wiebe mwwi...@gmail.com wrote:

 I wrote a new function, einsum, which implements Einstein summation
 notation, and I'd like comments/thoughts from people who might be interested
 in this kind of thing.

 In testing it, it is also faster than many of NumPy's built-in functions,
 except for dot and inner.  At the bottom of this email you can find the
 documentation blurb I wrote for it, and here are some timings:

 In [1]: import numpy as np
 In [2]: a = np.arange(25).reshape(5,5)

 In [3]: timeit np.einsum('ii', a)
 10 loops, best of 3: 3.45 us per loop
 In [4]: timeit np.trace(a)
 10 loops, best of 3: 9.8 us per loop

 In [5]: timeit np.einsum('ii-i', a)
 100 loops, best of 3: 1.19 us per loop
 In [6]: timeit np.diag(a)
 10 loops, best of 3: 7 us per loop

 In [7]: b = np.arange(30).reshape(5,6)

 In [8]: timeit np.einsum('ij,jk', a, b)
 1 loops, best of 3: 11.4 us per loop
 In [9]: timeit np.dot(a, b)
 10 loops, best of 3: 2.8 us per loop

 In [10]: a = np.arange(1.)

 In [11]: timeit np.einsum('i-', a)
 1 loops, best of 3: 22.1 us per loop
 In [12]: timeit np.sum(a)
 1 loops, best of 3: 25.5 us per loop

 -Mark

 The documentation:

 einsum(subscripts, *operands, out=None, dtype=None, order='K',
 casting='safe')

 Evaluates the Einstein summation convention on the operands.

 Using the Einstein summation convention, many common multi-dimensional
 array operations can be represented in a simple fashion.  This function
 provides a way compute such summations.

 The best way to understand this function is to try the examples below,
 which show how many common NumPy functions can be implemented as
 calls to einsum.

 The subscripts string is a comma-separated list of subscript labels,
 where each label refers to a dimension of the corresponding operand.
 Repeated subscripts labels in one operand take the diagonal.  For
 example,
 ``np.einsum('ii', a)`` is equivalent to ``np.trace(a)``.

 Whenever a label is repeated, it is summed, so ``np.einsum('i,i', a,
 b)``
 is equivalent to ``np.inner(a,b)``.  If a label appears only once,
 it is not summed, so ``np.einsum('i', a)`` produces a view of ``a``
 with no changes.

 The order of labels in the output is by default alphabetical.  This
 means that ``np.einsum('ij', a)`` doesn't affect a 2D array, while
 ``np.einsum('ji', a)`` takes its transpose.

 The output can be controlled by specifying output subscript labels
 as well.  This specifies the label order, and allows summing to be
 disallowed or forced when desired.  The call ``np.einsum('i-', a)``
 is equivalent to ``np.sum(a, axis=-1)``, and
 ``np.einsum('ii-i', a)`` is equivalent to ``np.diag(a)``.

 It is also possible to control how broadcasting occurs using
 an ellipsis.  To take the trace along the first and last axes,
 you can do ``np.einsum('i...i', a)``, or to do a matrix-matrix
 product with the left-most indices instead of rightmost, you can do
 ``np.einsum('ij...,jk...-ik...', a, b)``.

 When there is only one operand, no axes are summed, and no output
 parameter is provided, a view into the operand is returned instead
 of a new array.  Thus, taking the diagonal as ``np.einsum('ii-i', a)``
 produces a view.

 Parameters
 --
 subscripts : string
 Specifies the subscripts for summation.
 operands : list of array_like
 These are the arrays for the operation.
 out : None or array
 If provided, the calculation is done into this array.
 dtype : None or data type
 If provided, forces the calculation to use the data type specified.
 Note that you may have to also give a more liberal ``casting``
 parameter to allow the conversions.
 order : 'C', 'F', 'A', or 'K'
 Controls the memory layout of the output. 'C' means it should
 be Fortran contiguous. 'F' means it should be Fortran contiguous,
 'A' means it should be 'F' if the inputs are all 'F', 'C'
 otherwise.
 'K' means it should be as close to the layout as the inputs 

Re: [Numpy-discussion] Numpy 2.0 schedule

2011-01-26 Thread Charles R Harris
On Wed, Jan 26, 2011 at 1:10 PM, Mark Wiebe mwwi...@gmail.com wrote:

 On Wed, Jan 26, 2011 at 2:23 AM, Ralf Gommers ralf.gomm...@googlemail.com
  wrote:

 On Wed, Jan 26, 2011 at 12:28 PM, Mark Wiebe mwwi...@gmail.com wrote:

 On Tue, Jan 25, 2011 at 5:18 PM, Charles R Harris 
 charlesr.har...@gmail.com wrote:

 On Tue, Jan 25, 2011 at 1:13 PM, Travis Oliphant 
 oliph...@enthought.com wrote:


  It may make sense for a NumPy 1.6 to come out in March / April in the
 interim.


 Pulling out the changes to attain backward compatibility isn't getting
 any easier. I'd rather shoot for 2.0 in June. What can the rest of us do to
 help move things along?



 Focusing on 2.0 makes sense to me too. Besides that, March/April is bad
 timing for me so someone else should volunteer to be the release manager if
 we go for a 1.6.


 I think sooner than March/April might be a possibility.  I've gotten the
 ABI working so this succeeds on my machine:


If we go with a 1.6 I have some polynomial stuff I want to put in, probably
a weekend or two of work, and there are tickets and pull requests to look
through, so to me March-April looks like a good time. It sounds like Ralf
has stuff scheduled for the rest of the spring after the scipy release.
IIRC, there was at least one other person interested in managing a release
when David left for Silveregg, do we have any volunteers for a 1.6?

If we do go for 1.6 I would like to keep 2.0 in sight. If datetime, the new
iterator, einsum, and float16 are in 1.6 then 2.0 looks more like a cleanup
the library/inteface and support IronPython release and there isn't as much
pressure to get it out soon. Also it is important to get the ABI right so we
don't need to change it again soon and doing that might take a bit of trial
and error. Does September seem reasonable?

* Build SciPy against NumPy 1.5.1
 * Build NumPy trunk
 * Run NumPy trunk with the 1.5.1-built SciPy - all tests pass except for
 one (PIL image resize, which tests all float types and half lacks the
 precisions necessary)

 I took a shot at fixing the ABI compatibility, and if PyArray_ArrFunc was
 the main issue, then that might be done.  An ABI compatible 1.6 with the
 datetime and half types should be doable, just some extensions might get
 confused if they encounter arrays made with the new data types.

 Even if you fixed the ABI incompatibility (I don't know enough about the
 issue to confirm that), I'm not sure how much value there is in a release
 with as main new feature two dtypes that are not going to work well with
 scipy/other binaries compiled against 1.5.


 I've recently gotten the faster ufunc NEP implementation finished except
 for generalized ufuncs, and most things work the same or faster with
 it. Below are some timings of 1.5.1 vs the new_iterator branch.  In
 particular, the overhead on small arrays hasn't gotten worse, but the output
 memory layout speeds up some operations by a lot.


snip

Chuck



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 6:41 PM, Jonathan Rocher jroc...@enthought.comwrote:

 Nice function, and wonderful that it speeds some tasks up.

 some feedback: the following notation is a little counter intuitive to me:
  np.einsum('i...-', a)
 array([50, 55, 60, 65, 70])
  np.sum(a, axis=0)
 array([50, 55, 60, 65, 70])
  Since there is nothing after the -, I expected a scalar not a vector. I
 might suggest 'i...-...'


Hmm, the dimension that's left is a a broadcast dimension, and the dimension
labeled 'i' did go away.  I suppose disallowing the empty output string and
forcing a '...' is reasonable.  Would disallowing broadcasting by default be
a good approach?  Then,
einsum('ii-i', a) would only except two dimensional inputs, and you would
have to specify einsum('...ii-...i', a) to get the current default behavior
for it.

Just noticed also a typo in the doc:

  order : 'C', 'F', 'A', or 'K'
 Controls the memory layout of the output. 'C' means it should
 be Fortran contiguous. 'F' means it should be Fortran contiguous,
 should be changed to
 order : 'C', 'F', 'A', or 'K'
 Controls the memory layout of the output. 'C' means it should
 be C contiguous. 'F' means it should be Fortran contiguous,


Thanks,
Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 5:23 PM, josef.p...@gmail.com wrote:

 snip
 So, if I read the examples correctly we finally get dot along an axis

 np.einsum('ijk,ji-', a, b)
 np.einsum('ijk,jik-k', a, b)

 or something like this.

 the notation might require getting used to but it doesn't look worse
 than figuring out what tensordot does.


I thought of various extensions to the notation, but the idea is tricky
enough as is I think. Decoding a regex-like syntax probably wouldn't help.

The only disadvantage I see, is that choosing the axes to operate on
 in a program or function requires string manipulation.


One possibility would be for the Python exposure to accept lists or tuples
of integers.  The subscript 'ii' could be [(0,0)], and 'ij,jk-ik' could be
[(0,1), (1,2), (0,2)].  Internally it would convert this directly to a
C-string to pass to the API function.

-Mark
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Joshua Holbrook

 The only disadvantage I see, is that choosing the axes to operate on
 in a program or function requires string manipulation.


 One possibility would be for the Python exposure to accept lists or tuples
 of integers.  The subscript 'ii' could be [(0,0)], and 'ij,jk-ik' could be
 [(0,1), (1,2), (0,2)].  Internally it would convert this directly to a
 C-string to pass to the API function.
 -Mark


What if you made objects i, j, etc. such that i*j = (0, 1) and
etcetera? Maybe you could generate them with something like (i, j, k)
= einstein((1, 2, 3)) .

Feel free to disregard me since I haven't really thought too hard
about things and might not even really understand what the problem is
*anyway*. I'm just trying to help brainstorm. :)

--Josh
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] einsum

2011-01-26 Thread Mark Wiebe
On Wed, Jan 26, 2011 at 8:29 PM, Joshua Holbrook josh.holbr...@gmail.comwrote:

 
  The only disadvantage I see, is that choosing the axes to operate on
  in a program or function requires string manipulation.
 
 
  One possibility would be for the Python exposure to accept lists or
 tuples
  of integers.  The subscript 'ii' could be [(0,0)], and 'ij,jk-ik' could
 be
  [(0,1), (1,2), (0,2)].  Internally it would convert this directly to a
  C-string to pass to the API function.
  -Mark
 

 What if you made objects i, j, etc. such that i*j = (0, 1) and
 etcetera? Maybe you could generate them with something like (i, j, k)
 = einstein((1, 2, 3)) .

 Feel free to disregard me since I haven't really thought too hard
 about things and might not even really understand what the problem is
 *anyway*. I'm just trying to help brainstorm. :)


No worries. :)  The problem is that someone will probably want to
dynamically generate the axes to process in a loop, rather than having them
hardcoded beforehand.  For example, generalizing the diag function as
follows.  Within Python, creating lists and tuples is probably more natural.

-Mark

 def diagij(x, i, j):
... ss = 
... so = 
... # should error check i, j
... fill = ord('b')
... for k in range(x.ndim):
... if k in [i, j]:
... ss += 'a'
... else:
... ss += chr(fill)
... so += chr(fill)
... fill += 1
... ss += '-' + so + 'a'
... return np.einsum(ss, x)
...
 x = np.arange(3*3*3).reshape(3,3,3)
 diagij(x, 0, 1)
array([[ 0, 12, 24],
   [ 1, 13, 25],
   [ 2, 14, 26]])

 [np.diag(x[:,:,i]) for i in range(3)]
[array([ 0, 12, 24]), array([ 1, 13, 25]), array([ 2, 14, 26])]

 diagij(x, 1, 2)
array([[ 0,  4,  8],
   [ 9, 13, 17],
   [18, 22, 26]])
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] How to tell if I succeeded to build numpy with amd, umfpack and lapack

2011-01-26 Thread Paul Ivanov
Samuel John, on 2011-01-26 15:08,  wrote:
 Hi there!
 
 I have successfully built numpy 1.5 on ubuntu lucid (32 for now).
 I think I got ATLAS/lapack/BLAS support, and if I
   ldd linalg/lapack_lite.so
 I see that my libptf77blas.so etc. are successfully linked. :-)
 
 However, how to I find out, if (and where) libamd.a and libumfpack.a
 have been found and (statically) linked.
 As far as I understand, I they are not present, a fallback in pure
 python is used, right?
 
 Is there a recommended way, I can query against which libs numpy has
 been built?
 So I can be sure numpy uses my own compiled versions of libamd, lapack
 and so forth.

Hi Samuel,

take a look at numpy.show_config() and scipy.show_config() 

best,
-- 
Paul Ivanov
314 address only used for lists,  off-list direct email at:
http://pirsquared.org | GPG/PGP key id: 0x0F3E28F7 


signature.asc
Description: Digital signature
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion