Changeset 17215


Ignore:
Timestamp:
08/13/19 19:02:57 (10 days ago)
Author:
gkronber
Message:

#2994: use new code for NLOpt in the NLOpt evaluator

File:
1 edited

Legend:

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

    r17196 r17215  
    216216    }
    217217
    218     // for data exchange to/from optimizer in native code
    219     [StructLayout(LayoutKind.Sequential)]
    220     private struct ConstraintData {
    221       public ISymbolicExpressionTree Tree;
    222       public ISymbolicExpressionTreeNode[] ParameterNodes;
    223     }
    224218
    225219    public static double OptimizeConstants(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter,
     
    234228      if (!applyLinearScaling) throw new NotSupportedException("application without linear scaling is not supported");
    235229
    236       // we always update constants, so we don't need to calculate initial quality
    237       // double originalQuality = SymbolicRegressionSingleObjectiveMeanSquaredErrorEvaluator.Calculate(interpreter, tree, lowerEstimationLimit, upperEstimationLimit, problemData, rows, applyLinearScaling: false);
    238 
    239       if (counter == null) counter = new EvaluationsCounter();
    240       var rowEvaluationsCounter = new EvaluationsCounter();
    241 
    242       var intervalConstraints = problemData.IntervalConstraints;
    243       var dataIntervals = problemData.VariableRanges.GetIntervals();
    244       var trainingRows = problemData.TrainingIndices.ToArray();
    245       // buffers
    246       var target = problemData.TargetVariableTrainingValues.ToArray();
    247       var targetStDev = target.StandardDeviationPop();
    248       var targetVariance = targetStDev * targetStDev;
    249       var targetMean = target.Average();
    250       var pred = interpreter.GetSymbolicExpressionTreeValues(tree, problemData.Dataset, trainingRows).ToArray();
    251       if (pred.Any(pi => double.IsInfinity(pi) || double.IsNaN(pi))) return targetVariance;
    252 
    253       #region linear scaling
    254       var predStDev = pred.StandardDeviationPop();
    255       if (predStDev == 0) return targetVariance; // constant expression
    256       var predMean = pred.Average();
    257 
    258       var scalingFactor = targetStDev / predStDev;
    259       var offset = targetMean - predMean * scalingFactor;
    260 
    261       ISymbolicExpressionTree scaledTree = null;
    262       if (applyLinearScaling) scaledTree = CopyAndScaleTree(tree, scalingFactor, offset);
    263       #endregion
    264 
    265       // convert constants to variables named theta...
    266       var treeForDerivation = ReplaceConstWithVar(scaledTree, out List<string> thetaNames, out List<double> thetaValues); // copies the tree
    267 
    268       // create trees for relevant derivatives
    269       Dictionary<string, ISymbolicExpressionTree> derivatives = new Dictionary<string, ISymbolicExpressionTree>();
    270       var allThetaNodes = thetaNames.Select(_ => new List<ConstantTreeNode>()).ToArray();
    271       var constraintTrees = new List<ISymbolicExpressionTree>();
    272       foreach (var constraint in intervalConstraints.Constraints) {
    273         if (constraint.IsDerivation) {
    274           if (!problemData.AllowedInputVariables.Contains(constraint.Variable))
    275             throw new ArgumentException($"Invalid constraint: the variable {constraint.Variable} does not exist in the dataset.");
    276           var df = DerivativeCalculator.Derive(treeForDerivation, constraint.Variable);
    277 
    278           // alglib requires constraint expressions of the form c(x) <= 0
    279           // -> we make two expressions, one for the lower bound and one for the upper bound
    280 
    281           if (constraint.Interval.UpperBound < double.PositiveInfinity) {
    282             var df_smaller_upper = Subtract((ISymbolicExpressionTree)df.Clone(), CreateConstant(constraint.Interval.UpperBound));
    283             // convert variables named theta back to constants
    284             var df_prepared = ReplaceVarWithConst(df_smaller_upper, thetaNames, thetaValues, allThetaNodes);
    285             constraintTrees.Add(df_prepared);
    286           }
    287           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
    288             var df_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)df.Clone());
    289             // convert variables named theta back to constants
    290             var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);
    291             constraintTrees.Add(df_prepared);
    292           }
    293         } else {
    294           if (constraint.Interval.UpperBound < double.PositiveInfinity) {
    295             var f_smaller_upper = Subtract((ISymbolicExpressionTree)treeForDerivation.Clone(), CreateConstant(constraint.Interval.UpperBound));
    296             // convert variables named theta back to constants
    297             var df_prepared = ReplaceVarWithConst(f_smaller_upper, thetaNames, thetaValues, allThetaNodes);
    298             constraintTrees.Add(df_prepared);
    299           }
    300           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
    301             var f_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)treeForDerivation.Clone());
    302             // convert variables named theta back to constants
    303             var df_prepared = ReplaceVarWithConst(f_larger_lower, thetaNames, thetaValues, allThetaNodes);
    304             constraintTrees.Add(df_prepared);
    305           }
    306         }
     230
     231      using (var state = new ConstrainedNLSInternal(solver, tree, maxIterations, problemData, 0, 0, 0)) {
     232        state.Optimize();
     233        return state.BestError;
    307234      }
    308 
    309       var preparedTree = ReplaceVarWithConst(treeForDerivation, thetaNames, thetaValues, allThetaNodes);
    310       var preparedTreeParameterNodes = GetParameterNodes(preparedTree, allThetaNodes);
    311 
    312       // local function
    313       void UpdateThetaValues(double[] theta) {
    314         for (int i = 0; i < theta.Length; ++i) {
    315           foreach (var constNode in allThetaNodes[i]) constNode.Value = theta[i];
    316         }
    317       }
    318 
    319       var fi_eval = new double[target.Length];
    320       var jac_eval = new double[target.Length, thetaValues.Count];
    321 
    322       double calculate_obj(uint dim, double[] curX, double[] grad, IntPtr data) {
    323         UpdateThetaValues(curX);
    324 
    325         if (grad != null) {
    326           var autoDiffEval = new VectorAutoDiffEvaluator();
    327           autoDiffEval.Evaluate(preparedTree, problemData.Dataset, trainingRows,
    328             preparedTreeParameterNodes, fi_eval, jac_eval);
    329 
    330           // calc sum of squared errors and gradient
    331           var sse = 0.0;
    332           for (int j = 0; j < grad.Length; j++) grad[j] = 0;
    333           for (int i = 0; i < target.Length; i++) {
    334             var r = target[i] - fi_eval[i];
    335             sse += 0.5 * r * r;
    336             for (int j = 0; j < grad.Length; j++) {
    337               grad[j] -= r * jac_eval[i, j];
    338             }
    339           }
    340           if (double.IsNaN(sse)) return double.MaxValue;
    341           // average
    342           for (int j = 0; j < grad.Length; j++) { grad[j] /= target.Length; }
    343           return sse / target.Length;
    344         } else {
    345           var eval = new VectorEvaluator();
    346           var prediction = eval.Evaluate(preparedTree, problemData.Dataset, trainingRows);
    347 
    348           // calc sum of squared errors and gradient
    349           var sse = 0.0;
    350           for (int i = 0; i < target.Length; i++) {
    351             var r = target[i] - prediction[i];
    352             sse += 0.5 * r * r;
    353           }
    354           if (double.IsNaN(sse)) return double.MaxValue;
    355           // average
    356           return sse / target.Length;
    357         }
    358 
    359       }
    360 
    361       double calculate_constraint(uint dim, double[] curX, double[] grad, IntPtr data) {
    362         UpdateThetaValues(curX);
    363         var intervalEvaluator = new IntervalEvaluator();
    364         var constraintData = Marshal.PtrToStructure<ConstraintData>(data);
    365 
    366         if (grad != null) for (int j = 0; j < grad.Length; j++) grad[j] = 0; // clear grad
    367 
    368         var interval = intervalEvaluator.Evaluate(constraintData.Tree, dataIntervals, constraintData.ParameterNodes,
    369           out double[] lowerGradient, out double[] upperGradient);
    370 
    371         // we transformed this to a constraint c(x) <= 0, so only the upper bound is relevant for us
    372         if (grad != null) for (int j = 0; j < grad.Length; j++) { grad[j] = upperGradient[j]; }
    373         if (double.IsNaN(interval.UpperBound)) return double.MaxValue;
    374         else return interval.UpperBound;
    375       }
    376 
    377       var minVal = Math.Min(-1000.0, thetaValues.Min());
    378       var maxVal = Math.Max(1000.0, thetaValues.Max());
    379       var lb = Enumerable.Repeat(minVal, thetaValues.Count).ToArray();
    380       var up = Enumerable.Repeat(maxVal, thetaValues.Count).ToArray();
    381       IntPtr nlopt_opt = NLOpt.nlopt_create(GetAlgFromIdentifier(solver), (uint)thetaValues.Count); /* algorithm and dimensionality */
    382 
    383       NLOpt.nlopt_set_lower_bounds(nlopt_opt, lb);
    384       NLOpt.nlopt_set_upper_bounds(nlopt_opt, up);
    385       var calculateObjectiveDelegate = new NLOpt.nlopt_func(calculate_obj); // keep a reference to the delegate (see below)
    386       NLOpt.nlopt_set_min_objective(nlopt_opt, calculateObjectiveDelegate, IntPtr.Zero);
    387 
    388 
    389       var constraintDataPtr = new IntPtr[constraintTrees.Count];
    390       var 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)
    391       for (int i = 0; i < constraintTrees.Count; i++) {
    392         var constraintData = new ConstraintData() { Tree = constraintTrees[i], ParameterNodes = GetParameterNodes(constraintTrees[i], allThetaNodes) };
    393         constraintDataPtr[i] = Marshal.AllocHGlobal(Marshal.SizeOf<ConstraintData>());
    394         Marshal.StructureToPtr(constraintData, constraintDataPtr[i], fDeleteOld: false);
    395         calculateConstraintDelegates[i] = new NLOpt.nlopt_func(calculate_constraint);
    396         NLOpt.nlopt_add_inequality_constraint(nlopt_opt, calculateConstraintDelegates[i], constraintDataPtr[i], 1e-8);
    397       }
    398 
    399       NLOpt.nlopt_set_xtol_rel(nlopt_opt, 1e-4);
    400       NLOpt.nlopt_set_maxtime(nlopt_opt, 10.0); // 10 secs
    401       NLOpt.nlopt_set_maxeval(nlopt_opt, maxIterations);
    402 
    403       var x = thetaValues.ToArray();  /* initial guess */
    404       double minf = double.MaxValue; /* minimum objective value upon return */
    405       var res = NLOpt.nlopt_optimize(nlopt_opt, x, ref minf);
    406       if (res < 0) {
    407         throw new InvalidOperationException($"NLOpt failed {res} {NLOpt.nlopt_get_errmsg(nlopt_opt)}");
    408       } else {
    409         // update parameters in tree
    410         var pIdx = 0;
    411         // here we lose the two last parameters (for linear scaling)
    412         foreach (var node in tree.IterateNodesPostfix()) {
    413           if (node is ConstantTreeNode constTreeNode) {
    414             constTreeNode.Value = x[pIdx++];
    415           } else if (node is VariableTreeNode varTreeNode) {
    416             varTreeNode.Weight = x[pIdx++];
    417           }
    418         }
    419         // note: we keep the optimized constants even when the tree is worse.
    420         // assert that we lose the last two parameters
    421         if (pIdx != x.Length - 2) throw new InvalidProgramException();
    422       }
    423 
    424       NLOpt.nlopt_destroy(nlopt_opt);
    425       for (int i = 0; i < constraintDataPtr.Length; i++)
    426         Marshal.FreeHGlobal(constraintDataPtr[i]);
    427 
    428       counter.FunctionEvaluations += NLOpt.nlopt_get_numevals(nlopt_opt);
    429       // counter.GradientEvaluations += NLOpt.nlopt_get; // TODO
    430 
    431 
    432 
    433       return Math.Min(minf, targetVariance);
    434     }
    435 
    436     private static NLOpt.nlopt_algorithm GetAlgFromIdentifier(string solver) {
    437       if (solver.Contains("MMA")) return NLOpt.nlopt_algorithm.NLOPT_LD_MMA;
    438       if (solver.Contains("COBYLA")) return NLOpt.nlopt_algorithm.NLOPT_LN_COBYLA;
    439       if (solver.Contains("CCSAQ")) return NLOpt.nlopt_algorithm.NLOPT_LD_CCSAQ;
    440       if (solver.Contains("ISRES")) return NLOpt.nlopt_algorithm.NLOPT_GN_ISRES;
    441       throw new ArgumentException($"Unknown solver {solver}");
    442     }
    443 
    444     private static ISymbolicExpressionTree CopyAndScaleTree(ISymbolicExpressionTree tree, double scalingFactor, double offset) {
    445       var m = (ISymbolicExpressionTree)tree.Clone();
    446 
    447       var add = MakeNode<Addition>(MakeNode<Multiplication>(m.Root.GetSubtree(0).GetSubtree(0), CreateConstant(scalingFactor)), CreateConstant(offset));
    448       m.Root.GetSubtree(0).RemoveSubtree(0);
    449       m.Root.GetSubtree(0).AddSubtree(add);
    450       return m;
    451     }
    452 
    453     #region helper
    454     private static ISymbolicExpressionTreeNode[] GetParameterNodes(ISymbolicExpressionTree tree, List<ConstantTreeNode>[] allNodes) {
    455       // TODO better solution necessary
    456       var treeConstNodes = tree.IterateNodesPostfix().OfType<ConstantTreeNode>().ToArray();
    457       var paramNodes = new ISymbolicExpressionTreeNode[allNodes.Length];
    458       for (int i = 0; i < paramNodes.Length; i++) {
    459         paramNodes[i] = allNodes[i].SingleOrDefault(n => treeConstNodes.Contains(n));
    460       }
    461       return paramNodes;
    462     }
    463 
    464     private static ISymbolicExpressionTree ReplaceVarWithConst(ISymbolicExpressionTree tree, List<string> thetaNames, List<double> thetaValues, List<ConstantTreeNode>[] thetaNodes) {
    465       var copy = (ISymbolicExpressionTree)tree.Clone();
    466       var nodes = copy.IterateNodesPostfix().ToList();
    467       for (int i = 0; i < nodes.Count; i++) {
    468         var n = nodes[i] as VariableTreeNode;
    469         if (n != null) {
    470           var thetaIdx = thetaNames.IndexOf(n.VariableName);
    471           if (thetaIdx >= 0) {
    472             var parent = n.Parent;
    473             if (thetaNodes[thetaIdx].Any()) {
    474               // HACK: REUSE CONSTANT TREE NODE IN SEVERAL TREES
    475               // we use this trick to allow autodiff over thetas when thetas occurr multiple times in the tree (e.g. in derived trees)
    476               var constNode = thetaNodes[thetaIdx].First();
    477               var childIdx = parent.IndexOfSubtree(n);
    478               parent.RemoveSubtree(childIdx);
    479               parent.InsertSubtree(childIdx, constNode);
    480             } else {
    481               var constNode = (ConstantTreeNode)CreateConstant(thetaValues[thetaIdx]);
    482               var childIdx = parent.IndexOfSubtree(n);
    483               parent.RemoveSubtree(childIdx);
    484               parent.InsertSubtree(childIdx, constNode);
    485               thetaNodes[thetaIdx].Add(constNode);
    486             }
    487           }
    488         }
    489       }
    490       return copy;
    491     }
    492 
    493     private static ISymbolicExpressionTree ReplaceConstWithVar(ISymbolicExpressionTree tree, out List<string> thetaNames, out List<double> thetaValues) {
    494       thetaNames = new List<string>();
    495       thetaValues = new List<double>();
    496       var copy = (ISymbolicExpressionTree)tree.Clone();
    497       var nodes = copy.IterateNodesPostfix().ToList();
    498 
    499       int n = 1;
    500       for (int i = 0; i < nodes.Count; ++i) {
    501         var node = nodes[i];
    502         if (node is ConstantTreeNode constantTreeNode) {
    503           var thetaVar = (VariableTreeNode)new Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
    504           thetaVar.Weight = 1;
    505           thetaVar.VariableName = $"θ{n++}";
    506 
    507           thetaNames.Add(thetaVar.VariableName);
    508           thetaValues.Add(constantTreeNode.Value);
    509 
    510           var parent = constantTreeNode.Parent;
    511           if (parent != null) {
    512             var index = constantTreeNode.Parent.IndexOfSubtree(constantTreeNode);
    513             parent.RemoveSubtree(index);
    514             parent.InsertSubtree(index, thetaVar);
    515           }
    516         }
    517         if (node is VariableTreeNode varTreeNode) {
    518           var thetaVar = (VariableTreeNode)new Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
    519           thetaVar.Weight = 1;
    520           thetaVar.VariableName = $"θ{n++}";
    521 
    522           thetaNames.Add(thetaVar.VariableName);
    523           thetaValues.Add(varTreeNode.Weight);
    524 
    525           var parent = varTreeNode.Parent;
    526           if (parent != null) {
    527             var index = varTreeNode.Parent.IndexOfSubtree(varTreeNode);
    528             parent.RemoveSubtree(index);
    529             var prodNode = MakeNode<Multiplication>();
    530             varTreeNode.Weight = 1.0;
    531             prodNode.AddSubtree(varTreeNode);
    532             prodNode.AddSubtree(thetaVar);
    533             parent.InsertSubtree(index, prodNode);
    534           }
    535         }
    536       }
    537       return copy;
    538     }
    539 
    540     private static ISymbolicExpressionTreeNode CreateConstant(double value) {
    541       var constantNode = (ConstantTreeNode)new Constant().CreateTreeNode();
    542       constantNode.Value = value;
    543       return constantNode;
    544     }
    545 
    546     private static ISymbolicExpressionTree Subtract(ISymbolicExpressionTree t, ISymbolicExpressionTreeNode b) {
    547       var sub = MakeNode<Subtraction>(t.Root.GetSubtree(0).GetSubtree(0), b);
    548       t.Root.GetSubtree(0).RemoveSubtree(0);
    549       t.Root.GetSubtree(0).InsertSubtree(0, sub);
    550       return t;
    551     }
    552     private static ISymbolicExpressionTree Subtract(ISymbolicExpressionTreeNode b, ISymbolicExpressionTree t) {
    553       var sub = MakeNode<Subtraction>(b, t.Root.GetSubtree(0).GetSubtree(0));
    554       t.Root.GetSubtree(0).RemoveSubtree(0);
    555       t.Root.GetSubtree(0).InsertSubtree(0, sub);
    556       return t;
    557     }
    558 
    559     private static ISymbolicExpressionTreeNode MakeNode<T>(params ISymbolicExpressionTreeNode[] fs) where T : ISymbol, new() {
    560       var node = new T().CreateTreeNode();
    561       foreach (var f in fs) node.AddSubtree(f);
    562       return node;
    563     }
    564     #endregion
    565 
    566     private static void UpdateConstants(ISymbolicExpressionTreeNode[] nodes, double[] constants) {
    567       if (nodes.Length != constants.Length) throw new InvalidOperationException();
    568       for (int i = 0; i < nodes.Length; i++) {
    569         if (nodes[i] is VariableTreeNode varNode) varNode.Weight = constants[i];
    570         else if (nodes[i] is ConstantTreeNode constNode) constNode.Value = constants[i];
    571       }
    572     }
    573 
    574     private static alglib.ndimensional_fvec CreateFunc(ISymbolicExpressionTree tree, VectorEvaluator eval, ISymbolicExpressionTreeNode[] parameterNodes, IDataset ds, string targetVar, int[] rows) {
    575       var y = ds.GetDoubleValues(targetVar, rows).ToArray();
    576       return (double[] c, double[] fi, object o) => {
    577         UpdateConstants(parameterNodes, c);
    578         var pred = eval.Evaluate(tree, ds, rows);
    579         for (int i = 0; i < fi.Length; i++)
    580           fi[i] = pred[i] - y[i];
    581 
    582         var counter = (EvaluationsCounter)o;
    583         counter.FunctionEvaluations++;
    584       };
    585     }
    586 
    587     private static alglib.ndimensional_jac CreateJac(ISymbolicExpressionTree tree, VectorAutoDiffEvaluator eval, ISymbolicExpressionTreeNode[] parameterNodes, IDataset ds, string targetVar, int[] rows) {
    588       var y = ds.GetDoubleValues(targetVar, rows).ToArray();
    589       return (double[] c, double[] fi, double[,] jac, object o) => {
    590         UpdateConstants(parameterNodes, c);
    591         eval.Evaluate(tree, ds, rows, parameterNodes, fi, jac);
    592 
    593         for (int i = 0; i < fi.Length; i++)
    594           fi[i] -= y[i];
    595 
    596         var counter = (EvaluationsCounter)o;
    597         counter.GradientEvaluations++;
    598       };
    599     }
    600     public static bool CanOptimizeConstants(ISymbolicExpressionTree tree) {
    601       return TreeToAutoDiffTermConverter.IsCompatible(tree);
    602     }
     235    }
     236
     237
    603238  }
    604239}
Note: See TracChangeset for help on using the changeset viewer.