tulip_rs 0.1.15

//! # MESA Sine Wave — Optimised Scalar Implementation
//!
//! A clean reimplementation of [`super::msw`] that applies all three proposals
//! from `msw_optimization_proposal.md`:
//!
//! | Proposal | Change                               | Where applied                     |
//! |----------|--------------------------------------|-----------------------------------|
//! | **P3**   | `sin_cos` + angle-addition for lead  | `phase_from_rp_ip`                |
//! | **P1**   | Precomputed twiddle tables           | `precompute_twiddles`, `dot_product_simd` |
//! | **P2**   | Sliding DFT — O(1) per output bar    | `cycle_msw_sdft`, `batch_indicator` |
//!
//! This module is intentionally **independent of `msw`** so that `adaptivemsw`
//! (which calls `msw::calc` with a varying period each bar) is unaffected until
//! that migration is validated separately.
//!
//! ## Struct layout (mirrors the `stoch` pattern)
//!
//! | Struct           | Purpose                                                     |
//! |------------------|-------------------------------------------------------------|
//! | [`State`]        | Minimal per-bar update state — `(rp, ip, wr, wi)`          |
//! | [`IndicatorState`] | Cross-call state — price tail + cached twiddles + `State` |
//!
//! ## Algorithm summary
//!
//! **Initial call (`indicator`)**
//! 1. Precompute `period` twiddle factors with scalar `sin_cos` once → stored in
//!    `IndicatorState` as flat `Vec<f64>` (serializable, reusable).
//! 2. Seed the SDFT with one full dot-product over the first window.
//! 3. For every subsequent output bar, apply the O(1) SDFT recurrence via `calc`.
//!
//! **Streaming (`batch_indicator`)**
//! - Resume the SDFT from the stored `State` — zero trig per bar.
//! - Every `RE_ANCHOR_INTERVAL` bars a full dot-product is recomputed using the
//!   *cached twiddles* (no trig — pure FMA) to prevent floating-point drift.

use crate::common::{validate_inputs, validate_options};
pub use crate::indicator_types::TIndicatorState;
use crate::types::{DisplayGroup, DisplayType, IndicatorError, IndicatorType, Info};
use serde::{Deserialize, Serialize};
use std::f64::consts::PI;
use std::simd::{num::SimdFloat, Simd, StdFloat};

// ── Compile-time twiddle tables (generated by build.rs) ───────────────────────

include!(concat!(env!("OUT_DIR"), "/msw_twiddles.rs"));

/// Returns the precomputed cos/sin twiddle slices for a given `period` in `[6, 50]`.
///
/// Both slices have length `period`. The data lives in `.rodata` (no heap, no trig).
/// Used by `adaptivemsw` to make every per-bar DFT call allocation-free.
///
/// # Panics (debug only)
/// If `period` is outside `[6, 50]`.
#[inline(always)]
pub fn twiddles_for_period(period: usize) -> (&'static [f64], &'static [f64]) {
    debug_assert!((6..=50).contains(&period), "period {period} outside [6, 50]");
    let idx = period.saturating_sub(6).min(44);
    (&MSW_COS_TWIDDLES[idx][..period], &MSW_SIN_TWIDDLES[idx][..period])
}

/// Number of input price series required.
pub const INPUTS_WIDTH: usize = 1;

/// Number of option parameters required.
pub const OPTIONS_WIDTH: usize = 1;

/// SIMD-parallel variant that processes `N` assets with identical options simultaneously.
/// Requires the `simd_assets` Cargo feature. See [`by_assets`] for the module form.
#[cfg(feature = "simd_assets")]
pub use crate::indicators::simd_indicators::msw_simd::indicator_by_assets;

/// SIMD-parallel variant that processes a single asset with `N` different option
/// sets simultaneously. Requires the `simd_options` Cargo feature. See [`by_options`].
#[cfg(feature = "simd_options")]
pub use crate::indicators::simd_indicators::msw_simd::indicator_by_options;

