source: branches/3073_IA_constraint_splitting/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IABoundsEstimator.cs @ 18242

using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.Linq;
using HEAL.Attic;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Parameters;

namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
  [StorableType("C8539434-6FB0-47D0-9F5A-2CAE5D8B8B4F")]
  [Item("IA Bounds Estimator", "Interpreter for calculation of intervals of symbolic models.")]
  public sealed class IABoundsEstimator : ParameterizedNamedItem, IBoundsEstimator {
    #region Parameters

    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
    private const string UseIntervalSplittingParameterName = "Use Interval splitting";
    private const string SplittingIterationsParameterName = "Splitting Iterations";
    private const string SplittingWidthParameterName = "Splitting width";

    public IFixedValueParameter<IntValue> EvaluatedSolutionsParameter =>
      (IFixedValueParameter<IntValue>) Parameters[EvaluatedSolutionsParameterName];

    public IFixedValueParameter<BoolValue> UseIntervalSplittingParameter =>
      (IFixedValueParameter<BoolValue>) Parameters[UseIntervalSplittingParameterName];

    public IFixedValueParameter<IntValue> SplittingIterationsParameter =>
      (IFixedValueParameter<IntValue>) Parameters[SplittingIterationsParameterName];

    public IFixedValueParameter<DoubleValue> SplittingWidthParameter =>
      (IFixedValueParameter<DoubleValue>) Parameters[SplittingWidthParameterName];

    public int EvaluatedSolutions {
      get => EvaluatedSolutionsParameter.Value.Value;
      set => EvaluatedSolutionsParameter.Value.Value = value;
    }

    public bool UseIntervalSplitting {
      get => UseIntervalSplittingParameter.Value.Value;
      set => UseIntervalSplittingParameter.Value.Value = value;
    }

    public int SplittingIterations {
      get => SplittingIterationsParameter.Value.Value;
      set => SplittingIterationsParameter.Value.Value = value;
    }

    public double SplittingWidth {
      get => SplittingWidthParameter.Value.Value;
      set => SplittingWidthParameter.Value.Value = value;
    }

    #endregion

    #region Constructors

    [StorableConstructor]
    private IABoundsEstimator(StorableConstructorFlag _) : base(_) { }

    private IABoundsEstimator(IABoundsEstimator original, Cloner cloner) : base(original, cloner) { }

    public IABoundsEstimator() : base("IA Bounds Estimator",
      "Estimates the bounds of the model with interval arithmetic") {
      Parameters.Add(new FixedValueParameter<IntValue>(EvaluatedSolutionsParameterName,
        "A counter for the total number of solutions the estimator has evaluated.", new IntValue(0)));
      Parameters.Add(new FixedValueParameter<BoolValue>(UseIntervalSplittingParameterName,
        "Defines whether interval splitting is activated or not.", new BoolValue(false)));
      Parameters.Add(new FixedValueParameter<IntValue>(SplittingIterationsParameterName,
        "Defines the number of iterations of splitting.", new IntValue(200)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(SplittingWidthParameterName,
        "Width of interval, after the splitting should stop.", new DoubleValue(0.0)));
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new IABoundsEstimator(this, cloner);
    }

    #endregion

    #region IStatefulItem Members

    private readonly object syncRoot = new object();

    public void InitializeState() {
      EvaluatedSolutions = 0;
    }

    public void ClearState() { }

    #endregion

    #region Evaluation

    private static Instruction[] PrepareInterpreterState(
      ISymbolicExpressionTree tree,
      IDictionary<string, Interval> variableRanges) {
      if (variableRanges == null)
        throw new ArgumentNullException("No variablew ranges are present!", nameof(variableRanges));

      //Check if all variables used in the tree are present in the dataset
      foreach (var variable in tree.IterateNodesPrefix().OfType<VariableTreeNode>().Select(n => n.VariableName)
                                   .Distinct())
        if (!variableRanges.ContainsKey(variable))
          throw new InvalidOperationException($"No ranges for variable {variable} is present");

      var code = SymbolicExpressionTreeCompiler.Compile(tree, OpCodes.MapSymbolToOpCode);
      foreach (var instr in code.Where(i => i.opCode == OpCodes.Variable)) {
        var variableTreeNode = (VariableTreeNode) instr.dynamicNode;
        instr.data = variableRanges[variableTreeNode.VariableName];
      }

      return code;
    }

    public static Interval Evaluate(
      Instruction[] instructions, ref int instructionCounter,
      IDictionary<ISymbolicExpressionTreeNode, Interval> nodeIntervals = null,
      IDictionary<string, Interval> variableIntervals = null) {
      var currentInstr = instructions[instructionCounter];
      //Use ref parameter, because the tree will be iterated through recursively from the left-side branch to the right side
      //Update instructionCounter, whenever Evaluate is called
      instructionCounter++;
      Interval result = null;

      switch (currentInstr.opCode) {
        //Variables, Constants, ...
        case OpCodes.Variable: {
          var variableTreeNode = (VariableTreeNode) currentInstr.dynamicNode;
          var weightInterval = new Interval(variableTreeNode.Weight, variableTreeNode.Weight);

          Interval variableInterval;
          if (variableIntervals != null && variableIntervals.ContainsKey(variableTreeNode.VariableName))
            variableInterval = variableIntervals[variableTreeNode.VariableName];
          else
            variableInterval = (Interval) currentInstr.data;

          result = Interval.Multiply(variableInterval, weightInterval);
          break;
        }
        case OpCodes.Constant: {
          var constTreeNode = (ConstantTreeNode) currentInstr.dynamicNode;
          result = new Interval(constTreeNode.Value, constTreeNode.Value);
          break;
        }
        //Elementary arithmetic rules
        case OpCodes.Add: {
          result = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          for (var i = 1; i < currentInstr.nArguments; i++) {
            var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
            result = Interval.Add(result, argumentInterval);
          }

          break;
        }
        case OpCodes.Sub: {
          result = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          if (currentInstr.nArguments == 1)
            result = Interval.Multiply(new Interval(-1, -1), result);

          for (var i = 1; i < currentInstr.nArguments; i++) {
            var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
            result = Interval.Subtract(result, argumentInterval);
          }

          break;
        }
        case OpCodes.Mul: {
          result = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          for (var i = 1; i < currentInstr.nArguments; i++) {
            var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
            result = Interval.Multiply(result, argumentInterval);
          }

          break;
        }
        case OpCodes.Div: {
          result = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          if (currentInstr.nArguments == 1)
            result = Interval.Divide(new Interval(1, 1), result);

          for (var i = 1; i < currentInstr.nArguments; i++) {
            var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
            result = Interval.Divide(result, argumentInterval);
          }

          break;
        }
        //Trigonometric functions
        case OpCodes.Sin: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Sine(argumentInterval);
          break;
        }
        case OpCodes.Cos: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Cosine(argumentInterval);
          break;
        }
        case OpCodes.Tan: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Tangens(argumentInterval);
          break;
        }
        case OpCodes.Tanh: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.HyperbolicTangent(argumentInterval);
          break;
        }
        //Exponential functions
        case OpCodes.Log: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Logarithm(argumentInterval);
          break;
        }
        case OpCodes.Exp: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Exponential(argumentInterval);
          break;
        }
        case OpCodes.Square: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Square(argumentInterval);
          break;
        }
        case OpCodes.SquareRoot: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.SquareRoot(argumentInterval);
          break;
        }
        case OpCodes.Cube: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Cube(argumentInterval);
          break;
        }
        case OpCodes.CubeRoot: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.CubicRoot(argumentInterval);
          break;
        }
        case OpCodes.Absolute: {
          var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          result = Interval.Absolute(argumentInterval);
          break;
        }
        case OpCodes.AnalyticQuotient: {
          result = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
          for (var i = 1; i < currentInstr.nArguments; i++) {
            var argumentInterval = Evaluate(instructions, ref instructionCounter, nodeIntervals, variableIntervals);
            result = Interval.AnalyticalQuotient(result, argumentInterval);
          }

          break;
        }
        default:
          throw new NotSupportedException(
            $"The tree contains the unknown symbol {currentInstr.dynamicNode.Symbol}");
      }

      if (!(nodeIntervals == null || nodeIntervals.ContainsKey(currentInstr.dynamicNode)))
        nodeIntervals.Add(currentInstr.dynamicNode, result);

      return result;
    }

    #endregion

    #region Helpers

    private static IDictionary<string, Interval> GetOccurringVariableRanges(
      ISymbolicExpressionTree tree, IntervalCollection variableRanges) {
      var variables = tree.IterateNodesPrefix().OfType<VariableTreeNode>().Select(v => v.VariableName).Distinct()
                          .ToList();

      return variables.ToDictionary(x => x, x => variableRanges.GetReadonlyDictionary()[x]);
    }

    private static bool ContainsVariableMultipleTimes(ISymbolicExpressionTree tree, out List<String> variables) {
      variables = new List<string>();
      var varlist = tree.IterateNodesPrefix().OfType<VariableTreeNode>().GroupBy(x => x.VariableName);
      foreach (var group in varlist) {
        if (group.Count() > 1) {
          variables.Add(group.Key);
        }
      }

      return varlist.Any(group => group.Count() > 1);
    }

    // a multi-dimensional box with an associated bound
    // boxbounds are ordered first by bound (smaller first), then by size of box (larger first) then by distance of bottom left corner to origin
    private class BoxBound : IComparable<BoxBound> {
      public List<Interval> box;
      public double bound;

      public BoxBound(IEnumerable<Interval> box, double bound) {
        this.box = new List<Interval>(box);
        this.bound = bound;
      }

      public int CompareTo(BoxBound other) {
        if (bound != other.bound) return bound.CompareTo(other.bound);

        var thisSize = box.Aggregate(1.0, (current, dimExtent) => current * dimExtent.Width);
        var otherSize = other.box.Aggregate(1.0, (current, dimExtent) => current * dimExtent.Width);
        if (thisSize != otherSize) return -thisSize.CompareTo(otherSize);

        var thisDist = box.Sum(dimExtent => dimExtent.LowerBound * dimExtent.LowerBound);
        var otherDist = other.box.Sum(dimExtent => dimExtent.LowerBound * dimExtent.LowerBound);
        if (thisDist != otherDist) return thisDist.CompareTo(otherDist);

        // which is smaller first along the dimensions?
        for (int i = 0; i < box.Count; i++) {
          if (box[i].LowerBound != other.box[i].LowerBound) return box[i].LowerBound.CompareTo(other.box[i].LowerBound);
        }

        return 0;
      }
    }

    #endregion

    #region Splitting

    public static Interval EvaluateWithSplitting(Instruction[] instructions,
                                                 IDictionary<string, Interval> variableIntervals,
                                                 List<string> multipleOccurenceVariables, int splittingIterations,
                                                 double splittingWidth,
                                                 IDictionary<ISymbolicExpressionTreeNode, Interval> nodeIntervals =
                                                   null) {
      var min = FindBound(instructions, variableIntervals.ToDictionary(entry => entry.Key, entry => entry.Value),
        multipleOccurenceVariables, splittingIterations, splittingWidth, nodeIntervals,
        minimization: true);
      var max = FindBound(instructions, variableIntervals.ToDictionary(entry => entry.Key, entry => entry.Value),
        multipleOccurenceVariables, splittingIterations, splittingWidth, nodeIntervals,
        minimization: false);

      return new Interval(min, max);
    }

    private static double FindBound(Instruction[] instructions,
                                    IDictionary<string, Interval> variableIntervals,
                                    List<string> multipleOccurenceVariables, int splittingIterations,
                                    double splittingWidth,
                                    IDictionary<ISymbolicExpressionTreeNode, Interval> nodeIntervals = null,
                                    bool minimization = true, bool stopAtLimit = false, double limit = 0) {
      SortedSet<BoxBound> prioQ = new SortedSet<BoxBound>();
      var ic = 0;
      var stop = false;
      //Calculate full box
      var interval = Evaluate(instructions, ref ic, nodeIntervals, variableIntervals: variableIntervals);
      // the order of keys in a dictionary is guaranteed to be the same order as values in a dictionary
      // https://docs.microsoft.com/en-us/dotnet/api/system.collections.idictionary.keys?view=netcore-3.1#remarks
      //var box = variableIntervals.Values;
      //Box only contains intervals from multiple occurence variables
      var box = multipleOccurenceVariables.Select(k => variableIntervals[k]);
      if (minimization) {
        prioQ.Add(new BoxBound(box, interval.LowerBound));
        if (stopAtLimit && interval.LowerBound >= limit) stop = true;
      } else {
        prioQ.Add(new BoxBound(box, -interval.UpperBound));
        if (stopAtLimit && interval.UpperBound <= limit) stop = true;
      }

      var discardedBound = double.MaxValue;
      var runningBound = double.MaxValue;
      for (var depth = 0; depth < splittingIterations && prioQ.Count > 0 && !stop; ++depth) {
        var currentBound = prioQ.Min;
        prioQ.Remove(currentBound);

        if (currentBound.box.All(x => x.Width < splittingWidth)) {
          discardedBound = Math.Min(discardedBound, currentBound.bound);
          continue;
        }

        var newBoxes = Split(currentBound.box, splittingWidth);

        var innerBound = double.MaxValue;
        foreach (var newBox in newBoxes) {
          //var intervalEnum = newBox.GetEnumerator();
          //var keyEnum = readonlyRanges.Keys.GetEnumerator();
          //while (intervalEnum.MoveNext() & keyEnum.MoveNext()) {
          //  variableIntervals[keyEnum.Current] = intervalEnum.Current;
          //}
          //Set the splitted variables
          var intervalEnum = newBox.GetEnumerator();
          foreach (var key in multipleOccurenceVariables) {
            intervalEnum.MoveNext();
            variableIntervals[key] = intervalEnum.Current;
          }

          ic = 0;
          var res = Evaluate(instructions, ref ic, nodeIntervals,
            new ReadOnlyDictionary<string, Interval>(variableIntervals));

          var boxBound = new BoxBound(newBox, minimization ? res.LowerBound : -res.UpperBound);
          prioQ.Add(boxBound);
          innerBound = Math.Min(innerBound, boxBound.bound);
        }

        runningBound = innerBound;

        if (minimization) {
          if (stopAtLimit && innerBound >= limit)
            stop = true;
        } else {
          if (stopAtLimit && innerBound <= limit)
            stop = true;
        }
      }

      var bound = Math.Min(runningBound, discardedBound);
      if (bound == double.MaxValue)
        return minimization ? interval.LowerBound : interval.UpperBound;

      return minimization ? bound : -bound;
      //return minimization ? prioQ.First().bound : -prioQ.First().bound;
    }

    private static IEnumerable<IEnumerable<Interval>> Split(List<Interval> box) {
      var boxes = box.Select(region => region.Split())
                     .Select(split => new List<Interval> {split.Item1, split.Item2})
                     .ToList();

      return boxes.CartesianProduct();
    }

    private static IEnumerable<IEnumerable<Interval>> Split(List<Interval> box, double minWidth) {
      List<Interval> toList(Tuple<Interval, Interval> t) => new List<Interval> {t.Item1, t.Item2};
      var boxes = box.Select(region => region.Width > minWidth ? toList(region.Split()) : new List<Interval> {region})
                     .ToList();

      return boxes.CartesianProduct();
    }

    #endregion

    public Interval GetModelBound(ISymbolicExpressionTree tree, IntervalCollection variableRanges) {
      lock (syncRoot) {
        EvaluatedSolutions++;
      }

      var occuringVariableRanges = GetOccurringVariableRanges(tree, variableRanges);
      var instructions = PrepareInterpreterState(tree, occuringVariableRanges);
      Interval resultInterval;
      if (!UseIntervalSplitting) {
        var instructionCounter = 0;
        resultInterval = Evaluate(instructions, ref instructionCounter, variableIntervals: occuringVariableRanges);
      } else {
        var vars = ContainsVariableMultipleTimes(tree, out var variables);
        resultInterval = EvaluateWithSplitting(instructions, occuringVariableRanges, variables, SplittingIterations,
          SplittingWidth);
      }

      // because of numerical errors the bounds might be incorrect
      if (resultInterval.IsInfiniteOrUndefined || resultInterval.LowerBound <= resultInterval.UpperBound)
        return resultInterval;

      return new Interval(resultInterval.UpperBound, resultInterval.LowerBound);
    }

    public IDictionary<ISymbolicExpressionTreeNode, Interval> GetModelNodesBounds(
      ISymbolicExpressionTree tree, IntervalCollection variableRanges) {
      throw new NotImplementedException();
    }

    public double CheckConstraint(
      ISymbolicExpressionTree tree, IntervalCollection variableRanges, IntervalConstraint constraint) {
      var occuringVariableRanges = GetOccurringVariableRanges(tree, variableRanges);
      var instructions = PrepareInterpreterState(tree, occuringVariableRanges);
      if (!UseIntervalSplitting) {
        var instructionCounter = 0;
        var modelBound = Evaluate(instructions, ref instructionCounter, variableIntervals: occuringVariableRanges);
        if (constraint.Interval.Contains(modelBound)) return 0.0;


        var error = 0.0;

        if (!constraint.Interval.Contains(modelBound.LowerBound)) {
            error += Math.Abs(modelBound.LowerBound - constraint.Interval.LowerBound);
        }

        if (!constraint.Interval.Contains(modelBound.UpperBound)) {
            error += Math.Abs(modelBound.UpperBound - constraint.Interval.UpperBound);
        }

        return error;
        // return Math.Abs(modelBound.LowerBound - constraint.Interval.LowerBound) +
        //Math.Abs(modelBound.UpperBound - constraint.Interval.UpperBound);
      }

      if (double.IsNegativeInfinity(constraint.Interval.LowerBound) &&
          double.IsPositiveInfinity(constraint.Interval.UpperBound)) {
        return 0.0;
      }

      ContainsVariableMultipleTimes(tree, out var variables);

      var upperBound = 0.0;
      var lowerBound = 0.0;
      if (double.IsNegativeInfinity(constraint.Interval.LowerBound)) {
        upperBound = FindBound(instructions, occuringVariableRanges, variables, SplittingIterations, SplittingWidth,
          minimization: false, stopAtLimit: true, limit: constraint.Interval.UpperBound);

        return upperBound <= constraint.Interval.UpperBound
          ? 0.0
          : Math.Abs(upperBound - constraint.Interval.UpperBound);
      }

      if (double.IsPositiveInfinity(constraint.Interval.UpperBound)) {
        lowerBound = FindBound(instructions, occuringVariableRanges, variables, SplittingIterations, SplittingWidth,
          minimization: true, stopAtLimit: true, limit: constraint.Interval.LowerBound);

        return lowerBound >= constraint.Interval.LowerBound
          ? 0.0
          : Math.Abs(lowerBound - constraint.Interval.LowerBound);
      }

      upperBound = FindBound(instructions, occuringVariableRanges, variables, SplittingIterations, SplittingWidth,
        minimization: false, stopAtLimit: true, limit: constraint.Interval.UpperBound);
      lowerBound = FindBound(instructions, occuringVariableRanges, variables, SplittingIterations, SplittingWidth,
        minimization: true, stopAtLimit: true, limit: constraint.Interval.LowerBound);


      var res = 0.0;

      res += upperBound <= constraint.Interval.UpperBound ? 0.0 : Math.Abs(upperBound - constraint.Interval.UpperBound);
      res += lowerBound <= constraint.Interval.LowerBound ? 0.0 : Math.Abs(lowerBound - constraint.Interval.LowerBound);

      return res;
    }


    public bool IsCompatible(ISymbolicExpressionTree tree) {
      var containsUnknownSymbols = (
        from n in tree.Root.GetSubtree(0).IterateNodesPrefix()
        where
          !(n.Symbol is Variable) &&
          !(n.Symbol is Constant) &&
          !(n.Symbol is StartSymbol) &&
          !(n.Symbol is Addition) &&
          !(n.Symbol is Subtraction) &&
          !(n.Symbol is Multiplication) &&
          !(n.Symbol is Division) &&
          !(n.Symbol is Sine) &&
          !(n.Symbol is Cosine) &&
          !(n.Symbol is Tangent) &&
          !(n.Symbol is HyperbolicTangent) &&
          !(n.Symbol is Logarithm) &&
          !(n.Symbol is Exponential) &&
          !(n.Symbol is Square) &&
          !(n.Symbol is SquareRoot) &&
          !(n.Symbol is Cube) &&
          !(n.Symbol is CubeRoot) &&
          !(n.Symbol is Absolute) &&
          !(n.Symbol is AnalyticQuotient)
        select n).Any();
      return !containsUnknownSymbols;
    }
  }
}