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