/// Convenience module that re-exports [`indicator_by_assets`] as `indicator`,
/// allowing SIMD multi-asset computation to be used as a drop-in replacement
/// for the standard single-asset [`indicator`] function.
/// Requires the `simd_assets` Cargo feature.
#[cfg(feature = "simd_assets")]
pub mod by_assets {
    /// Processes `N` assets in parallel with shared options.
    /// See the parent module's [`super::indicator_by_assets`] for full documentation.
    pub use crate::indicators::simd_indicators::msw_simd::indicator_by_assets as indicator;
}

/// Convenience module that re-exports [`indicator_by_options`] as `indicator`,
/// allowing SIMD multi-option computation to be used as a drop-in replacement
/// for the standard single-asset [`indicator`] function.
/// Requires the `simd_options` Cargo feature.
#[cfg(feature = "simd_options")]
pub mod by_options {
    /// Processes a single asset with `N` different option sets in parallel.
    /// See the parent module's [`super::indicator_by_options`] for full documentation.
    pub use crate::indicators::simd_indicators::msw_simd::indicator_by_options as indicator;
}

// ── Scalar constants ──────────────────────────────────────────────────────────

const TPI: f64 = PI * 2.0;
const HPI: f64 = PI * 0.5;
/// 1/√2 — used in the angle-addition identity for sin(phase + π/4).
const INV_SQRT2: f64 = std::f64::consts::FRAC_1_SQRT_2;

// ── Indicator metadata ────────────────────────────────────────────────────────

pub const INFO: Info = Info {
    name: "msw",
    full_name: "Mesa Sine Wave",
    indicator_type: IndicatorType::Cycle,
    inputs: &["real"],
    options: &["period"],
    outputs: &["msw_sine", "msw_lead"],
    optional_outputs: &[],
    display_groups: &[DisplayGroup {
        offset: None,
        id: "msw",
        label: "MSW",
        display_type: DisplayType::Indicator,
        outputs: &["msw_sine", "msw_lead"],
    }],
};

// ── Per-bar state ─────────────────────────────────────────────────────────────

/// Minimal per-bar SDFT state — everything needed to advance one output bar.
///
/// Analogous to the `State` structs in other indicators (e.g. `stoch::State`).
/// Stored inside [`IndicatorState`] and passed to [`calc`] and [`cycle_sdft`].
#[derive(Serialize, Deserialize)]
pub struct State {
    /// SDFT real accumulator.
    pub rp: f64,
    /// SDFT imaginary accumulator.
    pub ip: f64,
    /// `cos(2π/period)` — rotation phasor real part (constant for a fixed period).
    pub wr: f64,
    /// `sin(2π/period)` — rotation phasor imaginary part.
    pub wi: f64,
}

impl State {
    /// Builds the rotation phasor from `multiplier = 1/period`; `rp` and `ip`
    /// are set to zero and must be initialised by a seeding DFT before use.
    pub fn new(multiplier: f64) -> Self {
        let angle = TPI * multiplier;
        let (wi, wr) = angle.sin_cos();
        Self {
            rp: 0.0,
            ip: 0.0,
            wr,
            wi,
        }
    }
}

/// Compatibility shim for `msw_simd` — provides the `TPI` SIMD constant used
/// in the SIMD by-assets DFT loop and the `Constants` trait surface in `msw_simd`.
pub struct MSWConstants<const N: usize>;
impl<const N: usize> MSWConstants<N> {
    pub const TPI: Simd<f64, N> = Simd::splat(TPI);
}

// ── Cross-call state ──────────────────────────────────────────────────────────

/// Cross-call streaming state for the optimised MSW indicator.
///
/// Holds everything required to resume computation across `batch_indicator` calls:
/// - `state`         — the 4-float SDFT accumulator (see [`State`])
/// - `real`          — the last `period` price bars for the SDFT `old_sample` term
/// - `cos_twiddles` / `sin_twiddles` — precomputed once, cached here so the
///   periodic re-anchor uses pure FMA (no trig) rather than recomputing `sin_cos`
#[derive(Serialize, Deserialize)]
pub struct IndicatorState {
    state: State,
    /// Last `period` input bars — sliding window tail.
    real: Vec<f64>,
    period: usize,
    /// Precomputed cos twiddle factors, length = `period`.
    cos_twiddles: Vec<f64>,
    /// Precomputed sin twiddle factors, length = `period`.
    sin_twiddles: Vec<f64>,
}

