QuanTAlib/lib/channels/bbands/Bbands.cs

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

namespace QuanTAlib;

/// <summary>
/// BBANDS: Bollinger Bands
/// A volatility-based channel indicator consisting of a moving average middle band
/// with upper and lower bands positioned at a specified number of standard deviations.
/// Developed by John Bollinger in the 1980s for volatility analysis.
/// </summary>
/// <remarks>
/// The BBands calculation process:
/// 1. Calculate SMA of price over the period
/// 2. Calculate standard deviation over the period
/// 3. Upper band = SMA + (multiplier × StdDev)
/// 4. Lower band = SMA - (multiplier × StdDev)
///
/// Key characteristics:
/// - Adapts dynamically to volatility changes
/// - Wider bands indicate higher volatility
/// - Narrower bands indicate lower volatility
/// - Price tends to oscillate between bands
/// - Can identify overbought/oversold conditions
///
/// Sources:
///     John Bollinger - "Bollinger on Bollinger Bands" (2001)
///     https://www.bollingerbands.com/
/// </remarks>
[SkipLocalsInit]
public sealed class Bbands : AbstractBase
{
    private readonly Sma _sma;
    private readonly StdDev _stdev;
    private readonly int _period;
    private readonly double _multiplier;
    private const int DefaultPeriod = 20;
    private const double DefaultMultiplier = 2.0;
    private const double MinMultiplier = 0.1;
    private const int MinPeriod = 2;

    public override bool IsHot => _index >= WarmupPeriod;
    private int _index;

    /// <summary>
    /// Middle band (SMA of price)
    /// </summary>
    public TValue Middle { get; private set; }

    /// <summary>
    /// Upper band (SMA + multiplier × StdDev)
    /// </summary>
    public TValue Upper { get; private set; }

    /// <summary>
    /// Lower band (SMA - multiplier × StdDev)
    /// </summary>
    public TValue Lower { get; private set; }

    /// <summary>
    /// Band width (Upper - Lower)
    /// </summary>
    public TValue Width { get; private set; }

