source: branches/3073_IA_constraint_splitting/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/IACompiledExpressionBoundsEstimator.cs @ 17886

using HEAL.Attic;

using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Parameters;

using System;
using System.Collections.Generic;
using System.Linq;
using System.Linq.Expressions;

namespace HeuristicLab.Problems.DataAnalysis.Symbolic
{
  [StorableType("60015D64-5D8B-408A-90A1-E4111BC114D4")]
  [Item("IA Compiled Expression Bounds Estimator", "Compile a symbolic model into a lambda and use it to evaluate model bounds.")]
  public class IACompiledExpressionBoundsEstimator : ParameterizedNamedItem, IBoundsEstimator
  {
    // interval method names
    private static readonly Dictionary<byte, string> methodName = new Dictionary<byte, string>() {
      { OpCodes.Add, "Add" },
      { OpCodes.Sub, "Subtract" },
      { OpCodes.Mul, "Multiply" },
      { OpCodes.Div, "Divide" },
      { OpCodes.Sin, "Sine" },
      { OpCodes.Cos, "Cosine" },
      { OpCodes.Tan, "Tangens" },
      { OpCodes.Tanh, "HyperbolicTangent" },
      { OpCodes.Log, "Logarithm" },
      { OpCodes.Exp, "Exponential" },
      { OpCodes.Square, "Square" },
      { OpCodes.Cube, "Cube" },
      { OpCodes.SquareRoot, "SquareRoot" },
      { OpCodes.CubeRoot, "CubicRoot" },
      { OpCodes.Absolute, "Absolute" },
      { OpCodes.AnalyticQuotient, "AnalyticalQuotient" },
    };

    private const string EvaluatedSolutionsParameterName = "EvaluatedSolutions";
    private const string UseIntervalSplittingParameterName = "Use Interval splitting";
    private const string MaxSplitParameterName = "MaxSplit";
    private const string MinWidthParameterName = "MinWidth";

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

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

    public IFixedValueParameter<IntValue> MaxSplitParameter {
      get => (IFixedValueParameter<IntValue>)Parameters[MaxSplitParameterName];
    }

    public IFixedValueParameter<DoubleValue> MinWidthParameter {
      get => (IFixedValueParameter<DoubleValue>)Parameters[MinWidthParameterName];
    }

    public int MaxSplit {
      get => MaxSplitParameter.Value.Value;
      set => MaxSplitParameter.Value.Value = value;
    }

    public double MinWidth {
      get => MinWidthParameter.Value.Value;
      set => MinWidthParameter.Value.Value = value;
    }

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

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

    private readonly object syncRoot = new object();

    public IACompiledExpressionBoundsEstimator() : base("IA Bounds Estimator",
      "Estimates the bounds of the model with interval arithmetic, by first compiling the model into a lambda.") {
      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>(MaxSplitParameterName,
        "Defines the number of iterations of splitting.", new IntValue(200)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(MinWidthParameterName,
        "Width of interval, after the splitting should stop.", new DoubleValue(0.0)));
    }

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

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

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



