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

// vector_space_utilities.hpp
#ifndef CALX_VECTOR_SPACE_UTILITIES_HPP
#define CALX_VECTOR_SPACE_UTILITIES_HPP

#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/concepts/algebraic_concepts.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <math/linalg/decomposition.hpp>

#include <vector>
#include <stdexcept>
#include <limits>
#include <cmath>
#include <cstdlib>

namespace calx {

    /**
     * @brief ベクトル空間の線形結合を計算
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param vectors ベクトルのリスト
     * @param coefficients 係数のリスト
     * @return 線形結合の結果
     * @throws std::invalid_argument ベクトル数と係数数が一致しない場合
     */
    template<typename T, typename V>
        requires concepts::Field<T>&& concepts::VectorOf<V, T>
    V linear_combination(
        const std::vector<V>& vectors,
        const std::vector<T>& coefficients) {

        if (vectors.size() != coefficients.size()) {
            throw std::invalid_argument("Number of vectors must match number of coefficients");
        }

        if (vectors.empty()) {
            throw std::invalid_argument("Empty vector list");
        }

        std::size_t size = vectors[0].size();
        V result(size);
        // ゼロベクトルで初期化
        result.zero();

        for (std::size_t i = 0; i < vectors.size(); ++i) {
            if (vectors[i].size() != size) {
                throw std::invalid_argument("All vectors must have the same size");
            }
            // result += vectors[i] * coefficients[i]
            V term = vectors[i];
            term *= coefficients[i];
            result += term;
        }

        return result;
    }

    /**
     * @brief グラム・シュミット直交化法
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param vectors 直交化するベクトルのリスト
     * @return 直交化されたベクトルのリスト
     */
    template<typename T, typename V>
        requires concepts::Field<T>&& concepts::VectorOf<V, T>&&
        requires(V a, V b) { { dot(a, b) } -> std::convertible_to<T>; }
    std::vector<V> gram_schmidt(std::vector<V> vectors) {
        if (vectors.empty()) {
            return {};
        }

        std::vector<V> orthogonal;
        orthogonal.reserve(vectors.size());

        for (const auto& v : vectors) {
            V u = v;

            // 既に直交化されたベクトルに対して直交成分を計算
            for (const auto& w : orthogonal) {
                T ww = dot(w, w);
                if (ww == T(0)) continue;  // ゼロ除算回避 (Rational 等で epsilon=0)
                T projection = dot(v, w) / ww;
                u -= w * projection;
            }

            // 非ゼロベクトルのみを追加
            if (norm(u) > std::numeric_limits<T>::epsilon()) {
                orthogonal.push_back(u);
            }
        }

        return orthogonal;
    }

    /**
     * @brief ベクトルを正規化（単位ベクトル化）する
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param v 正規化するベクトル
     * @return 正規化されたベクトル
     */
    template<typename T, typename V>
        requires concepts::Field<T>&& concepts::VectorOf<V, T>&&
        requires(V a) { { norm(a) } -> std::convertible_to<T>; }
    V normalize(const V& v) {
        T norm_val = norm(v);
        if (norm_val == T(0) || norm_val < std::numeric_limits<T>::epsilon()) {
            throw std::invalid_argument("Cannot normalize a zero vector");
        }
        return v / norm_val;
    }

    /**
     * @brief ベクトルの射影を計算
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param v 射影するベクトル
     * @param onto 射影先のベクトル
     * @return v の onto への射影
     */
    template<typename T, typename V>
        requires concepts::Field<T>&& concepts::VectorOf<V, T>&&
        requires(V a, V b) { { dot(a, b) } -> std::convertible_to<T>; }
    V projection(const V& v, const V& onto) {
        T onto_norm_squared = dot(onto, onto);
        if (onto_norm_squared == T(0) || onto_norm_squared < std::numeric_limits<T>::epsilon()) {
            throw std::invalid_argument("Cannot project onto a zero vector");
        }
        return onto * (dot(v, onto) / onto_norm_squared);
    }

