Add a guide on writing custom functions

[asmeurer]

This is still a work in progress, but early reviews are always welcome.

Reminder to reviewers that you can use the PR docs preview by clicking "Click here to see a preview of the documentation." in the checks list below.

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

• other
  • Add a guide to writing custom functions.

[sympy-bot]

✅

Hi, I am the SymPy bot (v167). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• other
  • Add a guide to writing custom functions. (#23488 by @asmeurer)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.11.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

This is still a work in progress, but early reviews are always welcome. 

Reminder to reviewers that you can use the PR docs preview by clicking "Click here to see a preview of the documentation." in the checks list below.

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- other
  - Add a guide to writing custom functions.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[sympy-bot]

🟠

Hi, I am the SymPy bot (v167). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• fb38248:
  • doc/src/guides/custom-functions.md

If these files were added/deleted on purpose, you can ignore this message.

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [77f1d79c]       [f4691d30]
     <sympy-1.10.1^0>                 
+      99.2±0.5ms        180±0.7ms     1.81  sum.TimeSum.time_doit

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

There is potential benefit if you can use evaluate=False and also have printing support. Another point to consider is that maybe it should be possible to call the function with symbols as arguments so the eval method should maybe handle that.

[asmeurer]

But do we really want to recommend that? evaluate=False is very fragile. It's a much better design to just avoid automatic evaluation and move it to some other method such as doit. Even for something like printing, it can be more complicated than it might seem, because that doesn't really tell the whole story. What if you want to manipulate the expression somewhat to get it in some form to print. What if you later decide to do some stuff that isn't printing. The more code you have build up around a function, the harder it is to remove automatic evaluation from it. That's my experience with SymPy functions and I believe it is also true of user defined functions.

And even if you are just doing printing and nothing else wouldn't it still be simpler to just have a function that doesn't evaluate (in fact, you could just use a Function('f') like I described earlier in the guide) than passing evaluate=False everywhere?

Another point to consider is that maybe it should be possible to call the function with symbols as arguments so the eval method should maybe handle that.

I'm not following this. Both a def function and an eval method of a Function subclass should handle symbols as arguments.

[asmeurer]

The only benefit I can see is maybe if we implement #4787, you could have a symbolic function object f which always evaluates. But even then, I don't know if having functions as objects that can be inside of other expressions is something that too many people really need. #4787 is important, but that's the most important part of it. Not many people request that feature.

The reason I am stressing the "don't just return an expression in eval" and also the "use a def function if you just want a simple shortcut" is because it's a common thing I see people do on StackOverflow and the mailing list. I think part of it is that people don't really sit down and think about what they actually want, but also part of it is that they don't completely understand how eval works and maybe don't even realize that just returning an expression in eval has this implication.

Finally, I would hope that everything in this guide should apply equally well to functions in library code as to user defined functions, and indeed I expect it to be useful for contributors. So I'm trying to apply the highest standard to the recommendations here. Users are of course free to do whatever they want, but we should not recommend something that we wouldn't do ourselves. The goal is to replace this wiki page, which isn't written very well https://github.com/sympy/sympy/wiki/Development-Tips-:-About-implementing-special-functions.

[oscarbenjamin]

I agree that automatic evaluation should be discouraged for symbolic functions. My suggestion though is that if someone is in the position of writing a Function subclass with an eval method that always returns maybe they should consider changing it so that the eval method does not always return so that then their function can also work with symbolic arguments.

[asmeurer]

I think which one you really want requires introspection as the whether you're looking for a symbolic function or just a shorthand. I can try to make it clearer here that there are two options.

[oscarbenjamin]

Also perhaps worth mentioning Lambda and maybe even lambdify here:

In [19]: f = Lambda(x, Piecewise((0, Eq(x, 0)), (x+1, True)))

In [20]: f
Out[20]: 
    ⎧  0    for x = 0
x ↦ ⎨                
    ⎩x + 1  otherwise

In [21]: f(0)
Out[21]: 0

[asmeurer]

Mentioning lambdify is a good idea. I thought about Lambda, but I'm a little unsure because I'm not sure just how well supported it is throughout SymPy. Very few things make use of it.

[oscarbenjamin]

The potential advantage of Lambda is that you get a symbolic object that can be used like a function with Python's call syntax (assuming that's what is wanted):