    public double CheckConstraint(ISymbolicExpressionTree tree, IntervalCollection variableRanges, IntervalConstraint constraint) {
      if (!UseIntervalSplitting) {
        var modelBound = GetModelBound(tree, variableRanges);
        if (constraint.Interval.Contains(modelBound)) return 0.0;
        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);

      lock (syncRoot) { EvaluatedSolutions++; }

      double upperBound;
      if (double.IsNegativeInfinity(constraint.Interval.LowerBound)) {
        upperBound = EstimateUpperBound(tree, variableRanges.GetReadonlyDictionary(), MaxSplit, MinWidth);

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

      double lowerBound;
      if (double.IsPositiveInfinity(constraint.Interval.UpperBound)) {
        lowerBound = EstimateLowerBound(tree, variableRanges.GetReadonlyDictionary(), MaxSplit, MinWidth);

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

      var ranges = variableRanges.GetReadonlyDictionary();
      lowerBound = EstimateLowerBound(tree, ranges, MaxSplit, MinWidth);
      upperBound = EstimateUpperBound(tree, ranges, MaxSplit, MinWidth);

      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 void ClearState() {
      EvaluatedSolutions = 0;
    }

    public Interval GetModelBound(ISymbolicExpressionTree tree, IntervalCollection variableRanges) {
      lock (syncRoot) { EvaluatedSolutions++; }
      var resultInterval = UseIntervalSplitting
        ? EstimateBounds(tree, variableRanges.GetReadonlyDictionary(), MaxSplit, MinWidth)
        : EstimateBounds(tree, variableRanges.GetReadonlyDictionary());

      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 NotSupportedException("Model nodes bounds are not supported.");
    }

    public void InitializeState() {
      EvaluatedSolutions = 0;
    }

    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;
    }

    #region compile a tree into a IA arithmetic lambda and estimate bounds
    static Expression MakeExpr(ISymbolicExpressionTreeNode node, IReadOnlyDictionary<string, Interval> variableRanges, IReadOnlyDictionary<string, int> variableIndices, Expression args) {
      Expression expr(ISymbolicExpressionTreeNode n) => MakeExpr(n, variableRanges, variableIndices, args);
      var opCode = OpCodes.MapSymbolToOpCode(node);

      switch (opCode) {
        case OpCodes.Variable: {
            var name = (node as VariableTreeNode).VariableName;
            var weight = (node as VariableTreeNode).Weight;
            var index = variableIndices[name];
            return Expression.Multiply(
              Expression.Constant(weight, typeof(double)),
              Expression.ArrayIndex(args, Expression.Constant(index, typeof(int)))
            );
          }
        case OpCodes.Constant: {
            var v = (node as ConstantTreeNode).Value;
            // we have to make an interval out of the constant because this may be the root of the tree (and we are expected to return an Interval)
            return Expression.Constant(new Interval(v, v), typeof(Interval));
          }
        case OpCodes.Add: {
            var e = expr(node.GetSubtree(0));
            foreach (var s in node.Subtrees.Skip(1)) {
              e = Expression.Add(e, expr(s));
            }
            return e;
          }
        case OpCodes.Sub: {
            var e = expr(node.GetSubtree(0));
            if (node.SubtreeCount == 1) {
              return Expression.Subtract(Expression.Constant(0.0, typeof(double)), e);
            }
            foreach (var s in node.Subtrees.Skip(1)) {
              e = Expression.Subtract(e, expr(s));
            }
            return e;
          }
        case OpCodes.Mul: {
            var e = expr(node.GetSubtree(0));
            foreach (var s in node.Subtrees.Skip(1)) {
              e = Expression.Multiply(e, expr(s));
            }
            return e;
          }
        case OpCodes.Div: {
            var e1 = expr(node.GetSubtree(0));
            if (node.SubtreeCount == 1) {
              return Expression.Divide(Expression.Constant(1.0, typeof(double)), e1);
            }
            // division is more expensive than multiplication so we use this construct
            var e2 = expr(node.GetSubtree(1));
            foreach (var s in node.Subtrees.Skip(2)) {
              e2 = Expression.Multiply(e2, expr(s));
            }
            return Expression.Divide(e1, e2);
          }
        // all these cases share the same code: get method info by name, emit call expression
        case OpCodes.Exp:
        case OpCodes.Log:
        case OpCodes.Sin:
        case OpCodes.Cos:
        case OpCodes.Tan:
        case OpCodes.Tanh:
        case OpCodes.Square:
        case OpCodes.Cube:
        case OpCodes.SquareRoot:
        case OpCodes.CubeRoot:
        case OpCodes.Absolute:
        case OpCodes.AnalyticQuotient: {
            var arg = expr(node.GetSubtree(0));
            var fun = typeof(Interval).GetMethod(methodName[opCode], new[] { arg.Type });
            return Expression.Call(fun, arg);
          }
        default: {
            throw new Exception($"Unsupported OpCode {opCode} encountered.");
          }
      }
    }

    public static IReadOnlyDictionary<string, int> GetVariableIndices(ISymbolicExpressionTree tree, IReadOnlyDictionary<string, Interval> variableIntervals, out Interval[] inputIntervals) {
      var variableIndices = new Dictionary<string, int>();
      var root = tree.Root;
      while (root.Symbol is ProgramRootSymbol || root.Symbol is StartSymbol) {
        root = root.GetSubtree(0);
      }
      inputIntervals = new Interval[variableIntervals.Count];
      int count = 0;
      foreach (var node in root.IterateNodesPrefix()) {
        if (node is VariableTreeNode varNode) {
          var name = varNode.VariableName;
          if (!variableIndices.ContainsKey(name)) {
            variableIndices[name] = count;
            inputIntervals[count] = variableIntervals[name];
            ++count;
          }
        }
      }
      Array.Resize<Interval>(ref inputIntervals, count);
      return variableIndices;
    }

    public static Func<Interval[], Interval> Compile(ISymbolicExpressionTree tree, IReadOnlyDictionary<string, Interval> variableRanges, IReadOnlyDictionary<string, int> variableIndices) {
      var root = tree.Root.GetSubtree(0).GetSubtree(0);
      var args = Expression.Parameter(typeof(Interval[]));
      var expr = MakeExpr(root, variableRanges, variableIndices, args);
      return Expression.Lambda<Func<Interval[], Interval>>(expr, args).Compile();
    }

    public static Interval EstimateBounds(ISymbolicExpressionTree tree, IReadOnlyDictionary<string, Interval> variableRanges, int n = 0, double w = 1e-5) {
      var variableIndices = GetVariableIndices(tree, variableRanges, out Interval[] x);
      var f = Compile(tree, variableRanges, variableIndices);
      if (n == 0) return f(x);
      var inf = EstimateBound(x, f, true, n, w);
      var sup = EstimateBound(x, f, false, n, w);
      return inf < sup ? new Interval(inf, sup) : new Interval(sup, inf);
    }
    public  double EstimateLowerBound(ISymbolicExpressionTree tree, IReadOnlyDictionary<string, Interval> variableRanges, int n = 1000, double w = 1e-5) {
      var variableIndices = GetVariableIndices(tree, variableRanges, out Interval[] x);
      var f = Compile(tree, variableRanges, variableIndices);
      var inf = EstimateBound(x, f, true, n, w);
      return inf;
    }

    public double EstimateUpperBound(ISymbolicExpressionTree tree, IReadOnlyDictionary<string, Interval> variableRanges, int n = 1000, double w = 1e-5) {
      var variableIndices = GetVariableIndices(tree, variableRanges, out Interval[] x);
      var f = Compile(tree, variableRanges, variableIndices);
      var sup = EstimateBound(x, f, false, n, w);
      return sup;
    }

    static double EstimateBound(Interval[] x, Func<Interval[], Interval> f, bool m = false, int n = 1000, double w = 1e-4) {
      double getBound(Interval iv) => m ? iv.LowerBound : -iv.UpperBound;
      double getVolume(Interval[] box) => box.Aggregate(1.0, (acc, iv) => acc * iv.Width);

      var splits = Enumerable.Range(0, x.Length).Select(_ => new List<Interval>()).ToArray();
      var newbox = new Interval[x.Length];

      int compare(Tuple<double, double, Interval[]> a, Tuple<double, double, Interval[]> b) {
        var res = a.Item1.CompareTo(b.Item1);
        if (res == 0) {
          res = b.Item2.CompareTo(a.Item2);
        }
        return res;
      }

      var q = new SortedSet<Tuple<double, double, Interval[]>>(Comparer<Tuple<double, double, Interval[]>>.Create(compare)) {
        Tuple.Create(getBound(f(x)), getVolume(x), x)
      };


      var bestBound = double.MaxValue;

      // examine all the ordered pairs in the cartesian product
      void next_pair(int i) {
        if (i == splits.Length) {
          var tmp = newbox.ToArray(); // make a copy to put in the queue
          q.Add(Tuple.Create(getBound(f(tmp)), getVolume(tmp), tmp));
          return;
        }

        foreach (var iv in splits[i]) {
          newbox[i] = iv;
          next_pair(i + 1);
        }
      }

      while (q.Count > 0 && n-- > 0) {
        var currentBound = q.Min; q.Remove(currentBound);
        var bound = currentBound.Item1;
        var box = currentBound.Item3;

        if (!box.Any(b => b.Width > w)) {
          bestBound = Math.Min(bestBound, bound);
          continue;
        }

        // do the splits
        for (int i = 0; i < box.Length; ++i) {
          splits[i].Clear();
          var iv = box[i];
          if (iv.Width > w) {
            var t = iv.Split();
            splits[i].Add(t.Item1);
            splits[i].Add(t.Item2);
          } else {
            splits[i].Add(iv);
          }
        }
        next_pair(0);
      }
      if (q.Count > 0) {
        bestBound = Math.Min(bestBound, q.First().Item1);
      }
      return m ? bestBound : -bestBound;
    }
    #endregion
  }
}