Symbolic higher order binary polynomials

[mhlr]

Application
Currently dimod has a BinaryPolynomial class that can be constructed from higher order mappings like {(1,2,3): 1}
which can be passed to make_quadratic to create BQMs. We also have symbolic experssions for constructing BQMs.
However it is not possible to construct higher order BinomialPolynomials symbolically

dimod.Binary(1)*dimod.Binary(2)*dimod.Binary(3)

throws

TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_1301/2702588800.py in <module>
----> 1 dimod.Binary(1)*dimod.Binary(2)*dimod.Binary(3)

~/anaconda3/lib/python3.8/site-packages/dimod/binary/binary_quadratic_model.py in __mul__(self, other)
    258         if isinstance(other, BinaryQuadraticModel):
    259             if not (self.is_linear() and other.is_linear()):
--> 260                 raise TypeError(
    261                     "cannot multiply BQMs with interactions")
    262             elif other.num_variables and other.vartype != self.vartype:

TypeError: cannot multiply BQMs with interactions

Proposed Solution
Extend symbolic expressions to BQMs with interactions and BinaryPolynomials:

BQM * BQM --> BinaryPolynomial, if either argument contains interactions
BQM BinaryPolynomial --> BinaryPolynomial
BinaryPolynomial BQM --> BinaryPolynomial
BinaryPolynomial BinaryPolynomial --> BinaryPolynomial

Also support non-negative integer powers of BQMs and BinaryPolynomials (see #864)

[arcondello]

I love this suggestion in principal but there are a few things that we need to do before we should allow this:

1. We need a performant implementation of BinaryPolynomial. Right not it's pure Python and not all that performant an implementation at that. In practice that means that we would see a big performance cliff as soon as the degree rises above two.
2. You cannot use a graph structure to "natively" encode a binary polynomial, so a lot of the core complexity assumptions start to fall apart. There are some encoding you can do (like a bipartite graph from terms to variables) but we couldn't re-use a lot of the existing infrastructure.
3. We need the implementation of BinaryPolynomial to be a closer match to QuadraticModel and BinaryQuadraticModel in terms of API. QM and BQM are deliberately similar in their methods/access, but BinaryPolynomial less so.
4. We probably need a Polynomial that supports integer variables while we're at it.
5. We need a solver for (binary) polynomials.

[mhlr]

Re 1) there is a lot of room for optimization within pure python.
I can start incrementally improving performance in the mean time
There are a lot of Python loops that can be optimised away.
The __eq__ method looks like particularly low hanging fruit.

[mhlr]

For 5) do you mean something beyond HigherOrderComposite?

[arcondello]

For those who want to do this in the mean time, you can actually use sympy to get something like the behaviour you want with

import itertools

import dimod

from sympy import Poly


def sympy_to_binary_polynomial(f, *, vartype):
    variables = tuple(f.free_symbols)
    f = Poly(f, variables)

    poly = dict()
    for powers, coeff in f.as_dict().items():
        term = tuple(itertools.chain(*(itertools.repeat(v, count) for v, count in zip(variables, powers))))
        poly[term] = coeff

    return dimod.BinaryPolynomial(poly, vartype=vartype)


if __name__ == '__main__':
    from sympy import symbols

    x, y, z = symbols('x y z')
    expr = x*(x - 1)*(y + 1) + z + 3
    print(sympy_to_binary_polynomial(expr, vartype="BINARY"))

I have not rigorously tested this though.