Changeset 17215
 Timestamp:
 08/13/19 19:02:57 (10 days ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/NLOptEvaluator.cs
r17196 r17215 216 216 } 217 217 218 // for data exchange to/from optimizer in native code219 [StructLayout(LayoutKind.Sequential)]220 private struct ConstraintData {221 public ISymbolicExpressionTree Tree;222 public ISymbolicExpressionTreeNode[] ParameterNodes;223 }224 218 225 219 public static double OptimizeConstants(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, … … 234 228 if (!applyLinearScaling) throw new NotSupportedException("application without linear scaling is not supported"); 235 229 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; 307 234 } 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/callbackdelegatesbeingcollected#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], 1e8); 397 } 398 399 NLOpt.nlopt_set_xtol_rel(nlopt_opt, 1e4); 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 603 238 } 604 239 }
Note: See TracChangeset
for help on using the changeset viewer.