Repository: sympy/sympy

Issue:
Lambdify slowed down by automated docstring
When experimenting with some mechanics problem I ran into the problem of a long runtime of `lambdify`. Using line profile I found that of the approximate 7 hours of running `lambdify` 84.8% was spend on [expr_str = str(expr)](https://github.com/sympy/sympy/blob/b0a27be36c69dc48d7b2145be82b77d5bffc56a0/sympy/utilities/lambdify.py#L898) (additionally 1.6% was spend on wrapping this string).

Would it be an option to either have a timeout on this line or a keyword argument in `lambdify` to not put the expression in the docstring? Another option would be to have some custom to string function, that just figures out the first 78 characters, as only the first 78 are put into the docstring.

Discussion:

Patch:
diff --git a/doc/src/modules/physics/control/lti.rst b/doc/src/modules/physics/control/lti.rst
index fcd94ef218..5c28064dce 100644
--- a/doc/src/modules/physics/control/lti.rst
+++ b/doc/src/modules/physics/control/lti.rst
@@ -32,5 +32,3 @@ lti
    :members:
 
 .. autofunction:: bilinear
-
-.. autofunction:: backward_diff
diff --git a/sympy/integrals/laplace.py b/sympy/integrals/laplace.py
index 6515db4673..040b4c9886 100644
--- a/sympy/integrals/laplace.py
+++ b/sympy/integrals/laplace.py
@@ -1539,10 +1539,10 @@ def _inverse_laplace_transform(fn, s_, t_, plane, simplify=True, dorational=True
             pass
         elif (r := _inverse_laplace_apply_simple_rules(f, s_, t_)) is not None:
             pass
-        elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
-            pass
         elif (r := _inverse_laplace_apply_prog_rules(f, s_, t_, plane)) is not None:
             pass
+        elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
+            pass
         elif any(undef.has(s_) for undef in f.atoms(AppliedUndef)):
             # If there are undefined functions f(t) then integration is
             # unlikely to do anything useful so we skip it and given an
diff --git a/sympy/physics/control/__init__.py b/sympy/physics/control/__init__.py
index 26d8a87de2..9bdddcb2e3 100644
--- a/sympy/physics/control/__init__.py
+++ b/sympy/physics/control/__init__.py
@@ -1,5 +1,5 @@
 from .lti import (TransferFunction, Series, MIMOSeries, Parallel, MIMOParallel,
-    Feedback, MIMOFeedback, TransferFunctionMatrix, bilinear, backward_diff)
+    Feedback, MIMOFeedback, TransferFunctionMatrix, bilinear)
 from .control_plots import (pole_zero_numerical_data, pole_zero_plot, step_response_numerical_data,
     step_response_plot, impulse_response_numerical_data, impulse_response_plot, ramp_response_numerical_data,
     ramp_response_plot, bode_magnitude_numerical_data, bode_phase_numerical_data, bode_magnitude_plot,
@@ -7,7 +7,7 @@
 
 __all__ = ['TransferFunction', 'Series', 'MIMOSeries', 'Parallel',
     'MIMOParallel', 'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix','bilinear',
-    'backward_diff', 'pole_zero_numerical_data',
+    'pole_zero_numerical_data',
     'pole_zero_plot', 'step_response_numerical_data', 'step_response_plot',
     'impulse_response_numerical_data', 'impulse_response_plot',
     'ramp_response_numerical_data', 'ramp_response_plot',
diff --git a/sympy/physics/control/lti.py b/sympy/physics/control/lti.py
index 71b13803e8..3b87d28c7a 100644
--- a/sympy/physics/control/lti.py
+++ b/sympy/physics/control/lti.py
@@ -22,7 +22,7 @@
 from mpmath.libmp.libmpf import prec_to_dps
 
 __all__ = ['TransferFunction', 'Series', 'MIMOSeries', 'Parallel', 'MIMOParallel',
-    'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix', 'bilinear', 'backward_diff']
+    'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix', 'bilinear']
 
 
 def _roots(poly, var):
@@ -35,8 +35,8 @@ def _roots(poly, var):
 
 def bilinear(tf, sample_per):
         """
-        Returns falling coefficients of H(z) from numerator and denominator.
-        Where H(z) is the corresponding discretized transfer function,
+        Returns falling coeffs of H(z) from numerator and denominator.
+        H(z) is the corresponding discretized transfer function,
         discretized with the bilinear transform method.
         H(z) is obtained from the continuous transfer function H(s)
         by substituting s(z) = 2/T * (z-1)/(z+1) into H(s), where T is the
@@ -44,6 +44,7 @@ def bilinear(tf, sample_per):
         Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned
         as [a, b], [c, d].
 
+
         Examples
         ========
 
@@ -57,10 +58,9 @@ def bilinear(tf, sample_per):
         [2*L + R*T, -2*L + R*T]
         """
 
-
+        z = Symbol('z') # discrete variable z
         T = sample_per  # and sample period T
         s = tf.var
-        z =  s # dummy discrete variable z
 
         np = tf.num.as_poly(s).all_coeffs()
         dp = tf.den.as_poly(s).all_coeffs()
@@ -76,49 +76,6 @@ def bilinear(tf, sample_per):
         return num_coefs, den_coefs
 
 
-def backward_diff(tf, sample_per):
-        """
-        Returns falling coefficients of H(z) from numerator and denominator.
-        Where H(z) is the corresponding discretized transfer function,
-        discretized with the backward difference transform method.
-        H(z) is obtained from the continuous transfer function H(s)
-        by substituting s(z) =  (z-1)/(T*z) into H(s), where T is the
-        sample period.
-        Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned
-        as [a, b], [c, d].
-
-        Examples
-        ========
-
-        >>> from sympy.physics.control.lti import TransferFunction, backward_diff
-        >>> from sympy.abc import s, L, R, T
-        >>> tf = TransferFunction(1, s*L + R, s)
-        >>> numZ, denZ = backward_diff(tf, T)
-        >>> numZ
-        [T, 0]
-        >>> denZ
-        [L + R*T, -L]
-        """
-
-        T = sample_per  # and sample period T
-        s = tf.var
-        z =  s         # dummy discrete variable z
-
-        np = tf.num.as_poly(s).all_coeffs()
-        dp = tf.den.as_poly(s).all_coeffs()
-
-        # The next line results from multiplying H(z) with z^N/z^N
-
-        N = max(len(np), len(dp)) - 1
-        num = Add(*[ T**(N-i)*c*(z-1)**i*(z)**(N-i) for c, i in zip(np[::-1], range(len(np))) ])
-        den = Add(*[ T**(N-i)*c*(z-1)**i*(z)**(N-i) for c, i in zip(dp[::-1], range(len(dp))) ])
-
-        num_coefs = num.as_poly(z).all_coeffs()
-        den_coefs = den.as_poly(z).all_coeffs()
-
-        return num_coefs, den_coefs
-
-
 class LinearTimeInvariant(Basic, EvalfMixin):
     """A common class for all the Linear Time-Invariant Dynamical Systems."""
 
diff --git a/sympy/utilities/lambdify.py b/sympy/utilities/lambdify.py
index 57665c1205..17941fc927 100644
--- a/sympy/utilities/lambdify.py
+++ b/sympy/utilities/lambdify.py
@@ -180,7 +180,7 @@ def _import(module, reload=False):
 
 @doctest_depends_on(modules=('numpy', 'scipy', 'tensorflow',), python_version=(3,))
 def lambdify(args, expr, modules=None, printer=None, use_imps=True,
-             dummify=False, cse=False):
+             dummify=False, cse=False, docstring_limit=1000):
     """Convert a SymPy expression into a function that allows for fast
     numeric evaluation.
 
@@ -350,6 +350,21 @@ def lambdify(args, expr, modules=None, printer=None, use_imps=True,
         When ``True``, ``sympy.simplify.cse`` is used, otherwise (the default)
         the user may pass a function matching the ``cse`` signature.
 
+    docstring_limit : int or None
+        When lambdifying large expressions, a significant proportion of the time
+        spent inside ``lambdify`` is spent producing a string representation of
+        the expression for use in the automatically generated docstring of the
+        returned function. For expressions containing hundreds or more nodes the
+        resulting docstring often becomes so long and dense that it is difficult
+        to read. To reduce the runtime of lambdify, the rendering of the full
+        expression inside the docstring can be disabled.
+
+        When ``None``, the full expression is rendered in the docstring. When
+        ``0`` or a negative ``int``, an ellipsis is rendering in the docstring
+        instead of the expression. When a strictly positive ``int``, if the
+        number of nodes in the expression exceeds ``docstring_limit`` an
+        ellipsis is rendered in the docstring, otherwise a string representation
+        of the expression is rendered as normal. The default is ``1000``.
 
     Examples
     ========
@@ -895,9 +910,14 @@ def _lambdifygenerated(x):
     # Apply the docstring
     sig = "func({})".format(", ".join(str(i) for i in names))
     sig = textwrap.fill(sig, subsequent_indent=' '*8)
-    expr_str = str(expr)
-    if len(expr_str) > 78:
-        expr_str = textwrap.wrap(expr_str, 75)[0] + '...'
+    if _too_large_for_docstring(expr, docstring_limit):
+        expr_str = 'EXPRESSION REDACTED DUE TO LENGTH'
+        src_str = 'SOURCE CODE REDACTED DUE TO LENGTH'
+    else:
+        expr_str = str(expr)
+        if len(expr_str) > 78:
+            expr_str = textwrap.wrap(expr_str, 75)[0] + '...'
+        src_str = funcstr
     func.__doc__ = (
         "Created with lambdify. Signature:\n\n"
         "{sig}\n\n"
@@ -907,7 +927,7 @@ def _lambdifygenerated(x):
         "{src}\n\n"
         "Imported modules:\n\n"
         "{imp_mods}"
-        ).format(sig=sig, expr=expr_str, src=funcstr, imp_mods='\n'.join(imp_mod_lines))
+        ).format(sig=sig, expr=expr_str, src=src_str, imp_mods='\n'.join(imp_mod_lines))
     return func
 
 def _module_present(modname, modlist):
@@ -1394,3 +1414,113 @@ def implemented_function(symfunc, implementation):
             symfunc should be either a string or
             an UndefinedFunction instance.'''))
     return symfunc
+
+
+def _too_large_for_docstring(expr, limit):
+    """Decide whether an ``Expr`` is too large to be fully rendered in a
+    ``lambdify`` docstring.
+
+    This is a fast alternative to ``count_ops``, which can become prohibitively
+    slow for large expressions, because in this instance we only care whether
+    ``limit`` is exceeded rather than counting the exact number of nodes in the
+    expression.
+
+    Parameters
+    ==========
+    expr : ``Expr``, (nested) ``list`` of ``Expr``, or ``Matrix``
+        The same objects that can be passed to the ``expr`` argument of
+        ``lambdify``.
+    limit : ``int`` or ``None``
+        The threshold above which an expression contains too many nodes to be
+        usefully rendered in the docstring. If ``None`` then there is no limit.
+
+    Returns
+    =======
+    bool
+        ``True`` if the number of nodes in the expression exceeds the limit,
+        ``False`` otherwise.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x, y, z
+    >>> from sympy.utilities.lambdify import _too_large_for_docstring
+    >>> expr = x
+    >>> _too_large_for_docstring(expr, None)
+    False
+    >>> _too_large_for_docstring(expr, 100)
+    False
+    >>> _too_large_for_docstring(expr, 1)
+    False
+    >>> _too_large_for_docstring(expr, 0)
+    True
+    >>> _too_large_for_docstring(expr, -1)
+    True
+
+    Does this split it?
+
+    >>> expr = [x, y, z]
+    >>> _too_large_for_docstring(expr, None)
+    False
+    >>> _too_large_for_docstring(expr, 100)
+    False
+    >>> _too_large_for_docstring(expr, 1)
+    True
+    >>> _too_large_for_docstring(expr, 0)
+    True
+    >>> _too_large_for_docstring(expr, -1)
+    True
+
+    >>> expr = [x, [y], z, [[x+y], [x*y*z, [x+y+z]]]]
+    >>> _too_large_for_docstring(expr, None)
+    False
+    >>> _too_large_for_docstring(expr, 100)
+    False
+    >>> _too_large_for_docstring(expr, 1)
+    True
+    >>> _too_large_for_docstring(expr, 0)
+    True
+    >>> _too_large_for_docstring(expr, -1)
+    True
+
+    >>> expr = ((x + y + z)**5).expand()
+    >>> _too_large_for_docstring(expr, None)
+    False
+    >>> _too_large_for_docstring(expr, 100)
+    True
+    >>> _too_large_for_docstring(expr, 1)
+    True
+    >>> _too_large_for_docstring(expr, 0)
+    True
+    >>> _too_large_for_docstring(expr, -1)
+    True
+
+    >>> from sympy import Matrix
+    >>> expr = Matrix([[(x + y + z), ((x + y + z)**2).expand(),
+    ...                 ((x + y + z)**3).expand(), ((x + y + z)**4).expand()]])
+    >>> _too_large_for_docstring(expr, None)
+    False
+    >>> _too_large_for_docstring(expr, 1000)
+    False
+    >>> _too_large_for_docstring(expr, 100)
+    True
+    >>> _too_large_for_docstring(expr, 1)
+    True
+    >>> _too_large_for_docstring(expr, 0)
+    True
+    >>> _too_large_for_docstring(expr, -1)
+    True
+
+    """
+    # Must be imported here to avoid a circular import error
+    from sympy.core.traversal import postorder_traversal
+
+    if limit is None:
+        return False
+
+    i = 0
+    for _ in postorder_traversal(expr):
+        i += 1
+        if i > limit:
+            return True
+    return False


Test patch:
diff --git a/sympy/integrals/tests/test_laplace.py b/sympy/integrals/tests/test_laplace.py
index 9f5a20d603..55b5185b33 100644
--- a/sympy/integrals/tests/test_laplace.py
+++ b/sympy/integrals/tests/test_laplace.py
@@ -483,12 +483,6 @@ def simp_hyp(expr):
     assert ILT(exp(2*s), s, t) == InverseLaplaceTransform(exp(2*s), s, t, None)
     r = Symbol('r', real=True)
     assert ILT(exp(r*s), s, t) == InverseLaplaceTransform(exp(r*s), s, t, None)
-    # Rules for testing whether Heaviside(t) is treated properly in diff rule
-    assert ILT(s**2/(a**2 + s**2), s, t) == (
-        -a*sin(a*t)*Heaviside(t) + DiracDelta(t))
-    assert ILT(s**2*(f(s) + 1/(a**2 + s**2)), s, t) == (
-        -a*sin(a*t)*Heaviside(t) + DiracDelta(t) +
-        Derivative(InverseLaplaceTransform(f(s), s, t, None), (t, 2)))
     # Rules from the previous test_inverse_laplace_transform_delta_cond():
     assert ILT(exp(r*s), s, t, noconds=False) ==\
         (InverseLaplaceTransform(exp(r*s), s, t, None), True)
diff --git a/sympy/physics/control/tests/test_lti.py b/sympy/physics/control/tests/test_lti.py
index 5d0f4b67e2..eeea0d4638 100644
--- a/sympy/physics/control/tests/test_lti.py
+++ b/sympy/physics/control/tests/test_lti.py
@@ -14,8 +14,7 @@
 from sympy.core.containers import Tuple
 from sympy.matrices import ImmutableMatrix, Matrix
 from sympy.physics.control import (TransferFunction, Series, Parallel,
-    Feedback, TransferFunctionMatrix, MIMOSeries, MIMOParallel, MIMOFeedback,
-    bilinear, backward_diff)
+    Feedback, TransferFunctionMatrix, MIMOSeries, MIMOParallel, MIMOFeedback, bilinear)
 from sympy.testing.pytest import raises
 
 a, x, b, s, g, d, p, k, a0, a1, a2, b0, b1, b2, tau, zeta, wn, T = symbols('a, x, b, s, g, d, p, k,\
@@ -1232,14 +1231,3 @@ def test_TransferFunction_bilinear():
     tf_test_manual = TransferFunction(s*T+T, s*(T*b+2*a)+T*b-2*a, s)
 
     assert S.Zero == (tf_test_bilinear-tf_test_manual).simplify().num
-
-def test_TransferFunction_backward_diff():
-    # simple transfer function, e.g. ohms law
-    tf = TransferFunction(1, a*s+b, s)
-    numZ, denZ = backward_diff(tf, T)
-    # discretized transfer function with coefs from tf.bilinear()
-    tf_test_bilinear = TransferFunction(s*numZ[0]+numZ[1], s*denZ[0]+denZ[1], s)
-    # corresponding tf with manually calculated coefs
-    tf_test_manual = TransferFunction(s*T, s*(T*b+a)-a, s)
-
-    assert S.Zero == (tf_test_bilinear-tf_test_manual).simplify().num
diff --git a/sympy/utilities/tests/test_lambdify.py b/sympy/utilities/tests/test_lambdify.py
index c5bad53e0b..a26bd7139b 100644
--- a/sympy/utilities/tests/test_lambdify.py
+++ b/sympy/utilities/tests/test_lambdify.py
@@ -1719,3 +1719,144 @@ def test_23536_lambdify_cse_dummy():
     eval_expr = lambdify(((f, g), z), expr, cse=True)
     ans = eval_expr((1.0, 2.0), 3.0)  # shouldn't raise NameError
     assert ans == 300.0  # not a list and value is 300
+
+
+class LambdifyDocstringTestCase:
+    SIGNATURE = None
+    EXPR = None
+    SRC = None
+
+    def __init__(self, docstring_limit, expected_redacted):
+        self.docstring_limit = docstring_limit
+        self.expected_redacted = expected_redacted
+
+    @property
+    def expected_expr(self):
+        expr_redacted_msg = 'EXPRESSION REDACTED DUE TO LENGTH'
+        return self.EXPR if not self.expected_redacted else expr_redacted_msg
+
+    @property
+    def expected_src(self):
+        src_redacted_msg = 'SOURCE CODE REDACTED DUE TO LENGTH'
+        return self.SRC if not self.expected_redacted else src_redacted_msg
+
+    @property
+    def expected_docstring(self):
+        expected_docstring = (
+            f'Created with lambdify. Signature:\n\n'
+            f'func({self.SIGNATURE})\n\n'
+            f'Expression:\n\n'
+            f'{self.expected_expr}\n\n'
+            f'Source code:\n\n'
+            f'{self.expected_src}\n\n'
+            f'Imported modules:\n\n'
+        )
+        return expected_docstring
+
+    def __len__(self):
+        return len(self.expected_docstring)
+
+    def __repr__(self):
+        return (
+            f'{self.__class__.__name__}('
+            f'docstring_limit={self.docstring_limit}, '
+            f'expected_redacted={self.expected_redacted})'
+        )
+
+
+def test_lambdify_docstring_size_limit_simple_symbol():
+
+    class SimpleSymbolTestCase(LambdifyDocstringTestCase):
+        SIGNATURE = 'x'
+        EXPR = 'x'
+        SRC = (
+            'def _lambdifygenerated(x):\n'
+            '    return x\n'
+        )
+
+    x = symbols('x')
+
+    test_cases = (
+        SimpleSymbolTestCase(docstring_limit=None, expected_redacted=False),
+        SimpleSymbolTestCase(docstring_limit=100, expected_redacted=False),
+        SimpleSymbolTestCase(docstring_limit=1, expected_redacted=False),
+        SimpleSymbolTestCase(docstring_limit=0, expected_redacted=True),
+        SimpleSymbolTestCase(docstring_limit=-1, expected_redacted=True),
+    )
+    for test_case in test_cases:
+        lambdified_expr = lambdify(
+            [x],
+            x,
+            'sympy',
+            docstring_limit=test_case.docstring_limit,
+        )
+        assert lambdified_expr.__doc__ == test_case.expected_docstring
+
+
+def test_lambdify_docstring_size_limit_nested_expr():
+
+    class ExprListTestCase(LambdifyDocstringTestCase):
+        SIGNATURE = 'x, y, z'
+        EXPR = (
+            '[x, [y], z, x**3 + 3*x**2*y + 3*x**2*z + 3*x*y**2 + 6*x*y*z '
+            '+ 3*x*z**2 +...'
+        )
+        SRC = (
+            'def _lambdifygenerated(x, y, z):\n'
+            '    return [x, [y], z, x**3 + 3*x**2*y + 3*x**2*z + 3*x*y**2 '
+            '+ 6*x*y*z + 3*x*z**2 + y**3 + 3*y**2*z + 3*y*z**2 + z**3]\n'
+        )
+
+    x, y, z = symbols('x, y, z')
+    expr = [x, [y], z, ((x + y + z)**3).expand()]
+
+    test_cases = (
+        ExprListTestCase(docstring_limit=None, expected_redacted=False),
+        ExprListTestCase(docstring_limit=200, expected_redacted=False),
+        ExprListTestCase(docstring_limit=50, expected_redacted=True),
+        ExprListTestCase(docstring_limit=0, expected_redacted=True),
+        ExprListTestCase(docstring_limit=-1, expected_redacted=True),
+    )
+    for test_case in test_cases:
+        lambdified_expr = lambdify(
+            [x, y, z],
+            expr,
+            'sympy',
+            docstring_limit=test_case.docstring_limit,
+        )
+        assert lambdified_expr.__doc__ == test_case.expected_docstring
+
+
+def test_lambdify_docstring_size_limit_matrix():
+
+    class MatrixTestCase(LambdifyDocstringTestCase):
+        SIGNATURE = 'x, y, z'
+        EXPR = (
+            'Matrix([[0, x], [x + y + z, x**3 + 3*x**2*y + 3*x**2*z + 3*x*y**2 '
+            '+ 6*x*y*z...'
+        )
+        SRC = (
+            'def _lambdifygenerated(x, y, z):\n'
+            '    return ImmutableDenseMatrix([[0, x], [x + y + z, x**3 '
+            '+ 3*x**2*y + 3*x**2*z + 3*x*y**2 + 6*x*y*z + 3*x*z**2 + y**3 '
+            '+ 3*y**2*z + 3*y*z**2 + z**3]])\n'
+        )
+
+    x, y, z = symbols('x, y, z')
+    expr = Matrix([[S.Zero, x], [x + y + z, ((x + y + z)**3).expand()]])
+
+    test_cases = (
+        MatrixTestCase(docstring_limit=None, expected_redacted=False),
+        MatrixTestCase(docstring_limit=200, expected_redacted=False),
+        MatrixTestCase(docstring_limit=50, expected_redacted=True),
+        MatrixTestCase(docstring_limit=0, expected_redacted=True),
+        MatrixTestCase(docstring_limit=-1, expected_redacted=True),
+    )
+    for test_case in test_cases:
+        lambdified_expr = lambdify(
+            [x, y, z],
+            expr,
+            'sympy',
+            docstring_limit=test_case.docstring_limit,
+        )
+        assert lambdified_expr.__doc__ == test_case.expected_docstring