From 9e70ce703dae44c91e7657350e2954f08401e707 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 10 Jul 2018 12:56:01 +0200 Subject: [PATCH] some adjustments, grooming and dependency fixes --- src/Numerics/Statistics/Correlation.cs | 53 +++++++++++++++----------- 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index a5389bde..dd183276 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -30,6 +30,8 @@ using System; using System.Collections.Generic; using System.Linq; +using System.Numerics; +using MathNet.Numerics.IntegralTransforms; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.Properties; @@ -41,7 +43,6 @@ namespace MathNet.Numerics.Statistics public static class Correlation { - /// /// autocorrelation function (ACF) based on fft (usually faster then direct brute force implementation) for all possible lags k /// First element is hidden since ACF(k = 0) = 1 @@ -49,7 +50,7 @@ namespace MathNet.Numerics.Statistics /// an array with the ACF as a function of the lags k public static double[] AutoCorrelation(IEnumerable x) { - return autoCorrFft(x, tmpk, 0, x.Count()-1); + return autoCorrelationFft(x, 0, x.Count() - 1); } /// @@ -58,14 +59,15 @@ namespace MathNet.Numerics.Statistics /// First element is hidden since ACF(k = 0) = 1 /// the data array to calculate auto correlation for /// max lag to calculate ACF for must be positive and smaller than x.Length-1 - /// min lag to calculate ACF for (0 = no shift with acf=1) must be zero or positive and smaller than x.Length-1 + /// min lag to calculate ACF for (0 = no shift with acf=1) must be zero or positive and smaller than x.Length-1 + /// an array with the ACF as a function of the lags k public static double[] AutoCorrelation(IEnumerable x, int kMax, int kMin = 0) { // assert max and min in proper order var kMax2 = Math.Max(kMax, kMin); var kMin2 = Math.Min(kMax, kMin); - return (CorrCov.autoCorrFft(x, kMin2, kMax2)); + return (autoCorrelationFft(x, kMin2, kMax2)); } @@ -74,6 +76,7 @@ namespace MathNet.Numerics.Statistics /// First element is skipped since ACF(k = 0) = 1 /// 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[] AutoCorrelation(IEnumerable x, int[] k) { if (k == null) @@ -83,7 +86,7 @@ namespace MathNet.Numerics.Statistics throw new ArgumentException("k"); // get acf between full range - var acf = autoCorrFft(x, k.Min(), k.Max()); + var acf = autoCorrelationFft(x, k.Min(), k.Max()); // map output by indexing var acfReturn = new double[k.Length]; @@ -93,9 +96,16 @@ namespace MathNet.Numerics.Statistics return acfReturn; } - private static double[] autoCorrFft(IEnumerable x, int k_low, int k_high) + /// + /// this is the internal core method for calculating the autocorrelation + /// + /// 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-1 + /// max lag to calculate ACF for must be positive and smaller than x.Length-1 + /// an array with the ACF as a function of the lags k + private static double[] autoCorrelationFft(IEnumerable x, int k_low, int k_high) { - if(x == null) + if (x == null) throw new ArgumentNullException("x"); if (k_low < 0 || k_low >= x.Count()) @@ -107,31 +117,34 @@ namespace MathNet.Numerics.Statistics return new double[0]; int N = x.Count(); // Sample size - int nFFT = (int)Math.Pow(2, Euclid.CeilingToPowerOfTwo(N)); + int nFFT = Euclid.CeilingToPowerOfTwo(N) * 2; Complex[] x_fft = new Complex[nFFT]; Complex[] x_fft2 = new Complex[nFFT]; double x_dash = Statistics.Mean(x); + double xArrNow = 0.0d; - using (IEnumerator ieX = x.GetEnumerator()) + using (IEnumerator iex = x.GetEnumerator()) { - for (int ii = 0; ii < x_fft.Length; ii++) + for (int ii = 0; ii < nFFT; ii++) { - if (!ieX.MoveNext()) - throw new ArgumentOutOfRangeException("x", Resources.ArgumentArraysSameLength); + if (!iex.MoveNext()) + { + throw new ArgumentOutOfRangeException("x"); + } + xArrNow = iex.Current; if (ii < N) - x_fft[ii] = new Complex(ieX.Current - x_dash, 0.0); // copy values in range and substract mean + 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 } - } + } Fourier.Forward(x_fft, FourierOptions.Matlab); - // maybe a Vector implementation here would be faster for (int ii = 0; ii < x_fft.Length; ii++) { @@ -145,15 +158,9 @@ namespace MathNet.Numerics.Statistics double[] acf_Val = new double[k_high + 1]; // normalize such that acf[0] would be 1.0 and drop the first element - for (int ii = 0; ii < k_high + 1; ii++) - { - acf_Val[ii] = x_fft2[ii + 1].Real / acf_Val1; - } - - // only return requested lags - for (int ii = 0; ii < (k_high - k_low + 1); ii++) + for (int ii = 0; ii < (k_high - k_low); ii++) { - acf_Vec[ii] = acf_Val[k_low + ii]; + acf_Vec[ii] = x_fft2[k_low + ii + 1].Real / acf_Val1; } return (acf_Vec);