viva_math

-module(viva_math@statistics).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/statistics.gleam").
-export([sum/1, mean/1, geometric_mean/1, harmonic_mean/1, weighted_mean/2, welford_empty/0, welford_update/2, welford/1, variance/1, sample_variance/1, stddev/1, sample_stddev/1, covariance/2, pearson/2, min/1, max/1, range/1, percentile/2, median/1, quartiles/1, iqr/1, z_score/1, min_max_normalize/1, moving_average/2, ema/2, ema_step/3, skewness/1, kurtosis/1]).
-export_type([welford/0]).

-if(?OTP_RELEASE >= 27).
-define(MODULEDOC(Str), -moduledoc(Str)).
-define(DOC(Str), -doc(Str)).
-else.
-define(MODULEDOC(Str), -compile([])).
-define(DOC(Str), -compile([])).
-endif.

?MODULEDOC(
    " Descriptive statistics over `List(Float)`.\n"
    "\n"
    " All functions are pure and operate on lists. For large datasets prefer\n"
    " `viva_tensor` and stream reductions; this module targets emotion-scale\n"
    " data (10²–10⁴ samples).\n"
    "\n"
    " ## Numerical notes\n"
    "\n"
    " - Variance uses Welford's online algorithm (`welford`) to avoid\n"
    "   catastrophic cancellation. `variance/1` materialises the result.\n"
    " - Population (N) and sample (N-1) variants are provided where the\n"
    "   distinction matters.\n"
    " - Percentile uses linear interpolation between adjacent order statistics\n"
    "   (NumPy-style \"linear\" method).\n"
).

-type welford() :: {welford, integer(), float(), float()}.

-file("src/viva_math/statistics.gleam", 29).
?DOC(
    " Sum of a list using Neumaier compensated summation. Returns 0.0 for empty.\n"
    "\n"
    " Recovers up to 16 extra bits of precision vs naive `fold (+.)`.\n"
).
-spec sum(list(float())) -> float().
sum(Xs) ->
    viva_math@precision:neumaier_sum(Xs).

-file("src/viva_math/statistics.gleam", 34).
?DOC(" Arithmetic mean using Neumaier-compensated sum.\n").
-spec mean(list(float())) -> {ok, float()} | {error, nil}.
mean(Xs) ->
    viva_math@precision:neumaier_mean(Xs).

-file("src/viva_math/statistics.gleam", 424).
-spec int_to_float(integer()) -> float().
int_to_float(N) ->
    erlang:float(N).

-file("src/viva_math/statistics.gleam", 39).
?DOC(" Geometric mean: ⁿ√(∏ xᵢ). Requires strictly positive inputs.\n").
-spec geometric_mean(list(float())) -> {ok, float()} | {error, nil}.
geometric_mean(Xs) ->
    case Xs of
        [] ->
            {error, nil};

        _ ->
            N = int_to_float(erlang:length(Xs)),
            All_positive = gleam@list:all(Xs, fun(X) -> X > +0.0 end),
            case All_positive of
                false ->
                    {error, nil};

                true ->
                    Log_sum = gleam@list:fold(
                        Xs,
                        +0.0,
                        fun(Acc, X@1) -> Acc + math:log(X@1) end
                    ),
                    {ok, math:exp(case N of
                                +0.0 -> +0.0;
                                -0.0 -> -0.0;
                                Gleam@denominator -> Log_sum / Gleam@denominator
                            end)}
            end
    end.

