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

// vector_operations.hpp
//
// ベクトル操作の実装
//
// このファイルでは、MKL代数ライブラリのベクトル操作関数を定義します。
// これらの関数は、ベクトルの算術演算、数学関数、およびユーティリティを提供します。
//
// 主な機能:
// - ベクトルの算術演算（加算、減算、スカラー乗算など）
// - ベクトルの内積、外積、ノルム計算
// - 要素ごとの数学関数（絶対値、平方根など）
// - 計算ポリシーによる実装の切り替え

#ifndef CALX_VECTOR_OPERATIONS_HPP
#define CALX_VECTOR_OPERATIONS_HPP

#include "math/core/common.hpp"
#include "math/core/traits.hpp"
#include "math/concepts/algebraic_concepts.hpp"
#include "computation_policy.hpp"

#include <algorithm>
#include <cmath>
#include <functional>
#include <numeric>
#include <type_traits>

namespace calx {
    namespace detail {
        namespace vector_operations {

            //-----------------------------------------------------------------------------
            // 基本ベクトル操作
            //-----------------------------------------------------------------------------

            // ベクトル加算: result = a + b
            template<typename VectorA, typename VectorB, typename VectorResult,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::AdditiveGroup<typename VectorA::value_type>&&
            concepts::AdditiveGroup<typename VectorB::value_type>&&
                concepts::AdditiveGroup<typename VectorResult::value_type>
                void add(const VectorA& a, const VectorB& b, VectorResult& result) {
                ComputationPolicy::add(a, b, result);
            }

            // ベクトル減算: result = a - b
            template<typename VectorA, typename VectorB, typename VectorResult,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::AdditiveGroup<typename VectorA::value_type>&&
            concepts::AdditiveGroup<typename VectorB::value_type>&&
                concepts::AdditiveGroup<typename VectorResult::value_type>
                void subtract(const VectorA& a, const VectorB& b, VectorResult& result) {
                ComputationPolicy::subtract(a, b, result);
            }

            // スカラー乗算: result = alpha * a
            template<typename VectorA, typename Scalar, typename VectorResult,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::Ring<typename VectorA::value_type>&&
            concepts::Ring<Scalar>&&
                concepts::Ring<typename VectorResult::value_type>
                void scale(const VectorA& a, const Scalar& alpha, VectorResult& result) {
                // サイズチェック
                if (result.size() != a.size()) {
                    throw DimensionError("Vector scaling: size mismatch");
                }

                // まず結果ベクトルにaをコピー
                if (&result != &a) {
                    for (std::size_t i = 0; i < a.size(); ++i) {
                        result[i] = a[i];
                    }
                }

                // 次にスケーリングを適用
                ComputationPolicy::scale(result, alpha);
            }

            // ベクトルの符号反転: result = -a
            template<typename VectorA, typename VectorResult,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::AdditiveGroup<typename VectorA::value_type>&&
            concepts::AdditiveGroup<typename VectorResult::value_type>
                void negate(const VectorA& a, VectorResult& result) {
                scale(a, typename VectorA::value_type(-1), result);
            }

            // 内積: return = a・b
            template<typename VectorA, typename VectorB,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::Ring<typename VectorA::value_type>&&
            concepts::Ring<typename VectorB::value_type>
                auto dot_product(const VectorA& a, const VectorB& b)
                -> decltype(a[0] * b[0]) {
                return ComputationPolicy::dot_product(a, b);
            }

            // 要素ごとの積: result = a .* b（要素ごとの乗算）
            template<typename VectorA, typename VectorB, typename VectorResult,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::Ring<typename VectorA::value_type>&&
            concepts::Ring<typename VectorB::value_type>&&
                concepts::Ring<typename VectorResult::value_type>
                void element_wise_multiply(const VectorA& a, const VectorB& b, VectorResult& result) {
                ComputationPolicy::element_wise_multiply(a, b, result);
            }

