Repository: sympy/sympy

Issue:
get_points() generates error when used with Heaviside() function
The following code used to work with older versions of sympy (e.g. 1.3) but generates an error with 1.7.1

```
import sympy as sym

t = sym.symbols('t')
p = sym.plot(sym.Heaviside(t), (t, -5, 5), show=False)
points = p[0].get_points()
```

results in

`AttributeError: 'ImmutableDenseNDimArray' object has no attribute 'as_coefficient'`

Discussion:

Patch:
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index 1c3b4e50c4..66a96335eb 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -442,7 +442,7 @@ def _eval_rewrite_as_Piecewise(self, arg, **kwargs):
     def _eval_rewrite_as_Heaviside(self, arg, **kwargs):
         from sympy.functions.special.delta_functions import Heaviside
         if arg.is_extended_real:
-            return Heaviside(arg, H0=S(1)/2) * 2 - 1
+            return Heaviside(arg) * 2 - 1
 
     def _eval_rewrite_as_Abs(self, arg, **kwargs):
         return Piecewise((0, Eq(arg, 0)), (arg / Abs(arg), True))
diff --git a/sympy/functions/special/delta_functions.py b/sympy/functions/special/delta_functions.py
index e5315e181d..e3ffc481f2 100644
--- a/sympy/functions/special/delta_functions.py
+++ b/sympy/functions/special/delta_functions.py
@@ -405,43 +405,44 @@ def _sage_(self):
 
 class Heaviside(Function):
     r"""
-    Heaviside Piecewise function.
+    Heaviside step function.
 
     Explanation
     ===========
 
-    Heaviside function has the following properties:
+    The Heaviside step function has the following properties:
 
     1) $\frac{d}{d x} \theta(x) = \delta(x)$
-    2) $\theta(x) = \begin{cases} 0 & \text{for}\: x < 0 \\ \text{undefined} &
+    2) $\theta(x) = \begin{cases} 0 & \text{for}\: x < 0 \\ \frac{1}{2} &
        \text{for}\: x = 0 \\1 & \text{for}\: x > 0 \end{cases}$
     3) $\frac{d}{d x} \max(x, 0) = \theta(x)$
 
     Heaviside(x) is printed as $\theta(x)$ with the SymPy LaTeX printer.
 
-    Regarding to the value at 0, Mathematica defines $\theta(0)=1$, but Maple
-    uses $\theta(0) = \text{undefined}$. Different application areas may have
-    specific conventions. For example, in control theory, it is common practice
-    to assume $\theta(0) = 0$ to match the Laplace transform of a DiracDelta
-    distribution.
+    The value at 0 is set differently in different fields. SymPy uses 1/2,
+    which is a convention from electronics and signal processing, and is
+    consistent with solving improper integrals by Fourier transform and
+    convolution.
 
-    To specify the value of Heaviside at ``x=0``, a second argument can be
-    given. Omit this 2nd argument or pass ``None`` to recover the default
-    behavior.
+    To specify a different value of Heaviside at ``x=0``, a second argument
+    can be given. Using ``Heaviside(x, nan)`` gives an expression that will
+    evaluate to nan for x=0.
+
+    .. versionchanged:: 1.9 ``Heaviside(0)`` now returns 1/2 (before: undefined)
 
     Examples
     ========
 
-    >>> from sympy import Heaviside, S
+    >>> from sympy import Heaviside, nan
     >>> from sympy.abc import x
     >>> Heaviside(9)
     1
     >>> Heaviside(-9)
     0
     >>> Heaviside(0)
-    Heaviside(0)
-    >>> Heaviside(0, S.Half)
     1/2
+    >>> Heaviside(0, nan)
+    nan
     >>> (Heaviside(x) + 1).replace(Heaviside(x), Heaviside(x, 1))
     Heaviside(x, 1) + 1
 
@@ -487,21 +488,17 @@ def fdiff(self, argindex=1):
 
         """
         if argindex == 1:
-            # property number 1
             return DiracDelta(self.args[0])
         else:
             raise ArgumentIndexError(self, argindex)
 
-    def __new__(cls, arg, H0=None, **options):
+    def __new__(cls, arg, H0=S.Half, **options):
         if isinstance(H0, Heaviside) and len(H0.args) == 1:
-            H0 = None
-
-        if H0 is None:
-            return super(cls, cls).__new__(cls, arg, **options)
+            H0 = S.Half
         return super(cls, cls).__new__(cls, arg, H0, **options)
 
     @classmethod
