1  using System;


2  using System.Collections.Generic;


3  using System.Linq;


4  using HeuristicLab.Common;


5  using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;


6  using HEAL.Attic;


7 


8  namespace HeuristicLab.Problems.DataAnalysis.Symbolic {


9  public abstract class Interpreter<T> where T : IAlgebraicType<T> {


10  public struct Instruction {


11  public byte opcode;


12  public ushort narg;


13  public int childIndex;


14  public double dblVal;


15  public object data; // any kind of data you want to store in instructions


16  public T value;


17  }


18 


19  public T Evaluate(Instruction[] code) {


20  for (int i = code.Length  1; i >= 0; i) {


21  var instr = code[i];


22  var c = instr.childIndex;


23  var n = instr.narg;


24 


25  switch (instr.opcode) {


26  case OpCodes.Variable: {


27  LoadVariable(instr);


28  break;


29  }


30  // case OpCodes.Constant: we initialize constants in Compile. The value never changes afterwards


31  case OpCodes.Add: {


32  instr.value.Assign(code[c].value);


33  for (int j = 1; j < n; ++j) {


34  instr.value.Add(code[c + j].value);


35  }


36  break;


37  }


38 


39  case OpCodes.Sub: {


40  if (n == 1) {


41  instr.value.AssignNeg(code[c].value);


42  } else {


43  instr.value.Assign(code[c].value);


44  for (int j = 1; j < n; ++j) {


45  instr.value.Sub(code[c + j].value);


46  }


47  }


48  break;


49  }


50 


51  case OpCodes.Mul: {


52  instr.value.Assign(code[c].value);


53  for (int j = 1; j < n; ++j) {


54  instr.value.Mul(code[c + j].value);


55  }


56  break;


57  }


58 


59  case OpCodes.Div: {


60  if (n == 1) {


61  instr.value.AssignInv(code[c].value);


62  } else {


63  instr.value.Assign(code[c].value);


64  for (int j = 1; j < n; ++j) {


65  instr.value.Div(code[c + j].value);


66  }


67  }


68  break;


69  }


70 


71  case OpCodes.Exp: {


72  instr.value.AssignExp(code[c].value);


73  break;


74  }


75 


76  case OpCodes.Log: {


77  instr.value.AssignLog(code[c].value);


78  break;


79  }


80  }


81  }


82  return code[0].value;


83  }


84 


85  protected Instruction[] Compile(ISymbolicExpressionTree tree) {


86  var root = tree.Root.GetSubtree(0).GetSubtree(0);


87  var code = new Instruction[root.GetLength()];


88  if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");


89  int c = 1, i = 0;


90  foreach (var node in root.IterateNodesBreadth()) {


91  if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");


92  code[i] = new Instruction {


93  opcode = OpCodes.MapSymbolToOpCode(node),


94  narg = (ushort)node.SubtreeCount,


95  childIndex = c


96  };


97  if (node is VariableTreeNode variable) {


98  InitializeTerminalInstruction(ref code[i], variable);


99  } else if (node is ConstantTreeNode constant) {


100  InitializeTerminalInstruction(ref code[i], constant);


101  } else {


102  InitializeInternalInstruction(ref code[i], node);


103  }


104  c += node.SubtreeCount;


105  ++i;


106  }


107  return code;


108  }


109 


110  protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant);


111  protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable);


112  protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node);


113 


114  protected abstract void LoadVariable(Instruction a);


115 


116  }


117 


118 


