make heaviside and step play nicely together.

[kcrisman]

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.

Depends on #22838

CC: @mforets @rwst

Component: symbolics

Author: Ralf Stephan

Branch/Commit: bbc921e

Reviewer: Travis Scrimshaw, Marcelo Forets

Issue created by migration from https://trac.sagemath.org/ticket/10070

[kcrisman]

comment:1

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

[mforets]

comment:7

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:

• heaviside has a similar behaviour to Mathematica's HeavisideTheta function.

• 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)

[tscrim]

comment:8

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.

[mforets]

comment:9

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

[mforets]

comment:11

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.

[tscrim]

comment:12

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).

[rwst]

comment:13

I opened an issue: pynac/pynac#243

[mforets]

comment:14

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

[rwst]

comment:15

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?

[mforets]

comment:16

Replying to @rwst:

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..

[tscrim]

comment:17

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.

[rwst]

comment:18

Replying to @mforets:

Replying to @rwst:

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.

[rwst]

comment:19

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).

[mforets]

comment:20

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 :/ )

[rwst]

comment:21

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.

[rwst]

comment:29

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

[rwst]

comment:30

Replying to @mforets:

Replying to @rwst:
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)

[egourgoulhon]

comment:31

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.

[tscrim]

comment:32

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.

[rwst]

Changed branch from u/mforets/make_heaviside_and_step_play_nicely_together_ to u/rws/10070-1

[rwst]

Dependencies: #22838

[rwst]

New commits:

3c848da	22838: pkg version/checksum
d89e6d0	22838: remove unnecessary patch
4a30119	22838: interface change
7d3fd60	22838: doctest adaptation to ticket 10070
79ab0ac	10070: heaviside, unit_step revamp

[rwst]

Author: Ralf Stephan

[rwst]

Changed commit from d7813f8 to 79ab0ac

[mforets]

comment:35

@rwst : cool, i was able to test it! (thanks to your answer in the other ticket). Wow! 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?

[rwst]

comment:36

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.

[tscrim]

comment:37

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.

[mforets]

Changed commit from 79ab0ac to 68e3a60

[mforets]

Changed branch from u/rws/10070-1 to u/mforets/10070

[mforets]

comment:38

update docstring of step, add see also blocks.

---

New commits:

68e3a60

add seealso blocks

[mforets]

comment:39

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.

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:40

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).

[sagetrac-git]

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

ca852a1	Some reviewer tweaks.
bbc921e	Merge branch 'u/tscrim/pynac_0_7_6-22838' of git://trac.sagemath.org/sage into u/mforets/10070

[sagetrac-git]

Changed commit from 68e3a60 to bbc921e

[mforets]

comment:42

@tscrim : merge done, thanks!

[mforets]

Changed reviewer from Travis Scrimshaw to Travis Scrimshaw, Marcelo Forets

[vbraun]

Changed branch from u/mforets/10070 to bbc921e