StatisticalAnalyzer.hpp

Go to the documentation of this file.

#pragma once
 
#include "MayaFlux/Transitive/Reflect/EnumReflect.hpp"
#include "MayaFlux/Yantra/OperationSpec/OperationHelper.hpp"
 
#include "UniversalAnalyzer.hpp"
 
#include "MayaFlux/Kinesis/Discrete/Analysis.hpp"
#include <Eigen/Dense>
 
/**
 * @file StatisticalAnalyzer_new.hpp
 * @brief Span-based statistical analysis for digital signals in Maya Flux
 *
 * Defines the StatisticalAnalyzer using the new UniversalAnalyzer framework with
 * zero-copy span processing and automatic structure handling via OperationHelper.
 * This analyzer extracts statistical features from digital data streams with multiple
 * computation methods and flexible output configurations.
 *
 * Key Features:
 * - **Zero-copy processing:** Uses spans for maximum efficiency
 * - **Template-flexible I/O:** Instance defines input/output types at construction
 * - **Multiple statistical methods:** Mean, variance, std dev, skewness, kurtosis, percentiles, etc.
 * - **Parallel processing:** Utilizes std::execution for performance
 * - **Cross-modal support:** Works on any numeric data stream - truly digital-first
 * - **Statistical classification:** Maps values to qualitative levels (outliers, normal, etc.)
 * - **Automatic data handling:** OperationHelper manages all extraction/conversion
 */
 
namespace MayaFlux::Yantra {
 
/**
 * @enum StatisticalMethod
 * @brief Supported statistical computation methods
 */
enum class StatisticalMethod : uint8_t {
    MEAN, ///< Arithmetic mean
    VARIANCE, ///< Population or sample variance
    STD_DEV, ///< Standard deviation
    SKEWNESS, ///< Third moment - asymmetry measure
    KURTOSIS, ///< Fourth moment - tail heaviness
    MIN, ///< Minimum value
    MAX, ///< Maximum value
    MEDIAN, ///< 50th percentile
    RANGE, ///< Max - min
    PERCENTILE, ///< Arbitrary percentile (requires parameter)
    MODE, ///< Most frequent value
    MAD, ///< Median Absolute Deviation
    CV, ///< Coefficient of Variation (std_dev/mean)
    SUM, ///< Sum of all values
    COUNT, ///< Number of values
    RMS, ///< Root Mean Square
    ENTROPY, ///< Shannon entropy for discrete data
    ZSCORE ///< Z-score normalization
};
 
/**
 * @enum StatisticalLevel
 * @brief Qualitative classification of statistical values
 */
enum class StatisticalLevel : uint8_t {
    EXTREME_LOW,
    LOW,
    NORMAL,
    HIGH,
    EXTREME_HIGH,
    OUTLIER
};
 
/**
 * @struct ChannelStatistics
 * @brief Statistical results for a single data channel
 */
struct MAYAFLUX_API ChannelStatistics {
    std::vector<double> statistical_values;
 
    double mean_stat {};
    double max_stat {};
    double min_stat {};
    double stat_variance {};
    double stat_std_dev {};
 
    double skewness {};
    double kurtosis {};
    double median {};
    std::vector<double> percentiles;
 
    std::vector<StatisticalLevel> stat_classifications;
    std::array<int, 6> level_counts {};
 
    std::vector<std::pair<size_t, size_t>> window_positions;
    std::map<std::string, std::any> method_specific_data;
};
 
/**
 * @struct StatisticalAnalysis
 * @brief Analysis result structure for statistical analysis
 */
struct MAYAFLUX_API StatisticalAnalysis {
    StatisticalMethod method_used { StatisticalMethod::MEAN };
    uint32_t window_size {};
    uint32_t hop_size {};
 
