using System.Buffers;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
namespace QuanTAlib;
///
/// TukeyBiweight: Tukey's Biweight (Bisquare) Loss
///
///
/// Tukey's Biweight is a robust loss function that completely rejects outliers
/// beyond a threshold c. Unlike Huber loss which downweights outliers, Tukey's
/// biweight assigns zero weight to extreme outliers, making it highly resistant
/// to contaminated data.
///
/// Formula:
/// ρ(x) = (c²/6) * (1 - (1 - (x/c)²)³) for |x| ≤ c
/// ρ(x) = c²/6 for |x| > c
///
/// Key properties:
/// - Completely rejects outliers beyond threshold c
/// - Redescending: influence function goes to zero for large errors
/// - Common c values: 4.685 (95% efficiency), 6.0 (more permissive)
/// - More robust than Huber for heavily contaminated data
/// - Smooth and differentiable everywhere
///
[SkipLocalsInit]
public sealed class TukeyBiweight : BiInputIndicatorBase
{
private readonly double _cSquaredOver6;
private const double DefaultC = 4.685; // 95% efficiency for normal distribution
public TukeyBiweight(int period, double c = DefaultC)
: base(period, $"TukeyBiweight({period},{c:F3})")
{
if (c <= 0)
{
throw new ArgumentException("Threshold c must be positive", nameof(c));
}
C = c;
_cSquaredOver6 = (c * c) / 6.0;
}
public double C { get; }
///
/// Computes Tukey's biweight loss for the error between actual and predicted values.
///
[MethodImpl(MethodImplOptions.AggressiveInlining)]
protected override double ComputeError(double actual, double predicted)
{
double error = actual - predicted;
double absError = Math.Abs(error);
if (absError > C)
{
return _cSquaredOver6;
}
double ratio = error / C;
double ratioSq = ratio * ratio;
double oneMinusRatioSq = 1.0 - ratioSq;
double cubed = oneMinusRatioSq * oneMinusRatioSq * oneMinusRatioSq;
return _cSquaredOver6 * (1.0 - cubed);
}
public static TSeries Batch(TSeries actual, TSeries predicted, int period, double c = DefaultC)
{
if (actual.Count != predicted.Count)
{
throw new ArgumentException("Actual and predicted series must have the same length", nameof(predicted));
}
int len = actual.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(actual.Values, predicted.Values, vSpan, period, c);
actual.Times.CopyTo(tSpan);
return new TSeries(t, v);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Batch(ReadOnlySpan actual, ReadOnlySpan predicted, Span output, int period, double c = DefaultC)
{
if (actual.Length != predicted.Length || actual.Length != output.Length)
{
throw new ArgumentException("All spans must have the same length", nameof(output));
}
if (period <= 0)
{
throw new ArgumentException("Period must be greater than 0", nameof(period));
}
if (c <= 0)
{
throw new ArgumentException("Threshold c must be positive", nameof(c));
}
int len = actual.Length;
if (len == 0)
{
return;
}
// Rent buffer for intermediate Tukey biweight errors
double[] rented = ArrayPool.Shared.Rent(len);
try
{
Span errors = rented.AsSpan(0, len);
// Step 1: Compute Tukey biweight errors using ErrorHelpers
ErrorHelpers.ComputeTukeyBiweightErrors(actual, predicted, errors, c);
// Step 2: Apply rolling mean
ErrorHelpers.ApplyRollingMean(errors, output, period);
}
finally
{
ArrayPool.Shared.Return(rented, clearArray: false);
}
}
public static (TSeries Results, TukeyBiweight Indicator) Calculate(TSeries actual, TSeries predicted, int period, double c = DefaultC)
{
var indicator = new TukeyBiweight(period, c);
TSeries results = Batch(actual, predicted, period, c);
return (results, indicator);
}
}