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

// ConvexOptimization.hpp — convex optimization solvers
//
// ADMM, proximal gradient methods (ISTA/FISTA), interior-point method (QP).

#ifndef SANGI_MATH_OPTIMIZATION_CONVEX_HPP
#define SANGI_MATH_OPTIMIZATION_CONVEX_HPP

#include <math/core/matrix.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <numeric>
#include <vector>

namespace sangi {

    // ================================================================
    // Proximal Operators
    // ================================================================

    /// Soft thresholding: prox_{λ||·||_1}(v) = sign(v) * max(|v| - λ, 0)
    template<typename T>
    Vector<T> proxL1(const BaseVector<T>& v, T lambda);

    /// Nonnegative projection: prox_{I_{x≥0}}(v) = max(v, 0)
    template<typename T>
    Vector<T> proxNonneg(const BaseVector<T>& v);

    // ================================================================
    // ISTA / FISTA (proximal gradient methods)
    // ================================================================

    /// FISTA result
    template<typename T>
    struct FISTAResult {
        Vector<T> x;        // solution
        T objective;         // final objective value
        size_t iterations;   // number of iterations
    };

    /// FISTA (Fast Iterative Shrinkage-Thresholding Algorithm)
    /// minimize f(x) + g(x)
    /// f: differentiable convex function, ∇f Lipschitz continuous (Lipschitz constant L)
    /// g: convex function with a computable proximal operator
    ///
    /// @param gradF   gradient function of f
    /// @param proxG   proximal operator of g
    /// @param x0      initial value
    /// @param L       Lipschitz constant (step size = 1/L)
    /// @param maxIter maximum number of iterations
    /// @param tol     convergence threshold
    template<typename T>
    FISTAResult<T> fista(
        std::function<Vector<T>(const Vector<T>&)> gradF,
        std::function<Vector<T>(const Vector<T>&, T)> proxG,
        const Vector<T>& x0,
        T L,
        size_t maxIter = 1000, T tol = T{1e-8});

    /// Lasso via FISTA: minimize 0.5 ||Ax - b||² + λ||x||₁
    template<typename T>
    FISTAResult<T> lassoFISTA(const BaseMatrix<T>& A, const BaseVector<T>& b,
                               T lambda, size_t maxIter = 1000, T tol = T{1e-8});

    // ================================================================
    // ADMM (Alternating Direction Method of Multipliers)
    // ================================================================

    /// ADMM result
    template<typename T>
    struct ADMMResult {
        Vector<T> x;         // solution
        size_t iterations;
    };

    /// ADMM for Lasso: minimize 0.5 ||Ax - b||² + λ||z||₁
    /// split: x = z, penalty ρ
    template<typename T>
    ADMMResult<T> admmLasso(const BaseMatrix<T>& A, const BaseVector<T>& b,
                             T lambda, T rho = T{1},
                             size_t maxIter = 1000, T tol = T{1e-6});

    // ================================================================
    // Interior-point method — quadratic programming (QP)
    // ================================================================

    /// QP result
    template<typename T>
    struct QPResult {
        Vector<T> x;        // solution
        T objective;         // objective value
        size_t iterations;
    };

    /// QP: minimize 0.5 x^T Q x + c^T x  subject to  Ax ≤ b
    /// Mehrotra predictor-corrector method (simplified version)
    template<typename T>
    QPResult<T> solveQP(const BaseMatrix<T>& Q, const BaseVector<T>& c,
                         const BaseMatrix<T>& A, const BaseVector<T>& b,
                         size_t maxIter = 100, T tol = T{1e-8});

    // Implementation is separated into ConvexOptimization_impl.hpp.

} // namespace sangi

#endif // SANGI_MATH_OPTIMIZATION_CONVEX_HPP