QuanTAlib/lib/trends_IIR/ema/Ema.cs

using System.Buffers;
using System.Diagnostics.Contracts;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;

namespace QuanTAlib;

/// <summary>
/// EMA: Exponential Moving Average
/// </summary>
/// <remarks>
/// Applies exponentially decreasing weights to give more importance to recent values.
/// Faster response to price changes than SMA; commonly used for trend identification.
///
/// Calculation: <c>EMA_t = α × Price_t + (1-α) × EMA_{t-1}</c>, where <c>α = 2/(period+1)</c>.
/// </remarks>
/// <seealso href="Ema.md">Detailed documentation</seealso>
/// <seealso href="ema.pine">Reference Pine Script implementation</seealso>
[SkipLocalsInit]
public sealed class Ema : AbstractBase
{
    [StructLayout(LayoutKind.Auto)]
    private record struct State(double Ema, double E, bool IsHot, bool IsCompensated)
    {
        public static State New() => new() { Ema = 0, E = 1.0, IsHot = false, IsCompensated = false };
    }

    private readonly double _alpha;
    private readonly double _decay;
    private State _state = State.New();
    private State _p_state = State.New();
    private double _lastValidValue;
    private double _p_lastValidValue;

    /// <summary>
    /// Creates EMA with specified period.
    /// Alpha = 2 / (period + 1)
    /// </summary>
    /// <param name="period">Period for EMA calculation (must be > 0)</param>
    public Ema(int period)
    {
        ArgumentOutOfRangeException.ThrowIfNegativeOrZero(period);

        _alpha = 2.0 / (period + 1);
        _decay = 1.0 - _alpha;
        Name = $"Ema({period})";
        WarmupPeriod = period;
    }

    /// <summary>
    /// Creates EMA with specified source and period.
    /// Subscribes to source.Pub event.
    /// </summary>
    /// <param name="source">Source to subscribe to</param>
    /// <param name="period">Period for EMA calculation</param>
    public Ema(ITValuePublisher source, int period) : this(period)
    {
        source.Pub += Handle;
    }

    public Ema(TSeries source, int period) : this(period)
    {
        Prime(source.Values);
        if (source.Count > 0)
        {
            Last = new TValue(source.LastTime, Last.Value);
        }
        source.Pub += Handle;
    }

    /// <summary>
    /// Creates EMA with specified alpha smoothing factor.
    /// </summary>
    /// <param name="alpha">Smoothing factor (0 < alpha <= 1)</param>
    public Ema(double alpha)
    {
        if (alpha <= 0 || alpha > 1)
        {
            throw new ArgumentException("Alpha must be greater than 0 and at most 1", nameof(alpha));
        }

        _alpha = alpha;
        _decay = 1.0 - alpha;
        Name = $"Ema(α={alpha:F4})";
        // Approximate period from alpha: alpha = 2/(N+1) => N = 2/alpha - 1
        WarmupPeriod = (int)(2.0 / alpha - 1.0);
    }

    /// <summary>
    /// True if the EMA has warmed up and is providing valid results.
    /// </summary>
    public override bool IsHot => _state.IsHot;

    /// <summary>
    /// Maximum size for stackalloc buffer in Prime().
    /// Larger datasets use ArrayPool to avoid stack overflow.
    /// </summary>
    private const int StackAllocThreshold = 512;

