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

Last change on this file since 11895 was 11895, checked in by gkronber, 7 years ago

#2283: constant opt, expressioncompiler, autodiff, fixes in GP solvers

File size: 8.0 KB
Line 
1using System;
2using System.CodeDom.Compiler;
3using System.Collections.Generic;
4using System.Diagnostics;
5using System.Linq;
6using System.Security;
7using System.Security.AccessControl;
8using System.Security.Authentication.ExtendedProtection.Configuration;
9using System.Text;
10using AutoDiff;
11using HeuristicLab.Common;
12using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
13using HeuristicLab.Problems.DataAnalysis;
14using HeuristicLab.Problems.Instances;
15using HeuristicLab.Problems.Instances.DataAnalysis;
16
17namespace HeuristicLab.Problems.GrammaticalOptimization.SymbReg {
18  // provides bridge to HL regression problem instances
19  public class SymbolicRegressionProblem : ISymbolicExpressionTreeProblem {
20
21    private const string grammarString = @"
22        G(E):
23        E -> V | C | V+E | V-E | V*E | V%E | (E) | C+E | C-E | C*E | C%E
24        C -> 0..9
25        V -> <variables>
26        ";
27    // C represents Koza-style ERC (the symbol is the index of the ERC), the values are initialized below
28
29    // S .. Sum (+), N .. Neg. sum (-), P .. Product (*), D .. Division (%)
30    private const string treeGrammarString = @"
31        G(E):
32        E -> V | C | S | N | P | D
33        S -> EE | EEE
34        N -> EE | EEE
35        P -> EE | EEE
36        D -> EE
37        C -> 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
38        V -> <variables>
39        ";
40
41    // when we use constants optimization we can completely ignore all constants by a simple strategy:
42    // introduce a constant factor for each complete term
43    // introduce a constant offset for each complete expression (including expressions in brackets)
44    // e.g. 1*(2*a + b - 3 + 4) is the same as c0*a + c1*b + c2
45
46    private readonly IGrammar grammar;
47
48    private readonly int N;
49    private readonly double[,] x;
50    private readonly double[] y;
51    private readonly int d;
52    private readonly double[] erc;
53
54
55    public SymbolicRegressionProblem(Random random, string partOfName) {
56      var instanceProvider = new RegressionRealWorldInstanceProvider();
57      var dds = instanceProvider.GetDataDescriptors().OfType<RegressionDataDescriptor>();
58
59      var problemData = instanceProvider.LoadData(dds.Single(ds => ds.Name.Contains(partOfName)));
60
61      this.N = problemData.TrainingIndices.Count();
62      this.d = problemData.AllowedInputVariables.Count();
63      if (d > 26) throw new NotSupportedException(); // we only allow single-character terminal symbols so far
64      this.x = new double[N, d];
65      this.y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
66
67      int i = 0;
68      foreach (var r in problemData.TrainingIndices) {
69        int j = 0;
70        foreach (var inputVariable in problemData.AllowedInputVariables) {
71          x[i, j++] = problemData.Dataset.GetDoubleValue(inputVariable, r);
72        }
73        i++;
74      }
75      // initialize ERC values
76      erc = Enumerable.Range(0, 10).Select(_ => Rand.RandNormal(random) * 10).ToArray();
77
78      char firstVar = 'a';
79      char lastVar = Convert.ToChar(Convert.ToByte('a') + d - 1);
80      this.grammar = new Grammar(grammarString.Replace("<variables>", firstVar + " .. " + lastVar));
81      this.TreeBasedGPGrammar = new Grammar(treeGrammarString.Replace("<variables>", firstVar + " .. " + lastVar));
82    }
83
84
85    public double BestKnownQuality(int maxLen) {
86      // for now only an upper bound is returned, ideally we have an R² of 1.0
87      return 1.0;
88    }
89
90    public IGrammar Grammar {
91      get { return grammar; }
92    }
93
94    public double Evaluate(string sentence) {
95      return OptimizeConstantsAndEvaluate(sentence);
96    }
97
98    public double SimpleEvaluate(string sentence) {
99      var interpreter = new ExpressionInterpreter();
100      var rowData = new double[d];
101      return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
102        for (int j = 0; j < d; j++) rowData[j] = x[i, j];
103        return interpreter.Interpret(sentence, rowData, erc);
104      }));
105    }
106
107
108    public string CanonicalRepresentation(string phrase) {
109      return phrase;
110      //var terms = phrase.Split('+');
111      //return string.Join("+", terms.Select(term => string.Join("", term.Replace("*", "").OrderBy(ch => ch)))
112      //  .OrderBy(term => term));
113    }
114
115    public IEnumerable<Feature> GetFeatures(string phrase) {
116      throw new NotImplementedException();
117    }
118
119
120
121    public double OptimizeConstantsAndEvaluate(string sentence) {
122      AutoDiff.Term func;
123      int pos = 0;
124
125      var compiler = new ExpressionCompiler();
126      Variable[] variables;
127      Variable[] constants;
128      compiler.Compile(sentence, out func, out variables, out constants);
129
130      // constants are manipulated
131
132      if (!constants.Any()) return SimpleEvaluate(sentence);
133
134      AutoDiff.IParametricCompiledTerm compiledFunc = func.Compile(constants, variables); // variate constants leave variables fixed to data
135
136      double[] c = constants.Select(_ => 1.0).ToArray(); // start with ones
137
138      alglib.lsfitstate state;
139      alglib.lsfitreport rep;
140      int info;
141
142      int n = x.GetLength(0);
143      int m = x.GetLength(1);
144      int k = c.Length;
145
146      alglib.ndimensional_pfunc function_cx_1_func = CreatePFunc(compiledFunc);
147      alglib.ndimensional_pgrad function_cx_1_grad = CreatePGrad(compiledFunc);
148
149      const int maxIterations = 10;
150      try {
151        alglib.lsfitcreatefg(x, y, c, n, m, k, false, out state);
152        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
153        //alglib.lsfitsetgradientcheck(state, 0.001);
154        alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
155        alglib.lsfitresults(state, out info, out c, out rep);
156      } catch (ArithmeticException) {
157        return 0.0;
158      } catch (alglib.alglibexception) {
159        return 0.0;
160      }
161
162      //info == -7  => constant optimization failed due to wrong gradient
163      if (info == -7) throw new ArgumentException();
164      {
165        var rowData = new double[d];
166        return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
167          for (int j = 0; j < d; j++) rowData[j] = x[i, j];
168          return compiledFunc.Evaluate(c, rowData);
169        }));
170      }
171    }
172
173
174    private static alglib.ndimensional_pfunc CreatePFunc(AutoDiff.IParametricCompiledTerm compiledFunc) {
175      return (double[] c, double[] x, ref double func, object o) => {
176        func = compiledFunc.Evaluate(c, x);
177      };
178    }
179
180    private static alglib.ndimensional_pgrad CreatePGrad(AutoDiff.IParametricCompiledTerm compiledFunc) {
181      return (double[] c, double[] x, ref double func, double[] grad, object o) => {
182        var tupel = compiledFunc.Differentiate(c, x);
183        func = tupel.Item2;
184        Array.Copy(tupel.Item1, grad, grad.Length);
185      };
186    }
187
188    public IGrammar TreeBasedGPGrammar { get; private set; }
189    public string ConvertTreeToSentence(ISymbolicExpressionTree tree) {
190      var sb = new StringBuilder();
191
192      TreeToSentence(tree.Root.GetSubtree(0).GetSubtree(0), sb);
193
194      return sb.ToString();
195    }
196
197    private void TreeToSentence(ISymbolicExpressionTreeNode treeNode, StringBuilder sb) {
198      if (treeNode.SubtreeCount == 0) {
199        // terminal
200        sb.Append(treeNode.Symbol.Name);
201      } else {
202        string op = string.Empty;
203        switch (treeNode.Symbol.Name) {
204          case "S": op = "+"; break;
205          case "N": op = "-"; break;
206          case "P": op = "*"; break;
207          case "D": op = "%"; break;
208          default: {
209              Debug.Assert(treeNode.SubtreeCount == 1);
210              break;
211            }
212        }
213        // nonterminal
214        if (op == "+" || op == "-") sb.Append("(");
215        TreeToSentence(treeNode.Subtrees.First(), sb);
216        foreach (var subTree in treeNode.Subtrees.Skip(1)) {
217          sb.Append(op);
218          TreeToSentence(subTree, sb);
219        }
220        if (op == "+" || op == "-") sb.Append(")");
221      }
222    }
223  }
224}
Note: See TracBrowser for help on using the repository browser.