[16674]  1  using System;


 2  using System.Collections.Generic;


[16695]  3  using System.Diagnostics;


[16674]  4  using System.Linq;


 5  using HeuristicLab.Common;


 6  using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;


 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  }


[16694]  30  case OpCodes.Constant: { break; } // we initialize constants in Compile. The value never changes afterwards


[16674]  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  }


[16694]  70  case OpCodes.Square: {


 71  instr.value.AssignIntPower(code[c].value, 2);


 72  break;


 73  }


[16727]  74  case OpCodes.SquareRoot: {


 75  instr.value.AssignIntRoot(code[c].value, 2);


 76  break;


 77  }


 78  case OpCodes.Cube: {


 79  instr.value.AssignIntPower(code[c].value, 3);


 80  break;


 81  }


 82  case OpCodes.CubeRoot: {


 83  instr.value.AssignIntRoot(code[c].value, 3);


 84  break;


 85  }


[16674]  86  case OpCodes.Exp: {


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


 88  break;


 89  }


 90  case OpCodes.Log: {


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


 92  break;


 93  }


[16727]  94  case OpCodes.Sin: {


 95  instr.value.AssignSin(code[c].value);


 96  break;


 97  }


 98  case OpCodes.Cos: {


 99  instr.value.AssignCos(code[c].value);


 100  break;


 101  }


 102  case OpCodes.Absolute: {


 103  instr.value.AssignAbs(code[c].value);


 104  break;


 105  }


 106  case OpCodes.AnalyticQuotient: {


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


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


 109  var t = instr.value.One;


 110  t.Add(code[c + j].value.Clone().IntPower(2));


 111  instr.value.Div(t.IntRoot(2));


 112  }


 113  break;


 114  }


[17131]  115  case OpCodes.Tanh: {


 116  instr.value.AssignTanh(code[c].value);


 117  break;


 118  }


 119 


 120  default: throw new ArgumentException($"Unknown opcode {instr.opcode}");


[16674]  121  }


 122  }


 123  return code[0].value;


 124  }


 125 


 126  protected Instruction[] Compile(ISymbolicExpressionTree tree) {


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


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


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


 130  int c = 1, i = 0;


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


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


 133  code[i] = new Instruction {


 134  opcode = OpCodes.MapSymbolToOpCode(node),


 135  narg = (ushort)node.SubtreeCount,


 136  childIndex = c


 137  };


 138  if (node is VariableTreeNode variable) {


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


 140  } else if (node is ConstantTreeNode constant) {


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


 142  } else {


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


 144  }


 145  c += node.SubtreeCount;


 146  ++i;


 147  }


 148  return code;


 149  }


 150 


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


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


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


 154 


 155  protected abstract void LoadVariable(Instruction a);


 156 


 157  }


 158 


 159 


[16695]  160  public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {


[16674]  161  private const int BATCHSIZE = 128;


 162  [ThreadStatic]


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


 164 


 165  [ThreadStatic]


 166  private IDataset dataset;


 167 


 168  [ThreadStatic]


 169  private int rowIndex;


 170 


 171  [ThreadStatic]


 172  private int[] rows;


 173 


 174  private void InitCache(IDataset dataset) {


 175  this.dataset = dataset;


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


 177  foreach (var v in dataset.DoubleVariables) {


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


 179  }


 180  }


 181 


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


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


 184  InitCache(dataset);


 185  }


 186 


 187  this.rows = rows;


 188  var code = Compile(tree);


 189  var remainingRows = rows.Length % BATCHSIZE;


 190  var roundedTotal = rows.Length  remainingRows;


 191 


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


 193 


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


 195  Evaluate(code);


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


 197  }


 198 


 199  if (remainingRows > 0) {


 200  Evaluate(code);


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


 202  }


 203 


 204  return result;


 205  }


 206 


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


 208  instruction.dblVal = constant.Value;


[16695]  209  instruction.value = new AlgebraicDoubleVector(BATCHSIZE);


[16674]  210  instruction.value.AssignConstant(instruction.dblVal);


 211  }


 212 


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


 214  instruction.dblVal = variable.Weight;


[16695]  215  instruction.value = new AlgebraicDoubleVector(BATCHSIZE);


