Free cookie consent management tool by TermsFeed Policy Generator

Changeset 16727


Ignore:
Timestamp:
03/29/19 15:01:47 (5 years ago)
Author:
gkronber
Message:

#2994: added a unit test and made some minor improvements to interpreters

Location:
branches/2994-AutoDiffForIntervals
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/2994-AutoDiffForIntervals/HeuristicLab.Algorithms.DataAnalysis.ConstrainedNonlinearRegression/3.4/ConstrainedNonlinearRegression.cs

    r16696 r16727  
    248248        }
    249249      }
    250 
    251       // var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
    252       //
    253       // SymbolicRegressionConstantOptimizationEvaluator.OptimizeConstants(interpreter, tree, problemData, problemData.TrainingIndices,
    254       //   applyLinearScaling: applyLinearScaling, maxIterations: maxIterations,
    255       //   updateVariableWeights: false, updateConstantsInTree: true);
    256 
    257 
    258250      var intervals = problemData.IntervalConstraints;
    259251      var constraintsParser = new IntervalConstraintsParser();
    … …  
    315307      }
    316308
     309      // local function
    317310      void UpdateThetaValues(double[] theta) {
    318311        for (int i = 0; i < theta.Length; ++i) {
    … …  
    321314      }
    322315
     316      // buffers for calculate_jacobian
     317      var target = problemData.TargetVariableTrainingValues.ToArray();
     318      var fi_eval = new double[target.Length];
     319      var jac_eval = new double[target.Length, thetaValues.Count];
     320
    323321      // define the callback used by the alglib optimizer
    324322      // the x argument for this callback represents our theta
     323      // local function
    325324      void calculate_jacobian(double[] x, double[] fi, double[,] jac, object obj) {
    326325        UpdateThetaValues(x);
    … …  
    328327        var autoDiffEval = new VectorAutoDiffEvaluator();
    329328        autoDiffEval.Evaluate(preparedTree, problemData.Dataset, problemData.TrainingIndices.ToArray(),
    330           GetParameterNodes(preparedTree, allThetaNodes), out double[] fi_eval, out double[,] jac_eval);
    331         var target = problemData.TargetVariableTrainingValues.ToArray();
     329          GetParameterNodes(preparedTree, allThetaNodes), fi_eval, jac_eval);
    332330
    333331        // calc sum of squared errors and gradient
    … …  
    336334        for (int i = 0; i < target.Length; i++) {
    337335          var res = target[i] - fi_eval[i];
    338           sse += res * res;
     336          sse += 0.5 * res * res;
    339337          for (int j = 0; j < g.Length; j++) {
    340             g[j] += -2.0 * res * jac_eval[i, j];
    341           }
    342         }
    343 
    344         fi[0] = sse;
    345         for (int j = 0; j < x.Length; j++) { jac[0, j] = g[j]; }
     338            g[j] -= res * jac_eval[i, j];
     339          }
     340        }
     341
     342        fi[0] = sse / target.Length;
     343        for (int j = 0; j < x.Length; j++) { jac[0, j] = g[j] / target.Length; }
    346344
    347345        var intervalEvaluator = new IntervalEvaluator();
    … …  
    358356      }
    359357
     358
     359
    360360      // prepare alglib
    361361      alglib.minnlcstate state;
    362362      alglib.minnlcreport rep;
     363      alglib.optguardreport optGuardRep;
    363364      var x0 = thetaValues.ToArray();
    364365
    365366      alglib.minnlccreate(x0.Length, x0, out state);
    366       double epsx = 1e-6;
    367       int maxits = 0;
    368       alglib.minnlcsetalgoslp(state);
    369       alglib.minnlcsetcond(state, 0, maxits);
     367      alglib.minnlcsetalgoslp(state);        // SLP is more robust but slower
     368      alglib.minnlcsetcond(state, 0, maxIterations);
    370369      var s = Enumerable.Repeat(1d, x0.Length).ToArray();  // scale is set to unit scale
    371370      alglib.minnlcsetscale(state, s);
    372371
    373       // set boundary constraints
    374       // var boundaryLower = Enumerable.Repeat(-10d, n).ToArray();
    375       // var boundaryUpper = Enumerable.Repeat(10d, n).ToArray();
    376       // alglib.minnlcsetbc(state, boundaryLower, boundaryUpper);
    377372      // set non-linear constraints: 0 equality constraints, 1 inequality constraint
    378373      alglib.minnlcsetnlc(state, 0, constraintTrees.Count);
    379374
     375      alglib.minnlcoptguardsmoothness(state);
     376      alglib.minnlcoptguardgradient(state, 0.001);
     377
    380378      alglib.minnlcoptimize(state, calculate_jacobian, null, null);
    381379      alglib.minnlcresults(state, out double[] xOpt, out rep);
     380      alglib.minnlcoptguardresults(state, out optGuardRep);
    382381
    383382      var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
    … …  
    416415            var parent = n.Parent;
    417416            if(thetaNodes[thetaIdx].Any()) {
    418               // HACKY: REUSE CONSTANT TREE NODE IN SEVERAL TREES
     417              // HACK: REUSE CONSTANT TREE NODE IN SEVERAL TREES
    419418              // we use this trick to allow autodiff over thetas when thetas occurr multiple times in the tree (e.g. in derived trees)
    420419              var constNode = thetaNodes[thetaIdx].First();
    … …  
    444443      for (int i = 0; i < nodes.Count; ++i) {
    445444        var node = nodes[i];
    446         /*if (node is VariableTreeNode variableTreeNode) {
    447           var thetaVar = (VariableTreeNode)new Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
    448           thetaVar.Weight = 1;
    449           thetaVar.VariableName = $"θ{n++}";
    450 
    451           thetaNames.Add(thetaVar.VariableName);
    452           thetaValues.Add(variableTreeNode.Weight);
    453           variableTreeNode.Weight = 1; // set to unit weight
    454 
    455           var parent = variableTreeNode.Parent;
    456           var prod = MakeNode<Multiplication>(thetaVar, variableTreeNode);
    457           if (parent != null) {
    458             var index = parent.IndexOfSubtree(variableTreeNode);
    459             parent.RemoveSubtree(index);
    460             parent.InsertSubtree(index, prod);
    461           }
    462         } else*/ if (node is ConstantTreeNode constantTreeNode) {
     445        if (node is ConstantTreeNode constantTreeNode) {
    463446          var thetaVar = (VariableTreeNode)new Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
    464447          thetaVar.Weight = 1;

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

    r16695 r16727  
    7272              break;
    7373            }
     74          case OpCodes.SquareRoot: {
     75              instr.value.AssignIntRoot(code[c].value, 2);
     76              break;
     77            }
     78          case OpCodes.Cube: {
     79              instr.value.AssignIntPower(code[c].value, 3);
     80              break;
     81            }
     82          case OpCodes.CubeRoot: {
     83              instr.value.AssignIntRoot(code[c].value, 3);
     84              break;
     85            }
    7486          case OpCodes.Exp: {
    7587              instr.value.AssignExp(code[c].value);
    7688              break;
    7789            }
    78 
    7990          case OpCodes.Log: {
    8091              instr.value.AssignLog(code[c].value);
     92              break;
     93            }
     94          case OpCodes.Sin: {
     95              instr.value.AssignSin(code[c].value);
     96              break;
     97            }
     98          case OpCodes.Cos: {
     99              instr.value.AssignCos(code[c].value);
     100              break;
     101            }
     102          case OpCodes.Absolute: {
     103              instr.value.AssignAbs(code[c].value);
     104              break;
     105            }
     106          case OpCodes.AnalyticQuotient: {
     107              instr.value.Assign(code[c].value);
     108              for (int j = 1; j < n; ++j) {
     109                var t = instr.value.One;
     110                t.Add(code[c + j].value.Clone().IntPower(2));
     111                instr.value.Div(t.IntRoot(2));
     112              }
    81113              break;
    82114            }
    … …  
    221253    }
    222254
    223     public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, out double[] fi, out double[,] jac) {
     255    /// <summary>
     256    ///
     257    /// </summary>
     258    /// <param name="tree"></param>
     259    /// <param name="dataset"></param>
     260    /// <param name="rows"></param>
     261    /// <param name="parameterNodes"></param>
     262    /// <param name="fi">Function output. Must be allocated by the caller.</param>
     263    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
     264    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
    224265      if (cachedData == null || this.dataset != dataset) {
    225266        InitCache(dataset);
    … …  
    235276      var roundedTotal = rows.Length - remainingRows;
    236277
    237       fi = new double[rows.Length];
    238       jac = new double[rows.Length, nParams];
    239 
    240278      this.rows = rows;
    241279
    … …  
    247285        var g = code[0].value.Gradient;
    248286        for (int j = 0; j < nParams; ++j) {
    249           g.Elements[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
     287          if(g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
     288            v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
     289          } else {
     290            for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
     291          }
    250292        }
    251293      }
    … …  
    257299        var g = code[0].value.Gradient;
    258300        for (int j = 0; j < nParams; ++j)
    259           g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
     301          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
     302            v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
     303          } else {
     304            for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
     305          }
    260306      }
    261307    }
    … …  
    315361        var g = a.value.Gradient.Elements[paramIdx];
    316362        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
    317           g[i] = data[rows[i]];
     363          g[i - rowIndex] = data[rows[i]];
    318364        }
    319365      }
    … …  
    377423  public interface IAlgebraicType<T> {
    378424    T Zero { get; }
     425    T One { get; }
    379426
    380427    T AssignAbs(T a); // set this to assign abs(a)
    … …  
    429476    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    430477    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
     478    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     479    public AlgebraicDouble One => new AlgebraicDouble(1.0);
    431480    public AlgebraicDouble() { }
    432481    public AlgebraicDouble(double value) { this.Value = value; }
    … …  
    441490    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
    442491    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
    443     public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
     492    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = a.Value <= 0?double.NegativeInfinity : Math.Log(a.Value); return this; } // alternative definiton of log to prevent NaN
    444493    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
    445494    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
    446495    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
    447496    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
    448     public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; }
     497    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value)? double.NaN : Math.Sign(a.Value); return this; }
    449498    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
    450499
    … …  
    473522    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    474523    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
     524    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     525    public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)
    475526    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
    476527    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
    … …  
    496547    }
    497548
    498     public void AssignConstant(double constantValue) {
     549    public AlgebraicDoubleVector AssignConstant(double constantValue) {
    499550      for (int i = 0; i < arr.Length; ++i) {
    500551        arr[i] = constantValue;
    501552      }
     553      return this;
    502554    }
    503555
    … …  
    522574  }
    523575
     576
     577  /*
    524578  // vectors of algebraic types
    525   public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T> {
     579  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
    526580    private T[] elems;
    527581
    … …  
    561615    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
    562616
    563     public AlgebraicVector<T> Concat(AlgebraicVector<T> other) {
    564       var oldLen = Length;
    565       Array.Resize(ref this.elems, oldLen + other.Length);
    566       for (int i = oldLen; i < Length; i++) {
    567         elems[i] = other.elems[i - oldLen].Clone();
    568       }
    569       return this;
    570     }
    571617
    572618    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    573619    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
     620    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     621    public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }
    574622    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
    575623    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
    … …  
    591639  }
    592640
     641  */
     642
    593643
    594644  /// <summary>
    … …  
    625675    }
    626676
    627 
    628 
    629     // combined elements from both vectors
    630     private void UnionAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    631       // elements from a
    632       foreach (var kvp in a.elems) {
    633         // this = f(a, this)
    634         if (elems.TryGetValue(kvp.Key, out T value))
    635           mapAssign(kvp.Value, value);
    636         else {
    637           // this = f(a, 0)
    638           var newValue = kvp.Value.Zero;
    639           elems.Add(kvp.Key, newValue);
    640           mapAssign(kvp.Value, newValue);
    641         }
    642       }
    643       // elements from this (without a)
    644       foreach (var kvp in elems) {
    645         if (a.elems.ContainsKey(kvp.Key)) continue; // already processed above
    646                                                     // this = f(0, this)
    647         mapAssign(kvp.Value.Zero, kvp.Value);
    648       }
    649     }
    650 
    651     // keep only elements in both vectors
    652     private void IntersectAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    653       List<K> keysToRemove = new List<K>();
    654       foreach (var kvp in elems) {
    655         if (a.elems.TryGetValue(kvp.Key, out T value))
    656           mapAssign(value, kvp.Value);
    657         else
    658           keysToRemove.Add(kvp.Key);
    659       }
    660       foreach (var o in keysToRemove) elems.Remove(o); // -> zero
    661     }
    662 
    663677    // keep only elements from a
    664678    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
    … …  
    680694    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    681695    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
     696    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     697    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
    682698
    683699    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
    … …  
    685701
    686702    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
    687     public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Add(src)); return this; }
    688     public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Sub(src)); return this; }
    689     public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) { IntersectAssign(a, (src, dest) => dest.Mul(src)); return this; }
    690     public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Div(src)); return this; }
    691703    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
    692704    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
    … …  
    699711    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
    700712    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
     713    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
     714      foreach (var kvp in a.elems) {
     715        if (elems.TryGetValue(kvp.Key, out T value))
     716          value.Add(kvp.Value);
     717        else
     718          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
     719      }
     720      return this;
     721    }
     722
     723    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
     724      foreach (var kvp in a.elems) {
     725        if (elems.TryGetValue(kvp.Key, out T value))
     726          value.Sub(kvp.Value);
     727        else
     728          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
     729      }
     730      return this;
     731    }
     732
     733    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
     734      var keys = elems.Keys.ToArray();
     735      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
     736      foreach (var kvp in a.elems) {
     737        if (elems.TryGetValue(kvp.Key, out T value))
     738          value.Mul(kvp.Value); // this * a
     739      }
     740      return this;
     741    }
     742
     743    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
     744      return Mul(a.Clone().Inv());
     745    }
    701746
    702747    public AlgebraicSparseVector<K, T> Clone() {
    … …  
    734779    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    735780    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
     781    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
    736782    public AlgebraicInterval Add(AlgebraicInterval a) {
    737783      low.Add(a.low);
    … …  
    799845          if (a.Contains(0.0)) {
    800846            low = new MultivariateDual<AlgebraicDouble>(0.0);
    801             high = Algebraic.Max(low.Clone().IntPower(p), high.Clone().IntPower(p));
     847            high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
    802848          } else {
    803             var lowPower = low.Clone().IntPower(p);
    804             var highPower = high.Clone().IntPower(p);
     849            var lowPower = a.low.IntPower(p);
     850            var highPower = a.high.IntPower(p);
    805851            low = Algebraic.Min(lowPower, highPower);
    806852            high = Algebraic.Max(lowPower, highPower);
    … …  
    808854        } else {
    809855          // p is uneven
    810           var lowPower = low.Clone().IntPower(p);
    811           var highPower = high.Clone().IntPower(p);
     856          var lowPower = a.low.IntPower(p);
     857          var highPower = a.high.IntPower(p);
    812858          low = Algebraic.Min(lowPower, highPower);
    813859          high = Algebraic.Max(lowPower, highPower);
    … …  
    951997    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    952998    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
     999    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
     1000    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
    9531001
    9541002    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
    … …  
    10341082
    10351083    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
     1084    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
    10361085
    10371086    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }

  • branches/2994-AutoDiffForIntervals/Tests/AutoDiffTest.cs

    r16696 r16727  
    7676        var ds = new Dataset(vars, values);
    7777        var problemData = new RegressionProblemData(ds, vars, "f(x)");
    78         evaluator.Evaluate(t, ds, problemData.TrainingIndices.ToArray(), paramNodes, out double[] train, out double[,] trainJac);
     78        var train = new double[problemData.TrainingIndices.Count()];
     79        var trainJac = new double[train.Length, 2];
     80        evaluator.Evaluate(t, ds, problemData.TrainingIndices.ToArray(), paramNodes, train, trainJac);
    7981        Assert.AreEqual(2, train.Length);
    8082        Assert.AreEqual(3, train[0]);
    … …  
    8688        Assert.AreEqual(1, trainJac[1, 1]);
    8789
    88         evaluator.Evaluate(t, ds, problemData.TestIndices.ToArray(), paramNodes, out double[] test, out double[,] testJac);
     90        var test = new double[problemData.TestIndices.Count()];
     91        var testJac = new double[test.Length, 2];
     92        evaluator.Evaluate(t, ds, problemData.TestIndices.ToArray(), paramNodes, test, testJac);
    8993        Assert.AreEqual(3, test.Length);
    9094        Assert.AreEqual(5, test[0]);
    … …  
    124128      }
    125129
     130      {
     131        // as discussed with Fabrício
     132        var intervals = new Dictionary<string, Interval>();
     133        intervals.Add("x1", new Interval(60.0, 65.0));
     134        intervals.Add("x2", new Interval(30.0, 40.0));
     135        intervals.Add("x3", new Interval(5.0, 10.0));
     136        intervals.Add("x4", new Interval(0.5, 0.8));
     137        intervals.Add("x5", new Interval(0.2, 0.5));
     138
     139        var parser = new InfixExpressionParser();
     140
     141        var t1 = parser.Parse("x5/x4");
     142        var t2 = parser.Parse("log(x5/x4)");
     143        var t3 = parser.Parse("x3 * log(x5/x4)");
     144        var t4 = parser.Parse("x1*x2*x5");
     145        var t5 = parser.Parse("x4/x5");
     146        var t6 = parser.Parse("sqr(x4/x5)");
     147        var t7 = parser.Parse("(1 - sqr(x4/x5)) ");
     148        var t8 = parser.Parse("x1*x2*x5 *(1 - sqr(x4/x5))");
     149        var t9 = parser.Parse("x1*x2*x5 *(1 - sqr(x4/x5)) + x3 * log(x5/x4)");
     150
     151        var evaluator = new IntervalEvaluator();
     152        var result = evaluator.Evaluate(t1, intervals);
     153        Assert.AreEqual(0.25, result.LowerBound);
     154        Assert.AreEqual(1, result.UpperBound);
     155
     156        result = evaluator.Evaluate(t2, intervals);
     157        Assert.AreEqual(-1.386294361, result.LowerBound, 1e-6);
     158        Assert.AreEqual(0, result.UpperBound);
     159
     160        result = evaluator.Evaluate(t3, intervals);
     161        Assert.AreEqual(-13.86294361, result.LowerBound, 1e-6);
     162        Assert.AreEqual(0, result.UpperBound);
     163
     164        result = evaluator.Evaluate(t4, intervals);
     165        Assert.AreEqual(360, result.LowerBound);
     166        Assert.AreEqual(1300, result.UpperBound);
     167
     168        result = evaluator.Evaluate(t5, intervals);
     169        Assert.AreEqual(1, result.LowerBound, 1e-6);
     170        Assert.AreEqual(4, result.UpperBound);
     171
     172        result = evaluator.Evaluate(t6, intervals);
     173        Assert.AreEqual(1, result.LowerBound);
     174        Assert.AreEqual(16, result.UpperBound);
     175
     176        result = evaluator.Evaluate(t7, intervals);
     177        Assert.AreEqual(-15, result.LowerBound);
     178        Assert.AreEqual(0, result.UpperBound);
     179
     180        result = evaluator.Evaluate(t8, intervals);
     181        Assert.AreEqual(-19500, result.LowerBound);
     182        Assert.AreEqual(0, result.UpperBound);
     183
     184        result = evaluator.Evaluate(t9, intervals);
     185        Assert.AreEqual(-19513.86294, result.LowerBound, 1e-3);
     186        Assert.AreEqual(0, result.UpperBound);
     187
     188
     189      }
     190
    126191    }
    127192  }
Note: See TracChangeset for help on using the changeset viewer.