-    def eval(cls, arg, H0=None):
+    def eval(cls, arg, H0=S.Half):
         """
         Returns a simplified form or a value of Heaviside depending on the
         argument passed by the Heaviside object.
@@ -522,13 +519,13 @@ class is about to be instantiated and it returns either some simplified
         >>> from sympy.abc import x
 
         >>> Heaviside(x)
-        Heaviside(x)
+        Heaviside(x, 1/2)
 
         >>> Heaviside(19)
         1
 
         >>> Heaviside(0)
-        Heaviside(0)
+        1/2
 
         >>> Heaviside(0, 1)
         1
@@ -539,7 +536,7 @@ class is about to be instantiated and it returns either some simplified
         >>> Heaviside(S.NaN)
         nan
 
-        >>> Heaviside(x).eval(100)
+        >>> Heaviside(x).eval(42)
         1
 
         >>> Heaviside(x - 100).subs(x, 5)
@@ -551,8 +548,9 @@ class is about to be instantiated and it returns either some simplified
         Parameters
         ==========
 
-        arg : argument passed by HeaviSide object
-        HO : boolean flag for HeaviSide Object is set to True
+        arg : argument passed by Heaviside object
+
+        H0 : value of Heaviside(0)
 
         """
         H0 = sympify(H0)
@@ -575,63 +573,66 @@ def _eval_rewrite_as_Piecewise(self, arg, H0=None, **kwargs):
         Examples
         ========
 
-        >>> from sympy import Heaviside, Piecewise, Symbol
+        >>> from sympy import Heaviside, Piecewise, Symbol, nan
         >>> x = Symbol('x')
 
         >>> Heaviside(x).rewrite(Piecewise)
-        Piecewise((0, x < 0), (Heaviside(0), Eq(x, 0)), (1, x > 0))
+        Piecewise((0, x < 0), (1/2, Eq(x, 0)), (1, x > 0))
+
+        >>> Heaviside(x,nan).rewrite(Piecewise)
+        Piecewise((0, x < 0), (nan, Eq(x, 0)), (1, x > 0))
 
         >>> Heaviside(x - 5).rewrite(Piecewise)
-        Piecewise((0, x - 5 < 0), (Heaviside(0), Eq(x - 5, 0)), (1, x - 5 > 0))
+        Piecewise((0, x - 5 < 0), (1/2, Eq(x - 5, 0)), (1, x - 5 > 0))
 
         >>> Heaviside(x**2 - 1).rewrite(Piecewise)
-        Piecewise((0, x**2 - 1 < 0), (Heaviside(0), Eq(x**2 - 1, 0)), (1, x**2 - 1 > 0))
+        Piecewise((0, x**2 - 1 < 0), (1/2, Eq(x**2 - 1, 0)), (1, x**2 - 1 > 0))
 
         """
-        if H0 is None:
-            return Piecewise((0, arg < 0), (Heaviside(0), Eq(arg, 0)), (1, arg > 0))
         if H0 == 0:
             return Piecewise((0, arg <= 0), (1, arg > 0))
         if H0 == 1:
             return Piecewise((0, arg < 0), (1, arg >= 0))
         return Piecewise((0, arg < 0), (H0, Eq(arg, 0)), (1, arg > 0))
 
-    def _eval_rewrite_as_sign(self, arg, H0=None, **kwargs):
+    def _eval_rewrite_as_sign(self, arg, H0=S.Half, **kwargs):
         """
         Represents the Heaviside function in the form of sign function.
 
         Explanation
         ===========
 
-        The value of the second argument of Heaviside must specify Heaviside(0)
-        = 1/2 for rewritting as sign to be strictly equivalent. For easier
-        usage, we also allow this rewriting when Heaviside(0) is undefined.
+        The value of Heaviside(0) must be 1/2 for rewritting as sign to be
+        strictly equivalent. For easier usage, we also allow this rewriting
+        when Heaviside(0) is undefined.
 
         Examples
         ========
 
-        >>> from sympy import Heaviside, Symbol, sign, S
+        >>> from sympy import Heaviside, Symbol, sign, nan
         >>> x = Symbol('x', real=True)
+        >>> y = Symbol('y')
 
-        >>> Heaviside(x, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x).rewrite(sign)
         sign(x)/2 + 1/2
 
         >>> Heaviside(x, 0).rewrite(sign)
         Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (0, True))
 