[16674]  216  if (cachedData.ContainsKey(variable.VariableName)) {


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


 218  } else {


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


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


 221  }


 222  }


 223 


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


[16695]  225  instruction.value = new AlgebraicDoubleVector(BATCHSIZE);


[16674]  226  }


 227 


 228  protected override void LoadVariable(Instruction a) {


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


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


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


 232  }


 233  }


 234 


[16695]  235  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {


[16674]  236  private const int BATCHSIZE = 128;


 237  [ThreadStatic]


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


 239 


 240  [ThreadStatic]


 241  private IDataset dataset;


 242 


 243  [ThreadStatic]


 244  private int rowIndex;


 245 


 246  [ThreadStatic]


 247  private int[] rows;


 248 


 249  [ThreadStatic]


 250  private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;


 251 


 252  private void InitCache(IDataset dataset) {


 253  this.dataset = dataset;


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


 255  foreach (var v in dataset.DoubleVariables) {


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


 257  }


 258  }


 259 


[16727]  260  /// <summary>


 261  ///


 262  /// </summary>


 263  /// <param name="tree"></param>


 264  /// <param name="dataset"></param>


 265  /// <param name="rows"></param>


 266  /// <param name="parameterNodes"></param>


 267  /// <param name="fi">Function output. Must be allocated by the caller.</param>


 268  /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>


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


[16674]  270  if (cachedData == null  this.dataset != dataset) {


 271  InitCache(dataset);


 272  }


 273 


 274  int nParams = parameterNodes.Length;


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


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


 277 


 278  var code = Compile(tree);


 279 


 280  var remainingRows = rows.Length % BATCHSIZE;


 281  var roundedTotal = rows.Length  remainingRows;


 282 


 283  this.rows = rows;


 284 


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


 286  Evaluate(code);


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


 288 


 289  // TRANSPOSE into JAC


[16682]  290  var g = code[0].value.Gradient;


 291  for (int j = 0; j < nParams; ++j) {


[16744]  292  if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {


[16727]  293  v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);


 294  } else {


 295  for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;


 296  }


[16682]  297  }


[16674]  298  }


 299 


 300  if (remainingRows > 0) {


 301  Evaluate(code);


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


[16682]  303 


 304  var g = code[0].value.Gradient;


[16674]  305  for (int j = 0; j < nParams; ++j)


[16727]  306  if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {


 307  v.CopyColumnTo(jac, j, roundedTotal, remainingRows);


 308  } else {


 309  for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;


 310  }


[16674]  311  }


 312  }


 313 


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


[16695]  315  var zero = new AlgebraicDoubleVector(BATCHSIZE);


 316  instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);


[16674]  317  }


 318 


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


 320  var g_arr = new double[BATCHSIZE];


[16682]  321  if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {


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


[16695]  323  var g = new AlgebraicDoubleVector(g_arr);


 324  instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient


[16682]  325  } else {


[16695]  326  instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));


[16674]  327  }


 328 


[16682]  329  instruction.dblVal = constant.Value;


[16674]  330  instruction.value.Value.AssignConstant(instruction.dblVal);


 331  }


 332 


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


 334  double[] data;


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


 336  data = cachedData[variable.VariableName];


 337  } else {


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


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


 340  }


 341 


 342  var paramIdx = 1;


 343  if (node2paramIdx.ContainsKey(variable)) {


 344  paramIdx = node2paramIdx[variable];


[16695]  345  var f = new AlgebraicDoubleVector(BATCHSIZE);


 346  var g = new AlgebraicDoubleVector(BATCHSIZE);


 347  instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);


[16682]  348  } else {


[16695]  349  var f = new AlgebraicDoubleVector(BATCHSIZE);


 350  instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);


[16674]  351  }


 352 


 353  instruction.dblVal = variable.Weight;


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


 355  }


 356 


 357  protected override void LoadVariable(Instruction a) {


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


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


 360 


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


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


 363 


 364  if (paramIdx >= 0) {


 365  // update gradient with variable values


[16682]  366  var g = a.value.Gradient.Elements[paramIdx];


[16674]  367  for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) {


[16727]  368  g[i  rowIndex] = data[rows[i]];


[16674]  369  }


 370  }


 371  }


 372  }


 373 


 374 


