viva_math

-module(viva_math@entropy).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/entropy.gleam").
-export([shannon/1, shannon_normalized/1, kl_divergence/2, kl_divergence_smoothed/3, symmetric_kl/2, jensen_shannon/2, cross_entropy/2, binary_cross_entropy/2, mutual_information/3, conditional_entropy/2, relative_entropy_rate/2, hybrid_shannon/3, kl_divergence_with_sensitivity/3, renyi/2, tsallis/2, fisher_information_gaussian/1, differential_entropy_gaussian/1]).
-export_type([kl_sensitivity/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(
    " Entropy and information theory functions.\n"
    "\n"
    " Based on Shannon (1948) and Kullback-Leibler (1951).\n"
    " Used for memory consolidation scoring and uncertainty quantification.\n"
    "\n"
    " References:\n"
    " - Shannon (1948) \"A Mathematical Theory of Communication\"\n"
    " - Cover & Thomas (2006) \"Elements of Information Theory\"\n"
).

-type kl_sensitivity() :: standard |
    {arousal_weighted, float()} |
    {custom_gamma, float()}.

-file("src/viva_math/entropy.gleam", 27).
?DOC(
    " Shannon entropy: H(X) = -Σ p(x) log₂ p(x)\n"
    "\n"
    " Measures uncertainty/information content of a distribution.\n"
    " Higher entropy = more uncertainty.\n"
    "\n"
    " ## Examples\n"
    "\n"
    " ```gleam\n"
    " shannon([0.5, 0.5])     // -> 1.0 (maximum for 2 outcomes)\n"
    " shannon([1.0, 0.0])     // -> 0.0 (no uncertainty)\n"
    " shannon([0.25, 0.25, 0.25, 0.25]) // -> 2.0\n"
    " ```\n"
).
-spec shannon(list(float())) -> float().
shannon(Probabilities) ->
    Terms = gleam@list:fold(Probabilities, [], fun(Acc, P) -> case P =< +0.0 of
                true ->
                    Acc;

                false ->
                    case gleam_community@maths:logarithm_2(P) of
                        {ok, Log_p} ->
                            [+0.0 - (P * Log_p) | Acc];

                        {error, _} ->
                            Acc
                    end
            end end),
    viva_math@precision:neumaier_sum(Terms).

-file("src/viva_math/entropy.gleam", 274).
-spec int_to_float(integer()) -> float().
int_to_float(N) ->
    case N of
        0 ->
            +0.0;

        1 ->
            1.0;

        2 ->
            2.0;

        _ ->
            Half = N div 2,
            Remainder = N - (Half * 2),
            (int_to_float(Half) * 2.0) + int_to_float(Remainder)
    end.

-file("src/viva_math/entropy.gleam", 47).
?DOC(
    " Normalized Shannon entropy (0 to 1 range).\n"
    "\n"
    " Divides by log₂(n) where n is number of outcomes.\n"
).
-spec shannon_normalized(list(float())) -> float().
shannon_normalized(Probabilities) ->
    N = erlang:length(Probabilities),
    case N =< 1 of
        true ->
            +0.0;

        false ->
            H = shannon(Probabilities),
            case gleam_community@maths:logarithm_2(int_to_float(N)) of
                {ok, Max_h} ->
                    case Max_h =:= +0.0 of
                        true ->
                            +0.0;

                        false ->
                            case Max_h of
                                +0.0 -> +0.0;
                                -0.0 -> -0.0;
                                Gleam@denominator -> H / Gleam@denominator
                            end
                    end;

                {error, _} ->
                    +0.0
            end
    end.

-file("src/viva_math/entropy.gleam", 71).
?DOC(
    " KL Divergence: D_KL(P || Q) = Σ p(x) log(p(x) / q(x))\n"
    "\n"
    " Measures how much P diverges from Q (not symmetric!).\n"
    " P is the \"true\" distribution, Q is the approximation.\n"
    "\n"
    " Returns Error if distributions have different lengths or Q has zeros where P is non-zero.\n"
).
-spec kl_divergence(list(float()), list(float())) -> {ok, float()} |
    {error, nil}.
kl_divergence(P, Q) ->
    case erlang:length(P) =:= erlang:length(Q) of
        false ->
            {error, nil};

        true ->
            Pairs = gleam@list:zip(P, Q),
            Result = gleam@list:fold(
                Pairs,
                {ok, +0.0},
                fun(Acc, Pair) -> case Acc of
                        {error, nil} ->
                            {error, nil};

                        {ok, Sum} ->
                            {Pi, Qi} = Pair,
                            case Pi =< +0.0 of
                                true ->
                                    {ok, Sum};

                                false ->
                                    case Qi =< +0.0 of
                                        true ->
                                            {error, nil};

                                        false ->
                                            case gleam_community@maths:natural_logarithm(
                                                case Qi of
                                                    +0.0 -> +0.0;
                                                    -0.0 -> -0.0;
                                                    Gleam@denominator -> Pi / Gleam@denominator
                                                end
                                            ) of
                                                {ok, Log_ratio} ->
                                                    {ok, Sum + (Pi * Log_ratio)};

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

-file("src/viva_math/entropy.gleam", 111).
?DOC(
    " KL divergence with additive smoothing.\n"
    "\n"
    " Adds `eps` to every q_i before normalising; useful when q contains zeros\n"
    " where p > 0 (which would otherwise return `Error`). Standard practice in\n"
    " language modelling and probabilistic ML.\n"
    "\n"
    " Recommended `eps ∈ [1e-12, 1e-6]`. Larger values bias the result more.\n"
).
-spec kl_divergence_smoothed(list(float()), list(float()), float()) -> {ok,
        float()} |
    {error, nil}.
kl_divergence_smoothed(P, Q, Eps) ->
    case {erlang:length(P) =:= erlang:length(Q), Eps < +0.0} of
        {false, _} ->
            {error, nil};

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

        {true, false} ->
            Smoothed_unnorm = gleam@list:map(Q, fun(Qi) -> Qi + Eps end),
            Total = gleam@list:fold(
                Smoothed_unnorm,
                +0.0,
                fun(Acc, X) -> Acc + X end
            ),
            case Total =< +0.0 of
                true ->
                    {error, nil};

                false ->
                    Q_norm = gleam@list:map(
                        Smoothed_unnorm,
                        fun(X@1) -> case Total of
                                +0.0 -> +0.0;
                                -0.0 -> -0.0;
                                Gleam@denominator -> X@1 / Gleam@denominator
                            end end
                    ),
                    kl_divergence(P, Q_norm)
            end
    end.

-file("src/viva_math/entropy.gleam", 137).
?DOC(
    " Symmetric KL Divergence (Jensen-Shannon divergence without the 1/2).\n"
    "\n"
    " D_sym(P, Q) = D_KL(P || Q) + D_KL(Q || P)\n"
).
-spec symmetric_kl(list(float()), list(float())) -> {ok, float()} | {error, nil}.
symmetric_kl(P, Q) ->
    case {kl_divergence(P, Q), kl_divergence(Q, P)} of
        {{ok, D1}, {ok, D2}} ->
            {ok, D1 + D2};

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

-file("src/viva_math/entropy.gleam", 148).
?DOC(
    " Jensen-Shannon Divergence: JS(P, Q) = (D_KL(P || M) + D_KL(Q || M)) / 2\n"
    " where M = (P + Q) / 2\n"
    "\n"
    " This is symmetric and bounded [0, 1] when using log₂.\n"
).
-spec jensen_shannon(list(float()), list(float())) -> {ok, float()} |
    {error, nil}.
jensen_shannon(P, Q) ->
    case erlang:length(P) =:= erlang:length(Q) of
        false ->
            {error, nil};

        true ->
            M = begin
                _pipe = gleam@list:zip(P, Q),
                gleam@list:map(
                    _pipe,
                    fun(Pair) ->
                        {Pi, Qi} = Pair,
                        (Pi + Qi) / 2.0
                    end
                )
            end,
            case {kl_divergence(P, M), kl_divergence(Q, M)} of
                {{ok, D_pm}, {ok, D_qm}} ->
                    {ok, (D_pm + D_qm) / 2.0};

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

-file("src/viva_math/entropy.gleam", 172).
?DOC(
    " Cross-entropy: H(P, Q) = -Σ p(x) log q(x)\n"
    "\n"
    " Used in machine learning loss functions.\n"
    " H(P, Q) = H(P) + D_KL(P || Q)\n"
).
-spec cross_entropy(list(float()), list(float())) -> {ok, float()} |
    {error, nil}.
cross_entropy(P, Q) ->
    case erlang:length(P) =:= erlang:length(Q) of
        false ->
            {error, nil};

        true ->
            Pairs = gleam@list:zip(P, Q),
            Result = gleam@list:fold(
                Pairs,
                {ok, +0.0},
                fun(Acc, Pair) -> case Acc of
                        {error, nil} ->
                            {error, nil};

                        {ok, Sum} ->
                            {Pi, Qi} = Pair,
                            case Pi =< +0.0 of
                                true ->
                                    {ok, Sum};

                                false ->
                                    case Qi =< +0.0 of
                                        true ->
                                            {error, nil};

                                        false ->
                                            case gleam_community@maths:natural_logarithm(
                                                Qi
                                            ) of
                                                {ok, Log_q} ->
                                                    {ok, Sum - (Pi * Log_q)};

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

-file("src/viva_math/entropy.gleam", 206).
?DOC(
    " Binary cross-entropy for single probability.\n"
    "\n"
    " H(p, q) = -[p log(q) + (1-p) log(1-q)]\n"
).
-spec binary_cross_entropy(float(), float()) -> {ok, float()} | {error, nil}.
binary_cross_entropy(P, Q) ->
    Q_clamped = gleam@float:max(gleam@float:min(Q, 0.999999), 0.000001),
    Q_inv = 1.0 - Q_clamped,
    P_inv = 1.0 - P,
    case {gleam_community@maths:natural_logarithm(Q_clamped),
        gleam_community@maths:natural_logarithm(Q_inv)} of
        {{ok, Log_q}, {ok, Log_q_inv}} ->
            Result = (P * Log_q) + (P_inv * Log_q_inv),
            {ok, +0.0 - Result};

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

-file("src/viva_math/entropy.gleam", 225).
?DOC(
    " Mutual Information: I(X; Y) = H(X) + H(Y) - H(X, Y)\n"
    "\n"
    " Measures shared information between two variables.\n"
    " Takes marginal distributions and joint distribution as input.\n"
).
-spec mutual_information(list(float()), list(float()), list(list(float()))) -> float().
mutual_information(Px, Py, Pxy) ->
    Hx = shannon(Px),
    Hy = shannon(Py),
    Pxy_flat = lists:append(Pxy),
    Hxy = shannon(Pxy_flat),
    (Hx + Hy) - Hxy.

-file("src/viva_math/entropy.gleam", 243).
?DOC(
    " Conditional entropy: H(X|Y) = H(X, Y) - H(Y)\n"
    "\n"
    " Uncertainty in X given knowledge of Y.\n"
).
-spec conditional_entropy(list(float()), list(list(float()))) -> float().
conditional_entropy(Px, Pxy) ->
    Hx = shannon(Px),
    Pxy_flat = lists:append(Pxy),
    Hxy = shannon(Pxy_flat),
    Hxy - Hx.

-file("src/viva_math/entropy.gleam", 253).
?DOC(
    " Relative entropy rate for sequences.\n"
    "\n"
    " Used for measuring \"surprise\" in temporal data.\n"
).
-spec relative_entropy_rate(list(float()), list(float())) -> {ok, float()} |
    {error, nil}.
relative_entropy_rate(Observed, Expected) ->
    case erlang:length(Observed) =:= erlang:length(Expected) of
        false ->
            {error, nil};

        true ->
            N = erlang:length(Observed),
            case N =:= 0 of
                true ->
                    {ok, +0.0};

                false ->
                    case kl_divergence(Observed, Expected) of
                        {ok, Kl} ->
                            {ok, case int_to_float(N) of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator -> Kl / Gleam@denominator
                                end};

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

-file("src/viva_math/entropy.gleam", 315).
?DOC(" Clamp value to [0, 1]\n").
-spec clamp_unit(float()) -> float().
clamp_unit(Value) ->
    case Value < +0.0 of
        true ->
            +0.0;

        false ->
            case Value > 1.0 of
                true ->
                    1.0;

                false ->
                    Value
            end
    end.

-file("src/viva_math/entropy.gleam", 303).
?DOC(
    " Hybrid entropy for mixed emotional states.\n"
    "\n"
    " H_hybrid(X) = α × H(X₁) + (1 - α) × H(X₂)\n"
    "\n"
    " Proposed by DeepSeek R1 for modeling hybrid emotions.\n"
    " α ∈ [0, 1] controls the blend between two emotional distributions.\n"
    "\n"
    " ## Examples\n"
    "\n"
    " ```gleam\n"
    " hybrid_shannon([0.5, 0.5], [0.7, 0.3], 0.5)  // Blend of two emotions\n"
    " ```\n"
).
-spec hybrid_shannon(list(float()), list(float()), float()) -> float().
hybrid_shannon(Probs1, Probs2, Alpha) ->
    H1 = shannon(Probs1),
    H2 = shannon(Probs2),
    Clamped_alpha = clamp_unit(Alpha),
    (Clamped_alpha * H1) + ((1.0 - Clamped_alpha) * H2).

-file("src/viva_math/entropy.gleam", 389).
-spec weighted_mean_helper(list(float()), integer(), float(), float()) -> float().
weighted_mean_helper(Probs, Index, Sum, Weight_sum) ->
    case Probs of
        [] ->
            case Weight_sum =:= +0.0 of
                true ->
                    +0.0;

                false ->
                    case Weight_sum of
                        +0.0 -> +0.0;
                        -0.0 -> -0.0;
                        Gleam@denominator -> Sum / Gleam@denominator
                    end
            end;

        [P | Rest] ->
            Idx_float = int_to_float(Index),
            weighted_mean_helper(
                Rest,
                Index + 1,
                Sum + (P * Idx_float),
                Weight_sum + P
            )
    end.

-file("src/viva_math/entropy.gleam", 385).
?DOC(
    " Calculate weighted mean of a distribution.\n"
    " Assumes indices 0, 1, 2, ... as values.\n"
).
-spec weighted_mean(list(float())) -> float().
weighted_mean(Probs) ->
    weighted_mean_helper(Probs, 0, +0.0, +0.0).

-file("src/viva_math/entropy.gleam", 376).
?DOC(" Calculate squared difference of distribution means.\n").
-spec mean_difference_squared(list(float()), list(float())) -> float().
mean_difference_squared(P, Q) ->
    Mean_p = weighted_mean(P),
    Mean_q = weighted_mean(Q),
    Diff = Mean_p - Mean_q,
    Diff * Diff.

-file("src/viva_math/entropy.gleam", 344).
?DOC(
    " KL divergence with sensitivity parameter.\n"
    "\n"
    " D_KL^γ(P || Q) = γ × (μ₁ - μ₂)² + D_KL(P || Q)\n"
    "\n"
    " Proposed by DeepSeek R1 for arousal-modulated divergence.\n"
    " Higher γ = more sensitive to mean differences.\n"
).
-spec kl_divergence_with_sensitivity(
    list(float()),
    list(float()),
    kl_sensitivity()
) -> {ok, float()} | {error, nil}.
kl_divergence_with_sensitivity(P, Q, Sensitivity) ->
    Gamma = case Sensitivity of
        standard ->
            +0.0;

        {arousal_weighted, Arousal} ->
            Abs_arousal = case Arousal >= +0.0 of
                true ->
                    Arousal;

                false ->
                    +0.0 - Arousal
            end,
            Abs_arousal;

        {custom_gamma, G} ->
            G
    end,
    case kl_divergence(P, Q) of
        {error, nil} ->
            {error, nil};

        {ok, Standard_kl} ->
            Mean_diff_sq = mean_difference_squared(P, Q),
            {ok, (Gamma * Mean_diff_sq) + Standard_kl}
    end.

-file("src/viva_math/entropy.gleam", 448).
?DOC(" Power function using exp(α × ln(x))\n").
-spec power(float(), float()) -> float().
power(Base, Exponent) ->
    case Base =< +0.0 of
        true ->
            +0.0;

        false ->
            case gleam_community@maths:natural_logarithm(Base) of
                {ok, Ln_base} ->
                    gleam_community@maths:exponential(Exponent * Ln_base);

                {error, _} ->
                    +0.0
            end
    end.

-file("src/viva_math/entropy.gleam", 421).
?DOC(
    " Renyi entropy of order α.\n"
    "\n"
    " H_α(X) = (1/(1-α)) × log(Σ p(x)^α)\n"
    "\n"
    " Generalizes Shannon entropy (α → 1 gives Shannon).\n"
    " α = 0: Hartley entropy (log of support size)\n"
    " α = 2: Collision entropy\n"
    " α → ∞: Min-entropy\n"
).
-spec renyi(list(float()), float()) -> {ok, float()} | {error, nil}.
renyi(Probabilities, Alpha) ->
    case Alpha =:= 1.0 of
        true ->
            {ok, shannon(Probabilities)};

        false ->
            Sum_p_alpha = gleam@list:fold(
                Probabilities,
                +0.0,
                fun(Acc, P) -> case P =< +0.0 of
                        true ->
                            Acc;

                        false ->
                            Acc + power(P, Alpha)
                    end end
            ),
            case Sum_p_alpha =< +0.0 of
                true ->
                    {error, nil};

                false ->
                    case gleam_community@maths:logarithm_2(Sum_p_alpha) of
                        {ok, Log_sum} ->
                            {ok, case (1.0 - Alpha) of
                                    +0.0 -> +0.0;
                                    -0.0 -> -0.0;
                                    Gleam@denominator -> Log_sum / Gleam@denominator
                                end};

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

-file("src/viva_math/entropy.gleam", 469).
?DOC(
    " Tsallis entropy S_q(X) = (1 - Σ pᵢ^q) / (q - 1).\n"
    "\n"
    " Generalises Shannon (q → 1) and Rényi. q < 1 emphasises rare events,\n"
    " q > 1 emphasises modes. Used in non-extensive statistical mechanics and\n"
    " to model heavy-tailed emotional distributions.\n"
).
-spec tsallis(list(float()), float()) -> {ok, float()} | {error, nil}.
tsallis(Probabilities, Q) ->
    Q_close_to_1 = case Q - 1.0 of
        D when D < +0.0 ->
            (+0.0 - D) < 1.0e-9;

        D@1 ->
            D@1 < 1.0e-9
    end,
    case Q_close_to_1 of
        true ->
            {ok, shannon(Probabilities)};

        false ->
            Sum_p_q = gleam@list:fold(
                Probabilities,
                +0.0,
                fun(Acc, P) -> case P =< +0.0 of
                        true ->
                            Acc;

                        false ->
                            Acc + power(P, Q)
                    end end
            ),
            {ok, case (Q - 1.0) of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> (1.0 - Sum_p_q) / Gleam@denominator
                end}
    end.

-file("src/viva_math/entropy.gleam", 494).
?DOC(
    " Fisher information for a Gaussian distribution: I(σ) = 1 / σ².\n"
    "\n"
    " Measures how much information an observation carries about the mean.\n"
).
-spec fisher_information_gaussian(float()) -> {ok, float()} | {error, nil}.
fisher_information_gaussian(Sigma) ->
    case Sigma =< +0.0 of
        true ->
            {error, nil};

        false ->
            {ok, case (Sigma * Sigma) of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> 1.0 / Gleam@denominator
                end}
    end.

-file("src/viva_math/entropy.gleam", 502).
?DOC(" Differential entropy of a Gaussian: h(N(μ, σ²)) = ½ ln(2πeσ²).\n").
-spec differential_entropy_gaussian(float()) -> {ok, float()} | {error, nil}.
differential_entropy_gaussian(Sigma) ->
    case Sigma =< +0.0 of
        true ->
            {error, nil};

        false ->
            Two_pi_e = 17.07946844534713,
            case gleam_community@maths:natural_logarithm(
                (Two_pi_e * Sigma) * Sigma
            ) of
                {ok, Ln_val} ->
                    {ok, 0.5 * Ln_val};

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