lambdify DiracDelta and Heaviside as Piecewise?

[G B]

As is probably becoming apparent, I've been busy lambdifying stuff...

My latest runtime problem is trying to convolve things with DiracDeltas.  The result comes back with a Heaviside function.  A quick look through the numpy docs doesn't show a numpy equivalent function for either DiracDelta or Heaviside-- does it makes sense to print these as Piecewise functions?

[G B]

Now having tried it, it looks like the current implementation of Piecewise using "a if test else b" expressions doesn't work with numpy.  The test clause returns an array result, which doesn't like being tested as a unit.

I've implemented _print_Heaviside like this:

    def _print_Heaviside(self,expr):

        a=expr.args

        retVal= "((%s)*(%s) + (%s)*(%s) + (%s)*(%s))" % (self._print(S(0)), self._print(Lt(a[0],0)),

                                                         self._print(div(1,2)), self._print(Eq(a[0],0)),

                                                         self._print(S(1)), self._print(Gt(a[0],0)))

        return retVal

I suspect that those self._print(S(constant)) calls are overkill, but this design seems to work.  I'd recommend recrafting Piecewise in a similar way-- though it might be tricky given that Piecewise as documented allows overlapping conditions.

I'm having a different problem with _print_DiracDelta, which I tried this way:

    def _print_DiracDelta(self,expr):

        a=expr.args

        retVal="((%s)*(%s) + (%s)*(%s))" % (self._print(oo), self._print(Eq(a[0],0)),

                                            self._print(S(0)), self._print(Ne(a[0],0)))

        return retVal

The problem being that oo*0 returns nan which, in turn, seizes the rest of my code by the brainstem.  Anybody have a clever idea of how to get around this?  

I could do some sort of clever indexing to solve the problem in a numpy way, I think, but I'd like to find a general solution.

Thanks-- 

 Greg

[Aaron Meurer]

On Mar 4, 2013, at 4:02 PM, G B <g.c.b....@gmail.com> wrote:

Now having tried it, it looks like the current implementation of Piecewise using "a if test else b" expressions doesn't work with numpy.  The test clause returns an array result, which doesn't like being tested as a unit.

I've implemented _print_Heaviside like this:

    def _print_Heaviside(self,expr):

        a=expr.args

        retVal= "((%s)*(%s) + (%s)*(%s) + (%s)*(%s))" % (self._print(S(0)), self._print(Lt(a[0],0)),

                                                         self._print(div(1,2)), self._print(Eq(a[0],0)),

                                                         self._print(S(1)), self._print(Gt(a[0],0)))

        return retVal

I suspect that those self._print(S(constant)) calls are overkill, but this design seems to work.  I'd recommend recrafting Piecewise in a similar way-- though it might be tricky given that Piecewise as documented allows overlapping conditions.

Ideally you could just call .rewrite(Piecewise). But it isn't implemented (issue 3478). It's easy to fix, though. 

I'm having a different problem with _print_DiracDelta, which I tried this way:

    def _print_DiracDelta(self,expr):

        a=expr.args

        retVal="((%s)*(%s) + (%s)*(%s))" % (self._print(oo), self._print(Eq(a[0],0)),

                                            self._print(S(0)), self._print(Ne(a[0],0)))

        return retVal

The problem being that oo*0 returns nan which, in turn, seizes the rest of my code by the brainstem.  Anybody have a clever idea of how to get around this?  

How exactly are you using it? DiracDelta is not really a function. It's a formalism that lets you use integration to get the value of a function at a point. That makes it inherently symbolic. I don't expect there will be an easy numeric version. You'll need to do some preprocessing to get what you really want later on when you (assumedly) integrate. SymPy can likely help you here, as it can symbolically manipulate Dirac deltas.

Aaron Meurer

I could do some sort of clever indexing to solve the problem in a numpy way, I think, but I'd like to find a general solution.

Thanks-- 

 Greg

On Monday, March 4, 2013 1:35:11 AM UTC-8, G B wrote:

As is probably becoming apparent, I've been busy lambdifying stuff...

My latest runtime problem is trying to convolve things with DiracDeltas.  The result comes back with a Heaviside function.  A quick look through the numpy docs doesn't show a numpy equivalent function for either DiracDelta or Heaviside-- does it makes sense to print these as Piecewise functions?

--
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.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[G B]

Sorry if this winds up being a duplicate post-- I thought I sent this earlier, but I haven't seen it come through.

How exactly are you using it? DiracDelta is not really a function. It's a formalism that lets you use integration to get the value of a function at a point. That makes it inherently symbolic. I don't expect there will be an easy numeric version. You'll need to do some preprocessing to get what you really want later on when you (assumedly) integrate. SymPy can likely help you here, as it can symbolically manipulate Dirac deltas.

I'm using the DiracDelta as one half of a convolution, so you're right, I don't really need a numeric representation.  Getting Heaviside to work is enough.  The likelihood of sampling the interesting part of the delta is numerically zero anyway.  As long as I was doing Heaviside though, DiracDelta seemed like an easy throw-in since it should have followed a similar pattern, but with inf being only a half step from nan, it's treacherous territory.

I'd put the priority on fixing Piecewise, I think.

Cheers--

 Greg

[Aaron Meurer]

On Mar 5, 2013, at 1:37 AM, G B <g.c.b....@gmail.com> wrote:

Sorry if this winds up being a duplicate post-- I thought I sent this earlier, but I haven't seen it come through.

How exactly are you using it? DiracDelta is not really a function. It's a formalism that lets you use integration to get the value of a function at a point. That makes it inherently symbolic. I don't expect there will be an easy numeric version. You'll need to do some preprocessing to get what you really want later on when you (assumedly) integrate. SymPy can likely help you here, as it can symbolically manipulate Dirac deltas.

I'm using the DiracDelta as one half of a convolution, so you're right, I don't really need a numeric representation.  Getting Heaviside to work is enough.  The likelihood of sampling the interesting part of the delta is numerically zero anyway.  As long as I was doing Heaviside though, DiracDelta seemed like an easy throw-in since it should have followed a similar pattern, but with inf being only a half step from nan, it's treacherous territory.

It's not so much inf, but the fact that delta(0) isn't *really* infinity at 0 (because it isn't really even a function). 

For this, you should just symbolically reduce the integral when you get a delta. SymPy should be able to handle this just fine. 

Aaron Meurer 

[Kelly Wilson]

I know this is an old post, but I found a way around it for the heaviside, at least.

lambdify allows you to specify different libraries to use (such as 'numpy') to look for functions, but you also give a dictionary specifying your own.  Just tried the following and it works:

Heavisidenew=lambda x: 1*(x>=0)

h=lambdify(x,sympy.special.delta_functions.Heaviside(x),({"Heaviside":Heavisidenew},'numpy')

and as a test:

h(-1)

h(1)

h(np.linspace(-2,2))

and no problems came up.