[16693]  375  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {


[16674]  376  [ThreadStatic]


[16831]  377  private IDictionary<string, Interval> intervals;


[16674]  378 


[16831]  379  public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals) {


[16674]  380  this.intervals = intervals;


 381  var code = Compile(tree);


 382  Evaluate(code);


[16682]  383  return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);


[16674]  384  }


 385 


[16831]  386  public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {


[16682]  387  this.intervals = intervals;


 388  var code = Compile(tree);


 389  Evaluate(code);


 390  lowerGradient = new double[paramNodes.Length];


 391  upperGradient = new double[paramNodes.Length];


 392  var l = code[0].value.LowerBound;


 393  var u = code[0].value.UpperBound;


 394  for (int i = 0; i < paramNodes.Length; ++i) {


[16694]  395  if (paramNodes[i] == null) continue;


[16738]  396  if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;


 397  if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;


[16682]  398  }


 399  return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);


 400  }


[16674]  401 


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


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


 404  }


 405 


 406 


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


 408  instruction.dblVal = constant.Value;


[16682]  409  instruction.value = new AlgebraicInterval(


[16693]  410  new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0),


 411  new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key


[16682]  412  );


[16674]  413  }


 414 


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


 416  instruction.dblVal = variable.Weight;


[16682]  417  instruction.value = new AlgebraicInterval(


[16693]  418  low: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound), // bounds change by variable value d/dc (c I(var)) = I(var)


 419  high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)


[16682]  420  );


[16674]  421  }


 422 


 423  protected override void LoadVariable(Instruction a) {


 424  // nothing to do


 425  }


 426  }


 427 


 428  public interface IAlgebraicType<T> {


[16693]  429  T Zero { get; }


[16727]  430  T One { get; }


[16693]  431 


[16682]  432  T AssignAbs(T a); // set this to assign abs(a)


[16674]  433  T Assign(T a); // assign this to same value as a (copy!)


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


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


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


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


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


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


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


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


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


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


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


[17131]  445  T AssignTanh(T a); // set this to tanh(a)


[16674]  446  T AssignIntPower(T a, int p);


 447  T AssignIntRoot(T a, int r);


[16682]  448  T AssignSgn(T a); // set this to sign(a)


[16674]  449  T Clone();


 450  }


 451 


[16682]  452  public static class Algebraic {


[16695]  453  public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }


 454  public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }


 455  public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }


 456  public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }


 457  public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }


 458  public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }


 459  public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }


 460  public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }


 461  public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }


 462  public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }


[16674]  463 


[16682]  464  public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {


 465  // ((a + b) + abs(b  a)) / 2


 466  return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);


 467  }


 468  public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {


 469  // ((a + b)  abs(a  b)) / 2


 470  return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);


 471  }


 472  }


 473 


 474 


 475  // algebraic type wrapper for a double value


[16695]  476  [DebuggerDisplay("{Value}")]


[16693]  477  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {


 478  public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }


 479  public static implicit operator double(AlgebraicDouble value) { return value.Value; }


[16682]  480  public double Value;


 481 


[16695]  482  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16693]  483  public AlgebraicDouble Zero => new AlgebraicDouble(0.0);


[16727]  484  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 485  public AlgebraicDouble One => new AlgebraicDouble(1.0);


[16693]  486  public AlgebraicDouble() { }


 487  public AlgebraicDouble(double value) { this.Value = value; }


 488  public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }


 489  public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }


 490  public AlgebraicDouble Sub(AlgebraicDouble a) { Value = a.Value; return this; }


 491  public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }


 492  public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }


 493  public AlgebraicDouble Scale(double s) { Value *= s; return this; }


 494  public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }


 495  public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = a.Value; return this; }


 496  public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }


 497  public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }


[17131]  498  public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }


[16744]  499  public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = a.Value <= 0 ? double.NegativeInfinity : Math.Log(a.Value); return this; } // alternative definiton of log to prevent NaN


[16693]  500  public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }


 501  public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }


 502  public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }


 503  public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }


[16744]  504  public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }


[16693]  505  public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }


[16695]  506 


 507  public override string ToString() {


 508  return Value.ToString();


 509  }


[16682]  510  }


 511 


[16674]  512  // a simple vector as an algebraic type


