using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
namespace QuanTAlib;
///
/// PTA: Ehlers Precision Trend Analysis
///
///
///
/// Dual highpass filter bandpass approach for near-zero-lag trend extraction.
/// Applies two 2-pole Butterworth highpass filters with different cutoff periods
/// to the same input, then subtracts: Trend = HP(longPeriod) - HP(shortPeriod).
/// This preserves cyclic components between shortPeriod and longPeriod bars,
/// producing a zero-centered trend indicator with near-zero lag.
///
///
/// Algorithm based on: John F. Ehlers, "Precision Trend Analysis," TASC September 2024.
///
///
/// Complexity: O(1) per bar — two IIR filter evaluations + subtraction.
///
///
[SkipLocalsInit]
public sealed class Pta : AbstractBase
{
// ── HP1 (long-period) coefficients ──
private readonly double _c1L, _c2L, _c3L;
// ── HP2 (short-period) coefficients ──
private readonly double _c1S, _c2S, _c3S;
[StructLayout(LayoutKind.Auto)]
private record struct State
{
// HP1 (long-period highpass)
public double Hp1;
public double Hp1_1;
// HP2 (short-period highpass)
public double Hp2;
public double Hp2_1;
// Source history (shared by both filters)
public double Src1;
public double Src2;
public int Count;
public static State New() => new()
{
Hp1 = 0, Hp1_1 = 0,
Hp2 = 0, Hp2_1 = 0,
Src1 = 0, Src2 = 0,
Count = 0
};
}
private State _state;
private State _p_state;
/// Long-period cutoff for the first highpass filter.
public int LongPeriod { get; }
/// Short-period cutoff for the second highpass filter.
public int ShortPeriod { get; }
///
/// Initializes a new instance of the class.
///
/// Long-period HP cutoff. Default is 250 (~1 year daily).
/// Short-period HP cutoff. Default is 40 (~2 months daily).
///
/// Thrown when longPeriod < 3, shortPeriod < 2, or longPeriod <= shortPeriod.
///
public Pta(int longPeriod = 250, int shortPeriod = 40)
{
ArgumentOutOfRangeException.ThrowIfLessThan(longPeriod, 3, nameof(longPeriod));
ArgumentOutOfRangeException.ThrowIfLessThan(shortPeriod, 2, nameof(shortPeriod));
if (longPeriod <= shortPeriod)
{
throw new ArgumentOutOfRangeException(nameof(longPeriod),
$"longPeriod ({longPeriod}) must be greater than shortPeriod ({shortPeriod}).");
}
LongPeriod = longPeriod;
ShortPeriod = shortPeriod;
// Precompute HP coefficients for long period
ComputeHpCoefficients(longPeriod, out _c1L, out _c2L, out _c3L);
// Precompute HP coefficients for short period
ComputeHpCoefficients(shortPeriod, out _c1S, out _c2S, out _c3S);
Name = $"PTA({longPeriod},{shortPeriod})";
WarmupPeriod = longPeriod;
_state = State.New();
_p_state = _state;
}
///
/// Initializes a new instance of the class with a publisher source.
///
public Pta(ITValuePublisher source, int longPeriod = 250, int shortPeriod = 40)
: this(longPeriod, shortPeriod)
{
source.Pub += (object? _, in TValueEventArgs args) => Update(args.Value, args.IsNew);
}
///
/// Computes 2-pole Butterworth highpass filter coefficients.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeHpCoefficients(int period, out double c1, out double c2, out double c3)
{
double f = 1.414 * Math.PI / period;
double a1 = Math.Exp(-f);
double b1 = 2.0 * a1 * Math.Cos(f);
c2 = b1;
c3 = -(a1 * a1);
c1 = (1.0 + c2 - c3) * 0.25;
}
public override bool IsHot => _state.Count >= 2;
/// Primes the indicator with historical data.
public override void Prime(ReadOnlySpan source, TimeSpan? step = null)
{
foreach (double v in source)
{
Update(new TValue(DateTime.MinValue, v), isNew: true);
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public override TValue Update(TValue input, bool isNew = true)
{
if (isNew)
{
_p_state = _state;
}
else
{
_state = _p_state;
}
double src = input.Value;
ref State s = ref _state;
double result;
if (s.Count < 2)
{
// Bootstrap: not enough bars for HP differentiation
if (s.Count == 0)
{
s.Src1 = src;
s.Src2 = src;
}
else
{
s.Src2 = s.Src1;
s.Src1 = src;
}
s.Count++;
result = 0.0;
}
else
{
// 2nd-order difference: src - 2*src1 + src2
double diff = src - 2.0 * s.Src1 + s.Src2;
// HP1 (long period): hp1 = c1L*diff + c2L*hp1 + c3L*hp1_1
double hp1 = Math.FusedMultiplyAdd(_c1L, diff,
Math.FusedMultiplyAdd(_c2L, s.Hp1, _c3L * s.Hp1_1));
// HP2 (short period): hp2 = c1S*diff + c2S*hp2 + c3S*hp2_1
double hp2 = Math.FusedMultiplyAdd(_c1S, diff,
Math.FusedMultiplyAdd(_c2S, s.Hp2, _c3S * s.Hp2_1));
// Trend = HP1 - HP2 (bandpass between shortPeriod and longPeriod)
result = hp1 - hp2;
// Update HP state
s.Hp1_1 = s.Hp1;
s.Hp1 = hp1;
s.Hp2_1 = s.Hp2;
s.Hp2 = hp2;
// Update source history
s.Src2 = s.Src1;
s.Src1 = src;
s.Count++;
}
Last = new TValue(input.Time, result);
PubEvent(Last, isNew);
return Last;
}
/// Updates with a full TSeries and returns results.
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
public override TSeries Update(TSeries source)
{
if (source.Count == 0)
{
return [];
}
var resultValues = new double[source.Count];
Batch(source.Values, resultValues, LongPeriod, ShortPeriod);
var result = new TSeries();
var times = source.Times;
for (int i = 0; i < source.Count; i++)
{
result.Add(new TValue(times[i], resultValues[i]));
}
// Sync internal state
int len = source.Count;
if (len >= 2)
{
var replay = new Pta(LongPeriod, ShortPeriod);
for (int i = 0; i < len; i++)
{
replay.Update(new TValue(times[i], source.Values[i]));
}
_state = replay._state;
}
_p_state = _state;
return result;
}
/// Static batch on TSeries.
public static TSeries Batch(TSeries source, int longPeriod = 250, int shortPeriod = 40)
{
var indicator = new Pta(longPeriod, shortPeriod);
return indicator.Update(source);
}
///
/// Static batch calculation on spans. Zero allocation on the hot path.
///
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Batch(ReadOnlySpan source, Span output,
int longPeriod = 250, int shortPeriod = 40)
{
if (source.Length != output.Length)
{
throw new ArgumentException("Source and output spans must be of equal length.", nameof(output));
}
if (source.Length == 0)
{
return;
}
ArgumentOutOfRangeException.ThrowIfLessThan(longPeriod, 3, nameof(longPeriod));
ArgumentOutOfRangeException.ThrowIfLessThan(shortPeriod, 2, nameof(shortPeriod));
if (longPeriod <= shortPeriod)
{
throw new ArgumentOutOfRangeException(nameof(longPeriod),
$"longPeriod ({longPeriod}) must be greater than shortPeriod ({shortPeriod}).");
}
// Precompute coefficients
ComputeHpCoefficients(longPeriod, out double c1L, out double c2L, out double c3L);
ComputeHpCoefficients(shortPeriod, out double c1S, out double c2S, out double c3S);
// Bar 0 and 1: output = 0 (not enough history for 2nd-order diff)
output[0] = 0.0;
if (source.Length < 2)
{
return;
}
output[1] = 0.0;
double hp1 = 0, hp1_1 = 0;
double hp2 = 0, hp2_1 = 0;
for (int i = 2; i < source.Length; i++)
{
double diff = source[i] - 2.0 * source[i - 1] + source[i - 2];
double newHp1 = Math.FusedMultiplyAdd(c1L, diff,
Math.FusedMultiplyAdd(c2L, hp1, c3L * hp1_1));
double newHp2 = Math.FusedMultiplyAdd(c1S, diff,
Math.FusedMultiplyAdd(c2S, hp2, c3S * hp2_1));
output[i] = newHp1 - newHp2;
hp1_1 = hp1; hp1 = newHp1;
hp2_1 = hp2; hp2 = newHp2;
}
}
/// Calculate factory returning results and indicator.
public static (TSeries Results, Pta Indicator) Calculate(TSeries source,
int longPeriod = 250, int shortPeriod = 40)
{
var indicator = new Pta(longPeriod, shortPeriod);
TSeries results = indicator.Update(source);
return (results, indicator);
}
public override void Reset()
{
_state = State.New();
_p_state = _state;
}
}