-module(viva_math@cusp).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/cusp.gleam").
-export([from_arousal_dominance/2, potential/2, gradient/2, discriminant/1, is_bistable/1, equilibria/1, nearest_equilibrium/2, would_jump/2, volatility/2]).
-export_type([cusp_params/0, cusp_result/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(
    " Cusp Catastrophe Theory implementation.\n"
    "\n"
    " Based on René Thom's catastrophe theory (1972).\n"
    " The cusp is the simplest catastrophe with two control parameters.\n"
    "\n"
    " Potential function: V(x) = x⁴/4 + αx²/2 + βx\n"
    " Gradient: dV/dx = x³ + αx + β\n"
    " Discriminant: Δ = -4α³ - 27β²\n"
    "\n"
    " When Δ > 0 and α < 0: bistable region (two stable states)\n"
    " This models emotional \"phase transitions\" - sudden mood shifts.\n"
    "\n"
    " References:\n"
    " - Grasman et al. (2009) \"Fitting the Cusp Catastrophe in R\"\n"
    " - Van der Maas et al. (2003) \"Sudden Transitions in Attitudes\"\n"
).

-type cusp_params() :: {cusp_params, float(), float()}.

-type cusp_result() :: {monostable, float()} |
    {bistable, float(), float(), float()}.

-file("src/viva_math/cusp.gleam", 43).
?DOC(
    " Create cusp parameters from PAD arousal and dominance.\n"
    "\n"
    " Mapping (from VIVA emotional dynamics):\n"
    " - alpha = -arousal (high arousal → negative alpha → bistability)\n"
    " - beta = dominance * 0.5 (dominance biases equilibrium)\n"
    "\n"
    " This means high arousal creates emotional volatility.\n"
).
-spec from_arousal_dominance(float(), float()) -> cusp_params().
from_arousal_dominance(Arousal, Dominance) ->
    Alpha = +0.0 - Arousal,
    Beta = Dominance * 0.5,
    {cusp_params, Alpha, Beta}.

-file("src/viva_math/cusp.gleam", 53).
?DOC(
    " Compute the cusp potential V(x) = x⁴/4 + αx²/2 + βx\n"
    "\n"
    " The potential represents emotional \"energy landscape\".\n"
    " Stable states are at local minima of this function.\n"
).
-spec potential(float(), cusp_params()) -> float().
potential(X, Params) ->
    X2 = X * X,
    X4 = X2 * X2,
    ((X4 / 4.0) + ((erlang:element(2, Params) * X2) / 2.0)) + (erlang:element(
        3,
        Params
    )
    * X).

-file("src/viva_math/cusp.gleam", 62).
?DOC(
    " Compute the gradient dV/dx = x³ + αx + β\n"
    "\n"
    " Equilibria are where gradient = 0.\n"
).
-spec gradient(float(), cusp_params()) -> float().
gradient(X, Params) ->
    X3 = (X * X) * X,
    (X3 + (erlang:element(2, Params) * X)) + erlang:element(3, Params).

-file("src/viva_math/cusp.gleam", 71).
?DOC(
    " Compute the discriminant Δ = -4α³ - 27β²\n"
    "\n"
    " Δ > 0 and α < 0: bistable region\n"
    " Δ ≤ 0 or α ≥ 0: monostable region\n"
).
-spec discriminant(cusp_params()) -> float().
discriminant(Params) ->
    Alpha3 = (erlang:element(2, Params) * erlang:element(2, Params)) * erlang:element(
        2,
        Params
    ),
    Beta2 = erlang:element(3, Params) * erlang:element(3, Params),
    (-4.0 * Alpha3) - (27.0 * Beta2).

-file("src/viva_math/cusp.gleam", 82).
?DOC(
    " Check if the system is in bistable region.\n"
    "\n"
    " Bistability requires:\n"
    " 1. α < 0 (necessary for two minima)\n"
    " 2. Δ > 0 (discriminant positive)\n"
).
-spec is_bistable(cusp_params()) -> boolean().
is_bistable(Params) ->
    (erlang:element(2, Params) < +0.0) andalso (discriminant(Params) > +0.0).

-file("src/viva_math/cusp.gleam", 217).
-spec trigonometric_roots(float(), float()) -> list(float()).
trigonometric_roots(Alpha, Beta) ->
    P = Alpha,
    Q = Beta,
    case P >= +0.0 of
        true ->
            [];

        false ->
            Neg_p = +0.0 - P,
            Neg_p_over_3 = Neg_p / 3.0,
            M = case gleam_community@maths:nth_root(Neg_p_over_3, 2) of
                {ok, V} ->
                    2.0 * V;

                {error, _} ->
                    +0.0
            end,
            Neg_p_cubed = (Neg_p_over_3 * Neg_p_over_3) * Neg_p_over_3,
            Cos_arg = case gleam_community@maths:nth_root(Neg_p_cubed, 2) of
                {ok, Denom} ->
                    case Denom =:= +0.0 of
                        true ->
                            +0.0;

                        false ->
                            Numerator = 3.0 * Q,
                            Denominator = (2.0 * P) * Denom,
                            Raw = case Denominator of
                                +0.0 -> +0.0;
                                -0.0 -> -0.0;
                                Gleam@denominator -> Numerator / Gleam@denominator
                            end,
                            +0.0 - Raw
                    end;

                {error, _} ->
                    +0.0
            end,
            Cos_arg_clamped = viva_math@common:clamp(Cos_arg, -1.0, 1.0),
            Theta = case gleam_community@maths:acos(Cos_arg_clamped) of
                {ok, V@1} ->
                    V@1 / 3.0;

                {error, _} ->
                    +0.0
            end,
            Pi = gleam_community@maths:pi(),
            Two_pi_over_3 = (2.0 * Pi) / 3.0,
            Four_pi_over_3 = (4.0 * Pi) / 3.0,
            R1 = M * gleam_community@maths:cos(Theta),
            R2 = M * gleam_community@maths:cos(Theta - Two_pi_over_3),
            R3 = M * gleam_community@maths:cos(Theta - Four_pi_over_3),
            [R1, R2, R3]
    end.

-file("src/viva_math/cusp.gleam", 273).
-spec cbrt(float()) -> float().
cbrt(X) ->
    case X >= +0.0 of
        true ->
            case gleam_community@maths:nth_root(X, 3) of
                {ok, V} ->
                    V;

                {error, _} ->
                    +0.0
            end;

        false ->
            Neg_x = +0.0 - X,
            case gleam_community@maths:nth_root(Neg_x, 3) of
                {ok, V@1} ->
                    +0.0 - V@1;

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

-file("src/viva_math/cusp.gleam", 178).
-spec cardano_single_root(float(), float()) -> float().
cardano_single_root(Alpha, Beta) ->
    P = Alpha,
    Q = Beta,
    P_over_3 = P / 3.0,
    Q_over_2 = Q / 2.0,
    P3 = (P_over_3 * P_over_3) * P_over_3,
    Q2 = Q_over_2 * Q_over_2,
    D = Q2 + P3,
    case D >= +0.0 of
        true ->
            Sqrt_d = case gleam_community@maths:nth_root(D, 2) of
                {ok, V} ->
                    V;

                {error, _} ->
                    +0.0
            end,
            Neg_q_over_2 = +0.0 - Q_over_2,
            U = cbrt(Neg_q_over_2 + Sqrt_d),
            V@1 = cbrt(Neg_q_over_2 - Sqrt_d),
            U + V@1;

        false ->
            case trigonometric_roots(Alpha, Beta) of
                [R | _] ->
                    R;

                [] ->
                    +0.0
            end
    end.

-file("src/viva_math/cusp.gleam", 291).
-spec sort_three(float(), float(), float()) -> {float(), float(), float()}.
sort_three(A, B, C) ->
    case A =< B of
        true ->
            case B =< C of
                true ->
                    {A, B, C};

                false ->
                    case A =< C of
                        true ->
                            {A, C, B};

                        false ->
                            {C, A, B}
                    end
            end;

        false ->
            case A =< C of
                true ->
                    {B, A, C};

                false ->
                    case B =< C of
                        true ->
                            {B, C, A};

                        false ->
                            {C, B, A}
                    end
            end
    end.

-file("src/viva_math/cusp.gleam", 90).
?DOC(
    " Calculate equilibria (roots of the cubic x³ + αx + β = 0).\n"
    "\n"
    " Uses Cardano's formula for the depressed cubic.\n"
    " Returns Monostable or Bistable depending on discriminant.\n"
).
-spec equilibria(cusp_params()) -> cusp_result().
equilibria(Params) ->
    _ = discriminant(Params),
    case is_bistable(Params) of
        false ->
            Root = cardano_single_root(
                erlang:element(2, Params),
                erlang:element(3, Params)
            ),
            {monostable, Root};

        true ->
            Roots = trigonometric_roots(
                erlang:element(2, Params),
                erlang:element(3, Params)
            ),
            case Roots of
                [R1, R2, R3] ->
                    Sorted = sort_three(R1, R2, R3),
                    {bistable,
                        erlang:element(1, Sorted),
                        erlang:element(2, Sorted),
                        erlang:element(3, Sorted)};

                _ ->
                    {monostable,
                        cardano_single_root(
                            erlang:element(2, Params),
                            erlang:element(3, Params)
                        )}
            end
    end.

-file("src/viva_math/cusp.gleam", 118).
?DOC(" Find nearest stable equilibrium to current state.\n").
-spec nearest_equilibrium(float(), cusp_params()) -> float().
nearest_equilibrium(X, Params) ->
    case equilibria(Params) of
        {monostable, Eq} ->
            Eq;

        {bistable, Lower, _, Upper} ->
            Dist_lower = gleam@float:absolute_value(X - Lower),
            Dist_upper = gleam@float:absolute_value(X - Upper),
            case Dist_lower < Dist_upper of
                true ->
                    Lower;

                false ->
                    Upper
            end
    end.

-file("src/viva_math/cusp.gleam", 136).
?DOC(
    " Check if state would \"jump\" to other attractor.\n"
    "\n"
    " In bistable regime, if state crosses the unstable equilibrium,\n"
    " it will rapidly transition to the opposite stable state.\n"
).
-spec would_jump(float(), cusp_params()) -> boolean().
would_jump(X, Params) ->
    case equilibria(Params) of
        {monostable, _} ->
            false;

        {bistable, Lower, Unstable, Upper} ->
            At_lower = gleam@float:absolute_value(X - Lower) < gleam@float:absolute_value(
                X - Upper
            ),
            case At_lower of
                true ->
                    X > Unstable;

                false ->
                    X < Unstable
            end
    end.

-file("src/viva_math/cusp.gleam", 159).
?DOC(
    " Compute emotional volatility based on cusp geometry.\n"
    "\n"
    " Volatility is high when:\n"
    " 1. In bistable region\n"
    " 2. Close to the unstable equilibrium (catastrophe manifold)\n"
).
-spec volatility(float(), cusp_params()) -> float().
volatility(X, Params) ->
    case is_bistable(Params) of
        false ->
            +0.0;

        true ->
            case equilibria(Params) of
                {monostable, _} ->
                    +0.0;

                {bistable, _, Unstable, _} ->
                    Dist = gleam@float:absolute_value(X - Unstable),
                    Neg_dist_sq = +0.0 - (Dist * Dist),
                    gleam_community@maths:exponential(Neg_dist_sq)
            end
    end.