[16695]  513  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]


 514  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {


[16674]  515  private double[] arr;


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


[16693]  517  public int Length => arr.Length;


 518 


[16695]  519  public AlgebraicDoubleVector(int length) { arr = new double[length]; }


[16674]  520 


[16695]  521  public AlgebraicDoubleVector() { }


[16682]  522 


[16674]  523  /// <summary>


 524  ///


 525  /// </summary>


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


[16695]  527  public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }


[16674]  528 


[16695]  529  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 530  public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)


[16727]  531  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 532  public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)


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


 534  public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }


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


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


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


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


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


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


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


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


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


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


[17131]  545  public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }


[16695]  546  public AlgebraicDoubleVector AssignIntPower(AlgebraicDoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }


 547  public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], 1.0 / r); } return this; }


 548  public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }


 549  public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }


[16674]  550 


[16695]  551  public AlgebraicDoubleVector Clone() {


 552  var v = new AlgebraicDoubleVector(this.arr.Length);


[16693]  553  Array.Copy(arr, v.arr, v.arr.Length);


 554  return v;


[16674]  555  }


 556 


[16727]  557  public AlgebraicDoubleVector AssignConstant(double constantValue) {


[16674]  558  for (int i = 0; i < arr.Length; ++i) {


 559  arr[i] = constantValue;


 560  }


[16727]  561  return this;


[16674]  562  }


 563 


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


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


[16674]  566  }


 567 


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


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


 570  }


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


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


 573  }


 574 


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


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


 577  }


[16695]  578 


 579  public override string ToString() {


 580  return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";


 581  }


[16674]  582  }


 583 


[16727]  584 


 585  /*


[16674]  586  // vectors of algebraic types


[16727]  587  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {


[16674]  588  private T[] elems;


 589 


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


 591 


 592  public int Length => elems.Length;


 593 


[16695]  594  private AlgebraicVector() { }


[16693]  595 


[16695]  596  public AlgebraicVector(int len) { elems = new T[len]; }


[16674]  597 


 598  /// <summary>


 599  ///


 600  /// </summary>


[16693]  601  /// <param name="elems">The array is copied (elementwise clone)</param>


[16695]  602  public AlgebraicVector(T[] elems) {


[16674]  603  this.elems = new T[elems.Length];


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


 605  }


 606 


 607  /// <summary>


 608  ///


 609  /// </summary>


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


 611  /// <returns></returns>


[16695]  612  public AlgebraicVector<T> FromArray(T[] elems) {


 613  var v = new AlgebraicVector<T>();


[16674]  614  v.elems = elems;


 615  return v;


 616  }


 617 


 618  public void CopyTo(T[] dest) {


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


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


 621  }


 622 


[16695]  623  public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }


[16674]  624 


 625 


[16695]  626  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 627  public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);


[16727]  628  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 629  public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


 645  public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }


 646  public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }


[16674]  647  }


 648 


[16727]  649  */


[16674]  650 


[16727]  651 


[16682]  652  /// <summary>


 653  /// A sparse vector of algebraic types. Elements are accessed via a key of type K


 654  /// </summary>


 655  /// <typeparam name="K">Key type</typeparam>


 656  /// <typeparam name="T">Element type</typeparam>


[16695]  657  [DebuggerDisplay("SparseVector: {ToString()}")]


 658  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {


 659  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16682]  660  private Dictionary<K, T> elems;


 661  public IReadOnlyDictionary<K, T> Elements => elems;


 662 


[16693]  663 


[16695]  664  public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {


[16682]  665  elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());


 666  }


 667 


[16693]  668  /// <summary>


 669  ///


 670  /// </summary>


 671  /// <param name="keys"></param>


 672  /// <param name="values">values are cloned</param>


[16695]  673  public AlgebraicSparseVector(K[] keys, T[] values) {


[16682]  674  if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");


 675  elems = new Dictionary<K, T>(keys.Length);


 676  for (int i = 0; i < keys.Length; ++i) {


[16693]  677  elems.Add(keys[i], values[i].Clone());


[16682]  678  }


 679  }


 680 


[16695]  681  public AlgebraicSparseVector() {


[16682]  682  this.elems = new Dictionary<K, T>();


 683  }


 684 


[16695]  685  // keep only elements from a


 686  private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {


 687  // remove elems from this which don't occur in a


 688  List<K> keysToRemove = new List<K>();


 689  foreach (var kvp in elems) {


 690  if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);


 691  }


 692  foreach (var o in keysToRemove) elems.Remove(o); // > zero


