//! HAAP — Hyperpolarization-Activated Adaptive Pulse-Bursting
//!
//! Uses `#` + `ComputeNeuron` for typed API.
//!
//! # Principles
//! - **Ohm's Law**: I = g·(V − E_rev), subtracted from I_inj
//! - **Forward Euler**: no gate_integrate, no hidden abstractions
//! - **Zero divisions**: all tau stored as inv_tau at init
//!
//! Equations:
//! ```text
//! I_leak = g_L·(V − E_L)
//! I_M = g_M·u·(V − E_K)
//! I_h = g_h·h·(V − E_h)
//! C_m·dV/dt = I_syn − I_leak − I_M − I_h
//! du/dt = (u_∞(V) − u) · inv_tau_u (Forward Euler)
//! dh/dt = (h_∞(V) − h) · inv_tau_h (Forward Euler)
//! Spike: V ≥ V_thresh → V = V_reset, u += delta_u
//! ```
//!
//! # Gate steady states (standard Boltzmann — `k` is the slope factor, a divisor)
//! ```text
//! u_∞(V) = 1 / (1 + exp(−(V − V_half_u) / k_u)) M-current (opens on depolarization)
//! h_∞(V) = 1 / (1 + exp(+(V − V_half_h) / k_h)) HCN (opens on hyperpolarization)
//! ```
//! This is the biologically correct convention, matching the model's name and
//! purpose: `h` gives hyperpolarization-activated rebound, `u` gives
//! depolarization-activated spike-frequency adaptation. Because `k` divides,
//! the exponent argument stays small for physiological `V` (|arg| ≲ 40) and
//! cannot overflow.

#! // deliberate: `pub mod haap` mirrors the file module

use pillars_macros::dual_neuron;

// ═══════════════════════════════════════════════════════════════
// ═══════════════════════════════════════════════════════════════

/// Hyperpolarization-activated adaptive bursting neuron.
///
/// # States (3)
/// - `v`: membrane potential (mV)
/// - `u`: M-current activation (0..1)
/// - `h`: h-current activation (0..1)
///
/// # Params (17)
/// - `dt`, `inv_cm`, `g_l`, `e_l`, `g_m`, `e_k`, `g_h`, `e_h`,
/// `inv_tau_u`, `inv_tau_h`, `v_half_u`, `k_u`, `v_half_h`,
/// `k_h`, `v_thresh`, `v_reset`, `delta_u`
#
pub mod haap {
use crate::core::compute::{ComputeNeuron, StepContext};
use crate::core::lane::SimdLane;
use pillars_macros::NeuronModel;

#
pub struct HAAPModel {
#
pub v: f32,
#
pub u: f32,
#
pub h: f32,

#
pub dt: f32,
#
pub inv_cm: f32,
#
pub g_l: f32,
#
pub e_l: f32,
#
pub g_m: f32,
#
pub e_k: f32,
#
pub g_h: f32,
#
pub e_h: f32,
#
pub inv_tau_u: f32,
#
pub inv_tau_h: f32,
#
pub v_half_u: f32,
#
pub k_u: f32,
#
pub v_half_h: f32,
#
pub k_h: f32,
#
pub v_thresh: f32,
#
pub v_reset: f32,
#
pub delta_u: f32,
}

impl ComputeNeuron for HAAPModel {
type State<S: SimdLane> = HAAPState<S>;
type Params<S: SimdLane> = HAAPParams<S>;

#
fn step<S: SimdLane>(
state: &mut HAAPState<S>,
params: &HAAPParams<S>,
ctx: &mut StepContext<S>,
) {
let one = S::one();
let dt = ctx.dt_simd();

// ─── Ohmic currents: I = g·(V − E), subtracted from I_inj ──
let i_leak = params.g_l * (state.v - params.e_l);
let i_m = params.g_m * state.u * (state.v - params.e_k);
let i_h = params.g_h * state.h * (state.v - params.e_h);

let i_total = i_leak + i_m + i_h;
let dv = (ctx.i_syn - i_total) * params.inv_cm;
let v_new = state.v + dv * dt;

// ─── Gate kinetics (standard Boltzmann: `k` = slope factor, a divisor) ──
// u_∞ = 1/(1+exp(−(V − V_half_u)/k_u)) M-current opens on depolarization
// h_∞ = 1/(1+exp(+(V − V_half_h)/k_h)) HCN opens on hyperpolarization
//
// With `k` dividing, |arg| stays ≲ 40 for physiological V, so exp
// cannot overflow — the earlier "steepness" rewrite (k as a
// multiplier) produced ±300 arguments, which was the real source of
// the overflow it then "fixed" with clamping. The ±50 clamp below
// is a pure safety net for pathological V (far outside biophysics).
let u_arg = (-(v_new - params.v_half_u) / params.k_u)
.max(S::splat(-50.0))
.min(S::splat(50.0));
let u_inf = one / (one + u_arg.exp());
let u_new = (state.u + dt * (u_inf - state.u) * params.inv_tau_u).clamp01();

let h_arg = ((v_new - params.v_half_h) / params.k_h)
.max(S::splat(-50.0))
.min(S::splat(50.0));
let h_inf = one / (one + h_arg.exp());
let h_new = (state.h + dt * (h_inf - state.h) * params.inv_tau_h).clamp01();

// ─── Spike detection & reset ───────────────────────────────
let spike = v_new.cmp_ge(params.v_thresh);
ctx.emit_spike(spike);

state.v = spike.cond().then(params.v_reset).else_(v_new);
let u_bump = (u_new + params.delta_u).clamp01();
state.u = spike.cond().then(u_bump).else_(u_new);
state.h = h_new;
}
}
}

// Re-export so existing paths keep working:
// `pillars::zoo::neurons::haap::HAAPModel` etc.
pub use haap::*;