    /**
     * @brief ベクトルセットが線形独立かどうかを判定
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @tparam M 行列型
     * @param vectors ベクトルのリスト
     * @return 線形独立ならtrue、そうでなければfalse
     */
    template<typename T, typename V, typename M>
        requires concepts::Field<T>&&
    concepts::VectorOf<V, T>&&
        concepts::MatrixOf<M, T>&&
        requires(V a, V b, M m) {
            { dot(a, b) } -> std::convertible_to<T>;
            { m.determinant() } -> std::convertible_to<T>;
    }
    bool is_linearly_independent(const std::vector<V>& vectors) {
        if (vectors.empty()) {
            return true;  // 空集合は線形独立
        }

        if (vectors.size() > vectors[0].size()) {
            return false;  // ベクトル数が次元を超える場合は線形従属
        }

        // グラム行列（内積行列）を作成
        std::size_t n = vectors.size();
        M gram(n, n);
        for (std::size_t i = 0; i < n; ++i) {
            for (std::size_t j = 0; j < n; ++j) {
                gram(i, j) = dot(vectors[i], vectors[j]);
            }
        }

        // 行列式の計算
        T det = gram.determinant();
        return std::abs(det) > std::numeric_limits<T>::epsilon();
    }

    /**
     * @brief 最小二乗法で線形回帰を解く
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @tparam M 行列型
     * @param X 説明変数の行列（各行がサンプル、各列が特徴量）
     * @param y 目的変数のベクトル
     * @return 回帰係数ベクトル
     * @throws DimensionError 行列とベクトルのサイズが一致しない場合
     * @throws LinearAlgebraError 行列が特異な場合
     */
    template<typename T, typename V, typename M>
        requires concepts::Field<T>&&
    concepts::VectorOf<V, T>&&
        concepts::MatrixOf<M, T>&&
        requires(M m, V v) {
            { m.transpose() } -> std::convertible_to<M>;
            { m* v } -> std::convertible_to<V>;
    }
    V linear_regression(const M& X, const V& y) {
        if (X.rows() != y.size()) {
            throw DimensionError("Number of samples (rows in X) must match length of y");
        }

        // 正規方程式を解く: beta = (X^T * X)^(-1) * X^T * y
        M X_transpose = X.transpose();
        M X_transpose_X = X_transpose * X;
        V X_transpose_y = X_transpose * y;

        // 逆行列を計算して係数を求める
        M X_transpose_X_inv = algorithms::lu_inverse(X_transpose_X);
        return X_transpose_X_inv * X_transpose_y;
    }

    /**
     * @brief ベクトル空間の基底を完成させる
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param partial_basis 部分基底となるベクトルのリスト
     * @param dimension 全空間の次元
     * @return 完成した基底ベクトルのリスト
     */
    template<typename T, typename V>
        requires concepts::Field<T>&&
    concepts::VectorOf<V, T>&&
        requires(V a, V b) {
            { dot(a, b) } -> std::convertible_to<T>;
            { norm(a) } -> std::convertible_to<T>;
    }
    std::vector<V> complete_basis(const std::vector<V>& partial_basis, std::size_t dimension) {
        if (partial_basis.empty()) {
            // 空の部分基底の場合は標準基底を返す
            std::vector<V> standard_basis;
            standard_basis.reserve(dimension);

            for (std::size_t i = 0; i < dimension; ++i) {
                V e(dimension, T(0));
                e[i] = T(1);
                standard_basis.push_back(e);
            }

            return standard_basis;
        }

        // 部分基底の正規直交化
        std::vector<V> orthogonal_basis = gram_schmidt(partial_basis);
        std::vector<V> result = orthogonal_basis;

        // すべてのベクトルのサイズが一致することを確認
        std::size_t vector_size = orthogonal_basis[0].size();
        if (vector_size < dimension) {
            throw std::invalid_argument("Vector dimension is smaller than the specified space dimension");
        }

        // 残りの基底ベクトルを生成
        while (result.size() < dimension) {
            // ランダムベクトルを生成
            V random_vector(vector_size);
            for (std::size_t i = 0; i < vector_size; ++i) {
                // 単純な擬似乱数（-1と1の間）
                random_vector[i] = T(2) * T(std::rand()) / T(RAND_MAX) - T(1);
            }

            // 既存の基底に対して直交化
            for (const auto& basis_vector : result) {
                T projection = dot(random_vector, basis_vector) / dot(basis_vector, basis_vector);
                random_vector -= basis_vector * projection;
            }

            // 正規化して基底に追加（非ゼロベクトルの場合のみ）
            T norm_val = norm(random_vector);
            if (norm_val > std::numeric_limits<T>::epsilon() * T(100)) {
                result.push_back(random_vector / norm_val);
            }
        }

        return result;
    }

