Movimentum - Backsubstitutions

Here is our next example:

The rotating bar

The idea is that a bar (a line) rotates about one of its endpoints—that's it.
I thought a little about the details and came up with the following:

• (a) The bar is defined in some "weird direction"—not horizontally, not vertically. The reason for this is that I want to see that the initial location is still as the time line defines it.
• (b) The length of the bar should be some nice number. Let's take 100.
• (c) There are eight images per full round. The first one shows the bar at 22.5° to the horizontal, each subsequent one rotates it by 360°/8 = 45°. The reason for this is that these positions can be verified visually easily, because the endpoints form an equilateral octagon. Moreover, there are solutions with equal x coordinates, but different y coordinates for the rotating endpoint. This could show problems with selecting the right solution—I expect that the bar would somehow "jump" up or down if there is such a problem.
• (d) For simplicity, the bar rotates once per second. Thus, we get images at .t=0, .t=0.125, .t=0.250, .t=0.375, .t=0.500 etc.

Here is the source for this test case.

    [Test]
    public void TestRotatingBar() {
        const string s =
            @".config(8, 600, 400);
            WH : .bar C = [0,0] P = [-60,80];

            @0  WH.C = [200,200];
                WH.P = WH.C + [10*x,0].r(360*.t + 22.5);
            @2";
        Script script = Program.Parse(s);

        string prefix = MkTestDir();
        Program.Interpret(script, prefix);
    }

When we run the test, it immediately returns the exception "No solution found for frame 1"—we are back in the business of working on our rewriting rules.

The debugging output of the test shows us eight dead nodes with very similar unmatched constraints. Here is the first such node:

SolverNode 1015 < root*1013 < root*1011 < root*1010 < 0=C-1009 < 0=V+C...
  ! {EqualsZeroConstraint}0 = 200+-(WH.P.X+sqrt(10000+-((200+-WH.P.Y)²+0)))
  ! {EqualsZeroConstraint}0 = 200+-(WH.P.Y+sqrt(10000+-((200+-WH.P.X)²+0)))
  ! {EqualsZeroConstraint}0 = 0+-(0+sqrt(10000+-((200+-WH.P.X)²+(200+-WH.P.Y)²)))
  ! {EqualsZeroConstraint}0 = WH.P.X+-(200+10*x*0.923879532511287+0)
  ! {EqualsZeroConstraint}0 = WH.P.Y+-(200+10*x*0.38268343236509+0)
  ! {AtLeastZeroConstraint}0 <= 10000+-((200+-WH.P.Y)²+0)
  ! {AtLeastZeroConstraint}0 <= 10000+-((200+-WH.P.X)²+0)
  ! {AtLeastZeroConstraint}0 <= 10000+-((200+-WH.P.X)²+(200+-WH.P.Y)²)
  WH.C.X = 200
  WH.C.Y = 200
  WH.C.Z = 0
  WH.P.Z = 0
  ZERO = 0

What can we do about these constraints to get nearer to a solution? The first three constraints look "complicated"—we'll handle them later. But the next two constraints are essentially assignments: A constraint of the form 0=V+E allows us to replace V with (–E) everywhere. And if V does not occur in E, we got rid of that variable entirely!
Of course, this does not give us a direct solution for V (as did 0=V+C). But when we later have solutions for the variables in E, their values can be "back-substituted" into E, which—after constant folding—will give us a solution for V. When we find a 0=V+E constraint, we therefore need to do two things:

• We have to replace all V's in all constraints by –E.
• We must remember somewhere that V is –E.

Handling 0=V+E

The first item is easy—we have the method SubstituteVariable for it.
For the second item, we need a data structure where we store that knowledge. There are at least the following three possibilities:

• a) We store it as a (special type of) constraint; nice effect: Rewrites will still work on it, especially replacements of variables with constants.
• b) We store that knowledge in a totally new data structure.
• c) We store it as a (new sort of) "VariableKnowledge."

