Wednesday, 27 March 2019

The Differentiator

Fun With C# Expressions

Previously:
Differentiating f(x)^g(x)
Note: The latest code can be found at https://github.com/dogbiscuituk/Sid/.

Take a seat with Visual Studio and create a new solution - Console, Winforms, whatever - your choice. Now add a new class library FormulaBuilder, and add to that a static class Differentiator. This is going to let us manipulate objects of type System.Linq.Expression, representing real-valued functions of a single real variable, and ultimately to plot the graph of any such function, as well as its nth derivatives. It will use an easy-to-read fluent interface, while illustrating some key design tenets I'd like to highlight.

Back to the class, and here come our first three requirements. We will need ways to:
  1. represent constants and variables, and combine these fluently into compound expressions;
  2. display any such expression as an easy-to-read formula, and evaluate it for given values of x; and
  3. verify correctness of these combinations and conversions.
Leaving aside the whole "combine these fluently into compound expressions" thing for just a moment, let's open with a ParameterExpression for x, and the following starter set of extension methods:

using System;
using System.Diagnostics;
using System.Linq.Expressions;

namespace FormulaBuilder
{
    public static class Differentiator
    {
        public static ParameterExpression x = Expression.Variable(typeof(double));

        public static ConstantExpression Constant(double c) => Expression.Constant(c);

        public static string AsString(this Expression e, string variableName = "x") =>
            e.ToString().Replace("Param_0", variableName);

        public static Func<double, double> AsFunction(this Expression e) =>
            Expression.Lambda<Func<double, double>>(e, x).Compile();

        public static double AsDouble(this Expression e, double x) => AsFunction(e)(x);

        public static void Check(double expected, double actual) =>
            Check(expected.ToString(), actual.ToString());

        public static void Check(string expected, string actual)
        {
            if (actual == expected)
                Debug.WriteLine("OK: " + actual);
            else
                Debug.WriteLine("*** Comparison failed. Expected: {0}, Actual: {1}.", expected, actual);
        }

On first look, the expression display requirement (1) seems to be met by the built-in ToString() method, available to every type; but that spits out an ugly Param_0 instead of a nice tidy x, so we supply a new AsString() method to tidy that up.

The expression evaluation/iteration requirement (2) means converting our expression into a LambdaExpression, compiling that into a function, and finally executing that function using our chosen values of x. For efficiency I've split this process into two steps. The conversion to a compiled LambdaExpression is performed by the ToFunction() method. Once we have this function, we can call it a million times with different x values, without having to perform a runtime lambda recompilation every time. The ToDouble() method is also provided for those out-of-loop test cases, where we just want to obtain a single value using a single method call.

As for testing requirement (3), the overloaded Check methods provide a line of debug output for each checked operation, beginning with either OK or *** Comparison failed.

Fluency First

Now back to those compounds. Actually, the System.Linq.Expression class already has some fine methods for compounding expressions into BinaryExpressions, but I'm going to deviate a little from these. The target is a fluent syntax for building expressions on the fly, and to that end, I find the names Plus, Minus, Times, Over a little easier to scan than Add, Subtract, Multiply, Divide. So let's add these to the Differentiator class:

        public static BinaryExpression Plus(this Expression f, Expression g) => Expression.Add(f, g);
        public static BinaryExpression Minus(this Expression f, Expression g) => Expression.Subtract(f, g);
        public static BinaryExpression Times(this Expression f, Expression g) => Expression.Multiply(f, g);
        public static BinaryExpression Over(this Expression f, Expression g) => Expression.Divide(f, g);
        public static BinaryExpression Power(this Expression f, Expression g) => Expression.Power(f, g);

Actually I also considered replacing Power with ToThe, but maybe that's just a colloquial idiom too far!

The Constant() syntax is quite clumsy, almost the opposite of the fluent design goal. It would be great to have an implicit operator for converting constants such as 2.0 or Math.PI into ConstantExpressions as required - just as it would be equally great to be able to overload the standard arithmetic operators for Expressions, and be finished with fluent method calls! Sadly, C# restrictions on overloaded operators make both of these dreams impossible, but at least in the Constant() case, what we can do instead is overload the calling methods at the point(s) of invocation. For example, we'd like to be able to write x.Power(2) or x.Power(2.0) instead of x.Power(Constant(2)). This can be achieved by overloading the Power() method, supplying a double in place of an Expression for its second argument, and similarly for the other BinaryExpression methods:

        public static BinaryExpression Plus(this Expression f, double g) => f.Plus(Constant(g));
        public static BinaryExpression Minus(this Expression f, double g) => f.Minus(Constant(g));
        public static BinaryExpression Times(this Expression f, double g) => f.Times(Constant(g));
        public static BinaryExpression Over(this Expression f, double g) => f.Over(Constant(g));
        public static BinaryExpression Power(this Expression f, double g) => f.Power(Constant(g));

This almost gets us where we want to go, but we still have to write e.g. Constant(3).Times(x), when we'd prefer just 3.Times(x). Here it's the first argument of the Times() extension method that needs overloading, and there's a spanner thrown in by C#'s rules of method resolution. Sure, we can supply alternate methods with a double in that position, but where previously we could rely on an int argument being promoted to double, no such luck this time. If we insist (and we do!) upon being able to write 3.Times(x), as opposed to workarounds like x.Times(3), 3.0.Times(x), or worst of all ((double)3).Times(x) - which seems hardly an improvement! - then we must supply separate double and int versions.

        public static BinaryExpression Plus(this double f, Expression g) => Constant(f).Plus(g);
        public static BinaryExpression Minus(this double f, Expression g) => Constant(f).Minus(g);
        public static BinaryExpression Times(this double f, Expression g) => Constant(f).Times(g);
        public static BinaryExpression Over(this double f, Expression g) => Constant(f).Over(g);
        public static BinaryExpression Power(this double f, Expression g) => Constant(f).Power(g);

