np.array symbol

[Nikolas Tezak]

Hi,

I am wondering if it were in principle possible to create subclass of
Symbol, that could handle array indices, so that at a later time
(after having performed symbolic calculations) I could just substitute
a numpy array in for that symbol.

e.g. (10 + a[:,0] *a[:,1]).subst(a, np.array([[...]])
Of course I could simply create symbols for individual array elements,
but I am working with a large number of coefficients (~1000+) for
which it would seem inefficient (at least to me) to create individual
symbols.

I am still trying to get a good feeling for the code and I will try to
implement it myself, but maybe some of you already know of some reason
why this would not work (efficiently).

Nikolas

[Aaron S. Meurer]

If I understand you correctly, it is related to this: http://code.google.com/p/sympy/issues/detail?id=16

Symbols with indices is something that would be nice to have implemented, but unfortunately isn't yet. So maybe the information on that page will help you.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>

[smichr]

On Jan 28, 3:44 am, Nikolas Tezak <Nikolas.Te...@gmx.de> wrote:
> Hi,
>
> I am wondering if it were in principle possible to create subclass of  
> Symbol, that could handle array indices, so that at a later time  
> (after having performed symbolic calculations) I could just substitute  
> a numpy array in for that symbol.
>
> e.g. (10 + a[:,0] *a[:,1]).subst(a, np.array([[...]])

Nikolas,

Not being familiar with numpy, I'm not sure what your expression
should expand to. Could you give me a concrete example with an array
of [1, 2, 3]?

/chris

[Nikolas Tezak]

Chris,

numpy has an Array class that adds to python multidimensional arrays
for various types and precisions.
It is suited for large (think 10^6+ entries) sets of data (python
lists aren't very efficient for this).

Numpy arrays support addition and multiplication, e.g.:

In [82]: a = np.array([0.,0.,1.,3.,2.])

adding a scalar to an array results in the addition of each array
element with the scalar

In [83]: a+1
Out[83]: array([ 1., 1., 2., 4., 3.])

adding or multiplying two arrays (necessarily of the same length)
result in element-wise addition or multiplication of the arrays'
contents.

For a project I am currently working on, I store a set of N dynamic,
i.e. time dependent variables in a numpy array where each row
corresponds to all variables at a given time (if you are familiar with
this: it is a quantum mechanical state vector, but that does not
really matter). I then calculate a time evolution for K timesteps for
my variables, resulting in an array of dimension (K,N) = (#rows. #cols).
I now would like to evaluate dynamical observables that are functions
of the state at each time. A very simple example might be the
occupation probability for a single state of the system which would be
given by the absolute square of one of the dynamic variables:

a = np.zeros((K,N), dtype = np.complex) # the array's shape is (K,N)
and its data type is complex

#calculate a[k,i] for all variables (i in range(0,N)) and times (k in
range(0,K))

p_0 = np.zeros(K) #empty 1d array for storage of the observable

for k in range(K):
p_0_t[k] = a[k,0]*a[k,0].conjugate()

#actually this may be written simpler (and calculated more
efficiently) as:

p_0_t = a[:,0]*a[:,0].conjugate()

The reason why I would like to do this is because I don't know the
right formula for my observables in terms of the dynamic state
variables in advance.
I am considering a variety of different systems for which even the
number of state variables (and more importantly their exact
interpretation) is not constant so I actually compute the right
expression of my observable in terms of the state variables on the fly
using a selfwritten and limited implementation of a symbolic algebra
(only addition and multiplication is allowed). This method is working
very well and nicely for small scale systems (N ~ 1000) but it does
get kind of hairy after that. I was therefore wondering if I could
compare its performance to something created with sympy, but so far I
am running into some problems.
I would therefore like to be able to work with expressions referencing
my array symbolically (Like the IndexedSymbol proposed under the link
Aaron posted) and only later actually substitute in the actual data.

The example you asked for would therefore be:

a = np.array([1, 2, 3])

x = IndexedSymbol('x')

(x**2 + 1).subs(x, a)

would result in an np.array containing:
np.array([2, 5, 10])

which would be consistent with the behaviour of numpy arrays (element-
wise powers and addition of an overall constant)
A bit more interesting would be:

(x[0:2]+10).sub(x,a)

which should now only result in a slice of the original shape:
-> np.array([11, 12])

Also if it could work such that the actual evaluation is performed
efficiently by numpy, then this expression already shows the potential
behind such a method.
The biggest problem I can currently see is that sympy requires
external objects that are substituted into an expression to be
sympified and this would be very impractical for a large array
containing millions of entries.
Any ideas concerning this would be greatly appreciated :)

Nikolas

[Nikolas Tezak]

This goes a lot further than what I would need. I would not even need
the array indices to be symbolic expressions (although it certainly
would be nice to have that option). I would only like to substitute a
very large numpy array into a symbolic expression that is then
completely determined and should be evaluated (as efficiently as
possible). I understand of course that this is a bit specialized.

Nikolas

[Vinzent Steinberg]

You can use sympy objects inside numpy arrays, see sympy/test_external/
test_numpy.py.

However, to use subs with numpy arrays you would have to implement
this operation first. Did you have actually a look at lambdify? It
supports numpy.

Vinzent