source: branches/2994-AutoDiffForIntervals/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis-3.4/AutoDiffInterpreterTest.cs @ 17308

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Problems.DataAnalysis.Symbolic;
using HeuristicLab.Random;
using Microsoft.VisualStudio.TestTools.UnitTesting;

namespace HeuristicLab.Problems.DataAnalysis.Tests {
  [TestClass]
  public class AutoDiffInterpreterTest {

    [TestMethod]
    [TestCategory("Problems.DataAnalysis")]
    [TestProperty("Time", "short")]
    public void TestAutoDiffUsingNumericDifferences() {

      // create random trees and evaluate on random data
      // calc gradient for all parameters
      // use numeric differences for approximate gradient calculation
      // compare gradients

      var grammar = new TypeCoherentExpressionGrammar();
      grammar.ConfigureAsDefaultRegressionGrammar();
      // activate supported symbols
      grammar.Symbols.First(s => s is Square).Enabled = true;
      grammar.Symbols.First(s => s is SquareRoot).Enabled = true;
      grammar.Symbols.First(s => s is Cube).Enabled = true;
      grammar.Symbols.First(s => s is CubeRoot).Enabled = true;
      grammar.Symbols.First(s => s is Sine).Enabled = true;
      grammar.Symbols.First(s => s is Cosine).Enabled = true;
      grammar.Symbols.First(s => s is Exponential).Enabled = true;
      grammar.Symbols.First(s => s is Logarithm).Enabled = true;
      grammar.Symbols.First(s => s is Absolute).Enabled = false; // XXX not yet supported by old interval calculator
      grammar.Symbols.First(s => s is AnalyticQuotient).Enabled = false; // not yet supported by old interval calculator

      var varSy = (Variable)grammar.Symbols.First(s => s is Variable);
      varSy.AllVariableNames = new string[] { "x", "y" };
      varSy.VariableNames = varSy.AllVariableNames;
      varSy.WeightMu = 1.0;
      varSy.WeightSigma = 1.0;
      var rand = new FastRandom(1234);

      // random data
      var values = new double[100, 2];
      for (int i = 0; i < 100; i++)
        for (int j = 0; j < 2; j++) {
          values[i, j] = rand.NextDouble() * 2 - 1;
        }
      var ds = new Dataset(varSy.AllVariableNames, values);
      // buffers
      var fi = new double[100];
      var rows = Enumerable.Range(0, 100).ToArray();

      var eval = new VectorAutoDiffEvaluator();
      var refEval = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();

      var formatter = new InfixExpressionFormatter();
      var sb = new StringBuilder();
      int N = 10000;
      int iter = 0;
      while (iter < N) {
        var t = ProbabilisticTreeCreator.Create(rand, grammar, maxTreeLength: 5, maxTreeDepth: 5);
        var parameterNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();

        var jac = new double[100, parameterNodes.Length];

        eval.Evaluate(t, ds, rows, parameterNodes, fi, jac);

        var refJac = ApproximateGradient(t, ds, rows, parameterNodes, refEval);

        for (int k = 0; k < rows.Length; k++) {
          if (double.IsNaN(fi[k]) || double.IsInfinity(fi[k])) continue; // skip outputs where we expect problematic gradients

          // check partial derivatives 
          for (int p = 0; p < parameterNodes.Length; p++) {
            if (double.IsNaN(jac[k, p]) && double.IsNaN(refJac[k, p])) continue; // both NaN
            if (jac[k, p] == refJac[k, p]) continue; // equal
            if (Math.Abs(jac[k, p]) <= 1e-12 && Math.Abs(refJac[k, p]) <= 1e-12) continue; // both very small

            // check relative error using the larger value as reference
            var refVal = Math.Max(Math.Abs(jac[k, p]), Math.Abs(refJac[k, p]));
            if (Math.Abs(jac[k, p] - refJac[k, p]) > refVal * 1e-4)
              sb.AppendLine($"{jac[k, p]} <> {refJac[k, p]} for {parameterNodes[p]} in {formatter.Format(t)} x={ds.GetDoubleValue("x", k)} y={ds.GetDoubleValue("y", k)}");
          }
        }

        iter++;
      }
      if (sb.Length > 0) {
        Console.WriteLine(sb.ToString());
        Assert.Fail("There were differences when validating AutoDiff using numeric differences");
      }
    }


    [TestMethod]
    [TestCategory("Problems.DataAnalysis")]
    [TestProperty("Time", "short")]
    public void TestVectorAutoDiffInterpreter() {
      var ds = new Dataset(new string[] { "x", "y" }, new double[,] { { 1, 0 }, { 2, 1 } });

      Assert.AreEqual(0.25, CalculateGradient("sqrt(4)", ds)[0]);
      Assert.AreEqual((1.0 / 3.0) * (1.0 / 4.0) , CalculateGradient("cuberoot(8)", ds)[0]);
      Assert.AreEqual((1.0 / 4.0), CalculateGradient("1.0 / 4.0", ds)[0]);
      Assert.AreEqual(-1.0 / 16.0, CalculateGradient("1.0 / 4.0", ds)[1]);
      Assert.AreEqual(1.0 / 16.0, CalculateGradient("1.0 / (-4.0)", ds)[1]);
    }

    #region helper

    private double[] CalculateGradient(string expr, IDataset ds) {
      var eval = new VectorAutoDiffEvaluator();
      var parser = new InfixExpressionParser();

      var rows = new int[1];
      var fi = new double[1];

      var t = parser.Parse(expr);
      var parameterNodes = t.IterateNodesPostfix().Where(n => n.SubtreeCount == 0).ToArray();
      var jac = new double[1, parameterNodes.Length];
      eval.Evaluate(t, ds, rows, parameterNodes, fi, jac);

      var g = new double[parameterNodes.Length];
      for (int i = 0; i < g.Length; i++) g[i] = jac[0, i];
      return g;
    }


    private double[,] ApproximateGradient(ISymbolicExpressionTree t, Dataset ds, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes,
      SymbolicDataAnalysisExpressionTreeLinearInterpreter eval) {
      var jac = new double[rows.Length, parameterNodes.Length];
      for (int p = 0; p < parameterNodes.Length; p++) {

        var x = GetValue(parameterNodes[p]);
        var x_diff = x * 1e-4; // relative change

        // calculate output for increased parameter value
        SetValue(parameterNodes[p], x + x_diff / 2);
        var f = eval.GetSymbolicExpressionTreeValues(t, ds, rows).ToArray();
        for (int i = 0; i < rows.Length; i++) {
          jac[i, p] = f[i];
        }

        // calculate output for decreased parameter value
        SetValue(parameterNodes[p], x - x_diff / 2);
        f = eval.GetSymbolicExpressionTreeValues(t, ds, rows).ToArray();
        for (int i = 0; i < rows.Length; i++) {
          jac[i, p] -= f[i]; // calc difference (and scale for x_diff)
          jac[i, p] /= x_diff;
        }

        // restore original value
        SetValue(parameterNodes[p], x);
      }
      return jac;
    }

    private void SetValue(ISymbolicExpressionTreeNode node, double v) {
      var varNode = node as VariableTreeNode;
      var constNode = node as ConstantTreeNode;
      if (varNode != null) varNode.Weight = v;
      else if (constNode != null) constNode.Value = v;
      else throw new InvalidProgramException();
    }

    private double GetValue(ISymbolicExpressionTreeNode node) {
      var varNode = node as VariableTreeNode;
      var constNode = node as ConstantTreeNode;
      if (varNode != null) return varNode.Weight;
      else if (constNode != null) return constNode.Value;
      throw new InvalidProgramException();
    }
    #endregion
  }
}