source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs @ 13645

Last change on this file since 13645 was 13645, checked in by gkronber, 3 years ago

#2581: added an MCTS for symbolic regression models

File size: 18.4 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22using System;
23using System.Collections.Generic;
24using System.Diagnostics.Contracts;
25using System.Linq;
26
27namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
28  // evalutes expressions (on vectors)
29  internal class ExpressionEvaluator {
30    // manages it's own vector buffers
31    private readonly List<double[]> vectorBuffers = new List<double[]>();
32    private readonly List<double[]> scalarBuffers = new List<double[]>(); // scalars are vectors of length 1 (to allow mixing scalars and vectors on the same stack)
33
34
35    private double[] GetVectorBuffer() {
36      var v = vectorBuffers[vectorBuffers.Count - 1];
37      vectorBuffers.RemoveAt(vectorBuffers.Count - 1);
38      return v;
39    }
40    private double[] GetScalarBuffer() {
41      var v = scalarBuffers[scalarBuffers.Count - 1];
42      scalarBuffers.RemoveAt(scalarBuffers.Count - 1);
43      return v;
44    }
45
46    private void ReleaseBuffer(double[] buf) {
47      (buf.Length == 1 ? scalarBuffers : vectorBuffers).Add(buf);
48    }
49
50    public const int MaxStackSize = 100;
51    public const int MaxParams = 50;
52    private readonly int vLen;
53    private readonly double lowerEstimationLimit;
54    private readonly double upperEstimationLimit;
55    private readonly double nanReplacementValue;
56
57    private readonly double[][] stack;
58    private readonly double[][][] gradientStack;
59
60    // preallocate stack and gradient stack
61    public ExpressionEvaluator(int vLen, double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue) {
62      if (vLen <= 1) throw new ArgumentException("number of elements in a variable must be > 1", "vlen");
63      this.vLen = vLen;
64      this.lowerEstimationLimit = lowerEstimationLimit;
65      this.upperEstimationLimit = upperEstimationLimit;
66      this.nanReplacementValue = (upperEstimationLimit - lowerEstimationLimit) / 2.0 + lowerEstimationLimit;
67
68      stack = new double[MaxStackSize][];
69      gradientStack = new double[MaxParams][][];
70
71      for (int k = 0; k < MaxParams; k++) {
72        gradientStack[k] = new double[MaxStackSize][];
73      }
74
75      // preallocate buffers
76      for (int i = 0; i < MaxStackSize; i++) {
77        ReleaseBuffer(new double[vLen]);
78        ReleaseBuffer(new double[1]);
79
80        for (int k = 0; k < MaxParams; k++) {
81          ReleaseBuffer(new double[vLen]);
82          ReleaseBuffer(new double[1]);
83        }
84      }
85    }
86
87    // pred must be allocated by the caller
88    // if adjustOffsetForLogAndExp is set to true we determine c in log(c + f(x)) to make sure that c + f(x) is positive
89    public void Exec(byte[] code, double[][] vars, double[] consts, double[] pred, bool adjustOffsetForLogAndExp = false) {
90      Contract.Assert(pred != null && pred.Length >= vLen);
91      int topOfStack = -1;
92      int pc = 0;
93      int curParamIdx = -1;
94      byte op;
95      short arg;
96      // checked at the end to make sure we do not leak buffers
97      int initialScalarCount = scalarBuffers.Count;
98      int initialVectorCount = vectorBuffers.Count;
99
100      while (true) {
101        ReadNext(code, ref pc, out op, out arg);
102        switch (op) {
103          case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed
104          case (byte)OpCodes.LoadConst0: {
105              ++topOfStack;
106              var z = GetScalarBuffer();
107              z[0] = 0;
108              stack[topOfStack] = z;
109              break;
110            }
111          case (byte)OpCodes.LoadConst1: {
112              ++topOfStack;
113              var z = GetScalarBuffer();
114              z[0] = 1.0;
115              stack[topOfStack] = z;
116              break;
117            }
118          case (byte)OpCodes.LoadParamN: {
119              ++topOfStack;
120              var c = consts[++curParamIdx];
121              var z = GetScalarBuffer();
122              z[0] = c;
123              stack[topOfStack] = z;
124              break;
125            }
126          case (byte)OpCodes.LoadVar: {
127              ++topOfStack;
128              var z = GetVectorBuffer();
129              Array.Copy(vars[arg], z, vars[arg].Length);
130              stack[topOfStack] = z;
131              break;
132            }
133          case (byte)OpCodes.Add: {
134              topOfStack--;
135              var a = stack[topOfStack + 1];
136              var b = stack[topOfStack];
137              stack[topOfStack] = Add(a, b);
138              ReleaseBuffer(a);
139              ReleaseBuffer(b);
140              break;
141            }
142          case (byte)OpCodes.Mul: {
143              topOfStack--;
144              var a = stack[topOfStack + 1];
145              var b = stack[topOfStack];
146              stack[topOfStack] = Mul(a, b);
147              ReleaseBuffer(a);
148              ReleaseBuffer(b);
149              break;
150            }
151          case (byte)OpCodes.Log: {
152              if (adjustOffsetForLogAndExp) {
153                // here we assume that the last used parameter is c in log(f(x) + c)
154                // this must match actions for producing code in the automaton!
155
156                // we can easily adjust c to make sure that f(x) + c is positive because at this point we all values for f(x)
157                var fxc = stack[topOfStack];
158                var minFx = fxc.Min() - consts[curParamIdx]; // stack[topOfStack] is f(x) + c
159
160                var delta = 1.0 - minFx - consts[curParamIdx];
161                // adjust c so that minFx + c = 1 ... log(minFx + c) = 0
162                consts[curParamIdx] += delta;
163
164                // also adjust values on stack
165                for (int i = 0; i < fxc.Length; i++) fxc[i] += delta;
166              }
167              var x = stack[topOfStack];
168              for (int i = 0; i < x.Length; i++)
169                x[i] = Math.Log(x[i]);
170              break;
171            }
172          case (byte)OpCodes.Exp: {
173              if (adjustOffsetForLogAndExp) {
174                // here we assume that the last used parameter is c in exp(f(x) * c)
175                // this must match actions for producing code in the automaton!
176
177                // adjust c to make sure that exp(f(x) * c) is not too large
178                var fxc = stack[topOfStack];
179                var maxFx = fxc.Max() / consts[curParamIdx]; // stack[topOfStack] is f(x) * c
180
181                var f = 1.0 / (maxFx * consts[curParamIdx]);
182                // adjust c so that maxFx*c = 1
183                consts[curParamIdx] *= f;
184
185                // also adjust values on stack
186                for (int i = 0; i < fxc.Length; i++) fxc[i] *= f;
187              }
188
189              var x = stack[topOfStack];
190              for (int i = 0; i < x.Length; i++)
191                x[i] = Math.Exp(x[i]);
192              break;
193            }
194          case (byte)OpCodes.Inv: {
195              var x = stack[topOfStack];
196              for (int i = 0; i < x.Length; i++)
197                x[i] = 1.0 / (x[i]);
198              break;
199            }
200          case (byte)OpCodes.Exit:
201            Contract.Assert(topOfStack == 0);
202            var r = stack[topOfStack];
203            if (r.Length == 1) {
204              var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
205              for (int i = 0; i < vLen; i++)
206                pred[i] = v;
207            } else {
208              for (int i = 0; i < vLen; i++) {
209                var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
210                pred[i] = v;
211              }
212            }
213            ReleaseBuffer(r);
214            Contract.Assert(vectorBuffers.Count == initialVectorCount);
215            Contract.Assert(scalarBuffers.Count == initialScalarCount);
216            return;
217        }
218      }
219    }
220
221
222    // evaluation with forward autodiff
223    // pred and gradients must be allocated by the caller
224    public void ExecGradient(byte[] code, double[][] vars, double[] consts, double[] pred, double[][] gradients) {
225      Contract.Assert(pred != null && pred.Length >= vLen);
226      int topOfStack = -1;
227      int pc = 0;
228      int curParamIdx = -1;
229      byte op;
230      short arg;
231      int nParams = consts.Length;
232      Contract.Assert(gradients != null && gradients.Length >= nParams && gradients.All(g => g.Length >= vLen));
233
234      // checked at the end to make sure we do not leak buffers
235      int initialScalarCount = scalarBuffers.Count;
236      int initialVectorCount = vectorBuffers.Count;
237
238      while (true) {
239        ReadNext(code, ref pc, out op, out arg);
240        switch (op) {
241          case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed
242          case (byte)OpCodes.LoadConst0: {
243              ++topOfStack;
244              var z = GetScalarBuffer();
245              z[0] = 0;
246              stack[topOfStack] = z;
247              for (int k = 0; k < nParams; ++k) {
248                var b = GetScalarBuffer();
249                b[0] = 0.0;
250                gradientStack[k][topOfStack] = b;
251              }
252              break;
253            }
254          case (byte)OpCodes.LoadConst1: {
255              ++topOfStack;
256              var z = GetScalarBuffer();
257              z[0] = 1.0;
258              stack[topOfStack] = z;
259              for (int k = 0; k < nParams; ++k) {
260                var b = GetScalarBuffer();
261                b[0] = 0.0;
262                gradientStack[k][topOfStack] = b;
263              }
264              break;
265            }
266          case (byte)OpCodes.LoadParamN: {
267              var c = consts[++curParamIdx];
268              ++topOfStack;
269              var z = GetScalarBuffer();
270              z[0] = c;
271              stack[topOfStack] = z;
272              for (int k = 0; k < nParams; ++k) {
273                var b = GetScalarBuffer();
274                b[0] = k == curParamIdx ? 1.0 : 0.0;
275                gradientStack[k][topOfStack] = b;
276              }
277              break;
278            }
279          case (byte)OpCodes.LoadVar: {
280              ++topOfStack;
281              var z = GetVectorBuffer();
282              Array.Copy(vars[arg], z, vars[arg].Length);
283              stack[topOfStack] = z;
284              for (int k = 0; k < nParams; ++k) {
285                var b = GetScalarBuffer();
286                b[0] = 0.0;
287                gradientStack[k][topOfStack] = b;
288              }
289            }
290            break;
291          case (byte)OpCodes.Add: {
292              topOfStack--;
293              var a = stack[topOfStack + 1];
294              var b = stack[topOfStack];
295              stack[topOfStack] = Add(a, b);
296              ReleaseBuffer(a);
297              ReleaseBuffer(b);
298
299              // same for gradient
300              for (int k = 0; k < nParams; ++k) {
301                var ag = gradientStack[k][topOfStack + 1];
302                var bg = gradientStack[k][topOfStack];
303                gradientStack[k][topOfStack] = Add(ag, bg);
304                ReleaseBuffer(ag);
305                ReleaseBuffer(bg);
306              }
307              break;
308            }
309          case (byte)OpCodes.Mul: {
310              topOfStack--;
311              var a = stack[topOfStack + 1];
312              var b = stack[topOfStack];
313              stack[topOfStack] = Mul(a, b);
314
315              // same for gradient
316              // f(x) g(x)  f '(x) g(x) + f(x) g'(x)
317              for (int k = 0; k < nParams; ++k) {
318                var ag = gradientStack[k][topOfStack + 1];
319                var bg = gradientStack[k][topOfStack];
320                var t1 = Mul(ag, b);
321                var t2 = Mul(a, bg);
322                gradientStack[k][topOfStack] = Add(t1, t2);
323                ReleaseBuffer(ag);
324                ReleaseBuffer(bg);
325                ReleaseBuffer(t1);
326                ReleaseBuffer(t2);
327              }
328
329              ReleaseBuffer(a);
330              ReleaseBuffer(b);
331
332              break;
333            }
334          case (byte)OpCodes.Log: {
335              var x = stack[topOfStack];
336              // calc gradients first before destroying x
337              // log(f(x))' = f(x)'/f(x)
338              for (int k = 0; k < nParams; k++) {
339                var xg = gradientStack[k][topOfStack];
340                gradientStack[k][topOfStack] = Frac(xg, x);
341                ReleaseBuffer(xg);
342              }
343
344              for (int i = 0; i < x.Length; i++)
345                x[i] = Math.Log(x[i]);
346
347              break;
348            }
349          case (byte)OpCodes.Exp: {
350              var x = stack[topOfStack];
351              for (int i = 0; i < x.Length; i++)
352                x[i] = Math.Exp(x[i]);
353
354              for (int k = 0; k < nParams; k++) {
355                var xg = gradientStack[k][topOfStack];
356                gradientStack[k][topOfStack] = Mul(x, xg);  // e(f(x))' = e(f(x)) * f(x)'
357                ReleaseBuffer(xg);
358              }
359              break;
360            }
361          case (byte)OpCodes.Inv: {
362              var x = stack[topOfStack];
363              for (int i = 0; i < x.Length; i++)
364                x[i] = 1.0 / x[i];
365
366              for (int k = 0; k < nParams; k++) {
367                var xg = gradientStack[k][topOfStack];
368                // x has already been inverted above
369                // (1/f)' = -f' / f²
370                var invF = Mul(xg, x);
371                gradientStack[k][topOfStack] = Mul(invF, x, factor: -1.0);
372                ReleaseBuffer(xg);
373                ReleaseBuffer(invF);
374              }
375              break;
376            }
377          case (byte)OpCodes.Exit:
378            Contract.Assert(topOfStack == 0);
379            var r = stack[topOfStack];
380            if (r.Length == 1) {
381              var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
382              for (int i = 0; i < vLen; i++)
383                pred[i] = v;
384            } else {
385              for (int i = 0; i < vLen; i++) {
386                var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
387                pred[i] = v;
388              }
389            }
390            ReleaseBuffer(r);
391
392            // same for gradients
393            for (int k = 0; k < nParams; k++) {
394              var g = gradientStack[k][topOfStack];
395              if (g.Length == 1) {
396                for (int i = 0; i < vLen; i++)
397                  gradients[k][i] = g[0];
398              } else
399                Array.Copy(g, gradients[k], vLen);
400              ReleaseBuffer(g);
401            }
402
403            Contract.Assert(vectorBuffers.Count == initialVectorCount);
404            Contract.Assert(scalarBuffers.Count == initialScalarCount);
405            return; // break loop
406        }
407      }
408    }
409
410    private double[] Add(double[] a, double[] b) {
411      double[] target = null;
412      if (a.Length > 1) {
413        target = GetVectorBuffer();
414        if (b.Length > 1) {
415          for (int i = 0; i < vLen; i++)
416            target[i] = a[i] + b[i];
417        } else {
418          // b == scalar
419          for (int i = 0; i < vLen; i++)
420            target[i] = a[i] + b[0];
421        }
422      } else {
423        // a == scalar
424        if (b.Length > 1) {
425          target = GetVectorBuffer();
426          for (int i = 0; i < vLen; i++)
427            target[i] = a[0] + b[i];
428        } else {
429          // b == scalar
430          target = GetScalarBuffer();
431          target[0] = a[0] + b[0];
432        }
433      }
434      return target;
435    }
436
437    private double[] Mul(double[] a, double[] b, double factor = 1.0) {
438      double[] target = null;
439      if (a.Length > 1) {
440        if (b.Length > 1) {
441          target = GetVectorBuffer();
442          for (int i = 0; i < vLen; i++)
443            target[i] = factor * a[i] * b[i];
444        } else {
445          // b == scalar
446          if (Math.Abs(b[0]) < 1E-12 /* == 0 */) {
447            target = GetScalarBuffer();
448            target[0] = 0.0;
449          } else {
450            target = GetVectorBuffer();
451            for (int i = 0; i < vLen; i++)
452              target[i] = factor * a[i] * b[0];
453          }
454        }
455      } else {
456        // a == scalar
457        if (b.Length > 1) {
458          if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
459            target = GetScalarBuffer();
460            target[0] = 0.0;
461          } else {
462            target = GetVectorBuffer();
463            for (int i = 0; i < vLen; i++)
464              target[i] = factor * a[0] * b[i];
465          }
466        } else {
467          // b == scalar
468          target = GetScalarBuffer();
469          target[0] = factor * a[0] * b[0];
470        }
471      }
472      return target;
473    }
474
475    private double[] Frac(double[] a, double[] b) {
476      double[] target = null;
477      if (a.Length > 1) {
478        target = GetVectorBuffer();
479        if (b.Length > 1) {
480          for (int i = 0; i < vLen; i++)
481            target[i] = a[i] / b[i];
482        } else {
483          // b == scalar
484          for (int i = 0; i < vLen; i++)
485            target[i] = a[i] / b[0];
486        }
487      } else {
488        // a == scalar
489        if (b.Length > 1) {
490          if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
491            target = GetScalarBuffer();
492            target[0] = 0.0;
493          } else {
494            target = GetVectorBuffer();
495            for (int i = 0; i < vLen; i++)
496              target[i] = a[0] / b[i];
497          }
498        } else {
499          // b == scalar
500          target = GetScalarBuffer();
501          target[0] = a[0] / b[0];
502        }
503      }
504      return target;
505    }
506
507    private void ReadNext(byte[] code, ref int pc, out byte op, out short s) {
508      op = code[pc++];
509      s = 0;
510      if (op == (byte)OpCodes.LoadVar) {
511        s = (short)(((short)code[pc] << 8) | (short)code[pc + 1]);
512        pc += 2;
513      }
514    }
515  }
516}
Note: See TracBrowser for help on using the repository browser.