Thursday, July 26, 2012

Movimentum - Polynomial folding creates our second clip

We happily replace the old ConstantFoldingVisitor with the new PolynomialFoldingVisitor everywhere and expect that everything still works—or rather, works even better. After all, the new visitor does all of what the old one did (constants are polynomials!), and more. But alas, a handful of tests fails. For example, one test tells us that
Expected: <{EqualsZeroConstraint}0 = y+3*-y+y+6>
But was:  <{EqualsZeroConstraint}0 = ((-y+6))>
However, this was to be expected: Where we formerly got an expression tree, we now get a cleanly folded polynomial. Therefore, we change this test, together with the few others where the behavior has changed.

But the crucial test is our rotating bar. Here is, once more, its movimentum script (I have replaced the former x variable with len, because this value should be equal to the length of the bar):

    .config(8, 600, 400);
    WH : .bar C = [0,0] P = [-60,80];
    // length of bar: 100 - Pythagoras!

    @0
        WH.C     = [200,200];
        WH.P     = WH.C + [len,0].r(360*.t + 22.5);
    @2;

As before, the solver fails already for the first frame. But this time, we have much more manageable constraints:

{EqualsZeroConstraint}0 = ((-0.923879532511287*len))+-sqrt((-0.146446609406726*len^2+10000))
{EqualsZeroConstraint}0 = ((-0.38268343236509*len))+-sqrt((-0.853553390593274*len^2+10000))
{EqualsZeroConstraint}0 = -sqrt((-len^2+10000))

Which rules should we define to solve these constraints? As always, there's more than one way to skin the cat—not that I would actually consider that!—, but I opt for very "atomic" rules. Here is my suggestion:
  • 0 = –E can be rewritten to 0 = E;
  • 0 = sqrt(E) can also be rewritten to 0 = E;
  • Finally, 0 = P[V] should solve the equation for V. Right now, we will only consider linear and quadratic equations here, because they suffice to solve our current rotating bar problem.

Here is the code for the first rule:

    var e = new TypeMatchTemplate<IAbstractExpr>();
    new RuleAction<ScalarConstraintMatcher>("0=-E",
        new EqualsZeroConstraintTemplate(-e).GetMatchDelegate(),
        matcher => matcher != null,
        (currNode, matcher, matchedConstraint) =>
            new SolverNode(currNode.Constraints
                                    .Except(matchedConstraint)
                                    .Union(new[] {
                                      new EqualsZeroConstraint(matcher & e),
                                    }),
                            currNode.ClosedVariables,
                            currNode));

(Compared to the rewrites in the PolynomialFoldingVisitor, this code seems long-winded. I'll leave it like this right now, but maybe I'll extract a few methods when the code gets on my nerves).

The second rule requires almost the same code, only the template and the actual match different:

    ...
    var t = new UnaryExpressionTemplate(new PositiveSquareroot(), e);
    new RuleAction<ScalarConstraintMatcher>("0=sqrt(E)",
        new EqualsZeroConstraintTemplate(t).GetMatchDelegate(),
        ...

Finally, we must find zeros for polynomials—a textbook job. Here is the rule—the crucial missing part is the method FindZeros, which must return the collection of all value where the poynomial intersects the x axis:

    var p = new TypeMatchTemplate<IPolynomial>();
    new RuleAction<ScalarConstraintMatcher>("0=P[V]",
        new EqualsZeroConstraintTemplate(p).GetMatchDelegate(),
        matcher => matcher != null,
        (currNode, matcher, matchedConstraint) =>
            FindZeros(matcher & p).Select(zero =>
                currNode.CloseVariable(
                    (matcher & p).Var,
                    Polynomial.CreateConstant(zero),
                    matchedConstraint)
                 )
            );

Here is the full code for FindZeros:

    private static IEnumerable<double> FindZeros(IPolynomial polynomial) {
        switch (polynomial.Degree) {
            case 0:
                throw new ArgumentException("Cannot find zeros for constant");
            case 1: {
                // 0 = ax + b --> x = -b/a                
                var a = polynomial.Coefficient(1);
                var b = polynomial.Coefficient(0);
                return new[] { -b / a };
            }
            case 2: {
                // 0 = ax² + bx + c --> with D = b² - 4ac,
                // we get x1,2 = (-b+-sqrt D)/2a.
                // Attention: These formulas are not good
                // from a numerical standpoint! - see Wikipedia.
                double a = polynomial.Coefficient(2);
                double b = polynomial.Coefficient(1);
                double c = polynomial.Coefficient(0);
                double d = b * b - 4 * a * c;
                double sqrtD = d.Near(0) ? 0 : Math.Sqrt(d));
                return new[] { (-b + sqrtD) / (2 * a),
                               (-b - sqrtD) / (2 * a) };
            }
            default:
                throw new NotImplementedException( "Cannot find zeros for " + 
                                "polynomial of degree " + polynomial.Degree);
        }
    }

And this time, all the constraints are solved! We have another clip:




But what's that? The bar does not rotate! Rather, it flaps around helplessly only on the left side. What is wrong? To find this out, I wrote a new test that asserts each endpoint coordinate against a hand-computed value. Here is its core:

