using System;
using System.Collections.Generic; // Required for List and Queue
using System.Linq; // Required for .Any() and .ToArray()
using System.Diagnostics; // For Debug.WriteLine, if needed for warnings/errors

namespace MS
{
    // ICalculation interface defines a contract for updatable and resettable calculation objects.
    // This allows different types of indicators or filters to be used interchangeably
    // if they adhere to this structure.
    public interface ICalculation
    {
        // Update processes a new data point (inputValue).
        // inputValue: The new data point to incorporate into the calculation.
        // commitState: If true, the internal state of the calculator is permanently updated.
        //              If false, the calculation is performed as a "peek" or "what-if" scenario
        //              without altering the persistent state of the calculator.
        // Returns: The calculated value (e.g., median, MMI, filtered output) as a double.
        //          May return double.NaN if the calculation cannot be performed (e.g., insufficient data).
        double Update(double inputValue, bool commitState = true);

        // Reset re-initializes the calculator to its default state.
        // This typically involves clearing any stored historical data or state variables,
        // making the calculator ready for a new sequence of data.
        void Reset();
    }

    // Median class calculates a rolling median over a specified period.
    // The median is the middle value of a data set sorted in ascending order.
    // This implementation allows for committing state changes or calculating a "what-if" median.
    public class Median : ICalculation
    {
        private readonly int _period; // The number of data points in the rolling window.
        private readonly List<double> _windowValues; // Stores the sorted values within the current window. Used for efficient median calculation.
        private readonly Queue<double> _windowOrder;  // Stores the values in insertion order to manage the rolling window (FIFO).

        // Constructor for the Median calculator.
        // period: The number of data points to include in the median calculation window. Must be positive.
        public Median(int period)
        {
            if (period <= 0)
                throw new ArgumentOutOfRangeException(nameof(period), "Period must be positive.");

            _period = period;
            _windowValues = new List<double>(_period); // Initialize with capacity for performance.
            _windowOrder = new Queue<double>(_period);  // Initialize with capacity.
        }

        // Updates the median calculation with a new data point (inputValue).
        // If commitState is true, the inputValue is added to the window, and the oldest value is removed if the window is full.
        // The median is then calculated from this updated (persisted) window.
        // If commitState is false, the median is calculated based on a temporary window:
        //   - If the persisted window is full, it simulates replacing the oldest values with inputValue (current logic uses P-1 smallest).
        //   - If the persisted window is not full, it adds inputValue to the current persisted values.
        // This temporary calculation does not alter the actual persisted state of _windowValues or _windowOrder.
        public double Update(double inputValue, bool commitState = true)
        {
            List<double> tempWindowForCalc = new List<double>(); // Temporary list for performing the median calculation.

            if (commitState) // Persist the new value and update the window.
            {
                _windowOrder.Enqueue(inputValue); // Add to FIFO queue to track age.
                int insertIndex = _windowValues.BinarySearch(inputValue); // Find position in sorted list.
                if (insertIndex < 0) insertIndex = ~insertIndex; // If not found, get complement for insertion index.
                _windowValues.Insert(insertIndex, inputValue); // Insert into sorted list.

                if (_windowOrder.Count > _period) // If window size exceeds period.
                {
                    double oldestValue = _windowOrder.Dequeue(); // Remove oldest from FIFO queue.
                    _windowValues.Remove(oldestValue); // Remove oldest from sorted list (can be slow if list is large).
                }
                tempWindowForCalc.AddRange(_windowValues); // Calculation based on the actual, updated window.
            }
            else // Calculate "what-if" median without persisting state.
            {
                if (_windowValues.Any()) // If there are existing values.
                {
                    if (_windowValues.Count == _period) // If the persisted window is full.
                    {
                        // Current "what-if" logic: takes the _period-1 smallest values from the current window
                        // and adds the new inputValue. This might not always reflect replacing the "oldest" actual item.
                        for (int i = 0; i < _period - 1; ++i) tempWindowForCalc.Add(_windowValues[i]);
                    }
                    else // Persisted window is not yet full.
                    {
                        tempWindowForCalc.AddRange(_windowValues); // Use all current values.
                    }
                }
                tempWindowForCalc.Add(inputValue); // Add the new input value to the temporary set.
            }

            if (!tempWindowForCalc.Any()) // If no data in the temporary window.
            {
                return double.NaN; // Cannot calculate median.
            }

            tempWindowForCalc.Sort(); // Median calculation requires a sorted list.

            int n = tempWindowForCalc.Count;
            double calculatedMedian;

            if (n == 0) return double.NaN; // Should be caught by .Any() check above.

            if (n % 2 == 0) // Even number of elements.
            {
                if (n >= 2) { // Median is the average of the two middle elements.
                    calculatedMedian = (tempWindowForCalc[n / 2 - 1] + tempWindowForCalc[n / 2]) / 2.0;
                } else {
                    return double.NaN; // Not enough elements for a meaningful even median.
                }
            }
            else // Odd number of elements.
            {
                calculatedMedian = tempWindowForCalc[n / 2]; // Median is the middle element.
            }
            return calculatedMedian;
        }