-        >>> Heaviside(x - 2, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x, nan).rewrite(sign)
+        Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (nan, True))
+
+        >>> Heaviside(x - 2).rewrite(sign)
         sign(x - 2)/2 + 1/2
 
-        >>> Heaviside(x**2 - 2*x + 1, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x**2 - 2*x + 1).rewrite(sign)
         sign(x**2 - 2*x + 1)/2 + 1/2
 
-        >>> y = Symbol('y')
-
         >>> Heaviside(y).rewrite(sign)
-        Heaviside(y)
+        Heaviside(y, 1/2)
 
         >>> Heaviside(y**2 - 2*y + 1).rewrite(sign)
-        Heaviside(y**2 - 2*y + 1)
+        Heaviside(y**2 - 2*y + 1, 1/2)
 
         See Also
         ========
@@ -648,7 +649,7 @@ def _eval_rewrite_as_sign(self, arg, H0=None, **kwargs):
                 (pw1, True))
             return pw2
 
-    def _eval_rewrite_as_SingularityFunction(self, args, **kwargs):
+    def _eval_rewrite_as_SingularityFunction(self, args, H0=S.Half, **kwargs):
         """
         Returns the Heaviside expression written in the form of Singularity
         Functions.
diff --git a/sympy/functions/special/singularity_functions.py b/sympy/functions/special/singularity_functions.py
index 69695de601..d568c28a1f 100644
--- a/sympy/functions/special/singularity_functions.py
+++ b/sympy/functions/special/singularity_functions.py
@@ -65,11 +65,11 @@ class SingularityFunction(Function):
 
     >>> expr = SingularityFunction(x, 4, 5) + SingularityFunction(x, -3, -1) - SingularityFunction(x, 0, -2)
     >>> expr.rewrite(Heaviside)
-    (x - 4)**5*Heaviside(x - 4) + DiracDelta(x + 3) - DiracDelta(x, 1)
+    (x - 4)**5*Heaviside(x - 4, 1/2) + DiracDelta(x + 3) - DiracDelta(x, 1)
     >>> expr.rewrite(DiracDelta)
-    (x - 4)**5*Heaviside(x - 4) + DiracDelta(x + 3) - DiracDelta(x, 1)
+    (x - 4)**5*Heaviside(x - 4, 1/2) + DiracDelta(x + 3) - DiracDelta(x, 1)
     >>> expr.rewrite('HeavisideDiracDelta')
-    (x - 4)**5*Heaviside(x - 4) + DiracDelta(x + 3) - DiracDelta(x, 1)
+    (x - 4)**5*Heaviside(x - 4, 1/2) + DiracDelta(x + 3) - DiracDelta(x, 1)
 
     See Also
     ========
diff --git a/sympy/integrals/deltafunctions.py b/sympy/integrals/deltafunctions.py
index e0acdab905..dd90867575 100644
--- a/sympy/integrals/deltafunctions.py
+++ b/sympy/integrals/deltafunctions.py
@@ -123,9 +123,9 @@ def deltaintegrate(f, x):
         >>> from sympy.integrals.deltafunctions import deltaintegrate
         >>> from sympy import sin, cos, DiracDelta
         >>> deltaintegrate(x*sin(x)*cos(x)*DiracDelta(x - 1), x)
-        sin(1)*cos(1)*Heaviside(x - 1)
+        sin(1)*cos(1)*Heaviside(x - 1, 1/2)
         >>> deltaintegrate(y**2*DiracDelta(x - z)*DiracDelta(y - z), y)
-        z**2*DiracDelta(x - z)*Heaviside(y - z)
+        z**2*DiracDelta(x - z)*Heaviside(y - z, 1/2)
 
     See Also
     ========
diff --git a/sympy/integrals/meijerint.py b/sympy/integrals/meijerint.py
index 86cecddb33..3acab5cac7 100644
--- a/sympy/integrals/meijerint.py
+++ b/sympy/integrals/meijerint.py
@@ -2087,7 +2087,7 @@ def meijerint_inversion(f, x, t):
     >>> from sympy.abc import x, t
     >>> from sympy.integrals.meijerint import meijerint_inversion
     >>> meijerint_inversion(1/x, x, t)
-    Heaviside(t)
+    Heaviside(t, 1/2)
     """
     from sympy import exp, expand, log, Add, Mul, Heaviside
     f_ = f
diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py
index 739450dd93..f1be3131dd 100644
--- a/sympy/integrals/transforms.py
+++ b/sympy/integrals/transforms.py
@@ -893,11 +893,11 @@ def inverse_mellin_transform(F, s, x, strip, **hints):
 
     >>> f = 1/(s**2 - 1)
     >>> inverse_mellin_transform(f, s, x, (-oo, -1))
-    x*(1 - 1/x**2)*Heaviside(x - 1)/2
+    x*(1 - 1/x**2)*Heaviside(x - 1, 1/2)/2
     >>> inverse_mellin_transform(f, s, x, (-1, 1))
-    -x*Heaviside(1 - x)/2 - Heaviside(x - 1)/(2*x)
+    -x*Heaviside(1 - x, 1/2)/2 - Heaviside(x - 1, 1/2)/(2*x)
     >>> inverse_mellin_transform(f, s, x, (1, oo))
-    (1/2 - x**2/2)*Heaviside(1 - x)/x
+    (1/2 - x**2/2)*Heaviside(1 - x, 1/2)/x
 
     See Also
     ========
@@ -1234,17 +1234,17 @@ def pw_simp(*args):
 
     u = Dummy('u')
 
-    def simp_heaviside(arg):
+    def simp_heaviside(arg, H0=S.Half):
         a = arg.subs(exp(-t), u)
         if a.has(t):
-            return Heaviside(arg)
+            return Heaviside(arg,H0)
         rel = _solve_inequality(a > 0, u)
         if rel.lts == u:
             k = log(rel.gts)
-            return Heaviside(t + k)
+            return Heaviside(t + k,H0)
         else:
             k = log(rel.lts)
-            return Heaviside(-(t + k))
+            return Heaviside(-(t + k),H0)
     f = f.replace(Heaviside, simp_heaviside)
 
     def simp_exp(arg):
@@ -1328,7 +1328,7 @@ def inverse_laplace_transform(F, s, t, plane=None, **hints):
     >>> from sympy.abc import s, t
     >>> a = Symbol('a', positive=True)
     >>> inverse_laplace_transform(exp(-a*s)/s, s, t)
-    Heaviside(-a + t)
+    Heaviside(-a + t, 1/2)
 
     See Also
     ========
diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index 917b703821..06bd234787 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -1573,7 +1573,12 @@ def _print_dirichlet_eta(self, e):
 
     def _print_Heaviside(self, e):
         func_name = greek_unicode['theta'] if self._use_unicode else 'Heaviside'
-        return self._print_Function(e, func_name=func_name)
+        if e.args[1]==1/2:
+            pform = prettyForm(*self._print(e.args[0]).parens())
+            pform = prettyForm(*pform.left(func_name))
+            return pform
+        else:
+            return self._print_Function(e, func_name=func_name)
 
     def _print_fresnels(self, e):
         return self._print_Function(e, func_name="S")


Test patch:
diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py
index 65b068349b..61b63a3198 100644
--- a/sympy/functions/elementary/tests/test_trigonometric.py
+++ b/sympy/functions/elementary/tests/test_trigonometric.py
@@ -1049,7 +1049,7 @@ def test_atan2():
     assert atan2(0, u) == pi
 
     assert atan2(y, oo) ==  0
-    assert atan2(y, -oo)==  2*pi*Heaviside(re(y)) - pi
+    assert atan2(y, -oo)==  2*pi*Heaviside(re(y), S.Half) - pi
 
     assert atan2(y, x).rewrite(log) == -I*log((x + I*y)/sqrt(x**2 + y**2))
     assert atan2(0, 0) is S.NaN
diff --git a/sympy/functions/special/tests/test_delta_functions.py b/sympy/functions/special/tests/test_delta_functions.py
index 72ba93e2c9..69722eac1b 100644
--- a/sympy/functions/special/tests/test_delta_functions.py
+++ b/sympy/functions/special/tests/test_delta_functions.py
@@ -6,6 +6,8 @@
 
 from sympy.testing.pytest import raises, warns_deprecated_sympy
 
+from sympy.core.expr import unchanged
+
 from sympy.core.function import ArgumentIndexError
 
 
@@ -80,19 +82,13 @@ def test_DiracDelta():
 
 
 def test_heaviside():
-    assert Heaviside(0).func == Heaviside
     assert Heaviside(-5) == 0
     assert Heaviside(1) == 1
-    assert Heaviside(nan) is nan
+    assert Heaviside(0) == S.Half
 
     assert Heaviside(0, x) == x
