Free cookie consent management tool by TermsFeed Policy Generator

Changeset 16693


Ignore:
Timestamp:
03/19/19 10:19:40 (5 years ago)
Author:
gkronber
Message:

#2994: worked on implementations of remaining methods for IAlgebraicType

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r16682 r16693  
    44using HeuristicLab.Common;
    55using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
    6 using HEAL.Attic;
    76
    87namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
    … …  
    117116
    118117
    119   public class VectorEvaluator : Interpreter<DoubleVector> {
     118  public sealed class VectorEvaluator : Interpreter<DoubleVector> {
    120119    private const int BATCHSIZE = 128;
    121120    [ThreadStatic]
    … …  
    192191  }
    193192
    194   public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
     193  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
    195194    private const int BATCHSIZE = 128;
    196195    [ThreadStatic]
    … …  
    318317
    319318
    320   public class IntervalEvaluator : Interpreter<AlgebraicInterval> {
     319  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
    321320    [ThreadStatic]
    322321    private Dictionary<string, Interval> intervals;
    … …  
    352351      instruction.dblVal = constant.Value;
    353352      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
     353        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0),
     354        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key
    356355        );
    357356    }
    … …  
    360359      instruction.dblVal = variable.Weight;
    361360      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)
     361        low: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound),  // bounds change by variable value d/dc (c I(var)) = I(var)
     362        high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
    364363        );
    365364    }
    … …  
    371370
    372371  public interface IAlgebraicType<T> {
     372    T Zero { get; }
     373
    373374    T AssignAbs(T a); // set this to assign abs(a)
    374375    T Assign(T a); // assign this to same value as a (copy!)
    … …  
    414415
    415416  // 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; }
     417  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
     418    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
     419    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
    419420    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 
     421
     422    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
     423    public AlgebraicDouble() { }
     424    public AlgebraicDouble(double value) { this.Value = value; }
     425    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
     426    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
     427    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
     428    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
     429    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
     430    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
     431    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
     432    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
     433    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
     434    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
     435    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
     436    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
     437    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
     438    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
     439    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
     440    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; }
     441    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
    441442  }
    442443
    … …  
    445446    private double[] arr;
    446447
     448
    447449    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
     450
     451    public int Length => arr.Length;
    448452
    449453    public DoubleVector(int length) {
    … …  
    461465    }
    462466
    463     public DoubleVector(int length, double constantValue) : this(length) {
    464       for (int i = 0; i < length; ++i) arr[i] = constantValue;
    465     }
    466 
    467     public void CopyTo(double[] dest, int idx, int length) {
    468       Array.Copy(arr, 0, dest, idx, length);
    469     }
    470 
    471     public void CopyRowTo(double[,] dest, int row) {
    472       for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
    473     }
    474     internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
    475       for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
    476     }
    477 
    478     public void AssignConstant(double constantValue) {
    479       for (int i = 0; i < arr.Length; ++i) {
    480         arr[i] = constantValue;
    481       }
    482     }
    483 
    484     public DoubleVector Assign(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
    485     public DoubleVector AssignCos(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
    486     public DoubleVector Div(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
    487     public DoubleVector AssignExp(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
    488     public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }
    489     public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }
    490     public DoubleVector AssignInv(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
    491     public DoubleVector AssignLog(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
     467    public DoubleVector Zero => new DoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
     468    public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
    492469    public DoubleVector Add(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
     470    public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
    493471    public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
    494     public DoubleVector AssignNeg(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
     472    public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
     473    public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
     474    public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
    495475    public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
    496     public DoubleVector AssignSin(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
    497     public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
    498     public DoubleVector AssignAbs(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
    499     public DoubleVector AssignSgn(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
     476    public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
     477    public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
     478    public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
     479    public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
     480    public DoubleVector AssignIntPower(DoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }
     481    public DoubleVector AssignIntRoot(DoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], 1.0 / r); } return this; }
     482    public DoubleVector AssignAbs(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
     483    public DoubleVector AssignSgn(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
     484
    500485    public DoubleVector Clone() {
    501486      var v = new DoubleVector(this.arr.Length);
    … …  
    504489    }
    505490
     491    public void AssignConstant(double constantValue) {
     492      for (int i = 0; i < arr.Length; ++i) {
     493        arr[i] = constantValue;
     494      }
     495    }
     496
     497    public void CopyTo(double[] dest, int idx, int length) {
     498      Array.Copy(arr, 0, dest, idx, length);
     499    }
     500
    506501    public void CopyFrom(double[] data, int rowIndex) {
    507502      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
    508503    }
    509504
     505    public void CopyRowTo(double[,] dest, int row) {
     506      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
     507    }
     508
     509    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
     510      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
     511    }
    510512  }
    511513
    512514  // vectors of algebraic types
    513   public class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {
     515  public sealed class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {
    514516    private T[] elems;
    515517
    … …  
    517519
    518520    public int Length => elems.Length;
     521
    519522
    520523    private Vector() { }
    … …  
    527530    ///
    528531    /// </summary>
    529     /// <param name="elems">The array is copied</param>
     532    /// <param name="elems">The array is copied (element-wise clone)</param>
    530533    public Vector(T[] elems) {
    531534      this.elems = new T[elems.Length];
    … …  
    562565    }
    563566
     567    public Vector<T> Zero => new Vector<T>(Length);
     568    public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
    564569    public Vector<T> Add(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
    565     public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
     570    public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
     571    public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
     572    public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
     573    public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
     574    public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
     575    public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
     576    public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
     577    public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
     578    public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
     579    public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
    566580    public Vector<T> AssignCos(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
    567     public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
    568581    public Vector<T> AssignIntPower(Vector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }
    569582    public Vector<T> AssignIntRoot(Vector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }
    570     public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
    571     public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
    572     public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
    573     public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
    574     public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
    575     public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
    576     public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
    577     public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
    578     public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
    579583    public Vector<T> AssignAbs(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
    580584    public Vector<T> AssignSgn(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
    … …  
    587591  /// <typeparam name="K">Key type</typeparam>
    588592  /// <typeparam name="T">Element type</typeparam>
    589   public class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {
     593  public sealed class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {
    590594
    591595    private Dictionary<K, T> elems;
    592 
    593596    public IReadOnlyDictionary<K, T> Elements => elems;
     597
    594598
    595599    public SparseVector(SparseVector<K, T> original) {
    … …  
    597601    }
    598602
     603    /// <summary>
     604    ///
     605    /// </summary>
     606    /// <param name="keys"></param>
     607    /// <param name="values">values are cloned</param>
    599608    public SparseVector(K[] keys, T[] values) {
    600609      if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
    601610      elems = new Dictionary<K, T>(keys.Length);
    602611      for (int i = 0; i < keys.Length; ++i) {
    603         elems.Add(keys[i], values[i]);
     612        elems.Add(keys[i], values[i].Clone());
    604613      }
    605614    }
    … …  
    609618    }
    610619
    611     public SparseVector<K, T> Add(SparseVector<K, T> a) {
     620
     621
     622    private void AssignTransformed(SparseVector<K, T> a, Func<T, T, T> mapAssign) {
    612623      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());
     624        if (elems.TryGetValue(kvp.Key, out T value))
     625          mapAssign(kvp.Value, value);
     626        else {
     627          var newValue = kvp.Value.Zero;
     628          elems.Add(kvp.Key, newValue);
     629          mapAssign(kvp.Value, newValue);
    617630        }
    618631      }
    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     }
     632    }
     633
     634    public SparseVector<K, T> Zero => new SparseVector<K, T>();
     635
     636    public SparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
     637    public SparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
     638
     639    public SparseVector<K, T> Assign(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Assign(src)); return this; }
     640    public SparseVector<K, T> Add(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Add(src)); return this; }
     641    public SparseVector<K, T> Mul(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Mul(src)); return this; }
     642    public SparseVector<K, T> Sub(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Sub(src)); return this; }
     643    public SparseVector<K, T> Div(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Div(src)); return this; }
     644    public SparseVector<K, T> AssignInv(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignInv(src)); return this; }
     645    public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignNeg(src)); return this; }
     646    public SparseVector<K, T> AssignLog(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignLog(src)); return this; }
     647    public SparseVector<K, T> AssignExp(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignExp(src)); return this; }
     648    public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) { AssignTransformed(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
     649    public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) { AssignTransformed(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
     650    public SparseVector<K, T> AssignSin(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSin(src)); return this; }
     651    public SparseVector<K, T> AssignCos(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignCos(src)); return this; }
     652    public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignAbs(src)); return this; }
     653    public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSgn(src)); return this; }
    674654
    675655    public SparseVector<K, T> Clone() {
    676656      return new SparseVector<K, T>(this);
    677657    }
    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     }
    713658  }
    714659
    715660
    716661  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
    717     private MultivariateDual<Double> low;
    718     private MultivariateDual<Double> high;
    719 
    720     public MultivariateDual<Double> LowerBound => low.Clone();
    721     public MultivariateDual<Double> UpperBound => high.Clone();
     662    private MultivariateDual<AlgebraicDouble> low;
     663    private MultivariateDual<AlgebraicDouble> high;
     664
     665    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
     666    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
     667
    722668
    723669    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
    724670
    725     public AlgebraicInterval(MultivariateDual<Double> low, MultivariateDual<Double> high) {
     671    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
    726672      this.low = low.Clone();
    727673      this.high = high.Clone();
    … …  
    729675
    730676    public AlgebraicInterval(double low, double high) {
    731       this.low = new MultivariateDual<Double>(new Double(low));
    732       this.high = new MultivariateDual<Double>(new Double(high));
    733     }
    734 
     677      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
     678      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
     679    }
     680
     681    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
    735682    public AlgebraicInterval Add(AlgebraicInterval a) {
    736683      low.Add(a.low);
    737684      high.Add(a.high);
    738       return this;
    739     }
    740 
    741     public AlgebraicInterval Assign(AlgebraicInterval a) {
    742       low = a.low;
    743       high = a.high;
    744       return this;
    745     }
    746 
    747     public AlgebraicInterval AssignCos(AlgebraicInterval a) {
    748       throw new NotImplementedException();
    749     }
    750 
    751     public AlgebraicInterval Div(AlgebraicInterval a) {
    752       if (a.Contains(0.0)) {
    753         if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
    754           low = new MultivariateDual<Double>(double.NegativeInfinity);
    755           high = new MultivariateDual<Double>(double.PositiveInfinity);
    756         } else if (a.low.Value.Value.IsAlmost(0.0))
    757           Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<Double>(double.PositiveInfinity)));
    758         else
    759           Mul(new AlgebraicInterval(new MultivariateDual<Double>(double.NegativeInfinity), a.low.Clone().Inv()));
    760       } else {
    761         Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
    762       }
    763       return this;
    764     }
    765 
    766     public AlgebraicInterval AssignExp(AlgebraicInterval a) {
    767       low.AssignExp(a.low);
    768       high.AssignExp(a.high);
    769       return this;
    770     }
    771 
    772     public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
    773       throw new NotImplementedException();
    774     }
    775 
    776     public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
    777       throw new NotImplementedException();
    778     }
    779 
    780     public AlgebraicInterval AssignInv(AlgebraicInterval a) {
    781       low = new MultivariateDual<Double>(1.0);
    782       high = new MultivariateDual<Double>(1.0);
    783       return Div(a);
    784     }
    785 
    786     public AlgebraicInterval AssignLog(AlgebraicInterval a) {
    787       low.AssignLog(a.low);
    788       high.AssignLog(a.high);
    789685      return this;
    790686    }
    … …  
    801697    }
    802698
     699    public AlgebraicInterval Assign(AlgebraicInterval a) {
     700      low = a.low;
     701      high = a.high;
     702      return this;
     703    }
     704
     705    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
     706      return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
     707    }
     708
     709    public AlgebraicInterval Div(AlgebraicInterval a) {
     710      if (a.Contains(0.0)) {
     711        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
     712          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
     713          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
     714        } else if (a.low.Value.Value.IsAlmost(0.0))
     715          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
     716        else
     717          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
     718      } else {
     719        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
     720      }
     721      return this;
     722    }
     723
     724    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
     725      low.AssignExp(a.low);
     726      high.AssignExp(a.high);
     727      return this;
     728    }
     729
     730    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
     731      if (p == 0) {
     732        // => 1
     733        low = new MultivariateDual<AlgebraicDouble>(1.0);
     734        high = new MultivariateDual<AlgebraicDouble>(1.0);
     735        return this;
     736      }
     737      if (p == 1) return this;
     738
     739      if (p < 0) {  // x^-3 == 1/(x^3)
     740        AssignIntPower(a, -p);
     741        return AssignInv(this);
     742      } else {
     743        // p is even => interval must be positive
     744        if (p % 2 == 0) {
     745          if (a.Contains(0.0)) {
     746            low = new MultivariateDual<AlgebraicDouble>(0.0);
     747            high = Algebraic.Max(low.Clone().IntPower(p), high.Clone().IntPower(p));
     748          } else {
     749            var lowPower = low.Clone().IntPower(p);
     750            var highPower = high.Clone().IntPower(p);
     751            low = Algebraic.Min(lowPower, highPower);
     752            high = Algebraic.Max(lowPower, highPower);
     753          }
     754        } else {
     755          // p is uneven
     756          var lowPower = low.Clone().IntPower(p);
     757          var highPower = high.Clone().IntPower(p);
     758          low = Algebraic.Min(lowPower, highPower);
     759          high = Algebraic.Max(lowPower, highPower);
     760        }
     761        return this;
     762      }
     763    }
     764
     765    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
     766      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
     767      if (r == 1) return this;
     768      if (r < 0) {
     769        // x^ (-1/2) = 1 / (x^(1/2))
     770        AssignIntRoot(a, -r);
     771        return AssignInv(this);
     772      } else {
     773        // root only defined for positive arguments
     774        if (a.LowerBound.Value.Value < 0) {
     775          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
     776          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
     777          return this;
     778        } else {
     779          low.AssignIntRoot(a.low, r);
     780          high.AssignIntRoot(a.high, r);
     781          return this;
     782        }
     783      }
     784    }
     785
     786    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
     787      low = new MultivariateDual<AlgebraicDouble>(1.0);
     788      high = new MultivariateDual<AlgebraicDouble>(1.0);
     789      return Div(a);
     790    }
     791
     792    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
     793      low.AssignLog(a.low);
     794      high.AssignLog(a.high);
     795      return this;
     796    }
     797
    803798    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
    804       throw new NotImplementedException();
     799      low.AssignNeg(a.high);
     800      high.AssignNeg(a.low);
     801      return this;
    805802    }
    806803
    … …  
    817814
    818815    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
    819       throw new NotImplementedException();
     816      if (Math.Abs(a.UpperBound.Value.Value - a.LowerBound.Value.Value) >= Math.PI * 2) {
     817        low = new MultivariateDual<AlgebraicDouble>(-1.0);
     818        high = new MultivariateDual<AlgebraicDouble>(1.0);
     819      }
     820
     821      //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)
     822      double Pihalf = Math.PI / 2;
     823      var scaled = this.Clone().Scale(1.0 / Pihalf);
     824      //move to positive scale
     825      if (scaled.LowerBound.Value.Value < 0) {
     826        int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1;
     827        scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4));
     828      }
     829
     830      double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0;
     831      double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0;
     832      if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0;
     833      List<double> sinValues = new List<double>();
     834      sinValues.Add(Math.Sin(scaledLowerBound * Pihalf));
     835      sinValues.Add(Math.Sin(scaledUpperBound * Pihalf));
     836
     837      int startValue = (int)Math.Ceiling(scaledLowerBound);
     838      while (startValue < scaledUpperBound) {
     839        sinValues.Add(Math.Sin(startValue * Pihalf));
     840        startValue += 1;
     841      }
     842
     843      low = new MultivariateDual<AlgebraicDouble>(sinValues.Min());
     844      high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());
     845      return this;
    820846    }
    821847
    … …  
    840866        var abshigh = a.high.Clone().Abs();
    841867        a.high.Assign(Algebraic.Max(abslow, abshigh));
    842         a.low.Assign(new MultivariateDual<Double>(0.0)); // lost gradient for lower bound
     868        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
    843869      } else {
    844870        var abslow = a.low.Clone().Abs();
    … …  
    851877
    852878    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
    853       throw new NotImplementedException();
     879      low.AssignSgn(a.low);
     880      high.AssignSgn(a.high);
     881      return this;
    854882    }
    855883  }
    … …  
    864892    public V Derivative => dv;
    865893
     894
    866895    public Dual(V v, V dv) {
    867896      this.v = v;
    … …  
    873902    }
    874903
     904    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
     905
    875906    public Dual<V> Add(Dual<V> a) {
    876907      v.Add(a.v);
    … …  
    892923
    893924    public Dual<V> Div(Dual<V> a) {
    894       throw new NotImplementedException();
     925      Mul(a.Inv());
     926      return this;
    895927    }
    896928
    … …  
    902934
    903935    public Dual<V> AssignIntPower(Dual<V> a, int p) {
    904       throw new NotImplementedException();
     936      v.AssignIntPower(a.v, p);
     937      dv.Assign(a.dv).Scale(p).Mul(a.v.IntPower(p - 1));
     938      return this;
    905939    }
    906940
    907941    public Dual<V> AssignIntRoot(Dual<V> a, int r) {
    908       throw new NotImplementedException();
     942      v.AssignIntRoot(a.v, r);
     943      dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r - 1));
     944      return this;
    909945    }
    910946
    911947    public Dual<V> AssignInv(Dual<V> a) {
    912       throw new NotImplementedException();
     948      v.AssignInv(a.v);
     949      dv.AssignNeg(a.dv).Mul(v).Mul(v); // (1/f(x))' = - f(x)' / f(x)^2
     950      return this;
    913951    }
    914952
    … …  
    935973
    936974    public Dual<V> AssignNeg(Dual<V> a) {
    937       throw new NotImplementedException();
     975      v.AssignNeg(a.v);
     976      dv.AssignNeg(a.dv);
     977      return this;
    938978    }
    939979
    … …  
    945985
    946986    public Dual<V> AssignSin(Dual<V> a) {
    947       throw new NotImplementedException();
     987      v.AssignSin(a.v);
     988      dv.Assign(a.dv).Mul(a.v.Clone().Cos());
     989      return this;
    948990    }
    949991
    950992    public Dual<V> Sub(Dual<V> a) {
    951       throw new NotImplementedException();
     993      v.Sub(a.v);
     994      dv.Sub(a.dv);
     995      return this;
    952996    }
    953997
    … …  
    10111055    }
    10121056
     1057    /// <summary>
     1058    ///
     1059    /// </summary>
     1060    /// <param name="v">not cloned</param>
     1061    /// <param name="gradient">not cloned</param>
     1062    internal MultivariateDual(V v, SparseVector<object, V> gradient) {
     1063      this.v = v;
     1064      this.dv = gradient;
     1065    }
     1066
    10131067    public MultivariateDual<V> Clone() {
    10141068      return new MultivariateDual<V>(this);
    10151069    }
     1070
     1071
     1072    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
     1073
    10161074
    10171075    public MultivariateDual<V> Add(MultivariateDual<V> a) {
    … …  
    10971155      v.AssignSgn(a.v);
    10981156      // sign(f(x)) = 0;
    1099       dv = new SparseVector<object, V>();
     1157      dv = a.dv.Zero;
    11001158      return this;
    11011159    }
Note: See TracChangeset for help on using the changeset viewer.