        // Resets the Median calculator to its initial state.
        // Clears all stored values from the calculation window.
        public void Reset()
        {
            _windowValues.Clear(); // Empties the sorted list of values.
            _windowOrder.Clear();  // Empties the queue tracking insertion order.
        }
    }

    // MMI (Market Meanness Index) class calculates a rolling index indicating if a market is trending or mean-reverting.
    // It uses a Median calculator internally.
    // The MMI value is typically between 0 and 100.
    public class MMI : ICalculation
    {
        private readonly int _period; // The period for MMI and its internal Median calculation.
        private readonly Median _medianCalculator; // Internal Median calculator instance.
        private readonly Queue<double> _priceWindow;  // Stores raw price values in insertion order for MMI logic.

        // Constructor for the MMI calculator.
        // period: The number of data points for the MMI window. Must be at least 2.
        public MMI(int period)
        {
            if (period < 2) // MMI calculation requires at least 2 data points to compare.
                throw new ArgumentOutOfRangeException(nameof(period), "Period must be at least 2 for MMI calculation.");

            _period = period;
            _medianCalculator = new Median(period); // Initialize the Median helper.
            _priceWindow = new Queue<double>(period); // Initialize price queue.
        }

        // Updates the MMI calculation with a new data point (inputValue).
        // It first updates its internal Median calculator.
        // Then, it updates its own price window based on commitState.
        // Finally, it calculates MMI based on the median and the price window.
        // MMI is calculated as 100 * (nl + nh) / (period - 1), where:
        //   nl: count of prices > median and > previous price.
        //   nh: count of prices < median and < previous price.
        // Returns NaN if the median is NaN or if the price window for MMI isn't full (_period size).
        public double Update(double inputValue, bool commitState = true)
        {
            // Update the median, passing through the commitState.
            double calculatedMedian = _medianCalculator.Update(inputValue, commitState);

            double[] pricesForMmiLoop; // Array to hold prices for the MMI loop.

            if (commitState) // Persist inputValue in MMI's price window.
            {
                _priceWindow.Enqueue(inputValue);
                if (_priceWindow.Count > _period)
                {
                    _priceWindow.Dequeue(); // Maintain window size.
                }
                pricesForMmiLoop = _priceWindow.ToArray(); // Use the actual, updated price window.
            }
            else // "What-if" scenario for MMI price window.
            {
                List<double> tempList = new List<double>();
                if (_priceWindow.Any())
                {
                    if (_priceWindow.Count == _period) // If full, simulate replacing the oldest.
                    {
                        var tempQueue = new Queue<double>(_priceWindow.Skip(1)); // Create a new queue, skipping the oldest.
                        tempQueue.Enqueue(inputValue); // Add the new value.
                        pricesForMmiLoop = tempQueue.ToArray();
                    }
                    else // Not full, just add to current MMI prices.
                    {
                        tempList.AddRange(_priceWindow);
                        tempList.Add(inputValue);
                        pricesForMmiLoop = tempList.ToArray();
                    }
                }
                else // Price window is empty.
                {
                     tempList.Add(inputValue);
                     pricesForMmiLoop = tempList.ToArray();
                }
            }

            // MMI can only be calculated if median is valid and the price window for the loop is full.
            if (double.IsNaN(calculatedMedian) || pricesForMmiLoop.Length != _period)
            {
                return double.NaN;
            }

            int nl = 0; // Counter for "trending up" conditions.
            int nh = 0; // Counter for "trending down" conditions.
            // Loop starts from 1 because it compares price[i] with price[i-1].
            for (int i = 1; i < _period; i++)
            {
                double priceInWindow = pricesForMmiLoop[i];
                double prevPriceInWindow = pricesForMmiLoop[i - 1];
                // Condition for nl: price is above median and current price is greater than previous price.
                if (priceInWindow > calculatedMedian && priceInWindow > prevPriceInWindow) nl++;
                // Condition for nh: price is below median and current price is less than previous price.
                else if (priceInWindow < calculatedMedian && priceInWindow < prevPriceInWindow) nh++;
            }
            // Denominator (_period - 1) is safe because constructor enforces _period >= 2.
            return 100.0 * (nl + nh) / (_period - 1);
        }

