Changeset 17295
 Timestamp:
 10/02/19 17:09:31 (12 days ago)
 Location:
 branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter
 Files:

 11 added
 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
r17294 r17295 1 1 using System; 2 using System.Collections.Generic;3 using System.Diagnostics;4 using System.Linq;5 2 using HeuristicLab.Common; 6 3 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding; … … 156 153 157 154 } 158 159 160 public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {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;209 instruction.value = new AlgebraicDoubleVector(BATCHSIZE);210 instruction.value.AssignConstant(instruction.dblVal);211 }212 213 protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {214 instruction.dblVal = variable.Weight;215 instruction.value = new AlgebraicDoubleVector(BATCHSIZE);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) {225 instruction.value = new AlgebraicDoubleVector(BATCHSIZE);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 235 public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {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 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) {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 JAC290 var g = code[0].value.Gradient;291 for (int j = 0; j < nParams; ++j) {292 if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {293 v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);294 } else {295 for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;296 }297 }298 }299 300 if (remainingRows > 0) {301 Evaluate(code);302 code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);303 304 var g = code[0].value.Gradient;305 for (int j = 0; j < nParams; ++j)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 }311 }312 }313 314 protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {315 var zero = new AlgebraicDoubleVector(BATCHSIZE);316 instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);317 }318 319 protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {320 var g_arr = new double[BATCHSIZE];321 if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {322 for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;323 var g = new AlgebraicDoubleVector(g_arr);324 instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient325 } else {326 instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));327 }328 329 instruction.dblVal = constant.Value;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];345 var f = new AlgebraicDoubleVector(BATCHSIZE);346 var g = new AlgebraicDoubleVector(BATCHSIZE);347 instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);348 } else {349 var f = new AlgebraicDoubleVector(BATCHSIZE);350 instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);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 values366 var g = a.value.Gradient.Elements[paramIdx];367 for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) {368 g[i  rowIndex] = data[rows[i]];369 }370 }371 }372 }373 374 375 public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {376 [ThreadStatic]377 private IDictionary<string, Interval> intervals;378 379 public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals) {380 this.intervals = intervals;381 var code = Compile(tree);382 Evaluate(code);383 if (code[0].value.LowerBound.Value.Value > code[0].value.UpperBound.Value.Value) throw new InvalidProgramException($"lower: {code[0].value.LowerBound.Value.Value} > upper: {code[0].value.UpperBound.Value.Value}");384 return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);385 }386 387 public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {388 this.intervals = intervals;389 var code = Compile(tree);390 Evaluate(code);391 lowerGradient = new double[paramNodes.Length];392 upperGradient = new double[paramNodes.Length];393 var l = code[0].value.LowerBound;394 var u = code[0].value.UpperBound;395 for (int i = 0; i < paramNodes.Length; ++i) {396 if (paramNodes[i] == null) continue;397 if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;398 if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;399 }400 return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);401 }402 403 protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {404 instruction.value = new AlgebraicInterval(0, 0);405 }406 407 408 protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {409 instruction.dblVal = constant.Value;410 instruction.value = new AlgebraicInterval(411 new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0),412 new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0) // use node as key413 );414 }415 416 protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {417 instruction.dblVal = variable.Weight;418 var v1 = instruction.dblVal * intervals[variable.VariableName].LowerBound;419 var v2 = instruction.dblVal * intervals[variable.VariableName].UpperBound;420 var lower = Math.Min(v1, v2);421 var upper = Math.Max(v1, v2);422 // we assume that the for variable nodes ( v(x,w) = w * x ) the gradient is returned for parameter w423 instruction.value = new AlgebraicInterval(424 low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),425 high: new MultivariateDual<AlgebraicDouble>(v: upper, key: variable, dv: intervals[variable.VariableName].UpperBound)426 );427 }428 429 protected override void LoadVariable(Instruction a) {430 // nothing to do431 }432 }433 434 public interface IAlgebraicType<T> {435 T Zero { get; } // Zero and One must create new objects436 T One { get; }437 438 T AssignAbs(T a); // set this to assign abs(a)439 T Assign(T a); // assign this to same value as a (copy!)440 T AssignNeg(T a); // set this to negative(a)441 T AssignInv(T a); // set this to inv(a);442 T Scale(double s); // scale this with s443 T Add(T a); // add a to this444 T Sub(T a); // subtract a from this445 T Mul(T a); // multiply this with a446 T Div(T a); // divide this by a447 T AssignLog(T a); // set this to log a448 T AssignExp(T a); // set this to exp(a)449 T AssignSin(T a); // set this to sin(a)450 T AssignCos(T a); // set this to cos(a)451 T AssignTanh(T a); // set this to tanh(a)452 T AssignIntPower(T a, int p);453 T AssignIntRoot(T a, int r);454 T AssignSgn(T a); // set this to sign(a)455 T AssignMin(T other); // set this min(this, other)456 T AssignMax(T other); // set this max(this, other)457 T Clone();458 }459 460 public static class Algebraic {461 public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }462 public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }463 public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }464 public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }465 public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }466 public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }467 public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }468 public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }469 public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }470 public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }471 472 internal static T Min<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMin(b); }473 internal static T Max<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMax(b); }474 475 // public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {476 // // ((a + b) + abs(b  a)) / 2477 // return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);478 // }479 // public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {480 // // ((a + b)  abs(a  b)) / 2481 // return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);482 // }483 }484 485 486 // algebraic type wrapper for a double value487 [DebuggerDisplay("{Value}")]488 public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {489 public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }490 public static implicit operator double(AlgebraicDouble value) { return value.Value; }491 public double Value;492 493 [DebuggerBrowsable(DebuggerBrowsableState.Never)]494 public AlgebraicDouble Zero => new AlgebraicDouble(0.0);495 [DebuggerBrowsable(DebuggerBrowsableState.Never)]496 public AlgebraicDouble One => new AlgebraicDouble(1.0);497 498 public bool IsInfinity => IsNegativeInfinity  IsPositiveInfinity;499 public bool IsNegativeInfinity => double.IsNegativeInfinity(Value);500 public bool IsPositiveInfinity => double.IsPositiveInfinity(Value);501 public AlgebraicDouble() { }502 public AlgebraicDouble(double value) { this.Value = value; }503 public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }504 public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }505 public AlgebraicDouble Sub(AlgebraicDouble a) { Value = a.Value; return this; }506 public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }507 public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }508 public AlgebraicDouble Scale(double s) { Value *= s; return this; }509 public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }510 public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = a.Value; return this; }511 public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }512 public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }513 public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }514 public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }515 public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }516 public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }517 public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = IntRoot(a.Value, r); return this; }518 public AlgebraicDouble AssignMin(AlgebraicDouble other) { Value = Math.Min(Value, other.Value); return this; }519 public AlgebraicDouble AssignMax(AlgebraicDouble other) { Value = Math.Max(Value, other.Value); return this; }520 521 // helper522 private static double IntRoot(double value, int r) {523 if (r % 2 == 0) return Math.Pow(value, 1.0 / r);524 else return value < 0 ? Math.Pow(value, 1.0 / r) : Math.Pow(value, 1.0 / r);525 }526 527 public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }528 public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }529 public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }530 531 public override string ToString() {532 return Value.ToString();533 }534 }535 536 // a simple vector as an algebraic type537 [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]538 public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {539 private double[] arr;540 public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }541 public int Length => arr.Length;542 543 public AlgebraicDoubleVector(int length) { arr = new double[length]; }544 545 public AlgebraicDoubleVector() { }546 547 /// <summary>548 ///549 /// </summary>550 /// <param name="arr">array is not copied</param>551 public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }552 553 [DebuggerBrowsable(DebuggerBrowsableState.Never)]554 public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)555 [DebuggerBrowsable(DebuggerBrowsableState.Never)]556 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)557 public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }558 public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }559 public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }560 public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }561 public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }562 public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }563 public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }564 public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }565 public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }566 public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }567 public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }568 public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }569 public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }570 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; }571 public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = IntRoot(a.arr[i], r); } return this; }572 public AlgebraicDoubleVector AssignMin(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Min(arr[i], other.arr[i]); } return this; }573 public AlgebraicDoubleVector AssignMax(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Max(arr[i], other.arr[i]); } return this; }574 575 // helper576 private double IntRoot(double v, int r) {577 if (r % 2 == 0) return Math.Pow(v, 1.0 / r);578 else return v < 0 ? Math.Pow(v, 1.0 / r) : Math.Pow(v, 1.0 / r);579 }580 581 public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }582 public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }583 584 public AlgebraicDoubleVector Clone() {585 var v = new AlgebraicDoubleVector(this.arr.Length);586 Array.Copy(arr, v.arr, v.arr.Length);587 return v;588 }589 590 public AlgebraicDoubleVector AssignConstant(double constantValue) {591 for (int i = 0; i < arr.Length; ++i) {592 arr[i] = constantValue;593 }594 return this;595 }596 597 public void CopyTo(double[] dest, int idx, int length) {598 Array.Copy(arr, 0, dest, idx, length);599 }600 601 public void CopyFrom(double[] data, int rowIndex) {602 Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length  rowIndex));603 }604 public void CopyRowTo(double[,] dest, int row) {605 for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];606 }607 608 internal void CopyColumnTo(double[,] dest, int column, int row, int len) {609 for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];610 }611 612 public override string ToString() {613 return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";614 }615 }616 617 618 /*619 // vectors of algebraic types620 public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {621 private T[] elems;622 623 public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }624 625 public int Length => elems.Length;626 627 private AlgebraicVector() { }628 629 public AlgebraicVector(int len) { elems = new T[len]; }630 631 /// <summary>632 ///633 /// </summary>634 /// <param name="elems">The array is copied (elementwise clone)</param>635 public AlgebraicVector(T[] elems) {636 this.elems = new T[elems.Length];637 for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }638 }639 640 /// <summary>641 ///642 /// </summary>643 /// <param name="elems">Array is not copied!</param>644 /// <returns></returns>645 public AlgebraicVector<T> FromArray(T[] elems) {646 var v = new AlgebraicVector<T>();647 v.elems = elems;648 return v;649 }650 651 public void CopyTo(T[] dest) {652 if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");653 Array.Copy(elems, dest, dest.Length);654 }655 656 public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }657 658 659 [DebuggerBrowsable(DebuggerBrowsableState.Never)]660 public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);661 [DebuggerBrowsable(DebuggerBrowsableState.Never)]662 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; } }663 public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }664 public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }665 public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }666 public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }667 public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }668 public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }669 public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }670 public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }671 public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }672 public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }673 public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }674 public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }675 public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }676 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; }677 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; }678 public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }679 public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }680 }681 682 */683 684 685 /// <summary>686 /// A sparse vector of algebraic types. Elements are accessed via a key of type K687 /// </summary>688 /// <typeparam name="K">Key type</typeparam>689 /// <typeparam name="T">Element type</typeparam>690 [DebuggerDisplay("SparseVector: {ToString()}")]691 public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {692 [DebuggerBrowsable(DebuggerBrowsableState.Never)]693 private Dictionary<K, T> elems;694 public IReadOnlyDictionary<K, T> Elements => elems;695 696 697 public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {698 elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());699 }700 701 /// <summary>702 ///703 /// </summary>704 /// <param name="keys"></param>705 /// <param name="values">values are cloned</param>706 public AlgebraicSparseVector(K[] keys, T[] values) {707 if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");708 elems = new Dictionary<K, T>(keys.Length);709 for (int i = 0; i < keys.Length; ++i) {710 elems.Add(keys[i], values[i].Clone());711 }712 }713 714 public AlgebraicSparseVector() {715 this.elems = new Dictionary<K, T>();716 }717 718 // keep only elements from a719 private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {720 // remove elems from this which don't occur in a721 List<K> keysToRemove = new List<K>();722 foreach (var kvp in elems) {723 if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);724 }725 foreach (var o in keysToRemove) elems.Remove(o); // > zero726 727 foreach (var kvp in a.elems) {728 if (elems.TryGetValue(kvp.Key, out T value))729 mapAssign(kvp.Value, value);730 else731 elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));732 }733 }734 735 [DebuggerBrowsable(DebuggerBrowsableState.Never)]736 public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();737 [DebuggerBrowsable(DebuggerBrowsableState.Never)]738 public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();739 740 public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }741 public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }742 743 public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }744 public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }745 public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }746 public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }747 public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }748 public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }749 public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }750 public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }751 public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }752 public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }753 public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }754 public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }755 public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {756 foreach (var kvp in a.elems) {757 if (elems.TryGetValue(kvp.Key, out T value))758 value.Add(kvp.Value);759 else760 elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a761 }762 return this;763 }764 765 public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {766 foreach (var kvp in a.elems) {767 if (elems.TryGetValue(kvp.Key, out T value))768 value.Sub(kvp.Value);769 else770 elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0  a771 }772 return this;773 }774 775 public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {776 var keys = elems.Keys.ToArray();777 foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0778 foreach (var kvp in a.elems) {779 if (elems.TryGetValue(kvp.Key, out T value))780 value.Mul(kvp.Value); // this * a781 }782 return this;783 }784 785 public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {786 return Mul(a.Clone().Inv());787 }788 789 public AlgebraicSparseVector<K, T> AssignMin(AlgebraicSparseVector<K, T> other) {790 // assumes that keys without a matching key in other are zero and vice versa791 foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMin(kvp.Value.Zero); // min(v, 0)792 foreach (var kvp in other.elems) {793 if (elems.TryGetValue(kvp.Key, out T value))794 value.AssignMin(kvp.Value);795 else796 elems.Add(kvp.Key, kvp.Value.Zero.AssignMin(kvp.Value));797 }798 return this;799 }800 801 public AlgebraicSparseVector<K, T> AssignMax(AlgebraicSparseVector<K, T> other) {802 // assumes that keys without a matching key in other are zero and vice versa803 foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMax(kvp.Value.Zero); // max(v, 0)804 foreach (var kvp in other.elems) {805 if (elems.TryGetValue(kvp.Key, out T value))806 value.AssignMax(kvp.Value);807 else808 elems.Add(kvp.Key, kvp.Value.Zero.AssignMax(kvp.Value));809 }810 return this;811 }812 813 814 public AlgebraicSparseVector<K, T> Clone() {815 return new AlgebraicSparseVector<K, T>(this);816 }817 818 public override string ToString() {819 return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";820 }821 }822 823 // this is our own implementation of interval arithmetic824 // for a well worked out definition of interval operations for IEEE reals see:825 // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995826 [DebuggerDisplay("[{low.Value}..{high.Value}]")]827 public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {828 [DebuggerBrowsable(DebuggerBrowsableState.Never)]829 private MultivariateDual<AlgebraicDouble> low;830 public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();831 832 [DebuggerBrowsable(DebuggerBrowsableState.Never)]833 private MultivariateDual<AlgebraicDouble> high;834 public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();835 836 837 public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }838 839 public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {840 this.low = low.Clone();841 this.high = high.Clone();842 }843 844 public AlgebraicInterval(double low, double high) {845 this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));846 this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));847 }848 849 [DebuggerBrowsable(DebuggerBrowsableState.Never)]850 public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);851 [DebuggerBrowsable(DebuggerBrowsableState.Never)]852 public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);853 854 public AlgebraicInterval Add(AlgebraicInterval a) {855 low.Add(a.low);856 high.Add(a.high);857 return this;858 }859 860 public AlgebraicInterval Mul(AlgebraicInterval a) {861 var v1 = low.Clone().Mul(a.low);862 var v2 = low.Clone().Mul(a.high);863 var v3 = high.Clone().Mul(a.low);864 var v4 = high.Clone().Mul(a.high);865 866 low = Min(Min(v1, v2), Min(v3, v4));867 high = Max(Max(v1, v2), Max(v3, v4));868 869 return this;870 }871 872 873 private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {874 return a.Value < b.Value ? a : b;875 }876 private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {877 return a.Value > b.Value ? a : b;878 }879 880 public AlgebraicInterval Assign(AlgebraicInterval a) {881 low = a.low;882 high = a.high;883 return this;884 }885 886 public AlgebraicInterval AssignCos(AlgebraicInterval a) {887 return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));888 }889 890 public AlgebraicInterval Div(AlgebraicInterval a) {891 if (a.Contains(0.0)) {892 if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {893 if (this.low.Value >= 0) {894 // pos / 0895 } else if (this.high.Value <= 0) {896 // neg / 0897 } else {898 low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);899 high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);900 }901 } else if (a.low.Value.Value >= 0) {902 // a is positive903 Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));904 } else if (a.high.Value <= 0) {905 // a is negative906 Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));907 } else {908 // a is interval over zero909 low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);910 high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);911 }912 } else {913 Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low914 }915 return this;916 }917 918 public AlgebraicInterval AssignExp(AlgebraicInterval a) {919 low.AssignExp(a.low);920 high.AssignExp(a.high);921 return this;922 }923 924 // tanh is a bijective function925 public AlgebraicInterval AssignTanh(AlgebraicInterval a) {926 low.AssignTanh(a.low);927 high.AssignTanh(a.high);928 return this;929 }930 931 public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {932 if (p < 0) { // x^3 == 1/(x^3)933 AssignIntPower(a, p);934 return AssignInv(this);935 } else if (p == 0) {936 if (a.Contains(0.0)) {937 // => 0^0 = 0 ; might not be relevant938 low = new MultivariateDual<AlgebraicDouble>(0.0);939 high = new MultivariateDual<AlgebraicDouble>(1.0);940 return this;941 } else {942 // => 1943 low = new MultivariateDual<AlgebraicDouble>(1.0);944 high = new MultivariateDual<AlgebraicDouble>(1.0);945 return this;946 }947 } else if (p == 1) return this;948 else if (p % 2 == 0) {949 // p is even => interval must be positive950 if (a.Contains(0.0)) {951 low = new MultivariateDual<AlgebraicDouble>(0.0);952 high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));953 } else {954 var lowPower = a.low.IntPower(p);955 var highPower = a.high.IntPower(p);956 low = lowPower.AssignMin(highPower);957 high = lowPower.AssignMax(highPower);958 }959 } else {960 // p is uneven961 if (a.Contains(0.0)) {962 low.AssignIntPower(a.low, p);963 high.AssignIntPower(a.high, p);964 } else {965 var lowPower = a.low.IntPower(p);966 var highPower = a.high.IntPower(p);967 low = lowPower.AssignMin(highPower);968 high = lowPower.AssignMax(highPower);969 }970 }971 return this;972 }973 974 public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {975 if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }976 if (r == 1) return this;977 if (r < 0) {978 // x^ (1/2) = 1 / (x^(1/2))979 AssignIntRoot(a, r);980 return AssignInv(this);981 } else {982 // root only defined for positive arguments for even roots983 if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {984 low = new MultivariateDual<AlgebraicDouble>(double.NaN);985 high = new MultivariateDual<AlgebraicDouble>(double.NaN);986 return this;987 } else {988 low.AssignIntRoot(a.low, r);989 high.AssignIntRoot(a.high, r);990 return this;991 }992 }993 }994 995 public AlgebraicInterval AssignInv(AlgebraicInterval a) {996 low = new MultivariateDual<AlgebraicDouble>(1.0);997 high = new MultivariateDual<AlgebraicDouble>(1.0);998 return Div(a);999 }1000 1001 public AlgebraicInterval AssignLog(AlgebraicInterval a) {1002 low.AssignLog(a.low);1003 high.AssignLog(a.high);1004 return this;1005 }1006 1007 public AlgebraicInterval AssignNeg(AlgebraicInterval a) {1008 low.AssignNeg(a.high);1009 high.AssignNeg(a.low);1010 return this;1011 }1012 1013 public AlgebraicInterval Scale(double s) {1014 low.Scale(s);1015 high.Scale(s);1016 if (s < 0) {1017 var t = low;1018 low = high;1019 high = t;1020 }1021 return this;1022 }1023 1024 public AlgebraicInterval AssignSin(AlgebraicInterval a) {1025 var lower = a.LowerBound.Value.Value;1026 var size = a.UpperBound.Value.Value  lower;1027 if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;1028 1029 if (size >= Math.PI * 2) {1030 low = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient1031 high = new MultivariateDual<AlgebraicDouble>(1.0);1032 return this;1033 }1034 1035 // assume low and high are in the same quadrant1036 low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());1037 high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());1038 1039 // override min and max if necessary1040 1041 // shift interval 'a' into the range [2pi .. 2pi] without changing the size of the interval to simplify the checks1042 lower = lower % (2 * Math.PI); // lower in [2pi .. 2pi]1043 1044 // handle min = 1 and max = 1 cases explicitly1045 var pi_2 = Math.PI / 2.0;1046 var maxima = new double[] { 3 * pi_2, pi_2 };1047 var minima = new double[] { pi_2, 3 * pi_2 };1048 1049 // override min and max if necessary1050 if (maxima.Any(m => lower < m && lower + size > m)) {1051 // max = 11052 high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient1053 }1054 1055 if (minima.Any(m => lower < m && lower + size > m)) {1056 // min = 1;1057 low = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient1058 }1059 return this;1060 }1061 1062 public AlgebraicInterval Sub(AlgebraicInterval a) {1063 // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]1064 low.Sub(a.high);1065 high.Sub(a.low);1066 return this;1067 }1068 1069 public AlgebraicInterval Clone() {1070 return new AlgebraicInterval(low, high);1071 }1072 1073 public bool Contains(double val) {1074 return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;1075 }1076 1077 public AlgebraicInterval AssignAbs(AlgebraicInterval a) {1078 if (a.Contains(0.0)) {1079 var abslow = a.low.Clone().Abs();1080 var abshigh = a.high.Clone().Abs();1081 a.high.Assign(Algebraic.Max(abslow, abshigh));1082 a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound1083 } else {1084 var abslow = a.low.Clone().Abs();1085 var abshigh = a.high.Clone().Abs();1086 a.low.Assign(Algebraic.Min(abslow, abshigh));1087 a.high.Assign(Algebraic.Max(abslow, abshigh));1088 }1089 return this;1090 }1091 1092 public AlgebraicInterval AssignSgn(AlgebraicInterval a) {1093 low.AssignSgn(a.low);1094 high.AssignSgn(a.high);1095 return this;1096 }1097 1098 public AlgebraicInterval AssignMin(AlgebraicInterval other) {1099 low.AssignMin(other.low);1100 high.AssignMin(other.high);1101 return this;1102 }1103 1104 public AlgebraicInterval AssignMax(AlgebraicInterval other) {1105 low.AssignMax(other.low);1106 high.AssignMax(other.high);1107 return this;1108 }1109 }1110 1111 public class Dual<V> : IAlgebraicType<Dual<V>>1112 where V : IAlgebraicType<V> {1113 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1114 private V v;1115 public V Value => v;1116 1117 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1118 private V dv;1119 public V Derivative => dv;1120 1121 public Dual(V v, V dv) { this.v = v; this.dv = dv; }1122 1123 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1124 public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);1125 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1126 public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);1127 1128 public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }1129 public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }1130 public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }1131 public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }1132 public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }1133 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)^21134 1135 // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);1136 public Dual<V> Mul(Dual<V> a) {1137 var t1 = a.dv.Clone().Mul(v);1138 var t2 = dv.Clone().Mul(a.v);1139 dv.Assign(t1).Add(t2);1140 1141 v.Mul(a.v);1142 return this;1143 }1144 public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }1145 1146 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)'1147 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)'1148 1149 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; }1150 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; }1151 1152 public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }1153 public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }1154 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; }1155 1156 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)1157 public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }1158 1159 public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }1160 1161 public Dual<V> AssignMin(Dual<V> other) {1162 throw new NotImplementedException();1163 }1164 1165 public Dual<V> AssignMax(Dual<V> other) {1166 throw new NotImplementedException();1167 }1168 }1169 1170 /// <summary>1171 /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.1172 /// </summary>1173 /// <typeparam name="V"></typeparam>1174 [DebuggerDisplay("v={Value}; dv={dv}")]1175 public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {1176 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1177 private V v;1178 public V Value => v;1179 1180 [DebuggerBrowsable(DebuggerBrowsableState.Never)]1181 private AlgebraicSparseVector<object, V> dv;1182 public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key1183 1184 private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }1185 1186 /// <summary>1187 /// Constructor without partial derivative1188 /// </summary>1189 /// <param name="v"></param>1190 public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }1191 1192 /// <summary>1193 /// Constructor for multiple partial derivatives1194 /// </summary>1195 /// <param name="v"></param>1196 /// <param name="keys"></param>1197 /// <param name="dv"></param>1198 public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }1199 1200 /// <summary>1201 /// Constructor for a single partial derivative1202 /// </summary>1203 /// <param name="v"></param>1204 /// <param name="key"></param>1205 /// <param name="dv"></param>1206 public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }1207 1208 /// <summary>1209 /// Constructor with a given value and gradient. For internal use.1210 /// </summary>1211 /// <param name="v">The value (not cloned).</param>1212 /// <param name="gradient">The gradient (not cloned).</param>1213 internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }1214 1215 public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }1216 1217 public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);1218 public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);1219 1220 public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }1221 1222 public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }1223 public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }1224 public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }1225 public MultivariateDual<V> Mul(MultivariateDual<V> a) {1226 // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);1227 var t1 = a.dv.Clone().Scale(v);1228 var t2 = dv.Clone().Scale(a.v);1229 dv.Assign(t1).Add(t2);1230 1231 v.Mul(a.v);1232 return this;1233 }1234 public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }1235 public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }1236 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)^21237 1238 public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }1239 public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }1240 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)))1241 1242 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; }1243 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; }1244 1245 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)'1246 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)'1247 1248 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 intervals1249 public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;1250 1251 public MultivariateDual<V> AssignMin(MultivariateDual<V> other) {1252 XXX1253 }1254 1255 public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {1256 XXX1257 }1258 }1259 155 }
Note: See TracChangeset
for help on using the changeset viewer.