Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3040_VectorBasedGP/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/SymbolicDataAnalysisExpressionTreeVectorInterpreter.cs

Last change on this file was 18234, checked in by pfleck, 2 years ago

#3040 Fixed vector-unrolling AutoDiff conversion.
File size: 53.7 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Linq;
25using HeuristicLab.Common;
26using HeuristicLab.Core;
27using HeuristicLab.Data;
28using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
29using HeuristicLab.Parameters;
30using HEAL.Attic;
31using MathNet.Numerics;
32using MathNet.Numerics.Statistics;
33using DoubleVector = MathNet.Numerics.LinearAlgebra.Vector<double>;
34
35namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
36  [StorableType("DE68A1D9-5AFC-4DDD-AB62-29F3B8FC28E0")]
37  [Item("SymbolicDataAnalysisExpressionTreeVectorInterpreter", "Interpreter for symbolic expression trees including vector arithmetic.")]
38  public class SymbolicDataAnalysisExpressionTreeVectorInterpreter : ParameterizedNamedItem, ISymbolicDataAnalysisExpressionTreeInterpreter {
39    [StorableType("2612504E-AD5F-4AE2-B60E-98A5AB59E164")]
40    public enum Aggregation {
41      Mean,
42      Median,
43      Sum,
44      First,
45      L1Norm,
46      L2Norm,
47      NaN,
48      Exception
49    }
50    public static double Aggregate(Aggregation aggregation, DoubleVector vector) {
51      switch (aggregation) {
52        case Aggregation.Mean: return Statistics.Mean(vector);
53        case Aggregation.Median: return Statistics.Median(vector);
54        case Aggregation.Sum: return vector.Sum();
55        case Aggregation.First: return vector.First();
56        case Aggregation.L1Norm: return vector.L1Norm();
57        case Aggregation.L2Norm: return vector.L2Norm();
58        case Aggregation.NaN: return double.NaN;
59        case Aggregation.Exception: throw new InvalidOperationException("Result of the tree is not a scalar.");
60        default: throw new ArgumentOutOfRangeException(nameof(aggregation), aggregation, null);
61      }
62    }
63
64    [StorableType("73DCBB45-916F-4139-8ADC-57BA610A1B66")]
65    public enum VectorLengthStrategy {
66      ExceptionIfDifferent,
67      FillShorterWithNaN,
68      FillShorterWithNeutralElement,
69      CutLonger,
70      ResampleToLonger,
71      ResampleToShorter,
72      CycleShorter
73    }
74
75    #region Implementation VectorLengthStrategy
76    public static (DoubleVector, DoubleVector) ExceptionIfDifferent(DoubleVector lhs, DoubleVector rhs) {
77      if (lhs.Count != rhs.Count)
78        throw new InvalidOperationException($"Vector Lengths incompatible ({lhs.Count} vs. {rhs.Count}");
79      return (lhs, rhs);
80    }
81
82    public static (DoubleVector, DoubleVector) FillShorter(DoubleVector lhs, DoubleVector rhs, double fillElement) {
83      var targetLength = Math.Max(lhs.Count, rhs.Count);
84
85      DoubleVector PadVector(DoubleVector v) {
86        if (v.Count == targetLength) return v;
87        var p = DoubleVector.Build.Dense(targetLength, fillElement);
88        v.CopySubVectorTo(p, 0, 0, v.Count);
89        return p;
90      }
91
92      return (PadVector(lhs), PadVector(rhs));
93    }
94
95    public static (DoubleVector, DoubleVector) CutLonger(DoubleVector lhs, DoubleVector rhs) {
96      var targetLength = Math.Min(lhs.Count, rhs.Count);
97
98      DoubleVector CutVector(DoubleVector v) {
99        if (v.Count == targetLength) return v;
100        return v.SubVector(0, targetLength);
101      }
102
103      return (CutVector(lhs), CutVector(rhs));
104    }
105
106    private static DoubleVector ResampleToLength(DoubleVector v, int targetLength) {
107      if (v.Count == targetLength) return v;
108
109      var indices = Enumerable.Range(0, v.Count).Select(x => (double)x);
110      var interpolation = Interpolate.Linear(indices, v);
111
112      var resampledIndices = Enumerable.Range(0, targetLength).Select(i => (double)i / targetLength * v.Count);
113      var interpolatedValues = resampledIndices.Select(interpolation.Interpolate);
114
115      return DoubleVector.Build.DenseOfEnumerable(interpolatedValues);
116    }
117    public static (DoubleVector, DoubleVector) ResampleToLonger(DoubleVector lhs, DoubleVector rhs) {
118      var maxLength = Math.Max(lhs.Count, rhs.Count);
119      return (ResampleToLength(lhs, maxLength), ResampleToLength(rhs, maxLength));
120    }
121    public static (DoubleVector, DoubleVector) ResampleToShorter(DoubleVector lhs, DoubleVector rhs) {
122      var minLength = Math.Min(lhs.Count, rhs.Count);
123      return (ResampleToLength(lhs, minLength), ResampleToLength(rhs, minLength));
124    }
125
126    public static (DoubleVector, DoubleVector) CycleShorter(DoubleVector lhs, DoubleVector rhs) {
127      var targetLength = Math.Max(lhs.Count, rhs.Count);
128
129      DoubleVector CycleVector(DoubleVector v) {
130        if (v.Count == targetLength) return v;
131        var cycledValues = Enumerable.Range(0, targetLength).Select(i => v[i % v.Count]);
132        return DoubleVector.Build.DenseOfEnumerable(cycledValues);
133      }
134
135      return (CycleVector(lhs), CycleVector(rhs));
136    }
137    #endregion
138
139    public static (DoubleVector lhs, DoubleVector rhs) ApplyVectorLengthStrategy(VectorLengthStrategy strategy, DoubleVector lhs, DoubleVector rhs,
140      double neutralElement = double.NaN) {
141
142      switch (strategy) {
143        case VectorLengthStrategy.ExceptionIfDifferent: return ExceptionIfDifferent(lhs, rhs);
144        case VectorLengthStrategy.FillShorterWithNaN: return FillShorter(lhs, rhs, double.NaN);
145        case VectorLengthStrategy.FillShorterWithNeutralElement: return FillShorter(lhs, rhs, neutralElement);
146        case VectorLengthStrategy.CutLonger: return CutLonger(lhs, rhs);
147        case VectorLengthStrategy.ResampleToLonger: return ResampleToLonger(lhs, rhs);
148        case VectorLengthStrategy.ResampleToShorter: return ResampleToShorter(lhs, rhs);
149        case VectorLengthStrategy.CycleShorter: return CycleShorter(lhs, rhs);
150        default: throw new ArgumentOutOfRangeException(nameof(strategy), strategy, null);
151      }
152    }
153
154    #region Aggregation Symbols
155    private static Type[] AggregationSymbols = new[] {
156      typeof(Sum), typeof(Mean), typeof(Length), typeof(StandardDeviation), typeof(Variance),
157      typeof(EuclideanDistance), typeof(Covariance)
158    };
159    #endregion
160
161    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
162    private const string FinalAggregationParameterName = "FinalAggregation";
163    private const string DifferentVectorLengthStrategyParameterName = "DifferentVectorLengthStrategy";
164
165    public override bool CanChangeName {
166      get { return false; }
167    }
168
169    public override bool CanChangeDescription {
170      get { return false; }
171    }
172
173    #region parameter properties
174    public IFixedValueParameter<IntValue> EvaluatedSolutionsParameter {
175      get { return (IFixedValueParameter<IntValue>)Parameters[EvaluatedSolutionsParameterName]; }
176    }
177    public IFixedValueParameter<EnumValue<Aggregation>> FinalAggregationParameter {
178      get { return (IFixedValueParameter<EnumValue<Aggregation>>)Parameters[FinalAggregationParameterName]; }
179    }
180    public IFixedValueParameter<EnumValue<VectorLengthStrategy>> DifferentVectorLengthStrategyParameter {
181      get { return (IFixedValueParameter<EnumValue<VectorLengthStrategy>>)Parameters[DifferentVectorLengthStrategyParameterName]; }
182    }
183    #endregion
184
185    #region properties
186    public int EvaluatedSolutions {
187      get { return EvaluatedSolutionsParameter.Value.Value; }
188      set { EvaluatedSolutionsParameter.Value.Value = value; }
189    }
190    public Aggregation FinalAggregation {
191      get { return FinalAggregationParameter.Value.Value; }
192      set { FinalAggregationParameter.Value.Value = value; }
193    }
194    public VectorLengthStrategy DifferentVectorLengthStrategy {
195      get { return DifferentVectorLengthStrategyParameter.Value.Value; }
196      set { DifferentVectorLengthStrategyParameter.Value.Value = value; }
197    }
198    #endregion
199
200    [StorableConstructor]
201    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(StorableConstructorFlag _) : base(_) { }
202
203    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(SymbolicDataAnalysisExpressionTreeVectorInterpreter original, Cloner cloner)
204      : base(original, cloner) { }
205
206    public override IDeepCloneable Clone(Cloner cloner) {
207      return new SymbolicDataAnalysisExpressionTreeVectorInterpreter(this, cloner);
208    }
209
210    public SymbolicDataAnalysisExpressionTreeVectorInterpreter()
211      : this("SymbolicDataAnalysisExpressionTreeVectorInterpreter", "Interpreter for symbolic expression trees including vector arithmetic.") { }
212
213    protected SymbolicDataAnalysisExpressionTreeVectorInterpreter(string name, string description)
214      : base(name, description) {
215      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName, "A counter for the total number of solutions the interpreter has evaluated", new IntValue(0)));
216      Parameters.Add(new FixedValueParameter<EnumValue<Aggregation>>(FinalAggregationParameterName, "If root node of the expression tree results in a Vector it is aggregated according to this parameter", new EnumValue<Aggregation>(Aggregation.Mean)));
217      Parameters.Add(new FixedValueParameter<EnumValue<VectorLengthStrategy>>(DifferentVectorLengthStrategyParameterName, "", new EnumValue<VectorLengthStrategy>(VectorLengthStrategy.ExceptionIfDifferent)));
218    }
219
220    [StorableHook(HookType.AfterDeserialization)]
221    private void AfterDeserialization() {
222      if (!Parameters.ContainsKey(FinalAggregationParameterName)) {
223        Parameters.Add(new FixedValueParameter<EnumValue<Aggregation>>(FinalAggregationParameterName, "If root node of the expression tree results in a Vector it is aggregated according to this parameter", new EnumValue<Aggregation>(Aggregation.Mean)));
224      }
225      if (!Parameters.ContainsKey(DifferentVectorLengthStrategyParameterName)) {
226        Parameters.Add(new FixedValueParameter<EnumValue<VectorLengthStrategy>>(DifferentVectorLengthStrategyParameterName, "", new EnumValue<VectorLengthStrategy>(VectorLengthStrategy.ExceptionIfDifferent)));
227      }
228    }
229
230    #region IStatefulItem
231    public void InitializeState() {
232      EvaluatedSolutions = 0;
233    }
234
235    public void ClearState() { }
236    #endregion
237
238    private readonly object syncRoot = new object();
239    public IEnumerable<double> GetSymbolicExpressionTreeValues(ISymbolicExpressionTree tree, IDataset dataset, IEnumerable<int> rows) {
240      lock (syncRoot) {
241        EvaluatedSolutions++; // increment the evaluated solutions counter
242      }
243      var state = PrepareInterpreterState(tree, dataset);
244
245      foreach (var rowEnum in rows) {
246        int row = rowEnum;
247        var result = Evaluate(dataset, ref row, state);
248        if (result.IsScalar)
249          yield return result.Scalar;
250        else if (result.IsVector) {
251          yield return Aggregate(FinalAggregation, result.Vector);
252        } else
253          yield return double.NaN;
254        state.Reset();
255      }
256    }
257
258    public IEnumerable<Dictionary<ISymbolicExpressionTreeNode, EvaluationResult>> GetIntermediateNodeValues(ISymbolicExpressionTree tree, IDataset dataset, IEnumerable<int> rows) {
259      var state = PrepareInterpreterState(tree, dataset);
260
261      foreach (var rowEnum in rows) {
262        int row = rowEnum;
263        var traceDict = new Dictionary<ISymbolicExpressionTreeNode, EvaluationResult>();
264        var result = Evaluate(dataset, ref row, state, traceDict);
265        traceDict.Add(tree.Root.GetSubtree(0), result); // Add StartSymbol
266        yield return traceDict;
267        state.Reset();
268      }
269    }
270
271    private static InterpreterState PrepareInterpreterState(ISymbolicExpressionTree tree, IDataset dataset) {
272      Instruction[] code = SymbolicExpressionTreeCompiler.Compile(tree, OpCodes.MapSymbolToOpCode);
273      int necessaryArgStackSize = 0;
274      foreach (Instruction instr in code) {
275        if (instr.opCode == OpCodes.Variable) {
276          var variableTreeNode = (VariableTreeNode)instr.dynamicNode;
277          if (dataset.VariableHasType<double>(variableTreeNode.VariableName))
278            instr.data = dataset.GetReadOnlyDoubleValues(variableTreeNode.VariableName);
279          else if (dataset.VariableHasType<DoubleVector>(variableTreeNode.VariableName))
280            instr.data = dataset.GetReadOnlyDoubleVectorValues(variableTreeNode.VariableName);
281          else throw new NotSupportedException($"Type of variable {variableTreeNode.VariableName} is not supported.");
282        } else if (instr.opCode == OpCodes.FactorVariable) {
283          var factorTreeNode = instr.dynamicNode as FactorVariableTreeNode;
284          instr.data = dataset.GetReadOnlyStringValues(factorTreeNode.VariableName);
285        } else if (instr.opCode == OpCodes.BinaryFactorVariable) {
286          var factorTreeNode = instr.dynamicNode as BinaryFactorVariableTreeNode;
287          instr.data = dataset.GetReadOnlyStringValues(factorTreeNode.VariableName);
288        } else if (instr.opCode == OpCodes.LagVariable) {
289          var laggedVariableTreeNode = (LaggedVariableTreeNode)instr.dynamicNode;
290          instr.data = dataset.GetReadOnlyDoubleValues(laggedVariableTreeNode.VariableName);
291        } else if (instr.opCode == OpCodes.VariableCondition) {
292          var variableConditionTreeNode = (VariableConditionTreeNode)instr.dynamicNode;
293          instr.data = dataset.GetReadOnlyDoubleValues(variableConditionTreeNode.VariableName);
294        } else if (instr.opCode == OpCodes.Call) {
295          necessaryArgStackSize += instr.nArguments + 1;
296        }
297      }
298      return new InterpreterState(code, necessaryArgStackSize);
299    }
300
301
302    public struct EvaluationResult {
303      public double Scalar { get; }
304      public bool IsScalar => !double.IsNaN(Scalar);
305
306      public DoubleVector Vector { get; }
307      public bool IsVector => !(Vector.Count == 1 && double.IsNaN(Vector[0]));
308
309      public bool IsNaN => !IsScalar && !IsVector;
310
311      public EvaluationResult(double scalar) {
312        Scalar = scalar;
313        Vector = NaNVector;
314      }
315      public EvaluationResult(DoubleVector vector) {
316        if (vector == null) throw new ArgumentNullException(nameof(vector));
317        Vector = vector;
318        Scalar = double.NaN;
319      }
320
321      public override string ToString() {
322        if (IsScalar) return Scalar.ToString();
323        if (IsVector) return Vector.ToVectorString();
324        return "NaN";
325      }
326
327      private static readonly DoubleVector NaNVector = DoubleVector.Build.Dense(1, double.NaN);
328      public static readonly EvaluationResult NaN = new EvaluationResult(double.NaN);
329    }
330
331    private static EvaluationResult ArithmeticApply(EvaluationResult lhs, EvaluationResult rhs,
332      Func<DoubleVector, DoubleVector, (DoubleVector, DoubleVector)> lengthStrategy,
333      Func<double, double, double> ssFunc = null,
334      Func<double, DoubleVector, DoubleVector> svFunc = null,
335      Func<DoubleVector, double, DoubleVector> vsFunc = null,
336      Func<DoubleVector, DoubleVector, DoubleVector> vvFunc = null) {
337
338      if (lhs.IsScalar && rhs.IsScalar && ssFunc != null) return new EvaluationResult(ssFunc(lhs.Scalar, rhs.Scalar));
339      if (lhs.IsScalar && rhs.IsVector && svFunc != null) return new EvaluationResult(svFunc(lhs.Scalar, rhs.Vector));
340      if (lhs.IsVector && rhs.IsScalar && vsFunc != null) return new EvaluationResult(vsFunc(lhs.Vector, rhs.Scalar));
341      if (lhs.IsVector && rhs.IsVector && vvFunc != null) {
342        if (lhs.Vector.Count == rhs.Vector.Count) {
343          return new EvaluationResult(vvFunc(lhs.Vector, rhs.Vector));
344        } else {
345          var (lhsVector, rhsVector) = lengthStrategy(lhs.Vector, rhs.Vector);
346          return new EvaluationResult(vvFunc(lhsVector, rhsVector));
347        }
348      }
349      return EvaluationResult.NaN;
350    }
351
352    private static EvaluationResult FunctionApply(EvaluationResult val,
353      Func<double, double> sFunc = null,
354      Func<DoubleVector, DoubleVector> vFunc = null) {
355      if (val.IsScalar && sFunc != null) return new EvaluationResult(sFunc(val.Scalar));
356      if (val.IsVector && vFunc != null) return new EvaluationResult(vFunc(val.Vector));
357      return EvaluationResult.NaN;
358    }
359    private static EvaluationResult AggregateApply(EvaluationResult val,
360      Func<double, double> sFunc = null,
361      Func<DoubleVector, double> vFunc = null) {
362      if (val.IsScalar && sFunc != null) return new EvaluationResult(sFunc(val.Scalar));
363      if (val.IsVector && vFunc != null) return new EvaluationResult(vFunc(val.Vector));
364      return EvaluationResult.NaN;
365    }
366
367    private static EvaluationResult WindowedAggregateApply(EvaluationResult val, WindowedSymbolTreeNode node,
368      Func<double, double> sFunc = null,
369      Func<DoubleVector, double> vFunc = null) {
370
371      // Parameters are interpreted as start and end with wrapping
372      var start = node.Offset;
373      var end = node.Length;
374
375      DoubleVector SubVector(DoubleVector v) {
376        int startIdx = (int)Math.Round(start * v.Count);
377        int endIdx = (int)Math.Round(end * v.Count);
378        int size = v.Count;
379        if (startIdx < endIdx) {
380          return v.SubVector(startIdx, count: endIdx - startIdx);
381        } else { // wrap around
382          var resultVector = DoubleVector.Build.Dense(size: size - (startIdx - endIdx));
383          v.CopySubVectorTo(resultVector, startIdx, 0, size - startIdx); // copy [startIdx:size] to [0:size-startIdx]
384          v.CopySubVectorTo(resultVector, 0, size - startIdx, endIdx);   // copy [0:endIdx]      to [size-startIdx:size]
385          return resultVector;
386        }
387      }
388
389      if (val.IsScalar && sFunc != null) return new EvaluationResult(sFunc(val.Scalar));
390      if (val.IsVector && vFunc != null) return new EvaluationResult(vFunc(SubVector(val.Vector)));
391      return EvaluationResult.NaN;
392    }
393    private static EvaluationResult WindowedFunctionApply(EvaluationResult val, IWindowedSymbolTreeNode node,
394      Func<double, double> sFunc = null,
395      Func<DoubleVector, DoubleVector> vFunc = null) {
396      // Parameters are interpreted as start and end with wrapping
397      var start = node.Offset;
398      var end = node.Length;
399
400      DoubleVector SubVector(DoubleVector v) {
401        int startIdx = (int)Math.Round(start * v.Count);
402        int endIdx = (int)Math.Round(end * v.Count);
403        int size = v.Count;
404        if (startIdx < endIdx) {
405          return v.SubVector(startIdx, count: endIdx - startIdx);
406        } else { // wrap around
407          var resultVector = DoubleVector.Build.Dense(size: size - (startIdx - endIdx));
408          v.CopySubVectorTo(resultVector, startIdx, 0, size - startIdx); // copy [startIdx:size] to [0:size-startIdx]
409          v.CopySubVectorTo(resultVector, 0, size - startIdx, endIdx);   // copy [0:endIdx]      to [size-startIdx:size]
410          return resultVector;
411        }
412      }
413
414      if (val.IsScalar && sFunc != null) return new EvaluationResult(sFunc(val.Scalar));
415      if (val.IsVector && vFunc != null) return new EvaluationResult(vFunc(SubVector(val.Vector)));
416      return EvaluationResult.NaN;
417    }
418
419    private static EvaluationResult AggregateMultipleApply(EvaluationResult lhs, EvaluationResult rhs,
420      Func<DoubleVector, DoubleVector, (DoubleVector, DoubleVector)> lengthStrategy,
421      Func<double, double, double> ssFunc = null,
422      Func<double, DoubleVector, double> svFunc = null,
423      Func<DoubleVector, double, double> vsFunc = null,
424      Func<DoubleVector, DoubleVector, double> vvFunc = null) {
425      if (lhs.IsScalar && rhs.IsScalar && ssFunc != null) return new EvaluationResult(ssFunc(lhs.Scalar, rhs.Scalar));
426      if (lhs.IsScalar && rhs.IsVector && svFunc != null) return new EvaluationResult(svFunc(lhs.Scalar, rhs.Vector));
427      if (lhs.IsVector && rhs.IsScalar && vsFunc != null) return new EvaluationResult(vsFunc(lhs.Vector, rhs.Scalar));
428      if (lhs.IsVector && rhs.IsVector && vvFunc != null) {
429        if (lhs.Vector.Count == rhs.Vector.Count) {
430          return new EvaluationResult(vvFunc(lhs.Vector, rhs.Vector));
431        } else {
432          var (lhsVector, rhsVector) = lengthStrategy(lhs.Vector, rhs.Vector);
433          return new EvaluationResult(vvFunc(lhsVector, rhsVector));
434        }
435      }
436      return EvaluationResult.NaN;
437    }
438
439    public virtual Type GetNodeType(ISymbolicExpressionTreeNode node) {
440      if (node.DataType != null)
441        return node.DataType;
442
443      if (AggregationSymbols.Contains(node.Symbol.GetType()))
444        return typeof(double);
445
446      var argumentTypes = node.Subtrees.Select(GetNodeType);
447      if (argumentTypes.Any(t => t == typeof(DoubleVector)))
448        return typeof(DoubleVector);
449
450      return typeof(double);
451    }
452
453
454    public virtual EvaluationResult Evaluate(IDataset dataset, ref int row, InterpreterState state,
455      IDictionary<ISymbolicExpressionTreeNode, EvaluationResult> traceDict = null) {
456
457      void TraceEvaluation(Instruction instr, EvaluationResult result) {
458        traceDict?.Add(instr.dynamicNode, result);
459      }
460
461      Instruction currentInstr = state.NextInstruction();
462      switch (currentInstr.opCode) {
463        case OpCodes.Add: {
464            var cur = Evaluate(dataset, ref row, state, traceDict);
465            for (int i = 1; i < currentInstr.nArguments; i++) {
466              var op = Evaluate(dataset, ref row, state, traceDict);
467              cur = ArithmeticApply(cur, op,
468                (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 0.0),
469                (s1, s2) => s1 + s2,
470                (s1, v2) => s1 + v2,
471                (v1, s2) => v1 + s2,
472                (v1, v2) => v1 + v2);
473            }
474            TraceEvaluation(currentInstr, cur);
475            return cur;
476          }
477        case OpCodes.Sub: {
478            var cur = Evaluate(dataset, ref row, state, traceDict);
479            for (int i = 1; i < currentInstr.nArguments; i++) {
480              var op = Evaluate(dataset, ref row, state, traceDict);
481              cur = ArithmeticApply(cur, op,
482                (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 0.0),
483                (s1, s2) => s1 - s2,
484                (s1, v2) => s1 - v2,
485                (v1, s2) => v1 - s2,
486                (v1, v2) => v1 - v2);
487            }
488            if (currentInstr.nArguments == 1)
489              cur = FunctionApply(cur,
490                s => -s,
491                v => -v);
492            TraceEvaluation(currentInstr, cur);
493            return cur;
494          }
495        case OpCodes.Mul: {
496            var cur = Evaluate(dataset, ref row, state, traceDict);
497            for (int i = 1; i < currentInstr.nArguments; i++) {
498              var op = Evaluate(dataset, ref row, state, traceDict);
499              cur = ArithmeticApply(cur, op,
500                (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 1.0),
501                (s1, s2) => s1 * s2,
502                (s1, v2) => s1 * v2,
503                (v1, s2) => v1 * s2,
504                (v1, v2) => v1.PointwiseMultiply(v2));
505            }
506            TraceEvaluation(currentInstr, cur);
507            return cur;
508          }
509        case OpCodes.Div: {
510            var cur = Evaluate(dataset, ref row, state, traceDict);
511            for (int i = 1; i < currentInstr.nArguments; i++) {
512              var op = Evaluate(dataset, ref row, state, traceDict);
513              cur = ArithmeticApply(cur, op,
514                (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 1.0),
515                (s1, s2) => s1 / s2,
516                (s1, v2) => s1 / v2,
517                (v1, s2) => v1 / s2,
518                (v1, v2) => v1 / v2);
519            }
520            if (currentInstr.nArguments == 1)
521              cur = FunctionApply(cur,
522                s => 1 / s,
523                v => 1 / v);
524            TraceEvaluation(currentInstr, cur);
525            return cur;
526          }
527        case OpCodes.Absolute: {
528            var cur = Evaluate(dataset, ref row, state, traceDict);
529            cur = FunctionApply(cur, Math.Abs, DoubleVector.Abs);
530            TraceEvaluation(currentInstr, cur);
531            return cur;
532          }
533        case OpCodes.Tanh: {
534            var cur = Evaluate(dataset, ref row, state, traceDict);
535            cur = FunctionApply(cur, Math.Tanh, DoubleVector.Tanh);
536            TraceEvaluation(currentInstr, cur);
537            return cur;
538          }
539        case OpCodes.Cos: {
540            var cur = Evaluate(dataset, ref row, state, traceDict);
541            cur = FunctionApply(cur, Math.Cos, DoubleVector.Cos);
542            TraceEvaluation(currentInstr, cur);
543            return cur;
544          }
545        case OpCodes.Sin: {
546            var cur = Evaluate(dataset, ref row, state, traceDict);
547            cur = FunctionApply(cur, Math.Sin, DoubleVector.Sin);
548            TraceEvaluation(currentInstr, cur);
549            return cur;
550          }
551        case OpCodes.Tan: {
552            var cur = Evaluate(dataset, ref row, state, traceDict);
553            cur = FunctionApply(cur, Math.Tan, DoubleVector.Tan);
554            TraceEvaluation(currentInstr, cur);
555            return cur;
556          }
557        case OpCodes.Square: {
558            var cur = Evaluate(dataset, ref row, state, traceDict);
559            cur = FunctionApply(cur,
560              s => Math.Pow(s, 2),
561              v => v.PointwisePower(2));
562            TraceEvaluation(currentInstr, cur);
563            return cur;
564          }
565        case OpCodes.Cube: {
566            var cur = Evaluate(dataset, ref row, state, traceDict);
567            cur = FunctionApply(cur,
568              s => Math.Pow(s, 3),
569              v => v.PointwisePower(3));
570            TraceEvaluation(currentInstr, cur);
571            return cur;
572          }
573        case OpCodes.Power: {
574            var x = Evaluate(dataset, ref row, state, traceDict);
575            var y = Evaluate(dataset, ref row, state, traceDict);
576            var cur = ArithmeticApply(x, y,
577              (lhs, rhs) => lhs.Count < rhs.Count
578                ? CutLonger(lhs, rhs)
579                : ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 1.0),
580              (s1, s2) => Math.Pow(s1, Math.Round(s2)),
581              (s1, v2) => DoubleVector.Build.Dense(v2.Count, s1).PointwisePower(DoubleVector.Round(v2)),
582              (v1, s2) => v1.PointwisePower(Math.Round(s2)),
583              (v1, v2) => v1.PointwisePower(DoubleVector.Round(v2)));
584            TraceEvaluation(currentInstr, cur);
585            return cur;
586          }
587        case OpCodes.SquareRoot: {
588            var cur = Evaluate(dataset, ref row, state, traceDict);
589            cur = FunctionApply(cur,
590              s => Math.Sqrt(s),
591              v => DoubleVector.Sqrt(v));
592            TraceEvaluation(currentInstr, cur);
593            return cur;
594          }
595        case OpCodes.CubeRoot: {
596            var cur = Evaluate(dataset, ref row, state, traceDict);
597            cur = FunctionApply(cur,
598              s => s < 0 ? -Math.Pow(-s, 1.0 / 3.0) : Math.Pow(s, 1.0 / 3.0),
599              v => v.Map(s => s < 0 ? -Math.Pow(-s, 1.0 / 3.0) : Math.Pow(s, 1.0 / 3.0)));
600            TraceEvaluation(currentInstr, cur);
601            return cur;
602          }
603        case OpCodes.Root: {
604            var x = Evaluate(dataset, ref row, state, traceDict);
605            var y = Evaluate(dataset, ref row, state, traceDict);
606            var cur = ArithmeticApply(x, y,
607              (lhs, rhs) => lhs.Count < rhs.Count
608                ? CutLonger(lhs, rhs)
609                : ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 1.0),
610              (s1, s2) => Math.Pow(s1, 1.0 / Math.Round(s2)),
611              (s1, v2) => DoubleVector.Build.Dense(v2.Count, s1).PointwisePower(1.0 / DoubleVector.Round(v2)),
612              (v1, s2) => v1.PointwisePower(1.0 / Math.Round(s2)),
613              (v1, v2) => v1.PointwisePower(1.0 / DoubleVector.Round(v2)));
614            TraceEvaluation(currentInstr, cur);
615            return cur;
616          }
617        case OpCodes.Exp: {
618            var cur = Evaluate(dataset, ref row, state, traceDict);
619            cur = FunctionApply(cur,
620              s => Math.Exp(s),
621              v => DoubleVector.Exp(v));
622            TraceEvaluation(currentInstr, cur);
623            return cur;
624          }
625        case OpCodes.Log: {
626            var cur = Evaluate(dataset, ref row, state, traceDict);
627            cur = FunctionApply(cur,
628              s => Math.Log(s),
629              v => DoubleVector.Log(v));
630            TraceEvaluation(currentInstr, cur);
631            return cur;
632          }
633        case OpCodes.Sum: {
634            var cur = Evaluate(dataset, ref row, state, traceDict);
635            cur = AggregateApply(cur,
636              s => s,
637              v => v.Sum());
638            TraceEvaluation(currentInstr, cur);
639            return cur;
640          }
641        case OpCodes.Mean: {
642            var cur = Evaluate(dataset, ref row, state, traceDict);
643            cur = AggregateApply(cur,
644              s => s,
645              v => Statistics.Mean(v));
646            TraceEvaluation(currentInstr, cur);
647            return cur;
648          }
649        case OpCodes.StandardDeviation: {
650            var cur = Evaluate(dataset, ref row, state, traceDict);
651            cur = AggregateApply(cur,
652              s => 0,
653              v => Statistics.PopulationStandardDeviation(v));
654            TraceEvaluation(currentInstr, cur);
655            return cur;
656          }
657        case OpCodes.Length: {
658            var cur = Evaluate(dataset, ref row, state, traceDict);
659            cur = AggregateApply(cur,
660              s => 1,
661              v => v.Count);
662            TraceEvaluation(currentInstr, cur);
663            return cur;
664          }
665        case OpCodes.Min: {
666            var cur = Evaluate(dataset, ref row, state, traceDict);
667            cur = AggregateApply(cur,
668              s => s,
669              v => Statistics.Minimum(v));
670            TraceEvaluation(currentInstr, cur);
671            return cur;
672          }
673        case OpCodes.Max: {
674            var cur = Evaluate(dataset, ref row, state, traceDict);
675            cur = AggregateApply(cur,
676              s => s,
677              v => Statistics.Maximum(v));
678            TraceEvaluation(currentInstr, cur);
679            return cur;
680          }
681        case OpCodes.Variance: {
682            var cur = Evaluate(dataset, ref row, state, traceDict);
683            cur = AggregateApply(cur,
684              s => 0,
685              v => Statistics.PopulationVariance(v));
686            TraceEvaluation(currentInstr, cur);
687            return cur;
688          }
689        case OpCodes.Skewness: {
690            var cur = Evaluate(dataset, ref row, state, traceDict);
691            cur = AggregateApply(cur,
692              s => double.NaN,
693              v => Statistics.PopulationSkewness(v));
694            TraceEvaluation(currentInstr, cur);
695            return cur;
696          }
697        case OpCodes.Kurtosis: {
698            var cur = Evaluate(dataset, ref row, state, traceDict);
699            cur = AggregateApply(cur,
700              s => double.NaN,
701              v => Statistics.PopulationKurtosis(v));
702            TraceEvaluation(currentInstr, cur);
703            return cur;
704          }
705        case OpCodes.EuclideanDistance: {
706            var x1 = Evaluate(dataset, ref row, state, traceDict);
707            var x2 = Evaluate(dataset, ref row, state, traceDict);
708            var cur = AggregateMultipleApply(x1, x2,
709              (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 0.0),
710              (s1, s2) => s1 - s2,
711              (s1, v2) => Math.Sqrt((s1 - v2).PointwisePower(2).Sum()),
712              (v1, s2) => Math.Sqrt((v1 - s2).PointwisePower(2).Sum()),
713              (v1, v2) => Math.Sqrt((v1 - v2).PointwisePower(2).Sum()));
714            TraceEvaluation(currentInstr, cur);
715            return cur;
716          }
717        case OpCodes.Covariance: {
718            var x1 = Evaluate(dataset, ref row, state, traceDict);
719            var x2 = Evaluate(dataset, ref row, state, traceDict);
720            var cur = AggregateMultipleApply(x1, x2,
721              (lhs, rhs) => ApplyVectorLengthStrategy(DifferentVectorLengthStrategy, lhs, rhs, 0.0),
722              (s1, s2) => 0,
723              (s1, v2) => 0,
724              (v1, s2) => 0,
725              (v1, v2) => Statistics.PopulationCovariance(v1, v2));
726            TraceEvaluation(currentInstr, cur);
727            return cur;
728          }
729        case OpCodes.SubVector: {
730            DoubleVector SubVector(DoubleVector v , double start, double end) {
731              int startIdx = ToVectorIdx(start, v.Count), endIdx = ToVectorIdx(end, v.Count);
732              return RoundTripSubVector(v, startIdx, endIdx);
733            }
734            var cur = Evaluate(dataset, ref row, state, traceDict);
735            TraceEvaluation(currentInstr, cur);
736            return FunctionApply(cur,
737              s => s,
738              v => {
739                var node = (IWindowedSymbolTreeNode)currentInstr.dynamicNode;
740                return SubVector(v, node.Offset, node.Length);
741              });
742          }
743        case OpCodes.SubVectorSubtree: {
744            DoubleVector SubVector(DoubleVector v, double start, double end) {
745              int startIdx = ToVectorIdx(start, v.Count), endIdx = ToVectorIdx(end, v.Count);
746              return RoundTripSubVector(v, startIdx, endIdx);
747            }
748            var cur = Evaluate(dataset, ref row, state, traceDict);
749            var offset = Evaluate(dataset, ref row, state, traceDict);
750            var length = Evaluate(dataset, ref row, state, traceDict);
751            TraceEvaluation(currentInstr, cur);
752            return FunctionApply(cur,
753              s => s,
754              v => SubVector(v, offset.Scalar, length.Scalar)
755            );
756          }
757        case OpCodes.Variable: {
758            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
759            var variableTreeNode = (VariableTreeNode)currentInstr.dynamicNode;
760            if (currentInstr.data is IList<double> doubleList) {
761              var cur = new EvaluationResult(doubleList[row] * variableTreeNode.Weight);
762              TraceEvaluation(currentInstr, cur);
763              return cur;
764            }
765            if (currentInstr.data is IList<DoubleVector> doubleVectorList) {
766              var cur = new EvaluationResult(doubleVectorList[row] * variableTreeNode.Weight);
767              TraceEvaluation(currentInstr, cur);
768              return cur;
769            }
770            throw new NotSupportedException($"Unsupported type of variable: {currentInstr.data.GetType().GetPrettyName()}");
771          }
772        case OpCodes.BinaryFactorVariable: {
773            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
774            var factorVarTreeNode = currentInstr.dynamicNode as BinaryFactorVariableTreeNode;
775            var cur = new EvaluationResult(((IList<string>)currentInstr.data)[row] == factorVarTreeNode.VariableValue ? factorVarTreeNode.Weight : 0);
776            TraceEvaluation(currentInstr, cur);
777            return cur;
778          }
779        case OpCodes.FactorVariable: {
780            if (row < 0 || row >= dataset.Rows) return EvaluationResult.NaN;
781            var factorVarTreeNode = currentInstr.dynamicNode as FactorVariableTreeNode;
782            var cur = new EvaluationResult(factorVarTreeNode.GetValue(((IList<string>)currentInstr.data)[row]));
783            TraceEvaluation(currentInstr, cur);
784            return cur;
785          }
786        case OpCodes.Constant: {
787            var constTreeNode = (ConstantTreeNode)currentInstr.dynamicNode;
788            var cur = new EvaluationResult(constTreeNode.Value);
789            TraceEvaluation(currentInstr, cur);
790            return cur;
791          }
792
793        #region Time Series Symbols
794        case OpCodes.Median: {
795            var cur = Evaluate(dataset, ref row, state, traceDict);
796            cur = AggregateApply(cur,
797              s => s,
798              v => Statistics.Median(v));
799            TraceEvaluation(currentInstr, cur);
800            return cur;
801          }
802        case OpCodes.Quantile: {
803            var cur = Evaluate(dataset, ref row, state, traceDict);
804            var q = Evaluate(dataset, ref row, state, traceDict);
805            cur = AggregateApply(cur,
806              s => s,
807              v => Statistics.Quantile(v, q.Scalar));
808            TraceEvaluation(currentInstr, cur);
809            return cur;
810          }
811
812        case OpCodes.AbsoluteEnergy: {
813            var cur = Evaluate(dataset, ref row, state, traceDict);
814            cur = AggregateApply(cur,
815              s => s * s,
816              v => v.PointwisePower(2.0).Sum());
817            TraceEvaluation(currentInstr, cur);
818            return cur;
819          }
820
821        case OpCodes.BinnedEntropy: {
822            var cur = Evaluate(dataset, ref row, state, traceDict);
823            var m = Evaluate(dataset, ref row, state, traceDict);
824            cur = AggregateApply(cur,
825              s => 0,
826              v => {
827                int bins = Math.Max((int)Math.Round(m.Scalar), 1);
828                double minValue = v.Minimum();
829                double maxValue = v.Maximum();
830                double intervalWidth = (maxValue - minValue) / bins;
831                int totalValues = v.Count;
832                double sum = 0;
833                for (int i = 0; i < Math.Max(bins, v.Count); i++) {
834                  double binMin = minValue * i;
835                  double binMax = binMin + intervalWidth;
836                  double countBin = v.Map(e => (e > binMin && e < binMax) ? 1.0 : 0.0).Sum();
837                  double percBin = countBin / totalValues;
838                  sum += percBin * Math.Log(percBin);
839                }
840
841                return sum;
842              });
843            TraceEvaluation(currentInstr, cur);
844            return cur;
845          }
846        case OpCodes.HasLargeStandardDeviation: {
847            var cur = Evaluate(dataset, ref row, state, traceDict);
848            cur = AggregateApply(cur,
849              s => 0,
850              v => Statistics.PopulationStandardDeviation(v) > (Statistics.Maximum(v) - Statistics.Minimum(v)) / 2 ? 1.0 : 0.0);
851            TraceEvaluation(currentInstr, cur);
852            return cur;
853          }
854        case OpCodes.HasVarianceLargerThanStd: {
855            var cur = Evaluate(dataset, ref row, state, traceDict);
856            cur = AggregateApply(cur,
857              s => 0,
858              v => Statistics.PopulationVariance(v) > Statistics.StandardDeviation(v) ? 1.0 : 0.0);
859            TraceEvaluation(currentInstr, cur);
860            return cur;
861          }
862        case OpCodes.IsSymmetricLooking: {
863            var cur = Evaluate(dataset, ref row, state, traceDict);
864            cur = AggregateApply(cur,
865              s => 0,
866              v => Math.Abs(Statistics.Mean(v) - Statistics.Median(v)) < (Statistics.Maximum(v) - Statistics.Minimum(v)) / 2 ? 1.0 : 0.0);
867            TraceEvaluation(currentInstr, cur);
868            return cur;
869          }
870        case OpCodes.NumberDataPointsAboveMean: {
871            var cur = Evaluate(dataset, ref row, state, traceDict);
872            cur = AggregateApply(cur,
873              s => 0,
874              v => {
875                double mean = Statistics.Mean(v);
876                return v.Map(e => e > mean ? 1.0 : 0.0).Sum();
877              });
878            TraceEvaluation(currentInstr, cur);
879            return cur;
880          }
881        case OpCodes.NumberDataPointsAboveMedian: {
882            var cur = Evaluate(dataset, ref row, state, traceDict);
883            cur = AggregateApply(cur,
884              s => 0,
885              v => {
886                double median = Statistics.Median(v);
887                return v.Map(e => e > median ? 1.0 : 0.0).Sum();
888              });
889            TraceEvaluation(currentInstr, cur);
890            return cur;
891          }
892        case OpCodes.NumberDataPointsBelowMean: {
893            var cur = Evaluate(dataset, ref row, state, traceDict);
894            cur = AggregateApply(cur,
895              s => 0,
896              v => {
897                double mean = Statistics.Mean(v);
898                return v.Map(e => e < mean ? 1.0 : 0.0).Sum();
899              });
900            TraceEvaluation(currentInstr, cur);
901            return cur;
902          }
903        case OpCodes.NumberDataPointsBelowMedian: {
904            var cur = Evaluate(dataset, ref row, state, traceDict);
905            cur = AggregateApply(cur,
906              s => 0,
907              v => {
908                double median = Statistics.Median(v);
909                return v.Map(e => e < median ? 1.0 : 0.0).Sum();
910              });
911            TraceEvaluation(currentInstr, cur);
912            return cur;
913          }
914
915        case OpCodes.ArimaModelCoefficients: {
916            var cur = Evaluate(dataset, ref row, state, traceDict);
917            var i = Evaluate(dataset, ref row, state, traceDict);
918            var k = Evaluate(dataset, ref row, state, traceDict);
919            cur = AggregateApply(cur,
920              s => 0,
921              v => throw new NotImplementedException(""));
922            TraceEvaluation(currentInstr, cur);
923            return cur;
924          }
925        case OpCodes.ContinuousWaveletTransformationCoefficients: {
926            var cur = Evaluate(dataset, ref row, state, traceDict);
927            var a = Evaluate(dataset, ref row, state, traceDict);
928            var b = Evaluate(dataset, ref row, state, traceDict);
929            cur = AggregateApply(cur,
930              s => 0,
931              v => throw new NotImplementedException(""));
932            TraceEvaluation(currentInstr, cur);
933            return cur;
934          }
935        case OpCodes.FastFourierTransformationCoefficient: {
936            var cur = Evaluate(dataset, ref row, state, traceDict);
937            var k = Evaluate(dataset, ref row, state, traceDict);
938            cur = AggregateApply(cur,
939              s => 0,
940              v => throw new NotImplementedException(""));
941            TraceEvaluation(currentInstr, cur);
942            return cur;
943          }
944        case OpCodes.FirstIndexMax: {
945            var cur = Evaluate(dataset, ref row, state, traceDict);
946            cur = AggregateApply(cur,
947              s => 0,
948              v => (double)v.MaximumIndex() / v.Count);
949            TraceEvaluation(currentInstr, cur);
950            return cur;
951          }
952        case OpCodes.FirstIndexMin: {
953            var cur = Evaluate(dataset, ref row, state, traceDict);
954            cur = AggregateApply(cur,
955              s => 0,
956              v => (double)v.MinimumIndex() / v.Count);
957            TraceEvaluation(currentInstr, cur);
958            return cur;
959          }
960        case OpCodes.LastIndexMax: {
961            var cur = Evaluate(dataset, ref row, state, traceDict);
962            cur = AggregateApply(cur,
963              s => 0,
964              v => (double)(v.Count - DoubleVector.Build.DenseOfEnumerable(v.Reverse()).MaximumIndex()) / v.Count);
965
966            TraceEvaluation(currentInstr, cur);
967            return cur;
968          }
969        case OpCodes.LastIndexMin: {
970            var cur = Evaluate(dataset, ref row, state, traceDict);
971            cur = AggregateApply(cur,
972              s => 0,
973              v => (double)(v.Count - DoubleVector.Build.DenseOfEnumerable(v.Reverse()).MinimumIndex()) / v.Count);
974            TraceEvaluation(currentInstr, cur);
975            return cur;
976          }
977        case OpCodes.LongestStrikeAboveMean: {
978            var cur = Evaluate(dataset, ref row, state, traceDict);
979            cur = AggregateApply(cur,
980              s => 0,
981              v => LongestStrikeAbove(v, Statistics.Mean(v)));
982            TraceEvaluation(currentInstr, cur);
983            return cur;
984          }
985        case OpCodes.LongestStrikeAboveMedian: {
986            var cur = Evaluate(dataset, ref row, state, traceDict);
987            cur = AggregateApply(cur,
988              s => 0,
989              v => LongestStrikeAbove(v, Statistics.Median(v)));
990            TraceEvaluation(currentInstr, cur);
991            return cur;
992          }
993        case OpCodes.LongestStrikeBelowMean: {
994            var cur = Evaluate(dataset, ref row, state, traceDict);
995            cur = AggregateApply(cur,
996              s => 0,
997              v => LongestStrikeBelow(v, Statistics.Mean(v)));
998            TraceEvaluation(currentInstr, cur);
999            return cur;
1000          }
1001        case OpCodes.LongestStrikeBelowMedian: {
1002            var cur = Evaluate(dataset, ref row, state, traceDict);
1003            cur = AggregateApply(cur,
1004              s => 0,
1005              v => LongestStrikeBelow(v, Statistics.Median(v)));
1006            TraceEvaluation(currentInstr, cur);
1007            return cur;
1008          }
1009        case OpCodes.LongestStrikePositive: {
1010            var cur = Evaluate(dataset, ref row, state, traceDict);
1011            cur = AggregateApply(cur,
1012              s => 0,
1013              v => LongestStrikeAbove(v, 0));
1014            TraceEvaluation(currentInstr, cur);
1015            return cur;
1016          }
1017        case OpCodes.LongestStrikeNegative: {
1018            var cur = Evaluate(dataset, ref row, state, traceDict);
1019            cur = AggregateApply(cur,
1020              s => 0,
1021              v => LongestStrikeAbove(v, 0));
1022            TraceEvaluation(currentInstr, cur);
1023            return cur;
1024          }
1025        case OpCodes.LongestStrikeZero: {
1026            var cur = Evaluate(dataset, ref row, state, traceDict);
1027            cur = AggregateApply(cur,
1028              s => 0,
1029              v => LongestStrikeEqual(v, 0));
1030            TraceEvaluation(currentInstr, cur);
1031            return cur;
1032          }
1033        case OpCodes.MeanAbsoluteChange: {
1034            var cur = Evaluate(dataset, ref row, state, traceDict);
1035            cur = AggregateApply(cur,
1036              s => 0,
1037              v => {
1038                double sum = 0.0;
1039                for (int i = 0; i < v.Count - 1; i++) {
1040                  sum += Math.Abs(v[i + 1] - v[i]);
1041                }
1042
1043                return sum / v.Count;
1044              });
1045            TraceEvaluation(currentInstr, cur);
1046            return cur;
1047          }
1048        case OpCodes.MeanAbsoluteChangeQuantiles: {
1049            var cur = Evaluate(dataset, ref row, state, traceDict);
1050            var ql = Evaluate(dataset, ref row, state, traceDict);
1051            var qu = Evaluate(dataset, ref row, state, traceDict);
1052            cur = AggregateApply(cur,
1053              s => 0,
1054              v => {
1055                var lowerBound = Statistics.Quantile(v, ql.Scalar);
1056                var upperBound = Statistics.Quantile(v, qu.Scalar);
1057                var inBounds = v.Select(e => e > lowerBound && e < upperBound).ToList();
1058                double sum = 0.0;
1059                int count = 0;
1060                for (int i = 0; i < v.Count - 1; i++) {
1061                  if (inBounds[i] && inBounds[i + 1]) {
1062                    sum += Math.Abs(v[i + 1] - v[i]);
1063                    count++;
1064                  }
1065                }
1066
1067                return sum / count;
1068              });
1069            TraceEvaluation(currentInstr, cur);
1070            return cur;
1071          }
1072        case OpCodes.MeanAutocorrelation: {
1073            var cur = Evaluate(dataset, ref row, state, traceDict);
1074            cur = AggregateApply(cur,
1075              s => 0,
1076              v => {
1077                double sum = 0.0;
1078                double mean = Statistics.Mean(v);
1079                for (int l = 0; l < v.Count; l++) {
1080                  for (int i = 0; i < v.Count - l; i++) {
1081                    sum += (v[i] - mean) * (v[i + l] - mean);
1082                  }
1083                }
1084
1085                return sum / (v.Count - 1) / Statistics.PopulationVariance(v);
1086              });
1087            TraceEvaluation(currentInstr, cur);
1088            return cur;
1089          }
1090        case OpCodes.LaggedAutocorrelation: {
1091            var cur = Evaluate(dataset, ref row, state, traceDict);
1092            var lVal = Evaluate(dataset, ref row, state, traceDict);
1093            cur = AggregateApply(cur,
1094              s => 0,
1095              v => {
1096                double sum = 0.0;
1097                int l = Math.Max((int)Math.Round(lVal.Scalar), 0);
1098                double mean = Statistics.Mean(v);
1099                for (int i = 0; i < v.Count - l; i++) {
1100                  sum += (v[i] - mean) * (v[i + l] - mean);
1101                }
1102
1103                return sum / Statistics.PopulationVariance(v);
1104              });
1105            TraceEvaluation(currentInstr, cur);
1106            return cur;
1107          }
1108        case OpCodes.MeanSecondDerivateCentral: {
1109            var cur = Evaluate(dataset, ref row, state, traceDict);
1110            cur = AggregateApply(cur,
1111              s => 0,
1112              v => {
1113                double sum = 0.0;
1114                for (int i = 1; i < v.Count - 1; i++) {
1115                  sum += (v[i - 1] - 2 * v[i] + v[i + 1]) / 2;
1116                }
1117
1118                return sum / (v.Count - 2);
1119              });
1120            TraceEvaluation(currentInstr, cur);
1121            return cur;
1122          }
1123        case OpCodes.NumberPeaksOfSize: {
1124            var cur = Evaluate(dataset, ref row, state, traceDict);
1125            var l = Evaluate(dataset, ref row, state, traceDict);
1126            cur = AggregateApply(cur,
1127              s => 0,
1128              v => CountNumberOfPeaks(v, l.Scalar));
1129            TraceEvaluation(currentInstr, cur);
1130            return cur;
1131          }
1132        case OpCodes.LargeNumberOfPeaks: {
1133            var cur = Evaluate(dataset, ref row, state, traceDict);
1134            var l = Evaluate(dataset, ref row, state, traceDict);
1135            var m = Evaluate(dataset, ref row, state, traceDict);
1136            cur = AggregateApply(cur,
1137              s => 0,
1138              v => CountNumberOfPeaks(v, l.Scalar) > m.Scalar ? 1.0 : 0.0);
1139            TraceEvaluation(currentInstr, cur);
1140            return cur;
1141          }
1142        case OpCodes.TimeReversalAsymmetryStatistic: {
1143            var cur = Evaluate(dataset, ref row, state, traceDict);
1144            var l = Evaluate(dataset, ref row, state, traceDict);
1145            cur = AggregateApply(cur,
1146              s => 0,
1147              v => {
1148                int lag = Math.Max((int)Math.Round(l.Scalar), 0);
1149                double sum = 0.0;
1150                for (int i = 0; i < v.Count - 2 * lag; i++) {
1151                  sum += Math.Pow(v[i + 2 * lag], 2) * v[i + lag] - v[i + lag] * Math.Pow(v[i], 2);
1152                }
1153
1154                return sum / (v.Count - 2 * lag);
1155              });
1156            TraceEvaluation(currentInstr, cur);
1157            return cur;
1158          }
1159        #endregion
1160
1161        default:
1162          throw new NotSupportedException($"Unsupported OpCode: {currentInstr.opCode}");
1163      }
1164    }
1165
1166    public static int ToVectorIdx(double relPosition, int size) {
1167      int idx = (int)Math.Round(relPosition * (size - 1));
1168      return (idx % size + size) % size; // positive mod
1169    }
1170
1171    public static IEnumerable<Tuple<int, int>>/*start, count*/ GetVectorSlices(int startIdx, int endIdx, int size) {
1172      if (startIdx <= endIdx) {
1173        yield return Tuple.Create(startIdx, endIdx - startIdx + 1); // incl end
1174      } else {
1175        yield return Tuple.Create(startIdx, size - startIdx); // startIdx to end of vector
1176        yield return Tuple.Create(0, endIdx);                 // start to endIdx of vector
1177      }
1178    }
1179   
1180    public static DoubleVector RoundTripSubVector(DoubleVector v, int startIdx, int endIdx) {
1181      var slices = GetVectorSlices(startIdx, endIdx, v.Count).ToList();
1182      var totalSize = slices.Sum(s => s.Item2);
1183      var resultVector = DoubleVector.Build.Dense(size: totalSize);
1184
1185      var curIdx = 0;
1186      foreach (var slice in slices) {
1187        int start = slice.Item1, count = slice.Item2;
1188        v.CopySubVectorTo(resultVector, sourceIndex: start, targetIndex: curIdx, count: count);
1189        curIdx += count;
1190      }
1191      return resultVector;
1192    }
1193
1194    private static int LongestStrikeAbove(DoubleVector v, double threshold) {
1195      int longestStrike = 0, currentStrike = 0;
1196      for (int i = 0; i < v.Count; i++) {
1197        if (v[i] > threshold) {
1198          currentStrike++;
1199          longestStrike = Math.Max(longestStrike, currentStrike);
1200        } else
1201          currentStrike = 0;
1202      }
1203      return longestStrike;
1204    }
1205    private static int LongestStrikeBelow(DoubleVector v, double threshold) {
1206      int longestStrike = 0, currentStrike = 0;
1207      for (int i = 0; i < v.Count; i++) {
1208        if (v[i] < threshold) {
1209          currentStrike++;
1210          longestStrike = Math.Max(longestStrike, currentStrike);
1211        } else
1212          currentStrike = 0;
1213      }
1214      return longestStrike;
1215    }
1216
1217    private static int LongestStrikeEqual(DoubleVector v, double value, double epsilon = double.Epsilon) {
1218      int longestStrike = 0, currentStrike = 0;
1219      for (int i = 0; i < v.Count; i++) {
1220        if (v[i].IsAlmost(epsilon)) {
1221          currentStrike++;
1222          longestStrike = Math.Max(longestStrike, currentStrike);
1223        } else
1224          currentStrike = 0;
1225      }
1226      return longestStrike;
1227    }
1228    private static int CountNumberOfPeaks(DoubleVector v, double heightDifference) {
1229      int count = 0;
1230      for (int i = 0; i < v.Count; i++) {
1231        bool largerThanPrev = i == 0 || v[i] > v[i - 1] + heightDifference;
1232        bool largerThanNext = i == v.Count - 1 || v[i] > v[i + 1] + heightDifference;
1233        if (largerThanPrev && largerThanNext)
1234          count++;
1235      }
1236      return count;
1237    }
1238  }
1239}
Note: See TracBrowser for help on using the repository browser.