119  public class VectorEvaluator : Interpreter<DoubleVector> {


120  private const int BATCHSIZE = 128;


121  [ThreadStatic]


122  private Dictionary<string, double[]> cachedData;


123 


124  [ThreadStatic]


125  private IDataset dataset;


126 


127  [ThreadStatic]


128  private int rowIndex;


129 


130  [ThreadStatic]


131  private int[] rows;


132 


133  private void InitCache(IDataset dataset) {


134  this.dataset = dataset;


135  cachedData = new Dictionary<string, double[]>();


136  foreach (var v in dataset.DoubleVariables) {


137  cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();


138  }


139  }


140 


141  public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {


142  if (cachedData == null  this.dataset != dataset) {


143  InitCache(dataset);


144  }


145 


146  this.rows = rows;


147  var code = Compile(tree);


148  var remainingRows = rows.Length % BATCHSIZE;


149  var roundedTotal = rows.Length  remainingRows;


150 


151  var result = new double[rows.Length];


152 


153  for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {


154  Evaluate(code);


155  code[0].value.CopyTo(result, rowIndex, BATCHSIZE);


156  }


157 


158  if (remainingRows > 0) {


159  Evaluate(code);


160  code[0].value.CopyTo(result, roundedTotal, remainingRows);


161  }


162 


163  return result;


164  }


165 


166  protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {


167  instruction.dblVal = constant.Value;


168  instruction.value = new DoubleVector(BATCHSIZE);


169  instruction.value.AssignConstant(instruction.dblVal);


170  }


171 


172  protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {


173  instruction.dblVal = variable.Weight;


174  instruction.value = new DoubleVector(BATCHSIZE);


175  if (cachedData.ContainsKey(variable.VariableName)) {


176  instruction.data = cachedData[variable.VariableName];


177  } else {


178  instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();


179  cachedData[variable.VariableName] = (double[])instruction.data;


180  }


181  }


182 


183  protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {


184  instruction.value = new DoubleVector(BATCHSIZE);


185  }


186 


187  protected override void LoadVariable(Instruction a) {


188  var data = (double[])a.data;


189  for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) a.value[i  rowIndex] = data[rows[i]];


190  a.value.Scale(a.dblVal);


191  }


192  }


193 


194  public class VectorAutoDiffEvaluator : Interpreter<VectorDual<DoubleVector>> {


195  private const int BATCHSIZE = 128;


196  [ThreadStatic]


197  private Dictionary<string, double[]> cachedData;


198 


199  [ThreadStatic]


200  private IDataset dataset;


201 


202  [ThreadStatic]


203  private int rowIndex;


204 


205  [ThreadStatic]


206  private int[] rows;


207 


208  [ThreadStatic]


209  private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;


210 


211  private void InitCache(IDataset dataset) {


212  this.dataset = dataset;


213  cachedData = new Dictionary<string, double[]>();


214  foreach (var v in dataset.DoubleVariables) {


215  cachedData[v] = dataset.GetDoubleValues(v).ToArray();


216  }


217  }


218 


219  public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, out double[] fi, out double[,] jac) {


220  if (cachedData == null  this.dataset != dataset) {


221  InitCache(dataset);


222  }


223 


224  int nParams = parameterNodes.Length;


225  node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();


226  for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);


227 


228  var code = Compile(tree);


229 


230  var remainingRows = rows.Length % BATCHSIZE;


231  var roundedTotal = rows.Length  remainingRows;


232 


233  fi = new double[rows.Length];


234  jac = new double[rows.Length, nParams];


235 


236  this.rows = rows;


237 


238  for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {


239  Evaluate(code);


240  code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);


241 


242  // TRANSPOSE into JAC


243  // here we assume that we usually calculate the gradient over _all_ parameters (otherwise it would be better to propagate only relevant vectors in evaluation)


244  for (int j = 0; j < nParams; ++j)


245  code[0].value.Gradient[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); XXX CONTINUE HERE


246  }


247 


248  if (remainingRows > 0) {


249  Evaluate(code);


250  code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);


251  for (int j = 0; j < nParams; ++j)


252  code[0].value.Gradient[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);


253  }


254  }


255 


256  protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {


257  var zero = new DoubleVector(BATCHSIZE);


258  // var g = new Vector<DoubleVector>(node2paramIdx.Count);


259  // for (int i = 0; i < BATCHSIZE; ++i) g[i] = new DoubleVector(BATCHSIZE);


260  instruction.value = new VectorDual<DoubleVector>(zero, null);


261  }


262 


