#!/usr/bin/env ruby

require_relative 'shared_constants'

# Statistical calculations module
module StatisticsCalculator
  extend self

  # Calculate percentile for a sorted array
  def percentile(sorted_array, percentile)
    return nil if sorted_array.empty?
    
    index = (percentile / 100.0) * (sorted_array.length - 1)
    lower_index = index.floor
    upper_index = index.ceil
    
    if lower_index == upper_index
      sorted_array[lower_index]
    else
      lower_value = sorted_array[lower_index]
      upper_value = sorted_array[upper_index]
      weight = index - lower_index
      lower_value + (upper_value - lower_value) * weight
    end
  end

  # Calculate basic statistical measures
  def basic_stats(values)
    return {} if values.empty?
    
    sorted_values = values.sort
    length = values.length
    
    sum = values.sum
    mean = sum.to_f / length
    
    # Calculate variance and standard deviation
    variance = values.map { |x| (x - mean) ** 2 }.sum / (length - 1)
    variance = 0 if variance.negative? || variance.nan?
    std_dev = Math.sqrt(variance)
    
    {
      count: length,
      sum: sum,
      mean: mean,
      median: percentile(sorted_values, 50),
      min: sorted_values.first,
      max: sorted_values.last,
      variance: variance,
      std_dev: std_dev,
      std_error: std_dev / Math.sqrt(length),
      q1: percentile(sorted_values, 25),
      q3: percentile(sorted_values, 75),
      p90: percentile(sorted_values, 90),
      p95: percentile(sorted_values, 95)
    }
  end

  # Detect outliers using IQR method
  def detect_outliers_iqr(values)
    return { outliers: [], clean_data: values, bounds: [nil, nil] } if values.length < 4
    
    sorted_values = values.sort
    q1 = percentile(sorted_values, 25)
    q3 = percentile(sorted_values, 75)
    iqr = q3 - q1
    
    multiplier = SharedConstants::Statistics::OUTLIER_IQR_MULTIPLIER
    lower_bound = q1 - multiplier * iqr
    upper_bound = q3 + multiplier * iqr
    
    outliers = values.select { |value| value < lower_bound || value > upper_bound }
    clean_data = values.select { |value| value >= lower_bound && value <= upper_bound }
    
    {
      outliers: outliers,
      clean_data: clean_data,
      bounds: [lower_bound, upper_bound],
      q1: q1,
      q3: q3,
      iqr: iqr
    }
  end

  # Calculate confidence interval
  def confidence_interval(values, confidence_level = SharedConstants::Defaults::CONFIDENCE_LEVEL)
    return [nil, nil] if values.length < 2
    
    stats = basic_stats(values)
    critical_value = get_critical_value(confidence_level, values.length)
    
    return [nil, nil] if critical_value.nil?
    
    margin_of_error = critical_value * stats[:std_error]
    
    [stats[:mean] - margin_of_error, stats[:mean] + margin_of_error]
  end

  # Calculate weighted mean
  def weighted_mean(values, weights = nil)
    return nil if values.empty?
    return values.sum / values.length if weights.nil? || weights.empty?
    
    weighted_sum = values.zip(weights).map { |value, weight| value * weight }.sum
    total_weight = weights.sum
    
    return nil if total_weight.zero?
    
    weighted_sum / total_weight
  end

  # Generate histogram data for visualization
  def generate_histogram(values, bins = SharedConstants::Defaults::HISTOGRAM_BINS)
    return { bins: [], counts: [], is_skewed: false } if values.empty?
    
    sorted_values = values.sort
    min_val = sorted_values.first
    max_val = sorted_values.last
    range = max_val - min_val
    
    # Handle case where all values are the same
    if range.zero?
      return {
        bins: [[min_val, max_val]],
        counts: [values.length],
        is_skewed: false,
        type: :uniform
      }
    end
    
    # Check for skewness
    skew_ratio = calculate_skewness(sorted_values)
    is_highly_skewed = skew_ratio > SharedConstants::Statistics::SKEW_THRESHOLD
    
    if is_highly_skewed && min_val > 0
      generate_log_histogram(sorted_values, min_val, max_val)
    else
      generate_linear_histogram(sorted_values, min_val, max_val, bins)
    end
  end

  # Calculate dose classification thresholds
  def dose_classification_thresholds(values)
    return {} if values.empty?
    
    stats = basic_stats(values)
    
    # Use percentile-based classification
    {
      very_light: { max: stats[:q1] * 0.5, color: :green },
      light: { max: stats[:q1], color: :cyan },
      common: { max: stats[:median], color: :yellow },
      strong: { max: stats[:q3], color: :magenta },
      heavy: { max: stats[:p90], color: :red },
      extreme: { max: Float::INFINITY, color: :red }
    }
  end

  # Calculate risk assessment thresholds
  def risk_thresholds(values)
    return {} if values.empty?
    
    stats = basic_stats(values)
    
    {
      low_risk: stats[:mean] - stats[:std_dev],
      moderate_risk: stats[:mean] + stats[:std_dev],
      high_risk: stats[:mean] + 2 * stats[:std_dev],
      extreme_risk: stats[:mean] + 3 * stats[:std_dev]
    }
  end

  # Calculate coefficient of variation
  def coefficient_of_variation(values)
    return nil if values.empty?
    
    stats = basic_stats(values)
    return nil if stats[:mean].zero?
    
    (stats[:std_dev] / stats[:mean]) * 100
  end

  # Classify data distribution type
  def classify_distribution(values)
    return :unknown if values.empty?
    
    stats = basic_stats(values)
    cv = coefficient_of_variation(values)
    
    case cv
    when 0..10
      :very_consistent
    when 10..25
      :consistent
    when 25..50
      :moderate_variation
    when 50..100
      :high_variation
    else
      :extreme_variation
    end
  end

  # Calculate z-score for a value
  def z_score(value, mean, std_dev)
    return nil if std_dev.zero?
    (value - mean) / std_dev
  end

  # Calculate percentile rank of a value
  def percentile_rank(value, values)
    return nil if values.empty?
    
    count_below = values.count { |v| v < value }
    count_equal = values.count { |v| v == value }
    
    ((count_below + count_equal * 0.5) / values.length) * 100
  end

  private

  # Get critical value for confidence interval
  def get_critical_value(confidence_level, sample_size)
    critical_values = SharedConstants::Statistics::CRITICAL_VALUES
    
    # Use t-distribution approximation for small samples
    if sample_size <= SharedConstants::Statistics::LARGE_SAMPLE_THRESHOLD
      case confidence_level
      when 0.95
        2.576  # Conservative approximation
      when 0.99
        3.291  # Conservative approximation
      else
        critical_values[confidence_level]
      end
    else
      critical_values[confidence_level]
    end
  end

  # Calculate skewness ratio
  def calculate_skewness(sorted_values)
    return 0 if sorted_values.length < 4
    
    min_val = sorted_values.first
    max_val = sorted_values.last
    range = max_val - min_val
    
    return 0 if range.zero?
    
    # Calculate how many values are in the first quarter of the range
    first_quarter_threshold = min_val + (range * 0.25)
    values_in_first_quarter = sorted_values.count { |v| v <= first_quarter_threshold }
    
    values_in_first_quarter.to_f / sorted_values.length
  end

  # Generate logarithmic histogram for skewed data
  def generate_log_histogram(sorted_values, min_val, max_val)
    log_min = Math.log10(min_val)
    log_max = Math.log10(max_val)
    log_range = log_max - log_min
    
    bins_count = SharedConstants::Statistics::LOG_BINS
    bin_edges = (0..bins_count).map { |i| 10**(log_min + i * log_range / bins_count) }
    
    # Ensure we capture all data points
    bin_edges[0] = [bin_edges[0], min_val].min
    bin_edges[-1] = [bin_edges[-1], max_val].max
    
    counts = calculate_bin_counts(sorted_values, bin_edges)
    bins = bin_edges.each_cons(2).to_a
    
    {
      bins: bins,
      counts: counts,
      is_skewed: true,
      type: :logarithmic
    }
  end

  # Generate linear histogram for normal data
  def generate_linear_histogram(sorted_values, min_val, max_val, bins_count)
    range = max_val - min_val
    bin_width = range / bins_count.to_f
    
    bin_edges = (0..bins_count).map { |i| min_val + i * bin_width }
    counts = calculate_bin_counts(sorted_values, bin_edges)
    bins = bin_edges.each_cons(2).to_a
    
    {
      bins: bins,
      counts: counts,
      is_skewed: false,
      type: :linear
    }
  end

  # Calculate counts for histogram bins
  def calculate_bin_counts(sorted_values, bin_edges)
    counts = Array.new(bin_edges.length - 1, 0)
    
    sorted_values.each do |value|
      bin_index = bin_edges.length - 2  # default to last bin
      
      (0...bin_edges.length - 1).each do |i|
        if value >= bin_edges[i] && value < bin_edges[i + 1]
          bin_index = i
          break
        elsif value == bin_edges[-1]  # handle max value
          bin_index = bin_edges.length - 2
          break
        end
      end
      
      counts[bin_index] += 1
    end
    
    counts
  end
end