-file("src/viva_math/statistics.gleam", 57).
?DOC(" Harmonic mean: n / Σ(1/xᵢ). Requires non-zero inputs.\n").
-spec harmonic_mean(list(float())) -> {ok, float()} | {error, nil}.
harmonic_mean(Xs) ->
    case Xs of
        [] ->
            {error, nil};

        _ ->
            N = int_to_float(erlang:length(Xs)),
            Inv_sum = gleam@list:fold(
                Xs,
                {ok, +0.0},
                fun(Acc, X) -> case {Acc, X =:= +0.0} of
                        {{error, _}, _} ->
                            {error, nil};

                        {_, true} ->
                            {error, nil};

                        {{ok, S}, false} ->
                            {ok, S + (case X of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator -> 1.0 / Gleam@denominator
                                end)}
                    end end
            ),
            case Inv_sum of
                {ok, S@1} ->
                    case S@1 =:= +0.0 of
                        true ->
                            {error, nil};

                        false ->
                            {ok, case S@1 of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator@1 -> N / Gleam@denominator@1
                                end}
                    end;

                {error, _} ->
                    {error, nil}
            end
    end.

-file("src/viva_math/statistics.gleam", 83).
?DOC(" Weighted mean: Σ(wᵢ · xᵢ) / Σwᵢ.\n").
-spec weighted_mean(list(float()), list(float())) -> {ok, float()} |
    {error, nil}.
weighted_mean(Xs, Ws) ->
    case {erlang:length(Xs) =:= erlang:length(Ws), Xs} of
        {false, _} ->
            {error, nil};

        {_, []} ->
            {error, nil};

        {true, _} ->
            Pairs = gleam@list:zip(Xs, Ws),
            Total_w = gleam@list:fold(
                Pairs,
                +0.0,
                fun(Acc, P) -> Acc + erlang:element(2, P) end
            ),
            case Total_w =:= +0.0 of
                true ->
                    {error, nil};

                false ->
                    Num = gleam@list:fold(
                        Pairs,
                        +0.0,
                        fun(Acc@1, P@1) ->
                            Acc@1 + (erlang:element(1, P@1) * erlang:element(
                                2,
                                P@1
                            ))
                        end
                    ),
                    {ok, case Total_w of
                            +0.0 -> +0.0;
                            -0.0 -> -0.0;
                            Gleam@denominator -> Num / Gleam@denominator
                        end}
            end
    end.

-file("src/viva_math/statistics.gleam", 111).
?DOC(" Empty Welford accumulator.\n").
-spec welford_empty() -> welford().
welford_empty() ->
    {welford, 0, +0.0, +0.0}.

-file("src/viva_math/statistics.gleam", 116).
?DOC(" Update a Welford accumulator with a new sample.\n").
-spec welford_update(welford(), float()) -> welford().
welford_update(W, X) ->
    Count = erlang:element(2, W) + 1,
    Delta = X - erlang:element(3, W),
    New_mean = erlang:element(3, W) + (case int_to_float(Count) of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> Delta / Gleam@denominator
    end),
    Delta2 = X - New_mean,
    {welford, Count, New_mean, erlang:element(4, W) + (Delta * Delta2)}.

-file("src/viva_math/statistics.gleam", 125).
?DOC(" Build a Welford accumulator from a list.\n").
-spec welford(list(float())) -> welford().
welford(Xs) ->
    gleam@list:fold(Xs, welford_empty(), fun welford_update/2).

-file("src/viva_math/statistics.gleam", 130).
?DOC(" Population variance σ² = M₂ / N.\n").
-spec variance(list(float())) -> {ok, float()} | {error, nil}.
variance(Xs) ->
    W = welford(Xs),
    case erlang:element(2, W) of
        0 ->
            {error, nil};

        N ->
            {ok, case int_to_float(N) of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> erlang:element(4, W) / Gleam@denominator
                end}
    end.

-file("src/viva_math/statistics.gleam", 139).
?DOC(" Sample variance s² = M₂ / (N - 1).\n").
-spec sample_variance(list(float())) -> {ok, float()} | {error, nil}.
sample_variance(Xs) ->
    W = welford(Xs),
    case erlang:element(2, W) of
        0 ->
            {error, nil};

        1 ->
            {error, nil};

        N ->
            {ok, case int_to_float(N - 1) of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> erlang:element(4, W) / Gleam@denominator
                end}
    end.

-file("src/viva_math/statistics.gleam", 148).
?DOC(" Population standard deviation.\n").
-spec stddev(list(float())) -> {ok, float()} | {error, nil}.
stddev(Xs) ->
    case variance(Xs) of
        {ok, V} ->
            {ok, math:sqrt(V)};

        {error, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 156).
