source: branches/2994-AutoDiffForIntervals/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis.Symbolic-3.4/VectorAutoDiffEvaluatorTest.cs @ 17582

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.Symbolic.Tests {
  [TestClass]
  public class VectorAutoDiffEvaluatorTest {

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

      // 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 = true; 
      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 VectorAutoDiffEvaluatorExamplesTest() {
      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]);
    }


    [TestMethod]
    [TestCategory("Problems.DataAnalysis.Symbolic")]
    [TestProperty("Time", "long")]
    public void VectorEvaluatorsEstimatedValuesConsistencyTest() {
      var twister = new MersenneTwister();
      twister.Seed(31415);
      const int numRows = 100;
      const int Columns = 50;
      const int N = 10000;
      var dataset = Util.CreateRandomDataset(twister, numRows, Columns);

      var grammar = new TypeCoherentExpressionGrammar();
      grammar.ConfigureAsDefaultRegressionGrammar();
      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 Exponential).Enabled = true;
      grammar.Symbols.First(s => s is Logarithm).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 Absolute).Enabled = true;
      grammar.Symbols.First(s => s is AnalyticQuotient).Enabled = true;

      var refInterpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
      var newEvaluator = new VectorEvaluator();
      var newAutoDiffEvaluator = new VectorAutoDiffEvaluator();

      var rows = Enumerable.Range(0, numRows).ToList();
      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 10, 0, 0);
      foreach (ISymbolicExpressionTree tree in randomTrees) {
        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
      }

      for (int i = 0; i < randomTrees.Length; ++i) {
        var tree = randomTrees[i];
        var refValues = refInterpreter.GetSymbolicExpressionTreeValues(tree, dataset, rows).ToArray();
        var newValues = newEvaluator.Evaluate(tree, dataset, rows.ToArray()).ToArray();
        var newAutoDiffValues = new double[numRows];
        newAutoDiffEvaluator.Evaluate(tree, dataset, rows.ToArray(), new ISymbolicExpressionTreeNode[0], newAutoDiffValues, null);

        for (int j = 0; j < rows.Count; j++) {
          if (double.IsNaN(refValues[j]) && double.IsNaN(newValues[j]) && double.IsNaN(newAutoDiffValues[j])) continue;
          string errorMessage = string.Format("Interpreters do not agree on tree {0} {1}.", i, (new InfixExpressionFormatter()).Format(tree));

          var relDelta = Math.Abs(refValues[j]) * 1e-5;
          Assert.AreEqual(refValues[j], newValues[j], relDelta, errorMessage);
          Assert.AreEqual(newValues[j], newAutoDiffValues[j], relDelta, errorMessage);
        }
      }
    }

    #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 ApproximateIntervalGradient(ISymbolicExpressionTree t, Dictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] parameterNodes, IntervalEvaluator eval, out double[] lowerGradient, out double[] upperGradient) {
      lowerGradient = new double[parameterNodes.Length];
      upperGradient = new double[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 r1 = eval.Evaluate(t, intervals);
        lowerGradient[p] = r1.LowerBound;
        upperGradient[p] = r1.UpperBound;

        // calculate output for decreased parameter value
        SetValue(parameterNodes[p], x - x_diff / 2);
        var r2 = eval.Evaluate(t, intervals);
        lowerGradient[p] -= r2.LowerBound;
        upperGradient[p] -= r2.UpperBound;

        lowerGradient[p] /= x_diff;
        upperGradient[p] /= x_diff;

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

    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
  }
}