// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_doc_optimization_samples.cpp
// Verification harness for samples in api/Optimization.html (ja+en).

#include <math/optimization/optimization_nd.hpp>
#include <math/optimization/optimization_nd_impl.hpp>
#include <math/optimization/linear_programming.hpp>
#include <math/optimization/linear_programming_impl.hpp>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

using namespace sangi;

int main() {
    std::cout << std::setprecision(15);

    // ---- Rosenbrock BFGS -----------------------------------------
    {
        std::cout << "[Rosenbrock BFGS]\n";
        auto f = [](const Vector<double>& x) {
            double a = 1.0 - x[0];
            double b = x[1] - x[0] * x[0];
            return a * a + 100.0 * b * b;
        };
        auto df = [](const Vector<double>& x) -> Vector<double> {
            double dx = -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]);
            double dy = 200.0 * (x[1] - x[0] * x[0]);
            return Vector<double>({dx, dy});
        };
        Vector<double> x0({-1.0, 1.0});
        auto result = bfgs_minimize<double, Vector<double>, Matrix<double>>(f, df, x0);
        std::cout << "  minimizer = (" << result.point[0] << ", " << result.point[1] << ")\n";
        std::cout << "  min value = " << result.value << "\n";
        std::cout << "  iters     = " << result.iterations << "\n";
        std::cout << "  success= " << std::boolalpha << result.success << "\n";
    }

    // ---- curveFit y = a*exp(-b*x) --------------------------------
    {
        std::cout << "[curveFit y = a*exp(-b*x)]\n";
        auto model = [](double x, const Vector<double>& p) {
            return p[0] * std::exp(-p[1] * x);
        };
        std::vector<double> xdata = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5};
        std::vector<double> ydata = {5.0, 3.1, 1.8, 1.2, 0.7, 0.4};
        Vector<double> p0({4.0, 0.5});
        auto result = curveFit<double>(model, xdata, ydata, p0);
        std::cout << "  a = " << result.parameters[0]
                  << " +/- " << result.standardErrors[0] << "\n";
        std::cout << "  b = " << result.parameters[1]
                  << " +/- " << result.standardErrors[1] << "\n";
        std::cout << "  chi^2/dof = " << result.reducedChiSquared << "\n";
        std::cout << "  iters     = " << result.iterations << "\n";
    }

    // ---- Linear programming --------------------------------------
    {
        std::cout << "[Linear programming]\n";
        LinearProgram<double> lp;
        lp.objective = {1.0, 2.0};
        lp.constraints = {
            LinearConstraint<double>({1.0, 1.0}, ConstraintType::LessEqual, 4.0),
            LinearConstraint<double>({1.0, 0.0}, ConstraintType::LessEqual, 3.0),
            LinearConstraint<double>({0.0, 1.0}, ConstraintType::LessEqual, 3.0),
        };
        auto result = simplexMaximize(lp);
        if (result.status == LPStatus::Optimal) {
            std::cout << "  optimum = " << result.objectiveValue << "\n";
            std::cout << "  x = " << result.solution[0]
                      << ", y = " << result.solution[1] << "\n";
        } else {
            std::cout << "  not optimal: status=" << static_cast<int>(result.status) << "\n";
        }
    }

    return 0;
}