From 543b37aab652b946bdc32179c4f921115c299836 Mon Sep 17 00:00:00 2001
From: Scott Stephens <stephens.js@gmail.com>
Date: Wed, 13 Feb 2013 22:44:08 -0600
Subject: [PATCH] Add most of GoldenSectionMinimizer implementation

---
 src/Numerics/Numerics.csproj                  |  3 +-
 .../Optimization/GoldenSectionMinimizer.cs    | 68 ++++++++++++-
 .../Optimization/ObjectiveFunction1D.cs       | 98 +++++++++++++++++++
 3 files changed, 166 insertions(+), 3 deletions(-)
 create mode 100644 src/Numerics/Optimization/ObjectiveFunction1D.cs

diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 2ee74519..7e7ca7e5 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -131,7 +131,6 @@
     <Compile Include="Optimization\IUnconstrainedMinimizer.cs" />
     <Compile Include="Optimization\MinimizationResult.cs" />
     <Compile Include="Optimization\MinimizationWithLineSearchResult.cs" />
-    <Compile Include="Optimization\NewtonMinimizer.cs" />
     <Compile Include="Precision.Comparison.cs" />
     <Compile Include="Precision.Equality.cs" />
     <Compile Include="Distributions\Bernoulli.cs" />
@@ -265,7 +264,9 @@
     <Compile Include="Optimization\BfgsMinimizer.cs" />
     <Compile Include="Optimization\ConjugateGradientMinimizer.cs" />
     <Compile Include="Optimization\GoldenSectionMinimizer.cs" />
+    <Compile Include="Optimization\NewtonMinimizer.cs" />
     <Compile Include="Optimization\ObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunction1D.cs" />
     <Compile Include="Optimization\OptimizationResult.cs" />
     <Compile Include="Optimization\LineSearch\StrongWolfeLineSearch.cs" />
     <Compile Include="SpecialFunctions\Evaluate.cs" />
diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs
index 516cceeb..6ce12c2d 100644
--- a/src/Numerics/Optimization/GoldenSectionMinimizer.cs
+++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs
@@ -1,6 +1,70 @@
-namespace MathNet.Numerics.Optimization
+using System;
+
+namespace MathNet.Numerics.Optimization
 {
-    class GoldenSectionMinimizer
+    public class GoldenSectionMinimizer
     {
+        public double XTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        public GoldenSectionMinimizer(double xTolerance=1e-5, int maxIterations=1000)
+        {
+            XTolerance = xTolerance;
+            MaximumIterations = maxIterations;
+        }
+
+        public MinimizationResult FindMinimum(IObjectiveFunction1D objective, double lowerBound, double upperBound)
+        {
+            double middlePointX = lowerBound + (upperBound - lowerBound) / (1 + Constants.GoldenRatio);
+            IEvaluation1D lower = objective.Evaluate(lowerBound);
+            IEvaluation1D middle = objective.Evaluate(middlePointX);
+            IEvaluation1D upper = objective.Evaluate(upperBound);
+
+            ValueChecker(lower.Value, lowerBound);
+            ValueChecker(middle.Value, middlePointX);
+            ValueChecker(upper.Value, upperBound);
+
+            if (upperBound <= lowerBound)
+                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) > XTolerance && iterations < MaximumIterations)
+            {
+                double testX = lower.Point + (upper.Point - middle.Point);
+                var test = objective.Evaluate(testX);
+                ValueChecker(test.Value, testX);
+
+                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 == MaximumIterations)
+                throw new MaximumIterationsException("Max iterations reached.");
+
+            return null;
+        }
+
+        private void ValueChecker(double value, double point)
+        {
+            if (Double.IsNaN(value) || Double.IsInfinity(value))
+                throw new Exception("Objective function returned non-finite value.");
+        }
     }
 }
diff --git a/src/Numerics/Optimization/ObjectiveFunction1D.cs b/src/Numerics/Optimization/ObjectiveFunction1D.cs
new file mode 100644
index 00000000..00a35d65
--- /dev/null
+++ b/src/Numerics/Optimization/ObjectiveFunction1D.cs
@@ -0,0 +1,98 @@
+using System;
+
+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? _secondDerivative;
+        private readonly SimpleObjectiveFunction1D _objectiveObject;
+        private readonly double _point;
+
+        public CachedEvaluation1D(SimpleObjectiveFunction1D f, double point)
+        {
+            _objectiveObject = f;
+            _point = point;
+        }
+        private double SetValue()
+        {
+            _value = _objectiveObject.Objective(_point);
+            return _value.Value;
+        }
+        private double SetDerivative()
+        {
+            _derivative = _objectiveObject.Derivative(_point);
+            return _derivative.Value;
+        }
+        private double SetSecondDerivative()
+        {
+            _secondDerivative = _objectiveObject.SecondDerivative(_point);
+            return _secondDerivative.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 _secondDerivative ?? SetSecondDerivative(); } }
+
+    }
+
+    public class SimpleObjectiveFunction1D : IObjectiveFunction1D
+    {
+        public Func<double, double> Objective { get; private set; }
+        public Func<double, double> Derivative { get; private set; }
+        public Func<double, double> SecondDerivative { get; private set; }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective)
+        {
+            Objective = objective;
+            Derivative = null;
+            SecondDerivative = null;
+        }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective, Func<double, double> derivative)
+        {
+            Objective = objective;
+            Derivative = derivative;
+            SecondDerivative = null;
+        }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective, Func<double, double> derivative, Func<double,double> secondDerivative)
+        {
+            Objective = objective;
+            Derivative = derivative;
+            SecondDerivative = secondDerivative;
+        }
+
+        public bool DerivativeSupported
+        {
+            get { return Derivative != null; }
+        }
+
+        public bool SecondDerivativeSupported
+        {
+            get { return SecondDerivative != null; }
+        }
+
+        public IEvaluation1D Evaluate(double point)
+        {
+            return new CachedEvaluation1D(this, point);
+        }
+    }
+}