        public static BinaryExpression Plus(this int f, Expression g) => Constant(f).Plus(g);
        public static BinaryExpression Minus(this int f, Expression g) => Constant(f).Minus(g);
        public static BinaryExpression Times(this int f, Expression g) => Constant(f).Times(g);
        public static BinaryExpression Over(this int f, Expression g) => Constant(f).Over(g);
        public static BinaryExpression Power(this int f, Expression g) => Constant(f).Power(g);

        public static UnaryExpression Negate(this Expression e) => Expression.Negate(e);
        public static BinaryExpression Reciprocal(this Expression e) => e.Power(-1);
        public static BinaryExpression Squared(this Expression e) => e.Power(2);
        public static BinaryExpression Cubed(this Expression e) => e.Power(3);

        public static void TestCompoundExpression()
        {
            var f = x.Squared().Plus(3.Times(x)).Minus(5);    // f(x) = x²+3x-5
            Check("(((x ^ 2) + (3 * x)) - 5)", f.AsString()); // Check the built expression formula
            Check(Math.Pow(7, 2) + 3 * 7 - 5, f.AsDouble(7)); // Check the expression value at x = 7 (should be 65)
        }

With the addition of these few minor utilities, we now have the means to build polynomials and other arithmetic expressions in x. Give your app a project reference to FormulaBuilder, add it to your using clauses, then add a call to Differentiator.TestCompoundExpression() somewhere in your main code and run it. You should see the following two lines of reassurance in your debug output:
OK: (((x ^ 2) + (3 * x)) - 5)
OK: 65
These show that the quadratic expression x² + 3x - 5 was built successfully, and that it correctly evaluates to 65 at x = 7.

Function Calls

Arithmetic operations are one way of building out our expressions, but what about invoking functions? What use a differentiator that can't tell you the derivative of sin is cos, while that of cos is -sin?

Elementary functions in the System.Linq.Expression hierarchy are represented by the MethodCallExpression class. The Function() code below generates such an expression from the correspondingly named method of the Math class, obtained via reflection. Below it are listed a generous dozen or so concrete functions, all which happen to be implemented in Math; and below these, a further group showing how functions "missing" from Math, such as the reciprocal trigonometrics and hyperbolics, as well as their inverses, can be built out of expressions in the elementary ones.

Naturally, the input to any functional expression can itself be an arbitrary expression tree, maybe containing further nested function calls, etc.

        public static MethodCallExpression Function(string functionName, Expression e) =>
            Expression.Call(typeof(Math).GetMethod(functionName, new[] { typeof(double) }), e);

        public static MethodCallExpression Abs(this Expression e) => Function("Abs", e);
        public static MethodCallExpression Sqrt(this Expression e) => Function("Sqrt", e);
        public static MethodCallExpression Exp(this Expression e) => Function("Exp", e);
        public static MethodCallExpression Log(this Expression e) => Function("Log", e);
        public static MethodCallExpression Log10(this Expression e) => Function("Log10", e);
        public static MethodCallExpression Sin(this Expression e) => Function("Sin", e);
        public static MethodCallExpression Cos(this Expression e) => Function("Cos", e);
        public static MethodCallExpression Tan(this Expression e) => Function("Tan", e);
        public static MethodCallExpression Asin(this Expression e) => Function("Asin", e);
        public static MethodCallExpression Acos(this Expression e) => Function("Acos", e);
        public static MethodCallExpression Atan(this Expression e) => Function("Atan", e);
        public static MethodCallExpression Sinh(this Expression e) => Function("Sinh", e);
        public static MethodCallExpression Cosh(this Expression e) => Function("Cosh", e);
        public static MethodCallExpression Tanh(this Expression e) => Function("Tanh", e);

        public static BinaryExpression Sec(this Expression e) => Reciprocal(Cos(e));
        public static BinaryExpression Csc(this Expression e) => Reciprocal(Sin(e));
        public static BinaryExpression Cot(this Expression e) => Reciprocal(Tan(e));
        public static MethodCallExpression Asec(this Expression e) => Acos(Reciprocal(e));
        public static MethodCallExpression Acsc(this Expression e) => Asin(Reciprocal(e));
        public static MethodCallExpression Acot(this Expression e) => Atan(Reciprocal(e));
        public static BinaryExpression Sech(this Expression e) => Reciprocal(Cosh(e));
        public static BinaryExpression Csch(this Expression e) => Reciprocal(Sinh(e));
        public static BinaryExpression Coth(this Expression e) => Reciprocal(Tanh(e));
        public static MethodCallExpression ArcSinh(this Expression e) => Log(e.Plus(Sqrt(e.Squared().Plus(1))));
        public static MethodCallExpression ArcCosh(this Expression e) => Log(e.Plus(Sqrt(e.Squared().Minus(1))));
        public static BinaryExpression ArcTanh(this Expression e) => Log(1.Plus(e).Over(1.Minus(e))).Over(2);
        public static MethodCallExpression ArcSech(this Expression e) => Log(1.Plus(Sqrt(1.Minus(e.Squared()))).Over(e));
        public static MethodCallExpression ArcCsch(this Expression e) => Log(1.Plus(Sqrt(1.Plus(e.Squared()))).Over(e));
        public static BinaryExpression ArcCoth(this Expression e) => Log(e.Plus(1).Over(e.Minus(1))).Over(2);

