// ch13_polynomial_models.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_polynomial_models.cpp -o ch13_polynomial_models
// 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';
}

struct Polynomial { Vec coefficients; double center,scale; };
// Householder QR least-squares solve; does not form normal equations.
Vec qr_least_squares(Mat A, Vec b) {
    const int m=static_cast<int>(A.size()), n=static_cast<int>(A[0].size());
    require(m>=n && b.size()==A.size(), "Invalid least-squares dimensions.");
    for (int k=0;k<n;++k) {
        double norm=0;
        for (int i=k;i<m;++i) norm=std::hypot(norm,A[i][k]);
        require(std::isfinite(norm) && norm>1e-12,"Rank-deficient or poorly scaled polynomial design.");
        double alpha=-std::copysign(norm,A[k][k]); Vec v(m-k);
        for (int i=k;i<m;++i) v[i-k]=A[i][k];
        v[0]-=alpha; double vv=std::inner_product(v.begin(),v.end(),v.begin(),0.0);
        for (int j=k;j<n;++j) {
            double dot=0;
            for (int i=k;i<m;++i) dot+=v[i-k]*A[i][j];
            for (int i=k;i<m;++i) A[i][j]-=2*v[i-k]*dot/vv;
        }
        double dot=0;
        for (int i=k;i<m;++i) dot+=v[i-k]*b[i];
        for (int i=k;i<m;++i) b[i]-=2*v[i-k]*dot/vv;
        A[k][k]=alpha;
    }
    Vec coef(n);
    for (int i=n-1;i>=0;--i) {
        double sum=b[i]; for (int j=i+1;j<n;++j) sum-=A[i][j]*coef[j]; coef[i]=sum/A[i][i];
    }
    finite_vector(coef); return coef;
}
Polynomial polynomial_fit(const Vec& x, const Vec& y, int degree) {
    finite_vector(x); finite_vector(y);
    require(x.size()==y.size() && x.size()>=2 && degree>=0 && static_cast<std::size_t>(degree)<x.size(),"Invalid polynomial data or degree.");
    Vec unique=x; std::sort(unique.begin(),unique.end());
    auto end=std::unique(unique.begin(),unique.end());
    require(static_cast<std::size_t>(degree)<static_cast<std::size_t>(end-unique.begin()),"Too few distinct x values.");
    Polynomial model; model.center=std::accumulate(x.begin(),x.end(),0.0)/x.size();
    double sum=0; for (double v:x) sum+=(v-model.center)*(v-model.center);
    model.scale=std::sqrt(sum/(x.size()-1));
    require(std::isfinite(model.scale) && model.scale>0,"No usable x spread.");
    Mat A(x.size(),Vec(degree+1,1.0));
    for (std::size_t i=0;i<x.size();++i) for (int j=1;j<=degree;++j) A[i][j]=A[i][j-1]*(x[i]-model.center)/model.scale;
    model.coefficients=qr_least_squares(A,y); return model;
}
double polynomial_predict(const Polynomial& model, double x) {
    double z=(x-model.center)/model.scale, value=0;
    for (auto i=model.coefficients.rbegin();i!=model.coefficients.rend();++i) value=value*z+*i;
    return value;
}

int main() {
  try {
    const unsigned int seed = 7;
    Random rng(seed);
    // Alternative: comment out the previous line and seed from system entropy.
    // Random rng(std::random_device{}());
    auto truth=[](double x) { return std::sin(2*pi*x)+.5*x; };
    double noise_sigma=.18;
    Vec x=linspace(0,1,18),y(x.size()),x_validation=linspace(0,1,101),y_validation(x_validation.size());
    for (std::size_t i=0;i<x.size();++i) y[i]=truth(x[i])+noise_sigma*rng.z();
    for (std::size_t i=0;i<x_validation.size();++i) y_validation[i]=truth(x_validation[i])+noise_sigma*rng.z();
    std::vector<int> degrees={1,5,14};
    // Substitute separate training/validation data. Do not fit validation data.
    // Reserve fresh test data for assessing a model after selection.
    std::cout << std::setprecision(10) << "Degree   Training RMSE    Validation RMSE\n";
    for (int degree:degrees) {
        Polynomial model=polynomial_fit(x,y,degree); double train=0,valid=0;
        for (std::size_t i=0;i<x.size();++i) { double e=y[i]-polynomial_predict(model,x[i]); train+=e*e; }
        for (std::size_t i=0;i<x_validation.size();++i) { double e=y_validation[i]-polynomial_predict(model,x_validation[i]); valid+=e*e; }
        std::cout << degree << "        " << std::sqrt(train/x.size()) << "        " << std::sqrt(valid/x_validation.size()) << '\n';
    }
    return 0;
  } catch (const std::exception& e) {
    std::cerr << "Error: " << e.what() << '\n'; return 1;
  }
}