[16682]  693 


[16695]  694  foreach (var kvp in a.elems) {


 695  if (elems.TryGetValue(kvp.Key, out T value))


 696  mapAssign(kvp.Value, value);


 697  else


 698  elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));


 699  }


 700  }


[16682]  701 


[16695]  702  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 703  public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();


[16727]  704  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 705  public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();


[16695]  706 


 707  public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }


 708  public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }


 709 


 710  public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }


 711  public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }


 712  public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }


 713  public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }


 714  public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }


 715  public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }


 716  public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }


 717  public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }


 718  public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }


[17131]  719  public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }


[16695]  720  public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }


 721  public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }


[16727]  722  public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {


 723  foreach (var kvp in a.elems) {


 724  if (elems.TryGetValue(kvp.Key, out T value))


 725  value.Add(kvp.Value);


 726  else


 727  elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a


 728  }


 729  return this;


 730  }


[16695]  731 


[16727]  732  public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {


 733  foreach (var kvp in a.elems) {


 734  if (elems.TryGetValue(kvp.Key, out T value))


 735  value.Sub(kvp.Value);


 736  else


 737  elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0  a


 738  }


 739  return this;


 740  }


 741 


 742  public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {


 743  var keys = elems.Keys.ToArray();


 744  foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0


 745  foreach (var kvp in a.elems) {


 746  if (elems.TryGetValue(kvp.Key, out T value))


 747  value.Mul(kvp.Value); // this * a


 748  }


 749  return this;


 750  }


 751 


 752  public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {


 753  return Mul(a.Clone().Inv());


 754  }


 755 


[16695]  756  public AlgebraicSparseVector<K, T> Clone() {


 757  return new AlgebraicSparseVector<K, T>(this);


[16682]  758  }


[16695]  759 


 760  public override string ToString() {


 761  return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";


 762  }


[16682]  763  }


 764 


[16695]  765  [DebuggerDisplay("[{low.Value}..{high.Value}]")]


[16674]  766  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {


[16695]  767  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16693]  768  private MultivariateDual<AlgebraicDouble> low;


[16695]  769  public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();


 770 


 771  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16693]  772  private MultivariateDual<AlgebraicDouble> high;


 773  public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();


[16674]  774 


[16693]  775 


[16674]  776  public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }


 777 


[16693]  778  public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {


[16682]  779  this.low = low.Clone();


 780  this.high = high.Clone();


 781  }


 782 


[16674]  783  public AlgebraicInterval(double low, double high) {


[16693]  784  this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));


 785  this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));


[16674]  786  }


 787 


[16695]  788  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16693]  789  public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);


[17200]  790  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16727]  791  public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);


[17200]  792 


[16674]  793  public AlgebraicInterval Add(AlgebraicInterval a) {


[16682]  794  low.Add(a.low);


 795  high.Add(a.high);


[16674]  796  return this;


 797  }


 798 


[16693]  799  public AlgebraicInterval Mul(AlgebraicInterval a) {


 800  var v1 = low.Clone().Mul(a.low);


 801  var v2 = low.Clone().Mul(a.high);


 802  var v3 = high.Clone().Mul(a.low);


 803  var v4 = high.Clone().Mul(a.high);


 804 


 805  low = Algebraic.Min(Algebraic.Min(v1, v2), Algebraic.Min(v3, v4));


 806  high = Algebraic.Max(Algebraic.Max(v1, v2), Algebraic.Max(v3, v4));


 807  return this;


 808  }


 809 


[16674]  810  public AlgebraicInterval Assign(AlgebraicInterval a) {


 811  low = a.low;


 812  high = a.high;


 813  return this;


 814  }


 815 


 816  public AlgebraicInterval AssignCos(AlgebraicInterval a) {


[16693]  817  return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));


[16674]  818  }


 819 


 820  public AlgebraicInterval Div(AlgebraicInterval a) {


[16682]  821  if (a.Contains(0.0)) {


[17200]  822  if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {


[16693]  823  low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);


 824  high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);


[17200]  825  } else if (a.low.Value.Value == 0)


[16693]  826  Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));


[16682]  827  else


