Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
03/13/19 23:55:38 (6 years ago)
Author:
gkronber
Message:

#2994: worked on auto diff for intervals and vectors

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r16674 r16682  
    192192  }
    193193
    194   public class VectorAutoDiffEvaluator : Interpreter<VectorDual<DoubleVector>> {
     194  public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
    195195    private const int BATCHSIZE = 128;
    196196    [ThreadStatic]
    … …  
    241241
    242242        // TRANSPOSE into JAC
    243         // here we assume that we usually calculate the gradient over _all_ parameters (otherwise it would be better to propagate only relevant vectors in evaluation)
    244         for (int j = 0; j < nParams; ++j)
    245           code[0].value.Gradient[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); XXX CONTINUE HERE
     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        }
    246247      }
    247248
    … …  
    249250        Evaluate(code);
    250251        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
     252
     253        var g = code[0].value.Gradient;
    251254        for (int j = 0; j < nParams; ++j)
    252           code[0].value.Gradient[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
     255          g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
    253256      }
    254257    }
    … …  
    256259    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
    257260      var zero = new DoubleVector(BATCHSIZE);
    258       //       var g = new Vector<DoubleVector>(node2paramIdx.Count);
    259       //       for (int i = 0; i < BATCHSIZE; ++i) g[i] = new DoubleVector(BATCHSIZE);
    260       instruction.value = new VectorDual<DoubleVector>(zero, null);
     261      instruction.value = new MultivariateDual<DoubleVector>(zero);
    261262    }
    262263
    263264    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
    264274      instruction.dblVal = constant.Value;
    265 
    266       var g_arr = new double[BATCHSIZE];
    267       if (node2paramIdx.ContainsKey(constant)) {
    268         for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
    269       }
    270       var g = new DoubleVector(g_arr);
    271 
    272       instruction.value = new VectorDual<DoubleVector>(new DoubleVector(BATCHSIZE), new Vector<DoubleVector>(new[] { g })); // only a single column for the gradient
    273275      instruction.value.Value.AssignConstant(instruction.dblVal);
    274276    }
    … …  
    286288      if (node2paramIdx.ContainsKey(variable)) {
    287289        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);
    288296      }
    289297
    290298      instruction.dblVal = variable.Weight;
    291299      instruction.data = new object[] { data, paramIdx };
    292 
    293       var f = new DoubleVector(BATCHSIZE);
    294       var g = new DoubleVector(BATCHSIZE);
    295       instruction.value = new VectorDual<DoubleVector>(f, new Vector<DoubleVector>(new[] { g }));
    296300    }
    297301
    … …  
    305309      if (paramIdx >= 0) {
    306310        // update gradient with variable values
    307         var g = a.value.Gradient[0];
     311        var g = a.value.Gradient.Elements[paramIdx];
    308312        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
    309313          g[i] = data[rows[i]];
    … …  
    316320  public class IntervalEvaluator : Interpreter<AlgebraicInterval> {
    317321    [ThreadStatic]
    318     private Dictionary<string, AlgebraicInterval> intervals;
    319 
    320     public AlgebraicInterval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, AlgebraicInterval> intervals) {
     322    private Dictionary<string, Interval> intervals;
     323
     324    public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals) {
    321325      this.intervals = intervals;
    322326      var code = Compile(tree);
    323327      Evaluate(code);
    324       return code[0].value;
    325     }
    326 
     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    }
    327345
    328346    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
    … …  
    333351    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
    334352      instruction.dblVal = constant.Value;
    335       instruction.value = new AlgebraicInterval(constant.Value, 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        );
    336357    }
    337358
    338359    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
    339360      instruction.dblVal = variable.Weight;
    340       instruction.value = intervals[variable.VariableName];
     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        );
    341365    }
    342366
    … …  
    347371
    348372  public interface IAlgebraicType<T> {
     373    T AssignAbs(T a); // set this to assign abs(a)
    349374    T Assign(T a); // assign this to same value as a (copy!)
    350375    T AssignNeg(T a); // set this to negative(a)
    … …  
    361386    T AssignIntPower(T a, int p);
    362387    T AssignIntRoot(T a, int r);
     388    T AssignSgn(T a); // set this to sign(a)
    363389    T Clone();
    364390  }
    365391
     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  }
    366442
    367443  // a simple vector as an algebraic type
    … …  
    374450      arr = new double[length];
    375451    }
     452
     453    public DoubleVector() { }
    376454
    377455    /// <summary>
    … …  
    398476    }
    399477
    400 
    401     public DoubleVector Add(DoubleVector a) {
    402       for (int i = 0; i < arr.Length; ++i) {
    403         arr[i] += a.arr[i];
    404       }
    405       return this;
    406     }
    407 
    408478    public void AssignConstant(double constantValue) {
    409479      for (int i = 0; i < arr.Length; ++i) {
    … …  
    412482    }
    413483
    414     public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
    415     public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
    416     public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
    417     public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
     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; }
    418488    public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }
    419489    public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }
    420     public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
    421     public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
     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; }
    422493    public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
    423     public DoubleVector AssignNeg(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; }
    424495    public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
    425     public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } 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; }
    426497    public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
    427 
     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; }
    428500    public DoubleVector Clone() {
    429501      var v = new DoubleVector(this.arr.Length);
    … …  
    437509
    438510  }
    439 
    440511
    441512  // vectors of algebraic types
    … …  
    506577    public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
    507578    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    }
    508713  }
    509714
    510715
    511716  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
    512     private double low;
    513     private double high;
    514 
    515     public double LowerBound => low;
    516     public double UpperBound => high;
     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();
    517722
    518723    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
    519724
     725    public AlgebraicInterval(MultivariateDual<Double> low, MultivariateDual<Double> high) {
     726      this.low = low.Clone();
     727      this.high = high.Clone();
     728    }
     729
    520730    public AlgebraicInterval(double low, double high) {
    521       this.low = low;
    522       this.high = high;
     731      this.low = new MultivariateDual<Double>(new Double(low));
     732      this.high = new MultivariateDual<Double>(new Double(high));
    523733    }
    524734
    525735    public AlgebraicInterval Add(AlgebraicInterval a) {
    526       low += a.low;
    527       high += a.high;
     736      low.Add(a.low);
     737      high.Add(a.high);
    528738      return this;
    529739    }
    … …  
    540750
    541751    public AlgebraicInterval Div(AlgebraicInterval a) {
    542       if (Contains(0.0)) {
    543         if (a.low.IsAlmost(0.0))
    544           Mul(new AlgebraicInterval(1.0 / a.high, double.PositiveInfinity));
    545         else if (a.high.IsAlmost(0.0))
    546           Mul(new AlgebraicInterval(double.NegativeInfinity, 1.0 / a.low));
    547         else {
    548           low = double.NegativeInfinity;
    549           high = double.PositiveInfinity;
    550         }
     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()));
    551760      } else {
    552         Mul(new AlgebraicInterval(1.0 / a.high, 1.0 / a.low));
     761        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
    553762      }
    554763      return this;
    … …  
    556765
    557766    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
    558       low = Math.Exp(a.low);
    559       high = Math.Exp(a.high);
     767      low.AssignExp(a.low);
     768      high.AssignExp(a.high);
    560769      return this;
    561770    }
    … …  
    570779
    571780    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
    572       low = 1.0;
    573       high = 1.0;
     781      low = new MultivariateDual<Double>(1.0);
     782      high = new MultivariateDual<Double>(1.0);
    574783      return Div(a);
    575784    }
    576785
    577786    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
    578       low = Math.Log(a.low);
    579       high = Math.Log(a.high);
     787      low.AssignLog(a.low);
     788      high.AssignLog(a.high);
    580789      return this;
    581790    }
    582791
    583792    public AlgebraicInterval Mul(AlgebraicInterval a) {
    584       double v1 = low * a.low;
    585       double v2 = low * a.high;
    586       double v3 = high * a.low;
    587       double v4 = high * a.high;
    588 
    589       low = Math.Min(Math.Min(v1, v2), Math.Min(v3, v4));
    590       high = Math.Max(Math.Max(v1, v2), Math.Max(v3, v4));
     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));
    591800      return this;
    592801    }
    593802
    594803    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
    595       low = 0;
    596       high = 0;
    597       return Sub(a);
     804      throw new NotImplementedException();
    598805    }
    599806
    600807    public AlgebraicInterval Scale(double s) {
     808      low.Scale(s);
     809      high.Scale(s);
    601810      if (s < 0) {
    602         low = s * high;
    603         high = s * low;
    604       } else {
    605         low = s * low;
    606         high = s * high;
     811        var t = low;
     812        low = high;
     813        high = t;
    607814      }
    608815      return this;
    … …  
    614821
    615822    public AlgebraicInterval Sub(AlgebraicInterval a) {
    616       low -= a.high;
    617       high -= a.low;
     823      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
     824      low.Sub(a.high);
     825      high.Sub(a.low);
    618826      return this;
    619827    }
    … …  
    624832
    625833    public bool Contains(double val) {
    626       return low <= val && val <= high;
    627     }
    628   }
     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
    629857
    630858  public class Dual<V> : IAlgebraicType<Dual<V>>
    … …  
    633861    private V dv;
    634862
     863    public V Value => v;
     864    public V Derivative => dv;
     865
    635866    public Dual(V v, V dv) {
    636867      this.v = v;
    … …  
    720951      throw new NotImplementedException();
    721952    }
    722   }
    723 
    724   public class VectorDual<V> : IAlgebraicType<VectorDual<V>> where V : IAlgebraicType<V> {
     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() {
    725971    private V v;
    726972    public V Value => v;
    727973
    728     private Vector<V> dv;
    729     public Vector<V> Gradient => dv; // column-oriented (dv[0] is the vector for the first parameter)
    730 
    731 
    732     public VectorDual(V v, Vector<V> dv) {
    733       this.v = v;
    734       this.dv = dv;
    735     }
    736 
    737     public VectorDual<V> Clone() {
    738       return new VectorDual<V>(v.Clone(), dv.Clone());
    739     }
    740 
    741     public VectorDual<V> Add(VectorDual<V> a) {
     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) {
    7421018      v.Add(a.v);
    743       dv.Concat(a.dv);
    744       return this;
    745     }
    746 
    747     public VectorDual<V> Assign(VectorDual<V> a) {
     1019      dv.Add(a.dv);
     1020      return this;
     1021    }
     1022
     1023    public MultivariateDual<V> Assign(MultivariateDual<V> a) {
    7481024      v.Assign(a.v);
    749       dv = a.dv.Clone();
    750       return this;
    751     }
    752 
    753     public VectorDual<V> AssignCos(VectorDual<V> a) {
    754       v.AssignCos(a.v);
    755       V t = default(V);
    756       t.AssignNeg(t.AssignSin(a.v));
    757       dv = a.dv.Clone().Scale(t);
    758       return this;
    759     }
    760 
    761     public VectorDual<V> AssignExp(VectorDual<V> a) {
    762       v.AssignExp(a.v);
    763       dv = a.dv.Clone().Scale(v);
    764       return this;
    765     }
    766 
    767     public VectorDual<V> AssignIntPower(VectorDual<V> a, int p) {
    768       throw new NotImplementedException();
    769     }
    770 
    771     public VectorDual<V> AssignIntRoot(VectorDual<V> a, int r) {
    772       throw new NotImplementedException();
    773     }
    774 
    775     public VectorDual<V> AssignInv(VectorDual<V> a) {
    776       throw new NotImplementedException();
    777     }
    778 
    779     public VectorDual<V> AssignLog(VectorDual<V> a) {
    780       throw new NotImplementedException();
    781     }
    782 
    783     public VectorDual<V> AssignNeg(VectorDual<V> a) {
    784       throw new NotImplementedException();
    785     }
    786 
    787     public VectorDual<V> AssignSin(VectorDual<V> a) {
    788       throw new NotImplementedException();
    789     }
    790 
    791     public VectorDual<V> Div(VectorDual<V> a) {
    792       throw new NotImplementedException();
    793     }
    794 
    795     public VectorDual<V> Mul(VectorDual<V> a) {
     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) {
    7961066      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
    7971067
    7981068      var t1 = a.dv.Clone().Scale(v);
    7991069      var t2 = dv.Clone().Scale(a.v);
    800       dv.Assign(t1).Concat(t2); // TODO: CHECK ORDER
     1070      dv.Assign(t1).Add(t2);
    8011071
    8021072      v.Mul(a.v);
    … …  
    8041074    }
    8051075
    806     public VectorDual<V> Scale(double s) {
     1076    public MultivariateDual<V> Scale(double s) {
    8071077      v.Scale(s);
    8081078      dv.Scale(s);
    … …  
    8101080    }
    8111081
    812     public VectorDual<V> Sub(VectorDual<V> a) {
     1082    public MultivariateDual<V> Sub(MultivariateDual<V> a) {
    8131083      v.Sub(a.v);
    814       dv.Concat(a.dv.AssignNeg(a.dv));
     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>();
    8151100      return this;
    8161101    }
Note: See TracChangeset for help on using the changeset viewer.