263  protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {


264  instruction.dblVal = constant.Value;


265 


266  var g_arr = new double[BATCHSIZE];


267  if (node2paramIdx.ContainsKey(constant)) {


268  for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;


269  }


270  var g = new DoubleVector(g_arr);


271 


272  instruction.value = new VectorDual<DoubleVector>(new DoubleVector(BATCHSIZE), new Vector<DoubleVector>(new[] { g })); // only a single column for the gradient


273  instruction.value.Value.AssignConstant(instruction.dblVal);


274  }


275 


276  protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {


277  double[] data;


278  if (cachedData.ContainsKey(variable.VariableName)) {


279  data = cachedData[variable.VariableName];


280  } else {


281  data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();


282  cachedData[variable.VariableName] = (double[])instruction.data;


283  }


284 


285  var paramIdx = 1;


286  if (node2paramIdx.ContainsKey(variable)) {


287  paramIdx = node2paramIdx[variable];


288  }


289 


290  instruction.dblVal = variable.Weight;


291  instruction.data = new object[] { data, paramIdx };


292 


293  var f = new DoubleVector(BATCHSIZE);


294  var g = new DoubleVector(BATCHSIZE);


295  instruction.value = new VectorDual<DoubleVector>(f, new Vector<DoubleVector>(new[] { g }));


296  }


297 


298  protected override void LoadVariable(Instruction a) {


299  var paramIdx = (int)((object[])a.data)[1];


300  var data = (double[])((object[])a.data)[0];


301 


302  for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) a.value.Value[i  rowIndex] = data[rows[i]];


303  a.value.Scale(a.dblVal);


304 


305  if (paramIdx >= 0) {


306  // update gradient with variable values


307  var g = a.value.Gradient[0];


308  for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) {


309  g[i] = data[rows[i]];


310  }


311  }


312  }


313  }


314 


315 


316  public class IntervalEvaluator : Interpreter<AlgebraicInterval> {


317  [ThreadStatic]


318  private Dictionary<string, AlgebraicInterval> intervals;


319 


320  public AlgebraicInterval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, AlgebraicInterval> intervals) {


321  this.intervals = intervals;


322  var code = Compile(tree);


323  Evaluate(code);


324  return code[0].value;


325  }


326 


327 


328  protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {


329  instruction.value = new AlgebraicInterval(0, 0);


330  }


331 


332 


333  protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {


334  instruction.dblVal = constant.Value;


335  instruction.value = new AlgebraicInterval(constant.Value, constant.Value);


336  }


337 


338  protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {


339  instruction.dblVal = variable.Weight;


340  instruction.value = intervals[variable.VariableName];


341  }


342 


343  protected override void LoadVariable(Instruction a) {


344  // nothing to do


345  }


346  }


347 


348  public interface IAlgebraicType<T> {


349  T Assign(T a); // assign this to same value as a (copy!)


350  T AssignNeg(T a); // set this to negative(a)


351  T AssignInv(T a); // set this to inv(a);


352  T Scale(double s); // scale this with s


353  T Add(T a); // add a to this


354  T Sub(T a); // subtract a from this


355  T Mul(T a); // multiply this with a


356  T Div(T a); // divide this by a


357  T AssignLog(T a); // set this to log a


358  T AssignExp(T a); // set this to exp(a)


359  T AssignSin(T a); // set this to sin(a)


360  T AssignCos(T a); // set this to cos(a)


361  T AssignIntPower(T a, int p);


362  T AssignIntRoot(T a, int r);


363  T Clone();


364  }


365 


366 


367  // a simple vector as an algebraic type


368  public class DoubleVector : IAlgebraicType<DoubleVector> {


369  private double[] arr;


370 


371  public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }


372 


373  public DoubleVector(int length) {


374  arr = new double[length];


375  }


376 


377  /// <summary>


378  ///


379  /// </summary>


380  /// <param name="arr">array is not copied</param>


381  public DoubleVector(double[] arr) {


382  this.arr = arr;


383  }


384 


385  public DoubleVector(int length, double constantValue) : this(length) {


386  for (int i = 0; i < length; ++i) arr[i] = constantValue;


387  }