In [1]: eq_y = Lambda(x, Piecewise((0, Eq(x, y)), (1, True)))

In [2]: eq_y
Out[2]: 
    ⎧0  for x = y
x ↦ ⎨            
    ⎩1  otherwise

In [3]: eq_y(y+1)
Out[3]: 1

In [4]: eq_y(y)
Out[4]: 0

Unlike other options the Lambda function can be used with subs to substitute something that is not the "argument" of the function:

In [5]: eq_4 = eq_y.subs(y, 4)

In [6]: eq_4
Out[6]: 
    ⎧0  for x = 4
x ↦ ⎨            
    ⎩1  otherwise

[asmeurer]

The thing I don't understand about Lambda is that it's an Expr. But does it actually make sense to include a Lambda inside of other expressions? You can't call it unless the top-level expression is a Lambda. Does something like Lambda(...) + 1 even make sense?

[oscarbenjamin]

I would prefer it if Lambda wasn't a subclass of Expr. There is a backwards way of calling it:

In [1]: f = Lambda(x, x**2)

In [2]: f
Out[2]: 
     2
x ↦ x 

In [3]: exp(f)
Out[3]: 
      2
 x ↦ x 
ℯ      

In [4]: exp(f).rcall(3)
Out[4]: 
 9
ℯ

I'm not sure I've ever seen that used in the wild though.

[asmeurer]