foreach (var f in frames) {
    IDictionary<string, IDictionary<string, ConstVector>>
        anchorLocations = f.SolveConstraints(10000, ref previousState);

    double alpha = f.FrameNo / 8.0 * 2 * Math.PI - Math.PI / 8;

    double expectedX = 200 + 100 * Math.Cos(alpha);
    double expectedY = 200 + 100 * Math.Sin(alpha);
    ConstVector whpLocation = anchorLocations["WH"]["P"];
    Assert.AreEqual(expectedX, whpLocation.X, 1e-2);
    Assert.AreEqual(expectedY, whpLocation.Y, 1e-2);
}

Here is its output before it fails:

++++ solution node ++++
SolverNode 1021 < 0<=C-1020 < 0<=C-1019 < 0<=C-1018 < 0<=C-1017 < 0=C-1015 < 0=P[V]...
  len:=-100
  WH.C.X:=200
  WH.C.Y:=200
  WH.C.Z:=0
  WH.P.X:=107,612046748871
  WH.P.Y:=161,731656763491
  WH.P.Z:=0
  ZERO:=0

  Expected: 292.38795325112869d +/- 0.01d
  But was:  107.61204674887132d

Looking at the variables in the last node, we see that len is minus 100—i.e., our solver found a solution where the bar has a negative length! That can only mean that the bar is oriented into the opposite direction of what we expect! We did not consider this when writing the script ... but the solution is easy: We add a constraint that restricts the length to positive values:

    ...;
    len       > 0; 
    ...;

And now the test that checks the coordinates happily runs green!

When we add this new constraint also to the clip-producing test, we finally get our second working clip out of Movimentum:



We are happy!

However, there is a usability issue lurking behind the problem with the negative length: How is the script writer ever going to find such missing constraints? It seems one needs a script debugging concept on the user level if it should be possible to create even moderately complex animation scripts. Mhm—one more item for the todo list.

Wednesday, July 25, 2012

Movimentum - Polynomial folding 3

5. Plus


Using "ASCII art," here is a table that shows how addition must work. Again, the goal is to merge polynomial subexpressions of a single variable at the left, and collect everything else in a single "rest" expression on the right. Be aware that E and F might also be polynomials, if their variables are different from V!

    //        | Q[V]+F              Q[V]            F
    // -------+------------------------------------------------
    // P[V]+E | {P[V]+Q[V]}+(E+F)   {P[V]+Q[V]}+E   P[V]+(E+F)
    //        |
    // P[V]   | {P[V]+Q[V]}+F       {P[V]+Q[V]}     P[V]+E
    //        |
    // E      | Q[V]+(E+F)          Q[V]+E          E+F

Why don't I consider entries with unary minus, or even square roots?
  • Square roots, of course, cannot be merged together by addition, so handling them separately will not create smaller, more compact expressions.
  • For unary minus, the argument is a little trickier: First, there cannot be a –P[V] subexpression anywhere—it would have been rewritten to {–P[V]}, i.e., a simple polynomial with inverted coefficients, on a lower level. The E's and F's could be expressions with a leading unary minus, but moving their minus upwards will not reduce the complexity of the expression in general. In specific cases, such a rewrite might get rid of expressions—e.g., if E happens to be –F or the like. But we do not optimize such "rare occurrences" here.
Polynomial folding is a heuristics. Like constant folding, it tries to standardize expressions quickly. However, if it oversees certain simplifications, this is not a big problem, as long as we do not rely on simplification guarantees.
Coding the table for the plus operator is easy—for example, here are the three entries for the first line of the table:

    new StandardExpressionRewrite("(P+E)+(Q+F)", _plusRewrites,
        (p + e) + (q + f),
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + ((m & e) + (m & f)));
    //                PO         AE         AE

    new StandardExpressionRewrite("(P+E)+Q", _plusRewrites,
        (p + e) + q,
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + (m & e));

    new StandardExpressionRewrite("(P+E)+F", _plusRewrites,
        (p + e) + f,
        m => (m & p) + ((m & e) + (m & f)));
    }

Below the first rewrite rule, I have marked the + operators with their "type":
  • PO is a polynomial operator, i.e., it will merge its two operands into a single polynomial.
  • AE, on the other hand, is an AbstractExpr operator, i.e., it will create an expression tree (a BinaryExpression).

The actual visiting code is similar to the square and unary minus operators:

    public IAbstractExpr Visit(Plus op, IAbstractExpr lhs,
                               IAbstractExpr rhs, Ignore p) {
        return Rewrite(new BinaryExpression(lhs, op, rhs),
                       "Cannot handle " + op, _plusRewrites);
    }


6. Times


My first attempt for multiplication was the same as for addition: Create a table that contains all combinations of "interesting expression structures," and explicitly define the result. From the square operator, we learn that we should handle square roots separately if they occur on both sides of the operator. In that first attempt, I identified 10 subexpression forms that I needed to distinguish: P+√R, P+–√R, P+√E, P+–√E, P+E, P+–E, P, C, –E, and general E. This gives us a table of 100 entries (or even more, because we have to distinguish (P+√R)·(Q+√R) from (P+√R)·(Q+√S), where all of P, Q, R, S are polynomials of the same variable). This is too much work—specifying, coding, and last, but not least, testing.

I thought a little bit and came up with the idea of stage-wise processing—it is what we do manually to "simplify" such expressions. For example, if we have (x + √4x) · (x - √x), we will
  • first multiply the two terms which gives us a sum of four terms;
  • then simplify each term;
  • and finally the whole sum.
