diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index aab028a8..71693a57 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -205,7 +205,9 @@ + + diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs index afa35d26..a43fda08 100644 --- a/src/Numerics/Optimization/GoldenSectionMinimizer.cs +++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs @@ -5,7 +5,69 @@ using System.Text; namespace MathNet.Numerics.Optimization { - class GoldenSectionMinimizer + public class GoldenSectionMinimizer { + public double XTolerance { get; set; } + public int MaximumIterations { get; set; } + + public GoldenSectionMinimizer(double x_tolerance=1e-5, int max_iterations=1000) + { + this.XTolerance = x_tolerance; + this.MaximumIterations = max_iterations; + } + + public MinimizationOutput FindMinimum(IObjectiveFunction1D objective, double lower_bound, double upper_bound) + { + if (!(objective is ObjectiveChecker1D)) + objective = new ObjectiveChecker1D(objective, this.ValueChecker, null, null); + + double middle_point_x = lower_bound + (upper_bound - lower_bound) / (1 + _golden_ratio); + IEvaluation1D lower = objective.Evaluate(lower_bound); + IEvaluation1D middle = objective.Evaluate(middle_point_x); + IEvaluation1D upper = objective.Evaluate(upper_bound); + + if (upper_bound <= lower_bound) + throw new OptimizationException("Lower bound must be lower than upper bound."); + + if (upper.Value < middle.Value || lower.Value < middle.Value) + throw new OptimizationException("Lower and upper bounds do not necessarily bound a minimum."); + + int iterations = 0; + while (Math.Abs(upper.Point - lower.Point) > this.XTolerance && iterations < this.MaximumIterations) + { + double test_x = lower.Point + (upper.Point - middle.Point); + var test = objective.Evaluate(test_x); + + if (test.Value > middle.Value) + { + if (test.Point < middle.Point) + lower = test; + else + upper = test; + } + else + { + if (test.Point < middle.Point) + upper = middle; + else + lower = middle; + } + + iterations += 1; + } + + if (iterations == this.MaximumIterations) + throw new MaximumIterationsException("Max iterations reached."); + else + return null; + + } + + private void ValueChecker(double value, double point) + { + if (Double.IsNaN(value) || Double.IsInfinity(value)) + throw new EvaluationException("Objective function returned non-finite value."); + } + private static double _golden_ratio = (1.0 + Math.Sqrt(5)) / 2.0; } } diff --git a/src/Numerics/Optimization/ObjectiveChecker1D.cs b/src/Numerics/Optimization/ObjectiveChecker1D.cs new file mode 100644 index 00000000..6ba0a7db --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveChecker1D.cs @@ -0,0 +1,131 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class CheckedEvaluation1D : IEvaluation1D + { + private ObjectiveChecker1D Checker; + private IEvaluation1D InnerEvaluation; + private bool ValueChecked; + private bool DerivativeChecked; + private bool SecondDerivativeChecked; + + public CheckedEvaluation1D(ObjectiveChecker1D checker, IEvaluation1D evaluation) + { + this.Checker = checker; + this.InnerEvaluation = evaluation; + } + + public double Point + { + get { return this.InnerEvaluation.Point; } + } + + public double Value + { + get + { + + if (!this.ValueChecked) + { + double tmp; + try + { + tmp = this.InnerEvaluation.Value; + } + catch (Exception e) + { + throw new EvaluationException("Objective function evaluation failed.", e); + } + this.Checker.ValueChecker(tmp, this.InnerEvaluation.Point); + } + return this.InnerEvaluation.Value; + } + } + + public double Derivative + { + get + { + + if (!this.DerivativeChecked) + { + double tmp; + try + { + tmp = this.InnerEvaluation.Derivative; + } + catch (Exception e) + { + throw new EvaluationException("Objective derivative evaluation failed.", e); + } + this.Checker.DerivativeChecker(tmp, this.InnerEvaluation.Point); + } + return this.InnerEvaluation.Derivative; + } + } + + public double SecondDerivative + { + get + { + + if (!this.SecondDerivativeChecked) + { + double tmp; + try + { + tmp = this.InnerEvaluation.SecondDerivative; + } + catch (Exception e) + { + throw new EvaluationException("Objective second derivative evaluation failed.", e); + } + this.Checker.SecondDerivativeChecker(tmp, this.InnerEvaluation.Point); + } + return this.InnerEvaluation.SecondDerivative; + } + } + } + + public class ObjectiveChecker1D : IObjectiveFunction1D + { + public IObjectiveFunction1D InnerObjective { get; private set; } + public Action ValueChecker { get; private set; } + public Action DerivativeChecker { get; private set; } + public Action SecondDerivativeChecker { get; private set; } + + public ObjectiveChecker1D(IObjectiveFunction1D objective, Action value_checker, Action gradient_checker, Action hessian_checker) + { + this.InnerObjective = objective; + this.ValueChecker = value_checker; + this.DerivativeChecker = gradient_checker; + this.SecondDerivativeChecker = hessian_checker; + } + + public bool DerivativeSupported + { + get { return this.InnerObjective.DerivativeSupported; } + } + + public bool SecondDerivativeSupported + { + get { return this.InnerObjective.SecondDerivativeSupported; } + } + + public IEvaluation1D Evaluate(double point) + { + try + { + return new CheckedEvaluation1D(this, this.InnerObjective.Evaluate(point)); + } + catch (Exception e) + { + throw new EvaluationException("Objective evaluation failed.", e); + } + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunction1D.cs b/src/Numerics/Optimization/ObjectiveFunction1D.cs new file mode 100644 index 00000000..8611eda4 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunction1D.cs @@ -0,0 +1,101 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public interface IEvaluation1D + { + double Point { get; } + double Value { get; } + double Derivative { get; } + double SecondDerivative { get; } + } + + public interface IObjectiveFunction1D + { + bool DerivativeSupported { get; } + bool SecondDerivativeSupported { get; } + IEvaluation1D Evaluate(double point); + } + + public class CachedEvaluation1D : IEvaluation1D + { + private double? _value; + private double? _derivative; + private double? _second_derivative; + private SimpleObjectiveFunction1D _objective_object; + private double _point; + + public CachedEvaluation1D(SimpleObjectiveFunction1D f, double point) + { + _objective_object = f; + _point = point; + } + private double setValue() + { + _value = _objective_object.Objective(_point); + return _value.Value; + } + private double setDerivative() + { + _derivative = _objective_object.Derivative(_point); + return _derivative.Value; + } + private double setSecondDerivative() + { + _second_derivative = _objective_object.SecondDerivative(_point); + return _second_derivative.Value; + } + + public double Point { get { return _point; } } + public double Value { get { return _value ?? setValue(); } } + public double Derivative { get { return _derivative ?? setDerivative(); } } + public double SecondDerivative { get { return _second_derivative ?? setSecondDerivative(); } } + + } + + public class SimpleObjectiveFunction1D : IObjectiveFunction1D + { + public Func Objective { get; private set; } + public Func Derivative { get; private set; } + public Func SecondDerivative { get; private set; } + + public SimpleObjectiveFunction1D(Func objective) + { + this.Objective = objective; + this.Derivative = null; + this.SecondDerivative = null; + } + + public SimpleObjectiveFunction1D(Func objective, Func derivative) + { + this.Objective = objective; + this.Derivative = derivative; + this.SecondDerivative = null; + } + + public SimpleObjectiveFunction1D(Func objective, Func derivative, Func second_derivative) + { + this.Objective = objective; + this.Derivative = derivative; + this.SecondDerivative = second_derivative; + } + + public bool DerivativeSupported + { + get { return this.Derivative != null; } + } + + public bool SecondDerivativeSupported + { + get { return this.SecondDerivative != null; } + } + + public IEvaluation1D Evaluate(double point) + { + return new CachedEvaluation1D(this, point); + } + } +}