Hopefully this will change in the future if we can integrate the kind system into Function (#4787, which might change other aspects of this guide as well).

[asmeurer]

GitHub has added $\mathrm{\LaTeX}$ math support to GitHub comments: $$\int_0^1\sin(x), dx.$$ Just use dollar signs, like $$\int_0^1\sin(x)\, dx.$$. Looks like they have a bug with \, and \:. Hopefully that will be fixed soon.

[oscarbenjamin]

GitHub has added
LATEX math

Nice!

(Although apparently it doesn't work with quoting...)

[agryman]

Very cool.

[agryman]

The new LaTeX math support if great!

[bertiewooster]

Suggest splitting this guide into several pages, to decrease page size and length.

[oscarbenjamin]

Maybe not worth mentioning here but perhaps check out the discussion in #20479 and #20493

[oscarbenjamin]

See also #20321

[oscarbenjamin]

It would be much better to have _eval_doit so that subclasses don't need to implement recursion (#17280)

[asmeurer]

Yes definitely. I noticed this problem for several methods. expand does it correctly, but most methods don't.

[oscarbenjamin]

Ideally we wouldn't use eval for type checking so that we can still check types when evaluate=False.

[oscarbenjamin]

Shouldn't these be sympfied

[asmeurer]

Oh I guess so. That's another problem with defining raw methods rather than _eval_* methods.

[oscarbenjamin]

Maybe this should use evaluate=False

[asmeurer]

I don't think there are any cases where this can make a difference. And even if there were, I'm unclear if that would be what you'd want.

[oscarbenjamin]

The idea should be that we don't need to check for evaluation again. I suppose that m and n might now be non-integer and divides could check that.

[asmeurer]

The m and n here are evaluated with doit (when deep=True), so they could be anything. In this example, the doit is just a more advanced version of the stuff in eval, but in general, they might not be related. I guess an optimization that can be made here is to just return self if deep=False since in that case the arguments are the same. But I also don't know if it's that important.

Really, I don't think evaluate=False should ever be used. If you're using it, that's a sign that you have something in eval that should be somewhere else like doit. If evaluate=False is being used to prevent some evaluation, you're creating an object that breaks the function's invariants, and will be fragile (due to the nature of evaluate=False). If evaluate=False is being used just to skip unnecessary work, that's maybe a sign that eval is doing too much. Ideally it should be so fast that it doesn't matter if it gets called again.

[oscarbenjamin]

Perhaps worth explicitly pointing out that argindex uses 1-based indexing rather than 0-based.

[asmeurer]

It's mentioned in the text above but I'll also note it here with a comment. I decided to ignore invalid argindex, since that should never happen from diff and users shouldn't really ever call fdiff (see #23490).

[oscarbenjamin]

Does returning None not work?

[asmeurer]

#23490

[oscarbenjamin]

how to define

[oscarbenjamin]

Creating a new expression like x/pi in an assumptions query is problematic. This kind of thing is bad for performance and also often leads to infinite recursion.

[asmeurer]

Hmm, OK. I thought this was a pretty common pattern, but I can see how it might cause problems.

So is a better pattern to parse the expression (like with the _peeloff_pi helper in the SymPy trig functions)? That seems more complex, especially for end users. What's the simplest way to write this?

[oscarbenjamin]

This kind of thing is more common in the codebase than it should be. There might not really be a better way but there could be some discussion about this in this guide.

The problem is that creating an expression leads to evaluation which leads to assumptions queries. So if assumptions queries create expressions then we have recursion. Whenever you have recursion you need some proof that your recursion terminates. If you know that you are always creating a smaller expression than the one that you started with then that is some proof. Likewise if your assumptions queries only ever query assumptions directly on their args then you know that you are always recursing to a query on a smaller expression tree.

In this case it is probable fine to do this but for a different reason: you are recursing from a query on a versin expression to a query on an expression from the set of expressions that do not involve versin. If you know that the set of expressions that do not involve versin will never create a versin then that is safe.

Getting this right is non-trivial and there should be clear warnings about it. The general rule is:

1. Ideally an assumptions query on x is resolved just by doing some direct assumptions queries on x.args.
2. Assumption queries should aboid creating any new expressions.
3. If it is necessary to create any new expression then you need to think very carefully about whether this recursion is guaranteed to terminate.

In this particular case you could use as_independent to do a structural check that you do actually have a Mul involving pi:

coeff, pi_ = x.as_independent(pi, as_Add=False)
if pi_ is pi:
    return fuzzy_and([x.is_real, fuzzy_not(coeff.is_even)])
elif x.is_real is False:
    return False

Here we always know that as_independent only creates smaller expressions. At worst it returns the original expression as the coeff and 1 for pi_.

[asmeurer]

I see. I do remember seeing issues in the past with that sort of thing. It's not something that we can really show with this example, because versin isn't used anywhere else in the codebase, so it's impossible for it to sneak in. But it can happen if someone is editing assumptions for a more core function in SymPy.

Actually I don't even think it's possible to "think very carefully" about whether the recursion will terminate. It really depends on what assumptions queries are being done elsewhere, and what expressions might create other expressions, which you cannot possibly know unless you've read the whole codebase. It's tough though because there's a lot of advanced logic in other class' assumptions handlers and you really want to be able to reuse that whenever possible.

This really is a major design flaw in the old assumptions. It makes the handlers implement both the specification of what the predicates are and the actual method by which they are resolved. That leads to problems like this where a handler would need to have global information to do that properly, but it can't. It also makes certain optimizations in the inference impossible, even basic ones like not evaluating the same inference twice. Contrast that with the new assumptions where the predicates are specified symbolically and a global solver determines how to resolve everything and which sub-predicates need to be evaluated. Global solvers (like SAT solvers) can handle cyclical inference chains just fine without infinite recursion:

>>> list(satisfiable(And(x >> y, y >> z, z >> x), all_models=True))
[{y: False, z: False, x: False}, {y: True, z: True, x: True}]

[oscarbenjamin]

Actually I don't even think it's possible to "think very carefully" about whether the recursion will terminate

In general no it isn't possible if arbitrary code can be in any of the handlers. However it is possible if we limit what can go into the handlers. Two rules are enough but they need to be applied universally:

1. Assumptions queries should not create new expressions.
2. A query on expr should be resolved only by querying assumptions on the args of expr (and not other queries on expr itself).

With these two rules we always know that any assumptions query only ever recurses to queries on smaller expressions and therefore the recursion must at some point terminate.

This really is a major design flaw in the old assumptions

The flaw is expecting the old assumptions to be something that they are not. They are relatively well defined for their purpose which is to be a limited inference system for use in evaluation. Remember that much of what they do is not really about "assumptions" in the sense of assuming things about symbols but actually just about computing basic properties often of numeric expressions.

[asmeurer]

Review comments have been addressed. I used your example for the _eval_is_positive method. Ad wrote an admonition against creating expressions in assumptions handlers.

It's a little unfortunate that we have to use if/elif instead of just a single fuzzy_* method with your version. Maybe it would help to add fuzzy_if (or fuzzy_implies). It wouldn't evaluate lazily (which is a problem with all the fuzzy_ methods), but it would be more robust to correct None handling.

[oscarbenjamin]

We need to be more careful here:

In [12]: (1 - cos(pi+I*pi))
Out[12]: 1 + cosh(π)

If x is not real then that does not imply that versin(x) is not nonnegative.

[asmeurer]

Ah good catch. Of course, the real range is just [0, 2], so there are complex values that map it to (2, oo).

Really I would like to just write (1 - cos(x)).is_real here. If you look in the guide, a lot of the methods do just delegate to 1 - cos(x). Partly that is just me being lazy/keeping the guide simple. But I also do think there are legitimate reasons why you'd want to do that in general. If something can be reduced to something already implemented somewhere else, you'd like to be able to use it. It's generally going to be more robust to reuse something already in SymPy as a user since it's going to be way better tested, and even for library code, it reduces duplication which also minimizes bugs. So it's really unfortunate to have a "no recursive expressions" rule, which I agree is better for performance and avoiding infinite recursion, but makes it harder to avoid rewriting hard mathematical logic.

But also, one of the use-cases for custom functions that this guide talks about is functions that are shorthands for another function (the see FMA example in the guide). In that case, you have a function that represents another expression and evaluates with doit.

Maybe we just need to have a specific way to create something like this, since the basic pattern is

def Unevaluated(Function):
    @classmethod
    def eval(cls, *args):
        # Evaluate to <expression> on explicit numerical inputs

    def doit(self, deep=True):
        return <expression>

    def _eval_whatever(self):
        return <expression>._eval_whatever()

where the last part is basically the same for every possible method that you would want to define. This whole thing can be wrapped up into a single Function subclass which you could subclass to specify the expression you want to represent, and only override methods if you don't want the default behavior (like defining a custom code printer).

[oscarbenjamin]

makes it harder to avoid rewriting hard mathematical logic.

It does but really I think the problem is the use of classes with automatic evaluation:

1. The mathematical logic is all tied to instances so you need an instance to be able to call the existing code.
2. The instances evaluate whenever you create them and actually most of the mathematical logic is in the evaluation code (so using evaluate=False probably doesn't work).
3. We don't want to run all of that that code every time there is an assumptions query because pretty much anything triggers an assumptions query.

[oscarbenjamin]

See also #23693 where creating expressions breaks assumptions queries under evaluate(False).

[asmeurer]

Made #23701 about the wrapper function superclass idea

[asmeurer]

Any more comments on this?

[oscarbenjamin]

The example above shows that this is not an if and only if.

[oscarbenjamin]

Maybe this is okay but in general I would prefer not to invoke symbolic evaluation as part of numeric evaluation.

[asmeurer]

I'm not following. What is your concern here?

[oscarbenjamin]

variables

[oscarbenjamin]

Maybe it's worth discussing the use of sympify somewhere. I think you've changed it so that both the inputs and outputs of this function are automatically sympified. That's very different from other Basic subclasses so it should probably be clarified somewhere.

[asmeurer]

With things being sympified automatically it's not something you really need to worry about. I'd love for all Basic subclasses and methods to work like that, and also for all classes to more generally use an eval pattern rather than overriding __new__. There are too many bugs caused by people forgetting to sympify inputs/outputs, plus people tend to use sympify instead of _sympify when they do it themselves.

[oscarbenjamin]

I had a couple of comments but I think this looks good in any case.

[oscarbenjamin]

Actually I meant to leave that for you to merge but hit the wrong button...

As I said I think it looks good!

[asmeurer]

Followup at #23732