?DOC(" Sample standard deviation.\n").
-spec sample_stddev(list(float())) -> {ok, float()} | {error, nil}.
sample_stddev(Xs) ->
    case sample_variance(Xs) of
        {ok, V} ->
            {ok, math:sqrt(V)};

        {error, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 168).
?DOC(" Population covariance Cov(X, Y) = E[(X - μₓ)(Y - μᵧ)].\n").
-spec covariance(list(float()), list(float())) -> {ok, float()} | {error, nil}.
covariance(Xs, Ys) ->
    case {erlang:length(Xs) =:= erlang:length(Ys), Xs, Ys} of
        {false, _, _} ->
            {error, nil};

        {_, [], _} ->
            {error, nil};

        {_, _, []} ->
            {error, nil};

        {true, _, _} ->
            N = int_to_float(erlang:length(Xs)),
            case {mean(Xs), mean(Ys)} of
                {{ok, Mx}, {ok, My}} ->
                    S = gleam@list:fold(
                        gleam@list:zip(Xs, Ys),
                        +0.0,
                        fun(Acc, P) ->
                            Acc + ((erlang:element(1, P) - Mx) * (erlang:element(
                                2,
                                P
                            )
                            - My))
                        end
                    ),
                    {ok, case N of
                            +0.0 -> +0.0;
                            -0.0 -> -0.0;
                            Gleam@denominator -> S / Gleam@denominator
                        end};

                {_, _} ->
                    {error, nil}
            end
    end.

-file("src/viva_math/statistics.gleam", 190).
?DOC(" Pearson correlation coefficient. Result in [-1, 1].\n").
-spec pearson(list(float()), list(float())) -> {ok, float()} | {error, nil}.
pearson(Xs, Ys) ->
    case {covariance(Xs, Ys), stddev(Xs), stddev(Ys)} of
        {{ok, C}, {ok, Sx}, {ok, Sy}} ->
            Denom = Sx * Sy,
            case Denom =:= +0.0 of
                true ->
                    {error, nil};

                false ->
                    {ok, case Denom of
                            +0.0 -> +0.0;
                            -0.0 -> -0.0;
                            Gleam@denominator -> C / Gleam@denominator
                        end}
            end;

        {_, _, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 208).
?DOC(" Minimum of a list.\n").
-spec min(list(float())) -> {ok, float()} | {error, nil}.
min(Xs) ->
    case Xs of
        [] ->
            {error, nil};

        [X | Rest] ->
            {ok, gleam@list:fold(Rest, X, fun gleam@float:min/2)}
    end.

-file("src/viva_math/statistics.gleam", 216).
?DOC(" Maximum of a list.\n").
-spec max(list(float())) -> {ok, float()} | {error, nil}.
max(Xs) ->
    case Xs of
        [] ->
            {error, nil};

        [X | Rest] ->
            {ok, gleam@list:fold(Rest, X, fun gleam@float:max/2)}
    end.

-file("src/viva_math/statistics.gleam", 224).
?DOC(" Range = max - min.\n").
-spec range(list(float())) -> {ok, float()} | {error, nil}.
range(Xs) ->
    case {min(Xs), max(Xs)} of
        {{ok, Lo}, {ok, Hi}} ->
            {ok, Hi - Lo};

        {_, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 431).
-spec float_to_int(float()) -> integer().
float_to_int(X) ->
    erlang:trunc(X).

-file("src/viva_math/statistics.gleam", 413).
-spec list_at(list(float()), integer()) -> {ok, float()} | {error, nil}.
list_at(Xs, Idx) ->
    case {Xs, Idx} of
        {[], _} ->
            {error, nil};

        {[X | _], 0} ->
            {ok, X};

        {[_ | Rest], N} ->
            list_at(Rest, N - 1)
    end.

