source: stable/HeuristicLab.Algorithms.DataAnalysis/3.4/MctsSymbolicRegression/ExpressionEvaluator.cs @ 15060

Last change on this file since 15060 was 15060, checked in by gkronber, 2 years ago

#2581: merged r13645,r13648,r13650,r13651,r13652,r13654,r13657,r13658,r13659,r13661,r13662,r13669,r13708,r14142 from trunk to stable (to be deleted in the next commit)

File size: 18.5 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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  // evalutes 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                var delta = 1.0 - minFx - consts[curParamIdx];
164                // adjust c so that minFx + c = 1 ... log(minFx + c) = 0
165                consts[curParamIdx] += delta;
166
167                // also adjust values on stack
168                for (int i = 0; i < fxc.Length; i++) fxc[i] += delta;
169              }
170              var x = stack[topOfStack];
171              for (int i = 0; i < x.Length; i++)
172                x[i] = Math.Log(x[i]);
173              break;
174            }
175          case (byte)OpCodes.Exp: {
176              if (adjustOffsetForLogAndExp) {
177                // here we assume that the last used parameter is c in exp(f(x) * c)
178                // this must match actions for producing code in the automaton!
179
180                // adjust c to make sure that exp(f(x) * c) is not too large
181                var fxc = stack[topOfStack];
182                var maxFx = fxc.Max() / consts[curParamIdx]; // stack[topOfStack] is f(x) * c
183
184                var f = 1.0 / (maxFx * consts[curParamIdx]);
185                // adjust c so that maxFx*c = 1 TODO: this is not ideal as it enforces positive arguments to exp()
186                consts[curParamIdx] *= f;
187
188                // also adjust values on stack
189                for (int i = 0; i < fxc.Length; i++) fxc[i] *= f;
190              }
191
192              var x = stack[topOfStack];
193              for (int i = 0; i < x.Length; i++)
194                x[i] = Math.Exp(x[i]);
195              break;
196            }
197          case (byte)OpCodes.Inv: {
198              var x = stack[topOfStack];
199              for (int i = 0; i < x.Length; i++)
200                x[i] = 1.0 / (x[i]);
201              break;
202            }
203          case (byte)OpCodes.Exit:
204            Contract.Assert(topOfStack == 0);
205            var r = stack[topOfStack];
206            if (r.Length == 1) {
207              var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
208              for (int i = 0; i < vLen; i++)
209                pred[i] = v;
210            } else {
211              for (int i = 0; i < vLen; i++) {
212                var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
213                pred[i] = v;
214              }
215            }
216            ReleaseBuffer(r);
217            Contract.Assert(lastVecBufIdx == initialVectorCount);
218            Contract.Assert(lastScalarBufIdx == initialScalarCount);
219            return;
220        }
221      }
222    }
223
224
225    // evaluation with forward autodiff
226    // pred and gradients must be allocated by the caller
227    public void ExecGradient(byte[] code, double[][] vars, double[] consts, double[] pred, double[][] gradients) {
228      Contract.Assert(pred != null && pred.Length >= vLen);
229      int topOfStack = -1;
230      int pc = 0;
231      int curParamIdx = -1;
232      byte op;
233      short arg;
234      int nParams = consts.Length;
235      Contract.Assert(gradients != null && gradients.Length >= nParams && gradients.All(g => g.Length >= vLen));
236
237      // checked at the end to make sure we do not leak buffers
238      int initialScalarCount = lastScalarBufIdx;
239      int initialVectorCount = lastVecBufIdx;
240
241      while (true) {
242        ReadNext(code, ref pc, out op, out arg);
243        switch (op) {
244          case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed
245          case (byte)OpCodes.LoadConst0: {
246              ++topOfStack;
247              var z = GetScalarBuffer();
248              z[0] = 0;
249              stack[topOfStack] = z;
250              for (int k = 0; k < nParams; ++k) {
251                var b = GetScalarBuffer();
252                b[0] = 0.0;
253                gradientStack[k][topOfStack] = b;
254              }
255              break;
256            }
257          case (byte)OpCodes.LoadConst1: {
258              ++topOfStack;
259              var z = GetScalarBuffer();
260              z[0] = 1.0;
261              stack[topOfStack] = z;
262              for (int k = 0; k < nParams; ++k) {
263                var b = GetScalarBuffer();
264                b[0] = 0.0;
265                gradientStack[k][topOfStack] = b;
266              }
267              break;
268            }
269          case (byte)OpCodes.LoadParamN: {
270              var c = consts[++curParamIdx];
271              ++topOfStack;
272              var z = GetScalarBuffer();
273              z[0] = c;
274              stack[topOfStack] = z;
275              for (int k = 0; k < nParams; ++k) {
276                var b = GetScalarBuffer();
277                b[0] = k == curParamIdx ? 1.0 : 0.0;
278                gradientStack[k][topOfStack] = b;
279              }
280              break;
281            }
282          case (byte)OpCodes.LoadVar: {
283              ++topOfStack;
284              var z = GetVectorBuffer();
285              Array.Copy(vars[arg], z, vars[arg].Length);
286              stack[topOfStack] = z;
287              for (int k = 0; k < nParams; ++k) {
288                var b = GetScalarBuffer();
289                b[0] = 0.0;
290                gradientStack[k][topOfStack] = b;
291              }
292            }
293            break;
294          case (byte)OpCodes.Add: {
295              topOfStack--;
296              var a = stack[topOfStack + 1];
297              var b = stack[topOfStack];
298              stack[topOfStack] = Add(a, b);
299              ReleaseBuffer(a);
300              ReleaseBuffer(b);
301
302              // same for gradient
303              for (int k = 0; k < nParams; ++k) {
304                var ag = gradientStack[k][topOfStack + 1];
305                var bg = gradientStack[k][topOfStack];
306                gradientStack[k][topOfStack] = Add(ag, bg);
307                ReleaseBuffer(ag);
308                ReleaseBuffer(bg);
309              }
310              break;
311            }
312          case (byte)OpCodes.Mul: {
313              topOfStack--;
314              var a = stack[topOfStack + 1];
315              var b = stack[topOfStack];
316              stack[topOfStack] = Mul(a, b);
317
318              // same for gradient
319              // f(x) g(x)  f '(x) g(x) + f(x) g'(x)
320              for (int k = 0; k < nParams; ++k) {
321                var ag = gradientStack[k][topOfStack + 1];
322                var bg = gradientStack[k][topOfStack];
323                var t1 = Mul(ag, b);
324                var t2 = Mul(a, bg);
325                gradientStack[k][topOfStack] = Add(t1, t2);
326                ReleaseBuffer(ag);
327                ReleaseBuffer(bg);
328                ReleaseBuffer(t1);
329                ReleaseBuffer(t2);
330              }
331
332              ReleaseBuffer(a);
333              ReleaseBuffer(b);
334
335              break;
336            }
337          case (byte)OpCodes.Log: {
338              var x = stack[topOfStack];
339              // calc gradients first before destroying x
340              // log(f(x))' = f(x)'/f(x)
341              for (int k = 0; k < nParams; k++) {
342                var xg = gradientStack[k][topOfStack];
343                gradientStack[k][topOfStack] = Frac(xg, x);
344                ReleaseBuffer(xg);
345              }
346
347              for (int i = 0; i < x.Length; i++)
348                x[i] = Math.Log(x[i]);
349
350              break;
351            }
352          case (byte)OpCodes.Exp: {
353              var x = stack[topOfStack];
354              for (int i = 0; i < x.Length; i++)
355                x[i] = Math.Exp(x[i]);
356
357              for (int k = 0; k < nParams; k++) {
358                var xg = gradientStack[k][topOfStack];
359                gradientStack[k][topOfStack] = Mul(x, xg);  // e(f(x))' = e(f(x)) * f(x)'
360                ReleaseBuffer(xg);
361              }
362              break;
363            }
364          case (byte)OpCodes.Inv: {
365              var x = stack[topOfStack];
366              for (int i = 0; i < x.Length; i++)
367                x[i] = 1.0 / x[i];
368
369              for (int k = 0; k < nParams; k++) {
370                var xg = gradientStack[k][topOfStack];
371                // x has already been inverted above
372                // (1/f)' = -f' / f²
373                var invF = Mul(xg, x);
374                gradientStack[k][topOfStack] = Mul(invF, x, factor: -1.0);
375                ReleaseBuffer(xg);
376                ReleaseBuffer(invF);
377              }
378              break;
379            }
380          case (byte)OpCodes.Exit:
381            Contract.Assert(topOfStack == 0);
382            var r = stack[topOfStack];
383            if (r.Length == 1) {
384              var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
385              for (int i = 0; i < vLen; i++)
386                pred[i] = v;
387            } else {
388              for (int i = 0; i < vLen; i++) {
389                var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
390                pred[i] = v;
391              }
392            }
393            ReleaseBuffer(r);
394
395            // same for gradients
396            for (int k = 0; k < nParams; k++) {
397              var g = gradientStack[k][topOfStack];
398              if (g.Length == 1) {
399                for (int i = 0; i < vLen; i++)
400                  gradients[k][i] = g[0];
401              } else
402                Array.Copy(g, gradients[k], vLen);
403              ReleaseBuffer(g);
404            }
405
406            Contract.Assert(lastVecBufIdx == initialVectorCount);
407            Contract.Assert(lastScalarBufIdx == initialScalarCount);
408            return; // break loop
409        }
410      }
411    }
412
413    private double[] Add(double[] a, double[] b) {
414      double[] target = null;
415      if (a.Length > 1) {
416        target = GetVectorBuffer();
417        if (b.Length > 1) {
418          for (int i = 0; i < vLen; i++)
419            target[i] = a[i] + b[i];
420        } else {
421          // b == scalar
422          for (int i = 0; i < vLen; i++)
423            target[i] = a[i] + b[0];
424        }
425      } else {
426        // a == scalar
427        if (b.Length > 1) {
428          target = GetVectorBuffer();
429          for (int i = 0; i < vLen; i++)
430            target[i] = a[0] + b[i];
431        } else {
432          // b == scalar
433          target = GetScalarBuffer();
434          target[0] = a[0] + b[0];
435        }
436      }
437      return target;
438    }
439
440    private double[] Mul(double[] a, double[] b, double factor = 1.0) {
441      double[] target = null;
442      if (a.Length > 1) {
443        if (b.Length > 1) {
444          target = GetVectorBuffer();
445          for (int i = 0; i < vLen; i++)
446            target[i] = factor * a[i] * b[i];
447        } else {
448          // b == scalar
449          if (Math.Abs(b[0]) < 1E-12 /* == 0 */) {
450            target = GetScalarBuffer();
451            target[0] = 0.0;
452          } else {
453            target = GetVectorBuffer();
454            for (int i = 0; i < vLen; i++)
455              target[i] = factor * a[i] * b[0];
456          }
457        }
458      } else {
459        // a == scalar
460        if (b.Length > 1) {
461          if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
462            target = GetScalarBuffer();
463            target[0] = 0.0;
464          } else {
465            target = GetVectorBuffer();
466            for (int i = 0; i < vLen; i++)
467              target[i] = factor * a[0] * b[i];
468          }
469        } else {
470          // b == scalar
471          target = GetScalarBuffer();
472          target[0] = factor * a[0] * b[0];
473        }
474      }
475      return target;
476    }
477
478    private double[] Frac(double[] a, double[] b) {
479      double[] target = null;
480      if (a.Length > 1) {
481        target = GetVectorBuffer();
482        if (b.Length > 1) {
483          for (int i = 0; i < vLen; i++)
484            target[i] = a[i] / b[i];
485        } else {
486          // b == scalar
487          for (int i = 0; i < vLen; i++)
488            target[i] = a[i] / b[0];
489        }
490      } else {
491        // a == scalar
492        if (b.Length > 1) {
493          if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
494            target = GetScalarBuffer();
495            target[0] = 0.0;
496          } else {
497            target = GetVectorBuffer();
498            for (int i = 0; i < vLen; i++)
499              target[i] = a[0] / b[i];
500          }
501        } else {
502          // b == scalar
503          target = GetScalarBuffer();
504          target[0] = a[0] / b[0];
505        }
506      }
507      return target;
508    }
509
510    private void ReadNext(byte[] code, ref int pc, out byte op, out short s) {
511      op = code[pc++];
512      s = 0;
513      if (op == (byte)OpCodes.LoadVar) {
514        s = (short)((code[pc] << 8) | code[pc + 1]);
515        pc += 2;
516      }
517    }
518  }
519}
Note: See TracBrowser for help on using the repository browser.