Opened 9 years ago

Closed 3 years ago

#10070 closed defect (fixed)

make heaviside and step play nicely together.

Reported by: kcrisman Owned by: burcin
Priority: minor Milestone: sage-8.0
Component: symbolics Keywords:
Cc: mforets, rws Merged in:
Authors: Ralf Stephan Reviewers: Travis Scrimshaw, Marcelo Forets
Report Upstream: N/A Work issues:
Branch: bbc921e (Commits) Commit: bbc921e655a6a914679f27919d11781f6eb5593c
Dependencies: #22838 Stopgaps:

Description

The heaviside built-in function and the step symbolic method are not the same at zero.

sage: heaviside(0)
heaviside(0)
sage: SR(0).step()
1/2

And the documentation seems to indicate that heaviside should be undefined (?) at zero, though it's not definitive.

In addition to reconciling these, we should probably unify notation or something.

Change History (43)

comment:1 Changed 9 years ago by kcrisman

See also #10071 for how this would allow step(hold=True) to be evaluated post hoc.

comment:2 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:3 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:4 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:5 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:6 Changed 3 years ago by mforets

  • Cc mforets added

comment:7 Changed 3 years ago by mforets

the problem is that in many contexts (that is, university courses and books), the heaviside step function is not assigned the value 1/2 at 0, as was the choice in SR(0).step(). if seen as a distribution, there is no need for such choice indeed, since this would be a difference by a set of measure 0.

let me mention in passing that a proper class of distributions with some flesh in it would be a (cool!) feature.

compare to:

  • unit_step has a similar behaviour to Mathematica's UnitStep function.

in sum, this is to suggest:

  • writing "See also" blocks, where it is missing
  • allowing the possibility to assign a value to heaviside at 0 (as in sympy)

comment:8 follow-up: Changed 3 years ago by tscrim

From the wikipedia page, it is mostly a convention choice and depends on the situation. The current behavior actually doesn't allow it to be evaluated at 0:

sage: heaviside(0)
heaviside(0)
sage: heaviside(0).n()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
...
TypeError: function returned None, expected return value of type sage.symbolic.expression.Expression

There seem to be some other differences with step:

sage: heaviside(I)
heaviside(I)
sage: SR(I).step()
1
sage: plot(x.step())
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
...
RuntimeError: cannot find SFunction in table

So both of these approaches might need some work.

comment:9 Changed 3 years ago by mforets

i think this ticket should also fix integration with the default interface (Maxima):

sage: integrate(heaviside(x), x, 0, 10)    # unevaluated?
integrate(heaviside(x), x, 0, 10)
sage: integrate(heaviside(x), x, 0, 10, algorithm='sympy')    # ok
10

comment:10 Changed 3 years ago by mforets

  • Milestone changed from sage-6.4 to sage-8.0

comment:11 in reply to: ↑ 8 Changed 3 years ago by mforets

Replying to tscrim:

From the wikipedia page, it is mostly a convention choice and depends on the situation. The current behavior actually doesn't allow it to be evaluated at 0:

sage: heaviside(0)
heaviside(0)
sage: heaviside(0).n()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
...
TypeError: function returned None, expected return value of type sage.symbolic.expression.Expression

ok. i think that it should rather break with: TypeError: cannot evaluate symbolic expression numerically (the problem is that currently it returns None). also in the help it should be clear that the value at 0 can be specified ad-hoc via an additional argument (as in sympy).

There seem to be some other differences with step:

sage: heaviside(I)
heaviside(I)
sage: SR(I).step()
1
sage: plot(x.step())
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
...
RuntimeError: cannot find SFunction in table

intereseting, i didn't yet understand how to "add new functions to the SFunction table", hints welcome. for instance, plot(x.abs()) does work, as well as a couple of other Pynac functions at /libs/pynac/pynac.pyd. however, in that file there exist all g_step, step_serial and py_step entries.