-file("src/viva_math/statistics.gleam", 439).
-spec float_ceil(float()) -> float().
float_ceil(X) ->
    math:ceil(X).

-file("src/viva_math/statistics.gleam", 435).
-spec float_floor(float()) -> float().
float_floor(X) ->
    math:floor(X).

-file("src/viva_math/statistics.gleam", 405).
-spec sort_floats(float(), float()) -> gleam@order:order().
sort_floats(A, B) ->
    case {A < B, A > B} of
        {true, _} ->
            lt;

        {_, true} ->
            gt;

        {_, _} ->
            eq
    end.

-file("src/viva_math/statistics.gleam", 237).
?DOC(" Percentile q ∈ [0, 1] using linear interpolation (NumPy default).\n").
-spec percentile(list(float()), float()) -> {ok, float()} | {error, nil}.
percentile(Xs, Q) ->
    case {Xs, (Q < +0.0) orelse (Q > 1.0)} of
        {[], _} ->
            {error, nil};

        {_, true} ->
            {error, nil};

        {_, false} ->
            Sorted = gleam@list:sort(Xs, fun sort_floats/2),
            N = erlang:length(Sorted),
            Pos = Q * int_to_float(N - 1),
            Lo = float_floor(Pos),
            Hi = float_ceil(Pos),
            Frac = Pos - Lo,
            case {list_at(Sorted, float_to_int(Lo)),
                list_at(Sorted, float_to_int(Hi))} of
                {{ok, A}, {ok, B}} ->
                    {ok, A + (Frac * (B - A))};

                {_, _} ->
                    {error, nil}
            end
    end.

-file("src/viva_math/statistics.gleam", 232).
?DOC(" Median (50th percentile). Linearly interpolates for even-length lists.\n").
-spec median(list(float())) -> {ok, float()} | {error, nil}.
median(Xs) ->
    percentile(Xs, 0.5).

-file("src/viva_math/statistics.gleam", 260).
?DOC(" Quartiles Q1, Q2 (median), Q3 as a tuple.\n").
-spec quartiles(list(float())) -> {ok, {float(), float(), float()}} |
    {error, nil}.
quartiles(Xs) ->
    case {percentile(Xs, 0.25), percentile(Xs, 0.5), percentile(Xs, 0.75)} of
        {{ok, Q1}, {ok, Q2}, {ok, Q3}} ->
            {ok, {Q1, Q2, Q3}};

        {_, _, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 268).
?DOC(" Interquartile range Q3 - Q1.\n").
-spec iqr(list(float())) -> {ok, float()} | {error, nil}.
iqr(Xs) ->
    case {percentile(Xs, 0.25), percentile(Xs, 0.75)} of
        {{ok, Q1}, {ok, Q3}} ->
            {ok, Q3 - Q1};

        {_, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 280).
?DOC(" Z-score: (x - μ) / σ. Returns an error if the input has zero variance.\n").
-spec z_score(list(float())) -> {ok, list(float())} | {error, nil}.
z_score(Xs) ->
    case {mean(Xs), stddev(Xs)} of
        {{ok, Mu}, {ok, Sigma}} ->
            case Sigma =:= +0.0 of
                true ->
                    {error, nil};

                false ->
                    {ok, gleam@list:map(Xs, fun(X) -> case Sigma of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator -> (X - Mu) / Gleam@denominator
                                end end)}
            end;

        {_, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 292).
?DOC(" Min-max normalisation to [0, 1].\n").
-spec min_max_normalize(list(float())) -> {ok, list(float())} | {error, nil}.
min_max_normalize(Xs) ->
    case {min(Xs), max(Xs)} of
        {{ok, Lo}, {ok, Hi}} ->
            Span = Hi - Lo,
            case Span =:= +0.0 of
                true ->
                    {error, nil};

                false ->
                    {ok, gleam@list:map(Xs, fun(X) -> case Span of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator -> (X - Lo) / Gleam@denominator
                                end end)}
            end;

        {_, _} ->
            {error, nil}
    end.