    /// <summary>
    /// Percent B: (Price - Lower) / (Upper - Lower)
    /// </summary>
    public TValue PercentB { get; private set; }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public Bbands(int period = DefaultPeriod, double multiplier = DefaultMultiplier)
    {
        if (period < MinPeriod)
        {
            throw new ArgumentOutOfRangeException(nameof(period),
                $"Period must be at least {MinPeriod}.");
        }
        if (multiplier < MinMultiplier)
        {
            throw new ArgumentOutOfRangeException(nameof(multiplier),
                $"Multiplier must be at least {MinMultiplier}.");
        }

        _period = period;
        _multiplier = multiplier;
        _sma = new Sma(period);
        _stdev = new StdDev(period, isPopulation: true);
        WarmupPeriod = period;
        Name = $"Bbands({period},{multiplier:F1})";
        Init();
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public Bbands(ITValuePublisher source, int period = DefaultPeriod, double multiplier = DefaultMultiplier)
        : this(period, multiplier)
    {
        source.Pub += Handle;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private void Init()
    {
        _index = 0;
        Middle = new TValue(DateTime.UtcNow, 0);
        Upper = new TValue(DateTime.UtcNow, 0);
        Lower = new TValue(DateTime.UtcNow, 0);
        Width = new TValue(DateTime.UtcNow, 0);
        PercentB = new TValue(DateTime.UtcNow, 0);
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
    private static double GetFiniteValue(double value, double fallback) =>
        double.IsFinite(value) ? value : fallback;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public override TValue Update(TValue input, bool isNew = true)
    {
        if (isNew)
        {
            _index++;
        }

        double finiteValue = GetFiniteValue(input.Value, Middle.Value);

        // Update SMA and StdDev
        TValue smaValue = _sma.Update(new TValue(input.Time, finiteValue), isNew);
        TValue stdevValue = _stdev.Update(new TValue(input.Time, finiteValue), isNew);

        double middle = smaValue.Value;
        double stdDev = stdevValue.Value;
        double offset = _multiplier * stdDev;

        double upper = middle + offset;
        double lower = middle - offset;
        double width = upper - lower;

        // Calculate Percent B
        double percentB = 0.0;
        if (width > double.Epsilon)
        {
            percentB = (finiteValue - lower) / width;
        }

        // Update all band values
        Middle = new TValue(input.Time, middle);
        Upper = new TValue(input.Time, upper);
        Lower = new TValue(input.Time, lower);
        Width = new TValue(input.Time, width);
        PercentB = new TValue(input.Time, percentB);
        Last = Middle;

        PubEvent(Middle, isNew);
        return Middle;
    }

    /// <summary>
    /// Updates the indicator with a new time series and returns the middle band series.
    /// </summary>
    public override TSeries Update(TSeries source)
    {
        if (source == null)
        {
            throw new ArgumentNullException(nameof(source));
        }

        ReadOnlySpan<double> sourceSpan = source.Values;
        ReadOnlySpan<long> timeSpan = source.Times;
        int len = sourceSpan.Length;

        TSeries middleSeries = new(capacity: len);

        // Use ArrayPool to avoid stack overflow for large series
        double[] middleRented = ArrayPool<double>.Shared.Rent(len);
        double[] upperRented = ArrayPool<double>.Shared.Rent(len);
        double[] lowerRented = ArrayPool<double>.Shared.Rent(len);

        try
        {
            Span<double> middleSpan = middleRented.AsSpan(0, len);
            Span<double> upperSpan = upperRented.AsSpan(0, len);
            Span<double> lowerSpan = lowerRented.AsSpan(0, len);

            Batch(sourceSpan, middleSpan, upperSpan, lowerSpan, _period, _multiplier);

            for (int i = 0; i < len; i++)
            {
                middleSeries.Add(timeSpan[i], middleSpan[i], isNew: true);
            }
        }
        finally
        {
            ArrayPool<double>.Shared.Return(middleRented);
            ArrayPool<double>.Shared.Return(upperRented);
            ArrayPool<double>.Shared.Return(lowerRented);
        }

        // Restore state from the last period values
        Reset();
        int startIdx = Math.Max(0, len - _period);
        for (int i = startIdx; i < len; i++)
        {
            Update(new TValue(timeSpan[i], sourceSpan[i]), isNew: true);
        }

        return middleSeries;
    }

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

    public override void Reset()
    {
        _sma.Reset();
        _stdev.Reset();
        Init();
    }

    public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
    {
        step ??= TimeSpan.FromSeconds(1);
        DateTime startTime = DateTime.UtcNow;

        for (int i = 0; i < source.Length; i++)
        {
            Update(new TValue(startTime + i * step.Value, source[i]), isNew: true);
        }
    }

    /// <summary>
    /// Calculates Bollinger Bands for the entire series and returns the middle band series.
    /// </summary>
    public static TSeries Batch(TSeries source, int period = DefaultPeriod, double multiplier = DefaultMultiplier)
    {
        Bbands bbands = new(period, multiplier);
        return bbands.Update(source);
    }

    /// <summary>
    /// Calculates Bollinger Bands across all input values using SIMD-optimized operations where possible.
    /// </summary>
    public static void Batch(
        ReadOnlySpan<double> source,
        Span<double> middle,
        Span<double> upper,
        Span<double> lower,
        int period = DefaultPeriod,
        double multiplier = DefaultMultiplier)
    {
        if (source.Length != middle.Length || source.Length != upper.Length || source.Length != lower.Length)
        {
            throw new ArgumentException("All spans must have the same length.", nameof(source));
        }
        if (period < MinPeriod)
        {
            throw new ArgumentOutOfRangeException(nameof(period),
                $"Period must be at least {MinPeriod}.");
        }
        if (multiplier < MinMultiplier)
        {
            throw new ArgumentOutOfRangeException(nameof(multiplier),
                $"Multiplier must be at least {MinMultiplier}.");
        }

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

        // Calculate SMA using static batch method
        Sma.Batch(source, middle, period);

        // Calculate standard deviation and bands using O(n) rolling sums
        // Instead of O(n²) nested loop, maintain running sum and sumSq
        // Track count of finite values to properly compute mean/variance
        double rollingSum = 0.0;
        double rollingSumSq = 0.0;
        int finiteCount = 0;

        // Initialize rolling sums for first window
        for (int i = 0; i < Math.Min(period, len); i++)
        {
            double val = source[i];
            if (double.IsFinite(val))
            {
                rollingSum += val;
                rollingSumSq += val * val;
                finiteCount++;
            }

            if (i < period - 1)
            {
                upper[i] = double.NaN;
                lower[i] = double.NaN;
            }
        }

        // Process first complete window
        if (len >= period)
        {
            if (finiteCount == period)
            {
                double mean = rollingSum / finiteCount;
                double variance = (rollingSumSq / finiteCount) - (mean * mean);
                variance = Math.Max(0.0, variance); // Guard against negative due to floating point
                double stdDev = Math.Sqrt(variance);
                double offset = multiplier * stdDev;
                upper[period - 1] = middle[period - 1] + offset;
                lower[period - 1] = middle[period - 1] - offset;
            }
            else
            {
                // Not all values in window are finite, emit NaN
                upper[period - 1] = double.NaN;
                lower[period - 1] = double.NaN;
            }
        }

        // Process remaining bars with O(1) rolling update
        for (int i = period; i < len; i++)
        {
            // Remove outgoing value (leftmost of previous window)
            double outgoing = source[i - period];
            if (double.IsFinite(outgoing))
            {
                rollingSum -= outgoing;
                rollingSumSq -= outgoing * outgoing;
                finiteCount--;
            }

            // Add incoming value (current)
            double incoming = source[i];
            if (double.IsFinite(incoming))
            {
                rollingSum += incoming;
                rollingSumSq += incoming * incoming;
                finiteCount++;
            }

            // Only compute bands when all values in window are finite
            if (finiteCount == period)
            {
                // Calculate variance from rolling sums: Var = E[X²] - E[X]²
                double mean = rollingSum / finiteCount;
                double variance = (rollingSumSq / finiteCount) - (mean * mean);
                variance = Math.Max(0.0, variance); // Guard against negative due to floating point

                double stdDev = Math.Sqrt(variance);
                double offset = multiplier * stdDev;

                upper[i] = middle[i] + offset;
                lower[i] = middle[i] - offset;
            }
            else
            {
                // Window contains non-finite values, emit NaN
                upper[i] = double.NaN;
                lower[i] = double.NaN;
            }
        }
    }

    public static (TSeries Results, Bbands Indicator) Calculate(TSeries source, int period = DefaultPeriod, double multiplier = DefaultMultiplier)
    {
        var indicator = new Bbands(period, multiplier);
        TSeries results = indicator.Update(source);
        return (results, indicator);
    }

}