        // Resets the MMI calculator to its initial state.
        // This involves resetting its internal Median calculator and clearing its own price window.
        public void Reset()
        {
            _medianCalculator.Reset(); // Reset the dependent Median calculator.
            _priceWindow.Clear();      // Clear the MMI's specific price data.
        }
    }

    // LowPassFilter class implements a 2nd order IIR (Infinite Impulse Response) digital filter.
    // It can be configured as various types of low-pass filters (Butterworth, Bessel, Chebyshev)
    // using static factory methods.
    // The filter equation is:
    // y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]
    public class LowPassFilter : ICalculation
    {
        // Filter coefficients for the numerator (feedforward)
        private readonly double _b0, _b1, _b2;
        // Filter coefficients for the denominator (feedback), signs are for the direct subtraction in formula.
        private readonly double _a1, _a2;

        // State variables storing previous input values.
        private double _x1; // x[n-1]
        private double _x2; // x[n-2]
        // State variables storing previous output values.
        private double _y1; // y[n-1]
        private double _y2; // y[n-2]

        // Constructor to initialize the filter with specific coefficients.
        // Used by the static factory methods or for custom filter designs.
        // b0_coeff, b1_coeff, b2_coeff: Numerator (feedforward) coefficients.
        // a1_coeff, a2_coeff: Denominator (feedback) coefficients for the form y[n] = ... - a1*y[n-1] - a2*y[n-2].
        public LowPassFilter(double b0_coeff, double b1_coeff, double b2_coeff, double a1_coeff, double a2_coeff)
        {
            _b0 = b0_coeff;
            _b1 = b1_coeff;
            _b2 = b2_coeff;
            _a1 = a1_coeff; // Corresponds to a1 in y[n] = ... - a1*y[n-1]
            _a2 = a2_coeff; // Corresponds to a2 in y[n] = ... - a2*y[n-2]
            Reset(); // Initialize state variables to zero.
        }

