Free cookie consent management tool by TermsFeed Policy Generator

source: branches/MCTS-SymbReg-2796/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs @ 15606

Last change on this file since 15606 was 15606, checked in by gkronber, 6 years ago

#2796: comments and typos

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