        public static void TestTrigonometricExpression()
        {
            var f = Sin(x).Squared().Plus(Cos(x).Squared());             // f(x) = sin²x+cos²x
            Check("((Sin(x) ^ 2) + (Cos(x) ^ 2))", f.AsString());        // Check the built expression formula
            double p3 = Math.PI / 3, s = Math.Sin(p3), c = Math.Cos(p3); // Check the expression value at x = π/3
            Check(Math.Pow(s, 2) + Math.Pow(c, 2), f.AsDouble(p3));      // (should be 1)
        }

In your app level code, add a call to Differentiator.TestTrigonometricExpression(), which builds the expression sin²x + cos²x: the LHS of the most famous trigonometric identity, whose RHS is 1 for all values of x. The new test method checks, first using AsString(), that the correct expression has in fact been built; then, using AsDouble(), that it evaluates to 1 as expected, in this case when supplied with π/3 as a parameter. At runtime, your debug window should display a further two comforting lines of verification:
OK: ((Sin(x) ^ 2) + (Cos(x) ^ 2))
OK: 1
Finally, A Derivative

Can't finish this article without a look at some actual differentiation, so let's add the function D() that differentiates expressions. It needs to handle just these four cases:
  • a constant;
  • x itself;
  • a function call;
  • an artithmetic operation (unary or binary).
The third and fourth cases here have potential for recursion. At long last, here is the code to differentiate any expression:

        public static Expression D(this Expression e)
        {
            if (e is ConstantExpression) return Constant(0);  // d(c)/dx = 0
            if (e is ParameterExpression) return Constant(1); // d(x)/dx = 1
            if (e is MethodCallExpression m)
            {
                var f = m.Arguments[0];
                return DifferentiateFunction(m.Method.Name, f).Times(D(f)); // Chain Rule
            }
            if (e is UnaryExpression u)
            {
                var f = D(u.Operand);
                return u.NodeType == ExpressionType.UnaryPlus ? f : Negate(f);
            }
            if (e is BinaryExpression b)
            {
                Expression f = b.Left, g = b.Right;
                switch (b.NodeType)
                {
                    case ExpressionType.Add:      // (f+g)' = f'+g'
                        return D(f).Plus(D(g));
                    case ExpressionType.Subtract: // (f-g)' = f'-g'
                        return D(f).Minus(D(g));
                    case ExpressionType.Multiply: // (fg)' = f'g+fg'
                        return D(f).Times(g).Plus(f.Times(D(g)));
                    case ExpressionType.Divide:   // (f÷g)' = (f'g-fg')÷(g^2) = f'÷g-fg'÷(g^2)
                        return D(f).Over(g).Minus(f.Times(D(g)).Over(g.Squared()));
                    case ExpressionType.Power:    // (f^g)' = (f^g)*(f'g÷f+g'Log(f)) = (f'g+fg'Log(f))f^(g-1)
                        return D(f).Times(g).Plus(f.Times(D(g)).Times(Log(f))).Times(f.Power(g.Minus(1)));
                }
            }
            throw new InvalidOperationException();
        }

        private static Expression DifferentiateFunction(string methodName, Expression x)
        {
            switch (methodName)
            {
                case "Abs": return x.Over(Abs(x));                            // d(|x|)/dx = x/|x|
                case "Sqrt": return x.Power(-0.5).Over(2);                    // d(√x)/dx = 1/(2√x)
                case "Exp": return Exp(x);                                    // d(eˣ)/dx = eˣ
                case "Log": return Reciprocal(x);                             // d(ln(x))/dx = 1/x
                case "Log10": return Reciprocal(x).Times(Math.Log10(Math.E)); // d(log₁₀(x))/dx = log₁₀(e)/x
                case "Sin": return Cos(x);                                    // d(sin(x))/dx = cos(x)
                case "Cos": return Negate(Sin(x));                            // d(cos(x))/dx = -sin(x)
                case "Tan": return Cos(x).Power(-2);                          // d(tan(x))/dx = sec²(x)
                case "Asin": return 1.Minus(x.Squared()).Power(-0.5);         // d(arcsin(x))/dx = 1/√(1-x²)
                case "Acos": return Negate(1.Minus(x.Squared()).Power(-0.5)); // d(arccos(x))/dx = -1/√(1-x²)
                case "Atan": return 1.Over(1.Plus(x.Squared()));              // d(arctan(x))/dx = 1/(1+x²)
                case "Sinh": return Cosh(x);                                  // d(sinh(x))/dx = cosh(x)
                case "Cosh": return Sinh(x);                                  // d(cosh(x))/dx = sinh(x)
                case "Tanh": return 1.Minus(Tanh(x).Squared());               // d(tanh(x))/dx = 1-tanh²(x)
            }
            throw new InvalidOperationException();
        }

DifferentiateFunction() provides the basic template for a function's derivative, without actually applying the chain rule to its argument Expression. That's why it has been implemented as a private helper function; calling it out of context could result in incorrect results if chain rule application was omitted by mistake.

Algebraic Simplification

Still we're not quite ready to supply the final batch of unit tests, stick a fork in it and call it done. While the D() method may be returning strictly correct results, good enough for numeric applications such as graphing, algebraically they are a bit of a mess. For example, differentiating any candidate containing a function call, even one as simple as Sin(x), results in a clumsy *1 becoming appended to the output expression. Some level of simplification is clearly needed to provide visually acceptable results. For this reason, D() gets called not directly, but from within a wrapper function Differentiate() which applies a layer of simplification before returning the resultant Expression:

        public static Expression Differentiate(this Expression e) => Simplify(D(e));

        public static Expression Simplify(this Expression e)
        {
            if (e is MethodCallExpression m) return SimplifyMethodCall(m);
            if (e is UnaryExpression u) return SimplifyUnary(u);
            if (e is BinaryExpression b) return SimplifyBinary(b);
            return e;
        }

        public static Expression SimplifyMethodCall(MethodCallExpression m)
        {
            var operand = Simplify(m.Arguments[0]);
            if (operand is ConstantExpression ce)
            {
                var c = (double)ce.Value;
                return Constant(Function(m.Method.Name, x).AsDouble(c));
            }
            return Function(m.Method.Name, operand);
        }

