1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022016 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


21  using System;


22  using System.Collections.Generic;


23  using System.Diagnostics.Contracts;


24  using System.Linq;


25 


26  namespace 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]) < 1E12 /* == 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]) < 1E12 /* == 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]) < 1E12 /* == 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  }

