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

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

#2283: worked on transformation of arithmetic expressions to canonical form

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