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

// ConvexOptimization.hpp — 凸最適化ソルバー
//
// ADMM, 近接勾配法 (ISTA/FISTA)。

#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 Vector<T>& v, T lambda);

    /// 非負射影: prox_{I_{x≥0}}(v) = max(v, 0)
    template<typename T>
    Vector<T> proxNonneg(const Vector<T>& v);

    // ================================================================
    // ISTA / FISTA (近接勾配法)
    // ================================================================

    /// FISTA の結果
    template<typename T>
    struct FISTAResult {
        Vector<T> x;        // 解
        T objective;         // 最終目的関数値
        size_t iterations;   // 反復回数
    };

    /// FISTA (Fast Iterative Shrinkage-Thresholding Algorithm)
    /// minimize f(x) + g(x)
    /// f: 微分可能な凸関数, ∇f が Lipschitz 連続 (Lipschitz 定数 L)
    /// g: 近接作用素が計算可能な凸関数
    ///
    /// @param gradF   f の勾配関数
    /// @param proxG   g の近接作用素
    /// @param x0      初期値
    /// @param L       Lipschitz 定数 (ステップサイズ = 1/L)
    /// @param maxIter 最大反復回数
    /// @param tol     収束判定閾値
    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 Matrix<T>& A, const Vector<T>& b,
                               T lambda, size_t maxIter = 1000, T tol = T{1e-8});

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

    /// ADMM の結果
    template<typename T>
    struct ADMMResult {
        Vector<T> x;         // 解
        size_t iterations;
    };

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

    // 実装は ConvexOptimization_impl.hpp に分離。

} // namespace sangi

#endif // SANGI_MATH_OPTIMIZATION_CONVEX_HPP