source: branches/HeuristicLab.Problems.GrammaticalOptimization/HeuristicLab.Problems.GrammaticalOptimization/SymbolicRegressionPoly10Problem.cs @ 11659

using System;
using System.Collections.Generic;
using System.Linq;
using System.Security;
using System.Security.AccessControl;
using System.Text;

namespace HeuristicLab.Problems.GrammaticalOptimization {
  public class SymbolicRegressionPoly10Problem : IProblem {
    // length of the longest palindrome in the sentence + number of different symbols occurring in the palindrome
    private const string grammarString = @"
G(E):
E -> V | V+E | V-E | V*E | V/E | (E)
V -> a .. j 
";

    private readonly IGrammar grammar;
    private readonly ExpressionInterpreter interpreter;

    private readonly int N;
    private readonly double[][] x;
    private readonly double[] y;

    public SymbolicRegressionPoly10Problem() {
      this.grammar = new Grammar(grammarString);
      this.interpreter = new ExpressionInterpreter();

      this.N = 500;
      this.x = new double[N][];
      this.y = new double[N];

      GenerateData();
    }

    private void GenerateData() {
      // generate data with fixed seed to make sure that data is always the same
      var rand = new Random(31415);
      for (int i = 0; i < N; i++) {
        x[i] = new double[10];
        for (int j = 0; j < 10; j++) {
          x[i][j] = rand.NextDouble() * 2 - 1;
        }
        // poly-10 no noise
        y[i] = x[i][0] * x[i][1] +
               x[i][2] * x[i][3] +
               x[i][4] * x[i][5] +
               x[i][0] * x[i][6] * x[i][8] +
               x[i][2] * x[i][5] * x[i][9];
      }
    }

    public double GetBestKnownQuality(int maxLen) {
      // for now only an upper bound is returned, ideally we have an R² of 1.0
      // the optimal R² can only be reached for sentences of at least 23 symbols
      return 1.0;
    }

    public IGrammar Grammar {
      get { return grammar; }
    }

    public double Evaluate(string sentence) {
      return RSq(y, Enumerable.Range(0, N).Select(i => interpreter.Interpret(sentence, x[i])).ToArray());
    }

    private double RSq(IEnumerable<double> xs, IEnumerable<double> ys) {
      // two pass implementation, but we don't care
      var meanX = xs.Average();
      var meanY = ys.Average();

      var s = 0.0;
      var ssX = 0.0;
      var ssY = 0.0;
      foreach (var p in xs.Zip(ys, (x, y) => new { x, y })) {
        s += (p.x - meanX) * (p.y - meanY);
        ssX += Math.Pow(p.x - meanX, 2);
        ssY += Math.Pow(p.y - meanY, 2);
      }

      if (s.IsAlmost(0)) return 0;
      if (ssX.IsAlmost(0) || ssY.IsAlmost(0)) return 0;
      return s * s / (ssX * ssY);
    }
  }
}