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.