388 


389  public void CopyTo(double[] dest, int idx, int length) {


390  Array.Copy(arr, 0, dest, idx, length);


391  }


392 


393  public void CopyRowTo(double[,] dest, int row) {


394  for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];


395  }


396  internal void CopyColumnTo(double[,] dest, int column, int row, int len) {


397  for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];


398  }


399 


400 


401  public DoubleVector Add(DoubleVector a) {


402  for (int i = 0; i < arr.Length; ++i) {


403  arr[i] += a.arr[i];


404  }


405  return this;


406  }


407 


408  public void AssignConstant(double constantValue) {


409  for (int i = 0; i < arr.Length; ++i) {


410  arr[i] = constantValue;


411  }


412  }


413 


414  public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }


415  public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }


416  public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }


417  public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }


418  public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }


419  public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }


420  public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }


421  public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }


422  public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }


423  public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }


424  public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }


425  public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }


426  public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }


427 


428  public DoubleVector Clone() {


429  var v = new DoubleVector(this.arr.Length);


430  Array.Copy(arr, v.arr, v.arr.Length);


431  return v;


432  }


433 


434  public void CopyFrom(double[] data, int rowIndex) {


435  Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length  rowIndex));


436  }


437 


438  }


439 


440 


441  // vectors of algebraic types


442  public class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {


443  private T[] elems;


444 


445  public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }


446 


447  public int Length => elems.Length;


448 


449  private Vector() { }


450 


451  public Vector(int len) {


452  elems = new T[len];


453  }


454 


455  /// <summary>


456  ///


457  /// </summary>


458  /// <param name="elems">The array is copied</param>


459  public Vector(T[] elems) {


460  this.elems = new T[elems.Length];


461  for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }


462  }


463 


464  /// <summary>


465  ///


466  /// </summary>


467  /// <param name="elems">Array is not copied!</param>


468  /// <returns></returns>


469  public Vector<T> FromArray(T[] elems) {


470  var v = new Vector<T>();


471  v.elems = elems;


472  return v;


473  }


474 


475  public void CopyTo(T[] dest) {


476  if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");


477  Array.Copy(elems, dest, dest.Length);


478  }


479 


480  public Vector<T> Clone() {


481  return new Vector<T>(elems);


482  }


483 


484  public Vector<T> Concat(Vector<T> other) {


485  var oldLen = Length;


486  Array.Resize(ref this.elems, oldLen + other.Length);


487  for (int i = oldLen; i < Length; i++) {


488  elems[i] = other.elems[i  oldLen].Clone();


489  }


490  return this;


491  }


492 


493  public Vector<T> Add(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }


494  public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }


495  public Vector<T> AssignCos(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }


496  public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }


497  public Vector<T> AssignIntPower(Vector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }


498  public Vector<T> AssignIntRoot(Vector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }


499  public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }


500  public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }


501  public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }


502  public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }


503  public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }


504  public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }


505  public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }


506  public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }


507  public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }


508  }


509 


510 


511  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {


512  private double low;


513  private double high;


514 


515  public double LowerBound => low;


516  public double UpperBound => high;


517 


518  public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }


519 


520  public AlgebraicInterval(double low, double high) {


521  this.low = low;


522  this.high = high;


523  }


524 


525  public AlgebraicInterval Add(AlgebraicInterval a) {


526  low += a.low;


527  high += a.high;


528  return this;


529  }


530 


531  public AlgebraicInterval Assign(AlgebraicInterval a) {


532  low = a.low;


533  high = a.high;


534  return this;


535  }


536 


537  public AlgebraicInterval AssignCos(AlgebraicInterval a) {


538  throw new NotImplementedException();


539  }


540 


541  public AlgebraicInterval Div(AlgebraicInterval a) {


542  if (Contains(0.0)) {


543  if (a.low.IsAlmost(0.0))


544  Mul(new AlgebraicInterval(1.0 / a.high, double.PositiveInfinity));


545  else if (a.high.IsAlmost(0.0))


546  Mul(new AlgebraicInterval(double.NegativeInfinity, 1.0 / a.low));


547  else {


548  low = double.NegativeInfinity;


549  high = double.PositiveInfinity;


550  }


551  } else {


552  Mul(new AlgebraicInterval(1.0 / a.high, 1.0 / a.low));


553  }


554  return this;


555  }