impl IndicatorState {
    /// Constructs streaming state from a completed indicator run.
    ///
    /// Used by the SIMD by-assets and by-options paths which compute outputs
    /// via their own drivers and then need a valid [`IndicatorState`] for
    /// subsequent `batch_indicator` calls.
    ///
    /// Computes `(rp, ip)` as the DFT of the last `period` bars — identical
    /// to the final SDFT accumulator value after a full `indicator()` run.
    pub fn new(real: &[f64], period: usize, multiplier: f64) -> Self {
        let (cos_twiddles, sin_twiddles) = precompute_twiddles(period, multiplier);
        let (rp, ip) = match period {
            0..=7 => {
                dot_product_simd::<4>(&real[real.len() - period..], &cos_twiddles, &sin_twiddles)
            }
            _ => dot_product_simd::<8>(&real[real.len() - period..], &cos_twiddles, &sin_twiddles),
        };
        Self {
            state: State {
                rp,
                ip,
                ..State::new(multiplier)
            },
            real: real[real.len() - period..].to_vec(),
            period,
            cos_twiddles,
            sin_twiddles,
        }
    }
}

impl TIndicatorState<1> for IndicatorState {
    fn batch_indicator(
        &mut self,
        inputs: &[&[f64]; INPUTS_WIDTH],
        _optional_outputs: Option<&[bool]>,
    ) -> Result<Vec<Vec<f64>>, IndicatorError> {
        validate_inputs(inputs, 1)?;

        let new_data = inputs[0];
        // self.real = [period old bars | new_data] after the extend.
        self.real.extend_from_slice(new_data);

        let n = new_data.len();
        let (mut sine_line, mut lead_line) =
            (crate::uninit_vec!(f64, n), crate::uninit_vec!(f64, n));

        cycle_sdft(
            &self.real,
            self.period,
            &mut self.state,
            &mut sine_line,
            &mut lead_line,
        );

        // Keep only the last `period` bars for the next batch's old_sample lookups.
        self.real.drain(..self.real.len() - self.period);

        Ok(vec![sine_line, lead_line])
    }
}

// ── Public API ────────────────────────────────────────────────────────────────


/// Returns the minimum number of input bars required.
pub fn min_data(options: &[f64]) -> usize {
    options[0] as usize + 1
}

/// Returns the number of output bars produced from `data_len` input bars.
pub fn output_length(data_len: usize, options: &[f64]) -> usize {
    data_len - min_data(options) + 1
}

/// Returns `1 / period` — the DFT frequency multiplier for a given period.
pub fn multiplier(period: usize) -> f64 {
    1.0 / period as f64
}

/// Calculates the optimised MESA Sine Wave over the full input dataset.
///
/// # Inputs
/// * `inputs[0]` — real (price series, e.g. close)
///
/// # Options
/// * `options[0]` — period
///
/// # Returns
/// `Ok((outputs, state))` where `outputs[0]` = `msw_sine`, `outputs[1]` = `msw_lead`,
/// and `state` can be passed to `IndicatorState::batch_indicator` for streaming.
pub fn indicator(
    inputs: &[&[f64]; INPUTS_WIDTH],
    options: &[f64; OPTIONS_WIDTH],
    _optional_outputs: Option<&[bool]>,
) -> Result<(Vec<Vec<f64>>, IndicatorState), IndicatorError> {
    validate_options(options)?;
    let period = options[0] as usize;
    let mult = multiplier(period);

    validate_inputs(inputs, min_data(options))?;
    let real = inputs[0];

    let capacity = output_length(real.len(), options);
    let mut sine_line = crate::uninit_vec!(f64, capacity);
    let mut lead_line = crate::uninit_vec!(f64, capacity);

    // Twiddles computed once here; stored in IndicatorState for re-use.
    let (cos_twiddles, sin_twiddles) = precompute_twiddles(period, mult);

    let state = match period {
        0..=7 => cycle_msw_sdft::<4>(
            real,
            period,
            mult,
            &cos_twiddles,
            &sin_twiddles,
            &mut sine_line,
            &mut lead_line,
        ),
        _ => cycle_msw_sdft::<8>(
            real,
            period,
            mult,
            &cos_twiddles,
            &sin_twiddles,
            &mut sine_line,
            &mut lead_line,
        ),
    };

    Ok((
        vec![sine_line, lead_line],
        IndicatorState {
            state,
            real: real[real.len() - period..].to_vec(),
            period,
            cos_twiddles,
            sin_twiddles,
        },
    ))
}

