source: branches/HeuristicLab.Problems.GrammaticalOptimization/HeuristicLab.Problems.GrammaticalOptimization.SymbReg/SymbolicRegressionProblem.cs @ 11851

using System;
using System.Collections.Generic;
using System.Linq;
using System.Security;
using System.Security.AccessControl;
using System.Text;
using HeuristicLab.Common;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Problems.Instances;
using HeuristicLab.Problems.Instances.DataAnalysis;

namespace HeuristicLab.Problems.GrammaticalOptimization.SymbReg {
  // provides bridge to HL regression problem instances
  public class SymbolicRegressionProblem : IProblem {

    private const string grammarString = @"
        G(E):
        E -> V | V+E | V-E | V*E | V/E | (E) | C | C+E | C-E | C*E | C/E 
        C -> 0..9
        V -> <variables>
        ";
    // C represents Koza-style ERC (the symbol is the index of the ERC), the values are initialized below

    private readonly IGrammar grammar;

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


    public SymbolicRegressionProblem(Random random, string partOfName) {
      var instanceProvider = new RegressionRealWorldInstanceProvider();
      var dds = instanceProvider.GetDataDescriptors().OfType<RegressionDataDescriptor>();

      var problemData = instanceProvider.LoadData(dds.Single(ds => ds.Name.Contains(partOfName)));

      this.N = problemData.TrainingIndices.Count();
      this.d = problemData.AllowedInputVariables.Count();
      if (d > 26) throw new NotSupportedException(); // we only allow single-character terminal symbols so far
      this.x = new double[N][];
      this.y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();

      int i = 0;
      foreach (var r in problemData.TrainingIndices) {
        x[i] = new double[d];
        int j = 0;
        foreach (var inputVariable in problemData.AllowedInputVariables) {
          x[i][j++] = problemData.Dataset.GetDoubleValue(inputVariable, r);
        }
        i++;
      }
      // initialize ERC values
      erc = Enumerable.Range(0, 10).Select(_ => Rand.RandNormal(random) * 10).ToArray();

      char firstVar = 'a';
      char lastVar = Convert.ToChar(Convert.ToByte('a') + d - 1);
      this.grammar = new Grammar(grammarString.Replace("<variables>", firstVar + " .. " + lastVar));
    }


    public double BestKnownQuality(int maxLen) {
      // for now only an upper bound is returned, ideally we have an R² of 1.0
      return 1.0;
    }

    public IGrammar Grammar {
      get { return grammar; }
    }