comment:12 follow-ups: Changed 3 years ago by tscrim

  • Cc rws added

Ralf, can you help us with the SFunction table (comment:11)? I haven't ever used it either.

I would also appreciate any thoughts you had about the Heaviside step function. On a sight tangent, I've been told it is somewhat slow for use in SageManifolds and if there was a way to make it faster (I can cc Eric too for more details).

comment:14 in reply to: ↑ 12 Changed 3 years ago by mforets

Replying to tscrim:

I would also appreciate any thoughts you had about the Heaviside step function. On a sight tangent, I've been told it is somewhat slow for use in SageManifolds and if there was a way to make it faster (I can cc Eric too for more details).

the module piecewise.py is using from sage.ext.fast_callable import ExpressionTreeBuilder, maybe this helps?

heaviside is just heaviside_piecewise = piecewise([((-oo, 0), 0), ((0, oo), 1)], var=x), if 0 is not in the domain

comment:15 in reply to: ↑ 12 ; follow-up: Changed 3 years ago by rws

Replying to tscrim:

Ralf, can you help us with the SFunction table (comment:11)? I haven't ever used it either.

Ok I have now looked a bit into this. The SFunctionTable problem comes from the fact that there is no step function class in sage/functions and you get the raw Pynac function object that is not registered with the (Python) table. The same BTW applies to x.csgn(). Since there is no Function object step(x) is also not available. So if you want a full feature step function different from heaviside and unit_step the best is to create such a GinacFunction in eg functions/generalized.py. Else I would argue to replace the code in Expression.step() with a call to heaviside or unit_step, and in the longer term, deprecate the member function in favor of an alias step = in functions/generalized.py.

I would also appreciate any thoughts you had about the Heaviside step function. On a sight tangent, I've been told it is somewhat slow for use in SageManifolds and if there was a way to make it faster (I can cc Eric too for more details).

Creation is not really slow but numerics is:

sage: %timeit _=heaviside(x)
The slowest run took 5.08 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 13.4 µs per loop
sage: %timeit _=heaviside(2.7)
The slowest run took 4.33 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 58.9 µs per loop
sage: %timeit _=heaviside(ZZ(8))
10000 loops, best of 3: 78.8 µs per loop
sage: %timeit _=heaviside(QQ(8))
10000 loops, best of 3: 82.5 µs per loop

The ZZ and QQ numbers can be easily and well improved by changing FunctionHeaviside to a GinacFunction and checking for these types in C++ because the GMP functions are available there and Python is not needed. The evalf part still needs Python. I guess this all applies to unit_step.

So, is a step() separate from heaviside and unit_step really needed? If not which is it?

comment:16 in reply to: ↑ 15 ; follow-up: Changed 3 years ago by mforets

Replying to rws:

The ZZ and QQ numbers can be easily and well improved by changing FunctionHeaviside to a GinacFunction and checking for these types in C++ [...]

thanks! so i tried BuiltinFunction? -> GinacFunction? in the def of FunctionHeaviside? (and at the init method), and sage crashes when we do ./sage -br (compiles but crashes at runtime). why?

then i defined a new class FunctionHeavisideGinac(GinacFunction): just calling it heaviside_ginac, and then sage doesn't crash and we can use it. i don't have a clue what may be the cause..

comment:17 follow-up: Changed 3 years ago by tscrim

I bet the slow part of _evalf_ comes from this line:

approx_x = ComplexIntervalField()(x)

which has a Python (i.e., not Cython) __call__ function (contrast this with Parent.__call__).

So, is a step() separate from heaviside and unit_step really needed? If not which is it?

This one is hard because they all have slightly different behaviors with dealing with 0. The natural inclination from me is to combine them all into 1 function, where we can give it an optional parameter for what to do when evaluating at 0. However, this has problems when converting to other programs and keeping with their conventions; e.g., mathematica for heaviside and unit_step are different functions. I also don't have good use cases, so I can't give definitive answers on this.

