Opened 12 years ago

Closed 6 years ago

# make heaviside and step play nicely together.

Reported by: Owned by: kcrisman burcin minor sage-8.0 symbolics mforets, rws Ralf Stephan Travis Scrimshaw, Marcelo Forets N/A bbc921e bbc921e655a6a914679f27919d11781f6eb5593c #22838

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

### comment:1 Changed 12 years ago by kcrisman

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

### comment:2 Changed 9 years ago by jdemeyer

Milestone: sage-5.11 → sage-5.12

### comment:3 Changed 9 years ago by vbraun_spam

Milestone: sage-6.1 → sage-6.2

### comment:4 Changed 9 years ago by vbraun_spam

Milestone: sage-6.2 → sage-6.3

### comment:5 Changed 9 years ago by vbraun_spam

Milestone: sage-6.3 → sage-6.4

### comment:7 Changed 6 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:

• `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:

• allowing the possibility to assign a value to heaviside at 0 (as in sympy)

### comment:8 follow-up:  11 Changed 6 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 6 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 6 years ago by mforets

Milestone: sage-6.4 → sage-8.0

### comment:11 in reply to:  8 Changed 6 years ago by mforets

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:  14  15 Changed 6 years ago by tscrim

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 6 years ago by mforets

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:  16 Changed 6 years ago by rws

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:  18 Changed 6 years ago by mforets

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:  20 Changed 6 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 6 years ago by 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 6 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:  21 Changed 6 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:  24 Changed 6 years ago by rws

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 `GinacFunction`s.

### comment:22 follow-up:  23 Changed 6 years ago by mforets

Branch: → u/mforets/make_heaviside_and_step_play_nicely_together_

### comment:23 in reply to:  22 ; follow-up:  27 Changed 6 years ago by mforets

Commit: → d7813f80f995233cb964d8bec9d39b50188d53a1

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:

 ​d7813f8 `call pynac's step from generalized.py`

### comment:24 in reply to:  21 ; follow-up:  30 Changed 6 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 `GinacFunction`s.

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 6 years ago by rws

Branch: u/mforets/make_heaviside_and_step_play_nicely_together_ → u/rws/10070

### comment:26 Changed 6 years ago by git

Commit: d7813f80f995233cb964d8bec9d39b50188d53a1 → 7211abe1d82d928a1b7547727b26e71a683672fc

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

 ​7211abe `10070: use Pynac's step() via GinacFunction`

### comment:27 in reply to:  23 Changed 6 years ago by rws

Branch: u/rws/10070 → u/mforets/make_heaviside_and_step_play_nicely_together_ 7211abe1d82d928a1b7547727b26e71a683672fc → d7813f80f995233cb964d8bec9d39b50188d53a1

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 6 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 6 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 6 years ago by rws

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

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 6 years ago by rws (previous) (diff)

### comment:31 Changed 6 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 6 years ago by egourgoulhon (previous) (diff)

### comment:32 Changed 6 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 6 years ago by rws

Branch: u/mforets/make_heaviside_and_step_play_nicely_together_ → u/rws/10070-1

### comment:34 Changed 6 years ago by rws

Authors: → Ralf Stephan d7813f80f995233cb964d8bec9d39b50188d53a1 → 79ab0ac26bc7dae8f3324b35520e239e1f14f774 → #22838 new → needs_review

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`

### comment:35 follow-up:  36 Changed 6 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 6 years ago by rws

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:  39 Changed 6 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 6 years ago by tscrim (previous) (diff)

### comment:38 Changed 6 years ago by mforets

Branch: u/rws/10070-1 → u/mforets/10070 79ab0ac26bc7dae8f3324b35520e239e1f14f774 → 68e3a60e4fd5ebef56437e4b9740fba993386ed4

update docstring of `step`, add see also blocks.

New commits:

 ​68e3a60 `add seealso blocks`

### comment:39 in reply to:  37 Changed 6 years ago by mforets

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 6 years ago by tscrim

Reviewers: → 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 6 years ago by git

Commit: 68e3a60e4fd5ebef56437e4b9740fba993386ed4 → bbc921e655a6a914679f27919d11781f6eb5593c

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`

### comment:42 Changed 6 years ago by mforets

Reviewers: Travis Scrimshaw → Travis Scrimshaw, Marcelo Forets needs_review → positive_review

@tscrim : merge done, thanks!

### comment:43 Changed 6 years ago by vbraun

Branch: u/mforets/10070 → bbc921e655a6a914679f27919d11781f6eb5593c → fixed positive_review → closed
Note: See TracTickets for help on using tickets.