# Source code for sympy.functions.special.singularity_functions

from __future__ import print_function, division

from sympy.core.function import Function, ArgumentIndexError
from sympy.core import S, sympify, oo, diff
from sympy.core.logic import fuzzy_not
from sympy.core.relational import Eq
from sympy.functions.elementary.complexes import im
from sympy.functions.elementary.piecewise import Piecewise
from sympy.functions.special.delta_functions import DiracDelta, Heaviside

###############################################################################
############################# SINGULARITY FUNCTION ############################
###############################################################################

[docs]class SingularityFunction(Function):
r"""
The Singularity functions are a class of discontinuous functions. They take a
variable, an offset and an exponent as arguments. These functions are
represented using Macaulay brackets as :

SingularityFunction(x, a, n) := <x - a>^n

The singularity function will automatically evaluate to
Derivative(DiracDelta(x - a), x, -n - 1) if n < 0
and (x - a)**n*Heaviside(x - a) if n >= 0.

Examples
========

>>> from sympy import SingularityFunction, diff, Piecewise, DiracDelta, Heaviside, Symbol
>>> from sympy.abc import x, a, n
>>> SingularityFunction(x, a, n)
SingularityFunction(x, a, n)
>>> y = Symbol('y', positive=True)
>>> n = Symbol('n', nonnegative=True)
>>> SingularityFunction(y, -10, n)
(y + 10)**n
>>> y = Symbol('y', negative=True)
>>> SingularityFunction(y, 10, n)
0
>>> SingularityFunction(x, 4, -1).subs(x, 4)
oo
>>> SingularityFunction(x, 10, -2).subs(x, 10)
oo
>>> SingularityFunction(4, 1, 5)
243
>>> diff(SingularityFunction(x, 1, 5) + SingularityFunction(x, 1, 4), x)
4*SingularityFunction(x, 1, 3) + 5*SingularityFunction(x, 1, 4)
>>> diff(SingularityFunction(x, 4, 0), x, 2)
SingularityFunction(x, 4, -2)
>>> SingularityFunction(x, 4, 5).rewrite(Piecewise)
Piecewise(((x - 4)**5, x - 4 > 0), (0, True))
>>> expr = SingularityFunction(x, a, n)
>>> y = Symbol('y', positive=True)
>>> n = Symbol('n', nonnegative=True)
>>> expr.subs({x: y, a: -10, n: n})
(y + 10)**n

The methods rewrite(DiracDelta), rewrite(Heaviside) and rewrite('HeavisideDiracDelta')
returns the same output. One can use any of these methods according to their choice.

>>> 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)
>>> expr.rewrite(DiracDelta)
(x - 4)**5*Heaviside(x - 4) + DiracDelta(x + 3) - DiracDelta(x, 1)
>>> expr.rewrite('HeavisideDiracDelta')
(x - 4)**5*Heaviside(x - 4) + DiracDelta(x + 3) - DiracDelta(x, 1)

========

DiracDelta, Heaviside

Reference
=========

.. [1] https://en.wikipedia.org/wiki/Singularity_function
"""

is_real = True

[docs]    def fdiff(self, argindex=1):
'''
Returns the first derivative of a DiracDelta Function.

The difference between diff() and fdiff() is:-
diff() is the user-level function and fdiff() is an object method.
fdiff() is just a convenience method available in the Function class.
It returns the derivative of the function without considering the chain rule.
diff(function, x) calls Function._eval_derivative which in turn calls
fdiff() internally to compute the derivative of the function.

'''
if argindex == 1:
x = sympify(self.args[0])
a = sympify(self.args[1])
n = sympify(self.args[2])
if n == 0 or n == -1:
return self.func(x, a, n-1)
elif n.is_positive:
return n*self.func(x, a, n-1)
else:
raise ArgumentIndexError(self, argindex)

@classmethod
[docs]    def eval(cls, variable, offset, exponent):
"""
Returns a simplified form or a value of Singularity Function depending on the
argument passed by the object.

The eval() method is automatically called when the SingularityFunction class
is about to be instantiated and it returns either some simplified instance
or the unevaluated instance depending on the argument passed. In other words,
eval() method is not needed to be called explicitly, it is being called
and evaluated once the object is called.

Examples
========
>>> from sympy import SingularityFunction, Symbol, nan
>>> from sympy.abc import x, a, n
>>> SingularityFunction(x, a, n)
SingularityFunction(x, a, n)
>>> SingularityFunction(5, 3, 2)
4
>>> SingularityFunction(x, a, nan)
nan
>>> SingularityFunction(x, 3, 0).subs(x, 3)
1
>>> SingularityFunction(x, a, n).eval(3, 5, 1)
0
>>> SingularityFunction(x, a, n).eval(4, 1, 5)
243
>>> x = Symbol('x', positive = True)
>>> a = Symbol('a', negative = True)
>>> n = Symbol('n', nonnegative = True)
>>> SingularityFunction(x, a, n)
(-a + x)**n
>>> x = Symbol('x', negative = True)
>>> a = Symbol('a', positive = True)
>>> SingularityFunction(x, a, n)
0

"""

x = sympify(variable)
a = sympify(offset)
n = sympify(exponent)
shift = (x - a)

if fuzzy_not(im(shift).is_zero):
raise ValueError("Singularity Functions are defined only for Real Numbers.")
if fuzzy_not(im(n).is_zero):
raise ValueError("Singularity Functions are not defined for imaginary exponents.")
if shift is S.NaN or n is S.NaN:
return S.NaN
if (n + 2).is_negative:
raise ValueError("Singularity Functions are not defined for exponents less than -2.")
if shift.is_negative:
return S.Zero
if n.is_nonnegative and shift.is_nonnegative:
return (x - a)**n
if n == -1 or n == -2:
if shift.is_negative or shift.is_positive:
return S.Zero
if shift.is_zero:
return S.Infinity

def _eval_rewrite_as_Piecewise(self, *args):
'''
Converts a Singularity Function expression into its Piecewise form.

'''
x = self.args[0]
a = self.args[1]
n = sympify(self.args[2])

if n == -1 or n == -2:
return Piecewise((oo, Eq((x - a), 0)), (0, True))
elif n.is_nonnegative:
return Piecewise(((x - a)**n, (x - a) > 0), (0, True))

def _eval_rewrite_as_Heaviside(self, *args):
'''
Rewrites a Singularity Function expression using Heavisides and DiracDeltas.

'''
x = self.args[0]
a = self.args[1]
n = sympify(self.args[2])

if n == -2:
return diff(Heaviside(x - a), x.free_symbols.pop(), 2)
if n == -1:
return diff(Heaviside(x - a), x.free_symbols.pop(), 1)
if n.is_nonnegative:
return (x - a)**n*Heaviside(x - a)

_eval_rewrite_as_DiracDelta = _eval_rewrite_as_Heaviside
_eval_rewrite_as_HeavisideDiracDelta = _eval_rewrite_as_Heaviside