comment:18 in reply to: ↑ 16 Changed 3 years ago by rws

Replying to mforets:

Replying to rws:

The ZZ and QQ numbers can be easily and well improved by changing FunctionHeaviside to a GinacFunction and checking for these types in C++ [...]

thanks! so i tried BuiltinFunction? -> GinacFunction? in the def of FunctionHeaviside? (and at the init method), and sage crashes when we do ./sage -br (compiles but crashes at runtime). why?

Probably because you still had "heaviside" as name argument in the GinacFunction.__init__ call but I didn't check (note that Pynac matches Pynac functions with Sage functions by name and the name of the only Pynac function of this kind is step). Rather, what I meant was then write a second Pynac function corresponding to heaviside.

then i defined a new class FunctionHeavisideGinac(GinacFunction): just calling it heaviside_ginac, and then sage doesn't crash and we can use it. i don't have a clue what may be the cause..

But then nothing is done in C++ so no speedup.

comment:19 Changed 3 years ago by rws

Note that if you allow a separate step() then a Function in src/sage/functions would be needed because atm only x.step() works not step(x).

comment:20 in reply to: ↑ 17 ; follow-up: Changed 3 years ago by mforets

i'd like to argue in favor of 1 heaviside 1 unit_step (not merged), that is to keep the current implementation, and to follow the recommendation of rws about pointing SR(x).step() into unit_step() (*)

for some people, the fact that heaviside is undefined at 0 is good and expected, because what it really defines is an equivalence class, and it is typically used in the context of distributions. in concrete, linear systems theory (to model a on/off switch) and for differential -> algebraic transformations in physics problems (applications of laplace transform).

but for other set of users (as in the OP description), this may be confusing, and in many contexts one wouldn't like to talk about equivalence class and the dirac_delta (more on this below).

if i had to change something it would be the derivative operator of unit_step: again, one wouldn't like to talk about dirac_delta in many contexts because that is a more advanced notion. instead i would give the answer piecewise([((-oo, 0), 0), ((0, oo), 0)]), which is nice because 0 is not in the domain when you ask f(0).. we can also add a message See heaviside for a generalized notion.

(*) argument: it's not possible to make a plot, plot(SR(x).step()), for such a long time, that i suspect it's not widely used. it should be merged with unit_step() or heaviside() with value 1/2 at 0. don't have any particular preference, probably the latter.