-    assert Heaviside(0, nan) is nan
-    assert Heaviside(x, None) == Heaviside(x)
-    assert Heaviside(0, None) == Heaviside(0)
-
-    # we do not want None and Heaviside(0) in the args:
-    assert Heaviside(x, H0=None).args == (x,)
-    assert Heaviside(x, H0=Heaviside(0)).args == (x,)
+    assert unchanged(Heaviside,x, nan)
+    assert Heaviside(0, nan) == nan
 
     assert adjoint(Heaviside(x)) == Heaviside(x)
     assert adjoint(Heaviside(x - y)) == Heaviside(x - y)
@@ -122,6 +118,8 @@ def test_rewrite():
         Piecewise((0, x <= 0), (1, x > 0)))
     assert Heaviside(x, 1).rewrite(Piecewise) == (
         Piecewise((0, x < 0), (1, x >= 0)))
+    assert Heaviside(x, nan).rewrite(Piecewise) == (
+        Piecewise((0, x < 0), (nan, Eq(x, 0)), (1, x > 0)))
 
     assert Heaviside(x).rewrite(sign) == \
         Heaviside(x, H0=Heaviside(0)).rewrite(sign) == \
@@ -153,13 +151,4 @@ def test_rewrite():
     assert Heaviside(x).rewrite(SingularityFunction) == SingularityFunction(x, 0, 0)
     assert 5*x*y*Heaviside(y + 1).rewrite(SingularityFunction) == 5*x*y*SingularityFunction(y, -1, 0)
     assert ((x - 3)**3*Heaviside(x - 3)).rewrite(SingularityFunction) == (x - 3)**3*SingularityFunction(x, 3, 0)
-    assert Heaviside(0).rewrite(SingularityFunction) == SingularityFunction(0, 0, 0)
-
-def test_issue_15923():
-    x = Symbol('x', real=True)
-    assert Heaviside(x).rewrite(Piecewise, H0=0) == (
-        Piecewise((0, x <= 0), (1, True)))
-    assert Heaviside(x).rewrite(Piecewise, H0=1) == (
-        Piecewise((0, x < 0), (1, True)))
-    assert Heaviside(x).rewrite(Piecewise, H0=S.Half) == (
-        Piecewise((0, x < 0), (S.Half, Eq(x, 0)), (1, x > 0)))
+    assert Heaviside(0).rewrite(SingularityFunction) == S.Half
diff --git a/sympy/printing/tests/test_mathml.py b/sympy/printing/tests/test_mathml.py
index ffd5982609..a2404718a9 100644
--- a/sympy/printing/tests/test_mathml.py
+++ b/sympy/printing/tests/test_mathml.py
@@ -1530,7 +1530,7 @@ def test_print_different_functions():
     assert mpp.doprint(primeomega(x)) == '<mrow><mi>&#x3A9;</mi><mfenced><mi>x</mi></mfenced></mrow>'
     assert mpp.doprint(fresnels(x)) == '<mrow><mi>S</mi><mfenced><mi>x</mi></mfenced></mrow>'
     assert mpp.doprint(fresnelc(x)) ==  '<mrow><mi>C</mi><mfenced><mi>x</mi></mfenced></mrow>'
-    assert mpp.doprint(Heaviside(x)) == '<mrow><mi>&#x398;</mi><mfenced><mi>x</mi></mfenced></mrow>'
+    assert mpp.doprint(Heaviside(x)) == '<mrow><mi>&#x398;</mi><mfenced><mi>x</mi><mfrac><mn>1</mn><mn>2</mn></mfrac></mfenced></mrow>'
 
 
 def test_mathml_builtins():
diff --git a/sympy/printing/tests/test_octave.py b/sympy/printing/tests/test_octave.py
index bdd0dca02e..b39116d460 100644
--- a/sympy/printing/tests/test_octave.py
+++ b/sympy/printing/tests/test_octave.py
@@ -90,7 +90,7 @@ def test_Function_change_name():
     assert mcode(RisingFactorial(x, y)) == "pochhammer(x, y)"
     assert mcode(DiracDelta(x)) == "dirac(x)"
     assert mcode(DiracDelta(x, 3)) == "dirac(3, x)"
-    assert mcode(Heaviside(x)) == "heaviside(x)"
+    assert mcode(Heaviside(x)) == "heaviside(x, 1/2)"
     assert mcode(Heaviside(x, y)) == "heaviside(x, y)"
     assert mcode(binomial(x, y)) == "bincoeff(x, y)"
     assert mcode(Mod(x, y)) == "mod(x, y)"