From a1492e23a83ae512a8ee29f1989db0947566c2ce Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 14 Oct 2018 16:17:09 +0200 Subject: [PATCH] Correlation: switch auto-correlation back to accept arrays to avoid repeated enumerable iteration --- src/Numerics/Fit.cs | 1 - src/Numerics/Statistics/Correlation.cs | 86 +++++++++++--------------- 2 files changed, 37 insertions(+), 50 deletions(-) diff --git a/src/Numerics/Fit.cs b/src/Numerics/Fit.cs index 8116a5f5..8d65c986 100644 --- a/src/Numerics/Fit.cs +++ b/src/Numerics/Fit.cs @@ -32,7 +32,6 @@ using System.Linq; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearRegression; using MathNet.Numerics.Providers.LinearAlgebra; -using MathNet.Numerics.Statistics; namespace MathNet.Numerics { diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index 2de8517e..de798291 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -43,23 +43,23 @@ namespace MathNet.Numerics.Statistics public static class Correlation { /// - /// Autocorrelation function (ACF) based on FFT for all possible lags k. + /// Auto-correlation function (ACF) based on FFT for all possible lags k. /// /// Data array to calculate auto correlation for. /// An array with the ACF as a function of the lags k. - public static double[] Auto(IEnumerable x) + public static double[] Auto(double[] x) { - return AutoCorrelationFft(x, 0, x.Count() - 1); + return AutoCorrelationFft(x, 0, x.Length - 1); } /// - /// Autocorrelation function (ACF) based on FFT for lags between kMin and kMax. + /// Auto-correlation function (ACF) based on FFT for lags between kMin and kMax. /// /// The data array to calculate auto correlation for. /// Max lag to calculate ACF for must be positive and smaller than x.Length. /// Min lag to calculate ACF for (0 = no shift with acf=1) must be zero or positive and smaller than x.Length. /// An array with the ACF as a function of the lags k. - public static double[] Auto(IEnumerable x, int kMax, int kMin = 0) + public static double[] Auto(double[] x, int kMax, int kMin = 0) { // assert max and min in proper order var kMax2 = Math.Max(kMax, kMin); @@ -69,12 +69,12 @@ namespace MathNet.Numerics.Statistics } /// - /// Autocorrelation function based on FFT for lags k. + /// Auto-correlation function based on FFT for lags k. /// /// The data array to calculate auto correlation for. /// Array with lags to calculate ACF for. /// An array with the ACF as a function of the lags k. - public static double[] Auto(IEnumerable x, int[] k) + public static double[] Auto(double[] x, int[] k) { if (k == null) { @@ -86,90 +86,78 @@ namespace MathNet.Numerics.Statistics throw new ArgumentException("k"); } - var k_min = k.Min(); - var k_max = k.Max(); + var kMin = k.Min(); + var kMax = k.Max(); + // get acf between full range - var acf = AutoCorrelationFft(x, k_min, k_max); + var acf = AutoCorrelationFft(x, kMin, kMax); // map output by indexing var result = new double[k.Length]; for (int i = 0; i < result.Length; i++) { - result[i] = acf[k[i] - k_min]; + result[i] = acf[k[i] - kMin]; } return result; } /// - /// The internal core method for calculating the autocorrelation. + /// The internal method for calculating the auto-correlation. /// - /// The data array to calculate auto correlation for - /// Min lag to calculate ACF for (0 = no shift with acf=1) must be zero or positive and smaller than x.Length - /// Max lag (EXCLUSIVE) to calculate ACF for must be positive and smaller than x.Length + /// The data array to calculate auto-correlation for + /// Min lag to calculate ACF for (0 = no shift with acf=1) must be zero or positive and smaller than x.Length + /// Max lag (EXCLUSIVE) to calculate ACF for must be positive and smaller than x.Length /// An array with the ACF as a function of the lags k. - private static double[] AutoCorrelationFft(IEnumerable x, int k_low, int k_high) + private static double[] AutoCorrelationFft(double[] x, int kLow, int kHigh) { if (x == null) throw new ArgumentNullException(nameof(x)); - int N = x.Count(); // Sample size - - if (k_low < 0 || k_low >= N) - throw new ArgumentOutOfRangeException(nameof(k_low), "kMin must be zero or positive and smaller than x.Length"); - if (k_high < 0 || k_high >= N) - throw new ArgumentOutOfRangeException(nameof(k_high), "kMax must be positive and smaller than x.Length"); + int N = x.Length; // Sample size + if (kLow < 0 || kLow >= N) + throw new ArgumentOutOfRangeException(nameof(kLow), "kMin must be zero or positive and smaller than x.Length"); + if (kHigh < 0 || kHigh >= N) + throw new ArgumentOutOfRangeException(nameof(kHigh), "kMax must be positive and smaller than x.Length"); if (N < 1) return new double[0]; int nFFT = Euclid.CeilingToPowerOfTwo(N) * 2; - Complex[] x_fft = new Complex[nFFT]; - Complex[] x_fft2 = new Complex[nFFT]; + Complex[] xFFT = new Complex[nFFT]; + Complex[] xFFT2 = new Complex[nFFT]; - double x_dash = Statistics.Mean(x); + double xDash = ArrayStatistics.Mean(x); double xArrNow = 0.0d; - using (IEnumerator iex = x.GetEnumerator()) + // copy values in range and substract mean - all the remaining parts are padded with zero. + for (int i = 0; i < x.Length; i++) { - for (int ii = 0; ii < nFFT; ii++) - { - - if (ii < N) - { - if (!iex.MoveNext()) - throw new ArgumentOutOfRangeException(nameof(x)); - xArrNow = iex.Current; - x_fft[ii] = new Complex(xArrNow - x_dash, 0.0); // copy values in range and substract mean - } - else - x_fft[ii] = new Complex(0.0, 0.0); // pad all remaining points - } - + xFFT[i] = new Complex(x[i] - xDash, 0.0); // copy values in range and substract mean } - Fourier.Forward(x_fft, FourierOptions.Matlab); + Fourier.Forward(xFFT, FourierOptions.Matlab); // maybe a Vector implementation here would be faster - for (int ii = 0; ii < x_fft.Length; ii++) + for (int i = 0; i < xFFT.Length; i++) { - x_fft2[ii] = Complex.Multiply(x_fft[ii], Complex.Conjugate(x_fft[ii])); + xFFT2[i] = Complex.Multiply(xFFT[i], Complex.Conjugate(xFFT[i])); } - Fourier.Inverse(x_fft2, FourierOptions.Matlab); + Fourier.Inverse(xFFT2, FourierOptions.Matlab); - double acf_Val1 = x_fft2[0].Real; + double dc = xFFT2[0].Real; - double[] acf_Vec = new double[k_high - k_low + 1]; + double[] result = new double[kHigh - kLow + 1]; // normalize such that acf[0] would be 1.0 - for (int ii = 0; ii < (k_high - k_low + 1); ii++) + for (int i = 0; i < (kHigh - kLow + 1); i++) { - acf_Vec[ii] = x_fft2[k_low + ii].Real / acf_Val1; + result[i] = xFFT2[kLow + i].Real / dc; } - return acf_Vec; + return result; } ///