1  using System;


2  using System.Diagnostics;


3  using System.Linq;


4 


5  namespace HeuristicLab.Problems.DataAnalysis.Symbolic {


6  // this is our own implementation of interval arithmetic


7  // for a well worked out definition of interval operations for IEEE reals see:


8  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995


9  [DebuggerDisplay("[{low.Value}..{high.Value}]")]


10  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {


11  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


12  private MultivariateDual<AlgebraicDouble> low;


13  public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();


14 


15  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


16  private MultivariateDual<AlgebraicDouble> high;


17  public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();


18 


19 


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


21 


22  public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {


23  this.low = low.Clone();


24  this.high = high.Clone();


25  }


26 


27  public AlgebraicInterval(double low, double high) {


28  this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));


29  this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));


30  }


31 


32  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


33  public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);


34  [DebuggerBrowsable(DebuggerBrowsableState.Never)]


35  public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);


36 


37  public AlgebraicInterval Add(AlgebraicInterval a) {


38  low.Add(a.low);


39  high.Add(a.high);


40  return this;


41  }


42 


43  public AlgebraicInterval Mul(AlgebraicInterval a) {


44  var v1 = low.Clone().Mul(a.low);


45  var v2 = low.Clone().Mul(a.high);


46  var v3 = high.Clone().Mul(a.low);


47  var v4 = high.Clone().Mul(a.high);


48 


49  low = Min(Min(v1, v2), Min(v3, v4));


50  high = Max(Max(v1, v2), Max(v3, v4));


51 


52  return this;


53  }


54 


55 


56  private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {


57  return a.Value < b.Value ? a : b;


58  }


59  private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {


60  return a.Value > b.Value ? a : b;


61  }


62 


63  public AlgebraicInterval Assign(AlgebraicInterval a) {


64  low = a.low;


65  high = a.high;


66  return this;


67  }


68 


69  public AlgebraicInterval AssignCos(AlgebraicInterval a) {


70  return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));


71  }


72 


73  public AlgebraicInterval Div(AlgebraicInterval a) {


74  if (a.Contains(0.0)) {


75  if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {


76  if (this.low.Value >= 0) {


77  // pos / 0


78  } else if (this.high.Value <= 0) {


79  // neg / 0


80  } else {


81  low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);


82  high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);


83  }


84  } else if (a.low.Value.Value >= 0) {


85  // a is positive


86  Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));


87  } else if (a.high.Value <= 0) {


88  // a is negative


89  Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));


90  } else {


91  // a is interval over zero


92  low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);


93  high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);


94  }


95  } else {


96  Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low


97  }


98  return this;


99  }


100 


101  public AlgebraicInterval AssignExp(AlgebraicInterval a) {


102  low.AssignExp(a.low);


103  high.AssignExp(a.high);


104  return this;


105  }


106 


107  // tanh is a bijective function


108  public AlgebraicInterval AssignTanh(AlgebraicInterval a) {


109  low.AssignTanh(a.low);


110  high.AssignTanh(a.high);


111  return this;


112  }


113 


114  public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {


115  if (p < 0) { // x^3 == 1/(x^3)


116  AssignIntPower(a, p);


117  return AssignInv(this);


118  } else if (p == 0) {


119  if (a.Contains(0.0)) {


120  // => 0^0 = 0 ; might not be relevant


121  low = new MultivariateDual<AlgebraicDouble>(0.0);


122  high = new MultivariateDual<AlgebraicDouble>(1.0);


123  return this;


124  } else {


125  // => 1


126  low = new MultivariateDual<AlgebraicDouble>(1.0);


127  high = new MultivariateDual<AlgebraicDouble>(1.0);


128  return this;


129  }


130  } else if (p == 1) return this;


131  else if (p % 2 == 0) {


132  // p is even => interval must be positive


133  if (a.Contains(0.0)) {


134  low = new MultivariateDual<AlgebraicDouble>(0.0);


135  high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));


136  } else {


137  var lowPower = a.low.IntPower(p);


138  var highPower = a.high.IntPower(p);


139  low = lowPower.AssignMin(highPower);


140  high = lowPower.AssignMax(highPower);


141  }


142  } else {


143  // p is uneven


144  if (a.Contains(0.0)) {


145  low.AssignIntPower(a.low, p);


146  high.AssignIntPower(a.high, p);


147  } else {


148  var lowPower = a.low.IntPower(p);


149  var highPower = a.high.IntPower(p);


150  low = lowPower.AssignMin(highPower);


151  high = lowPower.AssignMax(highPower);


152  }


153  }