    std::vector<ChannelStatistics> channel_statistics;
};
 
/**
 * @class StatisticalAnalyzer
 * @brief High-performance statistical analyzer with zero-copy processing
 *
 * The StatisticalAnalyzer provides comprehensive statistical analysis capabilities for
 * digital data streams using span-based processing for maximum efficiency.
 * All data extraction and conversion is handled automatically by OperationHelper.
 *
 * Example usage:
 * ```cpp
 * // DataVariant -> VectorXd analyzer
 * auto stat_analyzer = std::make_shared<StatisticalAnalyzer<Kakshya::DataVariant, Eigen::VectorXd>>();
 *
 * // User-facing analysis
 * auto analysis = stat_analyzer->analyze_data(numeric_data);
 * auto stat_result = safe_any_cast<StatisticalAnalysisResult>(analysis);
 *
 * // Pipeline usage
 * auto pipeline_output = stat_analyzer->apply_operation(IO{numeric_data});
 * ```
 */
template <ComputeData InputType = std::vector<Kakshya::DataVariant>, ComputeData OutputType = Eigen::VectorXd>
class MAYAFLUX_API StatisticalAnalyzer : public UniversalAnalyzer<InputType, OutputType> {
public:
    using input_type = Datum<InputType>;
    using output_type = Datum<OutputType>;
    using base_type = UniversalAnalyzer<InputType, OutputType>;
 
    /**
     * @brief Construct StatisticalAnalyzer with configurable window parameters
     * @param window_size Size of analysis window in samples (default: 512)
     * @param hop_size Step size between windows in samples (default: 256)
     */
    explicit StatisticalAnalyzer(uint32_t window_size = 512, uint32_t hop_size = 256)
        : m_window_size(window_size)
        , m_hop_size(hop_size)
    {
        validate_window_parameters();
    }
 
    /**
     * @brief Type-safe statistical analysis method
     * @param data Input data
     * @return StatisticalAnalysisResult directly
     */
    StatisticalAnalysis analyze_statistics(const InputType& data)
    {
        auto result = this->analyze_data(data);
        return safe_any_cast_or_throw<StatisticalAnalysis>(result);
    }
 
    /**
     * @brief Get last statistical analysis result (type-safe)
     * @return StatisticalAnalysisResult from last operation
     */
    [[nodiscard]] StatisticalAnalysis get_statistical_analysis() const
    {
        return safe_any_cast_or_throw<StatisticalAnalysis>(this->get_current_analysis());
    }
 
    /**
     * @brief Get analysis type category
     * @return AnalysisType::STATISTICAL
     */
    [[nodiscard]] AnalysisType get_analysis_type() const override
    {
        return AnalysisType::STATISTICAL;
    }
 
    /**
     * @brief Get available analysis methods
     * @return Vector of supported statistical method names
     */
    [[nodiscard]] std::vector<std::string> get_available_methods() const override
    {
        return Reflect::get_enum_names_lowercase<StatisticalMethod>();
    }
 
    /**
     * @brief Get supported methods for specific input type
     * @tparam T Input type to check
     * @return Vector of method names supported for this type
     */
    template <typename T>
    [[nodiscard]] std::vector<std::string> get_methods_for_type() const
    {
        return get_methods_for_type_impl(std::type_index(typeid(T)));
    }
 
    /**
     * @brief Check if analyzer supports given input type
     * @tparam T Input type to check
     * @return True if supported
     */
    template <typename T>
    [[nodiscard]] bool supports_input_type() const
    {
        return !get_methods_for_type<T>().empty();
    }
 
    /**
     * @brief Set statistical analysis method
     * @param method StatisticalMethod enum value
     */
    void set_method(StatisticalMethod method)
    {
        m_method = method;
        this->set_parameter("method", method_to_string(method));
    }
 
    /**
     * @brief Set method by string name
     * @param method_name String representation of method
     */
    void set_method(const std::string& method_name)
    {
        m_method = string_to_method(method_name);
        this->set_parameter("method", method_name);
    }
 
    /**
     * @brief Get current statistical method
     * @return StatisticalMethod enum value
     */
    [[nodiscard]] StatisticalMethod get_method() const
    {
        return m_method;
    }
 
    /**
     * @brief Set window size for windowed analysis
     * @param size Window size in samples
     */
    void set_window_size(uint32_t size)
    {
        m_window_size = size;
        validate_window_parameters();
    }
 
