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

// vector_space_utilities.hpp
#ifndef SANGI_VECTOR_SPACE_UTILITIES_HPP
#define SANGI_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 <cassert>

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

namespace sangi {

    /**
     * @brief Computes a linear combination in a vector space
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param vectors list of vectors
     * @param coefficients list of coefficients
     * @return the result of the linear combination
     * @throws std::invalid_argument when the number of vectors does not match the number of coefficients
     */
    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);
        // Initialize with the zero vector
        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 Gram-Schmidt orthogonalization
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param vectors list of vectors to orthogonalize
     * @return list of orthogonalized vectors
     */
    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;

            // Compute the orthogonal component against the already-orthogonalized vectors
            for (const auto& w : orthogonal) {
                T ww = dot(w, w);
                if (ww == T(0)) continue;  // avoid division by zero (epsilon=0 for Rational, etc.)
                T projection = dot(v, w) / ww;
                u -= w * projection;
            }

            // Add only nonzero vectors
            if (norm(u) > std::numeric_limits<T>::epsilon()) {
                orthogonal.push_back(u);
            }
        }

        return orthogonal;
    }

    /**
     * @brief Normalizes a vector (makes it a unit vector)
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param v the vector to normalize
     * @return the normalized vector
     */
    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 Computes the projection of a vector
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param v the vector to project
     * @param onto the vector to project onto
     * @return the projection of v onto 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 Determines whether a set of vectors is linearly independent
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @tparam M matrix type
     * @param vectors list of vectors
     * @return true if linearly independent, false otherwise
     */
    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;  // the empty set is linearly independent
        }

        if (vectors.size() > vectors[0].size()) {
            return false;  // linearly dependent if the number of vectors exceeds the dimension
        }

        // Build the Gram matrix (matrix of inner products)
        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]);
            }
        }

        // Compute the determinant
        T det = gram.determinant();
        return std::abs(det) > std::numeric_limits<T>::epsilon();
    }

    /**
     * @brief Solves linear regression via least squares
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @tparam M matrix type
     * @param X matrix of explanatory variables (each row is a sample, each column a feature)
     * @param y vector of response variables
     * @return vector of regression coefficients
     * @throws DimensionError when the sizes of the matrix and vector do not match
     * @throws LinearAlgebraError when the matrix is singular
     */
    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()) {
            assert(false && "DimensionError: Number of samples (rows in X) must match length of y");
            throw DimensionError("Number of samples (rows in X) must match length of y");
        }

        // Solve the normal equations: 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;

        // Compute the inverse to obtain the coefficients
        M X_transpose_X_inv = algorithms::lu_inverse(X_transpose_X);
        return X_transpose_X_inv * X_transpose_y;
    }

    /**
     * @brief Completes a basis of a vector space
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param partial_basis list of vectors forming a partial basis
     * @param dimension dimension of the full space
     * @return list of the completed basis vectors
     */
    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()) {
            // For an empty partial basis, return the standard basis
            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;
        }

        // Orthonormalize the partial basis
        std::vector<V> orthogonal_basis = gram_schmidt(partial_basis);
        std::vector<V> result = orthogonal_basis;

        // Verify that all vectors have the same size
        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");
        }

        // Generate the remaining basis vectors
        while (result.size() < dimension) {
            // Generate a random vector
            V random_vector(vector_size);
            for (std::size_t i = 0; i < vector_size; ++i) {
                // Simple pseudo-random value (between -1 and 1)
                random_vector[i] = T(2) * T(std::rand()) / T(RAND_MAX) - T(1);
            }

            // Orthogonalize against the existing basis
            for (const auto& basis_vector : result) {
                T projection = dot(random_vector, basis_vector) / dot(basis_vector, basis_vector);
                random_vector -= basis_vector * projection;
            }

            // Normalize and add to the basis (only for nonzero vectors)
            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 Orthonormalizes a basis using QR decomposition
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @tparam M matrix type
     * @param basis list of basis vectors to orthonormalize
     * @return a pair of the orthonormal basis and its transformation matrix
     */
    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) };
        }

        // Arrange the vectors as columns of a matrix
        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];
            }
        }

        // Perform the QR decomposition
        auto qr_result = algorithms::qr_decomposition(A);
        M Q = qr_result.first;
        M R = qr_result.second;

        // Extract the orthonormal basis
        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 Determines whether two vectors are parallel
     * @tparam T scalar type (field)
     * @tparam V vector type
     * @param v1 the first vector
     * @param v2 the second vector
     * @param tolerance tolerance for the test
     * @return true if parallel, false otherwise
     */
    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);

        // Treat a zero vector as parallel to any vector
        if (norm1 < tolerance || norm2 < tolerance) {
            return true;
        }

        // Normalize, then compute the dot product
        V unit1 = v1 / norm1;
        V unit2 = v2 / norm2;

        // The closer the absolute value of the dot product is to 1, the more parallel
        T dot_product = std::abs(dot(unit1, unit2));
        return std::abs(dot_product - T(1)) < tolerance;
    }

} // namespace sangi

#endif // SANGI_VECTOR_SPACE_UTILITIES_HPP