1  #region License Information


2  /* HeuristicLab


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


22  using System;


23  using System.Collections.Generic;


24  using System.Diagnostics.Contracts;


25  using System.Linq;


26 


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