[16693]  828  Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));


[16674]  829  } else {


[16682]  830  Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low


[16674]  831  }


 832  return this;


 833  }


 834 


 835  public AlgebraicInterval AssignExp(AlgebraicInterval a) {


[16682]  836  low.AssignExp(a.low);


 837  high.AssignExp(a.high);


[16674]  838  return this;


 839  }


 840 


[17131]  841  // tanh is a bijective function


 842  public AlgebraicInterval AssignTanh(AlgebraicInterval a) {


 843  low.AssignTanh(a.low);


 844  high.AssignTanh(a.high);


 845  return this;


 846  }


 847 


[16674]  848  public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {


[16693]  849  if (p < 0) { // x^3 == 1/(x^3)


 850  AssignIntPower(a, p);


 851  return AssignInv(this);


[16744]  852  } else if (p == 0) {


 853  if (a.Contains(0.0)) {


 854  // => 0^0 = 0 ; might not be relevant


 855  low = new MultivariateDual<AlgebraicDouble>(0.0);


 856  high = new MultivariateDual<AlgebraicDouble>(1.0);


 857  return this;


 858  } else {


 859  // => 1


 860  low = new MultivariateDual<AlgebraicDouble>(1.0);


 861  high = new MultivariateDual<AlgebraicDouble>(1.0);


 862  return this;


 863  }


 864  } else if (p == 1) return this;


 865  else if (p % 2 == 0) {


[16693]  866  // p is even => interval must be positive


[16744]  867  if (a.Contains(0.0)) {


 868  low = new MultivariateDual<AlgebraicDouble>(0.0);


 869  high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));


[16693]  870  } else {


[16727]  871  var lowPower = a.low.IntPower(p);


 872  var highPower = a.high.IntPower(p);


[16693]  873  low = Algebraic.Min(lowPower, highPower);


 874  high = Algebraic.Max(lowPower, highPower);


 875  }


[16744]  876  } else {


 877  // p is uneven


 878  if (a.Contains(0.0)) {


 879  low.AssignIntPower(a.low, p);


 880  high.AssignIntPower(a.high, p);


 881  } else {


 882  var lowPower = a.low.IntPower(p);


 883  var highPower = a.high.IntPower(p);


 884  low = Algebraic.Min(lowPower, highPower);


 885  high = Algebraic.Max(lowPower, highPower);


 886  }


[16693]  887  }


[16744]  888  return this;


[16674]  889  }


 890 


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


[16693]  892  if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }


 893  if (r == 1) return this;


 894  if (r < 0) {


 895  // x^ (1/2) = 1 / (x^(1/2))


 896  AssignIntRoot(a, r);


 897  return AssignInv(this);


 898  } else {


 899  // root only defined for positive arguments


 900  if (a.LowerBound.Value.Value < 0) {


 901  low = new MultivariateDual<AlgebraicDouble>(double.NaN);


 902  high = new MultivariateDual<AlgebraicDouble>(double.NaN);


 903  return this;


 904  } else {


 905  low.AssignIntRoot(a.low, r);


 906  high.AssignIntRoot(a.high, r);


 907  return this;


 908  }


 909  }


[16674]  910  }


 911 


 912  public AlgebraicInterval AssignInv(AlgebraicInterval a) {


[16693]  913  low = new MultivariateDual<AlgebraicDouble>(1.0);


 914  high = new MultivariateDual<AlgebraicDouble>(1.0);


[16674]  915  return Div(a);


 916  }


 917 


 918  public AlgebraicInterval AssignLog(AlgebraicInterval a) {


[16682]  919  low.AssignLog(a.low);


 920  high.AssignLog(a.high);


[16674]  921  return this;


 922  }


 923 


