Changeset 17328


Ignore:
Timestamp:
10/11/19 20:55:22 (3 years ago)
Author:
gkronber
Message:

#2994 made some fixes in the const-opt evaluator with constraints and added some debugging capabilities

Location:
branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs

    r17325 r17328  
    4141    private double[] bestConstraintValues;
    4242    public double[] BestConstraintValues => bestConstraintValues;
     43
     44    public NLOpt.nlopt_result OptResult { get; private set; }
    4345
    4446    private bool disposed = false;
     
    8486    private readonly VectorEvaluator eval;
    8587    private readonly bool invalidProblem = false;
     88    private readonly double targetVariance;
    8689
    8790    // end internal state
     
    112115      // buffers
    113116      target = problemData.TargetVariableTrainingValues.ToArray();
    114       var targetStDev = target.StandardDeviationPop();
    115       var targetVariance = targetStDev * targetStDev;
    116       var targetMean = target.Average();
     117      targetVariance = target.VariancePop();
    117118      var pred = interpreter.GetSymbolicExpressionTreeValues(expr, problemData.Dataset, trainingRows).ToArray();
    118119
    119       bestError = targetVariance;
    120 
    121       if (pred.Any(pi => double.IsInfinity(pi) || double.IsNaN(pi))) {
     120      // all trees are linearly scaled (to improve GP performance)
     121      if (pred.Any(pi => double.IsInfinity(pi) || double.IsNaN(pi)) || pred.StandardDeviationPop() == 0) {
    122122        invalidProblem = true;
    123       }
    124 
    125       // all trees are linearly scaled (to improve GP performance)
    126       #region linear scaling
    127       var predStDev = pred.StandardDeviationPop();
    128       if (predStDev == 0) {
    129         invalidProblem = true;
    130       }
    131       var predMean = pred.Average();
    132 
    133       var scalingFactor = targetStDev / predStDev;
    134       var offset = targetMean - predMean * scalingFactor;
    135 
    136       scaledTree = CopyAndScaleTree(expr, scalingFactor, offset);
    137       #endregion
    138 
    139       // convert constants to variables named theta...
    140       var treeForDerivation = ReplaceAndExtractParameters(scaledTree, out List<string> thetaNames, out thetaValues); // copies the tree
    141 
    142       // create trees for relevant derivatives
    143       Dictionary<string, ISymbolicExpressionTree> derivatives = new Dictionary<string, ISymbolicExpressionTree>();
    144       allThetaNodes = thetaNames.Select(_ => new List<ConstantTreeNode>()).ToArray();
    145       constraintTrees = new List<ISymbolicExpressionTree>();
    146       foreach (var constraint in intervalConstraints.Constraints) {
    147         if (!constraint.Enabled) continue;
    148         if (constraint.IsDerivation) {
    149           if (!problemData.AllowedInputVariables.Contains(constraint.Variable))
    150             throw new ArgumentException($"Invalid constraint: the variable {constraint.Variable} does not exist in the dataset.");
    151           var df = DerivativeCalculator.Derive(treeForDerivation, constraint.Variable);
    152 
    153           // NLOpt requires constraint expressions of the form c(x) <= 0
    154           // -> we make two expressions, one for the lower bound and one for the upper bound
    155 
    156           if (constraint.Interval.UpperBound < double.PositiveInfinity) {
    157             var df_smaller_upper = Subtract((ISymbolicExpressionTree)df.Clone(), CreateConstant(constraint.Interval.UpperBound));
    158             // convert variables named theta back to constants
    159             var df_prepared = ReplaceVarWithConst(df_smaller_upper, thetaNames, thetaValues, allThetaNodes);
    160             constraintTrees.Add(df_prepared);
    161           }
    162           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
    163             var df_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)df.Clone());
    164             // convert variables named theta back to constants
    165             var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);
    166             constraintTrees.Add(df_prepared);
    167           }
     123        bestError = targetVariance;
     124        bestSolution = new double[0];
     125        bestConstraintValues = new double[0];
     126      } else {
     127        #region linear scaling
     128        var cov = OnlineCovarianceCalculator.Calculate(pred, target, out OnlineCalculatorError covError);
     129        var scalingFactor = cov / pred.VariancePop();
     130        var offset = target.Average() - scalingFactor * pred.Average();
     131
     132        scaledTree = CopyAndScaleTree(expr, scalingFactor, offset);
     133        #endregion
     134
     135        // convert constants to variables named theta...
     136        var treeForDerivation = ReplaceAndExtractParameters(scaledTree, out List<string> thetaNames, out thetaValues); // copies the tree
     137
     138        // create trees for relevant derivatives
     139        Dictionary<string, ISymbolicExpressionTree> derivatives = new Dictionary<string, ISymbolicExpressionTree>();
     140        allThetaNodes = thetaNames.Select(_ => new List<ConstantTreeNode>()).ToArray();
     141        constraintTrees = new List<ISymbolicExpressionTree>();
     142        foreach (var constraint in intervalConstraints.Constraints) {
     143          if (!constraint.Enabled) continue;
     144          if (constraint.IsDerivation) {
     145            if (!problemData.AllowedInputVariables.Contains(constraint.Variable))
     146              throw new ArgumentException($"Invalid constraint: the variable {constraint.Variable} does not exist in the dataset.");
     147            var df = DerivativeCalculator.Derive(treeForDerivation, constraint.Variable);
     148
     149            // NLOpt requires constraint expressions of the form c(x) <= 0
     150            // -> we make two expressions, one for the lower bound and one for the upper bound
     151
     152            if (constraint.Interval.UpperBound < double.PositiveInfinity) {
     153              var df_smaller_upper = Subtract((ISymbolicExpressionTree)df.Clone(), CreateConstant(constraint.Interval.UpperBound));
     154              // convert variables named theta back to constants
     155              var df_prepared = ReplaceVarWithConst(df_smaller_upper, thetaNames, thetaValues, allThetaNodes);
     156              constraintTrees.Add(df_prepared);
     157            }
     158            if (constraint.Interval.LowerBound > double.NegativeInfinity) {
     159              var df_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)df.Clone());
     160              // convert variables named theta back to constants
     161              var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);
     162              constraintTrees.Add(df_prepared);
     163            }
     164          } else {
     165            if (constraint.Interval.UpperBound < double.PositiveInfinity) {
     166              var f_smaller_upper = Subtract((ISymbolicExpressionTree)treeForDerivation.Clone(), CreateConstant(constraint.Interval.UpperBound));
     167              // convert variables named theta back to constants
     168              var df_prepared = ReplaceVarWithConst(f_smaller_upper, thetaNames, thetaValues, allThetaNodes);
     169              constraintTrees.Add(df_prepared);
     170            }
     171            if (constraint.Interval.LowerBound > double.NegativeInfinity) {
     172              var f_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)treeForDerivation.Clone());
     173              // convert variables named theta back to constants
     174              var df_prepared = ReplaceVarWithConst(f_larger_lower, thetaNames, thetaValues, allThetaNodes);
     175              constraintTrees.Add(df_prepared);
     176            }
     177          }
     178        }
     179
     180        preparedTree = ReplaceVarWithConst(treeForDerivation, thetaNames, thetaValues, allThetaNodes);
     181        preparedTreeParameterNodes = GetParameterNodes(preparedTree, allThetaNodes);
     182
     183        var dim = thetaValues.Count;
     184        fi_eval = new double[target.Length]; // init buffer;
     185        jac_eval = new double[target.Length, dim]; // init buffer
     186
     187
     188        var minVal = Math.Min(-1000.0, thetaValues.Min());
     189        var maxVal = Math.Max(1000.0, thetaValues.Max());
     190        var lb = Enumerable.Repeat(minVal, thetaValues.Count).ToArray();
     191        var up = Enumerable.Repeat(maxVal, thetaValues.Count).ToArray();
     192        nlopt = NLOpt.nlopt_create(GetSolver(solver), (uint)dim);
     193
     194        NLOpt.nlopt_set_lower_bounds(nlopt, lb);
     195        NLOpt.nlopt_set_upper_bounds(nlopt, up);
     196        calculateObjectiveDelegate = new NLOpt.nlopt_func(CalculateObjective); // keep a reference to the delegate (see below)
     197        NLOpt.nlopt_set_min_objective(nlopt, calculateObjectiveDelegate, IntPtr.Zero); // --> without preconditioning
     198
     199        //preconditionDelegate = new NLOpt.nlopt_precond(PreconditionObjective);
     200        //NLOpt.nlopt_set_precond_min_objective(nlopt, calculateObjectiveDelegate, preconditionDelegate, IntPtr.Zero);
     201
     202
     203        constraintDataPtr = new IntPtr[constraintTrees.Count];
     204        calculateConstraintDelegates = new NLOpt.nlopt_func[constraintTrees.Count]; // make sure we keep a reference to the delegates (otherwise GC will free delegate objects see https://stackoverflow.com/questions/7302045/callback-delegates-being-collected#7302258)
     205        for (int i = 0; i < constraintTrees.Count; i++) {
     206          var constraintData = new ConstraintData() { Idx = i, Tree = constraintTrees[i], ParameterNodes = GetParameterNodes(constraintTrees[i], allThetaNodes) };
     207          constraintDataPtr[i] = Marshal.AllocHGlobal(Marshal.SizeOf<ConstraintData>());
     208          Marshal.StructureToPtr(constraintData, constraintDataPtr[i], fDeleteOld: false);
     209          calculateConstraintDelegates[i] = new NLOpt.nlopt_func(CalculateConstraint);
     210          NLOpt.nlopt_add_inequality_constraint(nlopt, calculateConstraintDelegates[i], constraintDataPtr[i], 1e-8);
     211          // NLOpt.nlopt_add_precond_inequality_constraint(nlopt, calculateConstraintDelegates[i], preconditionDelegate, constraintDataPtr[i], 1e-8);
     212        }
     213
     214
     215        var x = thetaValues.ToArray();  /* initial parameters */
     216                                        // calculate initial quality
     217        // calculate constraints of initial solution
     218        double[] constraintGrad = new double[x.Length];
     219        bestConstraintValues = new double[calculateConstraintDelegates.Length];
     220        for (int i = 0; i < calculateConstraintDelegates.Length; i++) {
     221          bestConstraintValues[i] = calculateConstraintDelegates[i].Invoke((uint)x.Length, x, constraintGrad, constraintDataPtr[i]);
     222        }
     223
     224        // if all constraints are OK then calculate the initial error (when there is any violation we use the variance)
     225        if (bestConstraintValues.Any(c => c > 0)) {
     226          bestError = targetVariance;
    168227        } else {
    169           if (constraint.Interval.UpperBound < double.PositiveInfinity) {
    170             var f_smaller_upper = Subtract((ISymbolicExpressionTree)treeForDerivation.Clone(), CreateConstant(constraint.Interval.UpperBound));
    171             // convert variables named theta back to constants
    172             var df_prepared = ReplaceVarWithConst(f_smaller_upper, thetaNames, thetaValues, allThetaNodes);
    173             constraintTrees.Add(df_prepared);
    174           }
    175           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
    176             var f_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)treeForDerivation.Clone());
    177             // convert variables named theta back to constants
    178             var df_prepared = ReplaceVarWithConst(f_larger_lower, thetaNames, thetaValues, allThetaNodes);
    179             constraintTrees.Add(df_prepared);
    180           }
    181         }
    182       }
    183 
    184       preparedTree = ReplaceVarWithConst(treeForDerivation, thetaNames, thetaValues, allThetaNodes);
    185       preparedTreeParameterNodes = GetParameterNodes(preparedTree, allThetaNodes);
    186 
    187       var dim = thetaValues.Count;
    188       fi_eval = new double[target.Length]; // init buffer;
    189       jac_eval = new double[target.Length, dim]; // init buffer
    190 
    191 
    192       var minVal = Math.Min(-1000.0, thetaValues.Min());
    193       var maxVal = Math.Max(1000.0, thetaValues.Max());
    194       var lb = Enumerable.Repeat(minVal, thetaValues.Count).ToArray();
    195       var up = Enumerable.Repeat(maxVal, thetaValues.Count).ToArray();
    196       nlopt = NLOpt.nlopt_create(GetSolver(solver), (uint)dim);
    197 
    198       NLOpt.nlopt_set_lower_bounds(nlopt, lb);
    199       NLOpt.nlopt_set_upper_bounds(nlopt, up);
    200       calculateObjectiveDelegate = new NLOpt.nlopt_func(CalculateObjective); // keep a reference to the delegate (see below)
    201       NLOpt.nlopt_set_min_objective(nlopt, calculateObjectiveDelegate, IntPtr.Zero); // --> without preconditioning
    202 
    203       //preconditionDelegate = new NLOpt.nlopt_precond(PreconditionObjective);
    204       //NLOpt.nlopt_set_precond_min_objective(nlopt, calculateObjectiveDelegate, preconditionDelegate, IntPtr.Zero);
    205 
    206 
    207       constraintDataPtr = new IntPtr[constraintTrees.Count];
    208       calculateConstraintDelegates = new NLOpt.nlopt_func[constraintTrees.Count]; // make sure we keep a reference to the delegates (otherwise GC will free delegate objects see https://stackoverflow.com/questions/7302045/callback-delegates-being-collected#7302258)
    209       for (int i = 0; i < constraintTrees.Count; i++) {
    210         var constraintData = new ConstraintData() { Idx = i, Tree = constraintTrees[i], ParameterNodes = GetParameterNodes(constraintTrees[i], allThetaNodes) };
    211         constraintDataPtr[i] = Marshal.AllocHGlobal(Marshal.SizeOf<ConstraintData>());
    212         Marshal.StructureToPtr(constraintData, constraintDataPtr[i], fDeleteOld: false);
    213         calculateConstraintDelegates[i] = new NLOpt.nlopt_func(CalculateConstraint);
    214         NLOpt.nlopt_add_inequality_constraint(nlopt, calculateConstraintDelegates[i], constraintDataPtr[i], 1e-8);
    215         // NLOpt.nlopt_add_precond_inequality_constraint(nlopt, calculateConstraintDelegates[i], preconditionDelegate, constraintDataPtr[i], 1e-8);
    216       }
    217 
    218       NLOpt.nlopt_set_ftol_rel(nlopt, ftol_rel);
    219       NLOpt.nlopt_set_ftol_abs(nlopt, ftol_abs);
    220       NLOpt.nlopt_set_maxtime(nlopt, maxTime);
    221       NLOpt.nlopt_set_maxeval(nlopt, maxIterations);
     228          bestError = CalculateObjective((uint)dim, x, null, IntPtr.Zero);
     229        }
     230
     231        NLOpt.nlopt_set_ftol_rel(nlopt, ftol_rel);
     232        NLOpt.nlopt_set_ftol_abs(nlopt, ftol_abs);
     233        NLOpt.nlopt_set_maxtime(nlopt, maxTime);
     234        NLOpt.nlopt_set_maxeval(nlopt, maxIterations);
     235      } // end if valid problem
    222236    }
    223237
     
    234248      var x = thetaValues.ToArray();  /* initial guess */
    235249      double minf = double.MaxValue; /* minimum objective value upon return */
    236       var res = NLOpt.nlopt_optimize(nlopt, x, ref minf);
    237 
    238       if (res < 0 && res != NLOpt.nlopt_result.NLOPT_FORCED_STOP) {
     250      OptResult = NLOpt.nlopt_optimize(nlopt, x, ref minf);
     251
     252      if (OptResult < 0 && OptResult != NLOpt.nlopt_result.NLOPT_FORCED_STOP) {
    239253        // throw new InvalidOperationException($"NLOpt failed {res} {NLOpt.nlopt_get_errmsg(nlopt)}");
    240254        return;
    241       } else /*if ( minf <= bestError ) */{
    242         bestSolution = x;
    243         bestError = minf;
     255      } else if (!double.IsNaN(minf) && minf < double.MaxValue) {
    244256
    245257        // calculate constraints of final solution
    246258        double[] _ = new double[x.Length];
    247         bestConstraintValues = new double[calculateConstraintDelegates.Length];
     259        var optimizedConstraintValues = new double[calculateConstraintDelegates.Length];
    248260        for (int i = 0; i < calculateConstraintDelegates.Length; i++) {
    249           bestConstraintValues[i] = calculateConstraintDelegates[i].Invoke((uint)x.Length, x, _, constraintDataPtr[i]);
    250         }
    251 
    252         // update parameters in tree
    253         UpdateParametersInTree(scaledTree, x);
    254 
    255         if (mode == OptimizationMode.UpdateParameters) {
    256           // update original expression (when called from evaluator we want to write back optimized parameters)
    257           expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
    258           expr.Root.GetSubtree(0).InsertSubtree(0,
    259             scaledTree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(0).GetSubtree(0) // insert the optimized sub-tree (without scaling nodes)
    260             );
    261         } else if (mode == OptimizationMode.UpdateParametersAndKeepLinearScaling) {
    262           expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
    263           expr.Root.GetSubtree(0).InsertSubtree(0, scaledTree.Root.GetSubtree(0).GetSubtree(0)); // insert the optimized sub-tree (including scaling nodes)
     261          optimizedConstraintValues[i] = calculateConstraintDelegates[i].Invoke((uint)x.Length, x, _, constraintDataPtr[i]);
     262        }
     263
     264        // we accept the optimized parameters when either the error has been reduced or at least one of the violated constraints was improved (all constraints < 0 are OK anyway)
     265        if (minf < bestError || bestConstraintValues.Zip(optimizedConstraintValues, Tuple.Create).Any(t => t.Item1 > 0 && t.Item1 > t.Item2)) {
     266
     267          bestConstraintValues = optimizedConstraintValues;
     268          bestSolution = x;
     269          bestError = Math.Min(minf, targetVariance);  // limit the error to variance
     270
     271          // update parameters in tree
     272          UpdateParametersInTree(scaledTree, x);
     273
     274          if (mode == OptimizationMode.UpdateParameters) {
     275            // update original expression (when called from evaluator we want to write back optimized parameters)
     276            var newChild = (ISymbolicExpressionTreeNode)scaledTree.Root.GetSubtree(0).GetSubtree(0).GetSubtree(0).GetSubtree(0).Clone(); // the optimized sub-tree (without scaling nodes), we need to clone again to remove the parameter nodes which are used in multiple trees (and have incorrect parent relations)
     277            var oldChild = expr.Root.GetSubtree(0).GetSubtree(0);
     278            expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
     279            expr.Root.GetSubtree(0).InsertSubtree(0, newChild);
     280
     281          } else if (mode == OptimizationMode.UpdateParametersAndKeepLinearScaling) {
     282            var oldChild = expr.Root.GetSubtree(0).GetSubtree(0);
     283            var newChild = (ISymbolicExpressionTreeNode)scaledTree.Root.GetSubtree(0).GetSubtree(0).Clone(); //  we need to clone again to remove the parameter nodes which are used in multiple trees(and have incorrect parent relations)
     284            expr.Root.GetSubtree(0).RemoveSubtree(0); // delete old tree
     285            expr.Root.GetSubtree(0).InsertSubtree(0, newChild); // insert the optimized sub-tree (including scaling nodes)
     286          }
    264287        }
    265288      }
     
    323346      }
    324347
    325       UpdateBestSolution(sse / target.Length, curX);
     348      // UpdateBestSolution(sse / target.Length, curX);
    326349      RaiseFunctionEvaluated();
    327350
     
    401424        out double[] lowerGradient, out double[] upperGradient);
    402425
    403       var refInterval = refIntervalEvaluator.GetSymbolicExpressionTreeInterval(constraintData.Tree, dataIntervals);
    404       if (Math.Abs(interval.LowerBound - refInterval.LowerBound) > Math.Abs(interval.LowerBound) * 1e-4) throw new InvalidProgramException($"Intervals don't match. {interval.LowerBound} <> {refInterval.LowerBound}");
    405       if (Math.Abs(interval.UpperBound - refInterval.UpperBound) > Math.Abs(interval.UpperBound) * 1e-4) throw new InvalidProgramException($"Intervals don't match. {interval.UpperBound} <> {refInterval.UpperBound}");
     426      // compare to reference interval calculation
     427      // var refInterval = refIntervalEvaluator.GetSymbolicExpressionTreeInterval(constraintData.Tree, dataIntervals);
     428      // if (Math.Abs(interval.LowerBound - refInterval.LowerBound) > Math.Abs(interval.LowerBound) * 1e-4) throw new InvalidProgramException($"Intervals don't match. {interval.LowerBound} <> {refInterval.LowerBound}");
     429      // if (Math.Abs(interval.UpperBound - refInterval.UpperBound) > Math.Abs(interval.UpperBound) * 1e-4) throw new InvalidProgramException($"Intervals don't match. {interval.UpperBound} <> {refInterval.UpperBound}");
    406430
    407431      // we transformed this to a constraint c(x) <= 0, so only the upper bound is relevant for us
     
    669693        for (int i = 0; i < constraintDataPtr.Length; i++)
    670694          if (constraintDataPtr[i] != IntPtr.Zero) {
     695            Marshal.DestroyStructure<ConstraintData>(constraintDataPtr[i]);
    671696            Marshal.FreeHGlobal(constraintDataPtr[i]);
    672697            constraintDataPtr[i] = IntPtr.Zero;
  • branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/NLOptEvaluator.cs

    r17325 r17328  
    3131using HEAL.Attic;
    3232using System.Runtime.InteropServices;
     33using System.Diagnostics;
    3334
    3435namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Regression {
     
    4748    private const string CountEvaluationsParameterName = "Count Function and Gradient Evaluations";
    4849
     50
     51    private const string AchievedQualityImprovementParameterName = "AchievedQualityImprovment";
     52    private const string NumberOfConstraintViolationsBeforeOptimizationParameterName = "NumberOfConstraintViolationsBeforeOptimization";
     53    private const string NumberOfConstraintViolationsAfterOptimizationParameterName = "NumberOfConstraintViolationsAfterOptimization";
     54    private const string ConstraintsBeforeOptimizationParameterName = "ConstraintsBeforeOptimization";
     55    private const string ViolationsAfterOptimizationParameterName = "ConstraintsAfterOptimization";
     56    private const string OptimizationDurationParameterName = "OptimizationDuration";
     57
    4958    public IFixedValueParameter<IntValue> ConstantOptimizationIterationsParameter {
    5059      get { return (IFixedValueParameter<IntValue>)Parameters[ConstantOptimizationIterationsParameterName]; }
     
    7786    public IConstrainedValueParameter<StringValue> SolverParameter {
    7887      get { return (IConstrainedValueParameter<StringValue>)Parameters["Solver"]; }
     88    }
     89
     90
     91    public ILookupParameter<DoubleValue> AchievedQualityImprovementParameter {
     92      get { return (ILookupParameter<DoubleValue>)Parameters[AchievedQualityImprovementParameterName]; }
     93    }
     94    public ILookupParameter<DoubleValue> NumberOfConstraintViolationsBeforeOptimizationParameter {
     95      get { return (ILookupParameter<DoubleValue>)Parameters[NumberOfConstraintViolationsBeforeOptimizationParameterName]; }
     96    }
     97    public ILookupParameter<DoubleValue> NumberOfConstraintViolationsAfterOptimizationParameter {
     98      get { return (ILookupParameter<DoubleValue>)Parameters[NumberOfConstraintViolationsAfterOptimizationParameterName]; }
     99    }
     100    public ILookupParameter<DoubleArray> ViolationsAfterOptimizationParameter {
     101      get { return (ILookupParameter<DoubleArray>)Parameters[ViolationsAfterOptimizationParameterName]; }
     102    }
     103    public ILookupParameter<DoubleArray> ViolationsBeforeOptimizationParameter {
     104      get { return (ILookupParameter<DoubleArray>)Parameters[ConstraintsBeforeOptimizationParameterName]; }
     105    }
     106    public ILookupParameter<DoubleValue> OptimizationDurationParameter {
     107      get { return (ILookupParameter<DoubleValue>)Parameters[OptimizationDurationParameterName]; }
    79108    }
    80109
     
    134163      Parameters.Add(new ResultParameter<IntValue>(FunctionEvaluationsResultParameterName, "The number of function evaluations performed by the constants optimization evaluator", "Results", new IntValue()));
    135164      Parameters.Add(new ResultParameter<IntValue>(GradientEvaluationsResultParameterName, "The number of gradient evaluations performed by the constants optimization evaluator", "Results", new IntValue()));
     165
     166
     167
     168      Parameters.Add(new LookupParameter<DoubleValue>(AchievedQualityImprovementParameterName));
     169      Parameters.Add(new LookupParameter<DoubleValue>(NumberOfConstraintViolationsBeforeOptimizationParameterName));
     170      Parameters.Add(new LookupParameter<DoubleValue>(NumberOfConstraintViolationsAfterOptimizationParameterName));
     171      Parameters.Add(new LookupParameter<DoubleArray>(ConstraintsBeforeOptimizationParameterName));
     172      Parameters.Add(new LookupParameter<DoubleArray>(ViolationsAfterOptimizationParameterName));
     173      Parameters.Add(new LookupParameter<DoubleValue>(OptimizationDurationParameterName));
    136174    }
    137175
     
    141179
    142180    [StorableHook(HookType.AfterDeserialization)]
    143     private void AfterDeserialization() { }
     181    private void AfterDeserialization() {
     182      if (!Parameters.ContainsKey(AchievedQualityImprovementParameterName)) {
     183        Parameters.Add(new LookupParameter<DoubleValue>(AchievedQualityImprovementParameterName));
     184        Parameters.Add(new LookupParameter<DoubleValue>(NumberOfConstraintViolationsBeforeOptimizationParameterName));
     185        Parameters.Add(new LookupParameter<DoubleValue>(NumberOfConstraintViolationsAfterOptimizationParameterName));
     186      }
     187      if(!Parameters.ContainsKey(ConstraintsBeforeOptimizationParameterName)) {
     188        Parameters.Add(new LookupParameter<DoubleArray>(ConstraintsBeforeOptimizationParameterName));
     189        Parameters.Add(new LookupParameter<DoubleArray>(ViolationsAfterOptimizationParameterName));
     190        Parameters.Add(new LookupParameter<DoubleValue>(OptimizationDurationParameterName));
     191
     192      }
     193    }
    144194
    145195    private static readonly object locker = new object();
     
    151201        IEnumerable<int> constantOptimizationRows = GenerateRowsToEvaluate(ConstantOptimizationRowsPercentage.Value);
    152202        var counter = new EvaluationsCounter();
    153         quality = OptimizeConstants(SymbolicDataAnalysisTreeInterpreterParameter.ActualValue, solution, ProblemDataParameter.ActualValue,
    154            constantOptimizationRows, ApplyLinearScalingParameter.ActualValue.Value, Solver, ConstantOptimizationIterations.Value, updateVariableWeights: UpdateVariableWeights, lowerEstimationLimit: EstimationLimitsParameter.ActualValue.Lower, upperEstimationLimit: EstimationLimitsParameter.ActualValue.Upper, updateConstantsInTree: UpdateConstantsInTree, counter: counter);
    155 
    156203        if (ConstantOptimizationRowsPercentage.Value != RelativeNumberOfEvaluatedSamplesParameter.ActualValue.Value) {
    157204          throw new NotSupportedException();
    158205        }
     206
     207        var sw = new Stopwatch();
     208        sw.Start();
     209        quality = OptimizeConstants(SymbolicDataAnalysisTreeInterpreterParameter.ActualValue, solution, ProblemDataParameter.ActualValue,
     210           constantOptimizationRows, ApplyLinearScalingParameter.ActualValue.Value, Solver,
     211           out double qDiff, out double[] constraintsBefore, out double[] constraintsAfter,
     212           ConstantOptimizationIterations.Value, updateVariableWeights: UpdateVariableWeights, lowerEstimationLimit: EstimationLimitsParameter.ActualValue.Lower, upperEstimationLimit: EstimationLimitsParameter.ActualValue.Upper, updateConstantsInTree: UpdateConstantsInTree, counter: counter);
     213
     214        AchievedQualityImprovementParameter.ActualValue = new DoubleValue(qDiff);
     215        NumberOfConstraintViolationsBeforeOptimizationParameter.ActualValue = new DoubleValue(constraintsBefore.Count(cv => cv > 0));
     216        NumberOfConstraintViolationsAfterOptimizationParameter.ActualValue = new DoubleValue(constraintsAfter.Count(cv => cv > 0));
     217        ViolationsBeforeOptimizationParameter.ActualValue = new DoubleArray(constraintsBefore);
     218        ViolationsAfterOptimizationParameter.ActualValue = new DoubleArray(constraintsAfter);
     219        OptimizationDurationParameter.ActualValue = new DoubleValue(sw.Elapsed.TotalSeconds);
    159220
    160221        if (CountEvaluations) {
     
    202263    public static double OptimizeConstants(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter,
    203264      ISymbolicExpressionTree tree, IRegressionProblemData problemData, IEnumerable<int> rows, bool applyLinearScaling,
    204       string solver,
     265      string solver, out double qDiff, out double[] constraintsBefore, out double[] constraintsAfter,
    205266      int maxIterations, bool updateVariableWeights = true,
    206267      double lowerEstimationLimit = double.MinValue, double upperEstimationLimit = double.MaxValue,
    207       bool updateConstantsInTree = true, Action<double[], double, object> iterationCallback = null, EvaluationsCounter counter = null) {
     268      bool updateConstantsInTree = true, Action<double[], double, object> iterationCallback = null, EvaluationsCounter counter = null
     269      ) {
    208270
    209271      if (!updateVariableWeights) throw new NotSupportedException("not updating variable weights is not supported");
     
    213275
    214276      using (var state = new ConstrainedNLSInternal(solver, tree, maxIterations, problemData, 0, 0, 0)) {
     277        constraintsBefore = state.BestConstraintValues;
     278        double qBefore = state.BestError;
    215279        state.Optimize(ConstrainedNLSInternal.OptimizationMode.UpdateParameters);
    216         return state.BestError;
     280        constraintsAfter = state.BestConstraintValues;
     281        var qOpt = state.BestError;
     282        if (constraintsAfter.Any(cv => cv > 1e-8)) qOpt = qBefore;
     283        qDiff = qOpt - qBefore;
     284        return qOpt;
    217285      }
    218286    }
Note: See TracChangeset for help on using the changeset viewer.