    /**
     * @brief Set hop size for windowed analysis
     * @param size Hop size in samples
     */
    void set_hop_size(uint32_t size)
    {
        m_hop_size = size;
        validate_window_parameters();
    }
 
    /**
     * @brief Get window size
     * @return Current window size
     */
    [[nodiscard]] uint32_t get_window_size() const { return m_window_size; }
 
    /**
     * @brief Get hop size
     * @return Current hop size
     */
    [[nodiscard]] uint32_t get_hop_size() const { return m_hop_size; }
 
    /**
     * @brief Enable/disable outlier classification
     * @param enabled Whether to classify outliers
     */
    void set_classification_enabled(bool enabled)
    {
        m_classification_enabled = enabled;
    }
 
    /**
     * @brief Check if classification is enabled
     * @return True if classification enabled
     */
    [[nodiscard]] bool is_classification_enabled() const { return m_classification_enabled; }
 
    /**
     * @brief Classify statistical value qualitatively
     * @param value Statistical value to classify
     * @return StatisticalLevel classification
     */
    [[nodiscard]] StatisticalLevel classify_statistical_level(double value) const
    {
        if (std::abs(value) > m_outlier_threshold)
            return StatisticalLevel::OUTLIER;
        if (value <= m_extreme_low_threshold)
            return StatisticalLevel::EXTREME_LOW;
        if (value <= m_low_threshold)
            return StatisticalLevel::LOW;
        if (value <= m_high_threshold)
            return StatisticalLevel::NORMAL;
        if (value <= m_extreme_high_threshold)
            return StatisticalLevel::HIGH;
        return StatisticalLevel::EXTREME_HIGH;
    }
 
    /**
     * @brief Convert statistical method enum to string
     * @param method StatisticalMethod value
     * @return String representation
     */
    static std::string method_to_string(StatisticalMethod method)
    {
        return Reflect::enum_to_lowercase_string(method);
    }
 
    /**
     * @brief Convert string to statistical method enum
     * @param str String representation
     * @return StatisticalMethod value
     */
    static StatisticalMethod string_to_method(const std::string& str)
    {
        if (str == "default")
            return StatisticalMethod::MEAN;
        return Reflect::string_to_enum_or_throw_case_insensitive<StatisticalMethod>(str, "StatisticalMethod");
    }
 
    /**
     * @brief Convert statistical level enum to string
     * @param level StatisticalLevel value
     * @return String representation
     */
    static std::string statistical_level_to_string(StatisticalLevel level)
    {
        return Reflect::enum_to_lowercase_string(level);
    }
 
protected:
    /**
     * @brief Get analyzer name
     * @return "StatisticalAnalyzer"
     */
    [[nodiscard]] std::string get_analyzer_name() const override
    {
        return "StatisticalAnalyzer";
    }
 
    /**
     * @brief Core analysis implementation - creates analysis result AND pipeline output
     * @param input Input data wrapped in Datum container
     * @return Pipeline output (data flow for chaining operations)
     */
    output_type analyze_implementation(const input_type& input) override
    {
        if (input.data.empty()) {
            error<std::runtime_error>(Journal::Component::Yantra, Journal::Context::ComputeMatrix, std::source_location::current(), "Input is empty");
        }
        try {
            auto [data_span, structure_info] = OperationHelper::extract_structured_double(
                const_cast<input_type&>(input));
 
            std::vector<std::span<const double>> channel_spans;
            for (const auto& span : data_span)
                channel_spans.emplace_back(span.data(), span.size());
 
            std::vector<std::vector<double>> stat_values;
            stat_values.reserve(channel_spans.size());
            for (const auto& ch_span : channel_spans) {
                stat_values.push_back(compute_statistical_values(ch_span, m_method));
            }
 
            StatisticalAnalysis analysis_result = create_analysis_result(
                stat_values, channel_spans, structure_info);
 
            this->store_current_analysis(analysis_result);
            return create_pipeline_output(input, analysis_result, structure_info);
        } catch (const std::exception& e) {
            MF_ERROR(Journal::Component::Yantra, Journal::Context::ComputeMatrix, "Statistical analysis failed: {}", e.what());
            output_type error_result;
            error_result.metadata = input.metadata;
            error_result.metadata["error"] = std::string("Analysis failed: ") + e.what();
            return error_result;
        }
    }
 
