Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
03/29/19 15:01:47 (6 years ago)
Author:
gkronber
Message:

#2994: added a unit test and made some minor improvements to interpreters

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs

    r16695 r16727  
    7272              break;
    7373            }
     74          case OpCodes.SquareRoot: {
     75              instr.value.AssignIntRoot(code[c].value, 2);
     76              break;
     77            }
     78          case OpCodes.Cube: {
     79              instr.value.AssignIntPower(code[c].value, 3);
     80              break;
     81            }
     82          case OpCodes.CubeRoot: {
     83              instr.value.AssignIntRoot(code[c].value, 3);
     84              break;
     85            }
    7486          case OpCodes.Exp: {
    7587              instr.value.AssignExp(code[c].value);
    7688              break;
    7789            }
    78 
    7990          case OpCodes.Log: {
    8091              instr.value.AssignLog(code[c].value);
     92              break;
     93            }
     94          case OpCodes.Sin: {
     95              instr.value.AssignSin(code[c].value);
     96              break;
     97            }
     98          case OpCodes.Cos: {
     99              instr.value.AssignCos(code[c].value);
     100              break;
     101            }
     102          case OpCodes.Absolute: {
     103              instr.value.AssignAbs(code[c].value);
     104              break;
     105            }
     106          case OpCodes.AnalyticQuotient: {
     107              instr.value.Assign(code[c].value);
     108              for (int j = 1; j < n; ++j) {
     109                var t = instr.value.One;
     110                t.Add(code[c + j].value.Clone().IntPower(2));
     111                instr.value.Div(t.IntRoot(2));
     112              }
    81113              break;
    82114            }
    … …  
    221253    }
    222254
    223     public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, out double[] fi, out double[,] jac) {
     255    /// <summary>
     256    ///
     257    /// </summary>
     258    /// <param name="tree"></param>
     259    /// <param name="dataset"></param>
     260    /// <param name="rows"></param>
     261    /// <param name="parameterNodes"></param>
     262    /// <param name="fi">Function output. Must be allocated by the caller.</param>
     263    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
     264    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
    224265      if (cachedData == null || this.dataset != dataset) {
    225266        InitCache(dataset);
    … …  
    235276      var roundedTotal = rows.Length - remainingRows;
    236277
    237       fi = new double[rows.Length];
    238       jac = new double[rows.Length, nParams];
    239 
    240278      this.rows = rows;
    241279
    … …  
    247285        var g = code[0].value.Gradient;
    248286        for (int j = 0; j < nParams; ++j) {
    249           g.Elements[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
     287          if(g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
     288            v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
     289          } else {
     290            for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
     291          }
    250292        }
    251293      }
    … …  
    257299        var g = code[0].value.Gradient;
    258300        for (int j = 0; j < nParams; ++j)
    259           g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
     301          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
     302            v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
     303          } else {
     304            for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
     305          }
    260306      }
    261307    }
    … …  
    315361        var g = a.value.Gradient.Elements[paramIdx];
    316362        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
    317           g[i] = data[rows[i]];
     363          g[i - rowIndex] = data[rows[i]];
    318364        }
    319365      }
    … …  
    377423  public interface IAlgebraicType<T> {
    378424    T Zero { get; }
     425    T One { get; }
    379426
    380427    T AssignAbs(T a); // set this to assign abs(a)
    … …  
    429476    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    430477    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
     478    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     479    public AlgebraicDouble One => new AlgebraicDouble(1.0);
    431480    public AlgebraicDouble() { }
    432481    public AlgebraicDouble(double value) { this.Value = value; }
    … …  
    441490    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
    442491    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
    443     public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
     492    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = a.Value <= 0?double.NegativeInfinity : Math.Log(a.Value); return this; } // alternative definiton of log to prevent NaN
    444493    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
    445494    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
    446495    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
    447496    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
    448     public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; }
     497    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value)? double.NaN : Math.Sign(a.Value); return this; }
    449498    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
    450499
    … …  
    473522    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    474523    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
     524    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     525    public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)
    475526    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
    476527    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
    … …  
    496547    }
    497548
    498     public void AssignConstant(double constantValue) {
     549    public AlgebraicDoubleVector AssignConstant(double constantValue) {
    499550      for (int i = 0; i < arr.Length; ++i) {
    500551        arr[i] = constantValue;
    501552      }
     553      return this;
    502554    }
    503555
    … …  
    522574  }
    523575
     576
     577  /*
    524578  // vectors of algebraic types
    525   public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T> {
     579  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
    526580    private T[] elems;
    527581
    … …  
    561615    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
    562616
    563     public AlgebraicVector<T> Concat(AlgebraicVector<T> other) {
    564       var oldLen = Length;
    565       Array.Resize(ref this.elems, oldLen + other.Length);
    566       for (int i = oldLen; i < Length; i++) {
    567         elems[i] = other.elems[i - oldLen].Clone();
    568       }
    569       return this;
    570     }
    571617
    572618    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    573619    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
     620    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     621    public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }
    574622    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
    575623    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
    … …  
    591639  }
    592640
     641  */
     642
    593643
    594644  /// <summary>
    … …  
    625675    }
    626676
    627 
    628 
    629     // combined elements from both vectors
    630     private void UnionAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    631       // elements from a
    632       foreach (var kvp in a.elems) {
    633         // this = f(a, this)
    634         if (elems.TryGetValue(kvp.Key, out T value))
    635           mapAssign(kvp.Value, value);
    636         else {
    637           // this = f(a, 0)
    638           var newValue = kvp.Value.Zero;
    639           elems.Add(kvp.Key, newValue);
    640           mapAssign(kvp.Value, newValue);
    641         }
    642       }
    643       // elements from this (without a)
    644       foreach (var kvp in elems) {
    645         if (a.elems.ContainsKey(kvp.Key)) continue; // already processed above
    646                                                     // this = f(0, this)
    647         mapAssign(kvp.Value.Zero, kvp.Value);
    648       }
    649     }
    650 
    651     // keep only elements in both vectors
    652     private void IntersectAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    653       List<K> keysToRemove = new List<K>();
    654       foreach (var kvp in elems) {
    655         if (a.elems.TryGetValue(kvp.Key, out T value))
    656           mapAssign(value, kvp.Value);
    657         else
    658           keysToRemove.Add(kvp.Key);
    659       }
    660       foreach (var o in keysToRemove) elems.Remove(o); // -> zero
    661     }
    662 
    663677    // keep only elements from a
    664678    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    … …  
    680694    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    681695    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
     696    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     697    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
    682698
    683699    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
    … …  
    685701
    686702    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
    687     public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Add(src)); return this; }
    688     public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Sub(src)); return this; }
    689     public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) { IntersectAssign(a, (src, dest) => dest.Mul(src)); return this; }
    690     public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Div(src)); return this; }
    691703    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
    692704    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
    … …  
    699711    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
    700712    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
     713    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
     714      foreach (var kvp in a.elems) {
     715        if (elems.TryGetValue(kvp.Key, out T value))
     716          value.Add(kvp.Value);
     717        else
     718          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
     719      }
     720      return this;
     721    }
     722
     723    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
     724      foreach (var kvp in a.elems) {
     725        if (elems.TryGetValue(kvp.Key, out T value))
     726          value.Sub(kvp.Value);
     727        else
     728          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
     729      }
     730      return this;
     731    }
     732
     733    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
     734      var keys = elems.Keys.ToArray();
     735      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
     736      foreach (var kvp in a.elems) {
     737        if (elems.TryGetValue(kvp.Key, out T value))
     738          value.Mul(kvp.Value); // this * a
     739      }
     740      return this;
     741    }
     742
     743    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
     744      return Mul(a.Clone().Inv());
     745    }
    701746
    702747    public AlgebraicSparseVector<K, T> Clone() {
    … …  
    734779    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    735780    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
     781    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
    736782    public AlgebraicInterval Add(AlgebraicInterval a) {
    737783      low.Add(a.low);
    … …  
    799845          if (a.Contains(0.0)) {
    800846            low = new MultivariateDual<AlgebraicDouble>(0.0);
    801             high = Algebraic.Max(low.Clone().IntPower(p), high.Clone().IntPower(p));
     847            high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
    802848          } else {
    803             var lowPower = low.Clone().IntPower(p);
    804             var highPower = high.Clone().IntPower(p);
     849            var lowPower = a.low.IntPower(p);
     850            var highPower = a.high.IntPower(p);
    805851            low = Algebraic.Min(lowPower, highPower);
    806852            high = Algebraic.Max(lowPower, highPower);
    … …  
    808854        } else {
    809855          // p is uneven
    810           var lowPower = low.Clone().IntPower(p);
    811           var highPower = high.Clone().IntPower(p);
     856          var lowPower = a.low.IntPower(p);
     857          var highPower = a.high.IntPower(p);
    812858          low = Algebraic.Min(lowPower, highPower);
    813859          high = Algebraic.Max(lowPower, highPower);
    … …  
    951997    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    952998    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
     999    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     1000    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
    9531001
    9541002    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
    … …  
    10341082
    10351083    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
     1084    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
    10361085
    10371086    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
Note: See TracChangeset for help on using the changeset viewer.