[16693]  924  public AlgebraicInterval AssignNeg(AlgebraicInterval a) {


 925  low.AssignNeg(a.high);


 926  high.AssignNeg(a.low);


[16674]  927  return this;


 928  }


 929 


 930  public AlgebraicInterval Scale(double s) {


[16682]  931  low.Scale(s);


 932  high.Scale(s);


[16674]  933  if (s < 0) {


[16682]  934  var t = low;


 935  low = high;


 936  high = t;


[16674]  937  }


 938  return this;


 939  }


 940 


 941  public AlgebraicInterval AssignSin(AlgebraicInterval a) {


[16693]  942  if (Math.Abs(a.UpperBound.Value.Value  a.LowerBound.Value.Value) >= Math.PI * 2) {


 943  low = new MultivariateDual<AlgebraicDouble>(1.0);


 944  high = new MultivariateDual<AlgebraicDouble>(1.0);


 945  }


 946 


 947  //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)


 948  double Pihalf = Math.PI / 2;


[17200]  949  var scaled = a.Clone().Scale(1.0 / Pihalf);


[16693]  950  //move to positive scale


 951  if (scaled.LowerBound.Value.Value < 0) {


 952  int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1;


 953  scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4));


 954  }


 955 


 956  double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0;


 957  double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0;


 958  if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0;


 959  List<double> sinValues = new List<double>();


 960  sinValues.Add(Math.Sin(scaledLowerBound * Pihalf));


 961  sinValues.Add(Math.Sin(scaledUpperBound * Pihalf));


 962 


 963  int startValue = (int)Math.Ceiling(scaledLowerBound);


 964  while (startValue < scaledUpperBound) {


 965  sinValues.Add(Math.Sin(startValue * Pihalf));


 966  startValue += 1;


 967  }


 968 


[17200]  969  low = new MultivariateDual<AlgebraicDouble>(sinValues.Min()); // TODO, XXX FIX CALCULATION TO SUPPORT GRADIENTS!


[16693]  970  high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());


 971  return this;


[16674]  972  }


 973 


 974  public AlgebraicInterval Sub(AlgebraicInterval a) {


[16682]  975  // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]


 976  low.Sub(a.high);


 977  high.Sub(a.low);


[16674]  978  return this;


 979  }


 980 


 981  public AlgebraicInterval Clone() {


 982  return new AlgebraicInterval(low, high);


 983  }


 984 


 985  public bool Contains(double val) {


[16682]  986  return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;


[16674]  987  }


[16682]  988 


 989  public AlgebraicInterval AssignAbs(AlgebraicInterval a) {


 990  if (a.Contains(0.0)) {


 991  var abslow = a.low.Clone().Abs();


 992  var abshigh = a.high.Clone().Abs();


 993  a.high.Assign(Algebraic.Max(abslow, abshigh));


[16693]  994  a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound


[16682]  995  } else {


 996  var abslow = a.low.Clone().Abs();


 997  var abshigh = a.high.Clone().Abs();


 998  a.low.Assign(Algebraic.Min(abslow, abshigh));


 999  a.high.Assign(Algebraic.Max(abslow, abshigh));


 1000  }


 1001  return this;


 1002  }


 1003 


 1004  public AlgebraicInterval AssignSgn(AlgebraicInterval a) {


[16693]  1005  low.AssignSgn(a.low);


 1006  high.AssignSgn(a.high);


 1007  return this;


[16682]  1008  }


[16674]  1009  }


 1010 


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


 1012  where V : IAlgebraicType<V> {


[16695]  1013  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16674]  1014  private V v;


[16694]  1015  public V Value => v;


 1016 


[16695]  1017  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16674]  1018  private V dv;


[16682]  1019  public V Derivative => dv;


 1020 


[16694]  1021  public Dual(V v, V dv) { this.v = v; this.dv = dv; }


[16693]  1022 


[16695]  1023  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16693]  1024  public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);


[16727]  1025  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 1026  public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);


[16693]  1027 


[16694]  1028  public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }


 1029  public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }


 1030  public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }


 1031  public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }


 1032  public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }


 1033  public Dual<V> AssignInv(Dual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Mul(v).Mul(v); return this; } // (1/f(x))' =  f(x)' / f(x)^2


[16674]  1034 


[16694]  1035  // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);


[16674]  1036  public Dual<V> Mul(Dual<V> a) {


[16694]  1037  var t1 = a.dv.Clone().Mul(v);


 1038  var t2 = dv.Clone().Mul(a.v);


[16674]  1039  dv.Assign(t1).Add(t2);


 1040 


 1041  v.Mul(a.v);


 1042  return this;


 1043  }


[16694]  1044  public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }


[16674]  1045 


