#region License Information
/* HeuristicLab
* Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
* and the BEACON Center for the Study of Evolution in Action.
*
* This file is part of HeuristicLab.
*
* HeuristicLab is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* HeuristicLab is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with HeuristicLab. If not, see .
*/
#endregion
using System;
using System.Collections.Generic;
using System.Diagnostics.Contracts;
using System.Linq;
namespace HeuristicLab.Algorithms.DataAnalysis.MctsSymbolicRegression {
// evalutes expressions (on vectors)
internal class ExpressionEvaluator {
// manages it's own vector buffers
private readonly List vectorBuffers = new List();
private readonly List scalarBuffers = new List(); // scalars are vectors of length 1 (to allow mixing scalars and vectors on the same stack)
private double[] GetVectorBuffer() {
var v = vectorBuffers[vectorBuffers.Count - 1];
vectorBuffers.RemoveAt(vectorBuffers.Count - 1);
return v;
}
private double[] GetScalarBuffer() {
var v = scalarBuffers[scalarBuffers.Count - 1];
scalarBuffers.RemoveAt(scalarBuffers.Count - 1);
return v;
}
private void ReleaseBuffer(double[] buf) {
(buf.Length == 1 ? scalarBuffers : vectorBuffers).Add(buf);
}
public const int MaxStackSize = 100;
public const int MaxParams = 50;
private readonly int vLen;
private readonly double lowerEstimationLimit;
private readonly double upperEstimationLimit;
private readonly double nanReplacementValue;
private readonly double[][] stack;
private readonly double[][][] gradientStack;
// preallocate stack and gradient stack
public ExpressionEvaluator(int vLen, double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue) {
if (vLen <= 1) throw new ArgumentException("number of elements in a variable must be > 1", "vlen");
this.vLen = vLen;
this.lowerEstimationLimit = lowerEstimationLimit;
this.upperEstimationLimit = upperEstimationLimit;
this.nanReplacementValue = (upperEstimationLimit - lowerEstimationLimit) / 2.0 + lowerEstimationLimit;
stack = new double[MaxStackSize][];
gradientStack = new double[MaxParams][][];
for (int k = 0; k < MaxParams; k++) {
gradientStack[k] = new double[MaxStackSize][];
}
// preallocate buffers
for (int i = 0; i < MaxStackSize; i++) {
ReleaseBuffer(new double[vLen]);
ReleaseBuffer(new double[1]);
for (int k = 0; k < MaxParams; k++) {
ReleaseBuffer(new double[vLen]);
ReleaseBuffer(new double[1]);
}
}
}
// pred must be allocated by the caller
// if adjustOffsetForLogAndExp is set to true we determine c in log(c + f(x)) to make sure that c + f(x) is positive
public void Exec(byte[] code, double[][] vars, double[] consts, double[] pred, bool adjustOffsetForLogAndExp = false) {
Contract.Assert(pred != null && pred.Length >= vLen);
int topOfStack = -1;
int pc = 0;
int curParamIdx = -1;
byte op;
short arg;
// checked at the end to make sure we do not leak buffers
int initialScalarCount = scalarBuffers.Count;
int initialVectorCount = vectorBuffers.Count;
while (true) {
ReadNext(code, ref pc, out op, out arg);
switch (op) {
case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed
case (byte)OpCodes.LoadConst0: {
++topOfStack;
var z = GetScalarBuffer();
z[0] = 0;
stack[topOfStack] = z;
break;
}
case (byte)OpCodes.LoadConst1: {
++topOfStack;
var z = GetScalarBuffer();
z[0] = 1.0;
stack[topOfStack] = z;
break;
}
case (byte)OpCodes.LoadParamN: {
++topOfStack;
var c = consts[++curParamIdx];
var z = GetScalarBuffer();
z[0] = c;
stack[topOfStack] = z;
break;
}
case (byte)OpCodes.LoadVar: {
++topOfStack;
var z = GetVectorBuffer();
Array.Copy(vars[arg], z, vars[arg].Length);
stack[topOfStack] = z;
break;
}
case (byte)OpCodes.Add: {
topOfStack--;
var a = stack[topOfStack + 1];
var b = stack[topOfStack];
stack[topOfStack] = Add(a, b);
ReleaseBuffer(a);
ReleaseBuffer(b);
break;
}
case (byte)OpCodes.Mul: {
topOfStack--;
var a = stack[topOfStack + 1];
var b = stack[topOfStack];
stack[topOfStack] = Mul(a, b);
ReleaseBuffer(a);
ReleaseBuffer(b);
break;
}
case (byte)OpCodes.Log: {
if (adjustOffsetForLogAndExp) {
// here we assume that the last used parameter is c in log(f(x) + c)
// this must match actions for producing code in the automaton!
// we can easily adjust c to make sure that f(x) + c is positive because at this point we all values for f(x)
var fxc = stack[topOfStack];
var minFx = fxc.Min() - consts[curParamIdx]; // stack[topOfStack] is f(x) + c
var delta = 1.0 - minFx - consts[curParamIdx];
// adjust c so that minFx + c = 1 ... log(minFx + c) = 0
consts[curParamIdx] += delta;
// also adjust values on stack
for (int i = 0; i < fxc.Length; i++) fxc[i] += delta;
}
var x = stack[topOfStack];
for (int i = 0; i < x.Length; i++)
x[i] = Math.Log(x[i]);
break;
}
case (byte)OpCodes.Exp: {
if (adjustOffsetForLogAndExp) {
// here we assume that the last used parameter is c in exp(f(x) * c)
// this must match actions for producing code in the automaton!
// adjust c to make sure that exp(f(x) * c) is not too large
var fxc = stack[topOfStack];
var maxFx = fxc.Max() / consts[curParamIdx]; // stack[topOfStack] is f(x) * c
var f = 1.0 / (maxFx * consts[curParamIdx]);
// adjust c so that maxFx*c = 1
consts[curParamIdx] *= f;
// also adjust values on stack
for (int i = 0; i < fxc.Length; i++) fxc[i] *= f;
}
var x = stack[topOfStack];
for (int i = 0; i < x.Length; i++)
x[i] = Math.Exp(x[i]);
break;
}
case (byte)OpCodes.Inv: {
var x = stack[topOfStack];
for (int i = 0; i < x.Length; i++)
x[i] = 1.0 / (x[i]);
break;
}
case (byte)OpCodes.Exit:
Contract.Assert(topOfStack == 0);
var r = stack[topOfStack];
if (r.Length == 1) {
var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
for (int i = 0; i < vLen; i++)
pred[i] = v;
} else {
for (int i = 0; i < vLen; i++) {
var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
pred[i] = v;
}
}
ReleaseBuffer(r);
Contract.Assert(vectorBuffers.Count == initialVectorCount);
Contract.Assert(scalarBuffers.Count == initialScalarCount);
return;
}
}
}
// evaluation with forward autodiff
// pred and gradients must be allocated by the caller
public void ExecGradient(byte[] code, double[][] vars, double[] consts, double[] pred, double[][] gradients) {
Contract.Assert(pred != null && pred.Length >= vLen);
int topOfStack = -1;
int pc = 0;
int curParamIdx = -1;
byte op;
short arg;
int nParams = consts.Length;
Contract.Assert(gradients != null && gradients.Length >= nParams && gradients.All(g => g.Length >= vLen));
// checked at the end to make sure we do not leak buffers
int initialScalarCount = scalarBuffers.Count;
int initialVectorCount = vectorBuffers.Count;
while (true) {
ReadNext(code, ref pc, out op, out arg);
switch (op) {
case (byte)OpCodes.Nop: throw new InvalidProgramException(); // not allowed
case (byte)OpCodes.LoadConst0: {
++topOfStack;
var z = GetScalarBuffer();
z[0] = 0;
stack[topOfStack] = z;
for (int k = 0; k < nParams; ++k) {
var b = GetScalarBuffer();
b[0] = 0.0;
gradientStack[k][topOfStack] = b;
}
break;
}
case (byte)OpCodes.LoadConst1: {
++topOfStack;
var z = GetScalarBuffer();
z[0] = 1.0;
stack[topOfStack] = z;
for (int k = 0; k < nParams; ++k) {
var b = GetScalarBuffer();
b[0] = 0.0;
gradientStack[k][topOfStack] = b;
}
break;
}
case (byte)OpCodes.LoadParamN: {
var c = consts[++curParamIdx];
++topOfStack;
var z = GetScalarBuffer();
z[0] = c;
stack[topOfStack] = z;
for (int k = 0; k < nParams; ++k) {
var b = GetScalarBuffer();
b[0] = k == curParamIdx ? 1.0 : 0.0;
gradientStack[k][topOfStack] = b;
}
break;
}
case (byte)OpCodes.LoadVar: {
++topOfStack;
var z = GetVectorBuffer();
Array.Copy(vars[arg], z, vars[arg].Length);
stack[topOfStack] = z;
for (int k = 0; k < nParams; ++k) {
var b = GetScalarBuffer();
b[0] = 0.0;
gradientStack[k][topOfStack] = b;
}
}
break;
case (byte)OpCodes.Add: {
topOfStack--;
var a = stack[topOfStack + 1];
var b = stack[topOfStack];
stack[topOfStack] = Add(a, b);
ReleaseBuffer(a);
ReleaseBuffer(b);
// same for gradient
for (int k = 0; k < nParams; ++k) {
var ag = gradientStack[k][topOfStack + 1];
var bg = gradientStack[k][topOfStack];
gradientStack[k][topOfStack] = Add(ag, bg);
ReleaseBuffer(ag);
ReleaseBuffer(bg);
}
break;
}
case (byte)OpCodes.Mul: {
topOfStack--;
var a = stack[topOfStack + 1];
var b = stack[topOfStack];
stack[topOfStack] = Mul(a, b);
// same for gradient
// f(x) g(x) f '(x) g(x) + f(x) g'(x)
for (int k = 0; k < nParams; ++k) {
var ag = gradientStack[k][topOfStack + 1];
var bg = gradientStack[k][topOfStack];
var t1 = Mul(ag, b);
var t2 = Mul(a, bg);
gradientStack[k][topOfStack] = Add(t1, t2);
ReleaseBuffer(ag);
ReleaseBuffer(bg);
ReleaseBuffer(t1);
ReleaseBuffer(t2);
}
ReleaseBuffer(a);
ReleaseBuffer(b);
break;
}
case (byte)OpCodes.Log: {
var x = stack[topOfStack];
// calc gradients first before destroying x
// log(f(x))' = f(x)'/f(x)
for (int k = 0; k < nParams; k++) {
var xg = gradientStack[k][topOfStack];
gradientStack[k][topOfStack] = Frac(xg, x);
ReleaseBuffer(xg);
}
for (int i = 0; i < x.Length; i++)
x[i] = Math.Log(x[i]);
break;
}
case (byte)OpCodes.Exp: {
var x = stack[topOfStack];
for (int i = 0; i < x.Length; i++)
x[i] = Math.Exp(x[i]);
for (int k = 0; k < nParams; k++) {
var xg = gradientStack[k][topOfStack];
gradientStack[k][topOfStack] = Mul(x, xg); // e(f(x))' = e(f(x)) * f(x)'
ReleaseBuffer(xg);
}
break;
}
case (byte)OpCodes.Inv: {
var x = stack[topOfStack];
for (int i = 0; i < x.Length; i++)
x[i] = 1.0 / x[i];
for (int k = 0; k < nParams; k++) {
var xg = gradientStack[k][topOfStack];
// x has already been inverted above
// (1/f)' = -f' / f²
var invF = Mul(xg, x);
gradientStack[k][topOfStack] = Mul(invF, x, factor: -1.0);
ReleaseBuffer(xg);
ReleaseBuffer(invF);
}
break;
}
case (byte)OpCodes.Exit:
Contract.Assert(topOfStack == 0);
var r = stack[topOfStack];
if (r.Length == 1) {
var v = double.IsNaN(r[0]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[0]));
for (int i = 0; i < vLen; i++)
pred[i] = v;
} else {
for (int i = 0; i < vLen; i++) {
var v = double.IsNaN(r[i]) ? nanReplacementValue : Math.Min(upperEstimationLimit, Math.Max(lowerEstimationLimit, r[i]));
pred[i] = v;
}
}
ReleaseBuffer(r);
// same for gradients
for (int k = 0; k < nParams; k++) {
var g = gradientStack[k][topOfStack];
if (g.Length == 1) {
for (int i = 0; i < vLen; i++)
gradients[k][i] = g[0];
} else
Array.Copy(g, gradients[k], vLen);
ReleaseBuffer(g);
}
Contract.Assert(vectorBuffers.Count == initialVectorCount);
Contract.Assert(scalarBuffers.Count == initialScalarCount);
return; // break loop
}
}
}
private double[] Add(double[] a, double[] b) {
double[] target = null;
if (a.Length > 1) {
target = GetVectorBuffer();
if (b.Length > 1) {
for (int i = 0; i < vLen; i++)
target[i] = a[i] + b[i];
} else {
// b == scalar
for (int i = 0; i < vLen; i++)
target[i] = a[i] + b[0];
}
} else {
// a == scalar
if (b.Length > 1) {
target = GetVectorBuffer();
for (int i = 0; i < vLen; i++)
target[i] = a[0] + b[i];
} else {
// b == scalar
target = GetScalarBuffer();
target[0] = a[0] + b[0];
}
}
return target;
}
private double[] Mul(double[] a, double[] b, double factor = 1.0) {
double[] target = null;
if (a.Length > 1) {
if (b.Length > 1) {
target = GetVectorBuffer();
for (int i = 0; i < vLen; i++)
target[i] = factor * a[i] * b[i];
} else {
// b == scalar
if (Math.Abs(b[0]) < 1E-12 /* == 0 */) {
target = GetScalarBuffer();
target[0] = 0.0;
} else {
target = GetVectorBuffer();
for (int i = 0; i < vLen; i++)
target[i] = factor * a[i] * b[0];
}
}
} else {
// a == scalar
if (b.Length > 1) {
if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
target = GetScalarBuffer();
target[0] = 0.0;
} else {
target = GetVectorBuffer();
for (int i = 0; i < vLen; i++)
target[i] = factor * a[0] * b[i];
}
} else {
// b == scalar
target = GetScalarBuffer();
target[0] = factor * a[0] * b[0];
}
}
return target;
}
private double[] Frac(double[] a, double[] b) {
double[] target = null;
if (a.Length > 1) {
target = GetVectorBuffer();
if (b.Length > 1) {
for (int i = 0; i < vLen; i++)
target[i] = a[i] / b[i];
} else {
// b == scalar
for (int i = 0; i < vLen; i++)
target[i] = a[i] / b[0];
}
} else {
// a == scalar
if (b.Length > 1) {
if (Math.Abs(a[0]) < 1E-12 /* == 0 */) {
target = GetScalarBuffer();
target[0] = 0.0;
} else {
target = GetVectorBuffer();
for (int i = 0; i < vLen; i++)
target[i] = a[0] / b[i];
}
} else {
// b == scalar
target = GetScalarBuffer();
target[0] = a[0] / b[0];
}
}
return target;
}
private void ReadNext(byte[] code, ref int pc, out byte op, out short s) {
op = code[pc++];
s = 0;
if (op == (byte)OpCodes.LoadVar) {
s = (short)(((short)code[pc] << 8) | (short)code[pc + 1]);
pc += 2;
}
}
}
}