            // 要素ごとの除算: result = a ./ b（要素ごとの除算）
            template<typename VectorA, typename VectorB, typename VectorResult>
                requires concepts::Field<typename VectorA::value_type>&&
            concepts::Field<typename VectorB::value_type>&&
                concepts::Field<typename VectorResult::value_type>
                void element_wise_divide(const VectorA& a, const VectorB& b, VectorResult& result) {
                const auto size = a.size();

                // サイズチェック
                if (b.size() != size || result.size() != size) {
                    throw DimensionError("Vector element-wise division: size mismatch");
                }

                // 要素ごとの除算（ゼロ除算チェック付き）
                for (std::size_t i = 0; i < size; ++i) {
                    if (approximately_equal(b[i], numeric_traits<typename VectorB::value_type>::zero())) {
                        throw MathError("Vector element-wise division: division by zero");
                    }
                    result[i] = a[i] / b[i];
                }
            }

            // axpy演算: y += alpha * x
            template<typename VectorY, typename Scalar, typename VectorX,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorY::value_type>>
                requires concepts::Ring<typename VectorY::value_type>&&
            concepts::Ring<Scalar>&&
                concepts::Ring<typename VectorX::value_type>
                void axpy(VectorY& y, const Scalar& alpha, const VectorX& x) {
                ComputationPolicy::axpy(y, alpha, x);
            }

            //-----------------------------------------------------------------------------
            // 数学関数
            //-----------------------------------------------------------------------------

            // ベクトルの各要素に対する関数適用: result = f(a)
            template<typename VectorA, typename VectorResult, typename Function>
            void apply_function(const VectorA& a, VectorResult& result, Function f) {
                const auto size = a.size();

                // サイズチェック
                if (result.size() != size) {
                    throw DimensionError("Vector function application: size mismatch");
                }

                // 各要素に関数を適用
                for (std::size_t i = 0; i < size; ++i) {
                    result[i] = f(a[i]);
                }
            }

            // ベクトルの絶対値（要素ごと）: result = |a|
            template<typename VectorA, typename VectorResult>
                requires requires(typename VectorA::value_type x) { { std::abs(x) } -> std::convertible_to<typename VectorResult::value_type>; }
            void abs(const VectorA& a, VectorResult& result) {
                apply_function(a, result, [](const auto& x) { return std::abs(x); });
            }

            // ベクトルの平方根（要素ごと）: result = √a
            template<typename VectorA, typename VectorResult>
                requires requires(typename VectorA::value_type x) { { std::sqrt(x) } -> std::convertible_to<typename VectorResult::value_type>; }
            void sqrt(const VectorA& a, VectorResult& result) {
                apply_function(a, result, [](const auto& x) {
                    if (x < numeric_traits<decltype(x)>::zero()) {
                        throw MathError("Vector sqrt: negative input");
                    }
                    return std::sqrt(x);
                    });
            }

            // ベクトルの二乗（要素ごと）: result = a²
            template<typename VectorA, typename VectorResult>
                requires concepts::Ring<typename VectorA::value_type>&&
            concepts::Ring<typename VectorResult::value_type>
                void square(const VectorA& a, VectorResult& result) {
                apply_function(a, result, [](const auto& x) { return x * x; });
            }

            // ベクトルの指数関数（要素ごと）: result = e^a
            template<typename VectorA, typename VectorResult>
                requires requires(typename VectorA::value_type x) { { std::exp(x) } -> std::convertible_to<typename VectorResult::value_type>; }
            void exp(const VectorA& a, VectorResult& result) {
                apply_function(a, result, [](const auto& x) { return std::exp(x); });
            }

            // ベクトルの対数関数（要素ごと）: result = log(a)
            template<typename VectorA, typename VectorResult>
                requires requires(typename VectorA::value_type x) { { std::log(x) } -> std::convertible_to<typename VectorResult::value_type>; }
            void log(const VectorA& a, VectorResult& result) {
                apply_function(a, result, [](const auto& x) {
                    if (x <= numeric_traits<decltype(x)>::zero()) {
                        throw MathError("Vector log: non-positive input");
                    }
                    return std::log(x);
                    });
            }

            // ベクトルのべき乗（要素ごと）: result = a^exponent
            template<typename VectorA, typename VectorResult, typename Exponent>
                requires requires(typename VectorA::value_type x, Exponent e) { { std::pow(x, e) } -> std::convertible_to<typename VectorResult::value_type>; }
            void pow(const VectorA& a, VectorResult& result, const Exponent& exponent) {
                apply_function(a, result, [&exponent](const auto& x) {
                    return std::pow(x, exponent);
                    });
            }

            //-----------------------------------------------------------------------------
            // ノルムと距離
            //-----------------------------------------------------------------------------