The result from the first step can be a sum of up to four terms, i.e., an expression tree of depth at most 3. Various factors in the sums might need simplification—e.g. in our example, where √4x · √x will have to be merged to √4x2, which means that we have to run the tree through folding again at least three levels down from the current root. Instead of introducing a special mechanism for this, we run the visitor over the complete result—this will uselessly visit expressions below depth 3, but we will have to optimize away useless visits anyway somewhen in the future, so we will ignore this overhead right now.

The three steps above are only a rough plan. Actually, we need to consider a little bit more what "simplifications" we have to apply to the result of the first step. Here is my suggestion for a three-step algorithm.

In the first version of this posting, I bragged "similarly to all previous algorithms, I do not give any proof why this accomplishes what we want; it is "more or less obvious," ...". Well, in addition to being obvious, it was also wrong. The tables for Step III claimed, among others, that P√Q could be transformed to √(P2Q). However, this cannot be true: The latter is never negative (as we define sqrt to be the positive square root), but the former can of course be negative—for example, if the polynomial P is x – 1, and x happens to be zero. Here is a better table; but still no proof—I risk keeping this a mere software project, not a scientific undertaking, although I'd like to do all these proofs. Maybe I find time for them somewhen later ...

    // Step I: Distribute.
    //
    //          Q+F                 F
    //
    //  P+E     {PQ}+Q*E+P*F+E*F    P*F+E*F
    //
    //  E       Q*E+E*F             =
    //

    // Step II: Lift unary minus to top.
    //
    //          -F       F
    //
    //  -E      E*F      -(E*F)
    //
    //  E       -(E*F)   =
    //

    // Step III: Pull up square roots as far as possible.
    //
    //          √Q              √F              Q       D               F
    //
    //  √P      P=Q:  P         √(P*F)          =       D<0: -√{PD²}    =
    //          else: √{PQ}                             else: √{PD²}
    //
    //  √E      √(Q*E)          E=F:  E         =       D<0: -√({D²}*E) =
    //                          else: √(E*F)            else: √({D²}*E)
    //
    //  P       =               =               {P*Q}   {P*D}           =
    //                                                X
    //  C       C<0: -√({C²Q}   C<0: -√({C²}*F) {C*Q}   {C*D}           =  
    //          else: √({C²Q}   else: √({C²}*F)
    //
    //  E       =             =                 =       =               =
    //

An equals sign, in the tables above, simple means that no rewrite is done for the indicated product. Together, these three steps require 3 + 3 + 15 = 21 rules—a lot less than the 100 or more needed in my first attempt!
The four rewritings around the X can be handled by the same rule, which simply matches the product of two polynomials—also constants are polynomials!