        public static Expression SimplifyUnary(UnaryExpression u)
        {
            var operand = Simplify(u.Operand);
            if (operand is ConstantExpression ce)
            {
                var c = (double)ce.Value;
                return Constant(u.NodeType == ExpressionType.UnaryPlus ? c : -c);
            }
            return u.NodeType == ExpressionType.UnaryPlus ? operand : Negate(operand);
        }

        public static Expression SimplifyBinary(BinaryExpression b)
        {
            b = Expression.MakeBinary(b.NodeType, Simplify(b.Left), Simplify(b.Right));
            if (b.Left is ConstantExpression ce1)
            {
                var c1 = (double)ce1.Value;
                if (b.Right is ConstantExpression ce2)
                    return Constant(b.NodeType, c1, (double)ce2.Value);
                switch (b.NodeType)
                {
                    case ExpressionType.Add:                   // c + x → x + c
                    case ExpressionType.Multiply:              // c * x → x * c
                        b = Expression.MakeBinary(b.NodeType, b.Right, b.Left);
                        break;
                    case ExpressionType.Subtract:
                        if (c1 == 0) return Negate(b.Right);   // 0 - x → -x
                        break;
                    case ExpressionType.Divide:
                        if (c1 == 0) return Constant(0);       // 0 ÷ x → 0
                        break;
                    case ExpressionType.Power:
                        if (c1 == 0) return Constant(0);       // 0 ^ x → 0
                        if (c1 == 1) return Constant(1);       // 1 ^ x → 1
                        break;
                }
            }
            if (b.Right is ConstantExpression ce3)
            {
                var c3 = (double)ce3.Value;
                switch (b.NodeType)
                {
                    case ExpressionType.Add:
                    case ExpressionType.Subtract:
                        if (c3 == 0) return b.Left;            // x ± 0 → x
                        break;
                    case ExpressionType.Multiply:
                        if (c3 == 0) return Constant(0);       // x * 0 → 0
                        if (c3 == 1) return b.Left;            // x * 1 → x
                        if (c3 == -1) return Negate(b.Left);   // x * -1 → -x
                        break;
                    case ExpressionType.Divide:
                        if (c3 == 0)
                            throw new DivideByZeroException();
                        if (c3 == 1) return b.Left;            // x ÷ 1 → x
                        if (c3 == -1) return Negate(b.Left);   // x ÷ -1 → -x
                        break;
                    case ExpressionType.Power:
                        if (c3 == 0) return Constant(1);       // x ^ 0 → 1
                        if (c3 == 1) return b.Left;            // x ^ 1 → x
                        break;
                }
                switch (b.NodeType)
                {
                    case ExpressionType.Add:                   // (x + c1) ± c2 → x + (c1 ± c2)
                    case ExpressionType.Subtract:              // (x - c1) ± c2 → x - (c1 ∓ c2)
                    case ExpressionType.Multiply:              // (x * c1) */ c2 → x * (c1 */ c2)
                    case ExpressionType.Divide:                // (x / c1) */ c2 → x / (c1 /* c2)
                        if (!(b.Left is BinaryExpression c && c.Right is ConstantExpression ce4))
                            break;
                        bool same = c.NodeType == b.NodeType, opps = c.NodeType == Invert(b.NodeType);
                        if (!(same || opps))
                            break;
                        var c4 = (double)ce4.Value;
                        var nodeType = c.NodeType == ExpressionType.Add || c.NodeType == ExpressionType.Subtract
                            ? ExpressionType.Add
                            : ExpressionType.Multiply;
                        if (opps) nodeType = Invert(nodeType);
                        return Expression.MakeBinary(c.NodeType, c.Left, Constant(nodeType, c4, c3));
                }
            }
            return b;
        }

        public static ExpressionType Invert(this ExpressionType nodeType)
        {
            switch (nodeType)
            {
                case ExpressionType.Add: return ExpressionType.Subtract;
                case ExpressionType.Subtract: return ExpressionType.Add;
                case ExpressionType.Multiply: return ExpressionType.Divide;
                case ExpressionType.Divide: return ExpressionType.Multiply;
            }
            throw new InvalidOperationException();
        }

        public static ConstantExpression Constant(ExpressionType nodeType, double c, double d) =>
            Constant(Apply(nodeType, c, d));

        public static double Apply(ExpressionType nodeType, double c, double d)
        {
            switch (nodeType)
            {
                case ExpressionType.Add: return c + d;
                case ExpressionType.Subtract: return c - d;
                case ExpressionType.Multiply: return c * d;
                case ExpressionType.Divide:
                    if (d == 0)
                        throw new DivideByZeroException();
                    return c / d;
                case ExpressionType.Power:
                    if (c == 0 && d == 0)
                        throw new InvalidOperationException();
                    return Math.Pow(c, d);
            }
            throw new InvalidOperationException();
        }

This Simplify() method does little more than apply some arithmetic identity, associativity and commutativity results, to remove some of the redundancy littering the differentiator's output. That's why I've implemented it here as simple, direct logic. There are many more simplification rules to consider, and as such rules are added, it becomes clear that a rule application engine is needed, allowing extension of the list without spaghettification of the logic.

Now we can at least start writing an acceptable suite of unit tests. Here are some examples:

        public static void TestSimplify(Expression e, string expected) => Check(expected, Simplify(e).AsString());