The nice effect of "store it as a constraint" idea a) is that variable substitutions by the 0=V+C rule will then automatically be performed. Together with constant folding, this will eventually lead to a constant in that constraint. The problem with the same idea is, of course, that we somehow must remove that constraint from the constraints under consideration—otherwise, we will go into an endless loop that finds that same constraint over and over. Ideas to solve this problem, all with their own problems:

• a.1) Add a flag to the constraint that says "don't select." However, this has the new problem that, when the constraint's expression has been rewritten to a constant, it will not be selected even though it matches the 0=V+C template.
• a.2) Remember by which rule that constraint was selected, and suppress selection in the corresponding rule action. This requires very special action, and possibly a mutable flag at the constraint.
• a.3) Add a new, fourth constraint class derived from ScalarConstraint. The problem with this approach is that all the visitors have to be extended.
• a.4) Add a new constraint class derived from EqualsZeroConstraint and suppress its selection in the same rule by a specific check. The problem with this is that it violates the Liskov substitution principle: Objects of the subclass do not behave like objects of the superclass (this is an effect often seen when there are non-abstract superclasses in an inheritance hierarchy). As a concrete risk, rewriting visitors that have to create a new EqualsZeroConstraint in some cases will ignore the information transported by the special class, and recreate them as EqualsZeroConstraint—and the special information is lost! So this is a very risky approach.
• a.5) Extend the class hierarchy: We make EqualsZeroConstraint abstract (which requires a modification to all creation locations), and add two new derived classes. If we update all visitors to use only concrete classes, this creates as much work as option c). If we let the visitors work on the abstract class, we will violate the principle that our visitors only visit leaf classes (the reason for this principle, which I forgot to state when I introduced the visitors here, is that it must always be clear which method is called for visiting so that one can override the right method. There are other, more complicated idioms for this, but I want to keep the design concepts of my solver as simple as possible).

Mhm. Actually, the "nice effect" that variables are rewritten can also be viewed as a problem: All the rules still run over the 0=V+E constraints, even though they already have been dealt with. The whole idea looks more and more suspicious.

How much work would option b) be?

• We could use a dictionary to store each variable with its expression in the SolverNode—simple.
• We need to pass in this new data structure into the SolverNode constructor.
• Somewhen after a variable is "solved," we need to rewrite the stored expressions. If this rewriting produces a constant, we add the variable-value pair to the existing variable knowledge dictionary. Option one: We do this in the central method Remember...Variable() that captures such solvings. Option two: We do it at the end of the whole solving process, i.e., only when we have found a solution node!

But that should be it. It seems that solution b) is much less problematic than a). Let us also think about option c) for a moment:

• We reuse the existing dictionary that right now stores variables with their values. However, the current homogeneous data structure becomes hetereogeneous, i.e., it contains objects of different types: Mappings from variables to doubles; and mappings from variables to expressions.
• No new data structures are passed around.
• The process for resolving the backsubstitutions is the same as for b).

Actually, I tried b) for an hour or so (and even checked a version into Github). But I had to pass in more and more emtpy dictionaries which I did not want to see. One remedy would have been to add a convenience constructor that handles this—but when writing new code, I am very wary of convenience methods that pass in default values for some parameters (and therefore also of default paramtetes): More than once this "convenience" has backfired, because I forgot to think about a critical parameter. Anyway, I ripped out that parameter again and decided to go with option c).

Why formal square roots must not be substituted

There is one last item we must be careful about: If the expression in 0=V+E contains a formal square root, we must not create a backsubstition. Why? The problem is that due to the substituion of V, the formal square root could emerge multiple times in the resulting constraints. Later, these square roots will be rewritten to positive and negative roots in different node branches—but in some of these branches, the roots will be rewritten differently, although they are part of the same variable definition! And this may lead to invalid "solutions."

Maybe an example shows the problem: Here are three equations:

a = root(b)
b = a
c = –a

Let us first follow the method to rewrite the formal square root beforehand. We end up with two branches:

a = sqrt(b)
b = a
c = –a

or

a = –sqrt(b)
b = a
c = –a