Here is the code for Step I: The rules are almost like earlier ones—but in addition, there is a call to RecursivelyFold on the result:

    new StandardExpressionRewrite("(P+E)*(Q+F)", _timesRewrites,
        (p + e) * (q + f),
        HaveSameVar(p, q),
        m => RecursivelyFold((m & p) * (m & q)
                            + (m & q) * (m & e)
                            + (m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("(P+E)*F", _timesRewrites,
        (p + e) * f,
        m => RecursivelyFold((m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("E*(Q+F)", _timesRewrites,
        e * (q + f),
        m => RecursivelyFold((m & q) * (m & e)
                            + (m & e) * (m & f)));

RecursivelyFold is a single line method:

    private static readonly PolynomialFoldingVisitor
        _recursiveFolder = new PolynomialFoldingVisitor();

    private static IAbstractExpr RecursivelyFold(AbstractExpr expr) {
        return expr.Accept(_recursiveFolder);
    }

Here is an excerpt from Step III:

    new StandardExpressionRewrite("√P*√P", _timesRewrites,
        sqrtP * sqrtP,
        m => RecursivelyFold(m & p));
    new StandardExpressionRewrite("√P*√Q", _timesRewrites,
        sqrtP * sqrtQ,
        HaveSameVar(p, q),
        m => PosSqrtAndRecursivelyFold((m & p) * (m & q)));
    new StandardExpressionRewrite("√P*√F", _timesRewrites,
        sqrtP * sqrtF,
        m => PosSqrtAndRecursivelyFold((m & p) * (m & f)));
    // √P*Q --> =
    new StandardExpressionRewrite("√P*D,D<0", _timesRewrites,
        sqrtP * d, m => (m & d).Value < 0,
        m => -PosSqrtAndRecursivelyFold(
                    (m & d).Value * (m & d).Value * (m & p)
        ));

PosSqrtAndRecursivelyFold wraps the parameter in a PositiveSquareroot after recursively folding it:

    private static IAbstractExpr PosSqrtAndRecursivelyFold(AbstractExpr expr) {
        return new UnaryExpression(RecursivelyFold(expr), new PositiveSquareroot());
    }

If you look closely, you will notice that I put all the rules for multiplication into a single list, called _timesRewrites. This means that later stages, i.e. the RecursiveFold calls, will run over the rules for earlier stages again. This is actually useless and could be optimized away, but I guess that that would require even more complex algorithms than the ones I invent here. So I leave it as it is.

Of course, we also need an implementation of the Visit method: But it is similar to the one for the + operator.

Together with more than one hundred unit tests, this code rewrites binary expressions with a Times operator to our intended P+E format!

6. Divide


Division, finally, is very simple:
  • If the division is by a constant, we replace it with a multiplication of the reciprocal value.
  • Otherwise, we leave the expression alone (we do not try polynomial division or the like!):

    new StandardExpressionRewrite("E/C", _divideRewrites,
        e / c,
        m => RecursivelyFold(
              (m & e) * Polynomial.CreateConstant(1 / (m & c).Value)
            )
    );
    new StandardExpressionRewrite("E", _divideRewrites,
        e,
        m => m & e);

The corresponding Visit method is, of course, like the ones for plus and times.

And this concludes the visitor for rewriting polynomials.

Movimentum - Polynomial folding 2

In the previous posting, we created a long list of rules how to rewrite square expressions in the polynomial folder. Now, we have to code this. Obviously, we want to reuse our template matching machinery, which we already use to find constraints that we can expand. Roughly, our code should look like this:

    // Static definition of match template for P[V] + sqrt(E)
    var p = new TypeMatchTemplate<IPolynomial>();
    var e = new TypeMatchTemplate<IAbstractExpr>();
    var sqrtE = new UnaryExpressionTemplate(new PositiveSquareroot(), e);

    var template = p + sqrtE;

    ...

    // Code for rewrite rule using template
    var m = template.CreateMatcher().TryMatch(innerOfSquare);
    if (m != null) {
        var newExpr = m.Match(p).Times(m.Match(p)) +
                        m.Match(p).Times(2) * m.Match(sqrtE) +
                        m.Match(e);
        return newExpr;
    }

However, this code is much too long, and quite unreadable. First of all, the control logic has to vanish into a class so that we can at least write:

    // Static definition of rewrite rule using template
    var rule = new StandardExpressionRewrite(null, null,
                    p + sqrtE,
                    m => m.Match(p).Times(m.Match(p)) +
                         m.Match(p).Times(2)*m.Match(sqrtE) +
                         m.Match(e));

    ...

    // Execution of rule, somewhere else
    var result = rule.SuccessfulMatch(innerOfSquare);
    if (result != null) {
        return result;
    }

Another change easily done is the defininition of a binary operator on a matcher to get rid of all these Match calls. I choose the operator & as a replacement for Match. The rule now suddenly gets much shorter:

    new StandardExpressionRewrite(p + -sqrtE,
        m => ((m & p) * (m & p))
             + (-2 * (m & p) * (m & sqrtE)
             + (m & e)));

I say that this code is short and readable enough (even though the many parentheses make this somewhat Lisp-like ...).
To also get rid of the many m & repetitions (and one level of parentheses), one would have to introduce a static context, probably in a ThreadStatic field. I'll not do this.
Using this notation, we can now more or less easily write our rewrite rules for the square operator. Here are the first seven rules from the previous posting reformulated as C# code:

    new StandardExpressionRewrite("P", _squareRewrites,
        p,
        m => (m & p) * (m & p));

    new StandardExpressionRewrite("sqrt(E)", _squareRewrites,
        sqrtE,
        m => m & e);
    new StandardExpressionRewrite("-sqrt(E)", _squareRewrites,
        -sqrtE,
        m => m & e);
    new StandardExpressionRewrite("-E", _squareRewrites,
        -e,
        m => new UnaryExpression(m & e, new Square()));

    new StandardExpressionRewrite("P[V]+-sqrt(Q[V])", _squareRewrites,
        p + -sqrtQ, HaveSameVar(p, q),
        m => ((m & p) * (m & p) + (m & q))
             + -2 * (m & p) * (m & sqrtQ));
    new StandardExpressionRewrite("P[V]+sqrt(Q[V])", _squareRewrites,
        p + sqrtQ, HaveSameVar(p, q),
        m => ((m & p) * (m & p) + (m & q))
             + 2 * (m & p) * (m & sqrtQ));
    new StandardExpressionRewrite("P+-sqrt(E)", _squareRewrites,
        p + -sqrtE,
        m => ((m & p) * (m & p))
             + (-2 * (m & p) * (m & sqrtE) + (m & e)));

I added three more features:
  • HaveSameVar(p, q) returns true if both polynomials use the same variable, or if at least one of them is a constant. This is necessary for rules that require that two matched polynomials coincide on their variable, so that they can be merged into a single polynomial.
  • Each rule gets a string to identify it during debugging.
  • Also, the rules are stored in a List (here, _squareRewrites) that is passed to a general rewrite loop. This, finally, can be used to implement the visiting method for the square operator:

    private static readonly List<ExpressionRewrite> _squareRewrites =  
                                        new List<ExpressionRewrite>();

    public IAbstractExpr Visit(Square op, IAbstractExpr inner, Ignore p) {
        return Rewrite(inner, "Cannot handle square", _squareRewrites);
    }

    private static IAbstractExpr Rewrite(IAbstractExpr expr,  
             string cannotHandleMessage,  
             IEnumerable<ExpressionRewrite> rewrites) {
        foreach (var r in rewrites) {
            IAbstractExpr result = r.SuccessfulMatch(expr);
            if (result != null) {
                return result;
            }
        }
        throw new InvalidOperationException(cannotHandleMessage);
    }

Of course, it is utterly necessary to write test cases for these rewrites. I have created ten test cases for the square operator. Coincidentally, there are ten rewrite rules for the square operator, but I doubt that my test cases check all rules. I'll add a few more tests ... later.

But already a compile attempt shows us that there is still a piece missing: I happily use as yet undefined operators for polynomials and doubles, e.g.

   new StandardExpressionRewrite(...
      ...
      m => ... -2 * (m & p) ...);

After adding an operator that allows us to multiply a double with a polynomial, the rules compile—but one of the first tests is now red. Why? Because I also use plus and times operators to connect polynomials, e.g. here:

   new StandardExpressionRewrite(...
      ...
      m => ((m & p) * (m & p) + ...);

Because I defined these operators already for the AbstractExpr class, they now create a small expression tree from the two polynomials, instead of merging them into one polynomial! But easily enough, redefining the operators in class Polynomial solves the problem. Here is, for example, the unary minus for polynomials:

    public static Polynomial operator -(Polynomial p) {
        return CreatePolynomial(p.Var, p.Coefficients.Select(c => -c)).P;
    }

Plus and times are quite a bit longer, but not really complicated. And after their definition, all the square tests happily run green!

4. Unary minus


Folding a unary minus is now easy. We need just five rules:

    // -(P + -E) --> {-P} + E
    new StandardExpressionRewrite("P+-E", _unaryMinusRewrites,
        p + -e, m => -(m & p) + (m & e));
    // -(P + E) --> {-P} + -E
    new StandardExpressionRewrite("P+E", _unaryMinusRewrites,
        p + e, m => -(m & p) + -(m & e));
    // -(-P) should not happen; the inner -P must have been
    // rewritten to a polynomial with inverted coefficients!
    // -(P) --> {-P}
    new StandardExpressionRewrite("P", _unaryMinusRewrites,
        p, m => -(m & p));
    // -(-E) --> E
    new StandardExpressionRewrite("-E", _unaryMinusRewrites,
        -e, m => m & e);
    // -(E) --> -(E)
    new StandardExpressionRewrite("E", _unaryMinusRewrites,
        e, m => -(m & e));

and a simple visiting method:

    public IAbstractExpr Visit(UnaryMinus op, IAbstractExpr inner, Ignore p) {
        return Rewrite(inner,  
                       "Cannot handle unary minus",
                       _unaryMinusRewrites);
    }

Maybe it is instructive to disect the operators used in the rewrite result of second rule. There, the code is

        m => -(m & p) + -(m & e)
  • (m & p) is a polynomial, therefore, the unary minus to the left of it is the operator we defined above. It creates a new polynomial, where all coefficients are inverted.
  • (m & e), on the other hand, is an AbstractExpr. Therefore, the preceding unary minus is the expression tree building operator that we defined many postings ago.
  • Finally, the + operator is applied to a polynomial and an AbstractExpr. Therefore, the compiler cannot select the polynomial addition, but must use the AbstractExpr operator for building a BinaryExpression.
And thus, this short piece of code actually computes {–P} + –E.

Because we rewrite – –E to E in the fourth rule, sequences of unary minuses are reduced to a single or no unary minus. Here is a test case for this behavior:

    [Test]
    public void TestManyMinuses1() {
        IAbstractExpr a = P("x", 1, 2);
        IAbstractExpr b = a;
        var foldingVisitor = new PolynomialFoldingVisitor();
        for (int i = 0; i < 4; i++) {
            // Apply double unary minus
            b = -(-b.C);
            Assert.AreEqual(a, b.Accept(foldingVisitor, Ig.nore));
        }
    }

This is the reason why we could ignore – –E inside square roots.

And after a mere pair of postings, we are done with folding unary operators around polynomials. Next on the agenda: Binary operators—including the horrible multiplication of P+E's.

Movimentum - Polynomial folding 1

As explained in two previous postings, we want to standardize expressions as much as possible so that rules can easily match them for rewriting. One case in point are expressions that are a sum of some powers of a single variable, and an additional square root term. We would like to see to it that such expressions are either of the form
P[x] + sqrt(E[x])
or of the form
P[x] + –sqrt(E[x])

where P is a polynomial in x, and E is some arbitrary expression. Is it necessary that E is also of that form? Is seems not, because we will need to match E only after rewriting the whole expression, and then, additional standardization possibilities might kick in. However, our folding visitor will eagerly rewrite also subexpressions—it is just easier to program it like that.

The core part of the PolynomialFoldingVisitor is very simple. It follows the lines of the ConstantFoldingVisitor, but does not even bother to optimize for non-rewritten expressions (I hope very much that we can do optimization later). We start with the standard visits for constraints ...

    public AbstractConstraint Visit(EqualsZeroConstraint equalsZero, Ignore p) {
        return new EqualsZeroConstraint(equalsZero.Expr.Accept(this, p));
    }

    public AbstractConstraint Visit(MoreThanZeroConstraint moreThanZero, Ignore p) {
        return new MoreThanZeroConstraint(moreThanZero.Expr.Accept(this, p));
    }

    public AbstractConstraint Visit(AtLeastZeroConstraint atLeastZero, Ignore p) {
        return new AtLeastZeroConstraint(atLeastZero.Expr.Accept(this, p));
    }

... as well as all the expressions: Constants, variables, and polynomials are already standardized to the P + Ei form (there is no Ei in all these cases!), therefore the code is trivial; unary and binary expressions delegate to the respective operators:

    public IAbstractExpr Visit(IConstant constant, Ignore p) {
        return constant;
    }

    public IAbstractExpr Visit(INamedVariable namedVar, Ignore p) {
        return namedVar;
    }

    public IAbstractExpr Visit(IAnchorVariable anchorVar, Ignore p) {
        return anchorVar;
    }

    public IAbstractExpr Visit(UnaryExpression unaryExpr, Ignore p) {
        IAbstractExpr newInner = unaryExpr.Inner.Accept(this, p);
        return unaryExpr.Op.Accept(this, newInner, p);
    }

    public IAbstractExpr Visit(BinaryExpression binaryExpr, Ignore p) {
        IAbstractExpr newLhs = binaryExpr.Lhs.Accept(this, p);
        IAbstractExpr newRhs = binaryExpr.Rhs.Accept(this, p);

        return binaryExpr.Op.Accept(this, newLhs, newRhs, p);
    }

    public IAbstractExpr Visit(IGeneralPolynomial polynomial, Ignore p) {
        return polynomial;
    }

We "only" have to complete the code for the operators, and then we are done. Of course, the real meat is to be found there—therefore, I'll explain that code in increasing order of complexity:
  1. Formal square root is easiest, because it cannot even be evaluated.
  2. Unary operators that cannot be rewritten to polynomials are als easy: sine, cosine, and positive square root. We only rewrite them if their argument is a constant, so that we do also constant folding.
  3. Next, we look at square (which rewrites polynomials and square roots to new polynomials).
  4. Unary minus (which rewrites a polynomial to another, quite similar polynomial) comes next. Why do we look at unary minus only after square? With the square operator, is is possible to see all the problems of polynomial rewriting "in a nutshell"; unary minus is actually too simple and would lead us to a design that we'd have to refactor later.
  5. Then, we tackle the plus operator.
  6. The times operator is the most challenging of all.
  7. As an aftershot, we handle division.

1. Formal square root rewriting


This rewriting "does nothing"—it recreates the input expression:
    public IAbstractExpr Visit(FormalSquareroot op, IAbstractExpr inner, Ignore p) {
         return new UnaryExpression(inner, new FormalSquareroot());
    }

2. Sine, cosine, and positive square root rewriting


The code for these rewritings essentially does the following:
  • If the inner expression is a constant, evaluate the operator on the constant and return the new constant.
  • Otherwise, wrap the inner expression in the operator.
The FoldIfConstant method gets inner as a parameter, but also a delegate that can do the evaluation if inner is a constant:

    public IAbstractExpr Visit(PositiveSquareroot op, IAbstractExpr inner, Ignore p) {
        return FoldIfConstant(op, inner, x => x.Near(0) ? 0 : Math.Sqrt(x));
    }
    public IAbstractExpr Visit(Sin op, IAbstractExpr inner, Ignore p) {
        return FoldIfConstant(op, inner, x => Math.Sin(x / 180 * Math.PI));
    }
    public IAbstractExpr Visit(Cos op, IAbstractExpr inner, Ignore p) {
        return FoldIfConstant(op, inner, x => Math.Cos(x / 180 * Math.PI));
    }
    private static IAbstractExpr FoldIfConstant(UnaryOperator op,
             IAbstractExpr inner, Func<double, double> eval) {
        IConstant innerAsConst = inner as IConstant;
        return innerAsConst == null
                    ? (IAbstractExpr)new UnaryExpression(inner, op)
                    : Polynomial.CreateConstant(eval(innerAsConst.Value));
    }

3. Square rewriting, part 1


Now we come to the crucial code: "Hard rewriting." We look at the square operator first, because its interaction with square roots will show more clearly what rewritings we want to do. We know that visitors of subexpressions will already have reformulated the inner expression to one of the following three forms:
  • P[V] + E
  • P[V]
  • E
Here, P[V] is a polynomial in the single variable V, and E is not—i.e., it is either not a polynomial at all, or a polynomial in some other variable. Two remarks are probably necessary:
  • P[V] might also be a constant, in which case we assume its variable is "some arbitrary or dummy variable for constants."
  • Because E may be a polynomial in some other variable, i.e. Q[W], it turns out that our standardization is not a normalization, as the following example shows: The expression 2x + 3y may be represented in two different ways—either the leading polynomial is 2x, and the "rest" E is 3y, or vice versa. One could probably repair this by requiring that the alphabetically lowest variable must be used to create the polynomial—but right now, I am not interested in this problem.
Now, we want to apply a square operator to such an expression, and again get an expression of on of these forms in as many cases as possible. Let us start with the simplest case, namely the second form P[V]: Of course, the square of a polynomial is again a polynomial in the same variable! We can sketch this rewriting as follows:
(P[V])2 → {P[V]2}
The curly braces mean that the expression inside actually has to be simplified to a polynomial—we will therefore need a method that squares a polynomial! But before we write this, let us look at the other possible inner expressions.

The next interesting case is the lonely E. One could believe that no rewritings are possible here except the trivial reconstruction of the input expression, but this is not true. Consider the case where E is of the form sqrt(E), where E is some other expression (sorry for using two different E's here—E and E; but in the code, I use E and e for the inner expressions, which I did not want to change).

Ok: What is the square of sqrt(E)? Of course, it is simply E, and we should capture this in a specific rewriting! And now it pays off that the visitor also folds inner expressions—we can now be sure that E is already standardized, and therefore we can happily sketch the rewriting as
(sqrt(E))2 → E
Next inside operator: If there is a unary minus operator inside the square, we can get rid of that operator, and we could be tempted to write
(–E)2 → E2
But this is not correct! Assume that E is sqrt(x). Then, (–sqrt(x))2 should not be rewritten to the unary expression (sqrt(x))2, but to the polynomial x. There are two ways to solve this problem:
  • One possibility is to recursively run it again through the standardization! We only have to take care that we do not produce an infinite recursion. But as long as our standardization is a "folding," i.e., a simplification, there is some metrics on the expression that gets strictly smaller with each rewriting—e.g., the number of operators in the whole expression. And when there is a lower limit on this number (zero, in the example metrics), we know that simplification must end somewhen. We will use this for the most difficult of all operators, namely multiplication.
  • Here, we have a simpler solution: We add a more specific rewriting for square roots:
(–sqrt(E))2 → E
Do we also have to consider "minus inside minus," i.e., (– – E)2? No, because – – E must already have been simplified by a lower visitor; we'll look to that in the section about the unary minus operator.

Finally, we look at the third, most general form of an inner expression: P[V] + E. If we square it, we get of course P[V]2 + 2 P[V] E + E2. Which parts of this can be simplified? Of course, P[V]2 and 2 P[V] are both polynomials.

But we must be more careful: Suppose that P[V] is x, and E is sqrt(x). Then, (x + sqrt(x))2 = x2 + 2x sqrt(x) + x, which we would like to have simplifed to the polynomial x2 + x plus the middle expression. And there might be other cases where also the middle term has to be folded into the result somehow ... What do we do? Again, we have the two possibilities:
  • Recursive standardization.
  • Explicitly add rules for square roots that handle special cases.
Again, I opt for special casing here. The cases are, from most specific to least specific:
  • inner expression is P[V] + –sqrt(Q[V]) → square is {P[V]2 + Q[V]} + {–2P[V]}sqrt(Q[V]).
  • inner expression is P[V] + sqrt(Q[V]) → square is {P[V]2 + Q[V]} + {2P[V]}sqrt(Q[V]).
  • inner expression is P[V] + –sqrt(E) → square is {P[V]2} + ({–2P[V]}sqrt(E) + E). We parenthesize the last two terms so that the complete result is of the form P + E.
  • inner expression is P[V] + sqrt(E) → square is {P[V]2} + ({2P[V]}sqrt(E) + E). Again, the last two terms are grouped together.
  • inner expression is P[V] + –E → square is {P[V]2} + ({–2P[V]}E + (E)2).
  • inner expression is P[V] + E → square is {P[V]2} + ({2P[V]}E + (E)2).

We now have a specification for folding all expressions with square operator. But how to code it?

Friday, July 20, 2012

Movimentum - Operators on interfaces do not work: A quick fix

The major change from the last posting was the introduction of interfaces for (most) expression types. This made it necessary to replace constructor calls for constants and variables with calls to factory methods whose return type is one of the new interfaces. The change from classes to interfaces ripples through the whole software, but thanks to ReSharper, it is quickly done.

Everything looks nice, until an intermediate compile suddenly turns up an bad problem: Operators cannot be defined solely on interfaces!

The problem is visible in many test cases, but also in some rewrite rules, where we could write

     new Constant(1) + new Constant(2);

until now, whereas

    Polynomial.CreateConstant(1) + Polynomial.CreateConstant(2)

is no longer possible: The + operator is defined for parameters of class AbstractExpr, but the results of the CreateConstant calls are IAbstractExprs! When we try to replace the parameters of operator+ with IAbstractExpr, the compiler tells us in no uncertain terms that at elast one of the parameters must be of the surrounding class.

What can we do?

Here is my quick fix: We add a property C to IAbstractExpr which returns the same object, but now as an AbstractExpr. The implementation, in class AbstractExpr, is trivial:

    public AbstractExpr C { get { return this; } }

The 1 + 2 example can now be written as follows, which happily compiles:

    Polynomial.CreateConstant(1).C + Polynomial.CreateConstant(2).C

To get rid of about half of these .Cs, C# allows us to pass in one interface parameter:

    public static AbstractExpr operator +(AbstractExpr lhs, IAbstractExpr rhs) {
        ...
    }

The example from above can now be written

    Polynomial.CreateConstant(1).C + Polynomial.CreateConstant(2)

As a side note, this is an example where interfaces are not introduced to be implemented by new classes. The design of the solver model is "closed" in the sense that no effort goes into considerations how to write new classes that implement one of the interfaces outside the controlled design of the solver engine. Therefore, there is no problem with the C property: This is not an abstract request to convert any sort of IAbstractExpr to an AbstractExpr; rather, it is simply the technically helpful statement that every IAbstractExpr object is an AbstractExpr object, anyway, and that we want to expose this simple fact because of C#'s operator limitation.
After this side note—that could be expanded to a whole article on how interfaces can be used; and what "abstraction" can mean—, we return to the mundane job of liberally sprinkling our code with .Cs, until it compiles again.

Tuesday, July 17, 2012

Movimentum - Creating polynomials

The whole goal of the solver model is the easy application of rewriting rules. And therefore, the main design force is standardization of expression: If a certain mathematical expression is always represented by the same expression tree in the solver model, it can reliably be matched against a rule. Consequences of this design rationale are:
  • The introduction of the solver model, in the first place: Expressions are now always represented as expression trees of scalar variables.
  • Constant folding: Constants are now always represented in the same way, i.e., they are normalized.
  • Polynomials: Chains of plus, times, and square operators containing only a single variable will always be represented by the same expression tree, namely a Polynomial object (we are not there yet).
However, in our current design, two types of polynomials can be represented by different expression trees: Constants, and variables.

A constant, right now, can either be represented by a Constant object or by a Polynomial object with degree zero, i.e., with only a single coefficient. Similarly, a variable can be represented either by an object whose type is assignable to Variable; or by a Polynomial object with coefficients [1, 0]. Obviously, this is not what we want.

What we rather would like to have, is twofold:
  • Rules matching a polynomial also handle constants and variables without any further ado.
  • Rules can reliably match constants and variables.
The solution for the first item is simple: Derive classes Constant and Variable from class Polynomial. However, if we allow Polynomial to remain a concrete class, this creates a subtle problem with visitors: We suddenly need visit methods that can override each other. In virtually all descriptions of the visitor pattern, this topic is ignored, because it is assumed—to my knowledge, always implicitly—that only leaf classes of a visited inheritance tree are concrete classes. This is not required by the visitor pattern, but it is good design practice, and we will follow it.

Therefore, we design our class structure for polynomials as follows:


What about the second item—reliably matching constants and variables?

This requires that constants are always represented by objects of type Constant; and likewise for variables. As long as the constructor for GeneralPolynomial is public, this cannot be enforced (or it could only be enforced by a runtime check, which would put the burden of using the correct constructors on the calling site).

The correct solution is to define a factory method, and to make the constructor of GeneralPolynomial private. However, the factory mathod must be able to access this constructor, and hence we must either put it into the GeneralPolynomial class, or we must make this class a private inner class of Polynomial. The first approach works, but it is (in my opinion, "strangely") asymmetric: The factory method can only access the private constructor of GeneralPolynomial, but must use the public constructors of Constant and Variable. I opt for a different design:
  • Instead of using classes for polynomials, we use interfaces everywhere.
  • All three concrete polynomial classes become private classes inside class Polynomial.
Here is a diagram showing this setup. The interfaces are colored yellow, the four classes comprising the polynomial inheritance hierarchy are green. The lines with the "plus circles" indicate "nesting," i.e., that the lower classes are nested inside the upper class (from which they are also derived—but these two relationships are orthogonal to each other!):


Why don't I introduce interfaces also for the remaining expression classes, namely UnaryExpression and BinaryExpression? Because ... mhm ... this would be YAGNI: There is, right now, no reason at all to change their design. Lame excuse?

Anyway, now we can write the factory method for polynomials:

    public static IPolynomial CreatePolynomial(IVariable variable,
                                    IEnumerable<double> coefficients) {
        // Skip all leading zero coefficients
        int zeroPrefixLength = 0;
        foreach (var c in coefficients) {
            if (!c.Near(0)) {
                break;
            }
            zeroPrefixLength++;
        }
        double[] normalizedCoefficients =
            coefficients.Skip(zeroPrefixLength).ToArray();

        // Create intermediate coefficient for zero polynomial.
        if (normalizedCoefficients.Length == 0) {
            normalizedCoefficients = new[] { 0.0 };
        }

        // Create result based on degree and coefficients.
        var deg = normalizedCoefficients.Length - 1;
        if (deg == 0) {
            return new Constant(normalizedCoefficients[0]);
        } else if (deg == 1
                   && normalizedCoefficients[0].Near(1)
                   && normalizedCoefficients[1].Near(0)) {
            return variable;
        } else {
            return new GeneralPolynomial(variable, normalizedCoefficients);
        }
    }

Because the constructors and variables are no longer visible, we must provide factory methods for variables (otherwise, we could never create a polynomial!). For simplicity, there is also a factory method for constants, so we do not have to create a polynomial with a throw-away variable and a single coefficient. By the way, these methods could, in the future, cache often-used variables and constants for performance reasons:

    public static IAnchorVariable CreateAnchorVariable(Anchor anchor,
                                    Anchor.Coordinate coordinate) {
        return new AnchorVariable(anchor, coordinate);
    }

    public static INamedVariable CreateNamedVariable(string name) {
        return new NamedVariable(name);
    }

    public static IConstant CreateConstant(double value) {
        return new Constant(value);
    }

Of course, we must also modify the interface for visiting expressions—instead of the three polynomial classes, it now uses interfaces. The same replacement must happen in all implementing classes. ReSharper can do this with a few mouse-clicks:

public interface ISolverModelExprVisitor<in TParameter, out TResult> {
    TResult Visit(IConstant constant, TParameter p);
    TResult Visit(INamedVariable namedVariable, TParameter p);
    TResult Visit(IAnchorVariable anchorVariable, TParameter p);
    TResult Visit(UnaryExpression unaryExpression, TParameter p);
    TResult Visit(BinaryExpression binaryExpression, TParameter p);
    TResult Visit(RangeExpr rangeExpr, TParameter p);
    TResult Visit(IGeneralPolynomial polynomial, TParameter p);
}

This concludes our revised solver model, which now contains polynomials as first class citizens. The next step is to rewrite all expressions to the form P + Ei. I call this "polynomial folding"—it is, so to speak, the "big cousin" of constant folding.