        // Factory method to create a 2nd order Butterworth low-pass filter.
        // Butterworth filters are known for their maximally flat passband response.
        // cutoffFrequency: The -3dB cutoff frequency of the filter.
        // samplingFrequency: The sampling frequency of the input signal.
        public static LowPassFilter CreateButterworthLowPass(double cutoffFrequency, double samplingFrequency)
        {
            if (cutoffFrequency <= 0)
                throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), "Cutoff frequency must be positive.");
            if (samplingFrequency <= 0)
                throw new ArgumentOutOfRangeException(nameof(samplingFrequency), "Sampling frequency must be positive.");
            // Cutoff frequency must be less than Nyquist frequency (samplingFrequency / 2.0).
            if (cutoffFrequency >= samplingFrequency / 2.0)
                throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), $"Cutoff frequency ({cutoffFrequency}) must be strictly less than half the sampling frequency ({samplingFrequency / 2.0}).");

            double b0_calc, b1_calc, b2_calc, a1_calc, a2_calc;
            // Prewarp the cutoff frequency for the bilinear transform.
            // C = tan(pi * fc / fs)
            double C = Math.Tan(Math.PI * cutoffFrequency / samplingFrequency);

            if (Math.Abs(C) < 1e-12) // Handle extremely low cutoff (approaching DC).
            {
                b0_calc = 0; b1_calc = 0; b2_calc = 0; // Effectively blocks AC signals.
                // Coefficients for y[n] = 2*y[n-1] - y[n-2] (integrator-like behavior for output).
                a1_calc = -2.0; // So that -a1_calc*y[n-1] becomes +2*y[n-1].
                a2_calc = 1.0;  // So that -a2_calc*y[n-2] becomes -1*y[n-2].
            }
            else
            {
                double sqrt2 = Math.Sqrt(2.0);
                // Denominator factor for Butterworth coefficients after bilinear transform.
                double D_butter = 1.0 / (1.0 + sqrt2 * C + C * C);

                // Numerator coefficients
                b0_calc = C * C * D_butter;
                b1_calc = 2.0 * b0_calc;
                b2_calc = b0_calc;

                // Denominator coefficients (for 1 + a1_std*z^-1 + a2_std*z^-2)
                // These are directly used as _a1, _a2 in the filter update equation.
                a1_calc = 2.0 * (C * C - 1.0) * D_butter;
                a2_calc = (1.0 - sqrt2 * C + C * C) * D_butter;
            }
            return new LowPassFilter(b0_calc, b1_calc, b2_calc, a1_calc, a2_calc);
        }

        // Factory method to create a 2nd order Bessel low-pass filter.
        // Bessel filters are known for their maximally flat group delay (linear phase response in passband).
        // cutoffFrequency: The -3dB cutoff frequency. Note: Bessel filter design often refers to other frequency points.
        //                  This implementation aims for -3dB at the specified cutoffFrequency.
        // samplingFrequency: The sampling frequency of the input signal.
        public static LowPassFilter CreateBesselLowPass(double cutoffFrequency, double samplingFrequency)
        {
            if (cutoffFrequency <= 0) throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), "Cutoff frequency must be positive.");
            if (samplingFrequency <= 0) throw new ArgumentOutOfRangeException(nameof(samplingFrequency), "Sampling frequency must be positive.");
            if (cutoffFrequency >= samplingFrequency / 2.0) throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), $"Cutoff frequency ({cutoffFrequency}) must be strictly less than half the sampling frequency ({samplingFrequency / 2.0}).");

            // C_bessel = tan(pi * fc / fs), prewarped frequency factor.
            double C_bessel = Math.Tan(Math.PI * cutoffFrequency / samplingFrequency);
            if (Math.Abs(C_bessel) < 1e-12) // Handle low cutoff.
            {
                return new LowPassFilter(0, 0, 0, -2.0, 1.0); // Similar to Butterworth low cutoff.
            }

            // Coefficients derived from analog Bessel prototype (s^2 + 3s + 3),
            // frequency scaled for -3dB at cutoffFrequency, then transformed using bilinear method.
            // omega_c_analog_warped = 2 * fs * tan(pi * fc / fs)
            double omega_c_analog_warped = 2.0 * samplingFrequency * C_bessel;
            // k_bessel_norm_3db is a scaling factor for the Bessel prototype's natural frequency
            // to achieve the -3dB point at the desired warped analog cutoff. Value is approx. 1.3617.
            double k_bessel_norm_3db = 1.3617070504; // Sqrt(3 * (Sqrt(5.0) - 1.0) / 2.0)
            double omega_0_scaled = omega_c_analog_warped / k_bessel_norm_3db; // Scaled analog natural frequency for prototype.

            // Analog prototype denominator: s^2 + A_s*s + B_s
            // For Bessel n=2, scaled: s^2 + (3*omega_0_scaled)*s + (3*omega_0_scaled^2)
            double as_coeff_analog = 3.0 * omega_0_scaled; // Coefficient for s-term.
            double bs_coeff_analog = 3.0 * omega_0_scaled * omega_0_scaled; // Constant term (and numerator for DC gain 1).

            // K_T_bilinear = 2 * fs, factor in bilinear transform s = K_T * (1-z^-1)/(1+z^-1) (if K_T is used for 2/T_sample)
            // Here, K_T_bilinear is 2*fs.
            double K_T_bilinear = 2.0 * samplingFrequency;

            // Denominator after bilinear transform of H(s) = bs_coeff_analog / (s^2 + as_coeff_analog*s + bs_coeff_analog)
            // gives A0_prime_denom = (K_T^2) + as_coeff_analog*K_T + bs_coeff_analog for normalization.
            double A0_prime_denom = K_T_bilinear * K_T_bilinear + as_coeff_analog * K_T_bilinear + bs_coeff_analog;
            if (Math.Abs(A0_prime_denom) < 1e-12) throw new InvalidOperationException("Denominator term A0_prime_denom is zero in Bessel filter calculation.");

            // Digital filter coefficients b_i = Bi_prime / A0_prime_denom, a_i = Ai_prime_std / A0_prime_denom
            double b0_calc = bs_coeff_analog / A0_prime_denom;
            double b1_calc = (2.0 * bs_coeff_analog) / A0_prime_denom; // From bilinear transform structure.
            double b2_calc = bs_coeff_analog / A0_prime_denom;

            // Denominator coefficients for 1 + a1_std*z^-1 + a2_std*z^-2 form.
            double a1_calc = (2.0 * bs_coeff_analog - 2.0 * K_T_bilinear * K_T_bilinear) / A0_prime_denom;
            double a2_calc = (K_T_bilinear * K_T_bilinear - as_coeff_analog * K_T_bilinear + bs_coeff_analog) / A0_prime_denom;

            return new LowPassFilter(b0_calc, b1_calc, b2_calc, a1_calc, a2_calc);
        }

        // Factory method to create a 2nd order Chebyshev Type I low-pass filter.
        // Chebyshev Type I filters have a steeper roll-off than Butterworth, but exhibit ripple in the passband.
        // cutoffFrequency: The frequency at which the passband ripple ends (e.g., gain drops by rippleDb).
        // samplingFrequency: The sampling frequency of the input signal.
        // rippleDb: The passband ripple in decibels (e.g., 0.5 dB, 1 dB). Must be positive.
        public static LowPassFilter CreateChebyshevLowPass(double cutoffFrequency, double samplingFrequency, double rippleDb)
        {
            if (cutoffFrequency <= 0) throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), "Cutoff frequency must be positive.");
            if (samplingFrequency <= 0) throw new ArgumentOutOfRangeException(nameof(samplingFrequency), "Sampling frequency must be positive.");
            if (cutoffFrequency >= samplingFrequency / 2.0) throw new ArgumentOutOfRangeException(nameof(cutoffFrequency), $"Cutoff frequency ({cutoffFrequency}) must be strictly less than half the sampling frequency ({samplingFrequency / 2.0}).");
            if (rippleDb <= 0) throw new ArgumentOutOfRangeException(nameof(rippleDb), "Passband ripple must be positive.");

            // C_cheby = tan(pi * fc / fs), prewarped frequency factor.
            double C_cheby = Math.Tan(Math.PI * cutoffFrequency / samplingFrequency);
            if (Math.Abs(C_cheby) < 1e-12) // Handle low cutoff.
            {
                 return new LowPassFilter(0,0,0, -2.0, 1.0); // Similar to Butterworth low cutoff.
            }

            // Epsilon is derived from the passband ripple.
            double epsilon = Math.Sqrt(Math.Pow(10, rippleDb / 10.0) - 1.0);
            int N_order = 2; // Filter order is fixed at 2.

            // Parameters for locating analog prototype poles on an ellipse.
            // alpha_val determines the eccentricity of the ellipse.
            double alpha_val = Math.Asinh(1.0 / epsilon) / N_order;
            // For N=2, poles are at s_k = -sigma_k +/- j*omega_k.
            // These are scaled by C_cheby (warped cutoff frequency) for the specific filter.
            double sin_pi_over_2N = Math.Sin(Math.PI / (2.0 * N_order)); // sin(pi/4) for N=2
            double cos_pi_over_2N = Math.Cos(Math.PI / (2.0 * N_order)); // cos(pi/4) for N=2

            double sigma_k_scaled = Math.Sinh(alpha_val) * sin_pi_over_2N * C_cheby;
            double omega_k_scaled = Math.Cosh(alpha_val) * cos_pi_over_2N * C_cheby;

            // Analog prototype denominator: s^2 + A_s*s + B_s
            // Derived from poles: (s + sigma_k_scaled)^2 + omega_k_scaled^2
            double as_coeff_analog = 2.0 * sigma_k_scaled; // s-term coefficient.
            double bs_coeff_analog = sigma_k_scaled * sigma_k_scaled + omega_k_scaled * omega_k_scaled; // Constant term.

            // Analog numerator: For Chebyshev Type I, N=2 (even order), DC gain H(0) = 1.
            // So, numerator constant of analog prototype is bs_coeff_analog.
            double num_s_coeff = bs_coeff_analog;

            double K_T_bilinear = 2.0 * samplingFrequency; // Factor for bilinear transform.

            // Denominator for normalization after bilinear transform.
            double A0_prime_denom = K_T_bilinear * K_T_bilinear + as_coeff_analog * K_T_bilinear + bs_coeff_analog;
            if (Math.Abs(A0_prime_denom) < 1e-12) throw new InvalidOperationException("Denominator term A0_prime_denom is zero in Chebyshev filter calculation.");

            // Digital filter coefficients.
            double b0_calc = num_s_coeff / A0_prime_denom;
            double b1_calc = (2.0 * num_s_coeff) / A0_prime_denom;
            double b2_calc = num_s_coeff / A0_prime_denom;

            double a1_calc = (2.0 * bs_coeff_analog - 2.0 * K_T_bilinear * K_T_bilinear) / A0_prime_denom;
            double a2_calc = (K_T_bilinear * K_T_bilinear - as_coeff_analog * K_T_bilinear + bs_coeff_analog) / A0_prime_denom;

            return new LowPassFilter(b0_calc, b1_calc, b2_calc, a1_calc, a2_calc);
        }

        // Updates the filter with a new input sample (inputValue).
        // Applies the difference equation:
        // y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]
        // If commitState is true, the internal state variables (previous inputs and outputs) are updated.
        // If false, the output is calculated using the current state, but the state is not advanced.
        public double Update(double inputValue, bool commitState = true)
        {
            // Calculate the current output y[n] based on current input x[n] and previous states.
            double output = _b0 * inputValue + _b1 * _x1 + _b2 * _x2 - _a1 * _y1 - _a2 * _y2;

            if (commitState) // If true, update the filter's internal state.
            {
                // Shift previous inputs: x[n-2] = x[n-1], x[n-1] = x[n] (current inputValue).
                _x2 = _x1;
                _x1 = inputValue;
                // Shift previous outputs: y[n-2] = y[n-1], y[n-1] = y[n] (current calculated output).
                _y2 = _y1;
                _y1 = output;
            }
            return output; // Return the filtered output.
        }

        // Resets the filter's internal state variables to zero.
        // This clears any memory of previous input or output values, useful for starting a new filtering sequence.
        public void Reset()
        {
            _x1 = 0; _x2 = 0; // Reset previous input states.
            _y1 = 0; _y2 = 0; // Reset previous output states.
        }
    }
}