We can solve the first branch by substituting sqrt(b) for a everywhere:

b = sqrt(b)
c = –sqrt(b)

from which we get the solutions

b = 0, c = 0, a = 0

and

b = 1, c = –1, a = –1

The second branch, on the other hand, yields

b = –sqrt(b)
c = sqrt(b)

with the single solution b = 0, c = 0, a = 0.

Now, let us look what happens if we immediately rewrite a with root(b) at the beginning. We get

b = root(b)
c = –root(b)

and as one of the four branches after rewriting the two roots, we get

b = sqrt(b)
c = –(–sqrt(b))

which has as a solution where b = 1, c = 1, and |a| = 1! Obviously, this contradicts the two original equations that say

b = a
c = –a

So, we must only rewrite 0=V+E constraints where E does not contain a formal square root. More generally, it must be possible to evaluate E uniquely. This is the second time we need to distinguish uniquely evaluatable expressions from others—the first occasion was in the constant folding visitor. We are not yet at the threshold of the "Rule of Three," and in order to not get sidetracked, we keep to our goal of rewriting 0=V+E—but next time we single out formal square roots, we will introduce some abstraction for them!

Two unit tests for backsubstitutions

Back to 0=V+E. We could simply implement it and try out our rotating bar. However, it seems that this rewriting is not as easy as it seemed. So we should write at least one focused unit test that tests the new data structures. Here is a first attempt—comments in the code explain the steps:

     [Test]
     public void Test0VERewritingOneStep() {
        // The 0=V+E constraint: 0 = x + y.
        EqualsZeroConstraint ve = new EqualsZeroConstraint(NV("x") + NV("y"));
        // Another constraint that contains x: 0 = y + (3x + y + 6)
        EqualsZeroConstraint other = new EqualsZeroConstraint(
            NV("y") + (new Constant(3) * NV("x") + NV("y") + new Constant(6)));
  
        var current = new SolverNode(new[] { other, ve }, null);
        // Do a step
        SolverNode solutionOrNull;
        IEnumerable<SolverNode> result =
            SolverNode.SolverStep(new[] { current }, out solutionOrNull);

        // Afterwards, we have a new node, but not yet a solution.
        Assert.AreEqual(1, result.Count());
        Assert.IsNull(solutionOrNull);
        SolverNode resultNode = result.ElementAt(0);
        Assert.AreEqual(1, resultNode.Constraints.Count());

        // The single constraint is other, with x replaced with -y.
        Assert.AreEqual(new EqualsZeroConstraint(
             NV("y") + (new Constant(3) * -NV("y") + NV("y") + new Constant(6))),
             resultNode.Constraints.ElementAt(0));

        // Moreover, we have a backsubstitution.
        Assert.AreEqual(-NV("y"), resultNode.Backsubstitution(NV("x")));
    }

Because we have not yet declared the backsubstitutions, this test does not even compile. We will add the necessary code so that it compiles and runs in a moment. However, this test obviously falls short of completely testing backsubstitutions:

• It does not check that backsubstitutions happen at all.
• It does not test that backsubstitutions that lead to solved variables also solve other dependent variables—i.e., test the "backsubstitution ripple effect."