556 


557  public AlgebraicInterval AssignExp(AlgebraicInterval a) {


558  low = Math.Exp(a.low);


559  high = Math.Exp(a.high);


560  return this;


561  }


562 


563  public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {


564  throw new NotImplementedException();


565  }


566 


567  public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {


568  throw new NotImplementedException();


569  }


570 


571  public AlgebraicInterval AssignInv(AlgebraicInterval a) {


572  low = 1.0;


573  high = 1.0;


574  return Div(a);


575  }


576 


577  public AlgebraicInterval AssignLog(AlgebraicInterval a) {


578  low = Math.Log(a.low);


579  high = Math.Log(a.high);


580  return this;


581  }


582 


583  public AlgebraicInterval Mul(AlgebraicInterval a) {


584  double v1 = low * a.low;


585  double v2 = low * a.high;


586  double v3 = high * a.low;


587  double v4 = high * a.high;


588 


589  low = Math.Min(Math.Min(v1, v2), Math.Min(v3, v4));


590  high = Math.Max(Math.Max(v1, v2), Math.Max(v3, v4));


591  return this;


592  }


593 


594  public AlgebraicInterval AssignNeg(AlgebraicInterval a) {


595  low = 0;


596  high = 0;


597  return Sub(a);


598  }


599 


600  public AlgebraicInterval Scale(double s) {


601  if (s < 0) {


602  low = s * high;


603  high = s * low;


604  } else {


605  low = s * low;


606  high = s * high;


607  }


608  return this;


609  }


610 


611  public AlgebraicInterval AssignSin(AlgebraicInterval a) {


612  throw new NotImplementedException();


613  }


614 


615  public AlgebraicInterval Sub(AlgebraicInterval a) {


616  low = a.high;


617  high = a.low;


618  return this;


619  }


620 


621  public AlgebraicInterval Clone() {


622  return new AlgebraicInterval(low, high);


623  }


624 


625  public bool Contains(double val) {


626  return low <= val && val <= high;


627  }


628  }


629 


630  public class Dual<V> : IAlgebraicType<Dual<V>>


631  where V : IAlgebraicType<V> {


632  private V v;


633  private V dv;


634 


635  public Dual(V v, V dv) {


636  this.v = v;


637  this.dv = dv;


638  }


639 


640  public Dual<V> Clone() {


641  return new Dual<V>(v.Clone(), dv.Clone());


642  }


643 


644  public Dual<V> Add(Dual<V> a) {


645  v.Add(a.v);


646  dv.Add(a.dv);


647  return this;


648  }


649 


650  public Dual<V> Assign(Dual<V> a) {


651  v.Assign(a.v);


652  dv.Assign(a.dv);


653  return this;


654  }


655 


656  public Dual<V> AssignCos(Dual<V> a) {


657  v.AssignCos(a.v);


658  dv.AssignNeg(dv.AssignSin(a.v));


659  return this;


660  }


661 


662  public Dual<V> Div(Dual<V> a) {


663  throw new NotImplementedException();


664  }


665 


666  public Dual<V> AssignExp(Dual<V> a) {


667  v.AssignExp(a.v);


668  dv.Assign(a.dv).Mul(v); // exp(f(x)) = exp(f(x))*f(x)'


669  return this;


670  }


671 


672  public Dual<V> AssignIntPower(Dual<V> a, int p) {


673  throw new NotImplementedException();


674  }


675 


676  public Dual<V> AssignIntRoot(Dual<V> a, int r) {


677  throw new NotImplementedException();


678  }


679 


680  public Dual<V> AssignInv(Dual<V> a) {


681  throw new NotImplementedException();


682  }


683 


684  public Dual<V> AssignLog(Dual<V> a) {


685  v.AssignLog(a.v);


686  dv.Assign(a.dv).Div(a.v); // log(x)' = 1/f(x) * f(x)'


687  return this;


688  }


689 


690  public Dual<V> Mul(Dual<V> a) {


691  // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);


692 


693  V t1 = default(V);


694  t1.Assign(a.dv).Mul(v);


695 


696  var t2 = default(V);


697  t2.Assign(dv).Mul(a.v);


698 


699  dv.Assign(t1).Add(t2);


700 


701  v.Mul(a.v);


702  return this;


703  }


704 


705  public Dual<V> AssignNeg(Dual<V> a) {


706  throw new NotImplementedException();


707  }


708 


709  public Dual<V> Scale(double s) {


710  v.Scale(s);


711  dv.Scale(s);


712  return this;


713  }


714 


715  public Dual<V> AssignSin(Dual<V> a) {


716  throw new NotImplementedException();


717  }


718 


719  public Dual<V> Sub(Dual<V> a) {


720  throw new NotImplementedException();


721  }


722  }