            // ベクトルのL1ノルム（絶対値の和）: ||a||₁
            template<typename VectorA>
                requires requires(typename VectorA::value_type x) { { std::abs(x) } -> std::convertible_to<typename VectorA::value_type>; }
            auto l1_norm(const VectorA& a) -> typename VectorA::value_type {
                using value_type = typename VectorA::value_type;
                value_type sum = numeric_traits<value_type>::zero();

                for (std::size_t i = 0; i < a.size(); ++i) {
                    sum += std::abs(a[i]);
                }

                return sum;
            }

            // ベクトルのL2ノルム（ユークリッドノルム）: ||a||₂
            template<typename VectorA,
                typename ComputationPolicy = detail::DefaultComputationPolicy<typename VectorA::value_type>>
                requires concepts::Ring<typename VectorA::value_type>&&
                requires(typename VectorA::value_type x) { { std::sqrt(x) } -> std::convertible_to<typename VectorA::value_type>; }
            auto l2_norm(const VectorA& a) -> typename VectorA::value_type {
                return ComputationPolicy::norm2(a);
            }

            // ベクトルのL∞ノルム（最大絶対値）: ||a||∞
            template<typename VectorA>
                requires requires(typename VectorA::value_type x) { { std::abs(x) } -> std::convertible_to<typename VectorA::value_type>; }
            auto linf_norm(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    return numeric_traits<typename VectorA::value_type>::zero();
                }

                using value_type = typename VectorA::value_type;
                value_type max_abs = std::abs(a[0]);

                for (std::size_t i = 1; i < a.size(); ++i) {
                    max_abs = std::max(max_abs, std::abs(a[i]));
                }

                return max_abs;
            }

            // ベクトル間のL2距離: ||a - b||₂
            template<typename VectorA, typename VectorB>
                requires concepts::AdditiveGroup<typename VectorA::value_type>&&
            concepts::AdditiveGroup<typename VectorB::value_type>&&
                requires(typename VectorA::value_type x) { { std::sqrt(x) } -> std::convertible_to<decltype(std::sqrt(std::declval<typename VectorA::value_type>()))>; }
            auto l2_distance(const VectorA& a, const VectorB& b)
                -> decltype(std::sqrt(std::declval<typename VectorA::value_type>())) {
                const auto size = a.size();

                // サイズチェック
                if (b.size() != size) {
                    throw DimensionError("Vector L2 distance: size mismatch");
                }

                using result_type = decltype(a[0] - b[0]);
                result_type sum_squares = numeric_traits<result_type>::zero();

                for (std::size_t i = 0; i < size; ++i) {
                    const result_type diff = a[i] - b[i];
                    sum_squares += diff * diff;
                }

                return std::sqrt(sum_squares);
            }

            //-----------------------------------------------------------------------------
            // 統計関数
            //-----------------------------------------------------------------------------

            // ベクトルの最小値
            template<typename VectorA>
                requires concepts::Comparable<typename VectorA::value_type>
            auto min(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    throw IndexError("Vector min: empty vector");
                }

                return *std::min_element(a.begin(), a.end());
            }

            // ベクトルの最大値
            template<typename VectorA>
                requires concepts::Comparable<typename VectorA::value_type>
            auto max(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    throw IndexError("Vector max: empty vector");
                }

