add findLocalMin

std/numeric.d

@@ -1382,6 +1382,234 @@ unittest
     static assert(__traits(compiles, findRoot((real x)=>cast(double)x, real.init, real.init)));
 }
 
+/++
+Find a real minimum of a real function $(D f(x)) via bracketing.
+Given a function $(D f) and a range $(D (ax..bx)),
+returns the value of $(D x) in the range which is closest to a minimum of $(D f(x)).
+$(D f) is never evaluted at the endpoints of $(D ax) and $(D bx).
+If $(D f(x)) has more than one minimum in the range, one will be chosen arbitrarily.
+If $(D f(x)) returns NaN or -Infinity, $(D (x, f(x), NaN)) will be returned;
+otherwise, this algorithm is guaranteed to succeed.
+
+Params:
+    f = Function to be analyzed
+    ax = Left bound of initial range of f known to contain the minimum.
+    bx = Right bound of initial range of f known to contain the minimum.
+    relTolerance = Relative tolerance.
+    absTolerance = Absolute tolerance.
+
+Preconditions:
+    $(D ax) and $(D bx) shall be finite reals. $(BR)
+    $(D relTolerance) shall be normal positive real. $(BR)
+    $(D absTolerance) shall be normal positive real no less then $(D T.epsilon*2).
+
+Returns:
+    A tuple consisting of $(D x), $(D y = f(x)) and $(D error = 3 * (absTolerance * fabs(x) + relTolerance)).
+
+    The method used is a combination of golden section search and
+successive parabolic interpolation. Convergence is never much slower
+than that for a Fibonacci search.
+
+References:
+    "Algorithms for Minimization without Derivatives", Richard Brent, Prentice-Hall, Inc. (1973)
+
+See_Also: $(LREF findRoot), $(XREF math, isNormal)
++/
+Tuple!(T, "x", Unqual!(ReturnType!DF), "y", T, "error")
+findLocalMin(T, DF)(
+        scope DF f,
+        in T ax,
+        in T bx,
+        in T relTolerance = sqrt(T.epsilon),
+        in T absTolerance = sqrt(T.epsilon),
+        )
+    if (isFloatingPoint!T
+        && __traits(compiles, {T _ = DF.init(T.init);}))
+in
+{
+    assert(isFinite(ax), "ax is not finite");
+    assert(isFinite(bx), "bx is not finite");
+    assert(isNormal(relTolerance), "relTolerance is not normal floating point number");
+    assert(isNormal(absTolerance), "absTolerance is not normal floating point number");
+    assert(relTolerance >= 0, "absTolerance is not positive");
+    assert(absTolerance >= T.epsilon*2, "absTolerance is not greater then `2*T.epsilon`");
+}
+out (result)
+{
+    assert(isFinite(result.x));
+}
+body
+{
+    alias R = Unqual!(CommonType!(ReturnType!DF, T));
+    // c is the squared inverse of the golden ratio
+    // (3 - sqrt(5))/2
+    // Value obtained from Wolfram Alpha.
+    enum T c = 0x0.61c8864680b583ea0c633f9fa31237p+0L;
+    enum T cm1 = 0x0.9e3779b97f4a7c15f39cc0605cedc8p+0L;
+    R tolerance;
+    T a = ax > bx ? bx : ax;
+    T b = ax > bx ? ax : bx;
+    // sequence of declarations suitable for SIMD instructions
+    T  v = a * cm1 + b * c;
+    assert(isFinite(v));
+    R fv = f(v);
+    if(isNaN(fv) || fv == -T.infinity)
+    {
+        return typeof(return)(v, fv, T.init);
+    }
+    T  w = v;
+    R fw = fv;
+    T  x = v;
+    R fx = fv;
+    size_t i;
+    for(R d = 0, e = 0;;)
+    {
+        i++;
+        T m = (a + b) / 2;
+        // This fix is not part of the original algorithm
+        if(!isFinite(m)) // fix infinity loop. Issue can be reproduced in R.
+        {
+            m = a / 2 + b / 2;
+            if(!isFinite(m)) // fast-math compiler switch is enabled
+            {
+                //SIMD instructions can be used by compiler, do not reduce declarations
+                int a_exp = void;
+                int b_exp = void;
+                immutable an = frexp(a, a_exp);
+                immutable bn = frexp(b, b_exp);
+                immutable am = ldexp(an, a_exp-1);
+                immutable bm = ldexp(bn, b_exp-1);
+                m = am + bm;
+                if(!isFinite(m)) // wrong input: constraints are disabled in release mode
+                {
+                    return typeof(return).init;
+                }
+            }
+        }
+        tolerance = absTolerance * fabs(x) + relTolerance;
+        immutable t2 = tolerance * 2;
+        // check stopping criterion
+        if (!(fabs(x - m) > t2 - (b - a) / 2))
+        {
+            break;
+        }
+        R p = 0;
+        R q = 0;
+        R r = 0;
+        // fit parabola
+        if (fabs(e) > tolerance)
+        {
+            immutable  xw =  x -  w;
+            immutable fxw = fx - fw;
+            immutable  xv =  x -  v;
+            immutable fxv = fx - fv;
+            immutable xwfxv = xw * fxv;
+            immutable xvfxw = xv * fxw;
+            p = xv * xvfxw - xw * xwfxv;
+            q = (xvfxw - xwfxv) * 2;
+            if (q > 0)
+                p = -p;
+            else
+                q = -q;
+            r = e;
+            e = d;
+        }
+        T u;
+        // a parabolic-interpolation step
+        if (fabs(p) < fabs(q * r / 2) && p > q * (a - x) && p < q * (b - x))
+        {
+            d = p / q;
+            u = x + d;
+            // f must not be evaluated too close to a or b
+            if (u - a < t2 || b - u < t2)
+                d = x < m ? tolerance : -tolerance;
+        }
+        // a golden-section step
+        else
+        {
+            e = (x < m ? b : a) - x;
+            d = c * e;
+        }
+        // f must not be evaluated too close to x
+        u = x + (fabs(d) >= tolerance ? d : d > 0 ? tolerance : -tolerance);
+        immutable fu = f(u);
+        if(isNaN(fu) || fu == -T.infinity)
+        {
+            return typeof(return)(u, fu, T.init);
+        }
+        //  update  a, b, v, w, and x
+        if (fu <= fx)
+        {
+            u < x ? b : a = x;
+            v = w; fv = fw;
+            w = x; fw = fx;
+            x = u; fx = fu;
+        }
+        else
+        {
+            u < x ? a : b = u;
+            if (fu <= fw || w == x)
+            {
+                v = w; fv = fw;
+                w = u; fw = fu;
+            }
+            else if (fu <= fv || v == x || v == w)
+            { // do not remove this braces
+                v = u; fv = fu;
+            }
+        }
+    }
+    return typeof(return)(x, fx, tolerance * 3);
+}
+
+///
+unittest {
+    auto ret = findLocalMin((double x) => (x-4)^^2, -1e7, 1e7);
+    assert(ret.x.approxEqual(4.0));
+    assert(ret.y.approxEqual(0.0));
+}
+
+unittest {
+    import std.meta: AliasSeq;
+    foreach(T; AliasSeq!(double, float, real))
+    {
+        {
+            auto ret = findLocalMin!T((T x) => (x-4)^^2, T.min_normal, 1e7);
+            assert(ret.x.approxEqual(T(4)));
+            assert(ret.y.approxEqual(T(0)));
+        }
+        {
+            auto ret = findLocalMin!T((T x) => fabs(x-1), -T.max/4, T.max/4, T.min_normal, 2*T.epsilon);
+            assert(approxEqual(ret.x, T(1)));
+            assert(approxEqual(ret.y, T(0)));
+            assert(ret.error <= 10 * T.epsilon);
+        }
+        {
+            auto ret = findLocalMin!T((T x) => T.init, 0, 1, T.min_normal, 2*T.epsilon);
+            assert(!ret.x.isNaN);
+            assert(ret.y.isNaN);
+            assert(ret.error.isNaN);
+        }
+        {
+            auto ret = findLocalMin!T((T x) => log(x), 0, 1, T.min_normal, 2*T.epsilon);
+            assert(ret.error < 3.00001 * ((2*T.epsilon)*fabs(ret.x)+ T.min_normal));
+            assert(ret.x >= 0 && ret.x <= ret.error);
+        }
+        {
+            auto ret = findLocalMin!T((T x) => log(x), 0, T.max, T.min_normal, 2*T.epsilon);
+            assert(ret.y < -18);
+            assert(ret.error < 5e-08);
+            assert(ret.x >= 0 && ret.x <= ret.error);
+        }
+        {
+            auto ret = findLocalMin!T((T x) => -fabs(x), -1, 1, T.min_normal, 2*T.epsilon);
+            assert(ret.x.fabs.approxEqual(T(1)));
+            assert(ret.y.fabs.approxEqual(T(1)));
+            assert(ret.error.approxEqual(T(0)));
+        }
+    }
+}
+
 /**
 Computes $(LUCKY Euclidean distance) between input ranges $(D a) and
 $(D b). The two ranges must have the same length. The three-parameter