diff --git a/.clinerules/AGENTS.md b/.clinerules/AGENTS.md index 2f9599df..87b4b328 100644 --- a/.clinerules/AGENTS.md +++ b/.clinerules/AGENTS.md @@ -44,6 +44,43 @@ We do not store objects in lists. We store primitive arrays. 3. **SIMD**: Batch operations (`Calculate`) should use `System.Runtime.Intrinsics` (AVX2) or `System.Numerics.Vector` where possible. Use `Vector.ConditionalSelect` to handle edge cases (e.g., division by zero) without branching. If SIMD is not possible due to recursive dependencies, use `stackalloc` for internal buffers to avoid heap allocations. 4. **Inlining**: Use `[MethodImpl(MethodImplOptions.AggressiveInlining)]` on hot methods. 5. **Locals**: Use `[SkipLocalsInit]` to avoid zero-init costs in tight loops. +6. **Fused Multiply-Add (FMA)**: Use `Math.FusedMultiplyAdd(a, b, c)` for `a*b+c` operations. FMA performs the operation with a single rounding step, improving both precision and potentially performance. + +### FMA Patterns + +Use `Math.FusedMultiplyAdd()` for these common patterns: + +| Pattern | Before | After | +| ------- | ------ | ----- | +| **EMA Smoothing** | `x + alpha * (y - x)` | `Math.FusedMultiplyAdd(x, decay, alpha * y)` where `decay = 1 - alpha` | +| **Weighted Sum** | `a * w1 + b * w2` | `Math.FusedMultiplyAdd(a, w1, b * w2)` | +| **Linear Combo** | `3.0 * a - b` | `Math.FusedMultiplyAdd(3.0, a, -b)` | +| **Cross Product** | `(a * b) + (c * d)` | `Math.FusedMultiplyAdd(a, b, c * d)` | +| **IIR Filter** | `coef * input + feedback * state` | `Math.FusedMultiplyAdd(coef, input, feedback * state)` | + +**When to use FMA:** + +* EMA-style smoothing operations (most moving averages) +* IIR filter calculations (Butterworth, Chebyshev, SSF) +* Homodyne discriminator calculations (HTIT, MAMA) +* Any `a*b+c` pattern in hot paths + +**When NOT to use FMA:** + +* Simple additions or multiplications (no benefit) +* When intermediate rounding is mathematically required +* In SIMD paths (use `Fma.MultiplyAdd`, `Avx512F.FusedMultiplyAdd`, or `AdvSimd.Arm64.FusedMultiplyAdd` instead) + +**Pre-compute decay constants:** + +```csharp +// In constructor or field initialization +private readonly double _alpha; +private readonly double _decay; // = 1 - _alpha + +// In hot path +result = Math.FusedMultiplyAdd(prevState, _decay, _alpha * newInput); +``` ## 3. Indicator Implementation Standards diff --git a/docs/_sidebar.md b/docs/_sidebar.md index 48e87573..b1edf560 100644 --- a/docs/_sidebar.md +++ b/docs/_sidebar.md @@ -69,6 +69,7 @@ - **Statistics** - [Overview](../lib/statistics/_index.md) + - [CMA - Cumulative MA](../lib/statistics/cma/Cma.md) - [COVARIANCE - Covariance](../lib/statistics/covariance/Covariance.md) - [LINREG - Linear Regression Curve](../lib/statistics/linreg/LinReg.md) - [MEDIAN - Rolling Median](../lib/statistics/median/Median.md) @@ -85,4 +86,4 @@ - [Overview](../lib/forecasts/_index.md) - **Cycles** - - [Overview](../lib/cycles/_index.md) + - [Overview](../lib/cycles/_index.md) \ No newline at end of file diff --git a/docs/indicators.md b/docs/indicators.md index 272f4411..74e60222 100644 --- a/docs/indicators.md +++ b/docs/indicators.md @@ -104,9 +104,10 @@ These measure the spread of data points around the mean. ### Statistics +- [**CMA**](../lib/statistics/cma/Cma.md) - Cumulative Moving Average - [**COVARIANCE**](../lib/statistics/covariance/Covariance.md) - Covariance - [**LINREG**](../lib/statistics/linreg/LinReg.md) - Linear Regression Curve - [**MEDIAN**](../lib/statistics/median/Median.md) - Rolling Median - [**SKEW**](../lib/statistics/skew/Skew.md) - Skewness - [**STDDEV**](../lib/statistics/stddev/StdDev.md) - Standard Deviation -- [**VARIANCE**](../lib/statistics/variance/Variance.md) - Population and Sample Variance +- [**VARIANCE**](../lib/statistics/variance/Variance.md) - Population and Sample Variance \ No newline at end of file diff --git a/docs/validation.md b/docs/validation.md index 43406fcf..92f47325 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -58,7 +58,7 @@ | **Conditional Volatility** | Cv | - | - | - | - | | **Convolution Moving Average** | [Conv](../lib/trends/conv/conv.md) | ✔️ | ✔️ | ✔️ | ✔️ | | **Correlation** | Correlation | CORREL | - | Correlation | - | -| **Cumulative Mean (Average)** | Cummean | - | - | - | - | +| **Cumulative Moving Average** | [Cma](../lib/statistics/cma/Cma.md) | - | - | - | - | | **Decay Min-Max Channel** | Decaychannel | - | - | - | - | | **DeMark Pivot Points** | Pivotdem | - | - | - | ❔ | | **Detrended Price Oscillator** | Dpo | - | dpo | Dpo | ❔ | @@ -283,4 +283,4 @@ | **Median (Statistical)** | [Median](../lib/statistics/median/Median.md) | ✔️ | - | - | - | | **Skewness** | [Skew](../lib/statistics/skew/Skew.md) | ✔️ | - | - | - | | **Standard Deviation** | [StdDev](../lib/statistics/stddev/StdDev.md) | ✔️ | ✔️ | ✔️ | ✔️ | -| **Variance** | [Variance](../lib/statistics/variance/Variance.md) | ✔️ | ✔️ | ✔️ | ✔️ | +| **Variance** | [Variance](../lib/statistics/variance/Variance.md) | ✔️ | ✔️ | ✔️ | ✔️ | \ No newline at end of file diff --git a/lib/core/tbarseries/tbarseries.cs b/lib/core/tbarseries/tbarseries.cs index 5b514c7b..0a773d3c 100644 --- a/lib/core/tbarseries/tbarseries.cs +++ b/lib/core/tbarseries/tbarseries.cs @@ -9,10 +9,26 @@ namespace QuanTAlib; /// Implemented as struct to avoid heap allocations in high-frequency event dispatch. /// [StructLayout(LayoutKind.Auto)] -public readonly struct TBarEventArgs +public readonly struct TBarEventArgs : IEquatable { public TBar Value { get; init; } public bool IsNew { get; init; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool Equals(TBarEventArgs other) => + Value.Equals(other.Value) && IsNew == other.IsNew; + + public override bool Equals(object? obj) => + obj is TBarEventArgs other && Equals(other); + + public override int GetHashCode() => + HashCode.Combine(Value, IsNew); + + public static bool operator ==(TBarEventArgs left, TBarEventArgs right) => + left.Equals(right); + + public static bool operator !=(TBarEventArgs left, TBarEventArgs right) => + !left.Equals(right); } /// @@ -277,4 +293,4 @@ public class TBarSeries : IReadOnlyList IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); -} \ No newline at end of file +} diff --git a/lib/core/tseries/ITValuePublisher.cs b/lib/core/tseries/ITValuePublisher.cs index b51518c2..a715ba07 100644 --- a/lib/core/tseries/ITValuePublisher.cs +++ b/lib/core/tseries/ITValuePublisher.cs @@ -1,3 +1,4 @@ +using System.Runtime.CompilerServices; using System.Runtime.InteropServices; namespace QuanTAlib; @@ -7,10 +8,26 @@ namespace QuanTAlib; /// Implemented as struct to avoid heap allocations in high-frequency event dispatch. /// [StructLayout(LayoutKind.Auto)] -public readonly struct TValueEventArgs +public readonly struct TValueEventArgs : IEquatable { public TValue Value { get; init; } public bool IsNew { get; init; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool Equals(TValueEventArgs other) => + Value.Equals(other.Value) && IsNew == other.IsNew; + + public override bool Equals(object? obj) => + obj is TValueEventArgs other && Equals(other); + + public override int GetHashCode() => + HashCode.Combine(Value, IsNew); + + public static bool operator ==(TValueEventArgs left, TValueEventArgs right) => + left.Equals(right); + + public static bool operator !=(TValueEventArgs left, TValueEventArgs right) => + !left.Equals(right); } // Performance-focused event args struct; not derived from EventArgs by design. diff --git a/lib/feeds/gbm/ValidationHelper.cs b/lib/feeds/gbm/ValidationHelper.cs index a3065c0d..510f594d 100644 --- a/lib/feeds/gbm/ValidationHelper.cs +++ b/lib/feeds/gbm/ValidationHelper.cs @@ -434,7 +434,7 @@ public static class ValidationHelper for (int i = 0; i < qSeries.Count; i++) { double? sValue = selector(sSeries[i]); - if (!sValue.HasValue || sValue.Value == 0) continue; + if (!sValue.HasValue || Math.Abs(sValue.Value) < double.Epsilon) continue; double relDiff = Math.Abs((qSeries[i].Value - sValue.Value) / sValue.Value); if (relDiff > maxDiff) diff --git a/lib/momentum/adx/Adx.cs b/lib/momentum/adx/Adx.cs index cfb9608f..b302015d 100644 --- a/lib/momentum/adx/Adx.cs +++ b/lib/momentum/adx/Adx.cs @@ -29,6 +29,8 @@ namespace QuanTAlib; public sealed class Adx : ITValuePublisher { private readonly int _period; + private readonly double _decay; // (period - 1) / period for RMA + private readonly double _invPeriod; // 1 / period private TBar _prevBar; private TBar _p_prevBar; private bool _isInitialized; @@ -93,6 +95,8 @@ public sealed class Adx : ITValuePublisher throw new ArgumentException("Period must be greater than 0", nameof(period)); _period = period; + _decay = (period - 1.0) / period; + _invPeriod = 1.0 / period; Name = $"Adx({period})"; WarmupPeriod = period * 2; // Needs period for TR/DM smoothing, then period for ADX smoothing _isInitialized = false; @@ -211,16 +215,11 @@ public sealed class Adx : ITValuePublisher } else { - // RMA: Previous + (Input - Previous) / Period - // Or: Previous * (1 - 1/Period) + Input * (1/Period) - // Or: (Previous * (Period - 1) + Input) / Period - // Wilder uses sums, but effectively it's RMA. - // Standard formula: - // Smooth = Smooth - (Smooth / Period) + Input - - _trSmooth = _trSmooth - (_trSmooth / _period) + tr; - _dmPlusSmooth = _dmPlusSmooth - (_dmPlusSmooth / _period) + dmPlus; - _dmMinusSmooth = _dmMinusSmooth - (_dmMinusSmooth / _period) + dmMinus; + // RMA: Smooth = Smooth * decay + Input * invPeriod + // Using FMA for precision + _trSmooth = Math.FusedMultiplyAdd(_trSmooth, _decay, tr * _invPeriod); + _dmPlusSmooth = Math.FusedMultiplyAdd(_dmPlusSmooth, _decay, dmPlus * _invPeriod); + _dmMinusSmooth = Math.FusedMultiplyAdd(_dmMinusSmooth, _decay, dmMinus * _invPeriod); } // Calculate DI and DX @@ -250,13 +249,13 @@ public sealed class Adx : ITValuePublisher if (_dxSamples == _period) { - _adx = _dxSum / _period; // First ADX is SMA of DX + _adx = _dxSum * _invPeriod; // First ADX is SMA of DX } } else { - // ADX = (Prior ADX * (Period - 1) + Current DX) / Period - _adx = ((_adx * (_period - 1)) + dx) / _period; + // ADX = Prior ADX * decay + DX * invPeriod (RMA smoothing) + _adx = Math.FusedMultiplyAdd(_adx, _decay, dx * _invPeriod); } } @@ -339,9 +338,10 @@ public sealed class Adx : ITValuePublisher } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void Smooth(double input, int period, ref double smoothed) + private static void Smooth(double input, double decay, double invPeriod, ref double smoothed) { - smoothed = smoothed - (smoothed / period) + input; + // RMA: smoothed = smoothed * decay + input * invPeriod + smoothed = Math.FusedMultiplyAdd(smoothed, decay, input * invPeriod); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -354,6 +354,9 @@ public sealed class Adx : ITValuePublisher return; } + double decay = (period - 1.0) / period; + double invPeriod = 1.0 / period; + // Phase 1: Accumulate TR, +DM, -DM for the first 'period' bars double trSum = 0; double dmPlusSum = 0; @@ -387,9 +390,9 @@ public sealed class Adx : ITValuePublisher { CalcTrDm(i, high, low, close, out double tr, out double dmPlus, out double dmMinus); - Smooth(tr, period, ref trSmooth); - Smooth(dmPlus, period, ref dmPlusSmooth); - Smooth(dmMinus, period, ref dmMinusSmooth); + Smooth(tr, decay, invPeriod, ref trSmooth); + Smooth(dmPlus, decay, invPeriod, ref dmPlusSmooth); + Smooth(dmMinus, decay, invPeriod, ref dmMinusSmooth); dx = CalcDx(trSmooth, dmPlusSmooth, dmMinusSmooth); dxSum += dx; @@ -397,7 +400,7 @@ public sealed class Adx : ITValuePublisher } // Initialize ADX (SMA of DX) - double adx = dxSum / period; + double adx = dxSum * invPeriod; destination[adxStart] = adx; // Phase 3: Calculate ADX for the rest of the series @@ -405,14 +408,14 @@ public sealed class Adx : ITValuePublisher { CalcTrDm(i, high, low, close, out double tr, out double dmPlus, out double dmMinus); - Smooth(tr, period, ref trSmooth); - Smooth(dmPlus, period, ref dmPlusSmooth); - Smooth(dmMinus, period, ref dmMinusSmooth); + Smooth(tr, decay, invPeriod, ref trSmooth); + Smooth(dmPlus, decay, invPeriod, ref dmPlusSmooth); + Smooth(dmMinus, decay, invPeriod, ref dmMinusSmooth); dx = CalcDx(trSmooth, dmPlusSmooth, dmMinusSmooth); // ADX Smoothing (RMA) - Smooth(dx / period, period, ref adx); + Smooth(dx, decay, invPeriod, ref adx); destination[i] = adx; } } @@ -434,4 +437,4 @@ public sealed class Adx : ITValuePublisher return new TSeries(tList, [.. v]); } -} \ No newline at end of file +} diff --git a/lib/momentum/rsx/Rsx.cs b/lib/momentum/rsx/Rsx.cs index 0e649ee5..b7cfcb35 100644 --- a/lib/momentum/rsx/Rsx.cs +++ b/lib/momentum/rsx/Rsx.cs @@ -25,6 +25,7 @@ public sealed class Rsx : ITValuePublisher { private readonly int _period; private readonly double _alpha; + private readonly double _decay; [StructLayout(LayoutKind.Auto)] private record struct State @@ -72,6 +73,7 @@ public sealed class Rsx : ITValuePublisher _period = period; WarmupPeriod = period; _alpha = 3.0 / (period + 2.0); + _decay = 1.0 - _alpha; Name = $"Rsx({period})"; _handler = Handle; } @@ -129,33 +131,34 @@ public sealed class Rsx : ITValuePublisher _state.LastPrice = price; } - // --- Momentum Smoothing --- - double m1_1 = _state.M1_1 + _alpha * (momentum - _state.M1_1); - double m1_2 = _state.M1_2 + _alpha * (m1_1 - _state.M1_2); - double m1_out = (3.0 * m1_1 - m1_2) * 0.5; + // --- Momentum Smoothing (using FMA for precision and performance) --- + // EMA update: new = old + alpha * (input - old) = old * (1-alpha) + alpha * input = old * decay + alpha * input + double m1_1 = Math.FusedMultiplyAdd(_state.M1_1, _decay, _alpha * momentum); + double m1_2 = Math.FusedMultiplyAdd(_state.M1_2, _decay, _alpha * m1_1); + double m1_out = Math.FusedMultiplyAdd(3.0, m1_1, -m1_2) * 0.5; - double m2_1 = _state.M2_1 + _alpha * (m1_out - _state.M2_1); - double m2_2 = _state.M2_2 + _alpha * (m2_1 - _state.M2_2); - double m2_out = (3.0 * m2_1 - m2_2) * 0.5; + double m2_1 = Math.FusedMultiplyAdd(_state.M2_1, _decay, _alpha * m1_out); + double m2_2 = Math.FusedMultiplyAdd(_state.M2_2, _decay, _alpha * m2_1); + double m2_out = Math.FusedMultiplyAdd(3.0, m2_1, -m2_2) * 0.5; - double m3_1 = _state.M3_1 + _alpha * (m2_out - _state.M3_1); - double m3_2 = _state.M3_2 + _alpha * (m3_1 - _state.M3_2); - double smoothedMomentum = (3.0 * m3_1 - m3_2) * 0.5; + double m3_1 = Math.FusedMultiplyAdd(_state.M3_1, _decay, _alpha * m2_out); + double m3_2 = Math.FusedMultiplyAdd(_state.M3_2, _decay, _alpha * m3_1); + double smoothedMomentum = Math.FusedMultiplyAdd(3.0, m3_1, -m3_2) * 0.5; - // --- Absolute Momentum Smoothing --- + // --- Absolute Momentum Smoothing (using FMA) --- double absMomentum = Math.Abs(momentum); - double a1_1 = _state.A1_1 + _alpha * (absMomentum - _state.A1_1); - double a1_2 = _state.A1_2 + _alpha * (a1_1 - _state.A1_2); - double a1_out = (3.0 * a1_1 - a1_2) * 0.5; + double a1_1 = Math.FusedMultiplyAdd(_state.A1_1, _decay, _alpha * absMomentum); + double a1_2 = Math.FusedMultiplyAdd(_state.A1_2, _decay, _alpha * a1_1); + double a1_out = Math.FusedMultiplyAdd(3.0, a1_1, -a1_2) * 0.5; - double a2_1 = _state.A2_1 + _alpha * (a1_out - _state.A2_1); - double a2_2 = _state.A2_2 + _alpha * (a2_1 - _state.A2_2); - double a2_out = (3.0 * a2_1 - a2_2) * 0.5; + double a2_1 = Math.FusedMultiplyAdd(_state.A2_1, _decay, _alpha * a1_out); + double a2_2 = Math.FusedMultiplyAdd(_state.A2_2, _decay, _alpha * a2_1); + double a2_out = Math.FusedMultiplyAdd(3.0, a2_1, -a2_2) * 0.5; - double a3_1 = _state.A3_1 + _alpha * (a2_out - _state.A3_1); - double a3_2 = _state.A3_2 + _alpha * (a3_1 - _state.A3_2); - double smoothedAbsMomentum = (3.0 * a3_1 - a3_2) * 0.5; + double a3_1 = Math.FusedMultiplyAdd(_state.A3_1, _decay, _alpha * a2_out); + double a3_2 = Math.FusedMultiplyAdd(_state.A3_2, _decay, _alpha * a3_1); + double smoothedAbsMomentum = Math.FusedMultiplyAdd(3.0, a3_1, -a3_2) * 0.5; if (isNew) { @@ -231,6 +234,7 @@ public sealed class Rsx : ITValuePublisher if (len == 0) return; double alpha = 3.0 / (period + 2.0); + double decay = 1.0 - alpha; // Momentum filters double m1_1 = 0, m1_2 = 0; @@ -267,33 +271,33 @@ public sealed class Rsx : ITValuePublisher double momentum = (price - lastPrice) * 100.0; lastPrice = price; - // Momentum Smoothing - m1_1 += alpha * (momentum - m1_1); - m1_2 += alpha * (m1_1 - m1_2); - double m1_out = (3.0 * m1_1 - m1_2) * 0.5; + // Momentum Smoothing (using FMA for precision and performance) + m1_1 = Math.FusedMultiplyAdd(m1_1, decay, alpha * momentum); + m1_2 = Math.FusedMultiplyAdd(m1_2, decay, alpha * m1_1); + double m1_out = Math.FusedMultiplyAdd(3.0, m1_1, -m1_2) * 0.5; - m2_1 += alpha * (m1_out - m2_1); - m2_2 += alpha * (m2_1 - m2_2); - double m2_out = (3.0 * m2_1 - m2_2) * 0.5; + m2_1 = Math.FusedMultiplyAdd(m2_1, decay, alpha * m1_out); + m2_2 = Math.FusedMultiplyAdd(m2_2, decay, alpha * m2_1); + double m2_out = Math.FusedMultiplyAdd(3.0, m2_1, -m2_2) * 0.5; - m3_1 += alpha * (m2_out - m3_1); - m3_2 += alpha * (m3_1 - m3_2); - double smoothedMomentum = (3.0 * m3_1 - m3_2) * 0.5; + m3_1 = Math.FusedMultiplyAdd(m3_1, decay, alpha * m2_out); + m3_2 = Math.FusedMultiplyAdd(m3_2, decay, alpha * m3_1); + double smoothedMomentum = Math.FusedMultiplyAdd(3.0, m3_1, -m3_2) * 0.5; - // Abs Momentum Smoothing + // Abs Momentum Smoothing (using FMA) double absMomentum = Math.Abs(momentum); - a1_1 += alpha * (absMomentum - a1_1); - a1_2 += alpha * (a1_1 - a1_2); - double a1_out = (3.0 * a1_1 - a1_2) * 0.5; + a1_1 = Math.FusedMultiplyAdd(a1_1, decay, alpha * absMomentum); + a1_2 = Math.FusedMultiplyAdd(a1_2, decay, alpha * a1_1); + double a1_out = Math.FusedMultiplyAdd(3.0, a1_1, -a1_2) * 0.5; - a2_1 += alpha * (a1_out - a2_1); - a2_2 += alpha * (a2_1 - a2_2); - double a2_out = (3.0 * a2_1 - a2_2) * 0.5; + a2_1 = Math.FusedMultiplyAdd(a2_1, decay, alpha * a1_out); + a2_2 = Math.FusedMultiplyAdd(a2_2, decay, alpha * a2_1); + double a2_out = Math.FusedMultiplyAdd(3.0, a2_1, -a2_2) * 0.5; - a3_1 += alpha * (a2_out - a3_1); - a3_2 += alpha * (a3_1 - a3_2); - double smoothedAbsMomentum = (3.0 * a3_1 - a3_2) * 0.5; + a3_1 = Math.FusedMultiplyAdd(a3_1, decay, alpha * a2_out); + a3_2 = Math.FusedMultiplyAdd(a3_2, decay, alpha * a3_1); + double smoothedAbsMomentum = Math.FusedMultiplyAdd(3.0, a3_1, -a3_2) * 0.5; // Final RSX double rsx; diff --git a/lib/statistics/_index.md b/lib/statistics/_index.md index 458f61d2..fec4fc08 100644 --- a/lib/statistics/_index.md +++ b/lib/statistics/_index.md @@ -9,7 +9,7 @@ Statistical analysis tools applied to price/returns. | COINTEGRATION | Cointegration | | | CORRELATION | Correlation (Pearson's) | | | [COVARIANCE](covariance/Covariance.md) | Covariance | | -| CUMMEAN | Cumulative Mean (Average) | | +| [CMA](cma/Cma.md) | Cumulative Moving Average | Running average of ALL values (Welford's algorithm). No window. | | ENTROPY | Normalized Shannon Entropy | | | GEOMEAN | Geometric Mean | | | GRANGER | Granger Causality Test | | @@ -30,4 +30,4 @@ Statistical analysis tools applied to price/returns. | THEIL | Theil Index | | | [VARIANCE](variance/Variance.md) | Variance | | | ZSCORE | Z-score standardization | | -| ZTEST | Z-Test | | +| ZTEST | Z-Test | | \ No newline at end of file diff --git a/lib/statistics/cma/Cma.Quantower.Tests.cs b/lib/statistics/cma/Cma.Quantower.Tests.cs new file mode 100644 index 00000000..9f0ccd43 --- /dev/null +++ b/lib/statistics/cma/Cma.Quantower.Tests.cs @@ -0,0 +1,169 @@ +using TradingPlatform.BusinessLayer; + +namespace QuanTAlib.Tests; + +public class CmaIndicatorTests +{ + [Fact] + public void CmaIndicator_Constructor_SetsDefaults() + { + var indicator = new CmaIndicator(); + + Assert.Equal(SourceType.Close, indicator.Source); + Assert.True(indicator.ShowColdValues); + Assert.Equal("CMA - Cumulative Moving Average", indicator.Name); + Assert.False(indicator.SeparateWindow); + Assert.True(indicator.OnBackGround); + } + + [Fact] + public void CmaIndicator_MinHistoryDepths_EqualsZero() + { + var indicator = new CmaIndicator(); + + Assert.Equal(0, CmaIndicator.MinHistoryDepths); + Assert.Equal(0, ((IWatchlistIndicator)indicator).MinHistoryDepths); + } + + [Fact] + public void CmaIndicator_ShortName_IncludesSource() + { + var indicator = new CmaIndicator(); + + Assert.Contains("CMA", indicator.ShortName, StringComparison.Ordinal); + } + + [Fact] + public void CmaIndicator_Initialize_CreatesInternalCma() + { + var indicator = new CmaIndicator(); + + // Initialize should not throw + indicator.Initialize(); + + // After init, line series should exist + Assert.Single(indicator.LinesSeries); + } + + [Fact] + public void CmaIndicator_ProcessUpdate_HistoricalBar_ComputesValue() + { + var indicator = new CmaIndicator(); + indicator.Initialize(); + + // Add historical data + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + + // Process update + var args = new UpdateArgs(UpdateReason.HistoricalBar); + indicator.ProcessUpdate(args); + + // Line series should have a value + Assert.Equal(1, indicator.LinesSeries[0].Count); + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0))); + } + + [Fact] + public void CmaIndicator_ProcessUpdate_NewBar_ComputesValue() + { + var indicator = new CmaIndicator(); + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + indicator.HistoricalData.AddBar(now.AddMinutes(1), 102, 108, 100, 106); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewBar)); + + Assert.Equal(2, indicator.LinesSeries[0].Count); + } + + [Fact] + public void CmaIndicator_ProcessUpdate_NewTick_ProcessesWithoutError() + { + var indicator = new CmaIndicator(); + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + double firstValue = indicator.LinesSeries[0].GetValue(0); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewTick)); + double secondValue = indicator.LinesSeries[0].GetValue(0); + + Assert.True(double.IsFinite(firstValue)); + Assert.True(double.IsFinite(secondValue)); + } + + [Fact] + public void CmaIndicator_MultipleUpdates_ProducesCorrectCmaSequence() + { + var indicator = new CmaIndicator(); + indicator.Initialize(); + + var now = DateTime.UtcNow; + double[] closes = { 100, 102, 104, 103, 105 }; + + foreach (var close in closes) + { + indicator.HistoricalData.AddBar(now, close, close + 2, close - 2, close); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + now = now.AddMinutes(1); + } + + // All values should be finite + for (int i = 0; i < closes.Length; i++) + { + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(closes.Length - 1 - i))); + } + + // Last CMA should be average of all values: (100 + 102 + 104 + 103 + 105) / 5 = 102.8 + double lastCma = indicator.LinesSeries[0].GetValue(0); + Assert.Equal(102.8, lastCma, 1e-10); + } + + [Fact] + public void CmaIndicator_DifferentSourceTypes_Work() + { + var sources = new[] { SourceType.Open, SourceType.High, SourceType.Low, SourceType.Close, SourceType.HL2, SourceType.HLC3 }; + + foreach (var source in sources) + { + var indicator = new CmaIndicator { Source = source }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 110, 90, 105); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0)), + $"Source {source} should produce finite value"); + } + } + + [Fact] + public void CmaIndicator_CalculatesRunningAverage() + { + var indicator = new CmaIndicator(); + indicator.Initialize(); + + var now = DateTime.UtcNow; + + // Add bars with known close prices: 10, 20, 30 + indicator.HistoricalData.AddBar(now, 10, 10, 10, 10); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + Assert.Equal(10.0, indicator.LinesSeries[0].GetValue(0), 1e-10); // CMA = 10 + + indicator.HistoricalData.AddBar(now.AddMinutes(1), 20, 20, 20, 20); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewBar)); + Assert.Equal(15.0, indicator.LinesSeries[0].GetValue(0), 1e-10); // CMA = (10+20)/2 = 15 + + indicator.HistoricalData.AddBar(now.AddMinutes(2), 30, 30, 30, 30); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewBar)); + Assert.Equal(20.0, indicator.LinesSeries[0].GetValue(0), 1e-10); // CMA = (10+20+30)/3 = 20 + } +} \ No newline at end of file diff --git a/lib/statistics/cma/Cma.Quantower.cs b/lib/statistics/cma/Cma.Quantower.cs new file mode 100644 index 00000000..d79ed2e8 --- /dev/null +++ b/lib/statistics/cma/Cma.Quantower.cs @@ -0,0 +1,52 @@ +using System.Drawing; +using System.Runtime.CompilerServices; +using TradingPlatform.BusinessLayer; + +namespace QuanTAlib; + +[SkipLocalsInit] +public sealed class CmaIndicator : Indicator, IWatchlistIndicator +{ + [IndicatorExtensions.DataSourceInput] + public SourceType Source { get; set; } = SourceType.Close; + + [InputParameter("Show cold values", sortIndex: 21)] + public bool ShowColdValues { get; set; } = true; + + private Cma? _cma; + private readonly LineSeries? _series; + private string? _sourceName; + private Func? _priceSelector; + + public static int MinHistoryDepths => 0; + int IWatchlistIndicator.MinHistoryDepths => MinHistoryDepths; + + public override string ShortName => $"CMA:{_sourceName}"; + + public CmaIndicator() + { + OnBackGround = true; + SeparateWindow = false; + Name = "CMA - Cumulative Moving Average"; + Description = "Cumulative Moving Average (Running Average)"; + _series = new(name: "CMA", color: IndicatorExtensions.Averages, width: 2, style: LineStyle.Solid); + AddLineSeries(_series); + } + + protected override void OnInit() + { + _priceSelector = Source.GetPriceSelector(); + _sourceName = Source.ToString(); + _cma = new Cma(); + base.OnInit(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + protected override void OnUpdate(UpdateArgs args) + { + bool isNew = args.IsNewBar(); + var item = HistoricalData[Count - 1, SeekOriginHistory.Begin]; + double value = _cma!.Update(new TValue(item.TimeLeft.Ticks, _priceSelector!(item)), isNew).Value; + _series!.SetValue(value, _cma.IsHot, ShowColdValues); + } +} \ No newline at end of file diff --git a/lib/statistics/cma/Cma.Tests.cs b/lib/statistics/cma/Cma.Tests.cs new file mode 100644 index 00000000..8aa2d6a6 --- /dev/null +++ b/lib/statistics/cma/Cma.Tests.cs @@ -0,0 +1,576 @@ +namespace QuanTAlib.Tests; + +public class CmaTests +{ + [Fact] + public void Cma_Calc_ReturnsValue() + { + var cma = new Cma(); + + Assert.Equal(0, cma.Last.Value); + + TValue result = cma.Update(new TValue(DateTime.UtcNow, 100)); + + Assert.True(result.Value > 0); + Assert.Equal(result.Value, cma.Last.Value); + } + + [Fact] + public void Cma_FirstValue_ReturnsItself() + { + var cma = new Cma(); + + TValue result = cma.Update(new TValue(DateTime.UtcNow, 100)); + + Assert.Equal(100.0, result.Value, 1e-10); + } + + [Fact] + public void Cma_Calc_IsNew_AcceptsParameter() + { + var cma = new Cma(); + + cma.Update(new TValue(DateTime.UtcNow, 100), isNew: true); + double value1 = cma.Last.Value; + + cma.Update(new TValue(DateTime.UtcNow, 200), isNew: true); + double value2 = cma.Last.Value; + + // Values should change with new bars + Assert.NotEqual(value1, value2); + } + + [Fact] + public void Cma_Calc_IsNew_False_UpdatesValue() + { + var cma = new Cma(); + + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, 110), isNew: true); + double beforeUpdate = cma.Last.Value; + + cma.Update(new TValue(DateTime.UtcNow, 120), isNew: false); + double afterUpdate = cma.Last.Value; + + // Update should change the value + Assert.NotEqual(beforeUpdate, afterUpdate); + } + + [Fact] + public void Cma_Reset_ClearsState() + { + var cma = new Cma(); + + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, 105)); + double valueBefore = cma.Last.Value; + + cma.Reset(); + + Assert.Equal(0, cma.Last.Value); + Assert.False(cma.IsHot); + + // After reset, should accept new values + cma.Update(new TValue(DateTime.UtcNow, 50)); + Assert.NotEqual(0, cma.Last.Value); + Assert.NotEqual(valueBefore, cma.Last.Value); + } + + [Fact] + public void Cma_Properties_Accessible() + { + var cma = new Cma(); + + Assert.Equal(0, cma.Last.Value); + Assert.False(cma.IsHot); + + cma.Update(new TValue(DateTime.UtcNow, 100)); + + Assert.NotEqual(0, cma.Last.Value); + Assert.True(cma.IsHot); + } + + [Fact] + public void Cma_IsHot_BecomesTrueAfterFirstValue() + { + var cma = new Cma(); + + Assert.False(cma.IsHot); + + cma.Update(new TValue(DateTime.UtcNow, 100)); + + Assert.True(cma.IsHot); + } + + [Fact] + public void Cma_CalculatesCorrectAverage() + { + var cma = new Cma(); + + cma.Update(new TValue(DateTime.UtcNow, 10)); + Assert.Equal(10.0, cma.Last.Value, 1e-10); // (10)/1 = 10 + + cma.Update(new TValue(DateTime.UtcNow, 20)); + Assert.Equal(15.0, cma.Last.Value, 1e-10); // (10+20)/2 = 15 + + cma.Update(new TValue(DateTime.UtcNow, 30)); + Assert.Equal(20.0, cma.Last.Value, 1e-10); // (10+20+30)/3 = 20 + + cma.Update(new TValue(DateTime.UtcNow, 40)); + Assert.Equal(25.0, cma.Last.Value, 1e-10); // (10+20+30+40)/4 = 25 + + cma.Update(new TValue(DateTime.UtcNow, 50)); + Assert.Equal(30.0, cma.Last.Value, 1e-10); // (10+20+30+40+50)/5 = 30 + } + + [Fact] + public void Cma_IncludesAllValues_NoSlidingWindow() + { + var cma = new Cma(); + + // Add 10 values: 10, 20, 30, ..., 100 + for (int i = 1; i <= 10; i++) + { + cma.Update(new TValue(DateTime.UtcNow, i * 10)); + } + + // CMA of 10,20,30,40,50,60,70,80,90,100 = 550/10 = 55 + Assert.Equal(55.0, cma.Last.Value, 1e-10); + + // Add one more value + cma.Update(new TValue(DateTime.UtcNow, 110)); + + // CMA now includes ALL 11 values: (550 + 110)/11 = 660/11 = 60 + Assert.Equal(60.0, cma.Last.Value, 1e-10); + } + + [Fact] + public void Cma_IterativeCorrections_RestoreToOriginalState() + { + var cma = new Cma(); + var gbm = new GBM(startPrice: 100.0, mu: 0.02, sigma: 0.1); + + // Feed 10 new values + TValue tenthInput = default; + for (int i = 0; i < 10; i++) + { + var bar = gbm.Next(isNew: true); + tenthInput = new TValue(bar.Time, bar.Close); + cma.Update(tenthInput, isNew: true); + } + + // Remember CMA state after 10 values + double cmaAfterTen = cma.Last.Value; + + // Generate 9 corrections with isNew=false (different values) + for (int i = 0; i < 9; i++) + { + var bar = gbm.Next(isNew: false); + cma.Update(new TValue(bar.Time, bar.Close), isNew: false); + } + + // Feed the remembered 10th input again with isNew=false + TValue finalCma = cma.Update(tenthInput, isNew: false); + + // CMA should match the original state after 10 values + Assert.Equal(cmaAfterTen, finalCma.Value, 1e-10); + } + + [Fact] + public void Cma_BatchCalc_MatchesIterativeCalc() + { + var cmaIterative = new Cma(); + var cmaBatch = new Cma(); + var gbm = new GBM(startPrice: 100.0, mu: 0.02, sigma: 0.1); + + // Generate data + var series = new TSeries(); + for (int i = 0; i < 100; i++) + { + var bar = gbm.Next(isNew: true); + series.Add(bar.Time, bar.Close); + } + + Assert.True(series.Count > 0); + + // Calculate iteratively + var iterativeResults = new TSeries(); + foreach (var item in series) + { + iterativeResults.Add(cmaIterative.Update(item)); + } + + // Calculate batch + var batchResults = cmaBatch.Update(series); + + // Compare + Assert.Equal(iterativeResults.Count, batchResults.Count); + for (int i = 0; i < iterativeResults.Count; i++) + { + Assert.Equal(iterativeResults[i].Value, batchResults[i].Value, 1e-10); + Assert.Equal(iterativeResults[i].Time, batchResults[i].Time); + } + } + + [Fact] + public void Cma_NaN_Input_UsesLastValidValue() + { + var cma = new Cma(); + + // Feed some valid values + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, 110)); + + // Feed NaN - should use last valid value (110) + var resultAfterNaN = cma.Update(new TValue(DateTime.UtcNow, double.NaN)); + + // Result should be finite (not NaN) + Assert.True(double.IsFinite(resultAfterNaN.Value)); + Assert.NotEqual(0, resultAfterNaN.Value); + } + + [Fact] + public void Cma_Infinity_Input_UsesLastValidValue() + { + var cma = new Cma(); + + // Feed some valid values + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, 110)); + + // Feed positive infinity - should use last valid value + var resultAfterPosInf = cma.Update(new TValue(DateTime.UtcNow, double.PositiveInfinity)); + Assert.True(double.IsFinite(resultAfterPosInf.Value)); + + // Feed negative infinity - should use last valid value + var resultAfterNegInf = cma.Update(new TValue(DateTime.UtcNow, double.NegativeInfinity)); + Assert.True(double.IsFinite(resultAfterNegInf.Value)); + } + + [Fact] + public void Cma_MultipleNaN_ContinuesWithLastValid() + { + var cma = new Cma(); + + // Feed valid values + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, 110)); + cma.Update(new TValue(DateTime.UtcNow, 120)); + + // Feed multiple NaN values + var r1 = cma.Update(new TValue(DateTime.UtcNow, double.NaN)); + var r2 = cma.Update(new TValue(DateTime.UtcNow, double.NaN)); + var r3 = cma.Update(new TValue(DateTime.UtcNow, double.NaN)); + + // All results should be finite + Assert.True(double.IsFinite(r1.Value)); + Assert.True(double.IsFinite(r2.Value)); + Assert.True(double.IsFinite(r3.Value)); + } + + [Fact] + public void Cma_BatchCalc_HandlesNaN() + { + var cma = new Cma(); + + // Create series with NaN values interspersed + var series = new TSeries(); + series.Add(DateTime.UtcNow.Ticks, 100); + series.Add(DateTime.UtcNow.Ticks + 1, 110); + series.Add(DateTime.UtcNow.Ticks + 2, double.NaN); + series.Add(DateTime.UtcNow.Ticks + 3, 120); + series.Add(DateTime.UtcNow.Ticks + 4, double.PositiveInfinity); + series.Add(DateTime.UtcNow.Ticks + 5, 130); + + var results = cma.Update(series); + + // All results should be finite + foreach (var result in results) + { + Assert.True(double.IsFinite(result.Value), $"Expected finite value but got {result.Value}"); + } + } + + [Fact] + public void Cma_Reset_ClearsLastValidValue() + { + var cma = new Cma(); + + // Feed values including NaN + cma.Update(new TValue(DateTime.UtcNow, 100)); + cma.Update(new TValue(DateTime.UtcNow, double.NaN)); + + // Reset + cma.Reset(); + + // After reset, first valid value should establish new baseline + var result = cma.Update(new TValue(DateTime.UtcNow, 50)); + Assert.Equal(50.0, result.Value, 1e-10); + } + + [Fact] + public void Cma_StaticBatch_Works() + { + var series = new TSeries(); + series.Add(DateTime.UtcNow.Ticks, 10); + series.Add(DateTime.UtcNow.Ticks + 1, 20); + series.Add(DateTime.UtcNow.Ticks + 2, 30); + series.Add(DateTime.UtcNow.Ticks + 3, 40); + series.Add(DateTime.UtcNow.Ticks + 4, 50); + + var results = Cma.Batch(series); + + Assert.Equal(5, results.Count); + // CMA for last value: (10+20+30+40+50)/5 = 30 + Assert.Equal(30.0, results.Last.Value, 1e-10); + } + + [Fact] + public void Cma_FlatLine_ReturnsSameValue() + { + var cma = new Cma(); + + for (int i = 0; i < 20; i++) + { + cma.Update(new TValue(DateTime.UtcNow, 100)); + } + + Assert.Equal(100.0, cma.Last.Value, 1e-10); + } + + // ============== Span API Tests ============== + + [Fact] + public void Cma_SpanBatch_ValidatesInput() + { + double[] source = [1, 2, 3, 4, 5]; + double[] wrongSizeOutput = new double[3]; + + // Output must be same length as source + Assert.Throws(() => Cma.Batch(source.AsSpan(), wrongSizeOutput.AsSpan())); + } + + [Fact] + public void Cma_SpanBatch_MatchesTSeriesBatch() + { + var series = new TSeries(); + double[] source = new double[100]; + double[] output = new double[100]; + + var gbm = new GBM(startPrice: 100.0, mu: 0.02, sigma: 0.1, seed: 42); + for (int i = 0; i < 100; i++) + { + var bar = gbm.Next(isNew: true); + source[i] = bar.Close; + series.Add(bar.Time, bar.Close); + } + + // Calculate with TSeries API + var tseriesResult = Cma.Batch(series); + + // Calculate with Span API + Cma.Batch(source.AsSpan(), output.AsSpan()); + + // Compare results + for (int i = 0; i < 100; i++) + { + Assert.Equal(tseriesResult[i].Value, output[i], 1e-10); + } + } + + [Fact] + public void Cma_SpanBatch_CalculatesCorrectly() + { + double[] source = [10, 20, 30, 40, 50]; + double[] output = new double[5]; + + Cma.Batch(source.AsSpan(), output.AsSpan()); + + Assert.Equal(10.0, output[0], 1e-10); // 10/1 = 10 + Assert.Equal(15.0, output[1], 1e-10); // (10+20)/2 = 15 + Assert.Equal(20.0, output[2], 1e-10); // (10+20+30)/3 = 20 + Assert.Equal(25.0, output[3], 1e-10); // (10+20+30+40)/4 = 25 + Assert.Equal(30.0, output[4], 1e-10); // (10+20+30+40+50)/5 = 30 + } + + [Fact] + public void Cma_SpanBatch_ZeroAllocation() + { + double[] source = new double[10000]; + double[] output = new double[10000]; + + var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 42); + for (int i = 0; i < source.Length; i++) + source[i] = gbm.Next().Close; + + // Warm up + Cma.Batch(source.AsSpan(), output.AsSpan()); + + // This test verifies the method runs without throwing + Assert.True(double.IsFinite(output[^1])); + } + + [Fact] + public void Cma_SpanBatch_HandlesNaN() + { + double[] source = [100, 110, double.NaN, 120, 130]; + double[] output = new double[5]; + + Cma.Batch(source.AsSpan(), output.AsSpan()); + + // All outputs should be finite + foreach (var val in output) + { + Assert.True(double.IsFinite(val), $"Expected finite value but got {val}"); + } + } + + [Fact] + public void Cma_AllModes_ProduceSameResult() + { + // Arrange + var gbm = new GBM(startPrice: 100, mu: 0.05, sigma: 0.2, seed: 123); + var bars = gbm.Fetch(1000, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + var series = bars.Close; + + // 1. Batch Mode + var batchSeries = Cma.Batch(series); + double expected = batchSeries.Last.Value; + + // 2. Span Mode + var tValues = series.Values.ToArray(); + var spanInput = new ReadOnlySpan(tValues); + var spanOutput = new double[tValues.Length]; + Cma.Batch(spanInput, spanOutput); + double spanResult = spanOutput[^1]; + + // 3. Streaming Mode + var streamingInd = new Cma(); + for (int i = 0; i < series.Count; i++) + { + streamingInd.Update(series[i]); + } + double streamingResult = streamingInd.Last.Value; + + // 4. Eventing Mode + var pubSource = new TSeries(); + var eventingInd = new Cma(pubSource); + for (int i = 0; i < series.Count; i++) + { + pubSource.Add(series[i]); + } + double eventingResult = eventingInd.Last.Value; + + // Assert + Assert.Equal(expected, spanResult, precision: 9); + Assert.Equal(expected, streamingResult, precision: 9); + Assert.Equal(expected, eventingResult, precision: 9); + } + + [Fact] + public void Chainability_Works() + { + var source = new TSeries(); + var cma = new Cma(source); + + source.Add(new TValue(DateTime.UtcNow, 100)); + Assert.Equal(100, cma.Last.Value); + } + + [Fact] + public void WarmupPeriod_IsSetCorrectly() + { + var cma = new Cma(); + Assert.Equal(1, cma.WarmupPeriod); + } + + [Fact] + public void Prime_SetsStateCorrectly() + { + var cma = new Cma(); + double[] history = [10, 20, 30, 40, 50]; // CMA = 30 + + cma.Prime(history); + + Assert.True(cma.IsHot); + Assert.Equal(30.0, cma.Last.Value, 1e-10); + + // Verify it continues correctly + cma.Update(new TValue(DateTime.UtcNow, 60)); // (10+20+30+40+50+60)/6 = 35 + Assert.Equal(35.0, cma.Last.Value, 1e-10); + } + + [Fact] + public void Prime_HandlesNaN_InHistory() + { + var cma = new Cma(); + double[] history = [10, 20, double.NaN, 40]; + // 10 -> 10 + // 10, 20 -> 15 + // 10, 20, 20 (NaN replaced by 20) -> 16.666... + // 10, 20, 20, 40 -> 22.5 + + cma.Prime(history); + + Assert.True(cma.IsHot); + Assert.Equal(22.5, cma.Last.Value, 1e-9); + } + + [Fact] + public void Calculate_ReturnsCorrectResultsAndHotIndicator() + { + var series = new TSeries(); + for (int i = 1; i <= 10; i++) series.Add(DateTime.UtcNow, i * 10); + // 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 + + var (results, indicator) = Cma.Calculate(series); + + // Check results + Assert.Equal(10, results.Count); + Assert.Equal(30.0, results[4].Value, 1e-10); // CMA after 5 values = 30 + Assert.Equal(55.0, results.Last.Value, 1e-10); // CMA of all 10 = 55 + + // Check indicator state + Assert.True(indicator.IsHot); + Assert.Equal(55.0, indicator.Last.Value, 1e-10); + Assert.Equal(1, indicator.WarmupPeriod); + + // Verify indicator continues correctly + indicator.Update(new TValue(DateTime.UtcNow, 110)); + // CMA now = (550 + 110)/11 = 60 + Assert.Equal(60.0, indicator.Last.Value, 1e-10); + } + + [Fact] + public void Cma_NumericalStability_LargeDataset() + { + // Test that CMA remains stable over a large number of values + var cma = new Cma(); + double expectedSum = 0; + + for (int i = 1; i <= 100000; i++) + { + cma.Update(new TValue(DateTime.UtcNow, 100.0)); // All same value + expectedSum += 100.0; + } + + // CMA of 100000 values all equal to 100 should be exactly 100 + Assert.Equal(100.0, cma.Last.Value, 1e-9); + } + + [Fact] + public void Cma_NumericalStability_VaryingValues() + { + // Test with alternating values + var cma = new Cma(); + + for (int i = 0; i < 10000; i++) + { + double value = (i % 2 == 0) ? 100.0 : 200.0; + cma.Update(new TValue(DateTime.UtcNow, value)); + } + + // CMA of alternating 100, 200 should converge to 150 + Assert.Equal(150.0, cma.Last.Value, 1e-9); + } +} \ No newline at end of file diff --git a/lib/statistics/cma/Cma.Validation.Tests.cs b/lib/statistics/cma/Cma.Validation.Tests.cs new file mode 100644 index 00000000..dd21f065 --- /dev/null +++ b/lib/statistics/cma/Cma.Validation.Tests.cs @@ -0,0 +1,318 @@ +using Xunit.Abstractions; + +namespace QuanTAlib.Tests; + +/// +/// Validation tests for CMA (Cumulative Moving Average). +/// CMA is not commonly found in standard TA libraries (like TA-Lib, Skender, etc.) +/// as it's a fundamental statistical concept rather than a trading indicator. +/// These tests validate against known mathematical results. +/// +public sealed class CmaValidationTests : IDisposable +{ + private readonly ValidationTestData _testData; + private readonly ITestOutputHelper _output; + private bool _disposed; + + public CmaValidationTests(ITestOutputHelper output) + { + _output = output; + _testData = new ValidationTestData(); + } + + public void Dispose() + { + Dispose(true); + } + + private void Dispose(bool disposing) + { + if (_disposed) + { + return; + } + + _disposed = true; + + if (disposing) + { + _testData?.Dispose(); + } + } + + [Fact] + public void Validate_MathematicalCorrectness_Batch() + { + // Calculate QuanTAlib CMA (batch TSeries) + var cma = new Cma(); + var qResult = cma.Update(_testData.Data); + + // Calculate expected CMA manually using running sum + double runningSum = 0; + int count = 0; + + foreach (var item in _testData.Data) + { + count++; + runningSum += item.Value; + double expectedCma = runningSum / count; + + // Get corresponding QuanTAlib result + double qValue = qResult[count - 1].Value; + + Assert.True( + Math.Abs(qValue - expectedCma) <= ValidationHelper.DefaultTolerance, + $"Mismatch at index {count - 1}: QuanTAlib={qValue:G17}, Expected={expectedCma:G17}"); + } + + _output.WriteLine("CMA Batch(TSeries) validated successfully against manual calculation"); + } + + [Fact] + public void Validate_MathematicalCorrectness_Streaming() + { + // Calculate QuanTAlib CMA (streaming) + var cma = new Cma(); + var qResults = new List(); + foreach (var item in _testData.Data) + { + qResults.Add(cma.Update(item).Value); + } + + // Calculate expected CMA manually + double runningSum = 0; + + for (int i = 0; i < _testData.Data.Count; i++) + { + runningSum += _testData.Data[i].Value; + double expectedCma = runningSum / (i + 1); + + Assert.True( + Math.Abs(qResults[i] - expectedCma) <= ValidationHelper.DefaultTolerance, + $"Mismatch at index {i}: QuanTAlib={qResults[i]:G17}, Expected={expectedCma:G17}"); + } + + _output.WriteLine("CMA Streaming validated successfully against manual calculation"); + } + + [Fact] + public void Validate_MathematicalCorrectness_Span() + { + // Prepare data for Span API + double[] sourceData = _testData.RawData.ToArray(); + double[] qOutput = new double[sourceData.Length]; + + // Calculate QuanTAlib CMA (Span API) + Cma.Batch(sourceData.AsSpan(), qOutput.AsSpan()); + + // Calculate expected CMA manually + double runningSum = 0; + + for (int i = 0; i < sourceData.Length; i++) + { + runningSum += sourceData[i]; + double expectedCma = runningSum / (i + 1); + + Assert.True( + Math.Abs(qOutput[i] - expectedCma) <= ValidationHelper.DefaultTolerance, + $"Mismatch at index {i}: QuanTAlib={qOutput[i]:G17}, Expected={expectedCma:G17}"); + } + + _output.WriteLine("CMA Span validated successfully against manual calculation"); + } + + [Fact] + public void Validate_WelfordAlgorithm_Stability() + { + // Test numerical stability with large values + // Welford's algorithm should handle this without overflow + var cma = new Cma(); + double[] largeValues = new double[1000]; + double baseValue = 1e10; + + for (int i = 0; i < largeValues.Length; i++) + { + largeValues[i] = baseValue + i; + } + + // Calculate CMA + foreach (var val in largeValues) + { + cma.Update(new TValue(DateTime.UtcNow, val)); + } + + // Expected: average of 1e10, 1e10+1, ..., 1e10+999 + // = 1e10 + average of 0,1,2,...,999 + // = 1e10 + 499.5 + double expectedMean = baseValue + 499.5; + + Assert.Equal(expectedMean, cma.Last.Value, 1e-6); + _output.WriteLine($"CMA Welford stability test passed: {cma.Last.Value:G17}"); + } + + [Fact] + public void Validate_WelfordAlgorithm_SmallDifferences() + { + // Test with values that have small differences (challenges precision) + var cma = new Cma(); + double[] values = new double[10000]; + double baseValue = 1e8; + + for (int i = 0; i < values.Length; i++) + { + values[i] = baseValue + (i % 2 == 0 ? 0.1 : -0.1); + } + + foreach (var val in values) + { + cma.Update(new TValue(DateTime.UtcNow, val)); + } + + // With alternating +0.1 and -0.1, the average offset is 0 + Assert.Equal(baseValue, cma.Last.Value, 1e-7); + _output.WriteLine($"CMA small differences test passed: {cma.Last.Value:G17}"); + } + + [Fact] + public void Validate_AgainstNaiveSum_ShortSequence() + { + // For short sequences, compare against naive sum/count + double[] values = [100, 200, 150, 175, 125, 180, 160, 140, 190, 170]; + var cma = new Cma(); + + double sum = 0; + for (int i = 0; i < values.Length; i++) + { + sum += values[i]; + cma.Update(new TValue(DateTime.UtcNow, values[i])); + + double naiveMean = sum / (i + 1); + Assert.Equal(naiveMean, cma.Last.Value, 1e-10); + } + + _output.WriteLine("CMA validated against naive sum for short sequence"); + } + + [Fact] + public void Validate_AgainstNaiveSum_LongSequence() + { + // For longer sequences, verify the final value + var gbm = new GBM(startPrice: 100, mu: 0.0, sigma: 0.1, seed: 42); + int count = 50000; + double sum = 0; + var cma = new Cma(); + + for (int i = 0; i < count; i++) + { + double value = gbm.Next().Close; + sum += value; + cma.Update(new TValue(DateTime.UtcNow, value)); + } + + double naiveMean = sum / count; + double welfordMean = cma.Last.Value; + + // Both should be very close + Assert.True( + Math.Abs(naiveMean - welfordMean) < 1e-8, + $"Naive={naiveMean:G17}, Welford={welfordMean:G17}, Diff={Math.Abs(naiveMean - welfordMean):G17}"); + + _output.WriteLine($"CMA long sequence: Naive={naiveMean:G10}, Welford={welfordMean:G10}"); + } + + [Fact] + public void Validate_KnownSequence_ArithmeticProgression() + { + // Arithmetic progression: 1, 2, 3, ..., n + // CMA at each point: 1, 1.5, 2, 2.5, 3, ... + // Formula: CMA_n = (n+1)/2 + + var cma = new Cma(); + + for (int n = 1; n <= 100; n++) + { + cma.Update(new TValue(DateTime.UtcNow, n)); + double expected = (n + 1.0) / 2.0; + Assert.Equal(expected, cma.Last.Value, 1e-10); + } + + _output.WriteLine("CMA validated for arithmetic progression"); + } + + [Fact] + public void Validate_KnownSequence_GeometricProgression() + { + // Geometric progression: r, r^2, r^3, ..., r^n + // Sum = r * (r^n - 1) / (r - 1) + // CMA = Sum / n + + double r = 1.1; + var cma = new Cma(); + + for (int n = 1; n <= 50; n++) + { + double value = Math.Pow(r, n); + cma.Update(new TValue(DateTime.UtcNow, value)); + + // Sum of geometric series: a * (r^n - 1) / (r - 1) where a = r + double sum = r * (Math.Pow(r, n) - 1) / (r - 1); + double expected = sum / n; + + Assert.Equal(expected, cma.Last.Value, 1e-9); + } + + _output.WriteLine("CMA validated for geometric progression"); + } + + [Fact] + public void Validate_ConstantSequence() + { + // CMA of constant sequence should be the constant + double constant = 42.5; + var cma = new Cma(); + + for (int i = 0; i < 10000; i++) + { + cma.Update(new TValue(DateTime.UtcNow, constant)); + } + + Assert.Equal(constant, cma.Last.Value, 1e-10); + _output.WriteLine("CMA validated for constant sequence"); + } + + [Fact] + public void Validate_AllModes_Consistency() + { + // Verify all three calculation modes produce identical results + var sourceData = _testData.RawData.ToArray(); + + // Mode 1: TSeries Batch + var cma1 = new Cma(); + var batchResult = cma1.Update(_testData.Data); + + // Mode 2: Streaming + var cma2 = new Cma(); + var streamingResults = new List(); + foreach (var item in _testData.Data) + { + streamingResults.Add(cma2.Update(item).Value); + } + + // Mode 3: Span + var spanOutput = new double[sourceData.Length]; + Cma.Batch(sourceData.AsSpan(), spanOutput.AsSpan()); + + // Compare all three + for (int i = 0; i < sourceData.Length; i++) + { + double batchVal = batchResult[i].Value; + double streamVal = streamingResults[i]; + double spanVal = spanOutput[i]; + + Assert.Equal(batchVal, streamVal, 1e-10); + Assert.Equal(batchVal, spanVal, 1e-10); + } + + _output.WriteLine("All CMA calculation modes produce consistent results"); + } +} \ No newline at end of file diff --git a/lib/statistics/cma/Cma.cs b/lib/statistics/cma/Cma.cs new file mode 100644 index 00000000..3a63b568 --- /dev/null +++ b/lib/statistics/cma/Cma.cs @@ -0,0 +1,272 @@ +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; + +namespace QuanTAlib; + +/// +/// CMA: Cumulative Moving Average (Running Average / Cumulative Mean) +/// +/// +/// CMA calculates the arithmetic mean of ALL data points seen so far, not just a fixed window. +/// Uses Welford's algorithm with FMA (Fused Multiply-Add) for maximum numerical precision. +/// +/// Calculation: +/// M_n = M_(n-1) + ± * (x_n - M_(n-1)) where ± = 1/n +/// +/// Implemented using FMA for single-rounding precision: +/// mean = FusedMultiplyAdd(alpha, delta, mean) +/// +/// This is equivalent to: +/// M_n = ((n-1) * M_(n-1) + x_n) / n +/// +/// Key Features: +/// - Zero window: includes ALL historical data with equal weight +/// - O(1) time complexity per update +/// - Maximum precision: FMA avoids intermediate rounding of alpha*delta +/// - Numerically stable: avoids overflow from summing large sequences +/// - No buffer required: only stores count and mean +/// +/// IsHot: +/// Always true after the first value (no warmup period needed). +/// +[SkipLocalsInit] +public sealed class Cma : AbstractBase +{ + [StructLayout(LayoutKind.Auto)] + private record struct State(double Mean, long Count, double LastValidValue); + private State _state; + private State _p_state; + + private readonly TValuePublishedHandler _handler; + + /// + /// Creates a new CMA indicator instance. + /// No period parameter required since CMA averages all values. + /// + public Cma() + { + Name = "Cma"; + WarmupPeriod = 1; + _handler = Handle; + } + + /// + /// Creates CMA with a source to subscribe to. + /// + /// Source to subscribe to + public Cma(ITValuePublisher source) : this() + { + source.Pub += _handler; + } + + /// + /// Creates CMA with a TSeries source to prime from and subscribe to. + /// + /// TSeries source + public Cma(TSeries source) : this() + { + Prime(source.Values); + if (source.Count > 0) + { + Last = new TValue(source.LastTime, Last.Value); + } + source.Pub += _handler; + } + + private void Handle(object? sender, in TValueEventArgs e) => Update(e.Value, e.IsNew); + + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Mode B: Streaming (Stateful) + ///////////////////////////////////////////////////////////////////////////////////////////////// + + /// + /// True if the CMA has enough data to produce valid results. + /// CMA is "hot" after the first value since no warmup is needed. + /// + public override bool IsHot => _state.Count > 0; + + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Mode C: Priming (The Bridge) + ///////////////////////////////////////////////////////////////////////////////////////////////// + + /// + /// Initializes the indicator state using the provided history. + /// + /// Historical data + /// Time interval between values (not used for CMA) + public override void Prime(ReadOnlySpan source, TimeSpan? step = null) + { + if (source.Length == 0) return; + + // Reset state + _state = default; + _p_state = default; + + // Find first valid value to seed lastValid + for (int i = 0; i < source.Length; i++) + { + if (double.IsFinite(source[i])) + { + _state.LastValidValue = source[i]; + break; + } + } + + // Process all values using Welford's algorithm with FMA + for (int i = 0; i < source.Length; i++) + { + double val = GetValidValue(source[i]); + _state.Count++; + double alpha = 1.0 / _state.Count; + double delta = val - _state.Mean; + _state.Mean = Math.FusedMultiplyAdd(alpha, delta, _state.Mean); + } + + Last = new TValue(DateTime.MinValue, _state.Mean); + _p_state = _state; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private double GetValidValue(double input) + { + if (double.IsFinite(input)) + { + _state.LastValidValue = input; + return input; + } + return _state.LastValidValue; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public override TValue Update(TValue input, bool isNew = true) + { + if (isNew) + { + _p_state = _state; + double val = GetValidValue(input.Value); + _state.Count++; + double alpha = 1.0 / _state.Count; + double delta = val - _state.Mean; + _state.Mean = Math.FusedMultiplyAdd(alpha, delta, _state.Mean); + } + else + { + _state = _p_state; + double val = GetValidValue(input.Value); + _state.Count++; + double alpha = 1.0 / _state.Count; + double delta = val - _state.Mean; + _state.Mean = Math.FusedMultiplyAdd(alpha, delta, _state.Mean); + } + + Last = new TValue(input.Time, _state.Mean); + 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); + source.Times.CopyTo(tSpan); + + Prime(source.Values); + + Last = new TValue(tSpan[len - 1], vSpan[len - 1]); + return new TSeries(t, v); + } + + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Mode A: Batch (Stateless) + ///////////////////////////////////////////////////////////////////////////////////////////////// + + /// + /// Calculates CMA for the entire series using a new instance. + /// + /// Input series + /// CMA series + public static TSeries Batch(TSeries source) + { + var cma = new Cma(); + return cma.Update(source); + } + + /// + /// Calculates CMA in-place, writing results to pre-allocated output span. + /// Zero-allocation method for maximum performance. + /// Uses Welford's algorithm for numerical stability. + /// + /// Input values + /// Output span (must be same length as source) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void Batch(ReadOnlySpan source, Span output) + { + 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; + + double mean = 0; + double lastValid = double.NaN; + + // Find first valid value to seed lastValid + for (int k = 0; k < len; k++) + { + if (double.IsFinite(source[k])) + { + lastValid = source[k]; + break; + } + } + + // Welford's algorithm for running mean with FMA + for (int i = 0; i < len; i++) + { + double val = source[i]; + if (double.IsFinite(val)) + lastValid = val; + else + val = lastValid; + + // M_n = M_(n-1) + alpha * delta using FMA for single-rounding precision + double alpha = 1.0 / (i + 1); + double delta = val - mean; + mean = Math.FusedMultiplyAdd(alpha, delta, mean); + output[i] = mean; + } + } + + /// + /// Runs a batch calculation on history and returns + /// a "Hot" Cma instance ready to process the next tick immediately. + /// + /// Historical time series + /// A tuple containing the full calculation results and the hot indicator instance + public static (TSeries Results, Cma Indicator) Calculate(TSeries source) + { + var cma = new Cma(); + TSeries results = cma.Update(source); + return (results, cma); + } + + /// + /// Resets the CMA state. + /// + public override void Reset() + { + _state = default; + _p_state = default; + Last = default; + } +} diff --git a/lib/statistics/cma/Cma.md b/lib/statistics/cma/Cma.md new file mode 100644 index 00000000..6f6bcf15 --- /dev/null +++ b/lib/statistics/cma/Cma.md @@ -0,0 +1,97 @@ +# CMA: Cumulative Moving Average + +> "The running average that never forgets. Every single tick you've ever fed it? Still in there, affecting the result. It's like the elephant of technical indicators." + +The Cumulative Moving Average (CMA) calculates the arithmetic mean of ALL data points seen so far, not just a fixed window. Unlike SMA or EMA which use a sliding window, CMA treats every historical value with equal weight. As the sample size grows, each new value has diminishing impact on the average. + +## Historical Context + +The concept of a running mean is fundamental to statistics and was formalized by B. P. Welford in 1962 for numerically stable computation. Donald Knuth popularized it in *The Art of Computer Programming*. While not a traditional trading indicator, CMA is essential for scenarios requiring the true average of all observed data: calculating session VWAP from scratch, averaging tick counts, or computing lifetime average fill prices. + +## Architecture & Physics + +The naive approach (sum all values, divide by count) works for small datasets but fails at scale. After millions of ticks, the running sum can overflow or lose precision. + +### Welford's Algorithm with FMA + +QuanTAlib uses Welford's numerically stable update, enhanced with Fused Multiply-Add (FMA) for maximum precision: + +$$ M_n = M_{n-1} + \alpha \cdot (x_n - M_{n-1}) \quad \text{where } \alpha = \frac{1}{n} $$ + +Implemented as: + +```csharp +double alpha = 1.0 / n; +double delta = x - mean; +mean = Math.FusedMultiplyAdd(alpha, delta, mean); +``` + +This formulation: + +1. Keeps intermediate values near the scale of the actual mean (no overflow) +2. Requires only O(1) memory (just count and mean) +3. Achieves O(1) time complexity per update +4. Uses FMA for single-rounding precision (avoids rounding `alpha * delta` before adding to `mean`) +5. Is mathematically equivalent to $M_n = \frac{(n-1) \cdot M_{n-1} + x_n}{n}$ + +### Why Not Just Sum? + +Consider averaging 10 million tick prices around 50,000 (a futures contract). The naive sum exceeds $5 \times 10^{11}$, approaching the precision limits of `double`. Welford's algorithm keeps the working value around 50,000 throughout, maintaining full precision. + +### The Diminishing Return Problem + +As $n$ grows large, each new value contributes only $\frac{1}{n}$ to the mean. After 1 million samples, a new tick moves the average by roughly 0.0001% of the difference from the current mean. This is mathematically correct but may not be what traders want for responsiveness (use EMA or SMA for that). + +## Mathematical Foundation + +### 1. Incremental Update (Welford) + +$$ M_n = M_{n-1} + \frac{x_n - M_{n-1}}{n} $$ + +Where: + +- $M_n$ = cumulative mean after $n$ values +- $M_{n-1}$ = previous cumulative mean +- $x_n$ = new value +- $n$ = total count of values + +### 2. Algebraic Equivalence + +$$ M_n = \frac{1}{n} \sum_{i=1}^{n} x_i = \frac{(n-1) \cdot M_{n-1} + x_n}{n} $$ + +## Performance Profile + +| Metric | Score | Notes | +| :--- | :--- | :--- | +| **Throughput** | ~5 ns/bar | Single division per update. | +| **Allocations** | 0 | Zero-allocation in hot paths. | +| **Complexity** | O(1) | Constant time regardless of history length. | +| **Accuracy** | 10 | Welford's algorithm ensures numerical stability. | +| **Timeliness** | 1 | Maximum lag; every historical value affects output. | +| **Overshoot** | 0 | Never overshoots the input data range. | +| **Smoothness** | 10 | Extremely smooth as $n$ grows (almost constant). | + +## Validation + +| Library | Status | Notes | +| :--- | :--- | :--- | +| **TA-Lib** | N/A | No CMA function. | +| **Skender** | N/A | No CMA function. | +| **Tulip** | N/A | No CMA function. | +| **Mathematical** | ✅ | Validated against known formulas. | + +CMA is a fundamental statistical operation rather than a standard TA library indicator. QuanTAlib validates against mathematical proofs: arithmetic progressions, geometric series, and direct sum/count calculations. + +## Use Cases + +1. **Session VWAP**: Calculate volume-weighted average price from session start +2. **Lifetime Averages**: Average fill price across all trades +3. **Quality Metrics**: Average latency, slippage, or fill rate over time +4. **Baseline Comparison**: Compare current price to "all-time average" + +## Common Pitfalls + +1. **Responsiveness**: CMA becomes nearly unresponsive after many values. For a reactive average, use SMA or EMA instead. +2. **Memory of Bad Data**: A single extreme outlier early in the stream permanently affects the average. Consider filtering before feeding CMA. +3. **No Period Parameter**: Unlike SMA/EMA, CMA has no period. It always includes all data. This is by design. +4. **Session Resets**: If you need per-session averages, call `Reset()` at session boundaries. \ No newline at end of file diff --git a/lib/trends/_index.md b/lib/trends/_index.md index 4e9eeee0..db53ac22 100644 --- a/lib/trends/_index.md +++ b/lib/trends/_index.md @@ -74,4 +74,4 @@ Trend indicators are the bread and butter of technical analysis—and often just | YZVAMA | Yang-Zhang Volatility Adjusted MA | | | ZLDEMA | Zero-Lag Double Exponential MA | | | ZLEMA | Zero-Lag Exponential MA | | -| ZLTEMA | Zero-Lag Triple Exponential MA | | +| ZLTEMA | Zero-Lag Triple Exponential MA | | \ No newline at end of file diff --git a/lib/trends/bessel/Bessel.cs b/lib/trends/bessel/Bessel.cs index de716300..2bef9b3c 100644 --- a/lib/trends/bessel/Bessel.cs +++ b/lib/trends/bessel/Bessel.cs @@ -135,7 +135,8 @@ public sealed class Bessel : AbstractBase, IDisposable double filt = _state.Count < 3 ? val - : (_c1 * val) + (_c2 * _state.F1) + (_c3 * _state.F2); + : Math.FusedMultiplyAdd(_c3, _state.F2, + Math.FusedMultiplyAdd(_c2, _state.F1, _c1 * val)); _state.F2 = _state.F1; _state.F1 = filt; @@ -184,7 +185,8 @@ public sealed class Bessel : AbstractBase, IDisposable double filt = _state.Count < 3 ? val - : (_c1 * val) + (_c2 * _state.F1) + (_c3 * _state.F2); + : Math.FusedMultiplyAdd(_c3, _state.F2, + Math.FusedMultiplyAdd(_c2, _state.F1, _c1 * val)); _state.F2 = _state.F1; _state.F1 = filt; @@ -276,7 +278,8 @@ public sealed class Bessel : AbstractBase, IDisposable double filt = state.Count < 3 ? val - : (c1 * val) + (c2 * state.F1) + (c3 * state.F2); + : Math.FusedMultiplyAdd(c3, state.F2, + Math.FusedMultiplyAdd(c2, state.F1, c1 * val)); state.F2 = state.F1; state.F1 = filt; diff --git a/lib/trends/butter/Butter.cs b/lib/trends/butter/Butter.cs index ee8ea65e..ba984108 100644 --- a/lib/trends/butter/Butter.cs +++ b/lib/trends/butter/Butter.cs @@ -113,9 +113,14 @@ public sealed class Butter : AbstractBase } double x = input.Value; + // IIR: y = (b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2) * invA0 + // Using chained FMA for precision double y = _state.Count < 2 ? x - : (_b0 * x + _b1 * _state.X1 + _b2 * _state.X2 - _a1 * _state.Y1 - _a2 * _state.Y2) * _invA0; + : Math.FusedMultiplyAdd(-_a2, _state.Y2, + Math.FusedMultiplyAdd(-_a1, _state.Y1, + Math.FusedMultiplyAdd(_b2, _state.X2, + Math.FusedMultiplyAdd(_b1, _state.X1, _b0 * x)))) * _invA0; // Update state _state.X2 = _state.X1; @@ -185,9 +190,14 @@ public sealed class Butter : AbstractBase destination[i] = i > 0 ? destination[i - 1] : initialLast; continue; } + // IIR: y = (b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2) * invA0 + // Using chained FMA for precision double y = i < 2 ? x - : (b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2) * invA0; + : Math.FusedMultiplyAdd(-a2, y2, + Math.FusedMultiplyAdd(-a1, y1, + Math.FusedMultiplyAdd(b2, x2, + Math.FusedMultiplyAdd(b1, x1, b0 * x)))) * invA0; x2 = x1; x1 = x; diff --git a/lib/trends/htit/Htit.cs b/lib/trends/htit/Htit.cs index 7724eec6..8374f0df 100644 --- a/lib/trends/htit/Htit.cs +++ b/lib/trends/htit/Htit.cs @@ -161,17 +161,17 @@ public sealed class Htit : AbstractBase double i2_val = i1 - jQ; double q2_val = q1 + jI; - // Smooth i2, q2 - _state.I2 = 0.2 * i2_val + 0.8 * _p_state.I2; - _state.Q2 = 0.2 * q2_val + 0.8 * _p_state.Q2; + // Smooth i2, q2 (using FMA for precision) + _state.I2 = Math.FusedMultiplyAdd(0.2, i2_val, 0.8 * _p_state.I2); + _state.Q2 = Math.FusedMultiplyAdd(0.2, q2_val, 0.8 * _p_state.Q2); // 6. Homodyne Discriminator - double re_val = (_state.I2 * _p_state.I2) + (_state.Q2 * _p_state.Q2); - double im_val = (_state.I2 * _p_state.Q2) - (_state.Q2 * _p_state.I2); + double re_val = Math.FusedMultiplyAdd(_state.I2, _p_state.I2, _state.Q2 * _p_state.Q2); + double im_val = Math.FusedMultiplyAdd(_state.I2, _p_state.Q2, -_state.Q2 * _p_state.I2); - // Smooth re, im - _state.Re = 0.2 * re_val + 0.8 * _p_state.Re; - _state.Im = 0.2 * im_val + 0.8 * _p_state.Im; + // Smooth re, im (using FMA) + _state.Re = Math.FusedMultiplyAdd(0.2, re_val, 0.8 * _p_state.Re); + _state.Im = Math.FusedMultiplyAdd(0.2, im_val, 0.8 * _p_state.Im); // 7. Calculate Period double angle = Math.Atan2(_state.Im, _state.Re); @@ -190,9 +190,9 @@ public sealed class Htit : AbstractBase if (period < 6) period = 6; if (period > 50) period = 50; - // Smooth the period - _state.Period = 0.2 * period + 0.8 * prevPeriod; - _state.SmoothPeriod = 0.33 * _state.Period + 0.67 * _p_state.SmoothPeriod; + // Smooth the period (using FMA) + _state.Period = Math.FusedMultiplyAdd(0.2, period, 0.8 * prevPeriod); + _state.SmoothPeriod = Math.FusedMultiplyAdd(0.33, _state.Period, 0.67 * _p_state.SmoothPeriod); // 8. Instantaneous Trend int dcPeriods = (int)(double.IsNaN(_state.SmoothPeriod) ? 0 : _state.SmoothPeriod + 0.5); @@ -376,15 +376,15 @@ public sealed class Htit : AbstractBase double i2_val = i1 - jQ; double q2_val = q1 + jI; - i2 = 0.2 * i2_val + 0.8 * p_i2; - q2 = 0.2 * q2_val + 0.8 * p_q2; + i2 = Math.FusedMultiplyAdd(0.2, i2_val, 0.8 * p_i2); + q2 = Math.FusedMultiplyAdd(0.2, q2_val, 0.8 * p_q2); // 6. Homodyne Discriminator - double re_val = (i2 * p_i2) + (q2 * p_q2); - double im_val = (i2 * p_q2) - (q2 * p_i2); + double re_val = Math.FusedMultiplyAdd(i2, p_i2, q2 * p_q2); + double im_val = Math.FusedMultiplyAdd(i2, p_q2, -q2 * p_i2); - re = 0.2 * re_val + 0.8 * p_re; - im = 0.2 * im_val + 0.8 * p_im; + re = Math.FusedMultiplyAdd(0.2, re_val, 0.8 * p_re); + im = Math.FusedMultiplyAdd(0.2, im_val, 0.8 * p_im); // 7. Calculate Period double angle = Math.Atan2(im, re); @@ -402,8 +402,8 @@ public sealed class Htit : AbstractBase if (newPeriod < 6) newPeriod = 6; if (newPeriod > 50) newPeriod = 50; - period = 0.2 * newPeriod + 0.8 * p_period; - smoothPeriod = 0.33 * period + 0.67 * p_smoothPeriod; + period = Math.FusedMultiplyAdd(0.2, newPeriod, 0.8 * p_period); + smoothPeriod = Math.FusedMultiplyAdd(0.33, period, 0.67 * p_smoothPeriod); // 8. Instantaneous Trend double safeSmooth = double.IsNaN(smoothPeriod) ? 0 : smoothPeriod; diff --git a/lib/trends/jma/Jma.cs b/lib/trends/jma/Jma.cs index 4954f17b..23c0c7cd 100644 --- a/lib/trends/jma/Jma.cs +++ b/lib/trends/jma/Jma.cs @@ -216,13 +216,17 @@ public sealed class Jma : AbstractBase prevJma = value; double alpha = Math.Exp(_logLengthDivider * d); + double decay = 1.0 - alpha; double alpha2 = alpha * alpha; - double c0 = (1.0 - alpha) * value + alpha * _state.LastC0; - double c8 = (value - c0) * (1.0 - _lengthDivider) + _lengthDivider * _state.LastC8; - double a8 = (_phaseParam * c8 + c0 - prevJma) * - (alpha * (-2.0) + alpha2 + 1.0) + - alpha2 * _state.LastA8; + // EMA smoothing: c0 = decay * value + alpha * LastC0 + double c0 = Math.FusedMultiplyAdd(_state.LastC0, alpha, decay * value); + // EMA smoothing: c8 = (value - c0) * (1 - lengthDivider) + lengthDivider * LastC8 + double lengthDecay = 1.0 - _lengthDivider; + double c8 = Math.FusedMultiplyAdd(_state.LastC8, _lengthDivider, lengthDecay * (value - c0)); + // IIR filter: a8 = (phase * c8 + c0 - prevJma) * coef + alpha2 * LastA8 + double coef = Math.FusedMultiplyAdd(alpha, -2.0, alpha2 + 1.0); + double a8 = Math.FusedMultiplyAdd(_state.LastA8, alpha2, Math.FusedMultiplyAdd(_phaseParam, c8, c0 - prevJma) * coef); double jma = prevJma + a8; diff --git a/lib/trends/mama/Mama.cs b/lib/trends/mama/Mama.cs index 701f0cb6..97531e82 100644 --- a/lib/trends/mama/Mama.cs +++ b/lib/trends/mama/Mama.cs @@ -170,17 +170,17 @@ public sealed class Mama : AbstractBase double i2_val = i1 - jQ; double q2_val = q1 + jI; - // Smooth i2, q2 - _state.I2 = SmoothCoef * i2_val + SmoothPrev * _p_state.I2; - _state.Q2 = SmoothCoef * q2_val + SmoothPrev * _p_state.Q2; + // Smooth i2, q2 (using FMA for precision) + _state.I2 = Math.FusedMultiplyAdd(SmoothCoef, i2_val, SmoothPrev * _p_state.I2); + _state.Q2 = Math.FusedMultiplyAdd(SmoothCoef, q2_val, SmoothPrev * _p_state.Q2); // Homodyne discriminator - double re_val = (_state.I2 * _p_state.I2) + (_state.Q2 * _p_state.Q2); - double im_val = (_state.I2 * _p_state.Q2) - (_state.Q2 * _p_state.I2); + double re_val = Math.FusedMultiplyAdd(_state.I2, _p_state.I2, _state.Q2 * _p_state.Q2); + double im_val = Math.FusedMultiplyAdd(_state.I2, _p_state.Q2, -_state.Q2 * _p_state.I2); - // Smooth re, im - _state.Re = SmoothCoef * re_val + SmoothPrev * _p_state.Re; - _state.Im = SmoothCoef * im_val + SmoothPrev * _p_state.Im; + // Smooth re, im (using FMA) + _state.Re = Math.FusedMultiplyAdd(SmoothCoef, re_val, SmoothPrev * _p_state.Re); + _state.Im = Math.FusedMultiplyAdd(SmoothCoef, im_val, SmoothPrev * _p_state.Im); // Calculate Period double angle = Math.Atan2(_state.Im, _state.Re); @@ -198,8 +198,8 @@ public sealed class Mama : AbstractBase if (period < MinPeriod) period = MinPeriod; if (period > MaxPeriod) period = MaxPeriod; - // Smooth Period - _state.Period = SmoothCoef * period + SmoothPrev * _p_state.Period; + // Smooth Period (using FMA) + _state.Period = Math.FusedMultiplyAdd(SmoothCoef, period, SmoothPrev * _p_state.Period); // Phase calculation _state.Phase = Math.Atan2(q1, i1); @@ -380,17 +380,17 @@ public sealed class Mama : AbstractBase double i2_val = i1 - jQ; double q2_val = q1 + jI; - // Smooth i2, q2 - i2 = SmoothCoef * i2_val + SmoothPrev * p_i2; - q2 = SmoothCoef * q2_val + SmoothPrev * p_q2; + // Smooth i2, q2 (using FMA for precision) + i2 = Math.FusedMultiplyAdd(SmoothCoef, i2_val, SmoothPrev * p_i2); + q2 = Math.FusedMultiplyAdd(SmoothCoef, q2_val, SmoothPrev * p_q2); // Homodyne discriminator - double re_val = (i2 * p_i2) + (q2 * p_q2); - double im_val = (i2 * p_q2) - (q2 * p_i2); + double re_val = Math.FusedMultiplyAdd(i2, p_i2, q2 * p_q2); + double im_val = Math.FusedMultiplyAdd(i2, p_q2, -q2 * p_i2); - // Smooth re, im - re = SmoothCoef * re_val + SmoothPrev * p_re; - im = SmoothCoef * im_val + SmoothPrev * p_im; + // Smooth re, im (using FMA) + re = Math.FusedMultiplyAdd(SmoothCoef, re_val, SmoothPrev * p_re); + im = Math.FusedMultiplyAdd(SmoothCoef, im_val, SmoothPrev * p_im); // Calculate Period double angle = Math.Atan2(im, re); @@ -408,8 +408,8 @@ public sealed class Mama : AbstractBase if (newPeriod < MinPeriod) newPeriod = MinPeriod; if (newPeriod > MaxPeriod) newPeriod = MaxPeriod; - // Smooth Period - period = SmoothCoef * newPeriod + SmoothPrev * p_period; + // Smooth Period (using FMA) + period = Math.FusedMultiplyAdd(SmoothCoef, newPeriod, SmoothPrev * p_period); // Phase calculation double phase = Math.Atan2(q1, i1); diff --git a/lib/trends/ssf/Ssf.cs b/lib/trends/ssf/Ssf.cs index 7931c878..97325d26 100644 --- a/lib/trends/ssf/Ssf.cs +++ b/lib/trends/ssf/Ssf.cs @@ -127,7 +127,8 @@ public sealed class Ssf : AbstractBase double ssf = (_state.Count < 4) ? val - : (_c1 * (val + _state.PrevInput) * 0.5) + (_c2 * _state.Ssf1) + (_c3 * _state.Ssf2); + : Math.FusedMultiplyAdd(_c3, _state.Ssf2, + Math.FusedMultiplyAdd(_c2, _state.Ssf1, _c1 * (val + _state.PrevInput) * 0.5)); _state.Ssf2 = _state.Ssf1; _state.Ssf1 = ssf; @@ -177,7 +178,8 @@ public sealed class Ssf : AbstractBase double ssf = (_state.Count < 4) ? val - : (_c1 * (val + _state.PrevInput) * 0.5) + (_c2 * _state.Ssf1) + (_c3 * _state.Ssf2); + : Math.FusedMultiplyAdd(_c3, _state.Ssf2, + Math.FusedMultiplyAdd(_c2, _state.Ssf1, _c1 * (val + _state.PrevInput) * 0.5)); _state.Ssf2 = _state.Ssf1; _state.Ssf1 = ssf; @@ -268,7 +270,8 @@ public sealed class Ssf : AbstractBase double ssf = (state.Count < 4) ? val - : (c1 * (val + state.PrevInput) * 0.5) + (c2 * state.Ssf1) + (c3 * state.Ssf2); + : Math.FusedMultiplyAdd(c3, state.Ssf2, + Math.FusedMultiplyAdd(c2, state.Ssf1, c1 * (val + state.PrevInput) * 0.5)); state.Ssf2 = state.Ssf1; state.Ssf1 = ssf; diff --git a/lib/trends/vidya/Vidya.cs b/lib/trends/vidya/Vidya.cs index c6a4cc44..449b30d6 100644 --- a/lib/trends/vidya/Vidya.cs +++ b/lib/trends/vidya/Vidya.cs @@ -131,7 +131,8 @@ public sealed class Vidya : AbstractBase, IDisposable } double dynamicAlpha = _alpha * vi; - _state.CurrentVidya = dynamicAlpha * price + (1.0 - dynamicAlpha) * _state.LastVidya; + double dynamicDecay = 1.0 - dynamicAlpha; + _state.CurrentVidya = Math.FusedMultiplyAdd(_state.LastVidya, dynamicDecay, dynamicAlpha * price); _state.CurrentClose = price; Last = new TValue(input.Time, _state.CurrentVidya); @@ -239,7 +240,8 @@ public sealed class Vidya : AbstractBase, IDisposable } double dynamicAlpha = _alpha * vi; - double currentVidya = dynamicAlpha * price + (1.0 - dynamicAlpha) * lastVidya; + double dynamicDecay = 1.0 - dynamicAlpha; + double currentVidya = Math.FusedMultiplyAdd(lastVidya, dynamicDecay, dynamicAlpha * price); _state.CurrentVidya = currentVidya; _state.CurrentClose = price; @@ -337,7 +339,8 @@ public sealed class Vidya : AbstractBase, IDisposable } double dynamicAlpha = alpha * vi; - double currentVidya = dynamicAlpha * price + (1.0 - dynamicAlpha) * lastVidya; + double dynamicDecay = 1.0 - dynamicAlpha; + double currentVidya = Math.FusedMultiplyAdd(lastVidya, dynamicDecay, dynamicAlpha * price); output[i] = currentVidya; diff --git a/lib/volume/adosc/Adosc.cs b/lib/volume/adosc/Adosc.cs index f88fabd7..4392275a 100644 --- a/lib/volume/adosc/Adosc.cs +++ b/lib/volume/adosc/Adosc.cs @@ -148,8 +148,13 @@ public sealed class Adosc : ITValuePublisher return adosc.Update(source); } + // EMA compensator threshold (same as in Ema.cs) + private const double COMPENSATOR_THRESHOLD = 1e-10; + /// - /// Calculates ADOSC for the entire span. + /// Calculates ADOSC for the entire span using a single-pass algorithm. + /// Zero allocation for maximum performance. + /// Uses compensator pattern from EMA for proper early-stage bias correction. /// /// High prices /// Low prices @@ -161,18 +166,101 @@ public sealed class Adosc : ITValuePublisher [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Calculate(ReadOnlySpan high, ReadOnlySpan low, ReadOnlySpan close, ReadOnlySpan volume, Span output, int fastPeriod = 3, int slowPeriod = 10) { - if (high.Length != output.Length) - throw new ArgumentException("Source and output spans must be of the same length.", nameof(output)); + if (high.Length != low.Length || high.Length != close.Length || + high.Length != volume.Length || high.Length != output.Length) + throw new ArgumentException("All spans must be of the same length.", nameof(output)); - Span adl = high.Length <= 1024 ? stackalloc double[high.Length] : new double[high.Length]; - Adl.Calculate(high, low, close, volume, adl); + if (fastPeriod <= 0) + throw new ArgumentException("Fast period must be greater than 0", nameof(fastPeriod)); + if (slowPeriod <= 0) + throw new ArgumentException("Slow period must be greater than 0", nameof(slowPeriod)); - Span fastEma = high.Length <= 1024 ? stackalloc double[high.Length] : new double[high.Length]; - Span slowEma = high.Length <= 1024 ? stackalloc double[high.Length] : new double[high.Length]; + int len = high.Length; + if (len == 0) return; - Ema.Batch(adl, fastEma, fastPeriod); - Ema.Batch(adl, slowEma, slowPeriod); + // EMA parameters (same formula as Ema.cs: alpha = 2 / (period + 1)) + double alphaFast = 2.0 / (fastPeriod + 1); + double alphaSlow = 2.0 / (slowPeriod + 1); + double decayFast = 1.0 - alphaFast; + double decaySlow = 1.0 - alphaSlow; - SimdExtensions.Subtract(fastEma, slowEma, output); + // State variables (no heap allocations) + double adl = 0; + double emaFast = 0; + double emaSlow = 0; + double eFast = 1.0; // Compensation factor for fast EMA (starts at 1, decays toward 0) + double eSlow = 1.0; // Compensation factor for slow EMA + bool fastCompensated = false; + bool slowCompensated = false; + + // Single pass: compute ADL, both EMAs, and output in one loop + for (int i = 0; i < len; i++) + { + double h = high[i]; + double l = low[i]; + double c = close[i]; + double vol = volume[i]; + + // 1. Compute Money Flow Multiplier and Volume + double hl = h - l; + double mfm = 0; + if (hl > double.Epsilon) + { + mfm = ((c - l) - (h - c)) / hl; + } + double mfv = mfm * vol; + + // 2. Update ADL (cumulative) + adl += mfv; + + // 3. Update Fast EMA with FMA (same pattern as Ema.cs Compute method) + // state.Ema = Math.FusedMultiplyAdd(state.Ema, decay, alpha * input) + emaFast = Math.FusedMultiplyAdd(emaFast, decayFast, alphaFast * adl); + + // 4. Update Slow EMA with FMA + emaSlow = Math.FusedMultiplyAdd(emaSlow, decaySlow, alphaSlow * adl); + + // 5. Compute compensated EMA values (same logic as Ema.cs Compute method) + // Compensator decays: e *= decay, then result = ema / (1 - e) until e <= threshold + double fastValue, slowValue; + + if (!fastCompensated) + { + eFast *= decayFast; + if (eFast <= COMPENSATOR_THRESHOLD) + { + fastCompensated = true; + fastValue = emaFast; + } + else + { + fastValue = emaFast / (1.0 - eFast); + } + } + else + { + fastValue = emaFast; + } + + if (!slowCompensated) + { + eSlow *= decaySlow; + if (eSlow <= COMPENSATOR_THRESHOLD) + { + slowCompensated = true; + slowValue = emaSlow; + } + else + { + slowValue = emaSlow / (1.0 - eSlow); + } + } + else + { + slowValue = emaSlow; + } + + output[i] = fastValue - slowValue; + } } } diff --git a/perf/Benchmark.cs b/perf/Benchmark.cs index 64c8d9bd..178b4c0c 100644 --- a/perf/Benchmark.cs +++ b/perf/Benchmark.cs @@ -100,7 +100,7 @@ public class IndicatorBenchmarks { _quotes.Add(new Quote { - Date = new DateTime(_closeTseries.Times[i]), + Date = new DateTime(_closeTseries.Times[i], DateTimeKind.Utc), Open = (decimal)bars.Open.Values[i], High = (decimal)bars.High.Values[i], Low = (decimal)bars.Low.Values[i], @@ -115,7 +115,7 @@ public class IndicatorBenchmarks { _ooplesData.Add(new TickerData { - Date = new DateTime(_closeTseries.Times[i]), + Date = new DateTime(_closeTseries.Times[i], DateTimeKind.Utc), Open = bars.Open.Values[i], High = bars.High.Values[i], Low = bars.Low.Values[i],