"""ch13_golden_section.py
Companion code for Chapter 13: Optimization and Data Fitting.

One-dimensional golden-section search.
Edit the example settings in main(); numerical routines are reusable on import.
Run: python3 ch13_golden_section.py [--no-plot]
Requires NumPy; optional plots require Matplotlib.
Fixed seeds reproduce runs with the same software; languages need not share sequences.
"""

import sys
import numpy as np


def main():
    make_plot = "--no-plot" not in sys.argv
    fun = lambda x: (x-2.0)**2+1.0
    a,b = -1.0,4.0
    tol_x,max_iter = 1e-6,200
    x,fx,info = golden_section(fun,a,b,tol_x,max_iter)
    report(x,fx,info)
    if make_plot: plot_scalar(fun,a,b,x,fx,"Golden-section search")


def vector(x, name="data"):
    if np.iscomplexobj(x):
        raise ValueError(name + " must be real")
    a = np.asarray(x, dtype=float)
    if a.ndim != 1 or a.size == 0 or not np.all(np.isfinite(a)):
        raise ValueError(name + " must be a nonempty finite vector")
    return a


def checked(fun, x):
    value = np.asarray(fun(x))
    if value.ndim != 0 or np.iscomplexobj(value) or not np.isfinite(value):
        raise ValueError("Objective must return a finite real scalar")
    return float(value)


def bounds(lower, upper):
    lo, hi = vector(lower, "lower"), vector(upper, "upper")
    if lo.shape != hi.shape or not np.all(np.isfinite(hi - lo)) or np.any(hi <= lo):
        raise ValueError("Require matching finite bounds with upper > lower")
    return lo, hi


def integer(value, minimum, name):
    if not np.isscalar(value) or not np.isfinite(value) or int(value) != value or value < minimum:
        raise ValueError(name + " must be an integer >= " + str(minimum))
    return int(value)


def positive(value, name):
    if not np.isscalar(value) or not np.isfinite(value) or value <= 0:
        raise ValueError(name + " must be positive and finite")
    return float(value)

def golden_section(fun, a, b, tol_x=1e-6, max_iter=200):
    """Minimize a unimodal function; tolerance controls final interval width.

    Rounded objective values can spoil localization near a minimum. The final
    interval is not a rigorous error certificate for rounded evaluations.
    """
    lo, hi = bounds([a], [b]); a, b = float(lo[0]), float(hi[0])
    positive(tol_x, "tol_x"); max_iter = integer(max_iter, 0, "max_iter")
    rho = (np.sqrt(5.0) - 1.0) / 2.0
    x1, x2 = b - rho*(b-a), a + rho*(b-a)
    f1, f2 = checked(fun, x1), checked(fun, x2)
    calls, iterations = 2, 0
    status = "Maximum iteration count reached."
    while b-a > tol_x and iterations < max_iter:
        if not a < x1 < x2 < b:
            status = "Floating-point interval stagnation."
            break
        if f1 > f2:
            a, x1, f1 = x1, x2, f2
            x2 = a + rho*(b-a); f2 = checked(fun, x2)
        else:
            b, x2, f2 = x2, x1, f1
            x1 = b - rho*(b-a); f1 = checked(fun, x1)
        calls += 1; iterations += 1
    if b-a <= tol_x:
        status = "Interval-width tolerance reached."
    x = a + 0.5*(b-a); fx = checked(fun, x)
    return x, fx, dict(a=a, b=b, iterations=iterations, evaluations=calls+1,
                       converged=b-a <= tol_x, status=status)

def report(p, value, info):
    for j,x in enumerate(np.atleast_1d(p),1):
        print(f"p({j}) = {x:.12g}")
    print(f"Objective value: {value:.12g}")
    for key,value in info.items():
        print(f"{key}: {value}")

def plot_scalar(fun, a, b, x, fx, title):
    import matplotlib.pyplot as plt
    xx=np.linspace(a,b,400)
    plt.plot(xx,[fun(z) for z in xx],"k-",label="Objective")
    plt.plot(x,fx,"ko",label="Best point");plt.xlabel("x");plt.ylabel("f(x)")
    plt.title(title,fontweight="normal");plt.xlim(a,b);plt.legend(loc="best");plt.grid(True);plt.tight_layout();plt.show()


if __name__ == "__main__":
    main()