154  return this;


155  }


156 


157  public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {


158  if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }


159  if (r == 1) return this;


160  if (r < 0) {


161  // x^ (1/2) = 1 / (x^(1/2))


162  AssignIntRoot(a, r);


163  return AssignInv(this);


164  } else {


165  // root only defined for positive arguments for even roots


166  if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {


167  low = new MultivariateDual<AlgebraicDouble>(double.NaN);


168  high = new MultivariateDual<AlgebraicDouble>(double.NaN);


169  return this;


170  } else {


171  low.AssignIntRoot(a.low, r);


172  high.AssignIntRoot(a.high, r);


173  return this;


174  }


175  }


176  }


177 


178  public AlgebraicInterval AssignInv(AlgebraicInterval a) {


179  low = new MultivariateDual<AlgebraicDouble>(1.0);


180  high = new MultivariateDual<AlgebraicDouble>(1.0);


181  return Div(a);


182  }


183 


184  public AlgebraicInterval AssignLog(AlgebraicInterval a) {


185  low.AssignLog(a.low);


186  high.AssignLog(a.high);


187  return this;


188  }


189 


190  public AlgebraicInterval AssignNeg(AlgebraicInterval a) {


191  low.AssignNeg(a.high);


192  high.AssignNeg(a.low);


193  return this;


194  }


195 


196  public AlgebraicInterval Scale(double s) {


197  low.Scale(s);


198  high.Scale(s);


199  if (s < 0) {


200  var t = low;


201  low = high;


202  high = t;


203  }


204  return this;


205  }


206 


207  public AlgebraicInterval AssignSin(AlgebraicInterval a) {


208  var lower = a.LowerBound.Value.Value;


209  var size = a.UpperBound.Value.Value  lower;


210  if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;


211 


212  if (size >= Math.PI * 2) {


213  low = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient


214  high = new MultivariateDual<AlgebraicDouble>(1.0);


215  return this;


216  }


217 


218  // assume low and high are in the same quadrant


219  low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());


220  high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());


221 


222  // override min and max if necessary


223 


224  // shift interval 'a' into the range [2pi .. 2pi] without changing the size of the interval to simplify the checks


225  lower = lower % (2 * Math.PI); // lower in [2pi .. 2pi]


226 


227  // handle min = 1 and max = 1 cases explicitly


228  var pi_2 = Math.PI / 2.0;


229  var maxima = new double[] { 3 * pi_2, pi_2 };


230  var minima = new double[] { pi_2, 3 * pi_2 };


231 


232  // override min and max if necessary


233  if (maxima.Any(m => lower < m && lower + size > m)) {


234  // max = 1


235  high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient


236  }


237 


238  if (minima.Any(m => lower < m && lower + size > m)) {


239  // min = 1;


240  low = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient


241  }


242  return this;


243  }


244 


245  public AlgebraicInterval Sub(AlgebraicInterval a) {


246  // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]


247  low.Sub(a.high);


248  high.Sub(a.low);


249  return this;


250  }


251 


252  public AlgebraicInterval Clone() {


253  return new AlgebraicInterval(low, high);


254  }


255 


256  public bool Contains(double val) {


257  return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;


258  }


259 


260  public AlgebraicInterval AssignAbs(AlgebraicInterval a) {


261  if (a.Contains(0.0)) {


262  var abslow = a.low.Clone().Abs();


263  var abshigh = a.high.Clone().Abs();


264  a.high.Assign(Algebraic.Max(abslow, abshigh));


265  a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound


266  } else {


267  var abslow = a.low.Clone().Abs();


268  var abshigh = a.high.Clone().Abs();


269  a.low.Assign(Algebraic.Min(abslow, abshigh));


270  a.high.Assign(Algebraic.Max(abslow, abshigh));


271  }


272  return this;


273  }


274 


275  public AlgebraicInterval AssignSgn(AlgebraicInterval a) {


276  low.AssignSgn(a.low);


277  high.AssignSgn(a.high);


278  return this;


279  }


280 


281  public AlgebraicInterval AssignMin(AlgebraicInterval other) {


282  low.AssignMin(other.low);


283  high.AssignMin(other.high);


284  return this;


285  }


286 


287  public AlgebraicInterval AssignMax(AlgebraicInterval other) {


288  low.AssignMax(other.low);


289  high.AssignMax(other.high);


290  return this;


291  }


292  }


293  } 