723 


724  public class VectorDual<V> : IAlgebraicType<VectorDual<V>> where V : IAlgebraicType<V> {


725  private V v;


726  public V Value => v;


727 


728  private Vector<V> dv;


729  public Vector<V> Gradient => dv; // columnoriented (dv[0] is the vector for the first parameter)


730 


731 


732  public VectorDual(V v, Vector<V> dv) {


733  this.v = v;


734  this.dv = dv;


735  }


736 


737  public VectorDual<V> Clone() {


738  return new VectorDual<V>(v.Clone(), dv.Clone());


739  }


740 


741  public VectorDual<V> Add(VectorDual<V> a) {


742  v.Add(a.v);


743  dv.Concat(a.dv);


744  return this;


745  }


746 


747  public VectorDual<V> Assign(VectorDual<V> a) {


748  v.Assign(a.v);


749  dv = a.dv.Clone();


750  return this;


751  }


752 


753  public VectorDual<V> AssignCos(VectorDual<V> a) {


754  v.AssignCos(a.v);


755  V t = default(V);


756  t.AssignNeg(t.AssignSin(a.v));


757  dv = a.dv.Clone().Scale(t);


758  return this;


759  }


760 


761  public VectorDual<V> AssignExp(VectorDual<V> a) {


762  v.AssignExp(a.v);


763  dv = a.dv.Clone().Scale(v);


764  return this;


765  }


766 


767  public VectorDual<V> AssignIntPower(VectorDual<V> a, int p) {


768  throw new NotImplementedException();


769  }


770 


771  public VectorDual<V> AssignIntRoot(VectorDual<V> a, int r) {


772  throw new NotImplementedException();


773  }


774 


775  public VectorDual<V> AssignInv(VectorDual<V> a) {


776  throw new NotImplementedException();


777  }


778 


779  public VectorDual<V> AssignLog(VectorDual<V> a) {


780  throw new NotImplementedException();


781  }


782 


783  public VectorDual<V> AssignNeg(VectorDual<V> a) {


784  throw new NotImplementedException();


785  }


786 


787  public VectorDual<V> AssignSin(VectorDual<V> a) {


788  throw new NotImplementedException();


789  }


790 


791  public VectorDual<V> Div(VectorDual<V> a) {


792  throw new NotImplementedException();


793  }


794 


795  public VectorDual<V> Mul(VectorDual<V> a) {


796  // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);


797 


798  var t1 = a.dv.Clone().Scale(v);


799  var t2 = dv.Clone().Scale(a.v);


800  dv.Assign(t1).Concat(t2); // TODO: CHECK ORDER


801 


802  v.Mul(a.v);


803  return this;


804  }


805 


806  public VectorDual<V> Scale(double s) {


807  v.Scale(s);


808  dv.Scale(s);


809  return this;


810  }


811 


812  public VectorDual<V> Sub(VectorDual<V> a) {


813  v.Sub(a.v);


814  dv.Concat(a.dv.AssignNeg(a.dv));


815  return this;


816  }


817  }


818  }