Even though our first test does not yet run, I'd like to keep on the "tester's hat" and write another test that checks these additional requirements. Unfortunately, it is somewhat long because it has to produce that ripple effect:

    [Test]
    public void Test0VERewritingTwoSteps() {
        EqualsZeroConstraint veConstraint1 =
            new EqualsZeroConstraint(NV("z") + (NV("x") + new Constant(4)));
        EqualsZeroConstraint veConstraint2 =
            new EqualsZeroConstraint(NV("x") + NV("y"));
        EqualsZeroConstraint constraint3 = new EqualsZeroConstraint(
            NV("y") + (new Constant(5) * NV("x") + NV("y") + NV("z")));

        var current = new SolverNode(new[] {
            veConstraint1, veConstraint2, constraint3}, null);

        // Do two steps that rewrite ve1 and ve2.
        SolverNode solutionOrNull;
        IEnumerable<SolverNode> expanded1 =
            SolverNode.SolverStep(new[] { current }, out solutionOrNull);
        Assert.AreEqual(1, expanded1.Count());
        Assert.IsNull(solutionOrNull);

        // After this first step, we should have the backsubstitution
        //   z &rarr; -(x+4)
        // and the constraints
        //   0 = x + y
        //   0 = y + (5x + y + -(x+4))
        IEnumerable<SolverNode> expanded2 =
            SolverNode.SolverStep(new[] {expanded1.First()}, out solutionOrNull);
        Assert.AreEqual(1, expanded2.Count());
        Assert.IsNull(solutionOrNull);

        // After this step, we should have the two backsubstitutions
        //   z &rarr; -(x+4)
        //   x &rarr; -y
        // and the single constraints
        //   0 = y + (5(-y) + y + -((-y)+4))
        SolverNode expanded2Node = expanded2.ElementAt(0);
        Assert.AreEqual(1, expanded2Node.Constraints.Count());

        // Now, we must solve the last single constraint. It can be
        // simplified to
        //   0 = y - 5y + y + y - 4
        // or
        //   0 = -2y - 4
        // which gives us
        //   y = -2
        // We simulate this solving by creating a new node with the
        // correct value for y:
        SolverNode result =
            expanded2Node.RememberAndSubstituteVariable(NV("y"), -2);

        // This now should put all solutions into the variable knowledge:
        Assert.AreEqual(-2, result.VariableInRangeKnowledges[NV("y")].Value);
        Assert.AreEqual(2, result.VariableInRangeKnowledges[NV("x")].Value);
        Assert.AreEqual(-6, result.VariableInRangeKnowledges[NV("z")].Value);
    }

Also this test only scratches the surface of backsubstitution testing: The backsubstituting algorithm will certainly contain loops (to handle more than one backsubstitution and/or more than one variable that has a value), and therefore the test should have data that will force such loops to run at least two rounds. However, in the blog, I do not show the even longer tests for such cases.

Implementing backsubstitutions—the rule

Here is the crucial change that should make the first test green: We add a rule that recognizes 0=V+E if E has no formal square root; and does the substitution V ? –E:

{
    // 5. Substitute variable with expression that does not
    //    contain the variable and has no formal square root.
    var v = new TypeMatchTemplate<Variable>();
    var e = new TypeMatchTemplate<AbstractExpr>();
    new RuleAction<ScalarConstraintMatcher>("V->-E",
        constraint => {
            // Check for v+e match.
            ScalarConstraintMatcher m =
                new ScalarConstraintMatcher(
                    new EqualsZeroConstraintTemplate(v + e))
                .TryMatch(constraint);
            if (m == null) {
                return null;
            }
            // Check that no formal square root exists.
            FindFormalSquarerootVisitor f =
                constraint.Expr.Accept(
                    new FindFormalSquarerootVisitor(), Ig.nore);
            if (f.SomeFormalSquareroot != null) {
                return null;
            }
            // Check that v is not in e.
            Dictionary<Variable, VariableDegree> degrees =
                m.Match(e).Accept(new VariableDegreeVisitor(), Ig.nore);
            Variable variable = m.Match(v);
            if (degrees.ContainsKey(variable)
                && degrees[variable] != VariableDegree.Zero) {
                return null;
            }
            return m;
        },
        m => m != null,
        (currNode, matcher, matchedConstraint) =>
            currNode.SubstituteVariable(matcher.Match(v),
                                       -matcher.Match(e))
    );
}

The first test should now fail only at its last line, as we do not yet fill the backsubstitutions. So let us run it! Yes, it is red—but the stack trace shows us that it fails much earlier! What went wrong?

Debugging shows us that there are still two constraints in the expanded node, instead of one—and now we remember: When we find a 0=V+C constraint and afterwards do a substitution V?–C, we do not remove the 0=V+C constraint. Our (my) reasoning was that the resulting 0=(–C)+C constraint would vanish after constant folding, anyway.

