[13645]  1  #region License Information


 2  /* HeuristicLab


 3  * Copyright (C) 20022015 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  // evalutes expressions (on vectors)


 28  internal class ExpressionEvaluator {


 29  // manages it's own vector buffers


[13651]  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;


[13645]  34 


 35 


 36  private double[] GetVectorBuffer() {


[13651]  37  return vectorBuffers[lastVecBufIdx];


[13645]  38  }


 39  private double[] GetScalarBuffer() {


[13651]  40  return scalarBuffers[lastScalarBufIdx];


[13645]  41  }


 42 


 43  private void ReleaseBuffer(double[] buf) {


[13651]  44  if (buf.Length == 1) {


 45  scalarBuffers[lastScalarBufIdx++] = buf;


 46  } else {


 47  vectorBuffers[lastVecBufIdx++] = buf;


 48  }


[13645]  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


[13651]  77  vectorBuffers = new double[MaxStackSize * (1 + MaxParams)][];


 78  scalarBuffers = new double[MaxStackSize * (1 + MaxParams)][];


[13645]  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


[14015]  91  // if adjustParameters is set to true we determine parameters in function to make sure that the function output is valid


 92  // e.g. in log(c + f(x)) to make sure that c + f(x) is positive


 93  public void Exec(byte[] code, double[][] vars, double[] consts, double[] pred, bool adjustParameters = false) {


[13645]  94  Contract.Assert(pred != null && pred.Length >= vLen);


 95  int topOfStack = 1;


 96  int pc = 0;


 97  int curParamIdx = 1;


 98  byte op;


 99  short arg;


 100  // checked at the end to make sure we do not leak buffers


[13651]  101  int initialScalarCount = lastScalarBufIdx;


 102  int initialVectorCount = lastVecBufIdx;


[13645]  103 


 104  while (true) {


 105  ReadNext(code, ref pc, out op, out arg);


 106  switch (op) {


 107  case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed


 108  case (byte)OpCodes.LoadConst0: {


 109  ++topOfStack;


 110  var z = GetScalarBuffer();


 111  z[0] = 0;


 112  stack[topOfStack] = z;


 113  break;


 114  }


 115  case (byte)OpCodes.LoadConst1: {


 116  ++topOfStack;


 117  var z = GetScalarBuffer();


 118  z[0] = 1.0;


 119  stack[topOfStack] = z;


 120  break;


 121  }


 122  case (byte)OpCodes.LoadParamN: {


 123  ++topOfStack;


 124  var c = consts[++curParamIdx];


 125  var z = GetScalarBuffer();


 126  z[0] = c;


 127  stack[topOfStack] = z;


 128  break;


 129  }


 130  case (byte)OpCodes.LoadVar: {


 131  ++topOfStack;


 132  var z = GetVectorBuffer();


 133  Array.Copy(vars[arg], z, vars[arg].Length);


 134  stack[topOfStack] = z;


 135  break;


 136  }


 137  case (byte)OpCodes.Add: {


 138  topOfStack;


 139  var a = stack[topOfStack + 1];


 140  var b = stack[topOfStack];


 141  stack[topOfStack] = Add(a, b);


 142  ReleaseBuffer(a);


 143  ReleaseBuffer(b);


 144  break;


 145  }


 146  case (byte)OpCodes.Mul: {


 147  topOfStack;


 148  var a = stack[topOfStack + 1];


 149  var b = stack[topOfStack];


 150  stack[topOfStack] = Mul(a, b);


 151  ReleaseBuffer(a);


 152  ReleaseBuffer(b);


 153  break;


 154  }


 155  case (byte)OpCodes.Log: {


[14015]  156  if (adjustParameters) {


[13645]  157  // here we assume that the last used parameter is c in log(f(x) + c)


 158  // this must match actions for producing code in the automaton!


 159 


 160  // we can easily adjust c to make sure that f(x) + c is positive because at this point we all values for f(x)


 161  var fxc = stack[topOfStack];


 162  var minFx = fxc.Min()  consts[curParamIdx]; // stack[topOfStack] is f(x) + c


 163 


 164  var delta = 1.0  minFx  consts[curParamIdx];


 165  // adjust c so that minFx + c = 1 ... log(minFx + c) = 0


 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: {


[14015]  177  if (adjustParameters) {


[13645]  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]);


[13652]  186  // adjust c so that maxFx*c = 1 TODO: this is not ideal as it enforces positive arguments to exp()


[13645]  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);


[13651]  218  Contract.Assert(lastVecBufIdx == initialVectorCount);


 219  Contract.Assert(lastScalarBufIdx == initialScalarCount);


[13645]  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


[13651]  239  int initialScalarCount = lastScalarBufIdx;


 240  int initialVectorCount = lastVecBufIdx;


[13645]  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 


[13651]  407  Contract.Assert(lastVecBufIdx == initialVectorCount);


 408  Contract.Assert(lastScalarBufIdx == initialScalarCount);


[13645]  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)(((short)code[pc] << 8)  (short)code[pc + 1]);


 516  pc += 2;


 517  }


 518  }


 519  }


 520  }

