Re: [sympy] Using a numpy time series to evalf a SymPy equation

[Aaron Meurer]

Use lambdify().

Aaron Meurer

On Wed, Jun 12, 2013 at 1:06 PM, Shawn Garbett <shawn....@gmail.com> wrote:
> I'm struggling with the most direct route to use a time series numpy array
> to get a time series back from a SymPy equation.
>
> Ideally I'd like to call something like sol.evalf(x), but that doesn't work.
> I have a loop structure that does it at present, but the continual
> reconstructing a dict causes the processor to heat up and the fans come on.
> I figure, there has to be a shorter / better / direct method. I'm at a loss
> to find it.
>
> Here's what the data looks like:
>
> In [212]: sol
> Out[212]: 1.0e-6*s4 + 1.0e-12*s6
>
> In [213]: x
> Out[213]:
> array([(1.0, 0.33, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0),
> (7.332991703456984e-07, 0.06344978343322456, 1.0,
> 1.9622076748233104e-10, 0.7332991703050197, 0.266668723224782,
> 3.1373201206573815e-05, 1.1995787763723961e-07),
> (6.853150238285165e-07, 0.015615017796943445, 1.0,
> 2.359623438104902e-09, 0.6853150237422331, 0.31460014764904676,
> 8.41433121893229e-05, 6.888955138586781e-07),
> ...,
> (7.447520789346351e-07, 0.00010069599712919727, 1.0,
> 1.5661950866577716, 0.7447520734699161, 0.2497844169697847,
> 0.005464700864489391, 0.00815509952413137),
> (7.446567577748946e-07, 0.0001007088880961995, 1.0,
> 1.5662440343737307, 0.7446567523064337, 0.24987599614231418,
> 0.005468443007273509, 0.008160817601891841),
> (7.445614740433947e-07, 0.00010072177731066373, 1.0,
> 1.5662930164082658, 0.7445614685711894, 0.24996753549603745,
> 0.005472187539255912, 0.008166538792436515)],
> dtype=[('s0', '<f8'), ('s1', '<f8'), ('s2', '<f8'), ('s3', '<f8'),
> ('s4', '<f8'), ('s5', '<f8'), ('s6', '<f8'), ('s7', '<f8')])
>
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy?hl=en-US.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Shawn Garbett]

On Wednesday, June 12, 2013 4:54:21 PM UTC-5, Aaron Meurer wrote:

Use lambdify().

I've just spent a couple hours fiddling with lambdify(). I can't seem to get it to work. First of all, lambdify requires that one specify the args for the first argument. I'm trying to write this to deal with sets of equations, that I don't know the args up front, and in fact they total number of variables will be large (10^26 is possible--but I've got a strategy to keep it under 10,000 if that happens). I could specify a single arg, but then that requires that the SymPy expression use numpy syntax for subscripting which I cannot get to work SymPy. So basically I'm lost on how lambdify helps me.

 

[Matthew Rocklin]

It may be that we don't understand your application given your example.  From your example I suspect you want to compute something like the following

In [1]: from sympy import *

In [2]: s4 = Symbol('s4')

In [3]: s6 = Symbol('s6')

In [4]: sol = 1.0e-6*s4 + 1.0e-12*s6

In [5]: sol

Out[5]: 1.0e-6*s4 + 1.0e-12*s6

In [6]: f = lambdify((s4, s6), sol)

In [7]: import numpy as np

In [8]: data_s4 = np.asarray([1,2,3,4])  # Some random data, presumably these are much longer

In [9]: data_s6 = np.asarray([1,2,3,4])

In [10]: f(data_s4, data_s6)

Out[10]: 

array([  1.00000100e-06,   2.00000200e-06,   3.00000300e-06,

         4.00000400e-06])

 

--

[Shawn Garbett]

Here's a simpler example:

from numpy import *

from sympy import *

x = array([(1.1, 2.2), (3.3, 4.4), (5.5, 6.6)], dtype=[('s0', '<f8'), ('s1', '<f8')])

eq = 2.0 * Symbol('s0') - Symbol('s1')

2.0*x['s0']-x['s1'] # Gives correct result at command prompt

eq.evalf(x) # Fails

The constraints are that the symbols in the equation are unknown, and there is a random number of them that can be huge, like thousands of variables as symbols. What is provided is the time course data x and the equation eq which contains matching variables. I.e., The code has no way of knowing beforehand what the variables are going to be, nor how many so a predefined lambda will not solve the problem.

[Aaron Meurer]

You can get all the variables in an expression using
expr.free_symbols. I hope that helps.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.

> Visit this group at http://groups.google.com/group/sympy.

[Shawn Garbett]

I figured out the problem. It's in evalf_symbol. 

def evalf_symbol(x, prec, options):

    val = options['subs'][x]

The x is the actual symbol and the numpy array provided is indexed by the string representation of the symbol.

Numpy can't index by a Symbol, ugh. 

In [23]: x = array([(1.1, 2.2), (3.3, 4.4), (5.5, 6.6)], dtype=[(Symbol('s0'), '<f8'), (Symbol('s1'), '<f8')])

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-23-3b59912775d8> in <module>()

----> 1 x = array([(1.1, 2.2), (3.3, 4.4), (5.5, 6.6)], dtype=[(Symbol('s0'), '<f8'), (Symbol('s1'), '<f8')])

TypeError: data type not understood

Could the evalf_symbol be modified to try the string representation as a fallback position, i.e. the path of least surprise?

def evalf_symbol(x, prec, options):

    try:

        val = options['subs'][x]

    except IndexError:

        val = options['subs'][x.name]

If this is an acceptable change, could someone point me in the direction of how to submit properly (e.g. ticket/ patch procedure)?

[Aaron Meurer]

https://github.com/sympy/sympy/wiki/development-workflow

Aaron Meurer

[Shawn Garbett]

Thanks Aaron, I'll follow that. However, that simple patch didn't work. It causes more problems up the evaluation stack for sympy. I'm going to see if I can monkey patch the numpy object, and if that works.

Shawn

[Shawn Garbett]

I ended up making a shim between the two.

from numpy import *

from sympy import *

x = array([(1.1, 2.2), (3.3, 4.4), (5.5, 6.6)], dtype=[('s0', '<f8'), ('s1', '<f8')])

eq = 2.0 * Symbol('s0') - sqrt(Symbol('s1'))

from collections import Mapping

class FilterNdarray(Mapping):

    def __init__(self, source, t):

        self.source = source

        self.t = t

    def __getitem__(self, key):

        return self.source[key.name][self.t]

    def __len__(self):

        return len(self.source)

    def __iter__(self):

        for key in self.source:

            yield key

    def set_time(self, t):

        self.t = t

        return self

y = FilterNdarray(x, 0)

[eq.evalf(subs=y.set_time(t)) for t in range(x.size)]

It's still horribly slow. I found someone who wrote an integrator, and they have an evalf that compiles it down to bytecode somehow. I'm going to pursue that, and it bears fruit, I'll post it as a sympy -> numpy evalf.