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:  sage8.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 builtin 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
comment:2 Changed 6 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:3 Changed 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:4 Changed 6 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:5 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:6 Changed 3 years ago by
 Cc mforets added
comment:7 Changed 3 years ago by
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)
comment:8 followup: ↓ 11 Changed 3 years ago by
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
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
 Milestone changed from sage6.4 to sage8.0
comment:11 in reply to: ↑ 8 Changed 3 years ago by
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 adhoc 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 followups: ↓ 14 ↓ 15 Changed 3 years ago by
 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:13 Changed 3 years ago by
I opened an issue: https://github.com/pynac/pynac/issues/243
comment:14 in reply to: ↑ 12 Changed 3 years ago by
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 ; followup: ↓ 16 Changed 3 years ago by
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 ; followup: ↓ 18 Changed 3 years ago by
Replying to rws:
The ZZ and QQ numbers can be easily and well improved by changing
FunctionHeaviside
to aGinacFunction
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 followup: ↓ 20 Changed 3 years ago by
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 fromheaviside
andunit_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
Replying to mforets:
Replying to rws:
The ZZ and QQ numbers can be easily and well improved by changing
FunctionHeaviside
to aGinacFunction
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 itheaviside_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
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 ; followup: ↓ 21 Changed 3 years ago by
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 ; followup: ↓ 24 Changed 3 years ago by
Replying to mforets:
i'd like to argue in favor of 1
heaviside
1unit_step
(not merged), that is to keep the current implementation, and to follow the recommendation of rws about pointingSR(x).step()
intounit_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 followup: ↓ 23 Changed 3 years ago by
 Branch set to u/mforets/make_heaviside_and_step_play_nicely_together_
comment:23 in reply to: ↑ 22 ; followup: ↓ 27 Changed 3 years ago by
 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:
d7813f8  call pynac's step from generalized.py

comment:24 in reply to: ↑ 21 ; followup: ↓ 30 Changed 3 years ago by
Replying to rws: and copy the code for
heaviside
in Pynac (again adapting the behaviour at 0), and after that make the Sage functionsGinacFunction
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 oneFunctionUnitStep(GinacFunction)
?
comment:25 Changed 3 years ago by
 Branch changed from u/mforets/make_heaviside_and_step_play_nicely_together_ to u/rws/10070
comment:26 Changed 3 years ago by
 Commit changed from d7813f80f995233cb964d8bec9d39b50188d53a1 to 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 3 years ago by
 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 tonew_Expression_from_GEx
requires acimport
?..
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
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
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
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 functionsGinacFunction
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 oneFunctionUnitStep(GinacFunction)
?
One FunctionHeaviside(GinacFunction)
plus one FunctionUnitStep(GinacFunction)
.
(EDITED)
comment:31 Changed 3 years ago by
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
.
comment:32 Changed 3 years ago by
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
 Branch changed from u/mforets/make_heaviside_and_step_play_nicely_together_ to u/rws/100701
comment:34 Changed 3 years ago by
 Commit changed from d7813f80f995233cb964d8bec9d39b50188d53a1 to 79ab0ac26bc7dae8f3324b35520e239e1f14f774
 Dependencies set to #22838
 Status changed from new to needs_review
comment:35 followup: ↓ 36 Changed 3 years ago by
@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
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 indesolve_rk4
fails)  but i guess that should be a new ticket forheaviside
only?
Agreed.
comment:37 followup: ↓ 39 Changed 3 years ago by
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.
comment:38 Changed 3 years ago by
 Branch changed from u/rws/100701 to u/mforets/10070
 Commit changed from 79ab0ac26bc7dae8f3324b35520e239e1f14f774 to 68e3a60e4fd5ebef56437e4b9740fba993386ed4
comment:39 in reply to: ↑ 37 Changed 3 years ago by
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
 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
 Commit changed from 68e3a60e4fd5ebef56437e4b9740fba993386ed4 to bbc921e655a6a914679f27919d11781f6eb5593c
comment:42 Changed 3 years ago by
 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
 Branch changed from u/mforets/10070 to bbc921e655a6a914679f27919d11781f6eb5593c
 Resolution set to fixed
 Status changed from positive_review to closed
See also #10071 for how this would allow
step(hold=True)
to be evaluated post hoc.