[16694]  1046  public Dual<V> AssignExp(Dual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Mul(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'


 1047  public Dual<V> AssignLog(Dual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Div(a.v); return this; } // log(x)' = 1/f(x) * f(x)'


[16674]  1048 


[16694]  1049  public Dual<V> AssignIntPower(Dual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Mul(a.v.Clone().IntPower(p  1)); return this; }


 1050  public Dual<V> AssignIntRoot(Dual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r  1)); return this; }


[16674]  1051 


[16694]  1052  public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }


 1053  public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }


[17131]  1054  public Dual<V> AssignTanh(Dual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Mul(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }


[16674]  1055 


[16694]  1056  public Dual<V> AssignAbs(Dual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Sgn()); return this; } // abs(f(x))' = f(x)*f'(x) / f(x)


 1057  public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }


[16682]  1058 


[16694]  1059  public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }


[16682]  1060 


[16674]  1061  }


 1062 


[16682]  1063  /// <summary>


 1064  /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.


 1065  /// </summary>


 1066  /// <typeparam name="V"></typeparam>


[16695]  1067  [DebuggerDisplay("v={Value}; dv={dv}")]


[16682]  1068  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {


[16695]  1069  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


[16674]  1070  private V v;


 1071  public V Value => v;


 1072 


[16695]  1073  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


 1074  private AlgebraicSparseVector<object, V> dv;


 1075  public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key


[16674]  1076 


[16694]  1077  private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }


[16674]  1078 


[16682]  1079  /// <summary>


 1080  /// Constructor without partial derivative


 1081  /// </summary>


 1082  /// <param name="v"></param>


[16695]  1083  public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }


[16674]  1084 


[16682]  1085  /// <summary>


 1086  /// Constructor for multiple partial derivatives


 1087  /// </summary>


 1088  /// <param name="v"></param>


 1089  /// <param name="keys"></param>


 1090  /// <param name="dv"></param>


[16695]  1091  public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }


[16674]  1092 


[16682]  1093  /// <summary>


 1094  /// Constructor for a single partial derivative


 1095  /// </summary>


 1096  /// <param name="v"></param>


 1097  /// <param name="key"></param>


 1098  /// <param name="dv"></param>


[16695]  1099  public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }


[16682]  1100 


[16693]  1101  /// <summary>


[16694]  1102  /// Constructor with a given value and gradient. For internal use.


[16693]  1103  /// </summary>


[16694]  1104  /// <param name="v">The value (not cloned).</param>


 1105  /// <param name="gradient">The gradient (not cloned).</param>


[16695]  1106  internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }


[16693]  1107 


[16694]  1108  public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }


[16682]  1109 


[16693]  1110  public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);


[16727]  1111  public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);


[16693]  1112 


[16694]  1113  public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }


[16693]  1114 


[16694]  1115  public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }


 1116  public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }


 1117  public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }


[16682]  1118  public MultivariateDual<V> Mul(MultivariateDual<V> a) {


[16674]  1119  // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);


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


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


[16682]  1122  dv.Assign(t1).Add(t2);


[16674]  1123 


 1124  v.Mul(a.v);


 1125  return this;


 1126  }


[16694]  1127  public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }


 1128  public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }


 1129  public MultivariateDual<V> AssignInv(MultivariateDual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Scale(v).Scale(v); return this; } // (1/f(x))' =  f(x)' / f(x)^2


[16674]  1130 


[16694]  1131  public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }


 1132  public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }


[17131]  1133  public MultivariateDual<V> AssignTanh(MultivariateDual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Scale(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; } // tanh(f(x))' = f(x)'sech²(f(x)) = f(x)'(1  tanh²(f(x)))


[16674]  1134 


[16694]  1135  public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Scale(a.v.Clone().IntPower(p  1)); return this; }


 1136  public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Scale(a.v.IntRoot(r  1)); return this; }


[16682]  1137 


[16694]  1138  public MultivariateDual<V> AssignExp(MultivariateDual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Scale(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'


 1139  public MultivariateDual<V> AssignLog(MultivariateDual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Inv()); return this; } // log(x)' = 1/f(x) * f(x)'


[16682]  1140 


[16694]  1141  public MultivariateDual<V> AssignAbs(MultivariateDual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Sgn()); return this; } // abs(f(x))' = f(x)*f'(x) / f(x) doesn't work for intervals


 1142  public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;


[16674]  1143  }


[16682]  1144  } 