// ── Proposal 2: Sliding DFT ───────────────────────────────────────────────────

/// Advances the Sliding DFT by one bar and writes the outputs into `state`.
///
/// This is the per-bar "calc" function, equivalent to `msw::calc` in the
/// original implementation.  Other indicators can call this directly once they
/// have seeded `state` (e.g. via `dot_product_simd` + `State::new`).
///
/// # Arguments
/// * `state`      — SDFT accumulator; updated in-place.
/// * `new_sample` — Price entering the window (newest bar).
/// * `old_sample` — Price leaving the window (oldest bar, `period` steps back).
///
/// # Returns
/// `(sine, lead_sine)` for this bar.
#[inline(always)]
pub fn calc(state: &mut State, new_sample: f64, old_sample: f64) -> (f64, f64) {
    let rp = state.wr.mul_add(state.rp, -(state.wi * state.ip)) + (new_sample - old_sample);
    let ip = state.wr.mul_add(state.ip, state.wi * state.rp);
    state.rp = rp;
    state.ip = ip;
    phase_from_rp_ip(rp, ip)
}

/// Full per-bar DFT returning `(rp, ip)` — for the SIMD by-options path where
/// each lane has a different period and SDFT state cannot be shared.
///
/// Precomputes twiddles on each call (O(period) trig); acceptable since this
/// function is called once per bar per option lane, not in a tight loop.
#[inline(always)]
pub fn calc_rp_ip<const N: usize>(window: &[f64], multiplier: f64) -> (f64, f64) {
    let (cos_twiddles, sin_twiddles) = precompute_twiddles(window.len(), multiplier);
    dot_product_simd::<N>(window, &cos_twiddles, &sin_twiddles)
}

/// Full per-bar DFT + phase — drop-in replacement for the old `msw::calc`.
///
/// Used by `adaptivemsw` where the period changes each bar and the SDFT
/// recurrence cannot be applied. O(period) per bar.
#[inline(always)]
pub fn calc_full<const N: usize>(window: &[f64], multiplier: f64) -> (f64, f64) {
    let (rp, ip) = calc_rp_ip::<N>(window, multiplier);
    phase_from_rp_ip(rp, ip)
}

/// SDFT hot loop — iterates over `real[period..]`, calling `calc` for each bar.
///
/// Used by both the initial full run (`cycle_msw_sdft`) and streaming
/// (`batch_indicator`).  `real` must be laid out as `[period old bars | new bars]`
/// so that `real[i - period]` is the sample leaving the window at step `i`.
#[inline(always)]
fn cycle_sdft(
    real: &[f64],
    period: usize,
    state: &mut State,
    sine_line: &mut [f64],
    lead_line: &mut [f64],
) {
    for i in period..real.len() {
        let j = i - period;
        let (sine, lead) = calc(state, real[i], real[i - period]);
        unsafe {
            *sine_line.get_unchecked_mut(j) = sine;
            *lead_line.get_unchecked_mut(j) = lead;
        }
    }
}