    public double Evaluate(string sentence) {
      var interpreter = new ExpressionInterpreter();
      return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => interpreter.Interpret(sentence, x[i], erc)));
    }


    public string CanonicalRepresentation(string phrase) {
      return phrase;
      //var terms = phrase.Split('+');
      //return string.Join("+", terms.Select(term => string.Join("", term.Replace("*", "").OrderBy(ch => ch)))
      //  .OrderBy(term => term));
    }

    public IEnumerable<Feature> GetFeatures(string phrase)
    {
      throw new NotImplementedException();
    }


    /*
    public static double OptimizeConstants(string sentence) {

      List<AutoDiff.Variable> variables = new List<AutoDiff.Variable>();
      List<AutoDiff.Variable> parameters = new List<AutoDiff.Variable>();
      List<string> variableNames = new List<string>();

      AutoDiff.Term func;
      if (!TryTransformToAutoDiff(sentence, 0, variables, parameters, out func))
        throw new NotSupportedException("Could not optimize constants of symbolic expression tree due to not supported symbols used in the tree.");
      if (variableNames.Count == 0) return 0.0;

      AutoDiff.IParametricCompiledTerm compiledFunc = AutoDiff.TermUtils.Compile(func, variables.ToArray(), parameters.ToArray());

      List<SymbolicExpressionTreeTerminalNode> terminalNodes = tree.Root.IterateNodesPrefix().OfType<SymbolicExpressionTreeTerminalNode>().ToList();
      double[] c = new double[variables.Count];

      {
        c[0] = 0.0;
        c[1] = 1.0;
        //extract inital constants
        int i = 2;
        foreach (var node in terminalNodes) {
          ConstantTreeNode constantTreeNode = node as ConstantTreeNode;
          VariableTreeNode variableTreeNode = node as VariableTreeNode;
          if (constantTreeNode != null)
            c[i++] = constantTreeNode.Value;
          else if (variableTreeNode != null)
            c[i++] = variableTreeNode.Weight;
        }
      }
      double[] originalConstants = (double[])c.Clone();
      double originalQuality = SymbolicRegressionSingleObjectivePearsonRSquaredEvaluator.Calculate(interpreter, tree, lowerEstimationLimit, upperEstimationLimit, problemData, rows, applyLinearScaling);

      alglib.lsfitstate state;
      alglib.lsfitreport rep;
      int info;

      Dataset ds = problemData.Dataset;
      double[,] x = new double[rows.Count(), variableNames.Count];
      int row = 0;
      foreach (var r in rows) {
        for (int col = 0; col < variableNames.Count; col++) {
          x[row, col] = ds.GetDoubleValue(variableNames[col], r);
        }
        row++;
      }
      double[] y = ds.GetDoubleValues(problemData.TargetVariable, rows).ToArray();
      int n = x.GetLength(0);
      int m = x.GetLength(1);
      int k = c.Length;

      alglib.ndimensional_pfunc function_cx_1_func = CreatePFunc(compiledFunc);
      alglib.ndimensional_pgrad function_cx_1_grad = CreatePGrad(compiledFunc);

      try {
        alglib.lsfitcreatefg(x, y, c, n, m, k, false, out state);
        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
        //alglib.lsfitsetgradientcheck(state, 0.001);
        alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
        alglib.lsfitresults(state, out info, out c, out rep);
      } catch (ArithmeticException) {
        return originalQuality;
      } catch (alglib.alglibexception) {
        return originalQuality;
      }

      //info == -7  => constant optimization failed due to wrong gradient
      if (info != -7) throw new ArgumentException();
    }


    private static alglib.ndimensional_pfunc CreatePFunc(AutoDiff.IParametricCompiledTerm compiledFunc) {
      return (double[] c, double[] x, ref double func, object o) => {
        func = compiledFunc.Evaluate(c, x);
      };
    }

    private static alglib.ndimensional_pgrad CreatePGrad(AutoDiff.IParametricCompiledTerm compiledFunc) {
      return (double[] c, double[] x, ref double func, double[] grad, object o) => {
        var tupel = compiledFunc.Differentiate(c, x);
        func = tupel.Item2;
        Array.Copy(tupel.Item1, grad, grad.Length);
      };
    }

    private static bool TryTransformToAutoDiff(string phrase, int symbolPos, List<AutoDiff.Variable> variables, List<AutoDiff.Variable> parameters, out AutoDiff.Term term)
    {
      var curSy = phrase[0];
      if () {
        var var = new AutoDiff.Variable();
        variables.Add(var);
        term = var;
        return true;
      }
      if (node.Symbol is Variable) {
        var varNode = node as VariableTreeNode;
        var par = new AutoDiff.Variable();
        parameters.Add(par);
        variableNames.Add(varNode.VariableName);
        var w = new AutoDiff.Variable();
        variables.Add(w);
        term = AutoDiff.TermBuilder.Product(w, par);
        return true;
      }
      if (node.Symbol is Addition) {
        List<AutoDiff.Term> terms = new List<Term>();
        foreach (var subTree in node.Subtrees) {
          AutoDiff.Term t;
          if (!TryTransformToAutoDiff(subTree, variables, parameters, variableNames, out t)) {
            term = null;
            return false;
          }
          terms.Add(t);
        }
        term = AutoDiff.TermBuilder.Sum(terms);
        return true;
      }
      if (node.Symbol is Subtraction) {
        List<AutoDiff.Term> terms = new List<Term>();
        for (int i = 0; i < node.SubtreeCount; i++) {
          AutoDiff.Term t;
          if (!TryTransformToAutoDiff(node.GetSubtree(i), variables, parameters, variableNames, out t)) {
            term = null;
            return false;
          }
          if (i > 0) t = -t;
          terms.Add(t);
        }
        term = AutoDiff.TermBuilder.Sum(terms);
        return true;
      }
      if (node.Symbol is Multiplication) {
        AutoDiff.Term a, b;
        if (!TryTransformToAutoDiff(node.GetSubtree(0), variables, parameters, variableNames, out a) ||
          !TryTransformToAutoDiff(node.GetSubtree(1), variables, parameters, variableNames, out b)) {
          term = null;
          return false;
        } else {
          List<AutoDiff.Term> factors = new List<Term>();
          foreach (var subTree in node.Subtrees.Skip(2)) {
            AutoDiff.Term f;
            if (!TryTransformToAutoDiff(subTree, variables, parameters, variableNames, out f)) {
              term = null;
              return false;
            }
            factors.Add(f);
          }
          term = AutoDiff.TermBuilder.Product(a, b, factors.ToArray());
          return true;
        }
      }
      if (node.Symbol is Division) {
        // only works for at least two subtrees
        AutoDiff.Term a, b;
        if (!TryTransformToAutoDiff(node.GetSubtree(0), variables, parameters, variableNames, out a) ||
          !TryTransformToAutoDiff(node.GetSubtree(1), variables, parameters, variableNames, out b)) {
          term = null;
          return false;
        } else {
          List<AutoDiff.Term> factors = new List<Term>();
          foreach (var subTree in node.Subtrees.Skip(2)) {
            AutoDiff.Term f;
            if (!TryTransformToAutoDiff(subTree, variables, parameters, variableNames, out f)) {
              term = null;
              return false;
            }
            factors.Add(1.0 / f);
          }
          term = AutoDiff.TermBuilder.Product(a, 1.0 / b, factors.ToArray());
          return true;
        }
      }
      
      if (node.Symbol is StartSymbol) {
        var alpha = new AutoDiff.Variable();
        var beta = new AutoDiff.Variable();
        variables.Add(beta);
        variables.Add(alpha);
        AutoDiff.Term branchTerm;
        if (TryTransformToAutoDiff(node.GetSubtree(0), variables, parameters, variableNames, out branchTerm)) {
          term = branchTerm * alpha + beta;
          return true;
        } else {
          term = null;
          return false;
        }
      }
      term = null;
      return false;
    }
     */
  }
}