using System.Buffers;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
namespace QuanTAlib;
///
/// PARZEN: Parzen (de la Vallée-Poussin) Window Moving Average
///
///
/// Symmetric FIR filter using the Parzen piecewise cubic window function.
/// The Parzen window is the self-convolution of two Bartlett (triangular) windows
/// at half-length, yielding continuous first and second derivatives and -24 dB/octave
/// sidelobe rolloff. All weights are non-negative.
///
/// Calculation: Precomputed piecewise cubic weights, applied as FIR convolution
/// over sliding window. O(period) per bar.
///
/// Detailed documentation
[SkipLocalsInit]
public sealed class Parzen : AbstractBase
{
private readonly int _period;
private readonly double[] _weights;
private readonly RingBuffer _buffer;
private readonly ITValuePublisher? _source;
private readonly TValuePublishedHandler? _pubHandler;
private bool _isNew = true;
private bool _disposed;
private double _lastValidValue = double.NaN;
private double _p_lastValidValue = double.NaN;
public bool IsNew => _isNew;
public override bool IsHot => _buffer.IsFull;
///
/// Creates PARZEN with specified period.
///
/// Lookback period (>= 2)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public Parzen(int period = 14)
{
if (period < 2)
{
throw new ArgumentException("Period must be at least 2", nameof(period));
}
_period = period;
Name = $"Parzen({_period.ToString(System.Globalization.CultureInfo.InvariantCulture)})";
WarmupPeriod = _period;
_buffer = new RingBuffer(_period);
_weights = new double[_period];
ComputeParzenWeights(_weights, _period);
}
///
/// Creates PARZEN connected to a data source for event-based updates.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public Parzen(ITValuePublisher source, int period = 14) : this(period)
{
_source = source;
_pubHandler = Handle;
_source.Pub += _pubHandler;
}
///
/// Computes Parzen (de la Vallée-Poussin) window weights and normalizes to sum=1.
/// Inner region (|u| <= 0.5): w = 1 - 6u² + 6|u|³
/// Outer region (0.5 < |u| <= 1.0): w = 2(1 - |u|)³
/// All weights are non-negative.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeParzenWeights(Span weights, int period)
{
double halfN = (period - 1) * 0.5;
double wsum = 0.0;
for (int k = 0; k < period; k++)
{
double u = halfN > 0 ? (k - halfN) / halfN : 0.0;
double absU = Math.Abs(u);
double w;
if (absU <= 0.5)
{
// Inner region: cubic spline
w = Math.FusedMultiplyAdd(6.0, absU * absU * absU, 1.0 - 6.0 * absU * absU);
}
else if (absU <= 1.0)
{
// Outer region: cubic taper to zero
double t = 1.0 - absU;
w = 2.0 * t * t * t;
}
else
{
w = 0.0;
}
weights[k] = w;
wsum += w;
}
if (Math.Abs(wsum) > double.Epsilon)
{
double inv = 1.0 / wsum;
for (int k = 0; k < period; k++)
{
weights[k] *= inv;
}
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public override TValue Update(TValue input, bool isNew = true)
{
_isNew = isNew;
return Update(input, isNew, publish: true);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private TValue Update(TValue input, bool isNew, bool publish)
{
if (isNew)
{
_p_lastValidValue = _lastValidValue;
}
else
{
_lastValidValue = _p_lastValidValue;
}
double val = GetValidValue(input.Value);
if (!double.IsFinite(val))
{
Last = new TValue(input.Time, double.NaN);
if (publish) { PubEvent(Last, isNew); }
return Last;
}
if (isNew)
{
_lastValidValue = val;
_buffer.Add(val);
int count = _buffer.Count;
double result;
if (count < _period)
{
result = val;
}
else
{
result = ConvolveFull(_buffer, _weights);
}
Last = new TValue(input.Time, result);
if (publish) { PubEvent(Last, isNew); }
return Last;
}
else
{
_buffer.Snapshot();
double prevLast = _lastValidValue;
double prevPLast = _p_lastValidValue;
_lastValidValue = val;
_buffer.UpdateNewest(val);
int count = _buffer.Count;
double result;
if (count < _period)
{
result = val;
}
else
{
result = ConvolveFull(_buffer, _weights);
}
Last = new TValue(input.Time, result);
_buffer.Restore();
_lastValidValue = prevLast;
_p_lastValidValue = prevPLast;
if (publish) { PubEvent(Last, isNew); }
return Last;
}
}
public override TSeries Update(TSeries source)
{
if (source.Count == 0)
{
return new TSeries([], []);
}
int len = source.Count;
var t = new List(len);
var v = new List(len);
CollectionsMarshal.SetCount(t, len);
CollectionsMarshal.SetCount(v, len);
var tSpan = CollectionsMarshal.AsSpan(t);
var vSpan = CollectionsMarshal.AsSpan(v);
Batch(source.Values, vSpan, _period);
source.Times.CopyTo(tSpan);
Reset();
int startIndex = Math.Max(0, len - _period);
for (int i = startIndex; i < len; i++)
{
Update(source[i], isNew: true, publish: false);
}
return new TSeries(t, v);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void Handle(object? sender, in TValueEventArgs e) => Update(e.Value, e.IsNew);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private double GetValidValue(double input)
{
if (double.IsFinite(input))
{
return input;
}
return double.IsFinite(_lastValidValue) ? _lastValidValue : double.NaN;
}
///
/// FIR convolution using SIMD DotProduct over circular buffer.
/// Weight[0] corresponds to oldest bar, Weight[period-1] to newest.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static double ConvolveFull(RingBuffer buffer, double[] weights)
{
ReadOnlySpan internalBuf = buffer.InternalBuffer;
int head = buffer.StartIndex;
int period = buffer.Capacity;
int part1Len = period - head;
double sum1 = internalBuf.Slice(head, part1Len).DotProduct(weights.AsSpan(0, part1Len));
double sum2 = internalBuf[..head].DotProduct(weights.AsSpan(part1Len));
return sum1 + sum2;
}
public override void Prime(ReadOnlySpan source, TimeSpan? step = null)
{
foreach (var value in source)
{
Update(new TValue(DateTime.MinValue, value));
}
}
public static TSeries Batch(TSeries source, int period = 14)
{
var parzen = new Parzen(period);
return parzen.Update(source);
}
///
/// Calculates Parzen Window MA over a span of values.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Batch(ReadOnlySpan source, Span output, int period = 14, double nanValue = double.NaN)
{
if (period < 2)
{
throw new ArgumentException("Period must be at least 2", nameof(period));
}
if (source.Length != output.Length)
{
throw new ArgumentException("Source and output must have the same length", nameof(output));
}
if (source.Length == 0)
{
return;
}
int len = source.Length;
const int StackallocThreshold = 256;
double[]? weightsRented = period > StackallocThreshold ? ArrayPool.Shared.Rent(period) : null;
Span weights = period <= StackallocThreshold
? stackalloc double[period]
: weightsRented!.AsSpan(0, period);
double[]? ringRented = period > StackallocThreshold ? ArrayPool.Shared.Rent(period) : null;
Span ring = period <= StackallocThreshold
? stackalloc double[period]
: ringRented!.AsSpan(0, period);
double[]? cleanRented = len > StackallocThreshold ? ArrayPool.Shared.Rent(len) : null;
Span clean = len <= StackallocThreshold
? stackalloc double[len]
: cleanRented!.AsSpan(0, len);
ComputeParzenWeights(weights, period);
try
{
double lastValid = nanValue;
for (int i = 0; i < len; i++)
{
double val = source[i];
if (double.IsFinite(val))
{
lastValid = val;
clean[i] = val;
}
else if (double.IsFinite(lastValid))
{
clean[i] = lastValid;
}
else
{
clean[i] = double.NaN;
}
}
int ringIdx = 0;
int count = 0;
for (int i = 0; i < len; i++)
{
double val = clean[i];
ring[ringIdx] = val;
ringIdx++;
if (ringIdx >= period)
{
ringIdx = 0;
}
if (count < period)
{
count++;
}
if (count < period)
{
output[i] = val;
continue;
}
int part1Len = period - ringIdx;
ReadOnlySpan ringRo = ring;
double sum = ringRo.Slice(ringIdx, part1Len).DotProduct(weights.Slice(0, part1Len))
+ ringRo[..ringIdx].DotProduct(weights.Slice(part1Len));
output[i] = sum;
}
}
finally
{
if (weightsRented != null)
{
ArrayPool.Shared.Return(weightsRented);
}
if (ringRented != null)
{
ArrayPool.Shared.Return(ringRented);
}
if (cleanRented != null)
{
ArrayPool.Shared.Return(cleanRented);
}
}
}
public static (TSeries Results, Parzen Indicator) Calculate(TSeries source, int period = 14)
{
var indicator = new Parzen(period);
TSeries results = indicator.Update(source);
return (results, indicator);
}
public override void Reset()
{
_buffer.Clear();
_lastValidValue = double.NaN;
_p_lastValidValue = double.NaN;
Last = default;
}
protected override void Dispose(bool disposing)
{
if (!_disposed)
{
if (disposing && _source != null && _pubHandler != null)
{
_source.Pub -= _pubHandler;
}
_disposed = true;
}
base.Dispose(disposing);
}
}