    /**
     * @brief QR分解を用いて基底の正規直交化を行う
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @tparam M 行列型
     * @param basis 正規直交化する基底ベクトルのリスト
     * @return 正規直交基底とその変換行列のペア
     */
    template<typename T, typename V, typename M>
        requires concepts::Field<T>&&
    concepts::VectorOf<V, T>&&
        concepts::MatrixOf<M, T>&&
        requires {
        typename M::size_type;
        typename V::size_type;
    }
    std::pair<std::vector<V>, M> orthonormalize_basis(const std::vector<V>& basis) {
        if (basis.empty()) {
            return { basis, M(0, 0) };
        }

        // ベクトルを行列の列として配置
        std::size_t n = basis.size();
        std::size_t m = basis[0].size();

        M A(m, n);
        for (std::size_t j = 0; j < n; ++j) {
            for (std::size_t i = 0; i < m; ++i) {
                A(i, j) = basis[j][i];
            }
        }

        // QR分解を実行
        auto qr_result = algorithms::qr_decomposition(A);
        M Q = qr_result.first;
        M R = qr_result.second;

        // 正規直交基底を抽出
        std::vector<V> orthonormal_basis;
        orthonormal_basis.reserve(n);

        for (std::size_t j = 0; j < n; ++j) {
            V q(m);
            for (std::size_t i = 0; i < m; ++i) {
                q[i] = Q(i, j);
            }
            orthonormal_basis.push_back(q);
        }

        return { orthonormal_basis, R };
    }

    /**
     * @brief ふたつのベクトルが平行かどうかを判定
     * @tparam T スカラー型（体）
     * @tparam V ベクトル型
     * @param v1 1つ目のベクトル
     * @param v2 2つ目のベクトル
     * @param tolerance 判定の許容誤差
     * @return 平行ならtrue、そうでなければfalse
     */
    template<typename T, typename V>
        requires concepts::Field<T>&&
    concepts::VectorOf<V, T>&&
        requires(V a, V b) {
            { norm(a) } -> std::convertible_to<T>;
            { dot(a, b) } -> std::convertible_to<T>;
    }
    bool are_parallel(
        const V& v1,
        const V& v2,
        T tolerance = std::numeric_limits<T>::epsilon() * T(100)) {

        T norm1 = norm(v1);
        T norm2 = norm(v2);

        // ゼロベクトルは任意のベクトルと平行と見なす
        if (norm1 < tolerance || norm2 < tolerance) {
            return true;
        }

        // 正規化してから内積を計算
        V unit1 = v1 / norm1;
        V unit2 = v2 / norm2;

        // 内積の絶対値が1に近いほど平行
        T dot_product = std::abs(dot(unit1, unit2));
        return std::abs(dot_product - T(1)) < tolerance;
    }

} // namespace calx

#endif // CALX_VECTOR_SPACE_UTILITIES_HPP