source: branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicInterval.cs @ 17296

using System;
using System.Diagnostics;
using System.Linq;

namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
  // this is our own implementation of interval arithmetic
  // for a well worked out definition of interval operations for IEEE reals see:
  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    private MultivariateDual<AlgebraicDouble> low;
    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();

    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    private MultivariateDual<AlgebraicDouble> high;
    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();


    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }

    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
      this.low = low.Clone();
      this.high = high.Clone();
    }

    public AlgebraicInterval(double low, double high) {
      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
    }

    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);

    public AlgebraicInterval Add(AlgebraicInterval a) {
      low.Add(a.low);
      high.Add(a.high);
      return this;
    }

    public AlgebraicInterval Mul(AlgebraicInterval a) {
      var v1 = low.Clone().Mul(a.low);
      var v2 = low.Clone().Mul(a.high);
      var v3 = high.Clone().Mul(a.low);
      var v4 = high.Clone().Mul(a.high);

      low = Min(Min(v1, v2), Min(v3, v4));
      high = Max(Max(v1, v2), Max(v3, v4));

      return this;
    }


    private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
      return a.Value < b.Value ? a : b;
    }
    private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
      return a.Value > b.Value ? a : b;
    }

    public AlgebraicInterval Assign(AlgebraicInterval a) {
      low = a.low;
      high = a.high;
      return this;
    }

    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
    }

    public AlgebraicInterval Div(AlgebraicInterval a) {
      if (a.Contains(0.0)) {
        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
          if (this.low.Value >= 0) {
            // pos / 0
          } else if (this.high.Value <= 0) {
            // neg / 0
          } else {
            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
          }
        } else if (a.low.Value.Value >= 0) {
          // a is positive 
          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
        } else if (a.high.Value <= 0) {
          // a is negative
          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
        } else {
          // a is interval over zero
          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
        }
      } else {
        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
      }
      return this;
    }

    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
      low.AssignExp(a.low);
      high.AssignExp(a.high);
      return this;
    }

    // tanh is a bijective function
    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
      low.AssignTanh(a.low);
      high.AssignTanh(a.high);
      return this;
    }

    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
      if (p < 0) {  // x^-3 == 1/(x^3)
        AssignIntPower(a, -p);
        return AssignInv(this);
      } else if (p == 0) {
        if (a.Contains(0.0)) {
          // => 0^0 = 0 ; might not be relevant
          low = new MultivariateDual<AlgebraicDouble>(0.0);
          high = new MultivariateDual<AlgebraicDouble>(1.0);
          return this;
        } else {
          // => 1
          low = new MultivariateDual<AlgebraicDouble>(1.0);
          high = new MultivariateDual<AlgebraicDouble>(1.0);
          return this;
        }
      } else if (p == 1) return this;
      else if (p % 2 == 0) {
        // p is even => interval must be positive
        if (a.Contains(0.0)) {
          low = new MultivariateDual<AlgebraicDouble>(0.0);
          high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));
        } else {
          var lowPower = a.low.IntPower(p);
          var highPower = a.high.IntPower(p);
          low = lowPower.AssignMin(highPower);
          high = lowPower.AssignMax(highPower);
        }
      } else {
        // p is uneven
        if (a.Contains(0.0)) {
          low.AssignIntPower(a.low, p);
          high.AssignIntPower(a.high, p);
        } else {
          var lowPower = a.low.IntPower(p);
          var highPower = a.high.IntPower(p);
          low = lowPower.AssignMin(highPower);
          high = lowPower.AssignMax(highPower);
        }
      }
      return this;
    }

    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
      if (r == 1) return this;
      if (r < 0) {
        // x^ (-1/2) = 1 / (x^(1/2))
        AssignIntRoot(a, -r);
        return AssignInv(this);
      } else {
        // root only defined for positive arguments for even roots
        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
          return this;
        } else {
          low.AssignIntRoot(a.low, r);
          high.AssignIntRoot(a.high, r);
          return this;
        }
      }
    }

    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
      low = new MultivariateDual<AlgebraicDouble>(1.0);
      high = new MultivariateDual<AlgebraicDouble>(1.0);
      return Div(a);
    }

    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
      low.AssignLog(a.low);
      high.AssignLog(a.high);
      return this;
    }

    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
      low.AssignNeg(a.high);
      high.AssignNeg(a.low);
      return this;
    }

    public AlgebraicInterval Scale(double s) {
      low.Scale(s);
      high.Scale(s);
      if (s < 0) {
        var t = low;
        low = high;
        high = t;
      }
      return this;
    }

    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
      var lower = a.LowerBound.Value.Value;
      var size = a.UpperBound.Value.Value - lower;
      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;

      if (size >= Math.PI * 2) {
        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
        high = new MultivariateDual<AlgebraicDouble>(1.0);
        return this;
      }

      // assume low and high are in the same quadrant
      low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
      high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());

      // override min and max if necessary

      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]      

      // handle min = -1 and max = 1 cases explicitly
      var pi_2 = Math.PI / 2.0;
      var maxima = new double[] { -3 * pi_2, pi_2 };
      var minima = new double[] { -pi_2, 3 * pi_2 };

      // override min and max if necessary
      if (maxima.Any(m => lower < m && lower + size > m)) {
        // max = 1
        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
      }

      if (minima.Any(m => lower < m && lower + size > m)) {
        // min = -1;
        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
      }
      return this;
    }

    public AlgebraicInterval Sub(AlgebraicInterval a) {
      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
      low.Sub(a.high);
      high.Sub(a.low);
      return this;
    }

    public AlgebraicInterval Clone() {
      return new AlgebraicInterval(low, high);
    }

    public bool Contains(double val) {
      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
    }

    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
      if (a.Contains(0.0)) {
        var abslow = a.low.Clone().Abs();
        var abshigh = a.high.Clone().Abs();
        a.high.Assign(Algebraic.Max(abslow, abshigh));
        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
      } else {
        var abslow = a.low.Clone().Abs();
        var abshigh = a.high.Clone().Abs();
        a.low.Assign(Algebraic.Min(abslow, abshigh));
        a.high.Assign(Algebraic.Max(abslow, abshigh));
      }
      return this;
    }

    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
      low.AssignSgn(a.low);
      high.AssignSgn(a.high);
      return this;
    }

    public AlgebraicInterval AssignMin(AlgebraicInterval other) {
      low.AssignMin(other.low);
      high.AssignMin(other.high);
      return this;
    }

    public AlgebraicInterval AssignMax(AlgebraicInterval other) {
      low.AssignMax(other.low);
      high.AssignMax(other.high);
      return this;
    }
  }
}