// ch13_gradient_descent.cpp
// Companion code for Chapter 13: Optimization and Data Fitting.
// Standalone C++11; standard library only; console output without plotting.
// Compile: c++ -std=c++11 -O2 -Wall -Wextra -pedantic ch13_gradient_descent.cpp -o ch13_gradient_descent
// Edit the example settings in main(); routines above main() can be reused.
// Fixed seeds reproduce runs with the same standard-library implementation.
// Random distributions need not give identical sequences across compilers/languages.

#include <algorithm>
#include <array>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

using Vec = std::vector<double>;
using Mat = std::vector<Vec>;
using Function = std::function<double(const Vec&)>;
using ScalarFunction = std::function<double(double)>;
const double pi = std::acos(-1.0);

void require(bool condition, const std::string& message) {
    if (!condition) throw std::invalid_argument(message);
}
void finite_vector(const Vec& x) {
    require(!x.empty(), "A nonempty vector is required.");
    for (double v : x) require(std::isfinite(v), "Data must be finite.");
}
void check_bounds(const Vec& lo, const Vec& hi) {
    finite_vector(lo); finite_vector(hi);
    require(lo.size() == hi.size(), "Bounds must have matching dimensions.");
    for (std::size_t j=0; j<lo.size(); ++j)
        require(hi[j]>lo[j] && std::isfinite(hi[j]-lo[j]), "Require finite upper > lower.");
}
double objective_value(double value) {
    require(std::isfinite(value), "Objective must return a finite real scalar.");
    return value;
}
double clip(double x, double lo, double hi) { return std::max(lo,std::min(hi,x)); }
Vec linspace(double a, double b, int n) {
    require(n>=2, "Grid needs at least two points.");
    Vec x(n);
    for (int i=0; i<n; ++i) x[i]=a+(b-a)*i/(n-1);
    return x;
}
struct Random {
    std::mt19937 engine;
    std::uniform_real_distribution<double> uniform{0.0,1.0};
    std::normal_distribution<double> normal{0.0,1.0};
    explicit Random(unsigned int seed) : engine(seed) {}
    double u() { return uniform(engine); }
    double z() { return normal(engine); }
    int index(int n) { return std::uniform_int_distribution<int>(0,n-1)(engine); }
};
struct Result {
    Vec point;
    double value=0.0;
    long long evaluations=0;
    int iterations=0;
    std::string status;
};
void print_result(const Result& r) {
    std::cout << std::setprecision(12);
    for (std::size_t j=0; j<r.point.size(); ++j)
        std::cout << "p(" << j+1 << ") = " << r.point[j] << '\n';
    std::cout << "Objective value: " << r.value << "\nIterations/updates: " << r.iterations
              << "\nFunction calls: " << r.evaluations << "\nStatus: " << r.status << '\n';
}

// Armijo backtracking; a small gradient establishes only stationarity.
Result gradient_descent(const Function& f, const std::function<Vec(const Vec&)>& grad,
                        Vec p, double step_initial, double tol_grad, int max_iter,
                        double& gradient_norm, long long& gradient_calls) {
    finite_vector(p);
    require(std::isfinite(step_initial) && step_initial>0 && std::isfinite(tol_grad)
            && tol_grad>0 && max_iter>=0, "Invalid gradient-descent settings.");
    Result r; r.value=objective_value(f(p)); r.evaluations=1;
    r.status="Maximum iteration count reached."; gradient_calls=0;
    while (true) {
        Vec g=grad(p); finite_vector(g); ++gradient_calls;
        require(g.size()==p.size(), "Gradient has the wrong dimension.");
        gradient_norm=0;
        for (double v:g) gradient_norm=std::hypot(gradient_norm,v);
        require(std::isfinite(gradient_norm), "Gradient norm overflow; rescale the problem.");
        if (gradient_norm<=tol_grad) { r.status="Gradient tolerance reached."; break; }
        if (r.iterations>=max_iter) break;
        double step=step_initial; bool accepted=false;
        for (int backtrack=0; backtrack<60; ++backtrack) {
            Vec trial=p; bool finite=true, changed=false;
            for (std::size_t j=0; j<p.size(); ++j) {
                trial[j]-=step*g[j]; finite=finite && std::isfinite(trial[j]);
                changed=changed || trial[j]!=p[j];
            }
            if (!changed) break;
            if (finite) {
                double ft=objective_value(f(trial)); ++r.evaluations;
                if (ft<=r.value-1e-4*step*gradient_norm*gradient_norm) {
                    p=trial; r.value=ft; accepted=true; break;
                }
            }
            step*=0.5;
        }
        if (!accepted) { r.status="Backtracking failed or the step stagnated."; break; }
        ++r.iterations;
    }
    r.point=p; return r;
}

int main() {
  try {
    Function f=[](const Vec& p) { return p[0]*p[0]+4*p[1]*p[1]+2*p[0]*p[1]; };
    auto grad=[](const Vec& p) { return Vec{2*p[0]+2*p[1],8*p[1]+2*p[0]}; };
    Vec p0={1.8,-1.6}; double step_initial=1,tol_grad=1e-6,gradient_norm;
    int max_iter=1000; long long gradient_calls;
    Result r=gradient_descent(f,grad,p0,step_initial,tol_grad,max_iter,gradient_norm,gradient_calls);
    print_result(r); std::cout << "Gradient norm: " << gradient_norm << "\nGradient calls: " << gradient_calls << '\n';
    return 0;
  } catch (const std::exception& e) {
    std::cerr << "Error: " << e.what() << '\n'; return 1;
  }
}