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

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

#2283: changed handling of inverse expressions in transformation of expressions to canonical form

File size: 9.8 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      var extender = new ExpressionExtender();
129      sentence = extender.CanonicalRepresentation(sentence);
130      if (useConstantOpt)
131        return OptimizeConstantsAndEvaluate(sentence);
132      else {
133
134        Debug.Assert(SimpleEvaluate(sentence) == SimpleEvaluate(extender.CanonicalRepresentation(sentence)));
135        return SimpleEvaluate(sentence);
136      }
137    }
138
139    public double SimpleEvaluate(string sentence) {
140      var interpreter = new ExpressionInterpreter();
141      var rowData = new double[d];
142      return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
143        for (int j = 0; j < d; j++) rowData[j] = x[i, j];
144        return interpreter.Interpret(sentence, rowData, erc);
145      }));
146    }
147
148
149    public string CanonicalRepresentation(string phrase) {
150      var extender = new ExpressionExtender();
151      return extender.CanonicalRepresentation(phrase);
152    }
153
154    public IEnumerable<Feature> GetFeatures(string phrase) {
155      // throw new NotImplementedException();
156      phrase = CanonicalRepresentation(phrase);
157      return phrase.Split('+').Distinct().Select(t => new Feature(t, 1.0));
158      // return new Feature[] { new Feature(phrase, 1.0) };
159    }
160
161
162
163    public double OptimizeConstantsAndEvaluate(string sentence) {
164      AutoDiff.Term func;
165      int n = x.GetLength(0);
166      int m = x.GetLength(1);
167
168      var compiler = new ExpressionCompiler();
169      Variable[] variables = Enumerable.Range(0, m).Select(_ => new Variable()).ToArray();
170      Variable[] constants;
171      compiler.Compile(sentence, out func, variables, out constants);
172
173      // constants are manipulated
174
175      if (!constants.Any()) return SimpleEvaluate(sentence);
176
177      AutoDiff.IParametricCompiledTerm compiledFunc = func.Compile(constants, variables); // variate constants leave variables fixed to data
178
179      double[] c = constants.Select(_ => 1.0).ToArray(); // start with ones
180
181      alglib.lsfitstate state;
182      alglib.lsfitreport rep;
183      int info;
184
185
186      int k = c.Length;
187
188      alglib.ndimensional_pfunc function_cx_1_func = CreatePFunc(compiledFunc);
189      alglib.ndimensional_pgrad function_cx_1_grad = CreatePGrad(compiledFunc);
190
191      const int maxIterations = 10;
192      try {
193        alglib.lsfitcreatefg(x, y, c, n, m, k, false, out state);
194        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
195        //alglib.lsfitsetgradientcheck(state, 0.001);
196        alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
197        alglib.lsfitresults(state, out info, out c, out rep);
198      } catch (ArithmeticException) {
199        return 0.0;
200      } catch (alglib.alglibexception) {
201        return 0.0;
202      }
203
204      //info == -7  => constant optimization failed due to wrong gradient
205      if (info == -7) throw new ArgumentException();
206      {
207        var rowData = new double[d];
208        return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
209          for (int j = 0; j < d; j++) rowData[j] = x[i, j];
210          return compiledFunc.Evaluate(c, rowData);
211        }));
212      }
213    }
214
215
216    private static alglib.ndimensional_pfunc CreatePFunc(AutoDiff.IParametricCompiledTerm compiledFunc) {
217      return (double[] c, double[] x, ref double func, object o) => {
218        func = compiledFunc.Evaluate(c, x);
219      };
220    }
221
222    private static alglib.ndimensional_pgrad CreatePGrad(AutoDiff.IParametricCompiledTerm compiledFunc) {
223      return (double[] c, double[] x, ref double func, double[] grad, object o) => {
224        var tupel = compiledFunc.Differentiate(c, x);
225        func = tupel.Item2;
226        Array.Copy(tupel.Item1, grad, grad.Length);
227      };
228    }
229
230    public IGrammar TreeBasedGPGrammar { get; private set; }
231    public string ConvertTreeToSentence(ISymbolicExpressionTree tree) {
232      var sb = new StringBuilder();
233
234      TreeToSentence(tree.Root.GetSubtree(0).GetSubtree(0), sb);
235
236      return sb.ToString();
237    }
238
239
240    private void TreeToSentence(ISymbolicExpressionTreeNode treeNode, StringBuilder sb) {
241      if (treeNode.SubtreeCount == 0) {
242        // terminal
243        sb.Append(treeNode.Symbol.Name);
244      } else {
245        string op = string.Empty;
246        switch (treeNode.Symbol.Name) {
247          case "S": op = "+"; break;
248          case "N": op = "-"; break;
249          case "P": op = "*"; break;
250          case "D": op = "%"; break;
251          default: {
252              Debug.Assert(treeNode.SubtreeCount == 1);
253              break;
254            }
255        }
256        // nonterminal
257        if (treeNode.SubtreeCount > 1) sb.Append("(");
258        TreeToSentence(treeNode.Subtrees.First(), sb);
259        foreach (var subTree in treeNode.Subtrees.Skip(1)) {
260          sb.Append(op);
261          TreeToSentence(subTree, sb);
262        }
263        if (treeNode.SubtreeCount > 1) sb.Append(")");
264      }
265    }
266  }
267}
Note: See TracBrowser for help on using the repository browser.