    /**
     * @brief Handle analysis-specific parameters
     */
    void set_analysis_parameter(const std::string& name, std::any value) override
    {
        try {
            if (name == "method") {
                try {
                    auto method_str = safe_any_cast_or_throw<std::string>(value);
                    m_method = string_to_method(method_str);
                } catch (const std::runtime_error&) {
                    auto method_enum = safe_any_cast_or_throw<StatisticalMethod>(value);
                    m_method = method_enum;
                }
            } else if (name == "window_size") {
                auto size = safe_any_cast_or_throw<uint32_t>(value);
                m_window_size = size;
                validate_window_parameters();
            } else if (name == "hop_size") {
                auto size = safe_any_cast_or_throw<uint32_t>(value);
                m_hop_size = size;
                validate_window_parameters();
            } else if (name == "classification_enabled") {
                auto enabled = safe_any_cast_or_throw<bool>(value);
                m_classification_enabled = enabled;
            } else if (name == "percentile") {
                auto percentile = safe_any_cast_or_throw<double>(value);
                if (percentile < 0.0 || percentile > 100.0) {
                    throw std::invalid_argument("Percentile must be between 0.0 and 100.0, got: " + std::to_string(percentile));
                }
                m_percentile_value = percentile;
            } else if (name == "sample_variance") {
                auto sample = safe_any_cast_or_throw<bool>(value);
                m_sample_variance = sample;
            } else {
                base_type::set_analysis_parameter(name, std::move(value));
            }
        } catch (const std::runtime_error& e) {
            error_rethrow(Journal::Component::Yantra, Journal::Context::ComputeMatrix, std::source_location::current(), "Failed to set parameter '{}': {}", name, e.what());
        }
    }
 
    /**
     * @brief Get analysis-specific parameter
     */
    [[nodiscard]] std::any get_analysis_parameter(const std::string& name) const override
    {
        if (name == "method")
            return std::any(method_to_string(m_method));
        if (name == "window_size")
            return std::any(m_window_size);
        if (name == "hop_size")
            return std::any(m_hop_size);
        if (name == "classification_enabled")
            return std::any(m_classification_enabled);
        if (name == "percentile")
            return std::any(m_percentile_value);
        if (name == "sample_variance")
            return std::any(m_sample_variance);
 
        return base_type::get_analysis_parameter(name);
    }
 
    /**
     * @brief Get supported methods for specific type index
     * @param type_info Type index to check
     * @return Vector of supported method names
     */
    [[nodiscard]] std::vector<std::string> get_methods_for_type_impl(std::type_index /*type_info*/) const
    {
        return get_available_methods();
    }
 
private:
    StatisticalMethod m_method { StatisticalMethod::MEAN };
    uint32_t m_window_size;
    uint32_t m_hop_size;
    bool m_classification_enabled { true };
    double m_percentile_value { 50.0 };
    bool m_sample_variance { true };
 
    double m_outlier_threshold { 3.0 };
    double m_extreme_low_threshold { -2.0 };
    double m_low_threshold { -1.0 };
    double m_high_threshold { 1.0 };
    double m_extreme_high_threshold { 2.0 };
 
    /**
     * @brief Validate window parameters
     */
    void validate_window_parameters() const
    {
        if (m_window_size == 0 || m_hop_size == 0) {
            error<std::invalid_argument>(Journal::Component::Yantra, Journal::Context::ComputeMatrix, std::source_location::current(), "Window size and hop size must be greater than 0");
        }
        if (m_hop_size > m_window_size) {
            error<std::invalid_argument>(Journal::Component::Yantra, Journal::Context::ComputeMatrix, std::source_location::current(), "Hop size should not exceed window size");
        }
    }
 
