source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Problems.GaussianProcessTuning/Interpreter.cs @ 9155

using System;
using System.Linq;
using System.Text;
using System.Threading;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Problems.GaussianProcessTuning {
  [StorableClass]
  [Item("Interpreter", "An interpreter for Gaussian process configurations represented as trees.")]
  public class Interpreter : Item {
    [StorableConstructor]
    protected Interpreter(bool deserializing) : base(deserializing) { }
    protected Interpreter(Interpreter original, Cloner cloner)
      : base(original, cloner) {
    }
    public Interpreter()
      : base() { }
    public override IDeepCloneable Clone(Cloner cloner) {
      return new Interpreter(this, cloner);
    }

    public void EvaluateGaussianProcessConfiguration(ISymbolicExpressionTree tree, IRegressionProblemData problemData, out double negLogLikelihood, out IGaussianProcessSolution solution) {
      var meanFunction = GetMeanFunction(tree);
      var covFunction = GetCovFunction(tree);

      var gprAlg = new GaussianProcessRegression();
      gprAlg.Problem.ProblemDataParameter.Value = problemData;
      gprAlg.CovarianceFunction = covFunction;
      gprAlg.MeanFunction = meanFunction;
      gprAlg.GaussianProcessModelCreatorParameter.Value =
        gprAlg.GaussianProcessModelCreatorParameter.ValidValues.First(
          v => v is TunedGaussianProcessRegressionModelCreator);
      gprAlg.MinimizationIterations = 50;

      var signal = new AutoResetEvent(false);
      double result = double.MaxValue;
      IGaussianProcessSolution regSolution = null;
      gprAlg.Stopped += (sender, args) => {
        result = ((DoubleValue)gprAlg.Results["NegativeLogLikelihood"].Value).Value;
        if (gprAlg.Results.ContainsKey("Solution"))
          regSolution = (IGaussianProcessSolution)gprAlg.Results["Solution"].Value;
        signal.Set();
      };
      Exception ex = null;
      gprAlg.ExceptionOccurred += (sender, args) => {
        result = double.MaxValue;
        regSolution = null;
        ex = args.Value;
        signal.Set();
      };

      gprAlg.Prepare();
      gprAlg.Start();

      signal.WaitOne();
      if (ex != null) throw ex;

      gprAlg.Prepare();
      gprAlg.Problem = null;
      solution = regSolution;
      negLogLikelihood = result;
    }

    /*
    public void EvaluateGaussianProcessConfiguration(ISymbolicExpressionTree tree, IRegressionProblemData trainingData, Dataset testData, IEnumerable<int> testRows,
      out double[] means, out double[] variances) {
      string meanExpression, meanHyperParameter;
      string covExpression, covHyperParameter;
      string likFunction, likHyperParameter;
      GetMeanFunction(tree, trainingData.AllowedInputVariables.Count(), out meanExpression, out meanHyperParameter);
      GetCovFunction(tree, trainingData.AllowedInputVariables.Count(), out covExpression, out covHyperParameter);
      GetLikelihoodFunction(tree, out likFunction, out likHyperParameter);

      double[,] y = new double[trainingData.TrainingIndices.Count(), 1];
      double[,] yImg = new double[trainingData.TrainingIndices.Count(), 1];

      int r, c;
      r = 0;
      foreach (var e in trainingData.Dataset.GetDoubleValues(trainingData.TargetVariable, trainingData.TrainingIndices)) {
        y[r++, 0] = e;
      }
      double[,] x = new double[y.Length, trainingData.AllowedInputVariables.Count()];
      double[,] xImg = new double[y.Length, trainingData.AllowedInputVariables.Count()];
      c = 0;
      foreach (var allowedInput in trainingData.AllowedInputVariables) {
        r = 0;
        foreach (var e in trainingData.Dataset.GetDoubleValues(allowedInput, trainingData.TrainingIndices)) {
          x[r++, c] = e;
        }
        c++;
      }

      double[,] xTest = new double[testRows.Count(), trainingData.AllowedInputVariables.Count()];
      double[,] xTestImg = new double[testRows.Count(), trainingData.AllowedInputVariables.Count()];
      c = 0;
      foreach (var allowedInput in trainingData.AllowedInputVariables) {
        r = 0;
        foreach (var e in testData.GetDoubleValues(allowedInput, testRows)) {
          xTest[r++, c] = e;
        }
        c++;
      }

      object oldX = null;
      try {
        oldX = MLApp.GetVariable("x", "base");
      }
      catch {
      }
      if (oldX == null || oldX is Missing || ((double[,])oldX).Length != x.Length) {
        MLApp.PutFullMatrix("y", "base", y, yImg);
        MLApp.PutFullMatrix("x", "base", x, xImg);
      }
      MLApp.PutFullMatrix("xTest", "base", xTest, xTestImg);
      ExecuteMatlab("hyp0 = " + GetHyperParameterString(meanHyperParameter, covHyperParameter, likHyperParameter) + ";");
      ExecuteMatlab("infFun =  " + GetInferenceMethodString(tree) + ";");
      ExecuteMatlab("meanExpr = " + meanExpression + ";");
      ExecuteMatlab("covExpr = " + covExpression + ";");
      ExecuteMatlab("likExp = " + likFunction + ";");
      ExecuteMatlab(
        "try " +
        "  hyp = minimize(hyp0,'gp', -50, infFun, meanExpr, covExpr, likExp, x, y);" +
        "  [ymu, ys2, fmu, fs2] = gp(hyp, infFun, meanExpr, covExpr, likExp, x, y, xTest); " +
        "catch " +
        "  ymu = zeros(size(xTest, 1), 1); " +
        "  ys2 = ones(size(xTest, 1), 1); " +
        "end");
      var meansMat = (double[,])MLApp.GetVariable("ymu", "base");
      var variancesMat = (double[,])MLApp.GetVariable("ys2", "base");
      means = new double[meansMat.GetLength(0)];
      for (int i = 0; i < means.Length; i++)
        means[i] = meansMat[i, 0];
      variances = new double[variancesMat.GetLength(0)];
      for (int i = 0; i < variances.Length; i++)
        variances[i] = variancesMat[i, 0];
    }
    */
    private string GetInferenceMethodString(ISymbolicExpressionTree tree) {
      return "'infEP'";
    }

    private IMeanFunction GetMeanFunction(ISymbolicExpressionTree tree) {
      return GetMeanFunction(tree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(0));
    }

    private ICovarianceFunction GetCovFunction(ISymbolicExpressionTree tree) {
      return GetCovFunction(tree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(1));
    }


    private void GetLikelihoodFunction(ISymbolicExpressionTree tree, out string expression, out string hyperParameter) {
      var expressionBuilder = new StringBuilder();
      var hyperParameterBuilder = new StringBuilder();
      hyperParameterBuilder.Append("[");
      GetLikelihoodFunction(tree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(2), expressionBuilder,
                                   hyperParameterBuilder);
      hyperParameterBuilder.Append("]");
      expression = expressionBuilder.ToString();
      hyperParameter = hyperParameterBuilder.ToString();
    }

    private void GetLikelihoodFunction(ISymbolicExpressionTreeNode node, StringBuilder expressionBuilder, StringBuilder hyperParameterBuilder) {
      if (node.Symbol is LikGauss) {
        var likNode = node as LikGaussTreeNode;
        expressionBuilder.Append("'likGauss'");
        hyperParameterBuilder.Append(likNode.Sigma);
      } else {
        throw new ArgumentException("unknown likelihood function " + node.Symbol);
      }
    }

    private string GetHyperParameterString(string meanHyp, string covHyp, string likHyp) {
      return "struct('mean', " + meanHyp +
             ", 'cov', " + covHyp +
             ", 'lik', " + likHyp +
             ")";
      ;
    }


    private ICovarianceFunction GetCovFunction(ISymbolicExpressionTreeNode node) {
      if (node.Symbol is CovConst) {
        return new CovarianceConst();
      } else if (node.Symbol is CovScale) {
        return new CovarianceScale();
      } else if (node.Symbol is CovMask) {
        var maskNode = node as CovMaskTreeNode;
        var covSymbol = node.Symbol as CovMask;
        var cov = new CovarianceMask();
        cov.SelectedDimensionsParameter.Value = new IntArray((from i in Enumerable.Range(0, covSymbol.Dimension)
                                                              where maskNode.Mask[i]
                                                              select i).ToArray());
        return cov;
      } else if (node.Symbol is CovLin) {
        return new CovarianceLinear();
      } else if (node.Symbol is CovLinArd) {
        return new CovarianceLinearArd();
      } else if (node.Symbol is CovMatern) {
        var covSymbol = node.Symbol as CovMatern;
        var cov = new CovarianceMaternIso();
        cov.DParameter.Value = cov.DParameter.ValidValues.Single(x => x.Value == covSymbol.D);
        return cov;
      } else if (node.Symbol is CovSeArd) {
        return new CovarianceSquaredExponentialArd();
      } else if (node.Symbol is CovSeIso) {
        return new CovarianceSquaredExponentialIso();
      } else if (node.Symbol is CovRQIso) {
        return new CovarianceRationalQuadraticIso();
      } else if (node.Symbol is CovRQArd) {
        return new CovarianceRationalQuadraticArd();
      } else if (node.Symbol is CovPeriodic) {
        return new CovariancePeriodic();
      } else if (node.Symbol is CovNoise) {
        return new CovarianceNoise();
      } else if (node.Symbol is CovSum) {
        var covSum = new Algorithms.DataAnalysis.CovarianceSum();
        covSum.Terms.Add(GetCovFunction(node.GetSubtree(0)));
        foreach (var subTree in node.Subtrees.Skip(1)) {
          covSum.Terms.Add(GetCovFunction(subTree));
        }
        return covSum;
      } else if (node.Symbol is CovProd) {
        var covProd = new Algorithms.DataAnalysis.CovarianceProduct();
        covProd.Factors.Add(GetCovFunction(node.GetSubtree(0)));
        foreach (var subTree in node.Subtrees.Skip(1)) {
          covProd.Factors.Add(GetCovFunction(subTree));
        }
        return covProd;
      } else {
        throw new ArgumentException("unknown symbol " + node.Symbol);
      }
    }


    private IMeanFunction GetMeanFunction(ISymbolicExpressionTreeNode node) {
      if (node.Symbol is MeanConst) {
        return new Algorithms.DataAnalysis.MeanConst();
      } else if (node.Symbol is MeanLinear) {
        return new Algorithms.DataAnalysis.MeanLinear();
      } else if (node.Symbol is MeanProd) {
        var meanProd = new Algorithms.DataAnalysis.MeanProduct();
        meanProd.Factors.Add(GetMeanFunction(node.GetSubtree(0)));
        foreach (var subTree in node.Subtrees.Skip(1)) {
          meanProd.Factors.Add(GetMeanFunction(subTree));
        }
        return meanProd;
      } else if (node.Symbol is MeanSum) {
        var meanSum = new Algorithms.DataAnalysis.MeanSum();
        meanSum.Terms.Add(GetMeanFunction(node.GetSubtree(0)));
        foreach (var subTree in node.Subtrees.Skip(1)) {
          meanSum.Terms.Add(GetMeanFunction(subTree));
        }
        return meanSum;
      } else if (node.Symbol is MeanZero) {
        return new Algorithms.DataAnalysis.MeanZero();
      } else {
        throw new ArgumentException("Unknown mean function" + node.Symbol);
      }
    }

    //private bool[] CalculateMask(ISymbolicExpressionTreeNode node) {
    //  var maskNode = node as MeanMaskTreeNode;
    //  if (maskNode != null) {
    //    bool[] newMask = CombineMasksProd(maskNode.Mask, CalculateMask(node.GetSubtree(0)));
    //    return newMask;
    //  } else if (node.Symbol is MeanProd) {
    //    bool[] newMask = CalculateMask(node.GetSubtree(0));
    //    foreach (var subTree in node.Subtrees.Skip(1)) {
    //      newMask = CombineMasksProd(newMask, CalculateMask(subTree));
    //    }
    //    return newMask;
    //  } else if (node.Symbol is MeanSum) {
    //    bool[] newMask = CalculateMask(node.GetSubtree(0));
    //    foreach (var subTree in node.Subtrees.Skip(1)) {
    //      newMask = CombineMasksSum(newMask, CalculateMask(subTree));
    //    }
    //    return newMask;
    //  } else if (node.SubtreeCount == 1) {
    //    return CalculateMask(node.GetSubtree(0));
    //  } else if (node is SymbolicExpressionTreeTerminalNode) {
    //    return null;
    //  } else {
    //    throw new NotImplementedException();
    //  }
    //}

    //private bool[] CombineMasksProd(bool[] m, bool[] n) {
    //  if (m == null) return n;
    //  if (n == null) return m;
    //  if (m.Length != n.Length) throw new ArgumentException();
    //  bool[] res = new bool[m.Length];
    //  for (int i = 0; i < res.Length; i++)
    //    res[i] = m[i] | n[i];
    //  return res;
    //}


    //private bool[] CombineMasksSum(bool[] m, bool[] n) {
    //  if (m == null) return n;
    //  if (n == null) return m;
    //  if (m.Length != n.Length) throw new ArgumentException();
    //  bool[] res = new bool[m.Length];
    //  for (int i = 0; i < res.Length; i++)
    //    res[i] = m[i] & n[i];
    //  return res;
    //}

    private string ToVectorString(bool[] b) {
      var strBuilder = new StringBuilder();
      strBuilder.Append("[");
      if (b.Length == 1) // workaround for bug in GPML
      {
        if (!b[0]) strBuilder.Append("1");
      } else {
        for (int i = 0; i < b.Length; i++) {
          if (i > 0) strBuilder.Append(", ");
          strBuilder.Append(b[i] ? "0" : "1");
        }
      }
      strBuilder.Append("]");
      return strBuilder.ToString();
    }
    private string ToVectorString(double[] xs, bool[] mask) {
      if (xs.Length != mask.Length) throw new ArgumentException();
      var strBuilder = new StringBuilder();
      strBuilder.Append("[");
      for (int i = 0; i < xs.Length; i++)
        if (!mask[i]) {
          if (i > 0) strBuilder.Append("; ");
          strBuilder.Append(xs[i]);
        }
      strBuilder.Append("]");
      return strBuilder.ToString();
    }
  }
}