-file("src/viva_math/statistics.gleam", 335).
-spec roll_loop(list(float()), list(float()), float(), float(), list(float())) -> list(float()).
roll_loop(Leaving, Entering, Running_sum, Window_float, Acc) ->
    case {Leaving, Entering} of
        {[Out | Lrest], [In_v | Erest]} ->
            New_sum = (Running_sum - Out) + In_v,
            roll_loop(Lrest, Erest, New_sum, Window_float, [case Window_float of
                        +0.0 -> +0.0;
                        -0.0 -> -0.0;
                        Gleam@denominator -> New_sum / Gleam@denominator
                    end | Acc]);

        {_, _} ->
            lists:reverse(Acc)
    end.

-file("src/viva_math/statistics.gleam", 312).
?DOC(
    " Simple moving average over a sliding window of size `window`.\n"
    "\n"
    " Output has length `len(xs) - window + 1`. Errors if window > len or < 1.\n"
).
-spec moving_average(list(float()), integer()) -> {ok, list(float())} |
    {error, nil}.
moving_average(Xs, Window) ->
    N = erlang:length(Xs),
    case (Window < 1) orelse (Window > N) of
        true ->
            {error, nil};

        false ->
            Head = gleam@list:take(Xs, Window),
            Body = gleam@list:drop(Xs, Window),
            Initial_sum = sum(Head),
            Window_float = int_to_float(Window),
            {ok,
                roll_loop(
                    Head,
                    Body,
                    Initial_sum,
                    Window_float,
                    [case Window_float of
                            +0.0 -> +0.0;
                            -0.0 -> -0.0;
                            Gleam@denominator -> Initial_sum / Gleam@denominator
                        end]
                )}
    end.

-file("src/viva_math/statistics.gleam", 365).
-spec ema_loop(list(float()), float(), float(), list(float())) -> list(float()).
ema_loop(Xs, Alpha, Prev, Acc) ->
    case Xs of
        [] ->
            lists:reverse(Acc);

        [X | Rest] ->
            Next = (Alpha * X) + ((1.0 - Alpha) * Prev),
            ema_loop(Rest, Alpha, Next, [Next | Acc])
    end.

-file("src/viva_math/statistics.gleam", 357).
?DOC(
    " Exponential moving average with smoothing factor α ∈ (0, 1].\n"
    "\n"
    " Recurrence: yₜ = α·xₜ + (1-α)·yₜ₋₁, y₀ = x₀.\n"
).
-spec ema(list(float()), float()) -> {ok, list(float())} | {error, nil}.
ema(Xs, Alpha) ->
    case {Xs, (Alpha =< +0.0) orelse (Alpha > 1.0)} of
        {[], _} ->
            {error, nil};

        {_, true} ->
            {error, nil};

        {[X0 | Rest], false} ->
            {ok, ema_loop(Rest, Alpha, X0, [X0])}
    end.

-file("src/viva_math/statistics.gleam", 381).
?DOC(" Single-step EMA update. Useful for streaming.\n").
-spec ema_step(float(), float(), float()) -> float().
ema_step(Previous, Observation, Alpha) ->
    (Alpha * Observation) + ((1.0 - Alpha) * Previous).

-file("src/viva_math/statistics.gleam", 392).
?DOC(
    " Skewness (Fisher-Pearson) via Pébay online accumulator.\n"
    "\n"
    " Numerically stable: avoids the cancellation of `Σ(x - μ)³` over a list.\n"
).
-spec skewness(list(float())) -> {ok, float()} | {error, nil}.
skewness(Xs) ->
    viva_math@precision:moments_skewness(
        viva_math@precision:moments_from_list(Xs)
    ).

-file("src/viva_math/statistics.gleam", 397).
?DOC(" Excess kurtosis via Pébay online accumulator. Zero for a Gaussian.\n").
-spec kurtosis(list(float())) -> {ok, float()} | {error, nil}.
kurtosis(Xs) ->
    viva_math@precision:moments_excess_kurtosis(
        viva_math@precision:moments_from_list(Xs)
    ).