                return *std::max_element(a.begin(), a.end());
            }

            // ベクトルの要素の和
            template<typename VectorA>
                requires concepts::AdditiveGroup<typename VectorA::value_type>
            auto sum(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    return numeric_traits<typename VectorA::value_type>::zero();
                }

                using value_type = typename VectorA::value_type;
                return std::accumulate(a.begin(), a.end(), numeric_traits<value_type>::zero());
            }

            // ベクトルの要素の積
            template<typename VectorA>
                requires concepts::MultiplicativeMonoid<typename VectorA::value_type>
            auto product(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    return numeric_traits<typename VectorA::value_type>::one();
                }

                using value_type = typename VectorA::value_type;
                return std::accumulate(a.begin(), a.end(), numeric_traits<value_type>::one(),
                    std::multiplies<value_type>());
            }

            // ベクトルの算術平均
            template<typename VectorA>
                requires concepts::Field<typename VectorA::value_type>
            auto mean(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty()) {
                    throw IndexError("Vector mean: empty vector");
                }

                using value_type = typename VectorA::value_type;
                const value_type total = sum(a);
                return total / static_cast<value_type>(a.size());
            }

            // ベクトルの分散
            template<typename VectorA>
                requires concepts::Field<typename VectorA::value_type>
            auto variance(const VectorA& a) -> typename VectorA::value_type {
                if (a.empty() || a.size() == 1) {
                    throw IndexError("Vector variance: insufficient data");
                }

                const auto avg = mean(a);
                using value_type = typename VectorA::value_type;
                value_type sum_sq_diff = numeric_traits<value_type>::zero();

                for (std::size_t i = 0; i < a.size(); ++i) {
                    const value_type diff = a[i] - avg;
                    sum_sq_diff += diff * diff;
                }

                return sum_sq_diff / static_cast<value_type>(a.size() - 1);
            }

            // ベクトルの標準偏差
            template<typename VectorA>
                requires concepts::Field<typename VectorA::value_type>&&
                requires(typename VectorA::value_type x) { { std::sqrt(x) } -> std::convertible_to<typename VectorA::value_type>; }
            auto standard_deviation(const VectorA& a) -> typename VectorA::value_type {
                return std::sqrt(variance(a));
            }

            // 数値状態を含むベクトルの内積計算
            template<typename VectorA, typename VectorB>
                requires HasNumericState<typename VectorA::value_type>&&
            HasNumericState<typename VectorB::value_type>
                auto dot_product_with_state(const VectorA& a, const VectorB& b)
                -> typename VectorA::value_type {
                const auto size = a.size();

                // サイズチェック
                if (b.size() != size) {
                    throw DimensionError("Vector dot product with state: size mismatch");
                }

                // 特殊状態のチェック
                for (std::size_t i = 0; i < size; ++i) {
                    if (numeric_state_traits<typename VectorA::value_type>::isNaN(a[i]) ||
                        numeric_state_traits<typename VectorB::value_type>::isNaN(b[i])) {
                        // NaNが含まれる場合はNaNを返す
                        return numeric_traits<typename VectorA::value_type>::quiet_NaN();
                    }
                    if (numeric_state_traits<typename VectorA::value_type>::isInfinite(a[i]) ||
                        numeric_state_traits<typename VectorB::value_type>::isInfinite(b[i])) {
                        // 無限大の場合は個別に処理
                        // 例: ∞ * 0 = NaN, ∞ * 非ゼロ = ∞
                        if ((numeric_state_traits<typename VectorA::value_type>::isInfinite(a[i]) &&
                            b[i] == numeric_traits<typename VectorB::value_type>::zero()) ||
                            (numeric_state_traits<typename VectorB::value_type>::isInfinite(b[i]) &&
                                a[i] == numeric_traits<typename VectorA::value_type>::zero())) {
                            return numeric_traits<typename VectorA::value_type>::quiet_NaN();
                        }
                        // その他の無限大の場合は無限大を返す
                        // 符号は計算結果によって異なる
                        return numeric_traits<typename VectorA::value_type>::infinity();
                    }
                }

                // 通常の内積計算
                using result_type = typename VectorA::value_type;
                result_type sum = numeric_traits<result_type>::zero();

                for (std::size_t i = 0; i < size; ++i) {
                    sum += a[i] * b[i];
                }

                return sum;
            }

            // 数値状態を含むベクトルのL2ノルム計算
            template<typename VectorA>
                requires HasNumericState<typename VectorA::value_type>
            auto l2_norm_with_state(const VectorA& a) -> typename VectorA::value_type {
                using value_type = typename VectorA::value_type;

                // 特殊状態のチェック
                for (std::size_t i = 0; i < a.size(); ++i) {
                    if (numeric_state_traits<value_type>::isNaN(a[i])) {
                        // NaNが含まれる場合はNaNを返す
                        return numeric_traits<value_type>::quiet_NaN();
                    }
                    if (numeric_state_traits<value_type>::isInfinite(a[i])) {
                        // 無限大の要素が含まれる場合は無限大を返す
                        return numeric_traits<value_type>::infinity();
                    }
                }

                // 通常のノルム計算
                return l2_norm(a);
            }

} // namespace vector_operations
} // namespace detail
} // namespace calx

#endif // CALX_VECTOR_OPERATIONS_HPP