(note: ..this comment could be more illustrative of different uses (or not) with proper references to books/courses but i don't have time to go into that now :/ )

comment:21 in reply to: ↑ 20 ; follow-up: Changed 3 years ago by rws

Replying to mforets:

i'd like to argue in favor of 1 heaviside 1 unit_step (not merged), that is to keep the current implementation, and to follow the recommendation of rws about pointing SR(x).step() into unit_step() (*)

step() timings:

sage: %timeit _=SR(1).step()
The slowest run took 11.30 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.77 µs per loop
sage: %timeit _=SR(1/2).step()
The slowest run took 7.51 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.92 µs per loop
sage: %timeit _=SR(1.5).step()
The slowest run took 6.56 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 8.21 µs per loop

So the Pynac's step can be used to speed up heaviside and unit_step. The procedure would be to rename step in Pynac to unit_step (adapting the behaviour at 0), and copy the code for heaviside in Pynac (again adapting the behaviour at 0), and after that make the Sage functions GinacFunctions.

comment:22 follow-up: Changed 3 years ago by mforets

  • Branch set to u/mforets/make_heaviside_and_step_play_nicely_together_

comment:23 in reply to: ↑ 22 ; follow-up: Changed 3 years ago by mforets

  • Commit set to d7813f80f995233cb964d8bec9d39b50188d53a1

Replying to mforets:

presuming that this commit is not correct, but is a starter. still i need some time to understand how things work.

about the timing, already it gets down to 22us with this little change. i suppose that calling symbolic expressions step is not good (?), but instead a call to new_Expression_from_GEx requires a cimport?..


New commits:

d7813f8call pynac's step from generalized.py

comment:24 in reply to: ↑ 21 ; follow-up: Changed 3 years ago by mforets

Replying to rws: and copy the code for heaviside in Pynac (again adapting the behaviour at 0), and after that make the Sage functions GinacFunctions.

Ok so in the end of the day, in generalized.py we shall be left with:

  • just one FunctionUnitStep(GinacFunction), or with
  • one FunctionUnitStep(BuiltinFunction) plus one FunctionUnitStep(GinacFunction)?

comment:25 Changed 3 years ago by rws

  • Branch changed from u/mforets/make_heaviside_and_step_play_nicely_together_ to u/rws/10070

comment:26 Changed 3 years ago by git

  • Commit changed from d7813f80f995233cb964d8bec9d39b50188d53a1 to 7211abe1d82d928a1b7547727b26e71a683672fc

Branch pushed to git repo; I updated commit sha1. New commits:

7211abe10070: use Pynac's step() via GinacFunction

comment:27 in reply to: ↑ 23 Changed 3 years ago by rws

  • Branch changed from u/rws/10070 to u/mforets/make_heaviside_and_step_play_nicely_together_
  • Commit changed from 7211abe1d82d928a1b7547727b26e71a683672fc to d7813f80f995233cb964d8bec9d39b50188d53a1

Replying to mforets:

Replying to mforets:

about the timing, already it gets down to 22us with this little change. i suppose that calling symbolic expressions step is not good (?), but instead a call to new_Expression_from_GEx requires a cimport?..

Yes but I think that's not the reason that it still needs 22us. The branch I uploaded leaves only the derivative function and uses Pynac's step for everything else. It gets down to 6us.

presuming that this commit is not correct, but is a starter. still i need some time to understand how things work.

comment:28 Changed 3 years ago by rws

Of course in the final version the name is set back to unit_step after step in Pynac is named unit_step. Also the doctests need to be moved to the class docstring.

comment:29 Changed 3 years ago by rws

I will try to include the relevant changes in the next Pynac version 0.7.6.

comment:30 in reply to: ↑ 24 Changed 3 years ago by rws

Replying to mforets:

Replying to rws: and copy the code for heaviside in Pynac (again adapting the behaviour at 0), and after that make the Sage functions GinacFunctions.

Ok so in the end of the day, in generalized.py we shall be left with:

  • just one FunctionUnitStep(GinacFunction), or with
  • one FunctionUnitStep(BuiltinFunction) plus one FunctionUnitStep(GinacFunction)?

One FunctionHeaviside(GinacFunction) plus one FunctionUnitStep(GinacFunction).

(EDITED)

Last edited 3 years ago by rws (previous) (diff)

comment:31 Changed 3 years ago by egourgoulhon

I don't know if this pertains to this ticket but, aside from the value at 0, another difference between heaviside and unit_step is that the former cannot be used in numerical resolutions:

sage: y = var('y')
sage: desolve_rk4(x,y,ics=[0,0],end_points=1,step=0.5) # solution is x^2/2
[[0, 0], [0.5, 0.125], [1.0, 0.4999999999999999]]
sage: desolve_rk4(x*unit_step(1+x),y,ics=[0,0],end_points=1,step=0.5) # OK
[[0, 0], [0.5, 0.125], [1.0, 0.4999999999999999]]
sage: desolve_rk4(x*heaviside(1+x),y,ics=[0,0],end_points=1,step=0.5) # ??
[[0, 0]]

The value at 0 cannot explain such difference because x varies in [0,1], so that 1+x>0 and heaviside(1+x) should always return 1.

Last edited 3 years ago by egourgoulhon (previous) (diff)

comment:32 Changed 3 years ago by tscrim

Note that

desolve_rk4(x*(1+x).step(),y,ics=[0,0],end_points=1,step=0.5)

results in an error, but that should be fixed by the proposed change.

comment:33 Changed 3 years ago by rws

  • Branch changed from u/mforets/make_heaviside_and_step_play_nicely_together_ to u/rws/10070-1

comment:34 Changed 3 years ago by rws

  • Authors set to Ralf Stephan
  • Commit changed from d7813f80f995233cb964d8bec9d39b50188d53a1 to 79ab0ac26bc7dae8f3324b35520e239e1f14f774
  • Dependencies set to #22838
  • Status changed from new to needs_review

New commits:

3c848da22838: pkg version/checksum
d89e6d022838: remove unnecessary patch
4a3011922838: interface change
7d3fd6022838: doctest adaptation to ticket 10070
79ab0ac10070: heaviside, unit_step revamp

comment:35 follow-up: Changed 3 years ago by mforets

@rws : cool, i was able to test it! (thanks to your answer in the other ticket). Wow! <surprise face> so heaviside is really defined in inifcns.cpp of the Ginac project.

we still have issues that have been raised in this thread that should be fixed (optional definition at 0, use in integrate fails, use in desolve_rk4 fails) -- but i guess that should be a new ticket for heaviside only?

comment:36 in reply to: ↑ 35 Changed 3 years ago by rws

Replying to mforets:

we still have issues that have been raised in this thread that should be fixed (optional definition at 0, use in integrate fails, use in desolve_rk4 fails) -- but i guess that should be a new ticket for heaviside only?

Agreed.

comment:37 follow-up: Changed 3 years ago by tscrim

There was a part of me that hoped this would fix comment:31 (it does fix comment:32 as I predicted), but it doesn't. However, I agree that this should be a separate ticket.

I think we need to also update the documentation for step to say it returns the unit_step, which has a different convention from heaviside. This should fix some of the issues with this ticket; the other would be to add .. SEEALSO:: links between the two with a mild note on the difference. I do not agree that they should have different derivatives as the Dirac delta behaves essentially in the same way.

Last edited 3 years ago by tscrim (previous) (diff)

comment:38 Changed 3 years ago by mforets

  • Branch changed from u/rws/10070-1 to u/mforets/10070
  • Commit changed from 79ab0ac26bc7dae8f3324b35520e239e1f14f774 to 68e3a60e4fd5ebef56437e4b9740fba993386ed4

update docstring of step, add see also blocks.


New commits:

68e3a60add seealso blocks

comment:39 in reply to: ↑ 37 Changed 3 years ago by mforets

Replying to tscrim:

There was a part of me that hoped this would fix comment:31 (it does fix comment:32 as I predicted), but it doesn't. However, I agree that this should be a separate ticket.

separate ticket: #22849, Heaviside in numerical resolutions.

comment:40 Changed 3 years ago by tscrim

  • Reviewers set to Travis Scrimshaw

Your last commit will likely conflict with a similar change I made on #22838. Once you merge in those changes, you can set a positive review on my behalf (and add your name to the reviewers list).

comment:41 Changed 3 years ago by git

  • Commit changed from 68e3a60e4fd5ebef56437e4b9740fba993386ed4 to bbc921e655a6a914679f27919d11781f6eb5593c

Branch pushed to git repo; I updated commit sha1. New commits:

ca852a1Some reviewer tweaks.
bbc921eMerge branch 'u/tscrim/pynac_0_7_6-22838' of git://trac.sagemath.org/sage into u/mforets/10070

comment:42 Changed 3 years ago by mforets

  • Reviewers changed from Travis Scrimshaw to Travis Scrimshaw, Marcelo Forets
  • Status changed from needs_review to positive_review

@tscrim : merge done, thanks!

comment:43 Changed 3 years ago by vbraun

  • Branch changed from u/mforets/10070 to bbc921e655a6a914679f27919d11781f6eb5593c
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.