    /// <summary>
    /// Initializes the indicator state using the provided history.
    /// Reuses CalculateCore with a temporary buffer to avoid code duplication.
    /// </summary>
    /// <param name="source">Historical data</param>
    public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
    {
        if (source.Length == 0)
        {
            return;
        }

        // Reset state
        _state = State.New();
        _p_state = State.New();
        _lastValidValue = 0;
        _p_lastValidValue = 0;

        int len = source.Length;

        // Find first valid value to seed lastValid
        bool foundValid = false;
        for (int k = 0; k < len; k++)
        {
            if (double.IsFinite(source[k]))
            {
                _lastValidValue = source[k];
                foundValid = true;
                break;
            }
        }

        if (!foundValid)
        {
            Last = new TValue(DateTime.MinValue, double.NaN);
            _p_state = _state;
            _p_lastValidValue = _lastValidValue;
            return;
        }

        // Use temporary buffer to run CalculateCore and extract final state
        // We only care about the state, not the output values
        double[]? rented = len > StackAllocThreshold ? ArrayPool<double>.Shared.Rent(len) : null;
        Span<double> tempOutput = rented != null
            ? rented.AsSpan(0, len)
            : stackalloc double[len];

        try
        {
            CalculateCore(source, tempOutput, _alpha, ref _state, ref _lastValidValue);

            // Extract the final result from the output
            double result = tempOutput[len - 1];
            Last = new TValue(DateTime.MinValue, result);

            // Backup state for the next update cycle
            _p_state = _state;
            _p_lastValidValue = _lastValidValue;
        }
        finally
        {
            if (rented != null)
            {
                ArrayPool<double>.Shared.Return(rented);
            }
        }
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private void Handle(object? sender, in TValueEventArgs e) => Update(e.Value, e.IsNew);

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private double GetValidValue(double input)
    {
        if (double.IsFinite(input))
        {
            _lastValidValue = input;
            return input;
        }
        return _lastValidValue;
    }

    private const double COVERAGE_THRESHOLD = 0.05;
    private const double COMPENSATOR_THRESHOLD = 1e-10;

    [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
    public override TValue Update(TValue input, bool isNew = true)
    {
        if (isNew)
        {
            _p_state = _state;
            _p_lastValidValue = _lastValidValue;
        }
        else
        {
            _state = _p_state;
            _lastValidValue = _p_lastValidValue;
        }

        double val = GetValidValue(input.Value);
        val = Compute(val, _alpha, _decay, ref _state);
        Last = new TValue(input.Time, val);
        PubEvent(Last, isNew);
        return Last;
    }

    [MethodImpl(MethodImplOptions.AggressiveOptimization)]
    public override TSeries Update(TSeries source)
    {
        if (source.Count == 0)
        {
            return [];
        }

        int len = source.Count;
        var t = new List<long>(len);
        var v = new List<double>(len);
        CollectionsMarshal.SetCount(t, len);
        CollectionsMarshal.SetCount(v, len);

        var tSpan = CollectionsMarshal.AsSpan(t);
        var vSpan = CollectionsMarshal.AsSpan(v);
        var sourceValues = source.Values;
        var sourceTimes = source.Times;

        State state = _state;
        double lastValidValue = _lastValidValue;

        CalculateCore(sourceValues, vSpan, _alpha, ref state, ref lastValidValue);

        _state = state;
        _lastValidValue = lastValidValue;

        sourceTimes.CopyTo(tSpan);

        _p_state = _state;
        _p_lastValidValue = _lastValidValue;
        Last = new TValue(tSpan[len - 1], vSpan[len - 1]);

        return new TSeries(t, v);
    }

    /// <summary>
    /// Core EMA computation with bias compensation.
    /// Computes the next EMA value and updates state in place.
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
    private static double Compute(double input, double alpha, double decay, ref State state)
    {
        // EMA update using FMA for precision:
        // state.Ema = state.Ema * decay + alpha * input
        state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * input);

        double result;
        if (!state.IsCompensated)
        {
            state.E *= decay;

            if (!state.IsHot && state.E <= COVERAGE_THRESHOLD)
            {
                state.IsHot = true;
            }

            if (state.E <= COMPENSATOR_THRESHOLD)
            {
                state.IsCompensated = true;
                result = state.Ema;
            }
            else
            {
                result = state.Ema / (1.0 - state.E);
            }
        }
        else
        {
            result = state.Ema;
        }

        return result;
    }

    /// <summary>
    /// Core EMA calculation with bias compensation and NaN handling.
    /// Uses FMA for precision. IIR filters are inherently self-correcting.
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveOptimization)]
    private static void CalculateCore(ReadOnlySpan<double> source, Span<double> output, double alpha, ref State state, ref double lastValidValue)
    {
        int len = source.Length;
        double decay = 1.0 - alpha;
        int i = 0;

        // Phase 1: Compensation phase (before warmup complete)
        if (!state.IsCompensated)
        {
            for (; i < len && state.E > COMPENSATOR_THRESHOLD; i++)
            {
                double val = source[i];
                if (double.IsFinite(val))
                {
                    lastValidValue = val;
                }
                else
                {
                    val = lastValidValue;
                }

                state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * val);
                state.E *= decay;

                if (!state.IsHot && state.E <= COVERAGE_THRESHOLD)
                {
                    state.IsHot = true;
                }

                output[i] = state.Ema / (1.0 - state.E);
            }
            if (state.E <= COMPENSATOR_THRESHOLD)
            {
                state.IsCompensated = true;
            }
        }

        // Phase 2: Post-compensation (hot path) - optimized with loop unrolling
        // Since EMA is inherently serial (each output depends on previous),
        // we optimize by minimizing branching and using FMA
        ref double srcRef = ref MemoryMarshal.GetReference(source);
        ref double outRef = ref MemoryMarshal.GetReference(output);

        // Unroll by 4 to reduce loop overhead and improve instruction-level parallelism
        int unrollEnd = i + ((len - i) / 4) * 4;
        for (; i < unrollEnd; i += 4)
        {
            double v0 = Unsafe.Add(ref srcRef, i);
            if (!double.IsFinite(v0))
            {
                v0 = lastValidValue;
            }
            else
            {
                lastValidValue = v0;
            }

            state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * v0);
            Unsafe.Add(ref outRef, i) = state.Ema;

            double v1 = Unsafe.Add(ref srcRef, i + 1);
            if (!double.IsFinite(v1))
            {
                v1 = lastValidValue;
            }
            else
            {
                lastValidValue = v1;
            }

            state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * v1);
            Unsafe.Add(ref outRef, i + 1) = state.Ema;

            double v2 = Unsafe.Add(ref srcRef, i + 2);
            if (!double.IsFinite(v2))
            {
                v2 = lastValidValue;
            }
            else
            {
                lastValidValue = v2;
            }

            state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * v2);
            Unsafe.Add(ref outRef, i + 2) = state.Ema;

            double v3 = Unsafe.Add(ref srcRef, i + 3);
            if (!double.IsFinite(v3))
            {
                v3 = lastValidValue;
            }
            else
            {
                lastValidValue = v3;
            }

            state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * v3);
            Unsafe.Add(ref outRef, i + 3) = state.Ema;

        }

        // Scalar remainder
        for (; i < len; i++)
        {
            double val = Unsafe.Add(ref srcRef, i);
            if (!double.IsFinite(val))
            {
                val = lastValidValue;
            }
            else
            {
                lastValidValue = val;
            }

            state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * val);
            Unsafe.Add(ref outRef, i) = state.Ema;
        }
    }

    /// <summary>
    /// SIMD-optimized EMA calculation for large, clean (NaN-free) datasets.
    /// Since EMA is inherently serial, this method uses SIMD for the input preprocessing
    /// and optimized scalar computation with loop unrolling.
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveOptimization)]
    private static void CalculateCleanCore(ReadOnlySpan<double> source, Span<double> output, double alpha)
    {
        int len = source.Length;
        double decay = 1.0 - alpha;

        ref double srcRef = ref MemoryMarshal.GetReference(source);
        ref double outRef = ref MemoryMarshal.GetReference(output);

        // Initialize with first value (no compensation since we assume clean data)
        double ema = Unsafe.Add(ref srcRef, 0);
        Unsafe.Add(ref outRef, 0) = ema;

        // Pre-multiply alpha for efficiency
        double alphaVal;

        // Unroll by 4 for better ILP
        int i = 1;
        int unrollEnd = 1 + ((len - 1) / 4) * 4;

        for (; i < unrollEnd; i += 4)
        {
            alphaVal = alpha * Unsafe.Add(ref srcRef, i);
            ema = Math.FusedMultiplyAdd(ema, decay, alphaVal);
            Unsafe.Add(ref outRef, i) = ema;

            alphaVal = alpha * Unsafe.Add(ref srcRef, i + 1);
            ema = Math.FusedMultiplyAdd(ema, decay, alphaVal);
            Unsafe.Add(ref outRef, i + 1) = ema;

            alphaVal = alpha * Unsafe.Add(ref srcRef, i + 2);
            ema = Math.FusedMultiplyAdd(ema, decay, alphaVal);
            Unsafe.Add(ref outRef, i + 2) = ema;

            alphaVal = alpha * Unsafe.Add(ref srcRef, i + 3);
            ema = Math.FusedMultiplyAdd(ema, decay, alphaVal);
            Unsafe.Add(ref outRef, i + 3) = ema;
        }

        // Scalar remainder
        for (; i < len; i++)
        {
            alphaVal = alpha * Unsafe.Add(ref srcRef, i);
            ema = Math.FusedMultiplyAdd(ema, decay, alphaVal);
            Unsafe.Add(ref outRef, i) = ema;
        }
    }

    /// <summary>
    /// Minimum dataset size to use optimized clean path.
    /// Below this threshold, the overhead of checking for NaN isn't worth it.
    /// </summary>
    private const int CleanPathThreshold = 256;

    /// <summary>
    /// Runs a high-performance batch calculation on history and returns
    /// a "Hot" Ema instance ready to process the next tick immediately.
    /// </summary>
    /// <param name="source">Historical time series</param>
    /// <param name="period">EMA Period</param>
    /// <returns>A tuple containing the full calculation results and the hot indicator instance</returns>

    /// <summary>
    /// Calculates EMA for the entire series using a new instance.
    /// </summary>
    /// <param name="source">Input series</param>
    /// <param name="period">EMA period</param>
    /// <returns>EMA series</returns>
    public static TSeries Batch(TSeries source, int period)
    {
        var ema = new Ema(period);
        return ema.Update(source);
    }

    /// <summary>
    /// Calculates EMA in-place using period, writing results to pre-allocated output span.
    /// Zero-allocation method for maximum performance.
    /// Alpha = 2 / (period + 1)
    /// </summary>
    /// <param name="source">Input values</param>
    /// <param name="output">Output span (must be same length as source)</param>
    /// <param name="period">EMA period (must be > 0)</param>
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static void Batch(ReadOnlySpan<double> source, Span<double> output, int period)
    {
        if (period <= 0)
        {
            throw new ArgumentException("Period must be greater than 0", nameof(period));
        }

        double alpha = 2.0 / (period + 1);
        Batch(source, output, alpha);
    }

    /// <summary>
    /// Calculates EMA in-place using alpha, writing results to pre-allocated output span.
    /// Automatically uses optimized path for large, NaN-free datasets.
    /// </summary>
    /// <param name="source">Input values</param>
    /// <param name="output">Output span (must be same length as source)</param>
    /// <param name="alpha">Smoothing factor (0 < alpha <= 1)</param>
    [MethodImpl(MethodImplOptions.AggressiveOptimization)]
    public static void Batch(ReadOnlySpan<double> source, Span<double> output, double alpha)
    {
        if (source.Length != output.Length)
        {
            throw new ArgumentException("Source and output must have the same length", nameof(output));
        }

        ArgumentOutOfRangeException.ThrowIfLessThanOrEqual(alpha, 0.0);
        ArgumentOutOfRangeException.ThrowIfGreaterThan(alpha, 1.0);

        if (source.Length == 0)
        {
            return;
        }

        // For large, clean datasets, use optimized path without NaN handling
        if (source.Length >= CleanPathThreshold && !source.ContainsNonFinite())
        {
            CalculateCleanCore(source, output, alpha);
            return;
        }

        // Standard path with NaN handling
        var state = State.New();
        double lastValid = 0;
        bool foundValid = false;

        // Find first valid value to seed lastValid
        for (int k = 0; k < source.Length; k++)
        {
            if (double.IsFinite(source[k]))
            {
                lastValid = source[k];
                foundValid = true;
                break;
            }
        }

        if (!foundValid)
        {
            output.Fill(double.NaN);
            return;
        }

        CalculateCore(source, output, alpha, ref state, ref lastValid);
    }
    public static (TSeries Results, Ema Indicator) Calculate(TSeries source, int period)
    {
        var ema = new Ema(period);
        TSeries results = ema.Update(source);
        return (results, ema);
    }

    /// <summary>
    /// Resets the EMA state.
    /// </summary>
    public override void Reset()
    {
        _state = State.New();
        _p_state = _state;
        _lastValidValue = 0;
        _p_lastValidValue = 0;
        Last = default;
    }
}