However, this is no longer true when we substitute a complex expression –E for V: The even more complex expression in 0=(–E)+E is not simplified to zero! What do we do?

I say we do what we should have done from the beginning: When we handle a constraint, we immediately remove it from the set of constraints. In concrete code, we pass a "source constraint" into SolverNode.SubstituteVariable; and that method will exclude this constraint when building a new SolverNode:

    private SolverNode SubstituteVariable(Variable variable,
              AbstractExpr expression,
              AbstractConstraint sourceConstraintToBeRemoved) {
         var rewriter = new RewritingVisitor(
            new Dictionary<AbstractExpr, AbstractExpr> {
                { variable, expression }
            });
        IEnumerable<AbstractConstraint> rewrittenConstraints =
            Constraints.Except(sourceConstraintToBeRemoved)
                        .Select(c2 => c2.Accept(rewriter, Ig.nore));
        return new SolverNode(rewrittenConstraints, Backsubstitutions, this);
    }

Two earlier tests that relied on the left-over 0=0 constraints now fail: Test2SimpleConstraintsAndOneStep and Test2SimpleConstraintsAndTwoSteps. But they are easily modified to reflect the new behavior—I simply delete half of the steps (and rename a few variables). When all the tests are green again, we can focus again on our first backsubstitution test. It is still red, but now it "fails correctly"—i.e., in the last line.

Implementing backsubstitutions—creating the backsubstitution

Here is a diagram of the data structure where we store knowledge about variables from now on, with the new VariableWithBacksubstitution class:

We call that knowledge now "closed variables." A variable is "closed" when it does no longer take part in the solving process—either because we know its value, or because we know how to compute its final value later by backsubstitutions. Creating such a backsubstition is easy: We just add it in SubstituteVariable if the substituted expression is not a Constant. Compare the following code with the one immediately above to see the modification:

    private SolverNode CloseVariable(Variable variable,
                                     AbstractExpr expression,
                                     AbstractConstraint sourceConstraintToBeRemoved) {
        expression = expression.Accept(new ConstantFoldingVisitor(), Ig.nore);

        // Rewrite all constraints
        var rewriter = new RewritingVisitor(
            new Dictionary<AbstractExpr, AbstractExpr> { { variable, expression } });
        IEnumerable<AbstractConstraint> rewrittenConstraints =
            Constraints.Except(sourceConstraintToBeRemoved)
                .Select(c => c.Accept(rewriter, Ig.nore));

        // Create new variable->value knowledge or new backsubstition.
        var newBacksubstitutions =  
            new Dictionary<Variable, AbstractClosedVariable> (_closedVariables) {{
              variable, expression is Constant
                ? (AbstractClosedVariable)
                    new VariableWithValue(variable, ((Constant) expression).Value)
                : new VariableWithBacksubstitution(variable, expression)
            }};

        // Create new node.
        return new SolverNode(rewrittenConstraints, newBacksubstitutions, this);
    }

After this change, our first test finally runs green!

Implementing backsubstitutions—rippling the backsubstitutions

The second test still fails, as expected. It does so in the last line but one, where it checks for the value of variable x. Debugging shows us that the final node still contains two backsubstitutions for x and z, as expected. Therefore, we must now add code to unwind the backsubstitutions remaining after actual solving the constraints.

The best place to do this is when we finally found a solution node, i.e., after the solver loop in Solve. The algorithm I chose is not that interesting—here is a rough outline:

• I keep two lists: One stores the variable-value pairs that we want to substitute, the other contains the open backsubstitutions.
• When we have substituted all values in all backsubstitutions, we are left over with two new lists: The values that now have emerged from completely subtituted backsubstitutions, and the backsubstitutions that are still expressions.
• We start the algorithm again with these two lists; and finish when all backsubstitutions have vanished.

For those who are interested, here is the algorithm as a pdf file.

And with this, both our backsubstitution tests are green! ... and we can, in the next posting, return to the rotating bar.