        public static void TestSimplifications()
        {
            TestSimplify(x.Plus(6).Plus(2), "(x + 8)");
            TestSimplify(6.Plus(x).Plus(2), "(x + 8)");
            TestSimplify(6.Plus(x.Plus(2)), "(x + 8)");
            TestSimplify(x.Plus(6).Minus(2), "(x + 4)");
            TestSimplify(6.Plus(x).Minus(2), "(x + 4)");
            TestSimplify(x.Minus(6).Plus(2), "(x - 4)");
            TestSimplify(x.Minus(6).Minus(2), "(x - 8)");
            TestSimplify(x.Times(6).Times(2), "(x * 12)");
            TestSimplify(6.Times(x).Times(2), "(x * 12)");
            TestSimplify(6.Times(x.Times(2)), "(x * 12)");
            TestSimplify(x.Times(6).Over(2), "(x * 3)");
            TestSimplify(6.Times(x).Over(2), "(x * 3)");
            TestSimplify(x.Over(2).Times(8), "(x / 0.25)");
            TestSimplify(x.Over(2).Over(8), "(x / 16)");
        }

        public static void TestDerivative(Expression e, string expected) => Check(expected, Differentiate(e).AsString());

        public static void TestFunctionDerivatives()
        {
            TestDerivative(Abs(x), "(x / Abs(x))");                       // d(|x|)/dx = x/|x|
            TestDerivative(Sqrt(x), "((x ^ -0.5) / 2)");                  // d(√x)/dx = 1/(2√x)
            TestDerivative(Exp(x), "Exp(x)");                             // d(eˣ)/dx = eˣ
            TestDerivative(Log(x), "(x ^ -1)");                           // d(ln(x))/dx = 1/x
            TestDerivative(Log10(x), "((x ^ -1) * 0.434294481903252)");   // d(log₁₀(x))/dx = log₁₀(e)/x
            TestDerivative(Sin(x), "Cos(x)");                             // d(sin(x))/dx = cos(x)
            TestDerivative(Cos(x), "-Sin(x)");                            // d(cos(x))/dx = -sin(x)
            TestDerivative(Tan(x), "(Cos(x) ^ -2)");                      // d(tan(x))/dx = sec²(x)
            TestDerivative(Asin(x), "((1 - (x ^ 2)) ^ -0.5)");            // d(arcsin(x))/dx = 1/√(1-x²)
            TestDerivative(Acos(x), "-((1 - (x ^ 2)) ^ -0.5)");           // d(arccos(x))/dx = -1/√(1-x²)
            TestDerivative(Atan(x), "(1 / ((x ^ 2) + 1))");               // d(arctan(x))/dx = 1/(x²+1)
            TestDerivative(Sinh(x), "Cosh(x)");                           // d(sinh(x))/dx = cosh(x)
            TestDerivative(Cosh(x), "Sinh(x)");                           // d(cosh(x))/dx = sinh(x)
            TestDerivative(Tanh(x), "(1 - (Tanh(x) ^ 2))");               // d(tanh(x))/dx = 1-tanh²(x)
        }

        public static void TestPolynomialDerivative()
        {
            TestDerivative(x.Power(4).Minus(3.Times(x.Cubed())).Plus(6.Times(x.Squared())).Minus(3.Times(x)).Plus(1),
                "(((((x ^ 3) * 4) - ((x ^ 2) * 9)) + (x * 12)) - 3)");    // d(x⁴-3x³+6x²-3x+1)/dx = 4x³-9x²+12x-3
        }

        public static void TestChainRule()
        {
            TestDerivative(Exp(x.Squared()), "(Exp((x ^ 2)) * (x * 2))"); // d(exp(x²))/dx = exp(x²)*2x
            TestDerivative(Log(Sin(x)), "((Sin(x) ^ -1) * Cos(x))");      // d(ln(sin(x)))/dx = cot(x)
            TestDerivative(Tan(x.Cubed().Plus(8.Times(x))),               // d(tan(x³+8x))/dx = sec²(x³+8x)*(3x²+8)
                "((Cos(((x ^ 3) + (x * 8))) ^ -2) * (((x ^ 2) * 3) + 8))");
        }

        public static void TestAll()
        {
            TestSimplifications();
            TestCompoundExpression();
            TestTrigonometricExpression();
            TestFunctionDerivatives();
            TestPolynomialDerivative();
            TestChainRule();
        }
    }
}

Next time, by overwhelming public demand (one request has flooded in), I'll be exercising these extension methods in a simple differential graphing application.

Acknowledgements
  • Thanks are due to my colleague Stuart Stein for his code review with useful comments, particularly the suggestion that I should embrace the expression body syntax.
  • Thanks also to Pavel Vladov for the free use of his C# Syntax Highlighter to format the code snippets in this article.

Thursday, 21 March 2019

Differentiating f(x)^g(x)

Black Pixel Red Pixel

My favourite mathematical "YouTutor" is Steve Chow, the guy behind the superlative blackpenredpen collection of mathematics videos on YouTube. Steve has developed a way of explaining proofs using two or more pen colours, hence the name of his YouTube channel. He holds multiple marker pens in one hand while writing on a whiteboard, deftly switching colours so as to highlight the essential differences between each pair of lines in the mathematical exposition. It's a very effective expository technique.

Until recently, my number one blackpenredpen video was this one, where Steve uses logarithmic differentiation and the chain rule to prove the power, product, and quotient rules of derivatives. The reason it's so high on my list is because of the additional bonus content at the end, kind of a "hidden track" type deal, where he goes on to solve for the derivative of a function raised to the power of another function.



Just recently Steve, a keen runner of physical marathons, published a truly marathon mathematical video - my new favourite - in which he performs a hundred and one integration proofs in a single take. I'll let you search for that 6 hour video yourself, and if you feel like sitting through it, may I suggest you limit your accompanying playlist to just these three songs: Steely Dan's "My Old School", Supertramp's "School", and Bowling For Soup's "High School Never Ends". This will help you digest the integration marathon in healthy 15 minute chunks.

Getting back to the bonus differentiation, I remember thinking, this procedure has the effect of subtracting one from the power of one of the functions. Yet it doesn't use the power rule in its own derivation. Aha, I bet this could be used to derive the power rule itself! And since I was looking for a suitable subject for my first foray into maths typesetting, using LaTeX in the Overleaf environment, that became my target.

