[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 | }
|
---|
[16694] | 115 | default: throw new ArgumentException($"Unknown opcode {code[c].opcode}");
|
---|
[16674] | 116 | }
|
---|
| 117 | }
|
---|
| 118 | return code[0].value;
|
---|
| 119 | }
|
---|
| 120 |
|
---|
| 121 | protected Instruction[] Compile(ISymbolicExpressionTree tree) {
|
---|
| 122 | var root = tree.Root.GetSubtree(0).GetSubtree(0);
|
---|
| 123 | var code = new Instruction[root.GetLength()];
|
---|
| 124 | if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
|
---|
| 125 | int c = 1, i = 0;
|
---|
| 126 | foreach (var node in root.IterateNodesBreadth()) {
|
---|
| 127 | if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
|
---|
| 128 | code[i] = new Instruction {
|
---|
| 129 | opcode = OpCodes.MapSymbolToOpCode(node),
|
---|
| 130 | narg = (ushort)node.SubtreeCount,
|
---|
| 131 | childIndex = c
|
---|
| 132 | };
|
---|
| 133 | if (node is VariableTreeNode variable) {
|
---|
| 134 | InitializeTerminalInstruction(ref code[i], variable);
|
---|
| 135 | } else if (node is ConstantTreeNode constant) {
|
---|
| 136 | InitializeTerminalInstruction(ref code[i], constant);
|
---|
| 137 | } else {
|
---|
| 138 | InitializeInternalInstruction(ref code[i], node);
|
---|
| 139 | }
|
---|
| 140 | c += node.SubtreeCount;
|
---|
| 141 | ++i;
|
---|
| 142 | }
|
---|
| 143 | return code;
|
---|
| 144 | }
|
---|
| 145 |
|
---|
| 146 | protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant);
|
---|
| 147 | protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable);
|
---|
| 148 | protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node);
|
---|
| 149 |
|
---|
| 150 | protected abstract void LoadVariable(Instruction a);
|
---|
| 151 |
|
---|
| 152 | }
|
---|
| 153 |
|
---|
| 154 |
|
---|
[16695] | 155 | public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {
|
---|
[16674] | 156 | private const int BATCHSIZE = 128;
|
---|
| 157 | [ThreadStatic]
|
---|
| 158 | private Dictionary<string, double[]> cachedData;
|
---|
| 159 |
|
---|
| 160 | [ThreadStatic]
|
---|
| 161 | private IDataset dataset;
|
---|
| 162 |
|
---|
| 163 | [ThreadStatic]
|
---|
| 164 | private int rowIndex;
|
---|
| 165 |
|
---|
| 166 | [ThreadStatic]
|
---|
| 167 | private int[] rows;
|
---|
| 168 |
|
---|
| 169 | private void InitCache(IDataset dataset) {
|
---|
| 170 | this.dataset = dataset;
|
---|
| 171 | cachedData = new Dictionary<string, double[]>();
|
---|
| 172 | foreach (var v in dataset.DoubleVariables) {
|
---|
| 173 | cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
|
---|
| 174 | }
|
---|
| 175 | }
|
---|
| 176 |
|
---|
| 177 | public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
|
---|
| 178 | if (cachedData == null || this.dataset != dataset) {
|
---|
| 179 | InitCache(dataset);
|
---|
| 180 | }
|
---|
| 181 |
|
---|
| 182 | this.rows = rows;
|
---|
| 183 | var code = Compile(tree);
|
---|
| 184 | var remainingRows = rows.Length % BATCHSIZE;
|
---|
| 185 | var roundedTotal = rows.Length - remainingRows;
|
---|
| 186 |
|
---|
| 187 | var result = new double[rows.Length];
|
---|
| 188 |
|
---|
| 189 | for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
|
---|
| 190 | Evaluate(code);
|
---|
| 191 | code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | if (remainingRows > 0) {
|
---|
| 195 | Evaluate(code);
|
---|
| 196 | code[0].value.CopyTo(result, roundedTotal, remainingRows);
|
---|
| 197 | }
|
---|
| 198 |
|
---|
| 199 | return result;
|
---|
| 200 | }
|
---|
| 201 |
|
---|
| 202 | protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
|
---|
| 203 | instruction.dblVal = constant.Value;
|
---|
[16695] | 204 | instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
[16674] | 205 | instruction.value.AssignConstant(instruction.dblVal);
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
|
---|
| 209 | instruction.dblVal = variable.Weight;
|
---|
[16695] | 210 | instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
[16674] | 211 | if (cachedData.ContainsKey(variable.VariableName)) {
|
---|
| 212 | instruction.data = cachedData[variable.VariableName];
|
---|
| 213 | } else {
|
---|
| 214 | instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
|
---|
| 215 | cachedData[variable.VariableName] = (double[])instruction.data;
|
---|
| 216 | }
|
---|
| 217 | }
|
---|
| 218 |
|
---|
| 219 | protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
|
---|
[16695] | 220 | instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
[16674] | 221 | }
|
---|
| 222 |
|
---|
| 223 | protected override void LoadVariable(Instruction a) {
|
---|
| 224 | var data = (double[])a.data;
|
---|
| 225 | for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
|
---|
| 226 | a.value.Scale(a.dblVal);
|
---|
| 227 | }
|
---|
| 228 | }
|
---|
| 229 |
|
---|
[16695] | 230 | public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {
|
---|
[16674] | 231 | private const int BATCHSIZE = 128;
|
---|
| 232 | [ThreadStatic]
|
---|
| 233 | private Dictionary<string, double[]> cachedData;
|
---|
| 234 |
|
---|
| 235 | [ThreadStatic]
|
---|
| 236 | private IDataset dataset;
|
---|
| 237 |
|
---|
| 238 | [ThreadStatic]
|
---|
| 239 | private int rowIndex;
|
---|
| 240 |
|
---|
| 241 | [ThreadStatic]
|
---|
| 242 | private int[] rows;
|
---|
| 243 |
|
---|
| 244 | [ThreadStatic]
|
---|
| 245 | private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
|
---|
| 246 |
|
---|
| 247 | private void InitCache(IDataset dataset) {
|
---|
| 248 | this.dataset = dataset;
|
---|
| 249 | cachedData = new Dictionary<string, double[]>();
|
---|
| 250 | foreach (var v in dataset.DoubleVariables) {
|
---|
| 251 | cachedData[v] = dataset.GetDoubleValues(v).ToArray();
|
---|
| 252 | }
|
---|
| 253 | }
|
---|
| 254 |
|
---|
[16727] | 255 | /// <summary>
|
---|
| 256 | ///
|
---|
| 257 | /// </summary>
|
---|
| 258 | /// <param name="tree"></param>
|
---|
| 259 | /// <param name="dataset"></param>
|
---|
| 260 | /// <param name="rows"></param>
|
---|
| 261 | /// <param name="parameterNodes"></param>
|
---|
| 262 | /// <param name="fi">Function output. Must be allocated by the caller.</param>
|
---|
| 263 | /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
|
---|
| 264 | public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
|
---|
[16674] | 265 | if (cachedData == null || this.dataset != dataset) {
|
---|
| 266 | InitCache(dataset);
|
---|
| 267 | }
|
---|
| 268 |
|
---|
| 269 | int nParams = parameterNodes.Length;
|
---|
| 270 | node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
|
---|
| 271 | for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
|
---|
| 272 |
|
---|
| 273 | var code = Compile(tree);
|
---|
| 274 |
|
---|
| 275 | var remainingRows = rows.Length % BATCHSIZE;
|
---|
| 276 | var roundedTotal = rows.Length - remainingRows;
|
---|
| 277 |
|
---|
| 278 | this.rows = rows;
|
---|
| 279 |
|
---|
| 280 | for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
|
---|
| 281 | Evaluate(code);
|
---|
| 282 | code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
|
---|
| 283 |
|
---|
| 284 | // TRANSPOSE into JAC
|
---|
[16682] | 285 | var g = code[0].value.Gradient;
|
---|
| 286 | for (int j = 0; j < nParams; ++j) {
|
---|
[16744] | 287 | if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
|
---|
[16727] | 288 | v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
|
---|
| 289 | } else {
|
---|
| 290 | for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
|
---|
| 291 | }
|
---|
[16682] | 292 | }
|
---|
[16674] | 293 | }
|
---|
| 294 |
|
---|
| 295 | if (remainingRows > 0) {
|
---|
| 296 | Evaluate(code);
|
---|
| 297 | code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
|
---|
[16682] | 298 |
|
---|
| 299 | var g = code[0].value.Gradient;
|
---|
[16674] | 300 | for (int j = 0; j < nParams; ++j)
|
---|
[16727] | 301 | if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
|
---|
| 302 | v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
|
---|
| 303 | } else {
|
---|
| 304 | for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
|
---|
| 305 | }
|
---|
[16674] | 306 | }
|
---|
| 307 | }
|
---|
| 308 |
|
---|
| 309 | protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
|
---|
[16695] | 310 | var zero = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
| 311 | instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);
|
---|
[16674] | 312 | }
|
---|
| 313 |
|
---|
| 314 | protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
|
---|
| 315 | var g_arr = new double[BATCHSIZE];
|
---|
[16682] | 316 | if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {
|
---|
[16674] | 317 | for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
|
---|
[16695] | 318 | var g = new AlgebraicDoubleVector(g_arr);
|
---|
| 319 | instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
|
---|
[16682] | 320 | } else {
|
---|
[16695] | 321 | instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
|
---|
[16674] | 322 | }
|
---|
| 323 |
|
---|
[16682] | 324 | instruction.dblVal = constant.Value;
|
---|
[16674] | 325 | instruction.value.Value.AssignConstant(instruction.dblVal);
|
---|
| 326 | }
|
---|
| 327 |
|
---|
| 328 | protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
|
---|
| 329 | double[] data;
|
---|
| 330 | if (cachedData.ContainsKey(variable.VariableName)) {
|
---|
| 331 | data = cachedData[variable.VariableName];
|
---|
| 332 | } else {
|
---|
| 333 | data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
|
---|
| 334 | cachedData[variable.VariableName] = (double[])instruction.data;
|
---|
| 335 | }
|
---|
| 336 |
|
---|
| 337 | var paramIdx = -1;
|
---|
| 338 | if (node2paramIdx.ContainsKey(variable)) {
|
---|
| 339 | paramIdx = node2paramIdx[variable];
|
---|
[16695] | 340 | var f = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
| 341 | var g = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
| 342 | instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);
|
---|
[16682] | 343 | } else {
|
---|
[16695] | 344 | var f = new AlgebraicDoubleVector(BATCHSIZE);
|
---|
| 345 | instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
|
---|
[16674] | 346 | }
|
---|
| 347 |
|
---|
| 348 | instruction.dblVal = variable.Weight;
|
---|
| 349 | instruction.data = new object[] { data, paramIdx };
|
---|
| 350 | }
|
---|
| 351 |
|
---|
| 352 | protected override void LoadVariable(Instruction a) {
|
---|
| 353 | var paramIdx = (int)((object[])a.data)[1];
|
---|
| 354 | var data = (double[])((object[])a.data)[0];
|
---|
| 355 |
|
---|
| 356 | for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value.Value[i - rowIndex] = data[rows[i]];
|
---|
| 357 | a.value.Scale(a.dblVal);
|
---|
| 358 |
|
---|
| 359 | if (paramIdx >= 0) {
|
---|
| 360 | // update gradient with variable values
|
---|
[16682] | 361 | var g = a.value.Gradient.Elements[paramIdx];
|
---|
[16674] | 362 | for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
|
---|
[16727] | 363 | g[i - rowIndex] = data[rows[i]];
|
---|
[16674] | 364 | }
|
---|
| 365 | }
|
---|
| 366 | }
|
---|
| 367 | }
|
---|
| 368 |
|
---|
| 369 |
|
---|
[16693] | 370 | public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
|
---|
[16674] | 371 | [ThreadStatic]
|
---|
[16682] | 372 | private Dictionary<string, Interval> intervals;
|
---|
[16674] | 373 |
|
---|
[16682] | 374 | public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals) {
|
---|
[16674] | 375 | this.intervals = intervals;
|
---|
| 376 | var code = Compile(tree);
|
---|
| 377 | Evaluate(code);
|
---|
[16682] | 378 | return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
|
---|
[16674] | 379 | }
|
---|
| 380 |
|
---|
[16682] | 381 | public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {
|
---|
| 382 | this.intervals = intervals;
|
---|
| 383 | var code = Compile(tree);
|
---|
| 384 | Evaluate(code);
|
---|
| 385 | lowerGradient = new double[paramNodes.Length];
|
---|
| 386 | upperGradient = new double[paramNodes.Length];
|
---|
| 387 | var l = code[0].value.LowerBound;
|
---|
| 388 | var u = code[0].value.UpperBound;
|
---|
| 389 | for (int i = 0; i < paramNodes.Length; ++i) {
|
---|
[16694] | 390 | if (paramNodes[i] == null) continue;
|
---|
[16738] | 391 | if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;
|
---|
| 392 | if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;
|
---|
[16682] | 393 | }
|
---|
| 394 | return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
|
---|
| 395 | }
|
---|
[16674] | 396 |
|
---|
| 397 | protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
|
---|
| 398 | instruction.value = new AlgebraicInterval(0, 0);
|
---|
| 399 | }
|
---|
| 400 |
|
---|
| 401 |
|
---|
| 402 | protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
|
---|
| 403 | instruction.dblVal = constant.Value;
|
---|
[16682] | 404 | instruction.value = new AlgebraicInterval(
|
---|
[16693] | 405 | new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0),
|
---|
| 406 | new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key
|
---|
[16682] | 407 | );
|
---|
[16674] | 408 | }
|
---|
| 409 |
|
---|
| 410 | protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
|
---|
| 411 | instruction.dblVal = variable.Weight;
|
---|
[16682] | 412 | instruction.value = new AlgebraicInterval(
|
---|
[16693] | 413 | 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)
|
---|
| 414 | high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
|
---|
[16682] | 415 | );
|
---|
[16674] | 416 | }
|
---|
| 417 |
|
---|
| 418 | protected override void LoadVariable(Instruction a) {
|
---|
| 419 | // nothing to do
|
---|
| 420 | }
|
---|
| 421 | }
|
---|
| 422 |
|
---|
| 423 | public interface IAlgebraicType<T> {
|
---|
[16693] | 424 | T Zero { get; }
|
---|
[16727] | 425 | T One { get; }
|
---|
[16693] | 426 |
|
---|
[16682] | 427 | T AssignAbs(T a); // set this to assign abs(a)
|
---|
[16674] | 428 | T Assign(T a); // assign this to same value as a (copy!)
|
---|
| 429 | T AssignNeg(T a); // set this to negative(a)
|
---|
| 430 | T AssignInv(T a); // set this to inv(a);
|
---|
| 431 | T Scale(double s); // scale this with s
|
---|
| 432 | T Add(T a); // add a to this
|
---|
| 433 | T Sub(T a); // subtract a from this
|
---|
| 434 | T Mul(T a); // multiply this with a
|
---|
| 435 | T Div(T a); // divide this by a
|
---|
| 436 | T AssignLog(T a); // set this to log a
|
---|
| 437 | T AssignExp(T a); // set this to exp(a)
|
---|
| 438 | T AssignSin(T a); // set this to sin(a)
|
---|
| 439 | T AssignCos(T a); // set this to cos(a)
|
---|
| 440 | T AssignIntPower(T a, int p);
|
---|
| 441 | T AssignIntRoot(T a, int r);
|
---|
[16682] | 442 | T AssignSgn(T a); // set this to sign(a)
|
---|
[16674] | 443 | T Clone();
|
---|
| 444 | }
|
---|
| 445 |
|
---|
[16682] | 446 | public static class Algebraic {
|
---|
[16695] | 447 | public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
|
---|
| 448 | public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
|
---|
| 449 | public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
|
---|
| 450 | public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
|
---|
| 451 | public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
|
---|
| 452 | public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
|
---|
| 453 | public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
|
---|
| 454 | public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
|
---|
| 455 | public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
|
---|
| 456 | public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
|
---|
[16674] | 457 |
|
---|
[16682] | 458 | public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
|
---|
| 459 | // ((a + b) + abs(b - a)) / 2
|
---|
| 460 | return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
|
---|
| 461 | }
|
---|
| 462 | public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
|
---|
| 463 | // ((a + b) - abs(a - b)) / 2
|
---|
| 464 | return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
|
---|
| 465 | }
|
---|
| 466 | }
|
---|
| 467 |
|
---|
| 468 |
|
---|
| 469 | // algebraic type wrapper for a double value
|
---|
[16695] | 470 | [DebuggerDisplay("{Value}")]
|
---|
[16693] | 471 | public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
|
---|
| 472 | public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
|
---|
| 473 | public static implicit operator double(AlgebraicDouble value) { return value.Value; }
|
---|
[16682] | 474 | public double Value;
|
---|
| 475 |
|
---|
[16695] | 476 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16693] | 477 | public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
|
---|
[16727] | 478 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 479 | public AlgebraicDouble One => new AlgebraicDouble(1.0);
|
---|
[16693] | 480 | public AlgebraicDouble() { }
|
---|
| 481 | public AlgebraicDouble(double value) { this.Value = value; }
|
---|
| 482 | public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
|
---|
| 483 | public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
|
---|
| 484 | public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
|
---|
| 485 | public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
|
---|
| 486 | public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
|
---|
| 487 | public AlgebraicDouble Scale(double s) { Value *= s; return this; }
|
---|
| 488 | public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
|
---|
| 489 | public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
|
---|
| 490 | public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
|
---|
| 491 | public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
|
---|
[16744] | 492 | 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] | 493 | public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
|
---|
| 494 | public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
|
---|
| 495 | public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
|
---|
| 496 | public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
|
---|
[16744] | 497 | public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }
|
---|
[16693] | 498 | public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
|
---|
[16695] | 499 |
|
---|
| 500 | public override string ToString() {
|
---|
| 501 | return Value.ToString();
|
---|
| 502 | }
|
---|
[16682] | 503 | }
|
---|
| 504 |
|
---|
[16674] | 505 | // a simple vector as an algebraic type
|
---|
[16695] | 506 | [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
|
---|
| 507 | public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
|
---|
[16674] | 508 | private double[] arr;
|
---|
| 509 | public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
|
---|
[16693] | 510 | public int Length => arr.Length;
|
---|
| 511 |
|
---|
[16695] | 512 | public AlgebraicDoubleVector(int length) { arr = new double[length]; }
|
---|
[16674] | 513 |
|
---|
[16695] | 514 | public AlgebraicDoubleVector() { }
|
---|
[16682] | 515 |
|
---|
[16674] | 516 | /// <summary>
|
---|
| 517 | ///
|
---|
| 518 | /// </summary>
|
---|
| 519 | /// <param name="arr">array is not copied</param>
|
---|
[16695] | 520 | public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
|
---|
[16674] | 521 |
|
---|
[16695] | 522 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 523 | public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
|
---|
[16727] | 524 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 525 | 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] | 526 | public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
|
---|
| 527 | public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
|
---|
| 528 | public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
|
---|
| 529 | public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
|
---|
| 530 | public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
|
---|
| 531 | public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
|
---|
| 532 | public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
|
---|
| 533 | public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
|
---|
| 534 | public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
|
---|
| 535 | public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
|
---|
| 536 | public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
|
---|
| 537 | public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
|
---|
| 538 | 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; }
|
---|
| 539 | 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; }
|
---|
| 540 | public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
|
---|
| 541 | public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
|
---|
[16674] | 542 |
|
---|
[16695] | 543 | public AlgebraicDoubleVector Clone() {
|
---|
| 544 | var v = new AlgebraicDoubleVector(this.arr.Length);
|
---|
[16693] | 545 | Array.Copy(arr, v.arr, v.arr.Length);
|
---|
| 546 | return v;
|
---|
[16674] | 547 | }
|
---|
| 548 |
|
---|
[16727] | 549 | public AlgebraicDoubleVector AssignConstant(double constantValue) {
|
---|
[16674] | 550 | for (int i = 0; i < arr.Length; ++i) {
|
---|
| 551 | arr[i] = constantValue;
|
---|
| 552 | }
|
---|
[16727] | 553 | return this;
|
---|
[16674] | 554 | }
|
---|
| 555 |
|
---|
[16693] | 556 | public void CopyTo(double[] dest, int idx, int length) {
|
---|
| 557 | Array.Copy(arr, 0, dest, idx, length);
|
---|
[16674] | 558 | }
|
---|
| 559 |
|
---|
| 560 | public void CopyFrom(double[] data, int rowIndex) {
|
---|
| 561 | Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
|
---|
| 562 | }
|
---|
[16693] | 563 | public void CopyRowTo(double[,] dest, int row) {
|
---|
| 564 | for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
|
---|
| 565 | }
|
---|
| 566 |
|
---|
| 567 | internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
|
---|
| 568 | for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
|
---|
| 569 | }
|
---|
[16695] | 570 |
|
---|
| 571 | public override string ToString() {
|
---|
| 572 | return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
|
---|
| 573 | }
|
---|
[16674] | 574 | }
|
---|
| 575 |
|
---|
[16727] | 576 |
|
---|
| 577 | /*
|
---|
[16674] | 578 | // vectors of algebraic types
|
---|
[16727] | 579 | public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
|
---|
[16674] | 580 | private T[] elems;
|
---|
| 581 |
|
---|
| 582 | public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
|
---|
| 583 |
|
---|
| 584 | public int Length => elems.Length;
|
---|
| 585 |
|
---|
[16695] | 586 | private AlgebraicVector() { }
|
---|
[16693] | 587 |
|
---|
[16695] | 588 | public AlgebraicVector(int len) { elems = new T[len]; }
|
---|
[16674] | 589 |
|
---|
| 590 | /// <summary>
|
---|
| 591 | ///
|
---|
| 592 | /// </summary>
|
---|
[16693] | 593 | /// <param name="elems">The array is copied (element-wise clone)</param>
|
---|
[16695] | 594 | public AlgebraicVector(T[] elems) {
|
---|
[16674] | 595 | this.elems = new T[elems.Length];
|
---|
| 596 | for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
|
---|
| 597 | }
|
---|
| 598 |
|
---|
| 599 | /// <summary>
|
---|
| 600 | ///
|
---|
| 601 | /// </summary>
|
---|
| 602 | /// <param name="elems">Array is not copied!</param>
|
---|
| 603 | /// <returns></returns>
|
---|
[16695] | 604 | public AlgebraicVector<T> FromArray(T[] elems) {
|
---|
| 605 | var v = new AlgebraicVector<T>();
|
---|
[16674] | 606 | v.elems = elems;
|
---|
| 607 | return v;
|
---|
| 608 | }
|
---|
| 609 |
|
---|
| 610 | public void CopyTo(T[] dest) {
|
---|
| 611 | if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
|
---|
| 612 | Array.Copy(elems, dest, dest.Length);
|
---|
| 613 | }
|
---|
| 614 |
|
---|
[16695] | 615 | public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
|
---|
[16674] | 616 |
|
---|
| 617 |
|
---|
[16695] | 618 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 619 | public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
|
---|
[16727] | 620 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 621 | 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] | 622 | public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
|
---|
| 623 | public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
|
---|
| 624 | public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
|
---|
| 625 | public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
|
---|
| 626 | public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
|
---|
| 627 | public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
|
---|
| 628 | public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
|
---|
| 629 | public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
|
---|
| 630 | public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
|
---|
| 631 | public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
|
---|
| 632 | public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
|
---|
| 633 | public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
|
---|
| 634 | public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
|
---|
| 635 | 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; }
|
---|
| 636 | 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; }
|
---|
| 637 | public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
|
---|
| 638 | public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
|
---|
[16674] | 639 | }
|
---|
| 640 |
|
---|
[16727] | 641 | */
|
---|
[16674] | 642 |
|
---|
[16727] | 643 |
|
---|
[16682] | 644 | /// <summary>
|
---|
| 645 | /// A sparse vector of algebraic types. Elements are accessed via a key of type K
|
---|
| 646 | /// </summary>
|
---|
| 647 | /// <typeparam name="K">Key type</typeparam>
|
---|
| 648 | /// <typeparam name="T">Element type</typeparam>
|
---|
[16695] | 649 | [DebuggerDisplay("SparseVector: {ToString()}")]
|
---|
| 650 | public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
|
---|
| 651 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16682] | 652 | private Dictionary<K, T> elems;
|
---|
| 653 | public IReadOnlyDictionary<K, T> Elements => elems;
|
---|
| 654 |
|
---|
[16693] | 655 |
|
---|
[16695] | 656 | public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
|
---|
[16682] | 657 | elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
|
---|
| 658 | }
|
---|
| 659 |
|
---|
[16693] | 660 | /// <summary>
|
---|
| 661 | ///
|
---|
| 662 | /// </summary>
|
---|
| 663 | /// <param name="keys"></param>
|
---|
| 664 | /// <param name="values">values are cloned</param>
|
---|
[16695] | 665 | public AlgebraicSparseVector(K[] keys, T[] values) {
|
---|
[16682] | 666 | if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
|
---|
| 667 | elems = new Dictionary<K, T>(keys.Length);
|
---|
| 668 | for (int i = 0; i < keys.Length; ++i) {
|
---|
[16693] | 669 | elems.Add(keys[i], values[i].Clone());
|
---|
[16682] | 670 | }
|
---|
| 671 | }
|
---|
| 672 |
|
---|
[16695] | 673 | public AlgebraicSparseVector() {
|
---|
[16682] | 674 | this.elems = new Dictionary<K, T>();
|
---|
| 675 | }
|
---|
| 676 |
|
---|
[16695] | 677 | // keep only elements from a
|
---|
| 678 | private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
|
---|
| 679 | // remove elems from this which don't occur in a
|
---|
| 680 | List<K> keysToRemove = new List<K>();
|
---|
| 681 | foreach (var kvp in elems) {
|
---|
| 682 | if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);
|
---|
| 683 | }
|
---|
| 684 | foreach (var o in keysToRemove) elems.Remove(o); // -> zero
|
---|
[16682] | 685 |
|
---|
[16695] | 686 | foreach (var kvp in a.elems) {
|
---|
| 687 | if (elems.TryGetValue(kvp.Key, out T value))
|
---|
| 688 | mapAssign(kvp.Value, value);
|
---|
| 689 | else
|
---|
| 690 | elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));
|
---|
| 691 | }
|
---|
| 692 | }
|
---|
[16682] | 693 |
|
---|
[16695] | 694 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 695 | public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
|
---|
[16727] | 696 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 697 | public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
|
---|
[16695] | 698 |
|
---|
| 699 | public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
|
---|
| 700 | public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
|
---|
| 701 |
|
---|
| 702 | public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
|
---|
| 703 | public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
|
---|
| 704 | public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
|
---|
| 705 | public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
|
---|
| 706 | public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
|
---|
| 707 | public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
|
---|
| 708 | public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
|
---|
| 709 | public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
|
---|
| 710 | public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
|
---|
| 711 | public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
|
---|
| 712 | public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
|
---|
[16727] | 713 | public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
|
---|
| 714 | foreach (var kvp in a.elems) {
|
---|
| 715 | if (elems.TryGetValue(kvp.Key, out T value))
|
---|
| 716 | value.Add(kvp.Value);
|
---|
| 717 | else
|
---|
| 718 | elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
|
---|
| 719 | }
|
---|
| 720 | return this;
|
---|
| 721 | }
|
---|
[16695] | 722 |
|
---|
[16727] | 723 | public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
|
---|
| 724 | foreach (var kvp in a.elems) {
|
---|
| 725 | if (elems.TryGetValue(kvp.Key, out T value))
|
---|
| 726 | value.Sub(kvp.Value);
|
---|
| 727 | else
|
---|
| 728 | elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
|
---|
| 729 | }
|
---|
| 730 | return this;
|
---|
| 731 | }
|
---|
| 732 |
|
---|
| 733 | public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
|
---|
| 734 | var keys = elems.Keys.ToArray();
|
---|
| 735 | foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
|
---|
| 736 | foreach (var kvp in a.elems) {
|
---|
| 737 | if (elems.TryGetValue(kvp.Key, out T value))
|
---|
| 738 | value.Mul(kvp.Value); // this * a
|
---|
| 739 | }
|
---|
| 740 | return this;
|
---|
| 741 | }
|
---|
| 742 |
|
---|
| 743 | public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
|
---|
| 744 | return Mul(a.Clone().Inv());
|
---|
| 745 | }
|
---|
| 746 |
|
---|
[16695] | 747 | public AlgebraicSparseVector<K, T> Clone() {
|
---|
| 748 | return new AlgebraicSparseVector<K, T>(this);
|
---|
[16682] | 749 | }
|
---|
[16695] | 750 |
|
---|
| 751 | public override string ToString() {
|
---|
| 752 | return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
|
---|
| 753 | }
|
---|
[16682] | 754 | }
|
---|
| 755 |
|
---|
[16695] | 756 | [DebuggerDisplay("[{low.Value}..{high.Value}]")]
|
---|
[16674] | 757 | public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
|
---|
[16695] | 758 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16693] | 759 | private MultivariateDual<AlgebraicDouble> low;
|
---|
[16695] | 760 | public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
|
---|
| 761 |
|
---|
| 762 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16693] | 763 | private MultivariateDual<AlgebraicDouble> high;
|
---|
| 764 | public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
|
---|
[16674] | 765 |
|
---|
[16693] | 766 |
|
---|
[16674] | 767 | public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
|
---|
| 768 |
|
---|
[16693] | 769 | public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
|
---|
[16682] | 770 | this.low = low.Clone();
|
---|
| 771 | this.high = high.Clone();
|
---|
| 772 | }
|
---|
| 773 |
|
---|
[16674] | 774 | public AlgebraicInterval(double low, double high) {
|
---|
[16693] | 775 | this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
|
---|
| 776 | this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
|
---|
[16674] | 777 | }
|
---|
| 778 |
|
---|
[16695] | 779 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16693] | 780 | public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
|
---|
[16727] | 781 | public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
|
---|
[16674] | 782 | public AlgebraicInterval Add(AlgebraicInterval a) {
|
---|
[16682] | 783 | low.Add(a.low);
|
---|
| 784 | high.Add(a.high);
|
---|
[16674] | 785 | return this;
|
---|
| 786 | }
|
---|
| 787 |
|
---|
[16693] | 788 | public AlgebraicInterval Mul(AlgebraicInterval a) {
|
---|
| 789 | var v1 = low.Clone().Mul(a.low);
|
---|
| 790 | var v2 = low.Clone().Mul(a.high);
|
---|
| 791 | var v3 = high.Clone().Mul(a.low);
|
---|
| 792 | var v4 = high.Clone().Mul(a.high);
|
---|
| 793 |
|
---|
| 794 | low = Algebraic.Min(Algebraic.Min(v1, v2), Algebraic.Min(v3, v4));
|
---|
| 795 | high = Algebraic.Max(Algebraic.Max(v1, v2), Algebraic.Max(v3, v4));
|
---|
| 796 | return this;
|
---|
| 797 | }
|
---|
| 798 |
|
---|
[16674] | 799 | public AlgebraicInterval Assign(AlgebraicInterval a) {
|
---|
| 800 | low = a.low;
|
---|
| 801 | high = a.high;
|
---|
| 802 | return this;
|
---|
| 803 | }
|
---|
| 804 |
|
---|
| 805 | public AlgebraicInterval AssignCos(AlgebraicInterval a) {
|
---|
[16693] | 806 | return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
|
---|
[16674] | 807 | }
|
---|
| 808 |
|
---|
| 809 | public AlgebraicInterval Div(AlgebraicInterval a) {
|
---|
[16682] | 810 | if (a.Contains(0.0)) {
|
---|
| 811 | if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
|
---|
[16693] | 812 | low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
|
---|
| 813 | high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
|
---|
[16682] | 814 | } else if (a.low.Value.Value.IsAlmost(0.0))
|
---|
[16693] | 815 | Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
|
---|
[16682] | 816 | else
|
---|
[16693] | 817 | Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
|
---|
[16674] | 818 | } else {
|
---|
[16682] | 819 | Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
|
---|
[16674] | 820 | }
|
---|
| 821 | return this;
|
---|
| 822 | }
|
---|
| 823 |
|
---|
| 824 | public AlgebraicInterval AssignExp(AlgebraicInterval a) {
|
---|
[16682] | 825 | low.AssignExp(a.low);
|
---|
| 826 | high.AssignExp(a.high);
|
---|
[16674] | 827 | return this;
|
---|
| 828 | }
|
---|
| 829 |
|
---|
| 830 | public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
|
---|
[16693] | 831 | if (p < 0) { // x^-3 == 1/(x^3)
|
---|
| 832 | AssignIntPower(a, -p);
|
---|
| 833 | return AssignInv(this);
|
---|
[16744] | 834 | } else if (p == 0) {
|
---|
| 835 | if (a.Contains(0.0)) {
|
---|
| 836 | // => 0^0 = 0 ; might not be relevant
|
---|
| 837 | low = new MultivariateDual<AlgebraicDouble>(0.0);
|
---|
| 838 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
| 839 | return this;
|
---|
| 840 | } else {
|
---|
| 841 | // => 1
|
---|
| 842 | low = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
| 843 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
| 844 | return this;
|
---|
| 845 | }
|
---|
| 846 | } else if (p == 1) return this;
|
---|
| 847 | else if (p % 2 == 0) {
|
---|
[16693] | 848 | // p is even => interval must be positive
|
---|
[16744] | 849 | if (a.Contains(0.0)) {
|
---|
| 850 | low = new MultivariateDual<AlgebraicDouble>(0.0);
|
---|
| 851 | high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
|
---|
[16693] | 852 | } else {
|
---|
[16727] | 853 | var lowPower = a.low.IntPower(p);
|
---|
| 854 | var highPower = a.high.IntPower(p);
|
---|
[16693] | 855 | low = Algebraic.Min(lowPower, highPower);
|
---|
| 856 | high = Algebraic.Max(lowPower, highPower);
|
---|
| 857 | }
|
---|
[16744] | 858 | } else {
|
---|
| 859 | // p is uneven
|
---|
| 860 | if (a.Contains(0.0)) {
|
---|
| 861 | low.AssignIntPower(a.low, p);
|
---|
| 862 | high.AssignIntPower(a.high, p);
|
---|
| 863 | } else {
|
---|
| 864 | var lowPower = a.low.IntPower(p);
|
---|
| 865 | var highPower = a.high.IntPower(p);
|
---|
| 866 | low = Algebraic.Min(lowPower, highPower);
|
---|
| 867 | high = Algebraic.Max(lowPower, highPower);
|
---|
| 868 | }
|
---|
[16693] | 869 | }
|
---|
[16744] | 870 | return this;
|
---|
[16674] | 871 | }
|
---|
| 872 |
|
---|
| 873 | public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
|
---|
[16693] | 874 | if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
|
---|
| 875 | if (r == 1) return this;
|
---|
| 876 | if (r < 0) {
|
---|
| 877 | // x^ (-1/2) = 1 / (x^(1/2))
|
---|
| 878 | AssignIntRoot(a, -r);
|
---|
| 879 | return AssignInv(this);
|
---|
| 880 | } else {
|
---|
| 881 | // root only defined for positive arguments
|
---|
| 882 | if (a.LowerBound.Value.Value < 0) {
|
---|
| 883 | low = new MultivariateDual<AlgebraicDouble>(double.NaN);
|
---|
| 884 | high = new MultivariateDual<AlgebraicDouble>(double.NaN);
|
---|
| 885 | return this;
|
---|
| 886 | } else {
|
---|
| 887 | low.AssignIntRoot(a.low, r);
|
---|
| 888 | high.AssignIntRoot(a.high, r);
|
---|
| 889 | return this;
|
---|
| 890 | }
|
---|
| 891 | }
|
---|
[16674] | 892 | }
|
---|
| 893 |
|
---|
| 894 | public AlgebraicInterval AssignInv(AlgebraicInterval a) {
|
---|
[16693] | 895 | low = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
| 896 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
[16674] | 897 | return Div(a);
|
---|
| 898 | }
|
---|
| 899 |
|
---|
| 900 | public AlgebraicInterval AssignLog(AlgebraicInterval a) {
|
---|
[16682] | 901 | low.AssignLog(a.low);
|
---|
| 902 | high.AssignLog(a.high);
|
---|
[16674] | 903 | return this;
|
---|
| 904 | }
|
---|
| 905 |
|
---|
[16693] | 906 | public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
|
---|
| 907 | low.AssignNeg(a.high);
|
---|
| 908 | high.AssignNeg(a.low);
|
---|
[16674] | 909 | return this;
|
---|
| 910 | }
|
---|
| 911 |
|
---|
| 912 | public AlgebraicInterval Scale(double s) {
|
---|
[16682] | 913 | low.Scale(s);
|
---|
| 914 | high.Scale(s);
|
---|
[16674] | 915 | if (s < 0) {
|
---|
[16682] | 916 | var t = low;
|
---|
| 917 | low = high;
|
---|
| 918 | high = t;
|
---|
[16674] | 919 | }
|
---|
| 920 | return this;
|
---|
| 921 | }
|
---|
| 922 |
|
---|
| 923 | public AlgebraicInterval AssignSin(AlgebraicInterval a) {
|
---|
[16693] | 924 | if (Math.Abs(a.UpperBound.Value.Value - a.LowerBound.Value.Value) >= Math.PI * 2) {
|
---|
| 925 | low = new MultivariateDual<AlgebraicDouble>(-1.0);
|
---|
| 926 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
| 927 | }
|
---|
| 928 |
|
---|
| 929 | //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)
|
---|
| 930 | double Pihalf = Math.PI / 2;
|
---|
| 931 | var scaled = this.Clone().Scale(1.0 / Pihalf);
|
---|
| 932 | //move to positive scale
|
---|
| 933 | if (scaled.LowerBound.Value.Value < 0) {
|
---|
| 934 | int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1;
|
---|
| 935 | scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4));
|
---|
| 936 | }
|
---|
| 937 |
|
---|
| 938 | double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0;
|
---|
| 939 | double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0;
|
---|
| 940 | if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0;
|
---|
| 941 | List<double> sinValues = new List<double>();
|
---|
| 942 | sinValues.Add(Math.Sin(scaledLowerBound * Pihalf));
|
---|
| 943 | sinValues.Add(Math.Sin(scaledUpperBound * Pihalf));
|
---|
| 944 |
|
---|
| 945 | int startValue = (int)Math.Ceiling(scaledLowerBound);
|
---|
| 946 | while (startValue < scaledUpperBound) {
|
---|
| 947 | sinValues.Add(Math.Sin(startValue * Pihalf));
|
---|
| 948 | startValue += 1;
|
---|
| 949 | }
|
---|
| 950 |
|
---|
| 951 | low = new MultivariateDual<AlgebraicDouble>(sinValues.Min());
|
---|
| 952 | high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());
|
---|
| 953 | return this;
|
---|
[16674] | 954 | }
|
---|
| 955 |
|
---|
| 956 | public AlgebraicInterval Sub(AlgebraicInterval a) {
|
---|
[16682] | 957 | // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
|
---|
| 958 | low.Sub(a.high);
|
---|
| 959 | high.Sub(a.low);
|
---|
[16674] | 960 | return this;
|
---|
| 961 | }
|
---|
| 962 |
|
---|
| 963 | public AlgebraicInterval Clone() {
|
---|
| 964 | return new AlgebraicInterval(low, high);
|
---|
| 965 | }
|
---|
| 966 |
|
---|
| 967 | public bool Contains(double val) {
|
---|
[16682] | 968 | return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
|
---|
[16674] | 969 | }
|
---|
[16682] | 970 |
|
---|
| 971 | public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
|
---|
| 972 | if (a.Contains(0.0)) {
|
---|
| 973 | var abslow = a.low.Clone().Abs();
|
---|
| 974 | var abshigh = a.high.Clone().Abs();
|
---|
| 975 | a.high.Assign(Algebraic.Max(abslow, abshigh));
|
---|
[16693] | 976 | a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
|
---|
[16682] | 977 | } else {
|
---|
| 978 | var abslow = a.low.Clone().Abs();
|
---|
| 979 | var abshigh = a.high.Clone().Abs();
|
---|
| 980 | a.low.Assign(Algebraic.Min(abslow, abshigh));
|
---|
| 981 | a.high.Assign(Algebraic.Max(abslow, abshigh));
|
---|
| 982 | }
|
---|
| 983 | return this;
|
---|
| 984 | }
|
---|
| 985 |
|
---|
| 986 | public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
|
---|
[16693] | 987 | low.AssignSgn(a.low);
|
---|
| 988 | high.AssignSgn(a.high);
|
---|
| 989 | return this;
|
---|
[16682] | 990 | }
|
---|
[16674] | 991 | }
|
---|
| 992 |
|
---|
| 993 | public class Dual<V> : IAlgebraicType<Dual<V>>
|
---|
| 994 | where V : IAlgebraicType<V> {
|
---|
[16695] | 995 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16674] | 996 | private V v;
|
---|
[16694] | 997 | public V Value => v;
|
---|
| 998 |
|
---|
[16695] | 999 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16674] | 1000 | private V dv;
|
---|
[16682] | 1001 | public V Derivative => dv;
|
---|
| 1002 |
|
---|
[16694] | 1003 | public Dual(V v, V dv) { this.v = v; this.dv = dv; }
|
---|
[16693] | 1004 |
|
---|
[16695] | 1005 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16693] | 1006 | public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
|
---|
[16727] | 1007 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 1008 | public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
|
---|
[16693] | 1009 |
|
---|
[16694] | 1010 | public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
|
---|
| 1011 | public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
|
---|
| 1012 | public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
|
---|
| 1013 | public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
|
---|
| 1014 | public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
|
---|
| 1015 | 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] | 1016 |
|
---|
[16694] | 1017 | // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
|
---|
[16674] | 1018 | public Dual<V> Mul(Dual<V> a) {
|
---|
[16694] | 1019 | var t1 = a.dv.Clone().Mul(v);
|
---|
| 1020 | var t2 = dv.Clone().Mul(a.v);
|
---|
[16674] | 1021 | dv.Assign(t1).Add(t2);
|
---|
| 1022 |
|
---|
| 1023 | v.Mul(a.v);
|
---|
| 1024 | return this;
|
---|
| 1025 | }
|
---|
[16694] | 1026 | public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
|
---|
[16674] | 1027 |
|
---|
[16694] | 1028 | 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)'
|
---|
| 1029 | 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] | 1030 |
|
---|
[16694] | 1031 | 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; }
|
---|
| 1032 | 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] | 1033 |
|
---|
[16694] | 1034 | public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
|
---|
| 1035 | public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
|
---|
[16674] | 1036 |
|
---|
[16694] | 1037 | 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)|
|
---|
| 1038 | public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
|
---|
[16682] | 1039 |
|
---|
[16694] | 1040 | public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
|
---|
[16682] | 1041 |
|
---|
[16674] | 1042 | }
|
---|
| 1043 |
|
---|
[16682] | 1044 | /// <summary>
|
---|
| 1045 | /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.
|
---|
| 1046 | /// </summary>
|
---|
| 1047 | /// <typeparam name="V"></typeparam>
|
---|
[16695] | 1048 | [DebuggerDisplay("v={Value}; dv={dv}")]
|
---|
[16682] | 1049 | public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
|
---|
[16695] | 1050 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
[16674] | 1051 | private V v;
|
---|
| 1052 | public V Value => v;
|
---|
| 1053 |
|
---|
[16695] | 1054 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
| 1055 | private AlgebraicSparseVector<object, V> dv;
|
---|
| 1056 | public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key
|
---|
[16674] | 1057 |
|
---|
[16694] | 1058 | private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }
|
---|
[16674] | 1059 |
|
---|
[16682] | 1060 | /// <summary>
|
---|
| 1061 | /// Constructor without partial derivative
|
---|
| 1062 | /// </summary>
|
---|
| 1063 | /// <param name="v"></param>
|
---|
[16695] | 1064 | public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }
|
---|
[16674] | 1065 |
|
---|
[16682] | 1066 | /// <summary>
|
---|
| 1067 | /// Constructor for multiple partial derivatives
|
---|
| 1068 | /// </summary>
|
---|
| 1069 | /// <param name="v"></param>
|
---|
| 1070 | /// <param name="keys"></param>
|
---|
| 1071 | /// <param name="dv"></param>
|
---|
[16695] | 1072 | public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }
|
---|
[16674] | 1073 |
|
---|
[16682] | 1074 | /// <summary>
|
---|
| 1075 | /// Constructor for a single partial derivative
|
---|
| 1076 | /// </summary>
|
---|
| 1077 | /// <param name="v"></param>
|
---|
| 1078 | /// <param name="key"></param>
|
---|
| 1079 | /// <param name="dv"></param>
|
---|
[16695] | 1080 | public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
|
---|
[16682] | 1081 |
|
---|
[16693] | 1082 | /// <summary>
|
---|
[16694] | 1083 | /// Constructor with a given value and gradient. For internal use.
|
---|
[16693] | 1084 | /// </summary>
|
---|
[16694] | 1085 | /// <param name="v">The value (not cloned).</param>
|
---|
| 1086 | /// <param name="gradient">The gradient (not cloned).</param>
|
---|
[16695] | 1087 | internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
|
---|
[16693] | 1088 |
|
---|
[16694] | 1089 | public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
|
---|
[16682] | 1090 |
|
---|
[16693] | 1091 | public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
|
---|
[16727] | 1092 | public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
|
---|
[16693] | 1093 |
|
---|
[16694] | 1094 | public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
|
---|
[16693] | 1095 |
|
---|
[16694] | 1096 | public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
|
---|
| 1097 | public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
|
---|
| 1098 | public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
|
---|
[16682] | 1099 | public MultivariateDual<V> Mul(MultivariateDual<V> a) {
|
---|
[16674] | 1100 | // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
|
---|
| 1101 | var t1 = a.dv.Clone().Scale(v);
|
---|
| 1102 | var t2 = dv.Clone().Scale(a.v);
|
---|
[16682] | 1103 | dv.Assign(t1).Add(t2);
|
---|
[16674] | 1104 |
|
---|
| 1105 | v.Mul(a.v);
|
---|
| 1106 | return this;
|
---|
| 1107 | }
|
---|
[16694] | 1108 | public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
|
---|
| 1109 | public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
|
---|
| 1110 | 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] | 1111 |
|
---|
[16694] | 1112 | public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
|
---|
| 1113 | public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
|
---|
[16674] | 1114 |
|
---|
[16694] | 1115 | 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; }
|
---|
| 1116 | 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] | 1117 |
|
---|
[16694] | 1118 | 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)'
|
---|
| 1119 | 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] | 1120 |
|
---|
[16694] | 1121 | 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
|
---|
| 1122 | public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;
|
---|
[16674] | 1123 | }
|
---|
[16682] | 1124 | } |
---|