    /**
     * @brief Compute statistical values using span (zero-copy processing)
     */
    [[nodiscard]] std::vector<double> compute_statistical_values(std::span<const double> data, StatisticalMethod method) const
    {
        const size_t num_windows = calculate_num_windows(data.size());
 
        namespace D = MayaFlux::Kinesis::Discrete;
 
        switch (method) {
        case StatisticalMethod::MEAN:
            return D::mean(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::VARIANCE:
            return D::variance(data, num_windows, m_hop_size, m_window_size, m_sample_variance);
        case StatisticalMethod::STD_DEV:
            return D::std_dev(data, num_windows, m_hop_size, m_window_size, m_sample_variance);
        case StatisticalMethod::SKEWNESS:
            return D::skewness(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::KURTOSIS:
            return D::kurtosis(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::MEDIAN:
            return D::median(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::PERCENTILE:
            return D::percentile(data, num_windows, m_hop_size, m_window_size, m_percentile_value);
        case StatisticalMethod::ENTROPY:
            return D::entropy(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::MIN:
            return D::min(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::MAX:
            return D::max(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::RANGE:
            return D::range(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::SUM:
            return D::sum(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::COUNT:
            return D::count(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::RMS:
            return D::rms(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::MAD:
            return D::mad(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::CV:
            return D::coefficient_of_variation(data, num_windows, m_hop_size, m_window_size, m_sample_variance);
        case StatisticalMethod::MODE:
            return D::mode(data, num_windows, m_hop_size, m_window_size);
        case StatisticalMethod::ZSCORE:
            return D::mean_zscore(data, num_windows, m_hop_size, m_window_size, m_sample_variance);
        default:
            return D::mean(data, num_windows, m_hop_size, m_window_size);
        }
    }
 
    /**
     * @brief Calculate number of windows for given data size
     */
    [[nodiscard]] size_t calculate_num_windows(size_t data_size) const
    {
        if (data_size < m_window_size)
            return 0;
        return (data_size - m_window_size) / m_hop_size + 1;
    }
 
    /**
     * @brief Create comprehensive analysis result
     */
    StatisticalAnalysis create_analysis_result(const std::vector<std::vector<double>>& stat_values,
        std::vector<std::span<const double>> original_data, const auto& /*structure_info*/) const
    {
        namespace D = MayaFlux::Kinesis::Discrete;
 
        StatisticalAnalysis result;
        result.method_used = m_method;
        result.window_size = m_window_size;
        result.hop_size = m_hop_size;
 
        if (stat_values.empty()) {
            return result;
        }
 
        result.channel_statistics.resize(stat_values.size());
 
        for (size_t ch = 0; ch < stat_values.size(); ++ch) {
            auto& channel_result = result.channel_statistics[ch];
            const auto& ch_stats = stat_values[ch];
 
            channel_result.statistical_values = ch_stats;
 
            if (ch_stats.empty())
                continue;
 
            const auto [min_it, max_it] = std::ranges::minmax_element(ch_stats);
            channel_result.min_stat = *min_it;
            channel_result.max_stat = *max_it;
 
            const std::span<const double> sp(ch_stats);
            const auto sz = static_cast<uint32_t>(ch_stats.size());
 
            const auto single = [&](auto fn) { return fn(sp, 1, 0, sz)[0]; };
 
            channel_result.mean_stat = single([](auto&&... a) { return D::mean(a...); });
            channel_result.stat_variance = single([&](auto&&... a) { return D::variance(a..., m_sample_variance); });
            channel_result.stat_std_dev = std::sqrt(channel_result.stat_variance);
            channel_result.skewness = single([](auto&&... a) { return D::skewness(a...); });
            channel_result.kurtosis = single([](auto&&... a) { return D::kurtosis(a...); });
            channel_result.median = single([](auto&&... a) { return D::median(a...); });
 
            channel_result.percentiles = {
                D::percentile(sp, 1, 0, sz, 25.0)[0],
                channel_result.median,
                D::percentile(sp, 1, 0, sz, 75.0)[0]
            };
 
            const size_t data_size = (ch < original_data.size()) ? original_data[ch].size() : 0;
            channel_result.window_positions.reserve(ch_stats.size());
            for (size_t i = 0; i < ch_stats.size(); ++i) {
                const size_t start = i * m_hop_size;
                const size_t end = std::min(start + m_window_size, data_size);
                channel_result.window_positions.emplace_back(start, end);
            }
 
            if (m_classification_enabled) {
                channel_result.stat_classifications.reserve(ch_stats.size());
                channel_result.level_counts.fill(0);
 
                for (double value : ch_stats) {
                    const StatisticalLevel level = classify_statistical_level(value);
                    channel_result.stat_classifications.push_back(level);
                    channel_result.level_counts[static_cast<size_t>(level)]++;
                }
            }
        }
 
        return result;
    }
 
    /**
     * @brief Create pipeline output for operation chaining
     */
    output_type create_pipeline_output(const input_type& input, const StatisticalAnalysis& analysis_result, DataStructureInfo& info)
    {
        std::vector<std::vector<double>> channel_stats;
        channel_stats.reserve(analysis_result.channel_statistics.size());
        for (const auto& ch : analysis_result.channel_statistics) {
            channel_stats.push_back(ch.statistical_values);
        }
 
        output_type output = this->convert_result(channel_stats, info);
 
        output.metadata = input.metadata;
 
        output.metadata["source_analyzer"] = "StatisticalAnalyzer";
        output.metadata["statistical_method"] = method_to_string(analysis_result.method_used);
        output.metadata["window_size"] = analysis_result.window_size;
        output.metadata["hop_size"] = analysis_result.hop_size;
        output.metadata["num_channels"] = analysis_result.channel_statistics.size();
 
        if (!analysis_result.channel_statistics.empty()) {
            std::vector<double> channel_means, channel_maxs, channel_mins, channel_variances, channel_stddevs, channel_skewness, channel_kurtosis, channel_medians;
            std::vector<size_t> channel_window_counts;
 
            for (const auto& ch : analysis_result.channel_statistics) {
                channel_means.push_back(ch.mean_stat);
                channel_maxs.push_back(ch.max_stat);
                channel_mins.push_back(ch.min_stat);
                channel_variances.push_back(ch.stat_variance);
                channel_stddevs.push_back(ch.stat_std_dev);
                channel_skewness.push_back(ch.skewness);
                channel_kurtosis.push_back(ch.kurtosis);
                channel_medians.push_back(ch.median);
                channel_window_counts.push_back(ch.statistical_values.size());
            }
 
            output.metadata["mean_per_channel"] = channel_means;
            output.metadata["max_per_channel"] = channel_maxs;
            output.metadata["min_per_channel"] = channel_mins;
            output.metadata["variance_per_channel"] = channel_variances;
            output.metadata["stddev_per_channel"] = channel_stddevs;
            output.metadata["skewness_per_channel"] = channel_skewness;
            output.metadata["kurtosis_per_channel"] = channel_kurtosis;
            output.metadata["median_per_channel"] = channel_medians;
            output.metadata["window_count_per_channel"] = channel_window_counts;
        }
 
        return output;
    }
};
 
/// Standard statistical analyzer: DataVariant -> MatrixXd
using StandardStatisticalAnalyzer = StatisticalAnalyzer<std::vector<Kakshya::DataVariant>, Eigen::MatrixXd>;
 
/// Container statistical analyzer: SignalContainer -> MatrixXd
using ContainerStatisticalAnalyzer = StatisticalAnalyzer<std::shared_ptr<Kakshya::SignalSourceContainer>, Eigen::MatrixXd>;
 
/// Region statistical analyzer: Region -> MatrixXd
using RegionStatisticalAnalyzer = StatisticalAnalyzer<Kakshya::Region, Eigen::VectorXd>;
 
/// Raw statistical analyzer: produces double vectors
template <ComputeData InputType = std::vector<Kakshya::DataVariant>>
using RawStatisticalAnalyzer = StatisticalAnalyzer<InputType, std::vector<std::vector<double>>>;
 
/// Variant statistical analyzer: produces DataVariant output
template <ComputeData InputType = std::vector<Kakshya::DataVariant>>
using VariantStatisticalAnalyzer = StatisticalAnalyzer<InputType, Kakshya::DataVariant>;
 
} // namespace MayaFlux::Yantra