Here is the resultant PDF. You'll notice the homage to blackpenredpen in some of the text colouring.



The actual Overleaf project which produced this PDF can be viewed here, if you're interested in examining the source (LaTeX is a markup language, just as HTML is a markup language).

This is the first in a series of articles looking at aspects of differentiation. Next time I'll be presenting a neat & simple differentiator, utilising the Expression Trees feature of Microsoft C#.

Saturday, 16 March 2019

Kedgeree

My Curry Here

Smoked haddock kedgeree, with added petits pois and raisins why. Because they seem to go nicely. This recipe has been refined for more than a quarter century (I made it once in 1993, then suddenly again today). Wife says it's terrific, and she's a fantastic cook, so there's always that.

Preparation time: 15 minutes.
Cook time: 40 minutes.
Serves 4.
Main Ingredients:
  • 300g smoked fillet of haddock
  • 4 eggs
  • 100g frozen petits pois
  • 300 ml milk
  • bunch of coriander
  • 2 bay leaves
Rice Ingredients:
  • 300g long grain easy cook rice (well rinsed)
  • 1 large onion (finely chopped)
  • 30g raisins
  • 1 tsp ground coriander
  • 1 tsp ground turmeric
  • 2 tsp roasted curry powder
  • 1 chicken or vegetable stock cube
  • 2 tbsp rapeseed oil