/// Seeds the SDFT from precomputed twiddles and runs the full output loop.
///
/// Returns the final [`State`] so `indicator` can move it into [`IndicatorState`].
fn cycle_msw_sdft<const N: usize>(
    real: &[f64],
    period: usize,
    multiplier: f64,
    cos_twiddles: &[f64],
    sin_twiddles: &[f64],
    sine_line: &mut [f64],
    lead_line: &mut [f64],
) -> State {
    let (rp, ip) = dot_product_simd::<N>(&real[..period], cos_twiddles, sin_twiddles);
    let mut state = State {
        rp,
        ip,
        ..State::new(multiplier)
    };
    cycle_sdft(real, period, &mut state, sine_line, lead_line);
    state
}

// ── Proposal 3: phase helper ──────────────────────────────────────────────────

/// Converts `(rp, ip)` DFT components to `(sine, lead_sine)`.
///
/// Uses a single `sin_cos` call plus the angle-addition identity
/// `sin(x + π/4) = (sin x + cos x) / √2` instead of two separate `sin` calls.
#[inline(always)]
pub fn phase_from_rp_ip(rp: f64, ip: f64) -> (f64, f64) {
    let mut p = if rp.abs() > 0.001 {
        (ip / rp).atan()
    } else {
        PI * if ip < 0.0 { -1.0 } else { 1.0 }
    };
    if rp < 0.0 {
        p += PI;
    }
    p += HPI;
    if p < 0.0 {
        p += TPI;
    }
    if p > TPI {
        p -= TPI;
    }
    let (sp, cp) = p.sin_cos();
    let lead = INV_SQRT2.mul_add(sp, INV_SQRT2 * cp);
    (sp, lead)
}

// ── Proposal 1: precomputed twiddle tables ────────────────────────────────────

/// Precomputes the DFT twiddle factors for a fixed `period` once, O(P).
///
/// Returns flat `(cos_twiddles, sin_twiddles)` of length `period`.
/// Flat `Vec<f64>` is directly serializable and is loaded into SIMD registers
/// on-the-fly inside `dot_product_simd` — no pre-packing needed.
///
/// Angle convention for position `k` in the price window:
/// `angle = 2π·(period−1−k) / period`
pub fn precompute_twiddles(period: usize, multiplier: f64) -> (Vec<f64>, Vec<f64>) {
    let mut cos_twiddles = Vec::with_capacity(period);
    let mut sin_twiddles = Vec::with_capacity(period);
    for k in 0..period {
        let angle = TPI * (period - 1 - k) as f64 * multiplier;
        let (s, c) = angle.sin_cos();
        cos_twiddles.push(c);
        sin_twiddles.push(s);
    }
    (cos_twiddles, sin_twiddles)
}

/// SIMD FMA dot product over a price window using flat twiddle slices.
///
/// Loads N-wide SIMD chunks from the flat `cos_twiddles` / `sin_twiddles` slices
/// directly — no pre-packing required.  Handles the scalar tail cleanly.
#[inline(always)]
pub(crate) fn dot_product_simd<const N: usize>(
    price_window: &[f64],
    cos_twiddles: &[f64],
    sin_twiddles: &[f64],
) -> (f64, f64) {
    let mut rp_acc = Simd::<f64, N>::splat(0.0);
    let mut ip_acc = Simd::<f64, N>::splat(0.0);

    let full_chunks = price_window.len() / N;
    for i in 0..full_chunks {
        let offset = i * N;
        let w = Simd::<f64, N>::from_slice(&price_window[offset..]);
        let c = Simd::<f64, N>::from_slice(&cos_twiddles[offset..]);
        let s = Simd::<f64, N>::from_slice(&sin_twiddles[offset..]);
        rp_acc = c.mul_add(w, rp_acc);
        ip_acc = s.mul_add(w, ip_acc);
    }

    let mut rp = rp_acc.reduce_sum();
    let mut ip = ip_acc.reduce_sum();

    // Scalar tail.
    for i in (full_chunks * N)..price_window.len() {
        rp = cos_twiddles[i].mul_add(price_window[i], rp);
        ip = sin_twiddles[i].mul_add(price_window[i], ip);
    }

    (rp, ip)
}