A tiny language for mechanical animations - the AWK script

The AWK compiler for my little animation language is quite simple. Right now, it has 104 lines, and I do not intend to expand it. You can download the 3kB file here: www.haraldmmueller.de/eb/animate.awk.

1. Overall structure

The overall structure of the script is simple: There is a single block for each of the 5 possible commands.

/^#/      { ...
          }
/^!/      { ...
          }
/^%/      { ...
          }
/^~/      { ...
          }
/^\$.../  { ...
          }

The state of the script is contained in the following eleven variables:

1. rem: the prefix for shell comments, to be set in the gawk call.
2. valuS[v]: value of scalar variable v.
3. stepS[v]: increment of scalar variable v.
4. moduS[v]: modulus of scalar variable v.
5. frmtS[v]: print format of scalar variable v (necessary because modulus controls leading zeros).
6. valuX[v]: x value of vector variable v.
7. valuY[v]: x value of vector variable v.
8. stepXOrM[v]: either a number → increment for x value of vector variable v; or a string starting with % → anchor for second format of vector increment.
9. stepYOrA[v]: increment for y value of vector variable v (if stepXOrM[v] is a number); or angle increment (if stepXOrM[v] starts with %).
10. currcmd: The previous command. 
11. out: The flag remembering the ~ state.

The following sections explain the five code pieces.

2. Explicit Comments

This is easy: We print rem, the line number, and the text; and then continue with the next line ("next"):

/^#/  { print rem " (" NR ":" $0 ")"
        next
      }

3.Scalar variable assignment

This straightforward code consists of four if statements:

• The first one checks that there are exactly three parameters present (a rudimentary attempt at error handling—more of this is needed ...). The print statements uses gawk's extension > "/dev/stderr" which works also under Windows.
• The other three store the respective values in valuS, stepS, and moduS, unless a dot was specified. frmtS is set when the value or the modulus is set.0 + ... is the AWK idiom for converting a string to a number:

 
/^!/  { if (NF != 4) {
          print NR ": ERROR - exactly 3 parameters expected in " $0 > "/dev/stderr"
          next
        }
        if ($2 != ".") {
          valuS[$1] = 0 + $2
          frmtS[$1] = "%0d"
        }
        if ($3 != ".") {
          stepS[$1] = 0 + $3
        }
        if ($4 != ".") {
          moduS[$1] = 0 + $4
          frmtS[$1] = "%0" int(log(moduS[$1]) / log(10) + 1) "d"
        }
        print rem " v=" valuS[$1] ...more debugging output...
        next
      }

4. Vector variable assignment

Also vector variable assignments just store the values:

/^%/  { if (NF != 5) {
          print NR ": ERROR - exactly 4 parameters expected in " $0 > "/dev/stderr"
          next
        }
        if ($2 != ".") {
          valuX[$1] = 0 + $2
        }
        if ($3 != ".") {
          valuY[$1] = 0 + $3
        }
        if ($4 != ".") {
          stepXOrM[$1] = $4
        }
        if ($5 != ".") {
          stepYOrA[$1] = 0 + $5
        }
        print rem " v=" valuS[$1] ":" valuY[$1] " ...more debugging output...
        next
      }

5. Handling ~

This is simple:

/^~/    { out = $2
          next
        }

6. Command handling

Obviously, this is the interesting part. A $ line must contain a subsequent number of steps. Therefore, the script checks not only for a $ at the beginning of the line, but also for a subsequent number. The steps are remembered, then we collect the command if it is present. At last, the command is emitted multiple times. Here is the overall structure:

/^\$[ \t]+[0-9]+[ \t]+/    {
        steps = 0 + $2               
        ...remove leading $ and $2
        if ($0 != "") {
          ...collect lines with trailing \ into currcmd
        }
        for (i = 0; i < steps; i++) {
          cmd = currcmd
          ...interpolate scalar variables in cmd
          ...interpolate vector variables in cmd
          print cmd
          ...increment scalar variables
          ...increment vector variables
        }
      }      

Removing the leading $ and count is easy in AWK:

        sub(/^\$[ \t]+[0-9]+[ \t]+/, "")

If there is a command after the count, we collect it in variable currcmd. Handling trailing \ for command concatenation requires five more lines:

        if ($0 != "") {
          currcmd = $0
          while (/[\\]$/) {
            sub(/[\\]$/, " ")     
            getline
            currcmd = currcmd $0
          }
        }

Inside the emitting loop, we handle variables, print the output line, and finally increment the variables.

Replacing scalars has a few special points:

• First, all values are rounded to the nearest integer. The reason is that the images work with pixels.
• If modulus is set, we take the variable value modulo that modulus
• And finally, we format the value using the current format.

          for (v in valuS) {
            vv = int(valuS[v]+0.5)
            if (moduS[v] > 0) vv %= moduS[v]
            vv = sprintf(frmtS[v], vv);
            gsub(v, vv, cmd)
          }         

Replacing vectors is actually easier than scalars, as there are no formats and modulus for vectors. Also here, we round to the nearest integer:
         
          for (v in valuX) {
            vx = int(valuX[v]+0.5);
            vy = int(valuY[v]+0.5);
            gsub(v ":x", vx, cmd)
            gsub(v ":y", vy, cmd)
          }

After (possibly) printing the command, we have to increment the values. For scalars, this is a "one-liner":

          for (v in valuS) {
            valuS[v] += stepS[v]
          }

For vectors, we have two variants. One is easy—we just add the increments. The 0 + is needed here because stepXOrM is not coerced to a number when reading because it might be a %name:

          for (v in valuX) {
            if (stepXOrM[v] ~ /^%/) {
              ...angle increment - see below...
            } else {
              valuX[v] += 0 + stepXOrM[v]
              valuY[v] += stepYOrA[v]
            }
          }

Finally, we need some simple trigonometry for the angle increment. Here is the algorithm:

• First, we need the coordinates of the rotation center (mx and my).
• Then, we compute the delta vector between the current location and the rotation center (dx and dy).
• We convert this to angular form (r and phi) and immediately add the angle increment, after converting it from degrees to radians.
• Finally, we compute the new location by transforming and adding the new delta vector to the rotation center:

              mx = valuX[stepXOrM[v]]
              my = valuY[stepXOrM[v]]
              dx = valuX[v] - mx
              dy = valuY[v] - my
              r = sqrt(dx*dx + dy*dy)
              phi = atan2(dy, dx) + (stepYOrA[v] / 180 * 3.1416)
              valuX[v] = mx + r * cos(phi)
              valuY[v] = my + r * sin(phi)

And this is it!

The next posting will contain two examples of simple animations and, additionally, some arguments on why cranks are difficult.

A tiny language for mechanical animations using ImageMagick

My Movimentum project is currently fast asleep—but in a discussion about railroad barriers, I felt that I needed to explain a certain mechanical linkage with animations. Christmas season is a good time to experiment with new ideas (at least new to me), so I created three "mechanical explanations" for an old Austrian gear that drove (and drives) railroad barriers. You can see the results in that posting (in German). If you are interested in what it explains, just ask here—then I'll give you a quick overview about what this is all about. However, this posting here is really about the tool chain I used to create these animations.

Here is the second animation I made:

Here is an overview of my process. Below, I will share my experiences for each step:

1. I create the parts of the animation with a CAD tool (I use an older version of Caddy++).
2. The result is exported ("printed") into a single PNG file with PDFCreator.
3. The PNG file is split into multiple PNGs, each of which gets background transparency (I use an old version of Paint Shop Pro here).
4. Next, the images need "location calibration": In step 3., they do not come out exactly overlapping. The small corrections to move them to the right position are done here, using ImageMagick.
5. Next, I write the animation script. For this, I defined a very simple language that is interpreted by a 100-line AWK program (I use Gnu's Gawk for this). The description of this language is the major part of this posting.
6. The result of this is a batch file that calls ImageMagick's convert program to assemble the images into frames.
7. At last, ffmpeg is used to create an animation file from the frames.

For the record, I work with Windows. However, I think most of this also applies e.g. to Linux.

1. Creating the parts

Most probably, the parts have to be aligned in the animation. For this, it is very useful to add cross hairs as anchor points at critical locations, e.g., at the center of rotating items.

Moreover, to place the anchors at a predictable point, I put each part into a box of equal size, where the anchors are at the same position relative to the left upper box corner. In addition, I create an empty box that contains only the cross hairs of an anchor.

Here is a part of such a drawing. The lower left box is the empty one.

2. Creating PNG files

Instead of snipping the parts images from the CAD program or a PDF rendering, PDFCreator can be installed as a virtual printer that directly creates PNG files (despite its name). This avoids color and other distortions. Currently, I create a single PNG file with all the parts—this is faster than creating a single PNG file per part by printing from the CAD program, because I have to adjust each parts file anyway in the next step.

Already here, I select the correct scaling of the PNG—at least with Caddy++, this is possible in the print dialog. This avoids any later scalings which will destroy the cross hairs of the anchor points.

PDFCreator creates the PNG in portrait orientation even if you "print" in landscape—I do not know whether this is a fault of PDFCreator or Caddy++ or whatever else. In Windows, the created PNG is immediately opened in the Photo Viewer, so it is only a single click to rotate it to the intended orientation.

3. Splitting the PNG image and adding global transparency

Using Paint Shop Pro, I select the content of each box in the result and paste each part into a new PNG file. Unfortunately, box selection uses screen coordinates in Paint Shop Pro. Therefore, in the new files, the pasted parts may not lie at perfectly the same relative spot—I consistently got offsets of ±1 for both coordinates. In step 4 I explain how I correct this—but maybe there is a better way ... Because of this correction, I save the parts in files called ...Approx.PNG; the corrected files will have names without Approx.

While I have the new part file open, I make the background transparent, because this is what is needed for stacking the parts in the animation. Also the empty box gets its own "part file" with transparent background.

4. "Location calibration"

As explained above, the anchor points of the parts files will not be exactly at the same relative position to the upper left corner. This is a nuisance, as ImageMagick uses the upper left corner as a reference point for its computations, whereas we want to use our anchors. I repair this problem manually as follows:

• With ImageMagick, I overlay the the empty box over each part.
• In the resulting image, I check the distance of the part anchors vs. the empty box anchors.
• Again using ImageMagick, I correct the location of the part.

Here are the concrete steps for this on the command line, using ImageMagick's convert:

• Overlay and view with convert Part1Approx.PNG Empty.PNG -geometry +0+0 -composite temp.png
• Viewing temp.png reveals that the anchors can e.g. be superimposed precisely with +1+-2 instead of +0+0. This is a trial-and-error process. Nicely enough, Windows's Photo Viewer will refresh the image immediately, so this is a very quick process.
• Knowing the size of the empty box (e.g. use ImageMagick's identify—here I use 565x612 as example), invert the correction from +1+-2 to +-1+2, because we now want to shift the part's image, not the empty box. Then create the perfect part file via convert -size 565x612 xc:none Part1Approx.PNG -geometry +-1+2 -composite Part1.PNG.

5. The animation script

Now we come to the central piece: The animation script. The basic idea is that we create a separate ImageMagick convert command for each frame. Each such command will place the parts at specific places. But of course, instead of writing single commands for thousands of frames, I write patterns that contain variables that are expanded many times. The script is then run through an Gnu AWK program that emits the convert commands that compute the frames. Here is the call I use:

  gawk -v rem=REM -f animate.awk script.txt > script.bat

I will explain the "script compiler" animate.awk in a separate posting.

Update: During documenting the small compiler, I found I could make the language a little bit terser: The = step definition has been merged into the $ and + commands.

5.1. Command patterns and variables - overview

Here is such a command pattern that overlays a Part1 at a fixed place and a Part2 after a rotation and a shift (you can look up convert's options in ImageMagick's documentation):

... convert -size 900x800 xc:white \
    Part1.PNG -geometry +50+600 -composite \
    ( Part2.PNG -distort SRT "333,329,-!A2" ) \
      -geometry +%2:x+%2:y -composite \
  out\f_!T.png

The variables in this pattern are the strings starting with ! or %:

• !... is a scalar variable, Very often this is used for angles (like !A2 above), but sometimes also for numbering of files (like the !T at the end of the command).
• %... is a vector variable. Its components are accessed with %...:x and %...:y. In the command above, a vector variable %2 is used to shift the rotated Part2.

5.2. Script segments

A complete script consists of script segments where each segment describes a homogeneous movement. Each segment consists of two parts:

• A (possibly empty) sequence of variable assignments.
• A command pattern (convert patterns usually become very long, therefore they can be broken over multiple lines with trailing backslashes, as shown above).

Here is an example of such a segment:

  !T 10000 1 .
  $ 4 convert -size 900x800 xc:white \
      Text1.PNG -geometry +50+600 -composite out\f_!T.png

• !T 10000 1 . is an assignment to the variable !T: Its initial value is 10000, and it is incremented by 1 for each emitted convert. The last parameter (which is a single dot here) will be explained later.
• $ 4 convert ... means that the convert command is to be emitted 4 times, with variables substituted in each command.

This script segment will emit four lines where the variable !T is replaced with 10000, 10001, 10002 and 10003 in turn:

  convert -size .... out\f_10000.png
  convert -size .... out\f_10001.png
  convert -size .... out\f_10002.png
  convert -size .... out\f_10003.png  

5.3. The details of scalar variables

The assignment to a scalar variable has the following form:

  !name  initial  increment  modulus

name can consist of any characters you like that are not special in regexps (reason: In the AWK script, regexp replacement is used to interpolate the variables).

The three parameters can either be numbers (decimals are ok) or single dots (.). The dot means that the value remains as it was before. This is a central concept that helps to write short scripts. E.g., here is how you move something linearly, then let it stop:

  # 15 frames that reduce !X from 200 to some value.
  !X 200 -10 . 
  $ 15 convert ...

  # 10 frames without movement
  !X . 0 .
  $ 10 convert ...

In the second assignment, we left !X at whatever value it had after the first segment completed by specifying initial as . (dot). However, we set the increment to 0, so that this variable does not change its value in the 10 frames of this segment.

The modulus is sometimes useful for angle variables: The replacement value for a pattern is the variable's value modulo the modulus. I needed this because I precomputed rotated images of a part into 360 files called something like Part1_000.PNG, Part1_001.PNG,..., Part1_359.PNG. Each file contained the part rotated by so-and-so many degrees. For rotating the part, I could now write something like

  !A 350 2 360
  $ 20 convert ... Part_!A.PNG -geometry +0+0 -composite ...

Starting at 350 (degrees), the variable !A will be incremented by 2 in each of 20 frames. After 5 frames, !A will be 360 and continue with 362, 364 etc. up to 390. Without 360, the calls using Part_360.PNG etc. will fail, as these files do not exist. However, with modules 360, the emitted commands will contain Part_000.PNG, Part_002.PNG etc. The modulus, if defined, also determines the number of digits in the output number: Because modulus is 360 here, we get three-digit numbers.

5.4. The details of vector variables

A vector variable consists of two components, x and y. There are two possible assignments to vector variables. The first form does what you expect:

  %name  initialX  initialY  incrementX  incrementY

The variable is set and/or incremented as specified. Instead of a number, each argument can again be a single dot, which keeps the current value.

The other form is necessary for rotational movements:

  %name  initialX  initialY  %other  degrees

This means that the change to %name is a rotation around %other by the angle specified as the last argument. %other cannot be a dot here (the % is the marker that this form is intended), but the other three values can be numbers or dots. If  %other itself changes, the results are somewhat unpredictable, as the order of variable computations is not defined. If someone has the idea to try this—well, it's on your own risk ...

5.5. The details of commands

There are two possible forms for commands. One is

  $ count command-text

This creates count instances of command-text. The variables are replaced in each instance, and the result is written to the output. As explained above, long commands can be wrapped with \ at the lines' ends (the \ must actually be the last character on the line—no trailing white space!).

text can be anything your operating systems understands. I have, up to now, only called ImageMagick's convert here. However, one could e.g. write intermediary files and use them later in the frame sequence, e.g. to improve the performance of frame generation (thousands of calls to convert are not that fast ...).

It makes much sense to write the output of commands—especially the created frame images—to a specific output directory, to keep things cleanly separated.

The other command form omits the command-text:

  $ count

This means that the previous command is to be repeated in this segment. This is a very important concept: It allows one to specify the "model" in a single convert command that places all the parts at various places; and then have this assembly move simply by maniulating variables. One could even argue that a good animation script sets up a model once at the beginning; and then only modifies variables.

In addition to the $ lines, it is possible to switch off the output with

  ~ -

and switch it on at a specific point with

  ~ +

This is very helpful for debugging: It creates a script that will start "later in the animation", so that the script will only compute certain frames that are to be fine-tuned. This is preferable to (and much quicker than) editing the script each time it has been written before running it.

5.6. Comments

There are two possible forms of comments. One is

  # text

The text of such a comment is copied into the output, with a leading comment marker that must be passed to the compiler. This is helpful for debugging, because you can line up the script with the output.

The other form is everything else which is not understood by animation.awk, i.e., every line that does not start with !, %, ~, $, or #. Such lines are simply ignored. I use this feature to embed the script into a text that explains the animations, which will later become a web page. However, one could call this an unfortunate idea (even though it is along the lines of Donald Knuth's Literate Programming).

6. Running the batch file

This is as easy as it gets. You just run it.

7. Creating the animation

At last, ffmpeg is used to create an animation file from the frames. Here is the command line I use:

ffmpeg -f image2 -r 12 -i f_1%04d.png 
             -vcodec libx264 -pix_fmt yuv420p -y out.mp4

Of course, you can use other formats than mp4 for your animation (especially if you do not like proprietory formats). For mp4, the -pix_fmt yuv420p is important. The default pixel format, whose name I forgot, is not understood by the majority of programs or web sites out there.

If you want to create multiple animations from a single script (which happens if you embed your script code into running text, as I do), it makes sense to start the frame files with new numbers for each script, using !T 10000 . ., !T 20000 . ., !T 30000 . . etc. Then, you can run ffmpeg with file patterns f_1%04d, f_2%04d, f_3%04d etc.

In my next posting, I explain the small compiler (which you can download here).
In a subsequent posting, I will give a few animation examples and some background on why cranks are difficult.

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.

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 √4x², 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 √(P²Q). 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.