using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Runtime.Intrinsics; using System.Runtime.Intrinsics.X86; namespace QuanTAlib; /// /// ALMA: Arnaud Legoux Moving Average /// /// /// Gaussian-weighted FIR filter with configurable offset and sigma parameters. /// The offset controls the peak of the Gaussian (0 = leftmost, 1 = rightmost). /// The sigma controls the width of the Gaussian curve. /// /// Calculation: ALMA = Σ(w_i × P_i) / Σ(w_i) where w_i = exp(-((i - m)²) / (2s²)), /// m = offset × (period - 1), s = period / sigma. /// /// Detailed documentation [SkipLocalsInit] public sealed class Alma : AbstractBase { private readonly int _period; private readonly double _offset; private readonly double _sigma; private readonly double[] _weights; private readonly RingBuffer _buffer; private readonly ITValuePublisher? _source; private readonly TValuePublishedHandler? _handler; private bool _disposed; [StructLayout(LayoutKind.Auto)] private record struct State(double LastInput, double LastValidValue, bool HasSeenValidData); private State _state; private State _pState; /// /// Default value to use for LastValidValue when no valid data has been seen yet. /// Defaults to double.NaN to avoid silently introducing zeros. /// public double DefaultLastValidValue { get; set; } = double.NaN; /// /// Initializes a new instance of the class. /// /// The lookback window size. Must be greater than 0. /// The Gaussian peak offset (0.0 to 1.0). Default: 0.85. /// The Gaussian width divisor. Default: 6.0. public Alma(int period, double offset = 0.85, double sigma = 6.0) { if (period <= 0) { throw new ArgumentException("Period must be greater than 0", nameof(period)); } if (offset < 0.0 || offset > 1.0) { throw new ArgumentOutOfRangeException(nameof(offset), offset, "Offset must be between 0.0 and 1.0"); } if (sigma <= 0.0) { throw new ArgumentException("Sigma must be greater than 0", nameof(sigma)); } _period = period; _offset = offset; _sigma = sigma; _buffer = new RingBuffer(period); _weights = ComputeNormalizedWeights(period, offset, sigma); Name = $"Alma({period},{offset:F2},{sigma:F1})"; WarmupPeriod = period; } /// /// Initializes a new instance of the class with a source publisher. /// public Alma(ITValuePublisher source, int period, double offset = 0.85, double sigma = 6.0) : this(period, offset, sigma) { _source = source; _handler = Handle; source.Pub += _handler; } protected override void Dispose(bool disposing) { if (!_disposed) { if (disposing && _source != null && _handler != null) { _source.Pub -= _handler; } _disposed = true; } base.Dispose(disposing); } public override bool IsHot => _buffer.IsFull; public bool IsNew { get; private set; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private double GetValidValue(double input) { if (double.IsFinite(input)) { _state.LastValidValue = input; _state.HasSeenValidData = true; return input; } return _state.HasSeenValidData ? _state.LastValidValue : DefaultLastValidValue; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private double CalculateWeightedSum() { double result = 0; int count = _buffer.Count; int weightOffset = _period - count; int idx = 0; foreach (double item in _buffer) { result = Math.FusedMultiplyAdd(_weights[weightOffset + idx], item, result); idx++; } // Normalize for partial windows if (count < _period) { double wSum = 0; for (int i = weightOffset; i < _period; i++) { wSum += _weights[i]; } return wSum > 0 ? result / wSum : result; } return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public override TValue Update(TValue input, bool isNew = true) { IsNew = isNew; if (isNew) { double val = GetValidValue(input.Value); _buffer.Add(val); _state.LastInput = val; _pState = _state; } else { if (_buffer.Count == 0) { throw new InvalidOperationException( "Cannot call Update with isNew=false when buffer is empty. " + "The first update must have isNew=true to initialize state."); } _state = _pState; double val = GetValidValue(input.Value); _buffer.UpdateNewest(val); } double result = CalculateWeightedSum(); Last = new TValue(input.Time, result); PubEvent(Last, isNew); return Last; } public override TSeries Update(TSeries source) { if (source.Count == 0) { return []; } int len = source.Count; var t = new List(len); var v = new List(len); CollectionsMarshal.SetCount(t, len); CollectionsMarshal.SetCount(v, len); var tSpan = CollectionsMarshal.AsSpan(t); var vSpan = CollectionsMarshal.AsSpan(v); Batch(source.Values, vSpan, _period, _offset, _sigma); source.Times.CopyTo(tSpan); Prime(source.Values); Last = new TValue(tSpan[len - 1], vSpan[len - 1]); return new TSeries(t, v); } private void Handle(object? sender, in TValueEventArgs e) => Update(e.Value, e.IsNew); public override void Prime(ReadOnlySpan source, TimeSpan? step = null) { if (source.Length == 0) { return; } int len = source.Length; int windowSize = Math.Min(len, _period); int startIndex = len - windowSize; // Seed LastValidValue _state.LastValidValue = DefaultLastValidValue; _state.HasSeenValidData = false; if (startIndex > 0) { for (int i = startIndex - 1; i >= 0; i--) { if (double.IsFinite(source[i])) { _state.LastValidValue = source[i]; _state.HasSeenValidData = true; break; } } } // Reset buffer and process window _buffer.Clear(); for (int i = startIndex; i < len; i++) { double val = GetValidValue(source[i]); _buffer.Add(val); _state.LastInput = val; } // Calculate Last double result = CalculateWeightedSum(); Last = new TValue(DateTime.MinValue, result); _pState = _state; } public override void Reset() { _buffer.Clear(); _state = default; _pState = default; Last = default; } /// /// Computes ALMA for a TSeries using batch processing. /// public static TSeries Batch(TSeries source, int period, double offset = 0.85, double sigma = 6.0) { var alma = new Alma(period, offset, sigma); return alma.Update(source); } /// /// Computes ALMA for raw spans using batch processing. /// [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Batch(ReadOnlySpan source, Span output, int period, double offset = 0.85, double sigma = 6.0) { if (period <= 0) { throw new ArgumentException("Period must be greater than 0", nameof(period)); } if (sigma <= 0.0) { throw new ArgumentException("Sigma must be greater than 0", nameof(sigma)); } if (offset < 0.0 || offset > 1.0) { throw new ArgumentOutOfRangeException(nameof(offset), offset, "Offset must be between 0.0 and 1.0"); } if (source.Length != output.Length) { throw new ArgumentException("Source and output must have the same length", nameof(output)); } int len = source.Length; if (len == 0) { return; } CalculateScalarCore(source, output, period, offset, sigma); } /// /// Computes ALMA and returns both the result series and a warmed-up indicator instance. /// public static (TSeries Results, Alma Indicator) Calculate(TSeries source, int period, double offset = 0.85, double sigma = 6.0) { var indicator = new Alma(period, offset, sigma); TSeries results = indicator.Update(source); return (results, indicator); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void CalculateScalarCore(ReadOnlySpan source, Span output, int period, double offset, double sigma) { int len = source.Length; double[] weights = ComputeNormalizedWeights(period, offset, sigma); double lastValid = double.NaN; Span buffer = period <= 512 ? stackalloc double[period] : new double[period]; int bufferCount = 0; int bufferIdx = 0; for (int i = 0; i < len; i++) { double val = source[i]; if (double.IsFinite(val)) { lastValid = val; } else { val = lastValid; } // Add to circular buffer buffer[bufferIdx] = val; bufferIdx++; if (bufferIdx >= period) { bufferIdx = 0; } if (bufferCount < period) { bufferCount++; } // Compute weighted sum double result = 0; int weightOffset = period - bufferCount; if (bufferCount == period) { // Full window — iterate from oldest to newest int readIdx = bufferIdx; // bufferIdx now points to oldest for (int k = 0; k < period; k++) { result = Math.FusedMultiplyAdd(weights[k], buffer[readIdx], result); readIdx++; if (readIdx >= period) { readIdx = 0; } } } else { // Partial window — use tail weights double wSum = 0; for (int k = 0; k < bufferCount; k++) { int wi = weightOffset + k; result = Math.FusedMultiplyAdd(weights[wi], buffer[k], result); wSum += weights[wi]; } result = wSum > 0 ? result / wSum : result; } output[i] = result; } } /// /// Pre-computes normalized Gaussian weights for the ALMA filter. /// Weights are normalized so that their sum equals 1.0. /// [MethodImpl(MethodImplOptions.AggressiveInlining)] private static double[] ComputeNormalizedWeights(int period, double offset, double sigma) { double[] w = new double[period]; double m = offset * (period - 1); double s = period / sigma; double s2 = 2.0 * s * s; double wSum = 0; for (int i = 0; i < period; i++) { double d = i - m; w[i] = Math.Exp(-(d * d) / s2); wSum += w[i]; } // Normalize weights to sum to 1.0 double invSum = 1.0 / wSum; for (int i = 0; i < period; i++) { w[i] *= invSum; } return w; } }