Multiple dispatch in SymPy

[Ondřej Čertík]

Hi,

Prasoon has described here precisely the problems with extending sympy
with user defined types:

https://groups.google.com/d/topic/sympy/QV6m9Nfp9kw/discussion

Currently the only way to robustly do that is through _op_priority as
Julien described.

The usual Python approach (NotImplemented) doesn't work, because if
for example Mul.__add__ thinks it knows how to work with user's type,
it will not return NotImplemented and thus the user has no way of
overriding it.

The other way to implement this is called "multiple dispatch"
(http://en.wikipedia.org/wiki/Multiple_dispatch). I mention it at the
end of my blog post [1]. Julia has multiple dispatch implemented, see
here how it works there:

http://docs.julialang.org/en/latest/manual/methods/

This is also related to Matthew's idea to use "mul" instead of "Mul":

http://code.google.com/p/sympy/issues/detail?id=3909

since the new function "mul" would (or could) use multiple dispatch to
decide what exactly it will do based on parameters, and I think that
it would just be a quick O(1) dictionary lookup based on argument
types. Then Symbol.__mul__, Integer.__mul__, Add.__mul__, Mul.__mul__,
..., would just call the function "mul".
So it should be fast.

The user would then register his own new types into SymPy and register
his own function for "mul", which would get called whenever he wants.

This would among other things allow to provide faster implementation
of SymPy core, as the faster core (implemented in Cython or C++ for
example) would simply register its own faster types and faster
implementation of "mul". We can also do multiple dispatch for things
like "series" etc., though for functions of just one argument the
default Python method resolution should work well (i.e. the faster
core simply implements its own "series" method). I think the problem
is for functions like "mul", which accept two arguments.

I suggest to read the Julia docs and learn from there all the details
how exactly they decide which method gets called. I think they have
thought about this thoroughly and so we should learn from their
experience.

Ondrej

[1] http://ondrejcertik.blogspot.com/2013/07/my-impressions-from-scipy-2013.html

[Aaron Meurer]

When we've talked about this in the past, we put the ideas at
https://github.com/sympy/sympy/wiki/Canonicalization. So take a look
at what is there, and add things if you have any more ideas.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 10:52 AM, Aaron Meurer <asme...@gmail.com> wrote:
> When we've talked about this in the past, we put the ideas at
> https://github.com/sympy/sympy/wiki/Canonicalization. So take a look
> at what is there, and add things if you have any more ideas.

I've added a link to this thread. Note that the page is about
canonical simplification, which happens after the expression is
constructed. The double dispatch (this thread) happens before any
canonical simplification happens. As correctly noted in the "Resuming
one possible solution taken from mail-list discussions" section in the
wiki.


Ondrej

[Aaron Meurer]

It's related, because currently all canonicalization happens in the
class constructor, and it's impossible for classes to define their own
canonicalization (i.e., what happens in the constructor) due to lack
of dispatch.

Aaron Meurer

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 11:04 AM, Aaron Meurer <asme...@gmail.com> wrote:
> It's related, because currently all canonicalization happens in the
> class constructor, and it's impossible for classes to define their own
> canonicalization (i.e., what happens in the constructor) due to lack
> of dispatch.

That's right. So if we implement multiple dispatch, which btw can be
implemented really easily (though I don't know if this implementation
would be sufficient for our purposes), see e.g.:

http://www.artima.com/weblogs/viewpost.jsp?thread=101605
http://blog.ianbicking.org/more-on-multimethods.html


then
x = Symbol("x")
add(x, x)

would call our standard Add which would do the usual simplification
like x+x -> 2*x.
But then we could implement SymbolHold as a subclass of Symbol, and then

x = SymbolHold("x")
add(x, x)

would dispatch (based on the argument types, which now is SymboHold
instead of Symbol) to a new class AddHold which would keep x+x as x+x.

And we can register AddHold to be called even for argument
combinations like (Symbol, SymbolHold)
or (SymbolHold, Symbol), so if you do:

x = SymbolHold("x")
x2 = Symbol("x")
add(x, x2)

it would still keep it as x+x.

I think this could really work well.

Ondrej

[Aaron Meurer]

I don't think it's so easy, because Add has *args.

Aaron Meurer

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 11:15 AM, Aaron Meurer <asme...@gmail.com> wrote:
> I don't think it's so easy, because Add has *args.

If you are talking about things like:

In [2]: Add(x, x, x, y, x, x)
Out[2]: 5⋅x + y

Then those are of course needed for fast construction many terms into
on Add, but this not really exposed to user's code and it is used
inside things like expand(). As such, we can provide a special static
method or function for doing the same. On the other hand, if you
write:

x + x + x + y + x + x

in Python, then Python itself will call Symbol.__add__ etc., so this
would work well with our new dispatch.

[Ondřej Čertík]

So we just keep our Add unmodified in SymPy. That way, code like expand()
can still call Add(*terms). No problem.

If you want to override the behavior of expand(), then AddHold can
provide a different implementation of expand().

Ondrej

[Aaron Meurer]

No, Add(x, y, z, t) is not just used for fast construction. Any
algorithm that recurses through an expression tree and rebuilds things
will rebuild an Add in that way, using expr.func(*expr.args).

Aaron Meurer

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 11:27 AM, Aaron Meurer <asme...@gmail.com> wrote:
> No, Add(x, y, z, t) is not just used for fast construction. Any
> algorithm that recurses through an expression tree and rebuilds things
> will rebuild an Add in that way, using expr.func(*expr.args).

Which will work with both Add and AddHold, whatever class happens to
be used. I think this is fine.

Ondrej

[Aaron Meurer]

So, going back to what we discussed the first time we met in Los
Alamos, how would you reimplement something like the oo logic so that
it lives entirely in the Infinity class, not in Add.flatten (say for
simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
3*I)?

Aaron Meurer

[Aaron Meurer]

Brian, I think you dropped off list. Forwarding to sympy.

Aaron Meurer

On Wed, Jul 3, 2013 at 2:29 PM, Brian Granger <elli...@gmail.com> wrote:
> I think that having multiple dispatch is something that is needed within
> SymPy. Multiple people have run into the difficulties with providing custom
> Add/Mul logic for their new types. In quantum, we just punted and used the
> builtin Add/Mul logic, but it was a bad compromise to make.
>
> For this, I think looking at Mathematica is really important. This is
> handled by Mathematica's Rules and Patterns which are documented here:
>
> http://reference.wolfram.com/mathematica/guide/RulesAndPatterns.html
> http://reference.wolfram.com/mathematica/tutorial/Introduction-Patterns.html
>
> If we don't have first class, nested and general pattern matching, we can't
> do multiple dispatch. The case in point is classes that accept *args such
> as Add and Mul. Good fast pattern matching would also allow us to do much
> of the logic in Add/Mul constructors using patterns:
>
> Add(x_, x_) -> 2*x_
> Mul(x_, x_) -> x**2
>
> etc.
>
> Summary - simple type based dispatch won't do - we need full blown pattern
> matching and transformation logic.
>
> Cheers,
>
> Brian

[Aaron Meurer]

So, how would you do my example, oo + 3 => oo, with pattern matching?

Aaron Meurer

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
> Brian, I think you dropped off list. Forwarding to sympy.
>
> Aaron Meurer
>
> On Wed, Jul 3, 2013 at 2:29 PM, Brian Granger <elli...@gmail.com> wrote:
>> I think that having multiple dispatch is something that is needed within
>> SymPy. Multiple people have run into the difficulties with providing custom
>> Add/Mul logic for their new types. In quantum, we just punted and used the
>> builtin Add/Mul logic, but it was a bad compromise to make.
>>
>> For this, I think looking at Mathematica is really important. This is
>> handled by Mathematica's Rules and Patterns which are documented here:
>>
>> http://reference.wolfram.com/mathematica/guide/RulesAndPatterns.html
>> http://reference.wolfram.com/mathematica/tutorial/Introduction-Patterns.html
>>
>> If we don't have first class, nested and general pattern matching, we can't
>> do multiple dispatch. The case in point is classes that accept *args such
>> as Add and Mul. Good fast pattern matching would also allow us to do much
>> of the logic in Add/Mul constructors using patterns:
>>
>> Add(x_, x_) -> 2*x_
>> Mul(x_, x_) -> x**2

Pattern matching is another way to implement CAS, but I am actually not
sure this is how Mathematica works inside. Also in pure Python, this might
be really slow (I don't know). I don't see all the little details for
pattern matching
approach. I think I can see more of the details for the sympy approach though.

>>
>> etc.
>>
>> Summary - simple type based dispatch won't do - we need full blown pattern
>> matching and transformation logic.

Why wouldn't simple type based dispatch work?
You might be right, I just want to understand the problem more.



To answer Aaron's question:

On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
> So, going back to what we discussed the first time we met in Los
> Alamos, how would you reimplement something like the oo logic so that
> it lives entirely in the Infinity class, not in Add.flatten (say for
> simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
> 3*I)?

This, and another example is x + O(x). Let's stick to oo + 3.

This is a very good question and it is one of the details that I don't
know the answer 100% yet.
But I feel it is solvable.

The general idea is that you create a new type, let's say Infinity and
register mul, add, div, sub, pow using dispatch, which specify how
your new type should be handled.

In particular, it would look something like:

@dispatch(Infinity, Basic)
def mul(a, b):
if isinstance(b, Integer):
return a
else:
return Add(a, b)

This covers most usecases for "oo", since when you construct it in
Python, you always use the two argument mul(a, b). Also when Add is
constructed, the algorithm simply takes one term at a time and "adds"
it to what is in Add already. So those are all binary operations that
could be overridden in principle.

For example a faster core might decide for speed reasons to not allow
"oo" and O(x) terms. Then multiple dispatch would pretty much make
sure that either the fast core will be used if you do add(x, y), but
slow and more general sympy core will be used if you do add(x, oo).


I think the best would be to create a demo (from scratch) where we can
play with these ideas. Hopefully I'll get to this eventually.

Ondrej

[Brian Granger]

I am not sure what our symbol syntax would look like but something like:

Infinity() + ZeroOrMore(Any()) -> Infinity()
Infinity() + ZeroOrMore(I*Any(x_)) -> Infinity() +ZeroOrMore(I*Any(x_))

Part of what is needed priorities in the matching rules...

--
Brian E. Granger
Cal Poly State University, San Luis Obispo
bgra...@calpoly.edu and elli...@gmail.com

[Brian Granger]

For example take a quantum time evolution operator exp(-I*t*H). What
if I wanted to develop custom logic for that operator. Simply
matching on type type (exp) won't work. You need to be able to say:

exp(-I*Numer(x_)*HermitianOperator(y_))

That requires pattern matching...

>
>
> To answer Aaron's question:
>
> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>> So, going back to what we discussed the first time we met in Los
>> Alamos, how would you reimplement something like the oo logic so that
>> it lives entirely in the Infinity class, not in Add.flatten (say for
>> simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>> 3*I)?
>
> This, and another example is x + O(x). Let's stick to oo + 3.
>
> This is a very good question and it is one of the details that I don't
> know the answer 100% yet.
> But I feel it is solvable.
>
> The general idea is that you create a new type, let's say Infinity and
> register mul, add, div, sub, pow using dispatch, which specify how
> your new type should be handled.
>
> In particular, it would look something like:
>
> @dispatch(Infinity, Basic)
> def mul(a, b):
> if isinstance(b, Integer):
> return a
> else:
> return Add(a, b)

But Mul is not a binary operator in SymPy - we could make it so and
build the overall *args version by repeatedly calling the binary
version. That would be nice... and would somewhat help our
situation...

[Ondřej Čertík]

So first of all, since H does not commute, H will not be a Symbol, but
rather an NCSymbol.
Then the "*" in there will not be Mul but NCMul, thanks to our binary dispatch.
exp can in principle be NCFunction or something.

What kind of custom logic are you thinking of? The automatic canonical
simplification should only be very simple things, most probably in
this case just things that NCMul can do.

Anything more advanced (like perturbation expansion of "exp") should
be implemented as an extra function.

Roughly that's how it works internally, that you take each term and
try to apply/combine it with all the previous ones. The actual
implementation in sympy has gotten quite complicated over the years.

Ondrej

[Brian Granger]

Here is another idea that might help:

What if we *always* mapped the *args logic for operators like Mul/Add
to the binary ones:

if len(args)==2:
self.dispatch(*args)
else:
return Mul(Mul(args[0],args[1]),args[2]) etc

Then we could focus on the rules and patterns for binary versions.
This doesn't solve the issue of needing pattern matching...

[Brian Granger]

But that is besides the point - there are infinitely other examples
where you need to apply pattern matching logic that is not merely type
based.

> What kind of custom logic are you thinking of? The automatic canonical
> simplification should only be very simple things, most probably in
> this case just things that NCMul can do.

I am thinking that the automatic canonical simplification logic
everywhere in sympy would be replaced by rules+patterns.

> Anything more advanced (like perturbation expansion of "exp") should
> be implemented as an extra function.

I completely agree...

[Ondřej Čertík]

Ok, can you give me another one?

>
>> What kind of custom logic are you thinking of? The automatic canonical
>> simplification should only be very simple things, most probably in
>> this case just things that NCMul can do.
>
> I am thinking that the automatic canonical simplification logic
> everywhere in sympy would be replaced by rules+patterns.
>
>> Anything more advanced (like perturbation expansion of "exp") should
>> be implemented as an extra function.
>
> I completely agree...

This extra function can use pattern matching, that's fine.

So we are talking about very minimal canonical simplification, like
H-H ->0, but H*A - A*H stays intact.
I think that two main cases are commutative and non-commutative
symbols. This covers 95% of cases.
Both should be in sympy itself. For the rest you need to create a
special type and dispatch.

Maybe it is possible to implement pattern matching efficiently, but I
am a bit skeptical. The dispatch based on just types on the other hand
should be fast.

Ondrej

[Ronan Lamy]

2013/7/3 Ondřej Čertík <ondrej...@gmail.com>

On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:

Why wouldn't simple type based dispatch work?
You might be right, I just want to understand the problem more.

To answer Aaron's question:

On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
> So, going back to what we discussed the first time we met in Los
> Alamos, how would you reimplement something like the oo logic so that
> it lives entirely in the Infinity class, not in Add.flatten (say for
> simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
> 3*I)?

This, and another example is x + O(x). Let's stick to oo + 3.

x + O(x) is a bad example, because it should really not be represented by an Add.

This is a very good question and it is one of the details that I don't
know the answer 100% yet.
But I feel it is solvable.

I think the best would be to create a demo (from scratch) where we can
play with these ideas. Hopefully I'll get to this eventually.

How about this: https://github.com/rlamy/sympy/commits/binop ?

[Aaron Meurer]

I don't think that magically solves the problem. It just changes the
form of it (and not in a particularly useful way in my opinion). The
issue is that you need for something like oo + I + 3 to work, where
the oo and the 3 are not next to one another. With your idea, it would
be represented as (oo + I) + 3 (or maybe oo + (I + 3)). In either
case, you have to add a complicated rule to dispatch against an Add,
which checks all the arguments.

And by the way, if you're shy on examples, just look at that wiki
page. There are even more that aren't on there.

Aaron Meurer

[Aaron Meurer]

It's tough to mull through a list of commits, so let me just ask you
some questions about it (I know you posted this branch a while ago,
but I forgot the details).

- Does it handle nary operations or just binary?

- What about *args like Mul.flatten?

- If two types register dispatchers against one another, what are the
precedence rules?

Aaron Meurer

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 4:40 PM, Ronan Lamy <ronan...@gmail.com> wrote:
> 2013/7/3 Ondřej Čertík <ondrej...@gmail.com>
>>
>> On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>
>> Why wouldn't simple type based dispatch work?
>> You might be right, I just want to understand the problem more.
>>
>> To answer Aaron's question:
>>
>> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>> > So, going back to what we discussed the first time we met in Los
>> > Alamos, how would you reimplement something like the oo logic so that
>> > it lives entirely in the Infinity class, not in Add.flatten (say for
>> > simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>> > 3*I)?
>>
>> This, and another example is x + O(x). Let's stick to oo + 3.
>
>
> x + O(x) is a bad example, because it should really not be represented by an
> Add.

So the Order class would simply contain both the expression and the
"x", so for example to put this into sympy:

x^2 + x + O(x)

the user would write:

Order(x^2 + x, x)

? I think that's a good idea.

>
>> This is a very good question and it is one of the details that I don't
>> know the answer 100% yet.
>> But I feel it is solvable.
>>
>> I think the best would be to create a demo (from scratch) where we can
>> play with these ideas. Hopefully I'll get to this eventually.
>
>
> How about this: https://github.com/rlamy/sympy/commits/binop ?

Yes! Thanks. Here is how to view changes once you are in this branch:

git diff c84e5df


So I can see that you defined the __pow__ operator in Expr to return
power(a, b) instead of the Power(a, b) class directly. The power(a, b)
is just a function, double dispatched. Then you change all Pow(a, b)
occurrences in sympy to a**b, which gets dispatched to power(a, b)
then. I assume you could have also just changed Power -> power?

Finally power() is then defined as follows:

+...@power.define(Expr, Expr)
+def _power_basecase(x, y):
+ return Pow(x, y)

+...@power.define(Expr, One)
+def _pow_Expr_One(x, one):
+ return x

+...@power.define(One, Expr)
+...@power.define(One, NaN)
+...@power.define(One, One)
+...@power.define(One, Zero)
+def _pow_One_Expr(one, x):
+ return one

etc. (there are some more rules, not important here)

So from this it is clear that power(Expr, One) is used first if
available, otherwise pow(Expr, Expr) is used as a backup plan.

Here are my questions:

* how is performance doing?

* currently your dispatch implementation uses issubclass(c_left, left)
etc., which potentially might be quite slow. Is there any way to just
check the types in a dictionary and only if it is not there, only then
do the slow dispatch that you implemented based on inheritance?

So for example if you put in (Add, One), then on the first run it
would figure out that it should call (Expr, One), and this first run
might be slower, that's ok. But on subsequent runs it would simply
return it from the dictionary directly, so this should be very fast?

Currently we can't use __class__ for it, because:

In [9]: Symbol.__class__
Out[9]: sympy.core.assumptions.ManagedProperties

In [10]: Zero.__class__
Out[10]: sympy.core.singleton.Singleton

In [11]: One.__class__
Out[11]: sympy.core.singleton.Singleton

Due to some sympy metaclasses machinery or something.
But we can create some attribute like _sympy_class_, which would be a
string or a number, unique for each sympy class. User defined types
would write there the name of the class, so that Expr, One, Zero, NaN
would all return unique name.


* What are your other conclusions or impressions from implementing it?

Ondrej

[Aaron Meurer]

On Wed, Jul 3, 2013 at 6:30 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> On Wed, Jul 3, 2013 at 4:40 PM, Ronan Lamy <ronan...@gmail.com> wrote:
>> 2013/7/3 Ondřej Čertík <ondrej...@gmail.com>
>>>
>>> On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>>
>>> Why wouldn't simple type based dispatch work?
>>> You might be right, I just want to understand the problem more.
>>>
>>> To answer Aaron's question:
>>>
>>> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>> > So, going back to what we discussed the first time we met in Los
>>> > Alamos, how would you reimplement something like the oo logic so that
>>> > it lives entirely in the Infinity class, not in Add.flatten (say for
>>> > simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>>> > 3*I)?
>>>
>>> This, and another example is x + O(x). Let's stick to oo + 3.
>>
>>
>> x + O(x) is a bad example, because it should really not be represented by an
>> Add.
>
> So the Order class would simply contain both the expression and the
> "x", so for example to put this into sympy:
>
> x^2 + x + O(x)
>
> the user would write:
>
> Order(x^2 + x, x)
>
> ? I think that's a good idea.

I think he meant it should be some kind of series class. But it
doesn't matter. You still need dispatch to let Add + Series => Series.

Be careful replacing isinstance with something else. I get the need
for speed, but we also want subclasses to be able to do the right
thing.

Aaron Meurer

>
>
> * What are your other conclusions or impressions from implementing it?
>
> Ondrej
>

[Ondřej Čertík]

On Wed, Jul 3, 2013 at 5:11 PM, Aaron Meurer <asme...@gmail.com> wrote:
> On Wed, Jul 3, 2013 at 5:40 PM, Ronan Lamy <ronan...@gmail.com> wrote:
>> 2013/7/3 Ondřej Čertík <ondrej...@gmail.com>
>>>
>>> On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>>
>>> Why wouldn't simple type based dispatch work?
>>> You might be right, I just want to understand the problem more.
>>>
>>> To answer Aaron's question:
>>>
>>> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>> > So, going back to what we discussed the first time we met in Los
>>> > Alamos, how would you reimplement something like the oo logic so that
>>> > it lives entirely in the Infinity class, not in Add.flatten (say for
>>> > simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>>> > 3*I)?
>>>
>>> This, and another example is x + O(x). Let's stick to oo + 3.
>>
>>
>> x + O(x) is a bad example, because it should really not be represented by an
>> Add.
>>
>>> This is a very good question and it is one of the details that I don't
>>> know the answer 100% yet.
>>> But I feel it is solvable.
>>>
>>> I think the best would be to create a demo (from scratch) where we can
>>> play with these ideas. Hopefully I'll get to this eventually.
>>
>>
>> How about this: https://github.com/rlamy/sympy/commits/binop ?
>
> It's tough to mull through a list of commits

See my previous email to easily see the changes.

>, so let me just ask you
> some questions about it (I know you posted this branch a while ago,
> but I forgot the details).
>
> - Does it handle nary operations or just binary?

It only handles Pow.

>
> - What about *args like Mul.flatten?

Pow doesn't have it, so nothing is implemented in this regard.

>
> - If two types register dispatchers against one another, what are the
> precedence rules?

If I understand the code correctly, it raises an exception if two
possibilities exist.
Which is the sane thing to do imho.

[Ronan Lamy]

2013/7/4 Aaron Meurer <asme...@gmail.com>

On Wed, Jul 3, 2013 at 5:40 PM, Ronan Lamy <ronan...@gmail.com> wrote:
> 2013/7/3 Ondřej Čertík <ondrej...@gmail.com>
>>
>> On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>
>> Why wouldn't simple type based dispatch work?
>> You might be right, I just want to understand the problem more.
>>
>> To answer Aaron's question:
>>
>> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>> > So, going back to what we discussed the first time we met in Los
>> > Alamos, how would you reimplement something like the oo logic so that
>> > it lives entirely in the Infinity class, not in Add.flatten (say for
>> > simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>> > 3*I)?
>>
>> This, and another example is x + O(x). Let's stick to oo + 3.
>
>
> x + O(x) is a bad example, because it should really not be represented by an
> Add.
>
>> This is a very good question and it is one of the details that I don't
>> know the answer 100% yet.
>> But I feel it is solvable.
>>
>> I think the best would be to create a demo (from scratch) where we can
>> play with these ideas. Hopefully I'll get to this eventually.
>
>
> How about this: https://github.com/rlamy/sympy/commits/binop ?

It's tough to mull through a list of commits,

You could just look at sympy/core/binop.py, it's quite short.

 

so let me just ask you
some questions about it (I know you posted this branch a while ago,
but I forgot the details).

- Does it handle nary operations or just binary?

Just binary.

- What about *args like Mul.flatten?

It doesn't do anything about it.

- If two types register dispatchers against one another, what are the
precedence rules?

There are no precedence rules. If the dispatcher doesn't find a unique most derived implementation (e.g. if it finds implementations for types (A1, A2) and (B1, B2) such that A1 strictly subclasses B1 and B2 strictly subclasses A2) then it raises an error.

[Aaron Meurer]

OK, this is what you really want to look at
https://github.com/rlamy/sympy/compare/binop.

>
>>
>> so let me just ask you
>> some questions about it (I know you posted this branch a while ago,
>> but I forgot the details).
>>
>> - Does it handle nary operations or just binary?
>
>
> Just binary.

Any thoughts how to extend it to nary?

>>
>>
>> - What about *args like Mul.flatten?
>
>
> It doesn't do anything about it.

Ditto here.

>>
>>
>> - If two types register dispatchers against one another, what are the
>> precedence rules?
>
>
> There are no precedence rules. If the dispatcher doesn't find a unique most
> derived implementation (e.g. if it finds implementations for types (A1, A2)
> and (B1, B2) such that A1 strictly subclasses B1 and B2 strictly subclasses
> A2) then it raises an error.

OK, but there is some kind of precedence on subclasses, right?

Aaron Meurer

[Joachim Durchholz]

Am 03.07.2013 19:15, schrieb Aaron Meurer:
> I don't think it's so easy, because Add has *args.

It should be possible to define the behaviour of the *args case by
declaring Add(a, b, c) to be equivalent to Add(Add(a, b), c).
With the obvious generalization for arbitrary parameter counts.

Doesn't sound difficult to me - what am I overlooking?

[Ronan Lamy]

The algorithm can trivially be extended to higher arities, by doing a triple (or quadruple, or ...) loop. Making it work for arbitrary (fixed) arity wouldn't be hard either. The only issue is that complexity is O(len(mro)**n).

>> - What about *args like Mul.flatten?
>
>
> It doesn't do anything about it.

Ditto here.

 I guess that Mul.flatten should just do the equivalent of reduce(mul, args).

>> - If two types register dispatchers against one another, what are the
>> precedence rules?
>
>
> There are no precedence rules. If the dispatcher doesn't find a unique most
> derived implementation (e.g. if it finds implementations for types (A1, A2)
> and (B1, B2) such that A1 strictly subclasses B1 and B2 strictly subclasses
> A2) then it raises an error.

OK, but there is some kind of precedence on subclasses, right?

No, it just tries all possible combinations of superclasses of the arguments, and finds the minimum of the set of registered implementations for the partial order "is a subclass of". It's hard to be more clever than that and still do the right thing wrt. multiple inheritance.

[Ronan Lamy]

2013/7/4 Joachim Durchholz <j...@durchholz.org>

We can certainly declare that, and it should already be true, but it doesn't help. If a, b, c are symbols, Add(Add(a, b), c) needs to be represented as something that looks like Add(a, b, c), so we're back to square one.

[Aaron Meurer]

It's not enough to just use the mro order? That's what everything else
uses (i.e., attribute lookup), so it seems to me that that should be
the correct way.

But you've probably thought about it more than I have, so I wouldn't
be surprised if I am missing something with it.

[Joachim Durchholz]

Am 03.07.2013 23:55, schrieb Ondřej Čertík:
> Maybe it is possible to implement pattern matching efficiently, but I
> am a bit skeptical.

From what I see in the linked web pages on Mathematica's rule engine,
these patterns are similar to grammar rules for compilers, plus a layer
of tests for applicability of a pattern.

It's possible to merge such patterns into a single super pattern. I.e.
don't test all patterns in sequence, use a state machine that's fed the
next input symbol and label each state with all the patterns that are
still valid at that point; apply a pattern as soon as a state is reached
where only one pattern is left. (Usually, you actually label only the
end states, but I found state machines easier to debug if one could see
the list of still-matching patterns on each node. This helps with
diagnosing ambiguities in the input and with diagnosing bugs in the
pattern matching code.)
This is essentially how lexical analyzers (scanners) are constructed.
Except we would feed tree nodes to the state machine, a scanner feeds
characters.

Setting up the state machine us usually O(number of tokens in rules),
though pathological rule sets can cause exponential size requirements.
Running a tree through such a pattern matcher, on the other hand, should
be linear with number of tree nodes and require no object creation,
which is as fast as it gets.

[Joachim Durchholz]

Am 04.07.2013 03:12, schrieb Ronan Lamy:
> We can certainly declare that, and it should already be true, but it
> doesn't help. If a, b, c are symbols, Add(Add(a, b), c) needs to be
> represented as something that looks like Add(a, b, c), so we're back to
> square one.

It helps in that it eliminates the need to have extra transformation
rules for the *args case.
If the result after transformation contains stuff like ((a*b)*c), it can
still be converted back to Mul(a,b,c).

[Ronan Lamy]

2013/7/4 Ondřej Čertík <ondrej...@gmail.com>

On Wed, Jul 3, 2013 at 4:40 PM, Ronan Lamy <ronan...@gmail.com> wrote:

> 2013/7/3 Ondřej Čertík <ondrej...@gmail.com>
>>
>> On Wed, Jul 3, 2013 at 1:48 PM, Aaron Meurer <asme...@gmail.com> wrote:
>>
>> Why wouldn't simple type based dispatch work?
>> You might be right, I just want to understand the problem more.
>>
>> To answer Aaron's question:
>>
>> On Wed, Jul 3, 2013 at 12:58 PM, Aaron Meurer <asme...@gmail.com> wrote:
>> > So, going back to what we discussed the first time we met in Los
>> > Alamos, how would you reimplement something like the oo logic so that
>> > it lives entirely in the Infinity class, not in Add.flatten (say for
>> > simplicity, oo + 3 should go to oo, but oo + 3*I should remain as oo +
>> > 3*I)?
>>
>> This, and another example is x + O(x). Let's stick to oo + 3.
>
>
> x + O(x) is a bad example, because it should really not be represented by an
> Add.

So the Order class would simply contain both the expression and the
"x", so for example to put this into sympy:

x^2 + x + O(x)

the user would write:

Order(x^2 + x, x)

? I think that's a good idea.

 

It would rather be something like AsymptoticExpansion(expr, order, variable), but the user would still write x**2 + x + O(x). ( O(x) would mean AsymptoticExpansion(0, x, x) )

>
>> This is a very good question and it is one of the details that I don't
>> know the answer 100% yet.
>> But I feel it is solvable.
>>
>> I think the best would be to create a demo (from scratch) where we can
>> play with these ideas. Hopefully I'll get to this eventually.
>
>
> How about this: https://github.com/rlamy/sympy/commits/binop ?

Yes! Thanks. Here is how to view changes once you are in this branch:

git diff c84e5df


So I can see that you defined the __pow__ operator in Expr to return
power(a, b) instead of the Power(a, b) class directly. The power(a, b)
is just a function, double dispatched. Then you change all Pow(a, b)
occurrences in sympy to a**b, which gets dispatched to power(a, b)
then. I assume you could have also just changed Power -> power?

No, because pow(x, 2) needs to return something, i.e. a Pow object. If Pow called pow(), we'd have a circular logic requiring a lot of work to avoid infinite recursion.

Finally power() is then defined as follows:

+...@power.define(Expr, Expr)
+def _power_basecase(x, y):
+    return Pow(x, y)

+...@power.define(Expr, One)
+def _pow_Expr_One(x, one):
+    return x

+...@power.define(One, Expr)
+...@power.define(One, NaN)
+...@power.define(One, One)
+...@power.define(One, Zero)
+def _pow_One_Expr(one, x):
+    return one

etc. (there are some more rules, not important here)

So from this it is clear that power(Expr, One) is used first if
available, otherwise pow(Expr, Expr) is used as a backup plan.

Here are my questions:

* how is performance doing?

I don't know. IIRC, I didn't find any benchmark where the speed of pow() itself made much of a difference.

 

* currently your dispatch implementation uses issubclass(c_left, left)
etc., which potentially might be quite slow. Is there any way to just
check the types in a dictionary and only if it is not there, only then
do the slow dispatch that you implemented based on inheritance?

So for example if you put in (Add, One), then on the first run it
would figure out that it should call (Expr, One), and this first run
might be slower, that's ok. But on subsequent runs it would simply
return it from the dictionary directly, so this should be very fast?

Yes, it would be possible to cache the dispatch. There are two problems with this:

* the cache needs to remain consistent when the dispatch dictionary is updated - the simplest solution being to clear it every time the dict is modified.

* the cache could grow quite big, because we have many classes that can be combined in many ways

The first issue is just a small matter of programming, but the second one may be more fundamental and would need to be tested (particularly with add and mul).

Currently we can't use __class__ for it, because:

In [9]: Symbol.__class__
Out[9]: sympy.core.assumptions.ManagedProperties

In [10]: Zero.__class__
Out[10]: sympy.core.singleton.Singleton

In [11]: One.__class__
Out[11]: sympy.core.singleton.Singleton

Due to some sympy metaclasses machinery or something.
But we can create some attribute like _sympy_class_, which would be a
string or a number, unique for each sympy class. User defined types
would write there the name of the class, so that Expr, One, Zero, NaN
would all return unique name.

What we need to use is the class of objects, not their metaclass. Expr, One, Zero, NaN are unique, hashable objects so we can use them as dictionary keys.

 

* What are your other conclusions or impressions from implementing it?

Ondrej

[Aaron Meurer]

To clarify something here, Symbol, One, and Zero already *are*
classes. You want something like x.__class__ or S(1).__class__.
Getting the class of the class does indeed give you the metaclass
(that's what metaclasses are).

Aaron Meurer

[F. B.]

Did this conversation die out? Multiple dispatch would be very useful...

I was pondering of using multiple dispatch to make the diffgeom and tensor modules easily interact, instead of adding if isinstance ... everywhere.

[Ondřej Čertík]

It looks like the Ronan's approach is the most promising one, so it
needs to be finished. This is quite a big task I think, though
maybe it can be done incrementally.

Ondrej

[Aaron Meurer]

Yes, we still need multiple dispatch, but it's a really hard problem,
and no one has been working on it (I'm not even sure what the best
solution is yet).

Anyway, I personally believe that if you want to refactor the core
with this, that you should clean up the assumptions first, so my plan
is to work on that first.

But feel free to try out some ideas. It might be useful to try some
things out in some of these submodules as a testing ground before
integrating something into the core.

Aaron Meurer

[Joachim Durchholz]

Am 07.11.2013 05:42, schrieb Aaron Meurer:
> Yes, we still need multiple dispatch,

I'm wondering - what do we need that for? What problem does it solve?

> but it's a really hard problem,

Actually, a wrapper-based solution shouldn't be too hard to do, so the
implementation side of things is just the amount of work.

At the conceptual level, it's indeed hard, it involves trade-offs with
nonobvious consequences.
I once tried to work it out for Eiffel, or for statically-typed
languages, and found (to my dismay) that it involves trading off
language properties that I all considered essential. (I do not think
that having a dynamically-typed language like Python helps here; in my
experience, dynamic typing just shifts the problem from the language
designer's shoulders to the library designer's shoulders, with the
library designers having less experience to deal successfully with the
issues.)

However, with a concrete goal in mind and when applying it to a
restricted problem domain, maybe a trade-off can be found that might work.
That's why I'm asking about the problem. Knowing the problem would help
me finding the path of least pain.

> no one has been working on it

Yeah, and I'm not volunteering either since my time is too limited.
What I can do is describe the issues and enable informed design decisions.

> I'm not even sure what the best solution is yet

What trade-offs are on your radar there?

[Aaron Meurer]

On Thu, Nov 7, 2013 at 1:22 AM, Joachim Durchholz <j...@durchholz.org> wrote:
> Am 07.11.2013 05:42, schrieb Aaron Meurer:
>
>> Yes, we still need multiple dispatch,
>
>
> I'm wondering - what do we need that for? What problem does it solve?

See the above discussion.

>
>
>> but it's a really hard problem,
>
> Actually, a wrapper-based solution shouldn't be too hard to do, so the
> implementation side of things is just the amount of work.

The implementation isn't what's hard. Python is dynamic enough that
anything you come up with can be implemented well enough.

The issue is how to actually do it, logically. How should the API
look? How should the logic of it work, especially regarding conflicts
and subclasses. How do we avoid having to define all n*n*k
possibilities (for n classes and k arguments).

>
> At the conceptual level, it's indeed hard, it involves trade-offs with
> nonobvious consequences.
> I once tried to work it out for Eiffel, or for statically-typed languages,
> and found (to my dismay) that it involves trading off language properties
> that I all considered essential. (I do not think that having a
> dynamically-typed language like Python helps here; in my experience, dynamic
> typing just shifts the problem from the language designer's shoulders to the
> library designer's shoulders, with the library designers having less
> experience to deal successfully with the issues.)
>
> However, with a concrete goal in mind and when applying it to a restricted
> problem domain, maybe a trade-off can be found that might work.
> That's why I'm asking about the problem. Knowing the problem would help me
> finding the path of least pain.

Well, F.B. may have his own view, but as far as I'm concerned, I only
care about one thing: I want custom classes to be able to define their
own behavior without having to modify the core. This is important for
two reasons: first, we want our users to be able to do things without
having to monkey-patch or subclass* or literally modify the core
classes, and second, we want to be able to do this within SymPy
itself, so that the core isn't such a huge mess, referencing other
parts of SymPy that it really shouldn't need to know directly about.

(*) - I'm actually not sure if subclassing the core classes is
completely wrong. It seems it actually can be useful in some cases.

>
>
>> no one has been working on it
>
> Yeah, and I'm not volunteering either since my time is too limited.
> What I can do is describe the issues and enable informed design decisions.
>
>
>> I'm not even sure what the best solution is yet
>
> What trade-offs are on your radar there?

At this point, I don't even know how to do it *at all*, so I'm not
even that far. Some things are not clear to me. For example, if you
want to say x defines x + y as one thing and y defines x + y as
another thing, who wins? Does it matter if y is a subclass of x?
Another question, if we have x + y + z and x, y, and z all define
addition behavior, how do we avoid performing O(n**2) computations
just to construct the Add object? Even for fixed argument dispatch,
this is still a concern.

Aaron Meurer

[Ondřej Čertík]

So this particulat case is clear to me: the point of a double dispatch (in
this case) is that you register "x+y" handler for types X + Y, so your handler
just gets called whenever x is a sublass od X and y is a subclass of Y.

> Another question, if we have x + y + z and x, y, and z all define
> addition behavior, how do we avoid performing O(n**2) computations
> just to construct the Add object?

Python decides that first (x+y) handler should by called (see my previous
note), which returns a result "r", and then a handler for "r+z" gets called.

I think that Ronan's code pretty much already implements this.

(Maybe your concern was about something else than I answered.)

Ondrej

[Joachim Durchholz]

Am 08.11.2013 01:40, schrieb Aaron Meurer:
> On Thu, Nov 7, 2013 at 1:22 AM, Joachim Durchholz <j...@durchholz.org> wrote:
>> Am 07.11.2013 05:42, schrieb Aaron Meurer:
>>
>>> Yes, we still need multiple dispatch,
>>
>>
>> I'm wondering - what do we need that for? What problem does it solve?
>
> See the above discussion.

I'm more after whether there's a consensus about what problem we're
actually trying to solve.

>>> but it's a really hard problem,
>>
>> Actually, a wrapper-based solution shouldn't be too hard to do, so the
>> implementation side of things is just the amount of work.
>
> The implementation isn't what's hard. Python is dynamic enough that
> anything you come up with can be implemented well enough.

Agreed.

> The issue is how to actually do it, logically. How should the API
> look? How should the logic of it work, especially regarding conflicts
> and subclasses. How do we avoid having to define all n*n*k
> possibilities (for n classes and k arguments).

Actually it's n^k on the zeroth approximation of n potential classes for
each of the k arguments.
More precisely, it's n1*n2*...*nk for k parameters with n1, n2, ... nk
candidate classes.

That's a hyperrectangle of class combinations that needs to be filled
with associations to code to run.
If you deal with each combination on a case-by-case basis, this becomes
unmaintainable due to sheer numbers.
If you use rules to map many combinations at once, depending on the
structure of the rules you risk overlaps (creating surprising behaviour
because you're usually not aware of the conflict) and oversights
(unmapped combinations). Standard MD rulesets tend to exclude
oversights; I'm somewhat sceptical that that's a good idea, an oversight
produces a clearly visible and easily fixable error message while an
overlap silently creates unexpected behaviour (which might not be wrong
but still using the wrong, too-slow algorithm or be otherwise undesirable).

>> However, with a concrete goal in mind and when applying it to a restricted
>> problem domain, maybe a trade-off can be found that might work.
>> That's why I'm asking about the problem. Knowing the problem would help me
>> finding the path of least pain.
>
> Well, F.B. may have his own view, but as far as I'm concerned, I only
> care about one thing: I want custom classes to be able to define their
> own behavior without having to modify the core. This is important for
> two reasons: first, we want our users to be able to do things without
> having to monkey-patch or subclass* or literally modify the core
> classes, and second, we want to be able to do this within SymPy
> itself, so that the core isn't such a huge mess, referencing other
> parts of SymPy that it really shouldn't need to know directly about.

As far as I know it's an unsolved problem.
And I have read a lot of research papers on the issue. I found that all
solved either the conflict or the oversight, but you'd need to solve
both, and I suspect it's just as unsolvable as the halting problem.

It's been a while since I last looked into the state of research (a
decade or so); maybe somebody found an entirely different angle to solve
that problem.

>> What trade-offs are on your radar there?
>
> At this point, I don't even know how to do it *at all*, so I'm not
> even that far. Some things are not clear to me. For example, if you
> want to say x defines x + y as one thing and y defines x + y as
> another thing, who wins? Does it matter if y is a subclass of x?
> Another question, if we have x + y + z and x, y, and z all define
> addition behavior, how do we avoid performing O(n**2) computations
> just to construct the Add object? Even for fixed argument dispatch,
> this is still a concern.

That's the combinatoric explosion of the hyperrectangle.

The design space is that you can have any three of the following properties:
- Allow end users to extend the hyperrectangle however they please
- Have no conflicts
- Have no oversights
- No combinatoric explosion

There might be a useful compromise between these goals somewhere, but I
don't know it.
Those compromises that I read about all came with more complex base
rules, which had the undesirable effect that it became difficult to
figure out which combination of parameters would result in the call to
what implementation.

[Joachim Durchholz]

Am 08.11.2013 01:40, schrieb Aaron Meurer:

> At this point, I don't even know how to do it *at all*, so I'm not
> even that far. Some things are not clear to me. For example, if you
> want to say x defines x + y as one thing and y defines x + y as
> another thing, who wins? Does it matter if y is a subclass of x?
> Another question, if we have x + y + z and x, y, and z all define
> addition behavior, how do we avoid performing O(n**2) computations
> just to construct the Add object? Even for fixed argument dispatch,
> this is still a concern.

I forgot to mention that commutativity and (in some sense) associativity
further complicate issues, making the effects of any ruleset even more
difficult to grasp.
I.e. I didn't ignore that, my previous message was just too long already.

Let me state my stance: MD is hard to do well, and therefore a bad idea
if there are alternatives.
BUT the logic issues are the same if it's inside a rule engine that does
structural pattern matching. So I'm taking back my knee-jerk reaction of
"extending a programming language with MD is always a bad idea" - it
still is, but in this particular case, the emphasis on "programming
language" doesn't solve any problem.
MD is still a bad idea, it's just that there seems to be no good
alternative.

Approaches that I see:
- Cut down on MD as far as possible. Usually by defining conversion
chains (e.g. Integer->Rational->Real->Complex, BUT ON A PER-FUNCTION
BASIS - a Real->Complex conversion would not preserve semantics for
computing whether a function has a derivative, for example), and
possibly by leveraging commutativity and associativity.
- Excel at diagnostics. Have a utility that spits out the full dispatch
table for any function. For functions with more than two parameters,
spit out slices, slicing the hyperrectangle at the dimensions that give
the fewest number of slices.
- Have each function define what parameters are subject to MD and what
aren't, cutting down on the number of dimensions.

More generally, keep the MD information (which functions have MD, and
the actual dispatch tables) in data structures that aren't intimately
tied to Python wrappers. We need that information for three cases:
- Pythonic multiple dispatch
- Diagnostics
- Rule-based tree transformations (once this actually happens)

[Aaron Meurer]

> On Nov 8, 2013, at 2:26 AM, Joachim Durchholz <j...@durchholz.org> wrote:
>
> Am 08.11.2013 01:40, schrieb Aaron Meurer:
>> On Thu, Nov 7, 2013 at 1:22 AM, Joachim Durchholz <j...@durchholz.org> wrote:
>>> Am 07.11.2013 05:42, schrieb Aaron Meurer:
>>>
>>>> Yes, we still need multiple dispatch,
>>>
>>>
>>> I'm wondering - what do we need that for? What problem does it solve?
>>
>> See the above discussion.
>
> I'm more after whether there's a consensus about what problem we're actually trying to solve.
>
>>>> but it's a really hard problem,
>>>
>>> Actually, a wrapper-based solution shouldn't be too hard to do, so the
>>> implementation side of things is just the amount of work.
>>
>> The implementation isn't what's hard. Python is dynamic enough that
>> anything you come up with can be implemented well enough.
>
> Agreed.
>
>> The issue is how to actually do it, logically. How should the API
>> look? How should the logic of it work, especially regarding conflicts
>> and subclasses. How do we avoid having to define all n*n*k
>> possibilities (for n classes and k arguments).
>
> Actually it's n^k on the zeroth approximation of n potential classes for each of the k arguments.
> More precisely, it's n1*n2*...*nk for k parameters with n1, n2, ... nk candidate classes.

Ah. For some reason I always think n*k when it should be n^k.

>
> That's a hyperrectangle of class combinations that needs to be filled with associations to code to run.
> If you deal with each combination on a case-by-case basis, this becomes unmaintainable due to sheer numbers.
> If you use rules to map many combinations at once, depending on the structure of the rules you risk overlaps (creating surprising behaviour because you're usually not aware of the conflict) and oversights (unmapped combinations). Standard MD rulesets tend to exclude oversights; I'm somewhat sceptical that that's a good idea, an oversight produces a clearly visible and easily fixable error message while an overlap silently creates unexpected behaviour (which might not be wrong but still using the wrong, too-slow algorithm or be otherwise undesirable).

Actually, for SymPy, if there's no rule, then the behavior is to just
create the unevaluated object. This will actually be the common case.
x + y will return Add(x, y) in almost every case. Only if x and y are
similar objects (like if they differ only by a number multiple) or
some special object like O(x**2) or oo will anything special happen.
So the hyperrectangle will be mostly sparse, and that's fine. It's not
like with numeric libraries where x + y has to return some kind of
number and so you have to catch every possible combination. We have a
default behavior: leave it unevaluated.

Aaron Meurer

[Joachim Durchholz]

Am 09.11.2013 06:15, schrieb Aaron Meurer:
>> On Nov 8, 2013, at 2:26 AM, Joachim Durchholz <j...@durchholz.org> wrote:
>> That's a hyperrectangle of class combinations that needs to be filled with associations to code to run.
>> If you deal with each combination on a case-by-case basis, this becomes unmaintainable due to sheer numbers.
>> If you use rules to map many combinations at once, depending on the structure of the rules you risk overlaps (creating surprising behaviour because you're usually not aware of the conflict) and oversights (unmapped combinations). Standard MD rulesets tend to exclude oversights; I'm somewhat sceptical that that's a good idea, an oversight produces a clearly visible and easily fixable error message while an overlap silently creates unexpected behaviour (which might not be wrong but still using the wrong, too-slow algorithm or be otherwise undesirable).
>
> Actually, for SymPy, if there's no rule, then the behavior is to just
> create the unevaluated object. This will actually be the common case.

Ah okay, so then we'll never have what I've been calling oversights.
That's helpful.

Hm. Well, on second thinking, there can still be oversights. Cases where
people find that case AxB works, case BxA doesn't but should because
it's symmetric. Or where Complex isn't covered but should be done by the
Real case.
But maybe it's good enough anyway. It's certainly what's already
happening, so we're not getting worse that way.

> x + y will return Add(x, y) in almost every case. Only if x and y are
> similar objects (like if they differ only by a number multiple) or
> some special object like O(x**2) or oo will anything special happen.
> So the hyperrectangle will be mostly sparse, and that's fine.

You do want some of the hyperrectangle entries to project.
E.g. you'll want to specify that oo * Number to return oo, and you don't
want to repeat yourself.
The issue to solve is two cases here:
1) Exceptional cases. One may be tempted to specify oo as the result of
oo*Number, but that's not the right one for oo*Complex, nor is it the
right one for oo*Zero. However, that's just "competence in math", I do
not think that any Python code can help with that (unless we rebuild
SymPy around a proof checker, which would be interesting but essentially
a full rewrite).
2) Conflicts. Somebody adds (say) Quaternions, somebody else adds a rule
for One, and suddenly we have both One*NumberLike and Number*Quaternion.
This can be diagnosed in unit tests if the association of type
combinations to functions to call is kept as a data structure (as
opposed to coding the decision tree into a piece of Python code).

> It's not
> like with numeric libraries where x + y has to return some kind of
> number and so you have to catch every possible combination. We have a
> default behavior: leave it unevaluated.

We have, but we still have the situation where the default behaviour is
inappropriate because it results in missed simplification opportunities.
And that's a biggie, because missed simplification means many more
integrals and integral equations that can't be solved, and such stuff.

Fortunately, that's not an architectural error that would be unfixable
without a complete rewrite, we can fix missed opportunities over time.

I like where this is heading.

[F. B.]

Concerning cases where there is no ambiguity in their usage, what about using this library?

https://bitbucket.org/coady/multimethod/src

I did not find the license, so I suppose it's public domain (correct me if I'm wrong).

There are some cases where the problems being discussed do not appear, I think this simple multimethod could be already used.

[Joachim Durchholz]

> think this simple *multimethod* could be already used.

Does it have a way to inspect the mappings?
That "create as needed" bit doesn't sound like it does, but I haven't
looked too deeply into the code to be sure.

I'm a bit concerned about that MRO usage.
MRO as a concept is fundamentally broken. That brokenness doesn't matter
as long as MRO doesn't matter, i.e. as long as no multiple inheritance
is in play - and I think (I hope!) that SymPy does not use it.

[Aaron Meurer]

On Sat, Nov 9, 2013 at 1:39 PM, F. B. <franz....@gmail.com> wrote:
> Concerning cases where there is no ambiguity in their usage, what about
> using this library?
>
> https://bitbucket.org/coady/multimethod/src
>
> I did not find the license, so I suppose it's public domain (correct me if
> I'm wrong).

No, the "default" license of things is no license. Unless someone
explicitly gives you permission to use their work, you are assumed not
to have it.

With that being said, almost all public code on GitHub and Bitbucket
is open source, as that's generally what those sites are used for, so
probably the author did intend for that code to be open source. Likely
you can just email him and ask him if he will provide it under
something BSD compatible if it's something that you want to use. I
didn't play with it, so I can't say if it actually is at this point.

Aaron Meurer

>
> There are some cases where the problems being discussed do not appear, I
> think this simple multimethod could be already used.
>

[Aaron Meurer]

Actually, according to the setup.py,
https://bitbucket.org/coady/multimethod/src/3179d8ffa4c6c920d7ec6acbcd125489bb03a359/setup.py?at=default#cl-16,
it uses PSF license, which is fine (I think there may be some
bookkeeping clause, but it's completely compatible with BSD).

Aaron Meurer

[Joachim Durchholz]

Am 09.11.2013 23:46, schrieb Aaron Meurer:
> On Sat, Nov 9, 2013 at 1:39 PM, F. B. <franz....@gmail.com> wrote:
>> Concerning cases where there is no ambiguity in their usage, what about
>> using this library?
>>
>> https://bitbucket.org/coady/multimethod/src
>>
>> I did not find the license, so I suppose it's public domain (correct me if
>> I'm wrong).
>
> No, the "default" license of things is no license. Unless someone
> explicitly gives you permission to use their work, you are assumed not
> to have it.

I just opened an issue so that at least that one does not get into the way.

Though I guess the code is small enough that it can be recreated.

[F. B.]

On Saturday, November 9, 2013 10:07:50 PM UTC+1, Joachim Durchholz wrote:

I'm a bit concerned about that MRO usage.
MRO as a concept is fundamentally broken. That brokenness doesn't matter
as long as MRO doesn't matter, i.e. as long as no multiple inheritance
is in play - and I think (I hope!) that SymPy does not use it.

I just checked, Add inherits both Expr and AssocOp, both of which on the other hand inherit Basic. So, unfortunately, SymPy does use multiple inheritance and presents the diamond problem, which is likely to be solved by Python's C3 algorithm.

I have some ideas in my mind about a possible suggestion. In the next days I'll write again.

[Joachim Durchholz]

Am 10.11.2013 21:34, schrieb F. B.:
>
> I just checked, *Add* inherits both *Expr* and *AssocOp*, both of which on
> the other hand inherit *Basic*. So, unfortunately, SymPy does use multiple

> inheritance and presents the diamond problem, which is likely to be solved
> by Python's C3 algorithm.

Hm. Not a problem if the same function isn't overridden in both Expr and
AssocOp, so maybe... hopefully...

... but yes, that's going to have to be considered for multiple dispatch.

[F. B.]

On Sunday, November 10, 2013 10:42:39 PM UTC+1, Joachim Durchholz wrote:

Hm. Not a problem if the same function isn't overridden in both Expr and
AssocOp, so maybe... hopefully...

... but yes, that's going to have to be considered for multiple dispatch.

I was thinking about that. That AssocOp is interesting, it defines the commutative property of Add.

The fact is, the Add object has two logical meanings:

• represent the plus "+" sign
• encode the properties of the addition

Actually, the properties of the addition can vary according to what objects are inside of Add. If you have just numbers, then you get a number, with symbols you get a summation of symbols, with infinities an complex numbers, more complex behavior.

In other parts of SymPy there are MatAdd, TensAdd to represent additions of matrices, tensors, etc..., which have their own properties. Developer deemed that it would be too complicated to add rules to Add, so they wrote other Add-like objects. In other cases, like the aforementioned infinities, Add is still kept, by adding behavioural rules to it.

But, isn't this breaking the logical correspondence of those two points in the bullet list? The plus sign is no more mapped to a single object, while the properties of addition are sometimes handled by Add, sometimes by other objects.

What about using instead of this an approach closer to C++ templates? I mean, the Add class would always represent a plus sign, for every expression, while its properties are managed by templates, so, keeping a C++ syntax:

• Add<AssocOpT>  ordinary
  
• Add<InfinitiesT>  for additions containing infinities.
  
• Add<MatrixT>  for additions containing matrices.
• Add<NonCommutativeT>  for additions containing non-commutative symbols.
• Add<NonCommutativeInfinitieT>  for both inf and non-commutative symbols.
• Add<QuantumExpressionT>  for the case mentioned by Brian, thus avoiding usage of pattern matching.
  

In this case, the properties and behaviour of the arguments of Add are managed by special classes.

The point is, the addition of new objects to an Add argument list, could change the template of the resulting new Add, e.g. Add(x, y) ==> Add<AssocOpT>(x, y), while Add(x, y) + oo ==> Add<InfinitiesT>(x, y, oo).

At this point, multiple dispatch could be defined so as to keep templates in regard, but I'm still unsure how this could work, but I made up my mind to post this anyway. I'm not even sure that properties should be regarded as templates, but, in any case, using such an external marker given immediately information about special symbols inside Add.

So, in quantum mechanics, the exp( ) function could immediately catch that its Mul/Add content contains quantum mechanical operators, thus acting appropriately.

[F. B.]

On Sunday, November 10, 2013 11:21:46 PM UTC+1, F. B. wrote:

I was thinking about that. That AssocOp is interesting, it defines the commutative property of Add.

Sorry, the associative, not commutative property.

[Aaron Meurer]

I think you've hit on the core of the issue. The problem is, if you
take just Symbol('x'), it can be interpreted as many things. x + 1
means the number x plus the number 1 (or does it mean the function
f(x) = x + 1?). x + O(x**2) means the set of functions whose series
expansion up to x is x. x + y could mean that we think of x and y as
matrices of the same shape.

The thing that really has bothered me is that it's hard to state
clearly, what exactly should be a subclass of Expr. The best
definition I've come up with is, "those things that make sense with
Add and Mul". I can tell you what's *not* Expr. Tuple is not Expr.
Boolean expressions (anything that can be thought of as having a truth
value) are not Expr.

Mathematically, + has a thousand meanings, and we blur the
distinctions because they all act similarly. But at some point in a
CAS, you have to unblur any notational convenience that you would
otherwise make on paper.

Aaron Meurer

[F. B.]

On Sunday, November 10, 2013 11:31:38 PM UTC+1, Aaron Meurer wrote:

Mathematically, + has a thousand meanings, and we blur the
distinctions because they all act similarly. But at some point in a
CAS, you have to unblur any notational convenience that you would
otherwise make on paper.

In Wolfram Mathematica, + is a Plus M-expression, always. They don't have any MatPlus, or whatever.

The behavior upon it is determined by pattern matching. I don't know how it works, but I guess that Mathematica builds up some indices to speed up the pattern search.

I don't think pattern matching is the only way. I would suggest to make Add, Mul and Power just containers for the corresponding operators (+, *, **), while creating parallel classes to represent the precise logical operation that those signs are meant to represent.

Otherwise, things could be done the opposite way, create an Add-like object for every logical operation, which means making Add abstract, and along with already existing classes like MatAdd, add stuff like: AssocAdd, InfinitiesAdd, NCAdd, InfinitiesNonCommutativeAdd, NumericAdd, QuantumOpAdd and so on. I don't think this is a good idea, too many classes to define.

[Aaron Meurer]

Yes, I agree this is the way to do it. It's essentially the idea from
https://code.google.com/p/sympy/issues/detail?id=2738. But again, I'm
really not sure exactly *how* to do it, down to the implementation
details. So I really feel (at least for myself) that I don't quite
understand the abstractions and the proper separation of ideas yet.
It's a hard problem, and as I mentioned earlier, not one that I plan
on working on myself in the foreseeable future (I'm already working on
another hard problem, namely the assumptions system).

Aaron Meurer

[Joachim Durchholz]

Am 11.11.2013 00:08, schrieb F. B.:
>
> In Wolfram Mathematica, + is a Plus M-expression, always. They don't have
> any MatPlus, or whatever.
>

> The behavior upon it is determined by pattern matching. [...]

>
> I don't think pattern matching is the only way.

Actually, it is.
*Somewhere*, you need to match expression structure with what to do. The
various *Add classes are essentially just precomputed pattern matches.

The more interesting question is how to organize the patterns.

Do we do it by application domain, i.e. one set of rules for infinities,
one set for numbers, one for functions, etc.? That's weak when dealing
with the cases that span application domains.

Or do we do it by operation? That's unsure footing because operations
tend to be defined in terms of application domain so it's a bit hard to
agree what "+" actually is; also, it tends to create huge, messy
specifications as the definition for "+" must cater for all combinations
but the kitchen sink.

I'm leaning towards a mixed approach: define groups of application
domains and do it on a per-function basis inside; however, now the
question is how to define groups, and the best answer I can come up for
that is "we'd need to experiment and see what works", which is hardly
satisfactory.

[Joachim Durchholz]

Am 10.11.2013 23:21, schrieb F. B.:
>
> But, isn't this breaking the logical correspondence of those two points in
> the bullet list? The plus sign is no more mapped to a single object, while
> the properties of addition are sometimes handled by Add, sometimes by other
> objects.

Not 100% sure that we're in the same boat here.
My perspective is that since we need to have a way to see all the cases
so that we can check for omissions and conflicts, "+" as such needs to
be represented in some way. (In theory, it doesn't even need to be
represented as data, it could be a piece of code.)
However, since this is getting large and unwieldy, I'd want to see the
dispatch table held somewhere centrally but the implementations could
(and probably should) be held in separate modules.

> What about using instead of this an approach closer to C++ templates?

The problem with C++ is that many people do not know C++, and those that
believe they do sometimes misinterpret the definitions.
For that reason, it is usually better to stick with Python nomenclature
and concepts.

> I mean, the Add class would always represent a plus sign, for every
> expression, while its properties are managed by templates,

Here, the question would be how you'd represent those templates would be
in Python.
I do not think that the proposal can be properly understood without that.

[F. B.]

On Monday, November 11, 2013 8:52:40 PM UTC+1, Joachim Durchholz wrote:

> What about using instead of this an approach closer to C++ templates?

The problem with C++ is that many people do not know C++, and those that
believe they do sometimes misinterpret the definitions.
For that reason, it is usually better to stick with Python nomenclature
and concepts.

OK, we don't have to use templates as they are defined in C++. It could be simply a link to an addition/multiplication rules object inside the class.

 > I mean, the Add class would always represent a plus sign, for every
> expression, while its properties are managed by templates,

Here, the question would be how you'd represent those templates would be
in Python.
I do not think that the proposal can be properly understood without that.

What about the __getitem__ of Add and Mul? So we could have

Add[AssocOp](x, y, z) ==> x + y + z
w = var('w', commutative=False)
Add[NonCommutative](x, w) ==> x + w

and the multiple dispatch could look like:

@dispatch(Add[NonCommutative], Add[InfinitiesT])
def f(x, y):
    # This function is called only when x is an Add with
    # non-commutative elements, while x is
    # an Add containing infinities.
    ...

Of course the hard part is not this, it's reorganizing the rules inside these template-like objects.

[F. B.]

Obviously, Add could change its template according to the cases:

Add(x, y) ===> Add[AssocOp](x, y)
Add(x, y) + Add(oo,1) ===> Add(Add[AssocOp](x,y), Add[InfinitiesT](oo))
    ===> Add[InfinitiesT](x,y,oo)

Add's __new__ constructor could be dispatched to handle specific cases for these templates.

[Aaron Meurer]

I'm not a fan of the syntax, but otherwise I think you have the right idea.

Aaron Meurer

[Ondřej Čertík]

Hi,

So there is this whole Parent/coercion machinery in Sage:

http://www.sagemath.org/doc/tutorial/tour_coercion.html
http://www.sagemath.org/doc/reference/coercion/

Isn't this what you are trying to design?

Ondrej

[F. B.]

I had a try at this with Julia... unfortunately there is a problem there because parametric type constructors are invariant, while they should be covariant for this particular case:

julia> immutable Point{T}
          x::T
          y::T
       end

julia> p = Point{Int64}(3, 5)
Point{Int64}(3,5)

julia> p2 = Point{Float64}(3.0, 5.0)
Point{Float64}(3.0,5.0)

julia> p.x
3

julia> p.y
5

julia> f(v::Point{Float64}) = "Point... parametric type Float64"
f (generic function with 3 methods)

julia> f(v::Point) = "Point... no parametric type"
f (generic function with 1 method)

julia> f(p2)
"Point... parametric type Float64"

julia> f(p)
"Point... no parametric type"

julia> f(v::Point{Int64}) = "Point... parametric type Int64"
f (generic function with 2 methods)

julia> f(p)
"Point... parametric type Int64"

If methods with detail are not matched, the more generic one is matched:

julia> p3 = Point{String}("hello", "world")
Point{String}("hello","world")

julia> f(p3)
"Point... no parametric type"

The problem is here, parametric types break all inheritance relations:

julia> g(v::Point{Integer}) = "Point... parametric type is Integer abstract class"
g (generic function with 1 method)

julia> g(p)
ERROR: no method g(Point{Int64})

In Julia may approach is not possible, or at least it has to be handled in another way.

Another language, Scala, has a similar behavior as Julia, but allows to put markers on parametric type to keep inheritance relations for derived types, + for covariant, - for invariant.

In my suggestion, this problem should be considered.

[Aaron Meurer]

That is reminiscent of how things are implemented in the polys module.
It's not clear to me that those ideas can extend to less rigorously
defined rings and fields (at least not without first rigorously
defining them).

Aaron Meurer

[F. B.]

Yet another alternative approach:

Instead of defining Add and Mul as dependent on a parametric type associated with the properties of the expression, just do the opposite.

Examples (Julia-like syntax, {T} is the parametric type, <: means inheritance):

type Expr{T} <: Basic
type NonCommutativeExpr{T} <: Expr{T}
type InfinitiesExpr{T} <: Expr{T}
type MatrixExpr{T} <: Expr{T}
type NonCommutativeInfinitiesExpr{T} <: Expr{T}  # or better use multiple inheritance?

Now, their corresponding Add, Mul, Pow classes are given by

NonCommutativeExpr{Add}
NonCommutativeExpr{Mul}
MatrixExpr{Add}  # instead of MatAdd
MatrixExpr{Mul}  # instead of MatMul

and so on.

This has the advantage to require developers to subclass the expression only once, while the corresponding Add, Mul, ... objects are automatically created by parametric typing. There is no covariance-invariance problem for the parametric type, because it is not using the template, it's using the main type in the dispatch, so this is compatible with lots of other programming languages.

[Ondřej Čertík]

I wonder how difficult would be to write just a very simple core in
Julia, just so that we can experiment with this.
See how extensible it is, and so on.

Ondrej

[F. B.]

On Wednesday, November 13, 2013 10:39:48 PM UTC+1, Ondřej Čertík wrote:

I wonder how difficult would be to write just a very simple core in
Julia, just so that we can experiment with this.
See how extensible it is, and so on.

The problem is how to call it from Python. For what I was able to understand, Python can currently call Julia from IPython or through Julia's C API. It's not that easy.

Another option would be to create a code generator to translate Julia back to Python, and gradually switch the development to Julia. This isn't that complicated, as Julia is homoiconic, i.e. it can load and represent its own code in its data structures, and also provides macros to alter its own code.

By the way, in Julia there are two possible approaches:

• Types and multiple dispatch.
• Assumptions and pattern matching (minimal use of types).

The former is natively supported, the latter requires a pattern matching library (there are 3 under development in Julia).

I think that using types properly would allow to avoid pattern matching.

One important thing concerning multiple dispatch: remove all multiple inheritances! With multiple inheritance there can be matching ambiguities, such as

class A
  pass

class B
  pass

class C(A, B)
  pass

@dispatch(A)
def func(x):
  pass

@dispatch(B)
def func(x)
  pass

c = C()
func(c)  # which one of the two dispatched methods should it match???

[Cristóvão Duarte Sousa]

I'm not a Julia specialist but I think that for the last case one must do:

julia> g{T<:Integer}(v::Point{T}) = "Point... parametric type is Integer abstract class"

g (generic function with 1 method)

julia> g(p)

"Point... parametric type is Integer abstract class"

Anyway, I think that there must be a reason to allow the function definition the way you have put, or else it should not be allowed in the first place.

We may ask that on Julia users group...

[Cristóvão Duarte Sousa]

Sorry guys for the OT, but let me conclude the answer to F.B.:

g(v::Point{Integer}) = "Point... parametric type is Integer abstract class"

can be used to do

g(Point{Integer}(int64(1),int32(2)))

"Point... parametric type is Integer abstract class"

It seems to me, who knows nothing of computer science,  that in fact Julia has in fact a powerful multiple dispatch system (though only on types not on values).

[F. B.]

On Sunday, November 17, 2013 11:29:30 AM UTC+1, Cristóvão Duarte Sousa wrote:

It seems to me, who knows nothing of computer science,  that in fact Julia has in fact a powerful multiple dispatch system (though only on types not on values)

Multiple dispatch is based on the types of the function parameters. It's similar to function overloading in C++, except that function overloading creates a static link at compile time based on the types the variables are declared. Function overloading, in short, matches the highest superclass, while pattern matching matches the lowest subclass.

If you want more advanced matching systems than multiple dispatch, you should look at pattern matching. There's a Julia package under development to add pattern matching capabilities to Julia:

https://github.com/kmsquire/Match.jl

[F. B.]

By the way, what about if I try to sketch a @dispatch decorator?

It could be used experimentally in a subset of sympy, such as the tensor module and the quantum module.

Things one should care of are:

• behavior of @dispatch with other sympy decorators
• how to display a @dispatch decorated function in the doctests?
• difference between dispatching a function and a class method.
• keep it as close as possible to Julia's dispatch mechanism.
  

[Ondřej Čertík]

By all means, please go ahead!

We will discuss this once you post a PR, see how it looks, see if the
speed is enough
for tensors/quantum etc. And we can play with it.

I will try to use it experimentally with csympy as well.

Ondrej

[F. B.]

On Monday, November 18, 2013 10:22:09 PM UTC+1, Ondřej Čertík wrote:

By all means, please go ahead!

We will discuss this once you post a PR, see how it looks, see if the
speed is enough
for tensors/quantum etc. And we can play with it.

 Here we go, though it does not work yet: https://github.com/sympy/sympy/pull/2680