Method:
1. Very important: partition your ingredients into rice and non-rice work areas!
2. Fry the onion in rapeseed oil until translucent (about 5 minutes).
3. Add the spices & salt seasoning; stir and continue frying for a couple of minutes.
4. Add the rice, raisins, crumbled stock cube, and 300ml of water.
5. Cover and simmer for 5 minutes, while...
6. ...poaching the haddock and bay in 300ml of milk, until flaky.
7. Drain the milk from the fish pan into the rice, then roughly flake the fish.
8. Cover the rice pan and remove from the heat.
9. Add the petits pois to the fish (they'll defrost almost instantly).
10. Boil the eggs for 4½ minutes, then plunge them into cold water.
11. After 5 minutes, peel and quarter the eggs, then add to the fish along with the rice and coriander.
12. Mix gently, season with black pepper, then serve with any additional herbage
(e.g. parsley, if you like that sort of thing).

Friday, 1 March 2019

Smiddy Shank Soup

A Photo Essay

A very hearty and rustic Ham Hough Red Lentil Soup. Seriously, if they had this recipe in Zelda: Breath of the Wild, it would break the cookery system. It brings health, stamina, attack and defence boosts, while simultaneously providing resistance to heat, cold and electricity. One plate of it would fetch at least two hundred and eighty rupees. Pro tip: for an additional speed boost, replace the carrots with "swift carrots".

Preparation time: 20 minutes.
Cook time: 3 hours 15 minutes.
Serves 8.
Ingredients:
  • 800g Blair Drummond Smiddy smoked ham hough (shank)
  • 2 baking potatoes
  • 1 large onion
  • 3 celery stalks with leaves
  • 6 carrots
  • 250g red split lentils (well rinsed)
  • 4 garlic cloves
  • 2 stock cubes (vegetable or chicken)
  • 6-8 black peppercorns
  • 1 teaspoon dried oregano
  • 1 bouquet garni
  • 2-3 bay leaves
Method:
1. To 3 litres of boiling water in a pot, add the ham, bouquet garni, peppercorns and bay.
2. Boil for 2½ hours.
3. Meanwhile, back at the vegetable station...
...peel the potatoes, onion, celery, carrots & garlic.
4. Dice the potatoes, onion, celery, and four of the carrots.
5.Crush the garlic and grate the remaining two carrots.
6. At the 2½ hour mark, remove the ham from the pot.
7. Skim the water, then add the potatoes...
...the chopped carrots...
...the celery...
...the onion...
...the grated carrots...
...the crushed garlic...
...the oregano and stock cubes...
...and finally the rinsed lentils. Simmer for 45 minutes.
8. Strip the ham off the shank and add back into the soup before serving.
9. It sticks tae yer ribs. Serve with a twist of pepper.
Acknowledgements

I am gratefully indebted to my beautiful assistant and arm model @catsjamas for her help with the visual production aspects of this project.

Monday, 18 June 2018

Japanese Numbers

The New Job

Nowadays I'm working in medical research. One recent task was to internationalise some UI, then localise it into Japanese. With little in the way of a clue how to do this, but having built up some Audible audiobook credits over the spring months, I decided to take advantage of my three hour daily commute to learn a little about the Japanese language and culture. After all, there's every chance I'll soon have the opportunity to visit Japan on the company ticket, and I'd like to be able to speak to people in their first language, without causing humongous offence. So I bought a couple of Teach Yourself Japanese audiobooks.

Linguistic Regularity

Boy, did they do a job on me. You would never think, looking at the Japanese writing system, reputedly one of the most horrendously complex in the world, that this language could be a clean and tidy, logically lovable, stone fox. But I was soon smitten. I mean, just look at some of these features. There are no singular or plural nouns, there's just the noun. The verb to be, one of the most irregular in my native tongue, doesn't conjugate at all; it's just like I be, you be, he she or it be. Month names are like first month, second month, and so on. Dates are in perfect ANSI order: the greater subdivision always comes first, so it's year-month-day, or yyyy-mm-dd. Ideal for sorting by computer! There are only 47 syllables in the language, and every vowel is always pronounced the same.

Of course that's just the first flush of romance; eventually reality and entropy make themselves heard. There might indeed be so many welcome regularities in Japanese, but then they go and spoil it all by saying something stupid like "there are three levels of formality for every sentence", or "men and women speak with subtly different inflections." This last, together with the fact that most teachers of Japanese are female, has the interesting consequence that male learners usually, unwittingly, start off sounding quite effeminate.

Numeric Logicality

But then there are the numbers. The Japanese counting system is beautifully simple and regular, with very few exceptions in structure and pronunciation. Again, the first idealistic impression is somewhat spoiled by reality, this time by the introduction of "counters". These are little particles which must be appended to numbers, when counting objects. Surprisingly, the particular counter used is often based upon the physical appearance of the thing being counted. Then again, certain numbers cannot be used in counting certain types of thing, but must be substituted with other special purpose counters...

Still, the basic counting system is quite easy to learn, interesting too, and if nothing else can at least be used while you count through the steps of kata in your chosen martial art (mine is press-ups-do). I considered it potentially a big win, both mathematically and scientifically, to count as high as I possibly could using this system. As a study aid I made a little text file, delighted to realise that all the kanji characters and superscript numerals available in unicode were supported in Notepad. Printed out and laminated in A6, it became a pocket reference card. Later in A4 format, it would become a placemat suitable, used in conjunction with a transparent plate, for learning your ichi-ni-san while eating your sushi.

Original Text Version

The version below can be copied to the clipboard as text, then saved in Notepad. Just remember to select one of the Unicode options, such as UTF-8, from the Encoding drop-down list.

0  〇/零  zero/rei  JAPANESE NUMBERS REFERENCE CARD  * Learn these exceptions!

       1            10               100              1,000           10,000
1   一 ichi       十 jū            百 hyaku         千 sen        一万 ichi-man
2   二 ni       二十 ni-jū       二百 ni-hyaku     二千 ni-sen     二万 ni-man
3   三 san      三十 san-jū      三百 san-byaku *  三千 san-zen *  三万 san-man
4   四 yon      四十 yon-jū      四百 yon-hyaku    四千 yon-sen    四万 yon-man
5   五 go       五十 go-jū       五百 go-hyaku     五千 go-sen     五万 go-man
6   六 roku     六十 roku-jū     六百 ro-ppyaku *  六千 roku-sen   六万 roku-man
7   七 nana     七十 nana-jū     七百 nana-hyaku   七千 nana-sen   七万 nana-man
8   八 hachi    八十 hachi-jū    八百 ha-ppyaku *  八千 ha-ssen *  八万 hachi-man
9   九 kyū      九十 kyū-jū      九百 kyū-hyaku    九千 kyū-sen    九万 kyū-man

10⁵ 十万 jū-man       10⁶ 百万 hyaku-man       10⁷ 千万 sen-man

10⁸    10¹²   10¹⁶   10²⁰   10²⁴     10²⁸   10³²   10³⁶   10⁴⁰   10⁴⁴   10⁴⁸
億     兆     京      垓     𥝱/秭    穣     溝      澗     正     載     極
oku    chō    kei    gai    jo/shi   jō     kō     kan    sei    sai   goku

10¹² 一兆 itchō      8x10¹² 八兆 hatchō        10¹³ 十兆 jutchō †
10¹⁶ 一京 ikkei      6x10¹⁶ 六京 rokkei      8x10¹⁶ 八京 hakkei

10¹⁷ 十京 jukkei †     † Multiples of 10:  change -jū to -jutchō or -jukkei
10¹⁸ 百京 hyakkei ‡    ‡ Multiples of 100: change -ku to -kkei

10⁵²⸍⁵⁶       10⁵⁶⸍⁶⁴       10⁶⁰⸍⁷²       10⁶⁴⸍⁸⁰       10⁶⁸⸍⁸⁸
恒河沙         阿僧祇     那由他/那由多     不可思議      無量大数
gōgasha       asōgi        nayuta        fukashigi    muryōtaisū

How It Works

The first thing to do is memorise the top left column, digits 1 to 9, and 10: ichi, ni, san, yon, go, roku, nana, hachi, kyū, jū. Note that the horizontal bar along the top of a vowel just makes it sound "long" - some guides use double vowels for this purpose.

Next, to form the numbers 11 to 19, just use literally ten-one, ten-two, ten-three, etc. So that's jū-ichi, jū-ni, jū-san, etc.

Multiples of ten up to ninety are found in the second column: ni-jū (20) is literally two tens, san-jū (30) three tens, and so on. Intermediate numbers are simply concatenated or added together, so 42 is "four tens two", or yon-jū-ni.

It's a similar story for multiples of 100, found in the third column, although there are three little exceptions to be aware of - more or less subtle variations in pronunciation. I've marked these with asterisks. So although 3 is san and 100 is hyaku, 300 is not san-hyaku but san-byaku; note the hard 'b' sound here, and also the hard 'pp' in 600 (ro-ppyaku instead of roku-hyaku) and 800 (ha-ppyaku rather than hachi-hyaku).

Multiples of 1,000, as seen in the fourth column, also have exceptions marked at 3,000 and 8,000, although there isn't an exception for 6,000, which is just business as usual: roku-sen.

Chinese Roots

Reaching 10,000 there's another little surprise in store, and quite an important one. It might help a bit to think first about English powers of 10. Notice that 100 (one hundred) and 1,000 (one thousand) both begin with the word "one", whereas 10 (ten) and 10,000 (ten thousand) don't. But actually, there is an old English number "myriad", one of whose older meanings is literally ten thousand. We might therefore think of 10,000 as "one myriad", and in fact that's similar to what the Japanese, following on from Chinese tradition, actually do. Ten thousand is best translated into Japanese via "one myriad", hence ichi-man.

With that in mind, column five otherwise proceeds much as before: 20,000 is ni-man (two myriad), 30,000 is san-man, and so on without further exceptions. But we have just passed an important milestone here. Just as our big numbers tend to group themselves into powers of 1,000 such as thousands, millions, billions, trillions, and so on, each time multiplying by 10³ and adding another comma, so the oriental tradition is to partition into powers of 10,000, and groups of four digits. The importance of this structure is reflected in the use of the ichi (one) prefix for 10,000: not just man (myriad), but ichi-man (one myriad).

The same caveat will apply to all subsequent powers of 10,000. But before getting to those, look at the next few powers which are not multiples of 4, and don't need the ichi prefix. As shown immediately below the main table, we already knew how to construct:
  • one hundred thousand (100,000  = 10⁵ = ten myriad = jū-man),
  • one million (1,000,000 = 10⁶ = hundred myriad = hyaku-man), and
  • ten million (10,000,000 = 10⁷ = thousand myriad = sen-man).
Next, we shall alight upon the Big Sequence, the powers of 10⁴ spanning 10⁸ to 10⁴⁸ (one quindecillion!), listed in the second horizontal table.

To Infinity And Beyond (almost)

So, now we avail ourselves of the long line of Powers of Myriad (that's definitely the name of my next witching and wizarding trilogy, watch out JK). We have reached 10⁸ which is labelled oku, but since 8 is a multiple of 4, we have to prepend this with ichi. Thus, 10⁸, or a hundred million, becomes ichi-oku. Two hundred million is ni-oku. And so on.

Things chug along quite nicely until we reach one trillion, or 10¹². Everything we've covered so far leads to the expectation ichi-chō, until we stop and consider the clumsiness of that utterance. Surely that is the reason the wise Japanese have decided to contract it to itchō. Similarly, eight trillion becomes not hachi-chō but hatchō, while ten trillion is not jū-chō but jutchō (with the shorter u).

Similar exceptions apply to the next big multiplier, 10¹⁶ (which becomes not ichi-kei but ikkei), as well as multiples of this by 6, 8, 10 and 100. There is a pattern to the modifications in these cases, codified by the daggered comments in the bottom half of the reference card (or placemat).

Finally, we reach the bottom 3 rows of the card, where the only reliable historical sources contradict each other. In a denouement reminiscent of the war between American and British billions (see Wikipedia for details, https://en.wikipedia.org/wiki/Names_of_large_numbers), there's a forking schism in which the same terms are used to count, simultaneously and incompatibly, using multiples of 10⁴ and/or 10⁸, into an uncaring infinity. Of limited scientific use admittedly, but why not take advantage of the Creative Commons licence, and print out and laminate your own set of Japanese Numbers placemats?

Thursday, 17 November 2016

My Collie Here

Meet our Archie

Every cat needs a collie, and this is Cleo's new wee brother, Archie.

"Wee" brother is right. Although he's grown to double her size, in the 4½ weeks he's been with us - compare the picture on the right with the one below, taken 4½ weeks later - he still gets "Wee Archie" due to his determination to cover the whole house in a uniform coat of wee!

Actually, that's a bit unfair. When he came to us at exactly 8 weeks old, he was already trained to go on the puppy training pads, and had more or less a 100% record on those. But that was only for wee; anything more demanding would be deposited instantly and immediately, just wherever he happened to be, as the mood took him. Now that we've started trying to take away the training pads, the result - with apologies to Kevin Bacon - has been a return to Everything Everywhere.

End of the Tunnel

Yesterday he attended his first puppy socialising class at The Dogs Trust, where he did us proud with his precocious knowledge of the Sit, Stay, and Down commands. All credit to Linda there, for her decades of experience training very intelligent border collies.

Socialising per se was less successful. Archie has had to be kept away from other dogs for longer than is usual, because his breeder failed to immunise him against deadly diseases such as Parvovirus. We didn't find out about this until we took Archie to our own vet, for what we thought would be his second and final inoculation, but which turned out to be a complete restart at day one.

So yesterday, two weeks later than he should have, he finally got to meet other pups. He was timid and reluctant, and basically failed to do so. However the instructor expressed his confidence that this delayed development can be put right over the coming month.

Archie and Cleo, everybody!

Thursday, 10 November 2016

FLAC Forensics

Frequency Fingerprints

I was asked by a friend if it was possible to check whether some downloaded files were actually faithful to the original WAV format entities, or whether they had in fact been rehydrated from an unfortunate, lossy intermediate MP3 excursion. The files in question were FLAC compressions of the original 3½ hour, 42 track Analord series of electronic pieces by Richard D. James, most of which were only ever issued on 12" vinyl by the artist's now defunct Rephlex Records outlet.

Update: My personal RD James expert informs me that these original 42 tracks were also released in lossless digital format on Rephlex, and that some (or all?) of the subsequently released additional, digital-only tracks were also added to these by Rephlex.

I loaded the first of these files, SteppingFilter 101, into Audacity, and took a look at the frequency domain graph (Analyze|Plot Spectrum...). After a brief complaint about only being able to analyze 237.8 seconds of audio at a time, the result was this:


Notice the slight uptick at the extreme right (high frequency) end, around 22kHz. This represents extraneous noise generated by the digital sampling process, which in this case appears to have been set naturally enough to the CD standard stereo setting of 44.1kHz. This effect will be present in any rip, at some frequency or other, and is of a different kind from the artefacts introduced by MP3 processing.

Next, I used FooBar2000/LAME to convert this file to MP3 format, using the highest available quality, Constant Bit Rate standard preset, namely 320kbps CBR (actually LAME can handle non-ISO bit rates of up to 640kbps via its freeformat option, but very few MP3 players can handle such files).

The result of this 320kbps conversion has a very obvious steep cutoff at about 20kHz:


Any attempt to convert this back to the WAV format will preserve this telltale high frequency cutoff. Does this mean we can be confident that the source file represents a good, high quality, uncompressed rip from the original vinyl? I would say confident, yes; certain, well, that's another bottle of kippers. We haven't ruled out nonstandard MP3 or other shenanigans with this test alone, but those compression antics are at worst extremely unlikely.