// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_optimization.cpp
// Optimization demo: Rosenbrock BFGS, Gaussian curve fitting, simplex LP
//
// Build:
//   cd build && cmake --build . --config Release --target example-optimization

#define _USE_MATH_DEFINES
#include <cmath>
#include <math/optimization/optimization_nd_impl.hpp>
#include <math/optimization/linear_programming_impl.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <iostream>
#include <iomanip>
#include <random>
#include <functional>

using namespace sangi;

int main() {
    // --- Demo 1: Rosenbrock BFGS ---
    {
        std::cout << "=== Rosenbrock minimization (BFGS) ===" << std::endl;

        auto rosenbrock = [](const Vector<double>& x) -> double {
            double a = 1.0 - x[0];
            double b = x[1] - x[0] * x[0];
            return a * a + 100.0 * b * b;
        };
        auto rosenbrock_grad = [](const Vector<double>& x) -> Vector<double> {
            Vector<double> g(2);
            g[0] = -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]);
            g[1] = 200.0 * (x[1] - x[0] * x[0]);
            return g;
        };

        Vector<double> x0(2);
        x0[0] = -1.0; x0[1] = -1.0;
        std::cout << "initial point: (" << x0[0] << ", " << x0[1] << ")" << std::endl;

        OptimizationOptions<double> opts;
        opts.tolerance = 1e-8;
        opts.max_iterations = 10000;

        auto result = bfgs_minimize<double, Vector<double>, Matrix<double>>(
            rosenbrock, rosenbrock_grad, x0, opts);

        std::cout << "iterations: " << result.iterations << std::endl;
        std::cout << std::fixed << std::setprecision(6);
        std::cout << "minimum at: (" << result.point[0] << ", " << result.point[1] << ")" << std::endl;
        std::cout << std::scientific << std::setprecision(2);
        std::cout << "f(x,y)    = " << result.value << std::endl;
        std::cout << "converged: " << (result.success ? "true" : "false") << std::endl;
    }

    // --- Demo 2: Gaussian curve fit ---
    {
        std::cout << "\n=== Gaussian curve fitting ===" << std::endl;

        double a_true = 3.0, b_true = 2.0, c_true = 0.5;

        std::mt19937 rng(42);
        std::normal_distribution<double> noise(0.0, 0.05);

        std::vector<double> xs, ys;
        for (int i = 0; i <= 40; ++i) {
            double x = i * 0.1;
            double y = a_true * std::exp(-(x - b_true) * (x - b_true)
                                          / (2.0 * c_true * c_true))
                       + noise(rng);
            xs.push_back(x);
            ys.push_back(y);
        }

        auto gaussian = [](double x, const Vector<double>& p) -> double {
            return p[0] * std::exp(-(x - p[1]) * (x - p[1])
                                    / (2.0 * p[2] * p[2]));
        };

        Vector<double> p0(3);
        p0[0] = 1.0; p0[1] = 1.0; p0[2] = 1.0;
        auto result = curveFit<double>(gaussian, xs, ys, p0);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "true parameters:    a=" << a_true
                  << ", b=" << b_true << ", c=" << c_true << std::endl;
        std::cout << "fitted parameters:  a=" << result.parameters[0]
                  << ", b=" << result.parameters[1]
                  << ", c=" << result.parameters[2] << std::endl;
        std::cout << "chi^2: " << result.chiSquared << std::endl;
    }

    // --- Demo 3: Linear programming (Simplex) ---
    {
        std::cout << "\n=== Linear programming (Simplex) ===" << std::endl;
        std::cout << "maximize: 3x + 5y" << std::endl;
        std::cout << "subject to:" << std::endl;
        std::cout << "  x +  y <= 4" << std::endl;
        std::cout << "  x + 3y <= 6" << std::endl;
        std::cout << "  x, y >= 0" << std::endl;

        LinearProgram<double> lp;
        lp.objective = {3.0, 5.0};
        lp.constraints.push_back(LinearConstraint<double>({1.0, 1.0}, ConstraintType::LessEqual, 4.0));
        lp.constraints.push_back(LinearConstraint<double>({1.0, 3.0}, ConstraintType::LessEqual, 6.0));

        auto result = simplexMaximize(lp);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "objective: " << result.objectiveValue << std::endl;
        std::cout << "x = " << result.solution[0] << ", y = " << result.solution[1] << std::endl;
    }

    return 0;
}