From 2a8a209202d43d9e926668f45f77080e33cb76d1 Mon Sep 17 00:00:00 2001 From: jmoo2880 Date: Sat, 25 Jul 2026 13:06:30 +1000 Subject: [PATCH 1/4] fix memory leak for early exit --- src/C/CO_AutoCorr.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/C/CO_AutoCorr.c b/src/C/CO_AutoCorr.c index dd475fb..5b20186 100644 --- a/src/C/CO_AutoCorr.c +++ b/src/C/CO_AutoCorr.c @@ -291,6 +291,7 @@ double CO_Embed2_Dist_tau_d_expfit_meandiff(const double y[], const int size) int nBins = num_bins_auto(d, size-tau-1); if (nBins == 0){ + free(d); return 0; } int * histCounts = malloc(nBins * sizeof(double)); From a7382f4f23e46393a189b42a229337e4e11c7f84 Mon Sep 17 00:00:00 2001 From: jmoo2880 Date: Sat, 25 Jul 2026 13:07:15 +1000 Subject: [PATCH 2/4] guard against size <= train_length in FC_LocalSimple --- src/C/FC_LocalSimple.c | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/C/FC_LocalSimple.c b/src/C/FC_LocalSimple.c index d88ead8..111db09 100644 --- a/src/C/FC_LocalSimple.c +++ b/src/C/FC_LocalSimple.c @@ -31,6 +31,11 @@ double FC_LocalSimple_mean_tauresrat(const double y[], const int size, const int } } + if(size <= train_length) + { + return NAN; + } + double * res = malloc((size - train_length) * sizeof *res); for (int i = 0; i < size - train_length; i++) @@ -66,6 +71,11 @@ double FC_LocalSimple_mean_stderr(const double y[], const int size, const int tr } } + if(size <= train_length) + { + return NAN; + } + double * res = malloc((size - train_length) * sizeof *res); for (int i = 0; i < size - train_length; i++) @@ -99,6 +109,11 @@ double FC_LocalSimple_mean1_tauresrat(const double y[], const int size){ double FC_LocalSimple_mean_taures(const double y[], const int size, const int train_length) { + if(size <= train_length) + { + return NAN; + } + double * res = malloc((size - train_length) * sizeof *res); // first z-score @@ -130,6 +145,11 @@ double FC_LocalSimple_lfit_taures(const double y[], const int size) // set tau from first AC zero crossing int train_length = co_firstzero(y, size, size); + if(size <= train_length) + { + return NAN; + } + double * xReg = malloc(train_length * sizeof * xReg); // double * yReg = malloc(train_length * sizeof * yReg); for(int i = 1; i < train_length+1; i++) From 42c2479d6991eb3cf5fe4049798caf86a2d7be80 Mon Sep 17 00:00:00 2001 From: jmoo2880 Date: Sat, 25 Jul 2026 13:09:41 +1000 Subject: [PATCH 3/4] use logical operators instead of bitwise in conditionals --- src/C/IN_AutoMutualInfoStats.c | 2 +- src/C/PD_PeriodicityWang.c | 4 ++-- src/C/SB_BinaryStats.c | 4 ++-- src/C/SP_Summaries.c | 4 ++-- src/C/splinefit.c | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/C/IN_AutoMutualInfoStats.c b/src/C/IN_AutoMutualInfoStats.c index e98de9e..7b09477 100644 --- a/src/C/IN_AutoMutualInfoStats.c +++ b/src/C/IN_AutoMutualInfoStats.c @@ -41,7 +41,7 @@ double IN_AutoMutualInfoStats_40_gaussian_fmmi(const double y[], const int size) // find first minimum of automutual information double fmmi = tau; for(int i = 1; i < tau-1; i++){ - if(ami[i] < ami[i-1] & ami[i] < ami[i+1]){ + if(ami[i] < ami[i-1] && ami[i] < ami[i+1]){ fmmi = i; // printf("found minimum at %i\n", i); break; diff --git a/src/C/PD_PeriodicityWang.c b/src/C/PD_PeriodicityWang.c index 75bfd9d..6ebf586 100644 --- a/src/C/PD_PeriodicityWang.c +++ b/src/C/PD_PeriodicityWang.c @@ -63,13 +63,13 @@ int PD_PeriodicityWang_th0_01(const double * y, const int size){ slopeIn = acf[i] - acf[i-1]; slopeOut = acf[i+1] - acf[i]; - if(slopeIn < 0 & slopeOut > 0) + if(slopeIn < 0 && slopeOut > 0) { // printf("trough at %i\n", i); troughs[nTroughs] = i; nTroughs += 1; } - else if(slopeIn > 0 & slopeOut < 0) + else if(slopeIn > 0 && slopeOut < 0) { // printf("peak at %i\n", i); peaks[nPeaks] = i; diff --git a/src/C/SB_BinaryStats.c b/src/C/SB_BinaryStats.c index 71fa3d8..58e2115 100644 --- a/src/C/SB_BinaryStats.c +++ b/src/C/SB_BinaryStats.c @@ -37,7 +37,7 @@ double SB_BinaryStats_diff_longstretch0(const double y[], const int size){ int maxstretch0 = 0; int last1 = 0; for(int i = 0; i < size-1; i++){ - if(yBin[i] == 1 | i == size-2){ + if(yBin[i] == 1 || i == size-2){ double stretch0 = i - last1; if(stretch0 > maxstretch0){ maxstretch0 = stretch0; @@ -75,7 +75,7 @@ double SB_BinaryStats_mean_longstretch1(const double y[], const int size){ int maxstretch1 = 0; int last1 = 0; for(int i = 0; i < size-1; i++){ - if(yBin[i] == 0 | i == size-2){ + if(yBin[i] == 0 || i == size-2){ double stretch1 = i - last1; if(stretch1 > maxstretch1){ maxstretch1 = stretch1; diff --git a/src/C/SP_Summaries.c b/src/C/SP_Summaries.c index 797d07b..b511df4 100644 --- a/src/C/SP_Summaries.c +++ b/src/C/SP_Summaries.c @@ -83,7 +83,7 @@ int welch(const double y[], const int size, const int NFFT, const double Fs, con *Pxx = malloc(Nout * sizeof(double)); for(int i = 0; i < Nout; i++){ (*Pxx)[i] = P[i]/KMU*dt; - if(i>0 & i < Nout-1){ + if(i>0 && i < Nout-1){ (*Pxx)[i] *= 2; } } @@ -148,7 +148,7 @@ double SP_Summaries_welch_rect(const double y[], const int size, const char what w[i] = 2*PI*f[i]; Sw[i] = S[i]/(2*PI); //printf("w[%i]=%1.3f, Sw[%i]=%1.3f\n", i, w[i], i, Sw[i]); - if(isinf(Sw[i]) | isinf(-Sw[i])){ + if(isinf(Sw[i]) || isinf(-Sw[i])){ return 0; } } diff --git a/src/C/splinefit.c b/src/C/splinefit.c index 9eb1712..802298c 100644 --- a/src/C/splinefit.c +++ b/src/C/splinefit.c @@ -603,7 +603,7 @@ int splinefit(const double *y, const int size, double *yOut) int breakInd = 1; for(int i = 0; i < size; i++){ - if(i >= breaks[breakInd] & breakInd= breaks[breakInd] && breakInd Date: Tue, 28 Jul 2026 05:39:48 +1000 Subject: [PATCH 4/4] replace vendored C code with git submodule drawing from catch22 --- .github/workflows/run_unit_tests.yaml | 2 + .gitmodules | 3 + MANIFEST.in | 1 + README.md | 8 + setup.py | 11 +- src/C/CO_AutoCorr.c | 552 ------------------ src/C/CO_AutoCorr.h | 35 -- src/C/DN_HistogramMode_10.c | 113 ---- src/C/DN_HistogramMode_10.h | 9 - src/C/DN_HistogramMode_5.c | 123 ---- src/C/DN_HistogramMode_5.h | 9 - src/C/DN_Mean.c | 11 - src/C/DN_Mean.h | 12 - src/C/DN_OutlierInclude.c | 213 ------- src/C/DN_OutlierInclude.h | 14 - src/C/DN_Spread_Std.c | 13 - src/C/DN_Spread_Std.h | 12 - src/C/FC_LocalSimple.c | 183 ------ src/C/FC_LocalSimple.h | 16 - src/C/IN_AutoMutualInfoStats.c | 54 -- src/C/IN_AutoMutualInfoStats.h | 16 - src/C/MD_hrv.c | 40 -- src/C/MD_hrv.h | 16 - src/C/PD_PeriodicityWang.c | 137 ----- src/C/PD_PeriodicityWang.h | 16 - src/C/SB_BinaryStats.c | 91 --- src/C/SB_BinaryStats.h | 17 - src/C/SB_CoarseGrain.c | 40 -- src/C/SB_CoarseGrain.h | 12 - src/C/SB_MotifThree.c | 375 ------------ src/C/SB_MotifThree.h | 12 - src/C/SB_TransitionMatrix.c | 163 ------ src/C/SB_TransitionMatrix.h | 15 - src/C/SC_FluctAnal.c | 350 ----------- src/C/SC_FluctAnal.h | 11 - src/C/SP_Summaries.c | 213 ------- src/C/SP_Summaries.h | 17 - src/C/butterworth.c | 298 ---------- src/C/butterworth.h | 15 - src/C/fft.c | 61 -- src/C/fft.h | 26 - src/C/helper_functions.c | 179 ------ src/C/helper_functions.h | 34 -- src/C/histcounts.c | 200 ------- src/C/histcounts.h | 22 - src/C/main.c | 420 -------------- src/C/main.h | 16 - src/C/runAllTS.sh | 71 --- src/C/splinefit.c | 801 -------------------------- src/C/splinefit.h | 16 - src/C/stats.c | 270 --------- src/C/stats.h | 28 - src/catch22 | 1 + src/{C => wrapper}/catch22_wrap.c | 0 54 files changed, 22 insertions(+), 5371 deletions(-) create mode 100644 .gitmodules delete mode 100644 src/C/CO_AutoCorr.c delete mode 100644 src/C/CO_AutoCorr.h delete mode 100644 src/C/DN_HistogramMode_10.c delete mode 100644 src/C/DN_HistogramMode_10.h delete mode 100644 src/C/DN_HistogramMode_5.c delete mode 100644 src/C/DN_HistogramMode_5.h delete mode 100644 src/C/DN_Mean.c delete mode 100644 src/C/DN_Mean.h delete mode 100644 src/C/DN_OutlierInclude.c delete mode 100644 src/C/DN_OutlierInclude.h delete mode 100644 src/C/DN_Spread_Std.c delete mode 100644 src/C/DN_Spread_Std.h delete mode 100644 src/C/FC_LocalSimple.c delete mode 100644 src/C/FC_LocalSimple.h delete mode 100644 src/C/IN_AutoMutualInfoStats.c delete mode 100644 src/C/IN_AutoMutualInfoStats.h delete mode 100644 src/C/MD_hrv.c delete mode 100644 src/C/MD_hrv.h delete mode 100644 src/C/PD_PeriodicityWang.c delete mode 100644 src/C/PD_PeriodicityWang.h delete mode 100644 src/C/SB_BinaryStats.c delete mode 100644 src/C/SB_BinaryStats.h delete mode 100644 src/C/SB_CoarseGrain.c delete mode 100644 src/C/SB_CoarseGrain.h delete mode 100644 src/C/SB_MotifThree.c delete mode 100644 src/C/SB_MotifThree.h delete mode 100644 src/C/SB_TransitionMatrix.c delete mode 100644 src/C/SB_TransitionMatrix.h delete mode 100644 src/C/SC_FluctAnal.c delete mode 100644 src/C/SC_FluctAnal.h delete mode 100644 src/C/SP_Summaries.c delete mode 100644 src/C/SP_Summaries.h delete mode 100644 src/C/butterworth.c delete mode 100644 src/C/butterworth.h delete mode 100644 src/C/fft.c delete mode 100644 src/C/fft.h delete mode 100644 src/C/helper_functions.c delete mode 100644 src/C/helper_functions.h delete mode 100644 src/C/histcounts.c delete mode 100644 src/C/histcounts.h delete mode 100644 src/C/main.c delete mode 100644 src/C/main.h delete mode 100755 src/C/runAllTS.sh delete mode 100644 src/C/splinefit.c delete mode 100644 src/C/splinefit.h delete mode 100644 src/C/stats.c delete mode 100644 src/C/stats.h create mode 160000 src/catch22 rename src/{C => wrapper}/catch22_wrap.c (100%) diff --git a/.github/workflows/run_unit_tests.yaml b/.github/workflows/run_unit_tests.yaml index 637ab41..ec05d46 100644 --- a/.github/workflows/run_unit_tests.yaml +++ b/.github/workflows/run_unit_tests.yaml @@ -12,6 +12,8 @@ jobs: python-version: ["3.8", "3.9", "3.10", "3.11"] steps: - uses: actions/checkout@v3 + with: + submodules: recursive - name: Setup python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..afd7b50 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "src/catch22"] + path = src/catch22 + url = https://github.com/DynamicsAndNeuralSystems/catch22.git diff --git a/MANIFEST.in b/MANIFEST.in index 11a5eb0..af6f3ed 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1 +1,2 @@ graft src +exclude src/catch22/.git diff --git a/README.md b/README.md index 077131d..85f5ca7 100644 --- a/README.md +++ b/README.md @@ -97,6 +97,14 @@ Thanks to [@jmoo2880](https://github.com/jmoo2880) for putting together a [demon ### Manual install +The upstream _catch22_ C sources are vendored as a git submodule, so clone with submodules before building from source: + +``` +git clone --recurse-submodules https://github.com/DynamicsAndNeuralSystems/pycatch22.git +``` + +(If you already cloned without `--recurse-submodules`, run `git submodule update --init` from the repository root.) + If you find issues with the `pip` install, you can also install using `setuptools`: ``` diff --git a/setup.py b/setup.py index fc2a169..d4bb667 100644 --- a/setup.py +++ b/setup.py @@ -2,11 +2,14 @@ import sysconfig import os -sourceDir = os.path.join("src", "C"); +# Upstream catch22 C sources live in a git submodule (src/catch22); the +# pycatch22-specific CPython wrapper lives alongside it in src/wrapper. +upstreamDir = os.path.join("src", "catch22", "C") +wrapDir = os.path.join("src", "wrapper") -sourceFileList = [os.path.join(sourceDir, file) for file in os.listdir(sourceDir) if file.endswith( +sourceFileList = [os.path.join(upstreamDir, file) for file in os.listdir(upstreamDir) if file.endswith( ".c") and not 'main' in file] - # and not (file == "sampen.c" or file == "run_features.c")] +sourceFileList.append(os.path.join(wrapDir, "catch22_wrap.c")) cflags = sysconfig.get_config_var('CFLAGS') if cflags is not None: @@ -19,7 +22,7 @@ # The c++ extension module: extension_mod = Extension(name = "catch22_C", sources = sourceFileList, - include_dirs = [sourceDir], + include_dirs = [upstreamDir], extra_compile_args = extra_compile_args) # Header files are here setup( diff --git a/src/C/CO_AutoCorr.c b/src/C/CO_AutoCorr.c deleted file mode 100644 index 5b20186..0000000 --- a/src/C/CO_AutoCorr.c +++ /dev/null @@ -1,552 +0,0 @@ -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -#include -#include -#include -#include - -#include "stats.h" -#include "fft.h" -#include "histcounts.h" - -#include "helper_functions.h" - -#ifndef CMPLX -#define CMPLX(x, y) ((cplx)((double)(x) + _Imaginary_I * (double)(y))) -#endif -#define pow2(x) (1 << x) - -int nextpow2(int n) -{ - n--; - n |= n >> 1; - n |= n >> 2; - n |= n >> 4; - n |= n >> 8; - n |= n >> 16; - n++; - return n; -} - -/* -static void apply_conj(cplx a[], int size, int normalize) -{ - switch(normalize) { - case(1): - for (int i = 0; i < size; i++) { - a[i] = conj(a[i]) / size; - } - break; - default: - for (int i = 0; i < size; i++) { - a[i] = conj(a[i]); - } - break; - } -} - */ - -void dot_multiply(cplx a[], cplx b[], int size) -{ - for (int i = 0; i < size; i++) { - a[i] = _Cmulcc(a[i], conj(b[i])); - } -} - -double * CO_AutoCorr(const double y[], const int size, const int tau[], const int tau_size) -{ - double m, nFFT; - m = mean(y, size); - nFFT = nextpow2(size) << 1; - - cplx * F = malloc(nFFT * sizeof *F); - cplx * tw = malloc(nFFT * sizeof *tw); - for (int i = 0; i < size; i++) { - - #if defined(__GNUC__) || defined(__GNUG__) - F[i] = CMPLX(y[i] - m, 0.0); - #elif defined(_MSC_VER) - cplx tmp = { y[i] - m, 0.0 }; - F[i] = tmp; - #endif - - } - for (int i = size; i < nFFT; i++) { - #if defined(__GNUC__) || defined(__GNUG__) - F[i] = CMPLX(0.0, 0.0); - #elif defined(_MSC_VER) - cplx tmp = { 0.0, 0.0 }; - F[i] = tmp; // CMPLX(0.0, 0.0); - #endif - - } - // size = nFFT; - - twiddles(tw, nFFT); - fft(F, nFFT, tw); - dot_multiply(F, F, nFFT); - fft(F, nFFT, tw); - cplx divisor = F[0]; - for (int i = 0; i < nFFT; i++) { - //F[i] = F[i] / divisor; - F[i] = _Cdivcc(F[i], divisor); - } - - double * out = malloc(tau_size * sizeof(out)); - for (int i = 0; i < tau_size; i++) { - out[i] = creal(F[tau[i]]); - } - free(F); - free(tw); - return out; -} - -double * co_autocorrs(const double y[], const int size) -{ - double m, nFFT; - m = mean(y, size); - nFFT = nextpow2(size) << 1; - - cplx * F = malloc(nFFT * 2 * sizeof *F); - cplx * tw = malloc(nFFT * 2 * sizeof *tw); - for (int i = 0; i < size; i++) { - - #if defined(__GNUC__) || defined(__GNUG__) - F[i] = CMPLX(y[i] - m, 0.0); - #elif defined(_MSC_VER) - cplx tmp = { y[i] - m, 0.0 }; - F[i] = tmp; - #endif - } - for (int i = size; i < nFFT; i++) { - - #if defined(__GNUC__) || defined(__GNUG__) - F[i] = CMPLX(0.0, 0.0); - #elif defined(_MSC_VER) - cplx tmp = { 0.0, 0.0 }; - F[i] = tmp; - #endif - } - //size = nFFT; - - twiddles(tw, nFFT); - fft(F, nFFT, tw); - dot_multiply(F, F, nFFT); - fft(F, nFFT, tw); - cplx divisor = F[0]; - for (int i = 0; i < nFFT; i++) { - F[i] = _Cdivcc(F[i], divisor); // F[i] / divisor; - } - - double * out = malloc(nFFT * 2 * sizeof(out)); - for (int i = 0; i < nFFT; i++) { - out[i] = creal(F[i]); - } - free(F); - free(tw); - return out; -} - -int co_firstzero(const double y[], const int size, const int maxtau) -{ - - //double * autocorrs = malloc(size * sizeof * autocorrs); - //autocorrs = co_autocorrs(y, size); - - double * autocorrs = co_autocorrs(y, size); - - int zerocrossind = 0; - while(autocorrs[zerocrossind] > 0 && zerocrossind < maxtau) - { - zerocrossind += 1; - } - - free(autocorrs); - return zerocrossind; - -} - -double CO_f1ecac(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return 0; - } - } - - // compute autocorrelations - double * autocorrs = co_autocorrs(y, size); - - // threshold to cross - double thresh = 1.0/exp(1); - - double out = (double)size; - for(int i = 0; i < size-2; i++){ - // printf("i=%d autocorrs_i=%1.3f\n", i, autocorrs[i]); - if ( autocorrs[i+1] < thresh ){ - double m = autocorrs[i+1] - autocorrs[i]; - double dy = thresh - autocorrs[i]; - double dx = dy/m; - out = ((double)i) + dx; - // printf("thresh=%1.3f AC(i)=%1.3f AC(i-1)=%1.3f m=%1.3f dy=%1.3f dx=%1.3f out=%1.3f\n", thresh, autocorrs[i], autocorrs[i-1], m, dy, dx, out); - free(autocorrs); - return out; - } - } - - free(autocorrs); - - return out; - -} - -double CO_Embed2_Basic_tau_incircle(const double y[], const int size, const double radius, const int tau) -{ - int tauIntern = 0; - - if(tau < 0) - { - tauIntern = co_firstzero(y, size, size); - } - else{ - tauIntern = tau; - } - - double insidecount = 0; - for(int i = 0; i < size-tauIntern; i++) - { - if(y[i]*y[i] + y[i+tauIntern]*y[i+tauIntern] < radius) - { - insidecount += 1; - } - } - - return insidecount/(size-tauIntern); -} - -double CO_Embed2_Dist_tau_d_expfit_meandiff(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - int tau = co_firstzero(y, size, size); - - //printf("co_firstzero ran\n"); - - if (tau > (double)size/10){ - tau = floor((double)size/10); - } - //printf("tau = %i\n", tau); - - double * d = malloc((size-tau) * sizeof(double)); - for(int i = 0; i < size-tau-1; i++) - { - - d[i] = sqrt((y[i+1]-y[i])*(y[i+1]-y[i]) + (y[i+tau]-y[i+tau+1])*(y[i+tau]-y[i+tau+1])); - - //printf("d[%i]: %1.3f\n", i, d[i]); - if (isnan(d[i])){ - free(d); - return NAN; - } - - /* - if(i<100) - printf("%i, y[i]=%1.3f, y[i+1]=%1.3f, y[i+tau]=%1.3f, y[i+tau+1]=%1.3f, d[i]: %1.3f\n", i, y[i], y[i+1], y[i+tau], y[i+tau+1], d[i]); - */ - } - - //printf("embedding finished\n"); - - // mean for exponential fit - double l = mean(d, size-tau-1); - - // count histogram bin contents - /* - int * histCounts; - double * binEdges; - int nBins = histcounts(d, size-tau-1, -1, &histCounts, &binEdges); - */ - - int nBins = num_bins_auto(d, size-tau-1); - if (nBins == 0){ - free(d); - return 0; - } - int * histCounts = malloc(nBins * sizeof(double)); - double * binEdges = malloc((nBins + 1) * sizeof(double)); - histcounts_preallocated(d, size-tau-1, nBins, histCounts, binEdges); - - //printf("histcount ran\n"); - - // normalise to probability - double * histCountsNorm = malloc(nBins * sizeof(double)); - for(int i = 0; i < nBins; i++){ - //printf("histCounts %i: %i\n", i, histCounts[i]); - histCountsNorm[i] = (double)histCounts[i]/(double)(size-tau-1); - //printf("histCounts norm %i: %1.3f\n", i, histCountsNorm[i]); - } - - /* - for(int i = 0; i < nBins; i++){ - printf("histCounts[%i] = %i\n", i, histCounts[i]); - } - for(int i = 0; i < nBins; i++){ - printf("histCountsNorm[%i] = %1.3f\n", i, histCountsNorm[i]); - } - for(int i = 0; i < nBins+1; i++){ - printf("binEdges[%i] = %1.3f\n", i, binEdges[i]); - } - */ - - - //printf("histcounts normed\n"); - - double * d_expfit_diff = malloc(nBins * sizeof(double)); - for(int i = 0; i < nBins; i++){ - double expf = exp(-(binEdges[i] + binEdges[i+1])*0.5/l)/l; - if (expf < 0){ - expf = 0; - } - d_expfit_diff[i] = fabs(histCountsNorm[i]-expf); - //printf("d_expfit_diff %i: %1.3f\n", i, d_expfit_diff[i]); - } - - double out = mean(d_expfit_diff, nBins); - - //printf("out = %1.6f\n", out); - //printf("reached free statements\n"); - - // arrays created dynamically in function histcounts - free(d); - free(d_expfit_diff); - free(binEdges); - free(histCountsNorm); - free(histCounts); - - return out; - -} - -int CO_FirstMin_ac(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return 0; - } - } - - double * autocorrs = co_autocorrs(y, size); - - int minInd = size; - for(int i = 1; i < size-1; i++) - { - if(autocorrs[i] < autocorrs[i-1] && autocorrs[i] < autocorrs[i+1]) - { - minInd = i; - break; - } - } - - free(autocorrs); - - return minInd; - -} - -double CO_trev_1_num(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - int tau = 1; - - double * diffTemp = malloc((size-1) * sizeof * diffTemp); - - for(int i = 0; i < size-tau; i++) - { - diffTemp[i] = pow(y[i+1] - y[i],3); - } - - double out; - - out = mean(diffTemp, size-tau); - - free(diffTemp); - - return out; -} - -#define tau 2 -#define numBins 5 - -double CO_HistogramAMI_even_2_5(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - //const int tau = 2; - //const int numBins = 5; - - double * y1 = malloc((size-tau) * sizeof(double)); - double * y2 = malloc((size-tau) * sizeof(double)); - - for(int i = 0; i < size-tau; i++){ - y1[i] = y[i]; - y2[i] = y[i+tau]; - } - - // set bin edges - const double maxValue = max_(y, size); - const double minValue = min_(y, size); - - double binStep = (maxValue - minValue + 0.2)/5; - //double binEdges[numBins+1] = {0}; - double binEdges[5+1] = {0}; - for(int i = 0; i < numBins+1; i++){ - binEdges[i] = minValue + binStep*i - 0.1; - // printf("binEdges[%i] = %1.3f\n", i, binEdges[i]); - } - - - // count histogram bin contents - int * bins1; - bins1 = histbinassign(y1, size-tau, binEdges, numBins+1); - - int * bins2; - bins2 = histbinassign(y2, size-tau, binEdges, numBins+1); - - /* - // debug - for(int i = 0; i < size-tau; i++){ - printf("bins1[%i] = %i, bins2[%i] = %i\n", i, bins1[i], i, bins2[i]); - } - */ - - // joint - double * bins12 = malloc((size-tau) * sizeof(double)); - //double binEdges12[(numBins + 1) * (numBins + 1)] = {0}; - double binEdges12[(5 + 1) * (5 + 1)] = {0}; - - for(int i = 0; i < size-tau; i++){ - bins12[i] = (bins1[i]-1)*(numBins+1) + bins2[i]; - // printf("bins12[%i] = %1.3f\n", i, bins12[i]); - } - - for(int i = 0; i < (numBins+1)*(numBins+1); i++){ - binEdges12[i] = i+1; - // printf("binEdges12[%i] = %1.3f\n", i, binEdges12[i]); - } - - // fancy solution for joint histogram here - int * jointHistLinear; - jointHistLinear = histcount_edges(bins12, size-tau, binEdges12, (numBins + 1) * (numBins + 1)); - - /* - // debug - for(int i = 0; i < (numBins+1)*(numBins+1); i++){ - printf("jointHistLinear[%i] = %i\n", i, jointHistLinear[i]); - } - */ - - // transfer to 2D histogram (no last bin, as in original implementation) - double pij[numBins][numBins]; - int sumBins = 0; - for(int i = 0; i < numBins; i++){ - for(int j = 0; j < numBins; j++){ - pij[j][i] = jointHistLinear[i*(numBins+1)+j]; - - // printf("pij[%i][%i]=%1.3f\n", i, j, pij[i][j]); - - sumBins += pij[j][i]; - } - } - - // normalise - for(int i = 0; i < numBins; i++){ - for(int j = 0; j < numBins; j++){ - pij[j][i] /= sumBins; - } - } - - // marginals - //double pi[numBins] = {0}; - double pi[5] = {0}; - //double pj[numBins] = {0}; - double pj[5] = {0}; - for(int i = 0; i < numBins; i++){ - for(int j = 0; j < numBins; j++){ - pi[i] += pij[i][j]; - pj[j] += pij[i][j]; - // printf("pij[%i][%i]=%1.3f, pi[%i]=%1.3f, pj[%i]=%1.3f\n", i, j, pij[i][j], i, pi[i], j, pj[j]); - } - } - - /* - // debug - for(int i = 0; i < numBins; i++){ - printf("pi[%i]=%1.3f, pj[%i]=%1.3f\n", i, pi[i], i, pj[i]); - } - */ - - // mutual information - double ami = 0; - for(int i = 0; i < numBins; i++){ - for(int j = 0; j < numBins; j++){ - if(pij[i][j] > 0){ - //printf("pij[%i][%i]=%1.3f, pi[%i]=%1.3f, pj[%i]=%1.3f, logarg=, %1.3f, log(...)=%1.3f\n", - // i, j, pij[i][j], i, pi[i], j, pj[j], pij[i][j]/(pi[i]*pj[j]), log(pij[i][j]/(pi[i]*pj[j]))); - ami += pij[i][j] * log(pij[i][j]/(pj[j]*pi[i])); - } - } - } - - free(bins1); - free(bins2); - free(jointHistLinear); - - free(y1); - free(y2); - free(bins12); - - return ami; -} diff --git a/src/C/CO_AutoCorr.h b/src/C/CO_AutoCorr.h deleted file mode 100644 index 34bbd28..0000000 --- a/src/C/CO_AutoCorr.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef CO_AUTOCORR_H -#define CO_AUTOCORR_H - -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -#include -#include -#include -#include -#include "stats.h" -#include "fft.h" - -extern int nextpow2(int n); -extern void dot_multiply(cplx a[], cplx b[], int size); -extern double * CO_AutoCorr(const double y[], const int size, const int tau[], const int tau_size); -extern double * co_autocorrs(const double y[], const int size); -extern int co_firstzero(const double y[], const int size, const int maxtau); -extern double CO_Embed2_Basic_tau_incircle(const double y[], const int size, const double radius, const int tau); -extern double CO_Embed2_Dist_tau_d_expfit_meandiff(const double y[], const int size); -extern int CO_FirstMin_ac(const double y[], const int size); -extern double CO_trev_1_num(const double y[], const int size); -extern double CO_f1ecac(const double y[], const int size); -extern double CO_HistogramAMI_even_2_5(const double y[], const int size); - -#endif diff --git a/src/C/DN_HistogramMode_10.c b/src/C/DN_HistogramMode_10.c deleted file mode 100644 index 717a394..0000000 --- a/src/C/DN_HistogramMode_10.c +++ /dev/null @@ -1,113 +0,0 @@ -#include -#include -#include -#include - -#include "stats.h" -#include "histcounts.h" - -double DN_HistogramMode_10(const double y[], const int size) -{ - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - const int nBins = 10; - - int * histCounts; - double * binEdges; - - histcounts(y, size, nBins, &histCounts, &binEdges); - - double maxCount = 0; - int numMaxs = 1; - double out = 0;; - for(int i = 0; i < nBins; i++) - { - // printf("binInd=%i, binCount=%i, binEdge=%1.3f \n", i, histCounts[i], binEdges[i]); - - if (histCounts[i] > maxCount) - { - maxCount = histCounts[i]; - numMaxs = 1; - out = (binEdges[i] + binEdges[i+1])*0.5; - } - else if (histCounts[i] == maxCount){ - - numMaxs += 1; - out += (binEdges[i] + binEdges[i+1])*0.5; - } - } - out = out/numMaxs; - - // arrays created dynamically in function histcounts - free(histCounts); - free(binEdges); - - return out; -} - -/* - double DN_HistogramMode_10(double y[], int size) - { - - double min = DBL_MAX, max=-DBL_MAX; - for(int i = 0; i < size; i++) - { - if (y[i] < min) - { - min = y[i]; - } - if (y[i] > max) - { - max = y[i]; - } - } - - double binStep = (max - min)/10; - - // fprintf(stdout, "min=%f, max=%f, binStep=%f \n", min, max, binStep); - - int histCounts[10] = {0}; - for(int i = 0; i < size; i++) - { - int binsLeft = 10; - int lowerInd = 0, upperInd = 10; - while(binsLeft > 1) - { - int limitInd = (upperInd - lowerInd)/2 + lowerInd; - double limit = limitInd * binStep + min; - - if (y[i] < limit) - { - upperInd = limitInd; - } - else - { - lowerInd = limitInd; - } - binsLeft = upperInd - lowerInd; - } - histCounts[lowerInd] += 1; - } - - double maxCount = 0; - int maxCountInd = 0; - for(int i = 0; i < 10; i++) - { - // fprintf(stdout, "binInd=%i, binCount=%i \n", i, histCounts[i]); - - if (histCounts[i] > maxCount) - { - maxCountInd = i; - maxCount = histCounts[i]; - } - } - return binStep*(maxCountInd+0.5) + min; - } - */ diff --git a/src/C/DN_HistogramMode_10.h b/src/C/DN_HistogramMode_10.h deleted file mode 100644 index 9107615..0000000 --- a/src/C/DN_HistogramMode_10.h +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef DN_HISTOGRAMMODE_10 -#define DN_HISTOGRAMMODE_10 -#include -#include -#include "stats.h" - -extern double DN_HistogramMode_10(const double y[], const int size); - -#endif diff --git a/src/C/DN_HistogramMode_5.c b/src/C/DN_HistogramMode_5.c deleted file mode 100644 index 0e60293..0000000 --- a/src/C/DN_HistogramMode_5.c +++ /dev/null @@ -1,123 +0,0 @@ -#include -#include -#include -#include -#include "stats.h" -#include "histcounts.h" - -double DN_HistogramMode_5(const double y[], const int size) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - const int nBins = 5; - - int * histCounts; - double * binEdges; - - histcounts(y, size, nBins, &histCounts, &binEdges); - - /* - for(int i = 0; i < nBins; i++){ - printf("histCounts[%i] = %i\n", i, histCounts[i]); - } - for(int i = 0; i < nBins+1; i++){ - printf("binEdges[%i] = %1.3f\n", i, binEdges[i]); - } - */ - - double maxCount = 0; - int numMaxs = 1; - double out = 0;; - for(int i = 0; i < nBins; i++) - { - // printf("binInd=%i, binCount=%i, binEdge=%1.3f \n", i, histCounts[i], binEdges[i]); - - if (histCounts[i] > maxCount) - { - maxCount = histCounts[i]; - numMaxs = 1; - out = (binEdges[i] + binEdges[i+1])*0.5; - } - else if (histCounts[i] == maxCount){ - - numMaxs += 1; - out += (binEdges[i] + binEdges[i+1])*0.5; - } - } - out = out/numMaxs; - - // arrays created dynamically in function histcounts - free(histCounts); - free(binEdges); - - return out; -} - -/* -double DN_HistogramMode_5(double y[], int size) -{ - - double min = DBL_MAX, max=-DBL_MAX; - for(int i = 0; i < size; i++) - { - if (y[i] < min) - { - min = y[i]; - } - if (y[i] > max) - { - max = y[i]; - } - } - - double binStep = (max - min)/5; - - // fprintf(stdout, "min=%f, max=%f, binStep=%f \n", min, max, binStep); - - int histCounts[5] = {0}; - for(int i = 0; i < size; i++) - { - int binsLeft = 5; - int lowerInd = 0, upperInd = 10; - while(binsLeft > 1) - { - int limitInd = (upperInd - lowerInd)/2 + lowerInd; - double limit = limitInd * binStep + min; - - if (y[i] < limit) - { - upperInd = limitInd; - } - else - { - lowerInd = limitInd; - } - binsLeft = upperInd - lowerInd; - } - histCounts[lowerInd] += 1; - } - - double maxCount = 0; - int maxCountInd = 0; - for(int i = 0; i < 5; i++) - { - // fprintf(stdout, "binInd=%i, binCount=%i \n", i, histCounts[i]); - - if (histCounts[i] > maxCount) - { - maxCountInd = i; - maxCount = histCounts[i]; - } - } - return binStep*(maxCountInd+0.5) + min; -} - - */ diff --git a/src/C/DN_HistogramMode_5.h b/src/C/DN_HistogramMode_5.h deleted file mode 100644 index 9f24c6a..0000000 --- a/src/C/DN_HistogramMode_5.h +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef DN_HISTOGRAMMODE_5 -#define DN_HISTOGRAMMODE_5 -#include -#include -#include "stats.h" - -extern double DN_HistogramMode_5(const double y[], const int size); - -#endif diff --git a/src/C/DN_Mean.c b/src/C/DN_Mean.c deleted file mode 100644 index f2e75ea..0000000 --- a/src/C/DN_Mean.c +++ /dev/null @@ -1,11 +0,0 @@ -#include - -double DN_Mean(const double a[], const int size) -{ - double m = 0.0; - for (int i = 0; i < size; i++) { - m += a[i]; - } - m /= size; - return m; -} diff --git a/src/C/DN_Mean.h b/src/C/DN_Mean.h deleted file mode 100644 index a39d135..0000000 --- a/src/C/DN_Mean.h +++ /dev/null @@ -1,12 +0,0 @@ -// -// Created by Trent Henderson 27 September 2021 -// - -#ifndef DN_MEAN -#define DN_MEAN - -#include - -extern double DN_Mean(const double a[], const int size); - -#endif /* DN_MEAN */ \ No newline at end of file diff --git a/src/C/DN_OutlierInclude.c b/src/C/DN_OutlierInclude.c deleted file mode 100644 index 1478eb9..0000000 --- a/src/C/DN_OutlierInclude.c +++ /dev/null @@ -1,213 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include "stats.h" - -double DN_OutlierInclude_np_001_mdrmd(const double y[], const int size, const int sign) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - double inc = 0.01; - int tot = 0; - double * yWork = malloc(size * sizeof(double)); - - // apply sign and check constant time series - int constantFlag = 1; - for(int i = 0; i < size; i++) - { - if(y[i] != y[0]) - { - constantFlag = 0; - } - - // apply sign, save in new variable - yWork[i] = sign*y[i]; - - // count pos/ negs - if(yWork[i] >= 0){ - tot += 1; - } - - } - if(constantFlag) return 0; // if constant, return 0 - - // find maximum (or minimum, depending on sign) - double maxVal = max_(yWork, size); - - // maximum value too small? return 0 - if(maxVal < inc){ - return 0; - } - - int nThresh = maxVal/inc + 1; - - // save the indices where y > threshold - double * r = malloc(size * sizeof * r); - - // save the median over indices with absolute value > threshold - double * msDti1 = malloc(nThresh * sizeof(double)); - double * msDti3 = malloc(nThresh * sizeof(double)); - double * msDti4 = malloc(nThresh * sizeof(double)); - - for(int j = 0; j < nThresh; j++) - { - //printf("j=%i, thr=%1.3f\n", j, j*inc); - - int highSize = 0; - - for(int i = 0; i < size; i++) - { - if(yWork[i] >= j*inc) - { - r[highSize] = i+1; - //printf("r[%i]=%1.f \n", highSize, r[highSize]); - highSize += 1; - } - } - - // intervals between high-values - double * Dt_exc = malloc(highSize * sizeof(double)); - - for(int i = 0; i < highSize-1; i++) - { - //printf("i=%i, r[i+1]=%1.f, r[i]=%1.f \n", i, r[i+1], r[i]); - Dt_exc[i] = r[i+1] - r[i]; - } - - /* - // median - double medianOut; - medianOut = median(r, highSize); - */ - - msDti1[j] = mean(Dt_exc, highSize-1); - msDti3[j] = (highSize-1)*100.0/tot; - msDti4[j] = median(r, highSize) / ((double)size/2) - 1; - - //printf("msDti1[%i] = %1.3f, msDti13[%i] = %1.3f, msDti4[%i] = %1.3f\n", - // j, msDti1[j], j, msDti3[j], j, msDti4[j]); - - free(Dt_exc); - - } - - int trimthr = 2; - int mj = 0; - int fbi = nThresh-1; - for(int i = 0; i < nThresh; i ++) - { - if (msDti3[i] > trimthr) - { - mj = i; - } - if (isnan(msDti1[nThresh-1-i])) - { - fbi = nThresh-1-i; - } - } - - double outputScalar; - int trimLimit = mj < fbi ? mj : fbi; - outputScalar = median(msDti4, trimLimit+1); - - free(r); - free(yWork); - free(msDti1); - free(msDti3); - free(msDti4); - - return outputScalar; -} - -double DN_OutlierInclude_p_001_mdrmd(const double y[], const int size) -{ - return DN_OutlierInclude_np_001_mdrmd(y, size, 1.0); -} - -double DN_OutlierInclude_n_001_mdrmd(const double y[], const int size) -{ - return DN_OutlierInclude_np_001_mdrmd(y, size, -1.0); -} - -double DN_OutlierInclude_abs_001(const double y[], const int size) -{ - double inc = 0.01; - double maxAbs = 0; - double * yAbs = malloc(size * sizeof * yAbs); - - for(int i = 0; i < size; i++) - { - // yAbs[i] = (y[i] > 0) ? y[i] : -y[i]; - yAbs[i] = (y[i] > 0) ? y[i] : -y[i]; - - if(yAbs[i] > maxAbs) - { - maxAbs = yAbs[i]; - } - } - - int nThresh = maxAbs/inc + 1; - - printf("nThresh = %i\n", nThresh); - - // save the indices where y > threshold - double * highInds = malloc(size * sizeof * highInds); - - // save the median over indices with absolute value > threshold - double * msDti3 = malloc(nThresh * sizeof * msDti3); - double * msDti4 = malloc(nThresh * sizeof * msDti4); - - for(int j = 0; j < nThresh; j++) - { - int highSize = 0; - - for(int i = 0; i < size; i++) - { - if(yAbs[i] >= j*inc) - { - // fprintf(stdout, "%i, ", i); - - highInds[highSize] = i; - highSize += 1; - } - } - - // median - double medianOut; - medianOut = median(highInds, highSize); - - msDti3[j] = (highSize-1)*100.0/size; - msDti4[j] = medianOut / (size/2) - 1; - - } - - int trimthr = 2; - int mj = 0; - for(int i = 0; i < nThresh; i ++) - { - if (msDti3[i] > trimthr) - { - mj = i; - } - } - - double outputScalar; - outputScalar = median(msDti4, mj); - - free(highInds); - free(yAbs); - free(msDti4); - - return outputScalar; -} diff --git a/src/C/DN_OutlierInclude.h b/src/C/DN_OutlierInclude.h deleted file mode 100644 index c5559a7..0000000 --- a/src/C/DN_OutlierInclude.h +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef DN_OUTLIERINCLUDE_ABS_001 -#define DN_OUTLIERINCLUDE_ABS_001 -#include -#include -#include -#include -#include "stats.h" - -extern double DN_OutlierInclude_abs_001(const double y[], const int size); -extern double DN_OutlierInclude_np_001_mdrmd(const double y[], const int size, const int sign); -extern double DN_OutlierInclude_p_001_mdrmd(const double y[], const int size); -extern double DN_OutlierInclude_n_001_mdrmd(const double y[], const int size); - -#endif diff --git a/src/C/DN_Spread_Std.c b/src/C/DN_Spread_Std.c deleted file mode 100644 index 5f5f035..0000000 --- a/src/C/DN_Spread_Std.c +++ /dev/null @@ -1,13 +0,0 @@ -#include -#include "stats.h" - -double DN_Spread_Std(const double a[], const int size) -{ - double m = mean(a, size); - double sd = 0.0; - for (int i = 0; i < size; i++) { - sd += pow(a[i] - m, 2); - } - sd = sqrt(sd / (size - 1)); - return sd; -} diff --git a/src/C/DN_Spread_Std.h b/src/C/DN_Spread_Std.h deleted file mode 100644 index c150e14..0000000 --- a/src/C/DN_Spread_Std.h +++ /dev/null @@ -1,12 +0,0 @@ -// -// Created by Trent Henderson 27 September 2021 -// - -#ifndef DN_SPREADSTD -#define DN_SPREADSTD - -#include - -extern double DN_Spread_Std(const double a[], const int size); - -#endif /* DN_SPREADSTD */ \ No newline at end of file diff --git a/src/C/FC_LocalSimple.c b/src/C/FC_LocalSimple.c deleted file mode 100644 index 111db09..0000000 --- a/src/C/FC_LocalSimple.c +++ /dev/null @@ -1,183 +0,0 @@ -#include -#include -#include "stats.h" -#include "CO_AutoCorr.h" - -static void abs_diff(const double a[], const int size, double b[]) -{ - for (int i = 1; i < size; i++) { - b[i - 1] = fabs(a[i] - a[i - 1]); - } -} - -double fc_local_simple(const double y[], const int size, const int train_length) -{ - double * y1 = malloc((size - 1) * sizeof *y1); - abs_diff(y, size, y1); - double m = mean(y1, size - 1); - free(y1); - return m; -} - -double FC_LocalSimple_mean_tauresrat(const double y[], const int size, const int train_length) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - if(size <= train_length) - { - return NAN; - } - - double * res = malloc((size - train_length) * sizeof *res); - - for (int i = 0; i < size - train_length; i++) - { - double yest = 0; - for (int j = 0; j < train_length; j++) - { - yest += y[i+j]; - - } - yest /= train_length; - - res[i] = y[i+train_length] - yest; - } - - double resAC1stZ = co_firstzero(res, size - train_length, size - train_length); - double yAC1stZ = co_firstzero(y, size, size); - double output = resAC1stZ/yAC1stZ; - - free(res); - return output; - -} - -double FC_LocalSimple_mean_stderr(const double y[], const int size, const int train_length) -{ - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - if(size <= train_length) - { - return NAN; - } - - double * res = malloc((size - train_length) * sizeof *res); - - for (int i = 0; i < size - train_length; i++) - { - double yest = 0; - for (int j = 0; j < train_length; j++) - { - yest += y[i+j]; - - } - yest /= train_length; - - res[i] = y[i+train_length] - yest; - } - - double output = stddev(res, size - train_length); - - free(res); - return output; - -} - -double FC_LocalSimple_mean3_stderr(const double y[], const int size) -{ - return FC_LocalSimple_mean_stderr(y, size, 3); -} - -double FC_LocalSimple_mean1_tauresrat(const double y[], const int size){ - return FC_LocalSimple_mean_tauresrat(y, size, 1); -} - -double FC_LocalSimple_mean_taures(const double y[], const int size, const int train_length) -{ - if(size <= train_length) - { - return NAN; - } - - double * res = malloc((size - train_length) * sizeof *res); - - // first z-score - // no, assume ts is z-scored!! - //zscore_norm(y, size); - - for (int i = 0; i < size - train_length; i++) - { - double yest = 0; - for (int j = 0; j < train_length; j++) - { - yest += y[i+j]; - - } - yest /= train_length; - - res[i] = y[i+train_length] - yest; - } - - int output = co_firstzero(res, size - train_length, size - train_length); - - free(res); - return output; - -} - -double FC_LocalSimple_lfit_taures(const double y[], const int size) -{ - // set tau from first AC zero crossing - int train_length = co_firstzero(y, size, size); - - if(size <= train_length) - { - return NAN; - } - - double * xReg = malloc(train_length * sizeof * xReg); - // double * yReg = malloc(train_length * sizeof * yReg); - for(int i = 1; i < train_length+1; i++) - { - xReg[i-1] = i; - } - - double * res = malloc((size - train_length) * sizeof *res); - - double m = 0.0, b = 0.0; - - for (int i = 0; i < size - train_length; i++) - { - linreg(train_length, xReg, y+i, &m, &b); - - // fprintf(stdout, "i=%i, m=%f, b=%f\n", i, m, b); - - res[i] = y[i+train_length] - (m * (train_length+1) + b); - } - - int output = co_firstzero(res, size - train_length, size - train_length); - - free(res); - free(xReg); - // free(yReg); - - return output; - -} - - diff --git a/src/C/FC_LocalSimple.h b/src/C/FC_LocalSimple.h deleted file mode 100644 index d821f79..0000000 --- a/src/C/FC_LocalSimple.h +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef FC_LOCALSIMPLE_H -#define FC_LOCALSIMPLE_H -#include -#include -#include "stats.h" -#include "CO_AutoCorr.h" - -extern double fc_local_simple(const double y[], const int size, const int train_length); -extern double FC_LocalSimple_mean_taures(const double y[], const int size, const int train_length); -extern double FC_LocalSimple_lfit_taures(const double y[], const int size); -extern double FC_LocalSimple_mean_tauresrat(const double y[], const int size, const int train_length); -extern double FC_LocalSimple_mean1_tauresrat(const double y[], const int size); -extern double FC_LocalSimple_mean_stderr(const double y[], const int size, const int train_length); -extern double FC_LocalSimple_mean3_stderr(const double y[], const int size); - -#endif diff --git a/src/C/IN_AutoMutualInfoStats.c b/src/C/IN_AutoMutualInfoStats.c deleted file mode 100644 index 7b09477..0000000 --- a/src/C/IN_AutoMutualInfoStats.c +++ /dev/null @@ -1,54 +0,0 @@ -// -// IN_AutoMutualInfoStats.c -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// -#include - -#include "IN_AutoMutualInfoStats.h" -#include "CO_AutoCorr.h" -#include "stats.h" - -double IN_AutoMutualInfoStats_40_gaussian_fmmi(const double y[], const int size) -{ - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - // maximum time delay - int tau = 40; - - // don't go above half the signal length - if(tau > ceil((double)size/2)){ - tau = ceil((double)size/2); - } - - // compute autocorrelations and compute automutual information - double * ami = malloc(size * sizeof(double)); - for(int i = 0; i < tau; i++){ - double ac = autocorr_lag(y,size, i+1); - ami[i] = -0.5 * log(1 - ac*ac); - // printf("ami[%i]=%1.7f\n", i, ami[i]); - } - - // find first minimum of automutual information - double fmmi = tau; - for(int i = 1; i < tau-1; i++){ - if(ami[i] < ami[i-1] && ami[i] < ami[i+1]){ - fmmi = i; - // printf("found minimum at %i\n", i); - break; - } - } - - free(ami); - - return fmmi; -} diff --git a/src/C/IN_AutoMutualInfoStats.h b/src/C/IN_AutoMutualInfoStats.h deleted file mode 100644 index def01ad..0000000 --- a/src/C/IN_AutoMutualInfoStats.h +++ /dev/null @@ -1,16 +0,0 @@ -// -// IN_AutoMutualInfoStats.h -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef IN_AutoMutualInfoStats_h -#define IN_AutoMutualInfoStats_h - -#include - -extern double IN_AutoMutualInfoStats_40_gaussian_fmmi(const double y[], const int size); - -#endif /* IN_AutoMutualInfoStats_h */ diff --git a/src/C/MD_hrv.c b/src/C/MD_hrv.c deleted file mode 100644 index dcf8424..0000000 --- a/src/C/MD_hrv.c +++ /dev/null @@ -1,40 +0,0 @@ -// -// MD_hrv.c -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#include "MD_hrv.h" -#include "stats.h" - -double MD_hrv_classic_pnn40(const double y[], const int size){ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - const int pNNx = 40; - - // compute diff - double * Dy = malloc((size-1) * sizeof(double)); - diff(y, size, Dy); - - double pnn40 = 0; - for(int i = 0; i < size-1; i++){ - if(fabs(Dy[i])*1000 > pNNx){ - pnn40 += 1; - } - } - - free(Dy); - - return pnn40/(size-1); -} - diff --git a/src/C/MD_hrv.h b/src/C/MD_hrv.h deleted file mode 100644 index 5588eb2..0000000 --- a/src/C/MD_hrv.h +++ /dev/null @@ -1,16 +0,0 @@ -// -// MD_hrv.h -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef MD_hrv_h -#define MD_hrv_h - -#include - -extern double MD_hrv_classic_pnn40(const double y[], const int size); - -#endif /* MD_hrv_h */ diff --git a/src/C/PD_PeriodicityWang.c b/src/C/PD_PeriodicityWang.c deleted file mode 100644 index 6ebf586..0000000 --- a/src/C/PD_PeriodicityWang.c +++ /dev/null @@ -1,137 +0,0 @@ -// -// PD_PeriodicityWang.c -// C_polished -// -// Created by Carl Henning Lubba on 28/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#include -#include - -#include "PD_PeriodicityWang.h" -#include "splinefit.h" -#include "stats.h" - -int PD_PeriodicityWang_th0_01(const double * y, const int size){ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return 0; - } - } - - const double th = 0.01; - - double * ySpline = malloc(size * sizeof(double)); - - // fit a spline with 3 nodes to the data - splinefit(y, size, ySpline); - - //printf("spline fit complete.\n"); - - // subtract spline from data to remove trend - double * ySub = malloc(size * sizeof(double)); - for(int i = 0; i < size; i++){ - ySub[i] = y[i] - ySpline[i]; - //printf("ySub[%i] = %1.5f\n", i, ySub[i]); - } - - // compute autocorrelations up to 1/3 of the length of the time series - int acmax = (int)ceil((double)size/3); - - double * acf = malloc(acmax*sizeof(double)); - for(int tau = 1; tau <= acmax; tau++){ - // correlation/ covariance the same, don't care for scaling (cov would be more efficient) - acf[tau-1] = autocov_lag(ySub, size, tau); - //printf("acf[%i] = %1.9f\n", tau-1, acf[tau-1]); - } - - //printf("ACF computed.\n"); - - // find troughts and peaks - double * troughs = malloc(acmax * sizeof(double)); - double * peaks = malloc(acmax * sizeof(double)); - int nTroughs = 0; - int nPeaks = 0; - double slopeIn = 0; - double slopeOut = 0; - for(int i = 1; i < acmax-1; i ++){ - slopeIn = acf[i] - acf[i-1]; - slopeOut = acf[i+1] - acf[i]; - - if(slopeIn < 0 && slopeOut > 0) - { - // printf("trough at %i\n", i); - troughs[nTroughs] = i; - nTroughs += 1; - } - else if(slopeIn > 0 && slopeOut < 0) - { - // printf("peak at %i\n", i); - peaks[nPeaks] = i; - nPeaks += 1; - } - } - - //printf("%i troughs and %i peaks found.\n", nTroughs, nPeaks); - - - // search through all peaks for one that meets the conditions: - // (a) a trough before it - // (b) difference between peak and trough is at least 0.01 - // (c) peak corresponds to positive correlation - int iPeak = 0; - double thePeak = 0; - int iTrough = 0; - double theTrough = 0; - - int out = 0; - - for(int i = 0; i < nPeaks; i++){ - iPeak = peaks[i]; - thePeak = acf[iPeak]; - - //printf("i=%i/%i, iPeak=%i, thePeak=%1.3f\n", i, nPeaks-1, iPeak, thePeak); - - // find trough before this peak - int j = -1; - while(troughs[j+1] < iPeak && j+1 < nTroughs){ - // printf("j=%i/%i, iTrough=%i, theTrough=%1.3f\n", j+1, nTroughs-1, (int)troughs[j+1], acf[(int)troughs[j+1]]); - j++; - } - if(j == -1) - continue; - - iTrough = troughs[j]; - theTrough = acf[iTrough]; - - // (a) should be implicit - - // (b) different between peak and trough it as least 0.01 - if(thePeak - theTrough < th) - continue; - - // (c) peak corresponds to positive correlation - if(thePeak < 0) - continue; - - // use this frequency that first fulfils all conditions. - out = iPeak; - break; - } - - //printf("Before freeing stuff.\n"); - - free(ySpline); - free(ySub); - free(acf); - free(troughs); - free(peaks); - - return out; - -} diff --git a/src/C/PD_PeriodicityWang.h b/src/C/PD_PeriodicityWang.h deleted file mode 100644 index 03d4dba..0000000 --- a/src/C/PD_PeriodicityWang.h +++ /dev/null @@ -1,16 +0,0 @@ -// -// PD_PeriodicityWang.h -// C_polished -// -// Created by Carl Henning Lubba on 28/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef PD_PeriodicityWang_h -#define PD_PeriodicityWang_h - -#include - -extern int PD_PeriodicityWang_th0_01(const double * y, const int size); - -#endif /* PD_PeriodicityWang_h */ diff --git a/src/C/SB_BinaryStats.c b/src/C/SB_BinaryStats.c deleted file mode 100644 index 58e2115..0000000 --- a/src/C/SB_BinaryStats.c +++ /dev/null @@ -1,91 +0,0 @@ -// -// SB_BinaryStats.c -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#include "SB_BinaryStats.h" -#include "stats.h" - -double SB_BinaryStats_diff_longstretch0(const double y[], const int size){ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - // binarize - int * yBin = malloc((size-1) * sizeof(int)); - for(int i = 0; i < size-1; i++){ - - double diffTemp = y[i+1] - y[i]; - yBin[i] = diffTemp < 0 ? 0 : 1; - - /* - if( i < 300) - printf("%i, y[i+1]=%1.3f, y[i]=%1.3f, yBin[i]=%i\n", i, y[i+1], y[i], yBin[i]); - */ - - } - - int maxstretch0 = 0; - int last1 = 0; - for(int i = 0; i < size-1; i++){ - if(yBin[i] == 1 || i == size-2){ - double stretch0 = i - last1; - if(stretch0 > maxstretch0){ - maxstretch0 = stretch0; - } - last1 = i; - } - } - - free(yBin); - - return maxstretch0; -} - -double SB_BinaryStats_mean_longstretch1(const double y[], const int size){ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - // binarize - int * yBin = malloc((size-1) * sizeof(int)); - double yMean = mean(y, size); - for(int i = 0; i < size-1; i++){ - - yBin[i] = (y[i] - yMean <= 0) ? 0 : 1; - //printf("yBin[%i]=%i\n", i, yBin[i]); - - } - - int maxstretch1 = 0; - int last1 = 0; - for(int i = 0; i < size-1; i++){ - if(yBin[i] == 0 || i == size-2){ - double stretch1 = i - last1; - if(stretch1 > maxstretch1){ - maxstretch1 = stretch1; - } - last1 = i; - } - - } - - free(yBin); - - return maxstretch1; -} diff --git a/src/C/SB_BinaryStats.h b/src/C/SB_BinaryStats.h deleted file mode 100644 index fb3ad88..0000000 --- a/src/C/SB_BinaryStats.h +++ /dev/null @@ -1,17 +0,0 @@ -// -// SB_BinaryStats.h -// C_polished -// -// Created by Carl Henning Lubba on 22/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef SB_BinaryStats_h -#define SB_BinaryStats_h - -#include - -extern double SB_BinaryStats_diff_longstretch0(const double y[], const int size); -extern double SB_BinaryStats_mean_longstretch1(const double y[], const int size); - -#endif /* SB_BinaryStats_h */ diff --git a/src/C/SB_CoarseGrain.c b/src/C/SB_CoarseGrain.c deleted file mode 100644 index 8805758..0000000 --- a/src/C/SB_CoarseGrain.c +++ /dev/null @@ -1,40 +0,0 @@ -#include -#include -#include -#include -#include "stats.h" -#include "helper_functions.h" - -void sb_coarsegrain(const double y[], const int size, const char how[], const int num_groups, int labels[]) -{ - int i, j; - if (strcmp(how, "quantile") == 1) { - fprintf(stdout, "ERROR in sb_coarsegrain: unknown coarse-graining method\n"); - exit(1); - } - - /* - for(int i = 0; i < size; i++){ - printf("yin coarsegrain[%i]=%1.4f\n", i, y[i]); - } - */ - - double * th = malloc((num_groups + 1) * 2 * sizeof(th)); - double * ls = malloc((num_groups + 1) * 2 * sizeof(th)); - linspace(0, 1, num_groups + 1, ls); - for (i = 0; i < num_groups + 1; i++) { - //double quant = quantile(y, size, ls[i]); - th[i] = quantile(y, size, ls[i]); - } - th[0] -= 1; - for (i = 0; i < num_groups; i++) { - for (j = 0; j < size; j++) { - if (y[j] > th[i] && y[j] <= th[i + 1]) { - labels[j] = i + 1; - } - } - } - - free(th); - free(ls); -} diff --git a/src/C/SB_CoarseGrain.h b/src/C/SB_CoarseGrain.h deleted file mode 100644 index cf9607e..0000000 --- a/src/C/SB_CoarseGrain.h +++ /dev/null @@ -1,12 +0,0 @@ -#ifndef SB_COARSEGRAIN_H -#define SB_COARSEGRAIN_H -#include -#include -#include -#include -#include "stats.h" -#include "helper_functions.h" - -extern void sb_coarsegrain(const double y[], const int size, const char how[], const int num_groups, int labels[]); - -#endif diff --git a/src/C/SB_MotifThree.c b/src/C/SB_MotifThree.c deleted file mode 100644 index 438a961..0000000 --- a/src/C/SB_MotifThree.c +++ /dev/null @@ -1,375 +0,0 @@ -#include -#include -#include -#include -#include "SB_CoarseGrain.h" -#include "helper_functions.h" - -double SB_MotifThree_quantile_hh(const double y[], const int size) -{ - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - int tmp_idx, r_idx; - int dynamic_idx; - int alphabet_size = 3; - int array_size; - int * yt = malloc(size * sizeof(yt)); // alphabetized array - double hh; // output - double * out = malloc(124 * sizeof(out)); // output array - - // transfer to alphabet - sb_coarsegrain(y, size, "quantile", 3, yt); - - // words of length 1 - array_size = alphabet_size; - int ** r1 = malloc(array_size * sizeof(*r1)); - int * sizes_r1 = malloc(array_size * sizeof(sizes_r1)); - double * out1 = malloc(array_size * sizeof(out1)); - for (int i = 0; i < alphabet_size; i++) { - r1[i] = malloc(size * sizeof(r1[i])); // probably can be rewritten - // using selfresizing array for memory efficiency. Time complexity - // should be comparable due to ammotization. - r_idx = 0; - sizes_r1[i] = 0; - for (int j = 0; j < size; j++) { - if (yt[j] == i + 1) { - r1[i][r_idx++] = j; - sizes_r1[i]++; - } - } - } - - // words of length 2 - array_size *= alphabet_size; - // removing last item if it is == max possible idx since later we are taking idx + 1 - // from yt - for (int i = 0; i < alphabet_size; i++) { - if (sizes_r1[i] != 0 && r1[i][sizes_r1[i] - 1] == size - 1) { - //int * tmp_ar = malloc((sizes_r1[i] - 1) * sizeof(tmp_ar)); - int* tmp_ar = malloc(sizes_r1[i] * sizeof(tmp_ar)); - subset(r1[i], tmp_ar, 0, sizes_r1[i]); - memcpy(r1[i], tmp_ar, (sizes_r1[i] - 1) * sizeof(tmp_ar)); - sizes_r1[i]--; - free(tmp_ar); - } - } - - /* - int *** r2 = malloc(array_size * sizeof(**r2)); - int ** sizes_r2 = malloc(array_size * sizeof(*sizes_r2)); - double ** out2 = malloc(array_size * sizeof(*out2)); - */ - int*** r2 = malloc(alphabet_size * sizeof(**r2)); - int** sizes_r2 = malloc(alphabet_size * sizeof(*sizes_r2)); - double** out2 = malloc(alphabet_size * sizeof(*out2)); - - - // allocate separately - for (int i = 0; i < alphabet_size; i++) { - r2[i] = malloc(alphabet_size * sizeof(*r2[i])); - sizes_r2[i] = malloc(alphabet_size * sizeof(*sizes_r2[i])); - //out2[i] = malloc(alphabet_size * sizeof(out2[i])); - out2[i] = malloc(alphabet_size * sizeof(**out2)); - for (int j = 0; j < alphabet_size; j++) { - r2[i][j] = malloc(size * sizeof(*r2[i][j])); - } - } - - // fill separately - for (int i = 0; i < alphabet_size; i++) { - // for (int i = 0; i < array_size; i++) { - //r2[i] = malloc(alphabet_size * sizeof(r2[i])); - //sizes_r2[i] = malloc(alphabet_size * sizeof(sizes_r2[i])); - //out2[i] = malloc(alphabet_size * sizeof(out2[i])); - for (int j = 0; j < alphabet_size; j++) { - //r2[i][j] = malloc(size * sizeof(r2[i][j])); - sizes_r2[i][j] = 0; - dynamic_idx = 0; //workaround as you can't just add elements to array - // like in python (list.append()) for example, so since for some k there will be no adding, - // you need to keep track of the idx at which elements will be inserted - for (int k = 0; k < sizes_r1[i]; k++) { - tmp_idx = yt[r1[i][k] + 1]; - if (tmp_idx == (j + 1)) { - r2[i][j][dynamic_idx++] = r1[i][k]; - sizes_r2[i][j]++; - // printf("dynamic_idx=%i, size = %i\n", dynamic_idx, size); - } - } - double tmp = (double)sizes_r2[i][j] / ((double)(size) - (double)(1.0)); - out2[i][j] = tmp; - } - } - - hh = 0.0; - for (int i = 0; i < alphabet_size; i++) { - hh += f_entropy(out2[i], alphabet_size); - } - - free(yt); - free(out); - free(out1); - - free(sizes_r1); - - // free nested array - for (int i = 0; i < alphabet_size; i++) { - free(r1[i]); - } - free(r1); - // free(sizes_r1); - - for (int i = 0; i < alphabet_size; i++) { - //for (int i = alphabet_size - 1; i >= 0; i--) { - - free(sizes_r2[i]); - free(out2[i]); - } - - //for (int i = alphabet_size-1; i >= 0 ; i--) { - for(int i = 0; i < alphabet_size; i++) { - for (int j = 0; j < alphabet_size; j++) { - free(r2[i][j]); - } - free(r2[i]); - } - - free(r2); - free(sizes_r2); - free(out2); - - - return hh; - -} - -double * sb_motifthree(const double y[], int size, const char how[]) -{ - int tmp_idx, r_idx, i, j, k, l, m, array_size; - int dynamic_idx; - int * tmp_ar; - int alphabet_size = 3; - int out_idx = 0; - int * yt = malloc(size * sizeof(yt)); - double tmp; - double * out = malloc(124 * sizeof(out)); // output array - if (strcmp(how, "quantile") == 0) { - sb_coarsegrain(y, size, how, alphabet_size, yt); - } else if (strcmp(how, "diffquant") == 0) { - double * diff_y = malloc((size - 1) * sizeof(diff_y)); - diff(y, size, diff_y); - sb_coarsegrain(diff_y, size, how, alphabet_size, yt); - size--; - } else { - fprintf(stdout, "ERROR in sb_motifthree: Unknown how method"); - exit(1); - } - - // words of length 1 - array_size = alphabet_size; - int ** r1 = malloc(array_size * sizeof(*r1)); - int * sizes_r1 = malloc(array_size * sizeof(sizes_r1)); - double * out1 = malloc(array_size * sizeof(out1)); - for (i = 0; i < array_size; i++) { - r1[i] = malloc(size * sizeof(r1[i])); // probably can be rewritten - // using selfresizing array for memory efficiency. Time complexity - // should be comparable due to ammotization. - r_idx = 0; - sizes_r1[i] = 0; - for (j = 0; j < size; j++) { - if (yt[j] == i + 1) { - r1[i][r_idx++] = j; - sizes_r1[i]++; - } - } - tmp = (double)sizes_r1[i] / size; - - out1[i] = tmp; - out[out_idx++] = tmp; - } - out[out_idx++] = f_entropy(out1, array_size); - - // words of length 2 - array_size *= alphabet_size; - // removing last item if it is == max possible idx since later we are taking idx + 1 - // from yt - for (i = 0; i < alphabet_size; i++) { - if (sizes_r1[i] != 0 && r1[i][sizes_r1[i] - 1] == size - 1) { - tmp_ar = malloc((sizes_r1[i] - 1) * sizeof(tmp_ar)); - subset(r1[i], tmp_ar, 0, sizes_r1[i]); - memcpy(r1[i], tmp_ar, (sizes_r1[i] - 1) * sizeof(tmp_ar)); - sizes_r1[i]--; - } - } - - int *** r2 = malloc(array_size * sizeof(**r2)); - int ** sizes_r2 = malloc(array_size * sizeof(*sizes_r2)); - double ** out2 = malloc(array_size * sizeof(*out2)); - for (i = 0; i < alphabet_size; i++) { - r2[i] = malloc(alphabet_size * sizeof(r2[i])); - sizes_r2[i] = malloc(alphabet_size * sizeof(sizes_r2[i])); - out2[i] = malloc(alphabet_size * sizeof(out2[i])); - for (j = 0; j < alphabet_size; j++) { - r2[i][j] = malloc(size * sizeof(r2[i][j])); - sizes_r2[i][j] = 0; - dynamic_idx = 0; //workaround as you can't just add elements to array - // like in python (list.append()) for example, so since for some k there will be no adding, - // you need to keep track of the idx at which elements will be inserted - for (k = 0; k < sizes_r1[i]; k++) { - tmp_idx = yt[r1[i][k] + 1]; - if (tmp_idx == (j + 1)) { - r2[i][j][dynamic_idx++] = r1[i][k]; - sizes_r2[i][j]++; - } - } - tmp = (double)sizes_r2[i][j] / (size - 1); - out2[i][j] = tmp; - out[out_idx++] = tmp; - } - } - tmp = 0.0; - for (i = 0; i < alphabet_size; i++) { - tmp += f_entropy(out2[i], alphabet_size); - } - out[out_idx++] = tmp; - - // words of length 3 - array_size *= alphabet_size; - for (i = 0; i < alphabet_size; i++) { - for (j = 0; j < alphabet_size; j++) { - if (sizes_r2[i][j] != 0 && r2[i][j][sizes_r2[i][j] - 1] == size - 2) { - subset(r2[i][j], tmp_ar, 0, sizes_r2[i][j]); - memcpy(r2[i][j], tmp_ar, (sizes_r2[i][j] - 1) * sizeof(tmp_ar)); - sizes_r2[i][j]--; - } - } - } - - int **** r3 = malloc(array_size * sizeof(***r3)); - int *** sizes_r3 = malloc(array_size * sizeof(**sizes_r3)); - double *** out3 = malloc(array_size * sizeof(**out3)); - for (i = 0; i < alphabet_size; i++) { - r3[i] = malloc(alphabet_size * sizeof(r3[i])); - sizes_r3[i] = malloc(alphabet_size * sizeof(sizes_r3[i])); - out3[i] = malloc(alphabet_size * sizeof(out3[i])); - for (j = 0; j < alphabet_size; j++) { - r3[i][j] = malloc(alphabet_size * sizeof(r3[i][j])); - sizes_r3[i][j] = malloc(alphabet_size * sizeof(sizes_r3[i][j])); - out3[i][j] = malloc(alphabet_size * sizeof(out3[i][j])); - for (k = 0; k < alphabet_size; k++) { - r3[i][j][k] = malloc(size * sizeof(r3[i][j][k])); - sizes_r3[i][j][k] = 0; - dynamic_idx = 0; - for (l = 0; l < sizes_r2[i][j]; l++) { - tmp_idx = yt[r2[i][j][l] + 2]; - if (tmp_idx == (k + 1)) { - r3[i][j][k][dynamic_idx++] = r2[i][j][l]; - sizes_r3[i][j][k]++; - } - } - tmp = (double)sizes_r3[i][j][k] / (size - 2); - out3[i][j][k] = tmp; - out[out_idx++] = tmp; - } - } - } - tmp = 0.0; - for (i = 0; i < alphabet_size; i++) { - for (j = 0; j < alphabet_size; j++) { - tmp += f_entropy(out3[i][j], alphabet_size); - } - } - out[out_idx++] = tmp; - - // words of length 4 - array_size *= alphabet_size; - for (i = 0; i < alphabet_size; i++) { - for (j = 0; j < alphabet_size; j++) { - for (k = 0; k < alphabet_size; k++) { - if (sizes_r3[i][j][k] != 0 && r3[i][j][k][sizes_r3[i][j][k] - 1] == size - 3) { - subset(r3[i][j][k], tmp_ar, 0, sizes_r3[i][j][k]); - memcpy(r3[i][j][k], tmp_ar, (sizes_r3[i][j][k] - 1) * sizeof(tmp_ar)); - sizes_r3[i][j][k]--; - } - } - } - } - - int ***** r4 = malloc(array_size * sizeof(****r4)); - // just an array of pointers of array of pointers of array of pointers - // of array of pointers of array of ints... We need to go deeper (c) - int **** sizes_r4 = malloc(array_size * sizeof(***sizes_r3)); - double **** out4 = malloc(array_size * sizeof(***out4)); - for (i = 0; i < alphabet_size; i++) { - r4[i] = malloc(alphabet_size * sizeof(r4[i])); - sizes_r4[i] = malloc(alphabet_size * sizeof(sizes_r4[i])); - out4[i] = malloc(alphabet_size * sizeof(out4[i])); - for (j = 0; j < alphabet_size; j++) { - r4[i][j] = malloc(alphabet_size * sizeof(r4[i][j])); - sizes_r4[i][j] = malloc(alphabet_size * sizeof(sizes_r4[i][j])); - out4[i][j] = malloc(alphabet_size * sizeof(out4[i][j])); - for (k = 0; k < alphabet_size; k++) { - r4[i][j][k] = malloc(alphabet_size * sizeof(r4[i][j][k])); - sizes_r4[i][j][k] = malloc(alphabet_size * sizeof(sizes_r4[i][j][k])); - out4[i][j][k] = malloc(alphabet_size * sizeof(out4[i][j][k])); - for (l = 0; l < alphabet_size; l++) { - r4[i][j][k][l] = malloc(size * sizeof(r4[i][j][k][l])); - sizes_r4[i][j][k][l] = 0; - dynamic_idx = 0; - for (m = 0; m < sizes_r3[i][j][k]; m++) { - tmp_idx = yt[r3[i][j][k][m] + 3]; - if (tmp_idx == l + 1) { - r4[i][j][k][l][dynamic_idx++] = r3[i][j][k][m]; - sizes_r4[i][j][k][l]++; - } - } - tmp = (double)sizes_r4[i][j][k][l] / (size - 3); - out4[i][j][k][l] = tmp; - out[out_idx++] = tmp; - } - } - } - } - tmp = 0.0; - for (i = 0; i < alphabet_size; i++) { - for (j = 0; j < alphabet_size; j++) { - for (k = 0; k < alphabet_size; k++) { - tmp += f_entropy(out4[i][j][k], alphabet_size); - } - } - } - out[out_idx++] = tmp; - - return out; -} - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/C/SB_MotifThree.h b/src/C/SB_MotifThree.h deleted file mode 100644 index 2a471d7..0000000 --- a/src/C/SB_MotifThree.h +++ /dev/null @@ -1,12 +0,0 @@ -#ifndef SB_MOTIFTHREE_H -#define SB_MOTIFTHREE_H -#include -#include -#include -#include "SB_CoarseGrain.h" -#include "helper_functions.h" - -extern double SB_MotifThree_quantile_hh(const double y[], const int size); -extern double * sb_motifthree(const double y[], int size, const char how[]); - -#endif diff --git a/src/C/SB_TransitionMatrix.c b/src/C/SB_TransitionMatrix.c deleted file mode 100644 index 037df09..0000000 --- a/src/C/SB_TransitionMatrix.c +++ /dev/null @@ -1,163 +0,0 @@ -// -// SB_TransitionMatrix.c -// -// -// Created by Carl Henning Lubba on 23/09/2018. -// - -#include "SB_TransitionMatrix.h" -#include "butterworth.h" -#include "CO_AutoCorr.h" -#include "SB_CoarseGrain.h" -#include "stats.h" - -double SB_TransitionMatrix_3ac_sumdiagcov(const double y[], const int size) -{ - - // NaN and const check - int constant = 1; - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - if(y[i] != y[0]){ - constant = 0; - } - } - if (constant){ - return NAN; - } - - const int numGroups = 3; - - int tau = co_firstzero(y, size, size); - - double * yFilt = malloc(size * sizeof(double)); - - // sometimes causes problems in filt!!! needs fixing. - /* - if(tau > 1){ - butterworthFilter(y, size, 4, 0.8/tau, yFilt); - } - */ - - for(int i = 0; i < size; i++){ - yFilt[i] = y[i]; - } - - /* - for(int i = 0; i < size; i++){ - printf("yFilt[%i]=%1.4f\n", i, yFilt[i]); - } - */ - - int nDown = (size-1)/tau+1; - double * yDown = malloc(nDown * sizeof(double)); - - for(int i = 0; i < nDown; i++){ - yDown[i] = yFilt[i*tau]; - } - - /* - for(int i = 0; i < nDown; i++){ - printf("yDown[%i]=%1.4f\n", i, yDown[i]); - } - */ - - - // transfer to alphabet - int * yCG = malloc(nDown * sizeof(double)); - sb_coarsegrain(yDown, nDown, "quantile", numGroups, yCG); - - /* - for(int i = 0; i < nDown; i++){ - printf("yCG[%i]=%i\n", i, yCG[i]); - } - */ - - - double T[3][3]; - for(int i = 0; i < numGroups; i++){ - for(int j = 0; j < numGroups; j++){ - T[i][j] = 0; - } - } - - // more efficient way of doing the below - for(int j = 0; j < nDown-1; j++){ - T[yCG[j]-1][yCG[j+1]-1] += 1; - } - - /* - for(int i = 0; i < numGroups; i++){ - for(int j = 0; j < numGroups; j++){ - printf("%1.f, ", T[i][j]); - } - printf("\n"); - } - */ - - /* - for(int i = 0; i < numGroups; i++){ - for(int j = 0; j < nDown-1; j++){ - if(yCG[j] == i+1){ - T[i][yCG[j+1]-1] += 1; - } - } - } - */ - - for(int i = 0; i < numGroups; i++){ - for(int j = 0; j < numGroups; j++){ - T[i][j] /= (nDown-1); - // printf("T(%i, %i) = %1.3f\n", i, j, T[i][j]); - - } - } - - double column1[3] = {0}; - double column2[3] = {0}; - double column3[3] = {0}; - - for(int i = 0; i < numGroups; i++){ - column1[i] = T[i][0]; - column2[i] = T[i][1]; - column3[i] = T[i][2]; - // printf("column3(%i) = %1.3f\n", i, column3[i]); - } - - double *columns[3]; - columns[0] = &(column1[0]); - columns[1] = &(column2[0]); - columns[2] = &(column3[0]); - - - double COV[3][3]; - double covTemp = 0; - for(int i = 0; i < numGroups; i++){ - for(int j = i; j < numGroups; j++){ - - covTemp = cov(columns[i], columns[j], 3); - - COV[i][j] = covTemp; - COV[j][i] = covTemp; - - // printf("COV(%i , %i) = %1.3f\n", i, j, COV[i][j]); - } - } - - double sumdiagcov = 0; - for(int i = 0; i < numGroups; i++){ - sumdiagcov += COV[i][i]; - } - - free(yFilt); - free(yDown); - free(yCG); - - return sumdiagcov; - - -} diff --git a/src/C/SB_TransitionMatrix.h b/src/C/SB_TransitionMatrix.h deleted file mode 100644 index 06d71e4..0000000 --- a/src/C/SB_TransitionMatrix.h +++ /dev/null @@ -1,15 +0,0 @@ -// -// SB_TransitionMatrix.h -// -// -// Created by Carl Henning Lubba on 23/09/2018. -// - -#ifndef SB_TransitionMatrix_h -#define SB_TransitionMatrix_h - -#include - -extern double SB_TransitionMatrix_3ac_sumdiagcov(const double y[], const int size); - -#endif /* SB_TransitionMatrix_h */ diff --git a/src/C/SC_FluctAnal.c b/src/C/SC_FluctAnal.c deleted file mode 100644 index 938f8e3..0000000 --- a/src/C/SC_FluctAnal.c +++ /dev/null @@ -1,350 +0,0 @@ -#include -#include -#include -#include -#include -#include "stats.h" -#include "CO_AutoCorr.h" - -double SC_FluctAnal_2_50_1_logi_prop_r1(const double y[], const int size, const int lag, const char how[]) -{ - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - // generate log spaced tau vector - double linLow = log(5); - double linHigh = log(size/2); - - int nTauSteps = 50; - double tauStep = (linHigh - linLow) / (nTauSteps-1); - - int tau[50]; - for(int i = 0; i < nTauSteps; i++) - { - tau[i] = round(exp(linLow + i*tauStep)); - } - - // check for uniqueness, use ascending order - int nTau = nTauSteps; - for(int i = 0; i < nTauSteps-1; i++) - { - - while (tau[i] == tau[i+1] && i < nTau-1) - { - for(int j = i+1; j < nTauSteps-1; j++) - { - tau[j] = tau[j+1]; - } - // lost one - nTau -= 1; - } - } - - // fewer than 12 points -> leave. - if(nTau < 12){ - return 0; - } - - int sizeCS = size/lag; - double * yCS = malloc(sizeCS * sizeof(double)); - - /* - for(int i = 0; i < 50; i++) - { - printf("y[%i]=%1.3f\n", i, y[i]); - } - */ - - // transform input vector to cumsum - yCS[0] = y[0]; - for(int i = 0; i < sizeCS-1; i++) - { - yCS[i+1] = yCS[i] + y[(i+1)*lag]; - - /* - if(i<300) - printf("yCS[%i]=%1.3f\n", i, yCS[i]); - */ - } - - //for each value of tau, cut signal into snippets of length tau, detrend and - - // first generate a support for regression (detrending) - double * xReg = malloc(tau[nTau-1] * sizeof * xReg); - for(int i = 0; i < tau[nTau-1]; i++) - { - xReg[i] = i+1; - } - - // iterate over taus, cut signal, detrend and save amplitude of remaining signal - double * F = malloc(nTau * sizeof * F); - for(int i = 0; i < nTau; i++) - { - int nBuffer = sizeCS/tau[i]; - double * buffer = malloc(tau[i] * sizeof * buffer); - double m = 0.0, b = 0.0; - - //printf("tau[%i]=%i\n", i, tau[i]); - - F[i] = 0; - for(int j = 0; j < nBuffer; j++) - { - - //printf("%i th buffer\n", j); - - linreg(tau[i], xReg, yCS+j*tau[i], &m, &b); - - - for(int k = 0; k < tau[i]; k++) - { - buffer[k] = yCS[j*tau[i]+k] - (m * (k+1) + b); - //printf("buffer[%i]=%1.3f\n", k, buffer[k]); - } - - if (strcmp(how, "rsrangefit") == 0) { - F[i] += pow(max_(buffer, tau[i]) - min_(buffer, tau[i]), 2); - } - else if (strcmp(how, "dfa") == 0) { - for(int k = 0; k leave. - if(nTau < 8){ - return 0; - } - - // transform input vector to cumsum - for(int i = 0; i < size-1; i++) - { - y[i+1] = y[i] + y[i+1]; - } - - //for each value of tau, cut signal into snippets of length tau, detrend and - - // first generate a support for regression (detrending) - double * xReg = malloc(tau[nTau-1] * sizeof * xReg); - for(int i = 0; i < tau[nTau-1]; i++) - { - xReg[i] = i+1; - } - - // iterate over taus, cut signal, detrend and save amplitude of remaining signal - double * F = malloc(nTau * sizeof * F); - for(int i = 0; i < nTau; i++) - { - int nBuffer = size/tau[i]; - double * buffer = malloc(tau[i] * sizeof * buffer); - double m = 0.0, b = 0.0; - - F[i] = 0; - for(int j = 0; j < nBuffer; j++) - { - - linreg(tau[i], xReg, y+j*tau[i], &m, &b); - - for(int k = 0; k < tau[i]; k++) - { - buffer[k] = y[j*tau[i]+k] - (m * (k+1) + b); - } - - F[i] += pow(max(buffer, tau[i]) - min(buffer, tau[i]), 2); - } - - F[i] = sqrt(F[i]/nBuffer); - - free(buffer); - - } - - double * logtt = malloc(nTau * sizeof * logtt); - double * logFF = malloc(nTau * sizeof * logFF); - int ntt = nTau; - - for (int i = 0; i < nTau; i++) - { - logtt[i] = log(tau[i]); - logFF[i] = log(F[i]); - } - - int minPoints = 6; - int nsserr = (ntt - 2*minPoints + 1); - double * sserr = malloc(nsserr * sizeof * sserr); - double * buffer = malloc((ntt - minPoints + 1) * sizeof * buffer); - for (int i = minPoints; i < ntt - minPoints + 1; i++) - { - // this could be done with less variables of course - double m1 = 0.0, b1 = 0.0; - double m2 = 0.0, b2 = 0.0; - - sserr[i - minPoints] = 0.0; - - linreg(i, logtt, logFF, &m1, &b1); - linreg(ntt-i+1, logtt+i-1, logFF+i-1, &m2, &b2); - - for(int j = 0; j < i; j ++) - { - buffer[j] = logtt[j] * m1 + b1 - logFF[j]; - } - - sserr[i - minPoints] += norm(buffer, i); - - for(int j = 0; j < ntt-i+1; j++) - { - buffer[j] = logtt[j+i-1] * m2 + b2 - logFF[j+i-1]; - } - - sserr[i - minPoints] += norm(buffer, ntt-i+1); - - } - - double firstMinInd = 0.0; - double minimum = min(sserr, nsserr); - for(int i = 0; i < nsserr; i++) - { - if(sserr[i] == minimum) - { - firstMinInd = i + minPoints - 1; - break; - } - } - - free(xReg); - free(F); - free(logtt); - free(logFF); - free(sserr); - free(buffer); - - return (firstMinInd+1)/ntt; - -} - */ diff --git a/src/C/SC_FluctAnal.h b/src/C/SC_FluctAnal.h deleted file mode 100644 index ea2b648..0000000 --- a/src/C/SC_FluctAnal.h +++ /dev/null @@ -1,11 +0,0 @@ -#ifndef SC_FLUCTANAL -#define SC_FLUCTANAL -#include -#include -#include "stats.h" -#include "CO_AutoCorr.h" - -extern double SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1(const double y[], const int size); -extern double SC_FluctAnal_2_50_1_logi_prop_r1(const double y[], const int size, const char how[]); -extern double SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1(const double y[], const int size); -#endif diff --git a/src/C/SP_Summaries.c b/src/C/SP_Summaries.c deleted file mode 100644 index b511df4..0000000 --- a/src/C/SP_Summaries.c +++ /dev/null @@ -1,213 +0,0 @@ -// -// SP_Summaries.c -// -// -// Created by Carl Henning Lubba on 23/09/2018. -// - -#include "SP_Summaries.h" -#include "CO_AutoCorr.h" - -int welch(const double y[], const int size, const int NFFT, const double Fs, const double window[], const int windowWidth, double ** Pxx, double ** f){ - - double dt = 1.0/Fs; - double df = 1.0/(nextpow2(windowWidth))/dt; - double m = mean(y, size); - - // number of windows, should be 1 - int k = floor((double)size/((double)windowWidth/2.0))-1; - - // normalising scale factor - double KMU = k * pow(norm_(window, windowWidth),2); - - double * P = malloc(NFFT * sizeof(double)); - for(int i = 0; i < NFFT; i++){ - P[i] = 0; - } - - // fft variables - cplx * F = malloc(NFFT * sizeof *F); - cplx * tw = malloc(NFFT * sizeof *tw); - twiddles(tw, NFFT); - - double * xw = malloc(windowWidth * sizeof(double)); - for(int i = 0; i0 && i < Nout-1){ - (*Pxx)[i] *= 2; - } - } - /* - for(int i = 0; i < Nout; i++){ - printf("Pxx[%i]: %1.3f\n", i, Pxx[i]); - } - */ - - *f = malloc(Nout * sizeof(double)); - for(int i = 0; i < Nout; i++){ - (*f)[i] = (double)i*df; - } - /* - for(int i = 0; i < Nout; i++){ - printf("f[%i]: %1.3f\n", i, (*f)[i]); - } - */ - - free(P); - free(F); - free(tw); - free(xw); - - return Nout; -} - -double SP_Summaries_welch_rect(const double y[], const int size, const char what[]) -{ - - // NaN check - for(int i = 0; i < size; i++) - { - if(isnan(y[i])) - { - return NAN; - } - } - - // rectangular window for Welch-spectrum - double * window = malloc(size * sizeof(double)); - for(int i = 0; i < size; i++){ - window[i] = 1; - } - - double Fs = 1.0; // sampling frequency - int N = nextpow2(size); - - double * S; - double * f; - - // compute Welch-power - int nWelch = welch(y, size, N, Fs, window, size, &S, &f); - free(window); - - // angualr frequency and spectrum on that - double * w = malloc(nWelch * sizeof(double)); - double * Sw = malloc(nWelch * sizeof(double)); - - double PI = 3.14159265359; - for(int i = 0; i < nWelch; i++){ - w[i] = 2*PI*f[i]; - Sw[i] = S[i]/(2*PI); - //printf("w[%i]=%1.3f, Sw[%i]=%1.3f\n", i, w[i], i, Sw[i]); - if(isinf(Sw[i]) || isinf(-Sw[i])){ - return 0; - } - } - - double dw = w[1] - w[0]; - - double * csS = malloc(nWelch * sizeof(double)); - cumsum(Sw, nWelch, csS); - /* - for(int i=0; i csSThres){ - centroid = w[i]; - break; - } - } - - output = centroid; - - } - else if(strcmp(what, "area_5_1") == 0){ - double area_5_1 = 0;; - for(int i=0; i - -extern double SP_Summaries_welch_rect(const double y[], const int size, const char what[]); -extern double SP_Summaries_welch_rect_area_5_1(const double y[], const int size); -extern double SP_Summaries_welch_rect_centroid(const double y[], const int size); - -#endif /* SP_Summaries_h */ diff --git a/src/C/butterworth.c b/src/C/butterworth.c deleted file mode 100644 index 62cba09..0000000 --- a/src/C/butterworth.c +++ /dev/null @@ -1,298 +0,0 @@ -// -// butterworth.c -// -// -// Created by Carl Henning Lubba on 23/09/2018. -// - -#include -#include - -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) - typedef double complex cplx; -#elif defined(_MSC_VER) - typedef _Dcomplex cplx; -#endif -#endif - -#include "helper_functions.h" -#include "butterworth.h" - -#ifndef CMPLX -#define CMPLX(x, y) ((cplx)((double)(x) + _Imaginary_I * (double)(y))) -#endif - -void poly(cplx x[], int size, cplx out[]) -{ - /* Convert roots x to polynomial coefficients */ - - // initialise - #if defined(__GNUC__) || defined(__GNUG__) - out[0] = 1; - for(int i=1; i= 0) - { - out[i] += b[j]*(y[i-j]-offset); - out[i] -= a[j]*out[i-j]; - } - else{ - out[i] += 0; //b[j]*offset; // 'padding' - out[i] -= 0; //a[j]*offset; - } - } - } - - for(int i = 0; i < size; i++){ - out[i] += offset; - } -} - -void reverse_array(double a[], int size){ - - /* Reverse the order of the elements in an array. Write back into the input array.*/ - - double temp; - for(int i = 0; i < size/2; i++){ - temp = a[i]; - a[i] = a[size-i-1]; - a[size-1-i] = temp; - /* - printf("indFrom = %i, indTo = %i\n", i, size-1-i); - for(int i=0; i < size; i++){ - printf("reversed[%i]=%1.3f\n", i, a[i]); - } - */ - } -} - -void filt_reverse(double y[], int size, double a[], double b[], int nCoeffs, double out[]){ - - /* Filter a signal y with the filter coefficients a and b _in reverse order_, output to array out.*/ - - double * yTemp = malloc(size * sizeof(double)); - for(int i = 0; i < size; i++){ - yTemp[i] = y[i]; - } - - /* - for(int i=0; i < size; i++){ - printf("yTemp[%i]=%1.3f\n", i, yTemp[i]); - } - */ - - reverse_array(yTemp, size); - - /* - for(int i=0; i < size; i++){ - printf("reversed[%i]=%1.3f\n", i, yTemp[i]); - } - */ - - double offset = yTemp[0]; - - for(int i = 0; i < size; i++){ - out[i] = 0; - for(int j = 0; j < nCoeffs; j++){ - if(i - j >= 0) - { - out[i] += b[j]*(yTemp[i-j]-offset); - out[i] -= a[j]*out[i-j]; - } - else{ - out[i] += 0; //b[j]*offset; // 'padding' - out[i] -= 0; //a[j]*offset; - } - } - } - - for(int i = 0; i < size; i++){ - out[i] += offset; - } - - reverse_array(out, size); - - free(yTemp); - -} - -/* -void butterworthFilter(const double y[], int size, const int nPoles, const double W, double out[]){ - - double PI = 3.14159265359; - - double V = tan(W * PI/2); - cplx * Q = malloc(nPoles * sizeof(cplx)); - - for(int i = 0; i 0){ - // prod1mSp *= (1 - Sp[i]); - //} - - } - - double G = creal(Sg / prod1mSp); - - cplx * Zpoly = malloc((nPoles+1) * sizeof(cplx)); - cplx * Ppoly = malloc((nPoles+1) * sizeof(cplx)); - - // polynomial coefficients from poles and zeros for filtering - poly(Z, nPoles, Zpoly); - - //for(int i = 0; i < nPoles+1; i++){ - // printf("Zpoly[%i]= %1.3f + %1.3f i\n", i, creal(Zpoly[i]), cimag(Zpoly[i])); - //} - - poly(P, nPoles, Ppoly); - - //for(int i = 0; i < nPoles+1; i++){ - // printf("Ppoly[%i]= %1.3f + %1.3f i\n", i, creal(Ppoly[i]), cimag(Ppoly[i])); - //} - - - // coeffs for filtering - double * b = malloc((nPoles+1) * sizeof(double)); // zeros - double * a = malloc((nPoles+1) * sizeof(double)); // poles - - for(int i = 0; i - -extern void butterworthFilter(const double y[], const int size, const int nPoles, const double W, double out[]); - -#endif /* butterworth_h */ diff --git a/src/C/fft.c b/src/C/fft.c deleted file mode 100644 index 31f5037..0000000 --- a/src/C/fft.c +++ /dev/null @@ -1,61 +0,0 @@ -#include -#include -#include - -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -#ifndef CMPLX -#define CMPLX(x, y) ((cplx)((double)(x) + _Imaginary_I * (double)(y))) -#endif - -#include "helper_functions.h" - -void twiddles(cplx a[], int size) -{ - - double PI = 3.14159265359; - - for (int i = 0; i < size; i++) { - // cplx tmp = { 0, -PI * i / size }; - #if defined(__GNUC__) || defined(__GNUG__) - cplx tmp = 0.0 - PI * i / size * I; - #elif defined(_MSC_VER) - cplx tmp = {0.0, -PI * i / size }; - #endif - a[i] = cexp(tmp); - //a[i] = cexp(-I * M_PI * i / size); - } -} - -static void _fft(cplx a[], cplx out[], int size, int step, cplx tw[]) -{ - if (step < size) { - _fft(out, a, size, step * 2, tw); - _fft(out + step, a + step, size, step * 2, tw); - - for (int i = 0; i < size; i += 2 * step) { - //cplx t = tw[i] * out[i + step]; - cplx t = _Cmulcc(tw[i], out[i + step]); - a[i / 2] = _Caddcc(out[i], t); - a[(i + size) / 2] = _Cminuscc(out[i], t); - } - } -} - -void fft(cplx a[], int size, cplx tw[]) -{ - cplx * out = malloc(size * sizeof(cplx)); - memcpy(out, a, size * sizeof(cplx)); - _fft(a, out, size, 1, tw); - free(out); -} diff --git a/src/C/fft.h b/src/C/fft.h deleted file mode 100644 index b26a0ae..0000000 --- a/src/C/fft.h +++ /dev/null @@ -1,26 +0,0 @@ -#ifndef FFT_H -#define FFT_H -//#include - -#if __cplusplus -# include - typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -#include -#include -#ifndef CMPLX -#define CMPLX(x, y) ((cplx)((double)(x) + _Complex_I * (double)(y))) -#endif -extern void twiddles(cplx a[], int size); -// extern void _fft(cplx a[], cplx out[], int size, int step, cplx tw[]); -extern void fft(cplx a[], int size, cplx tw[]); -extern void ifft(cplx a[], int size, cplx tw[]); -#endif diff --git a/src/C/helper_functions.c b/src/C/helper_functions.c deleted file mode 100644 index fb9c08d..0000000 --- a/src/C/helper_functions.c +++ /dev/null @@ -1,179 +0,0 @@ -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -#include -#include -#include -#include -#include "stats.h" - -// compare function for qsort, for array of doubles -static int compare (const void * a, const void * b) -{ - if (*(double*)a < *(double*)b) { - return -1; - } else if (*(double*)a > *(double*)b) { - return 1; - } else { - return 0; - } -} - -// wrapper for qsort for array of doubles. Sorts in-place -void sort(double y[], int size) -{ - qsort(y, size, sizeof(*y), compare); -} - -// linearly spaced vector -void linspace(double start, double end, int num_groups, double out[]) -{ - double step_size = (end - start) / (num_groups - 1); - for (int i = 0; i < num_groups; i++) { - out[i] = start; - start += step_size; - } - return; -} - -double quantile(const double y[], const int size, const double quant) -{ - double quant_idx, q, value; - int idx_left, idx_right; - double * tmp = malloc(size * sizeof(*y)); - memcpy(tmp, y, size * sizeof(*y)); - sort(tmp, size); - - /* - for(int i=0; i < size; i++){ - printf("y[%i]=%1.4f\n", i, y[i]); - } - for(int i=0; i < size; i++){ - printf("sorted[%i]=%1.4f\n", i, tmp[i]); - } - */ - - // out of range limit? - q = 0.5 / size; - if (quant < q) { - value = tmp[0]; // min value - free(tmp); - return value; - } else if (quant > (1 - q)) { - value = tmp[size - 1]; // max value - free(tmp); - return value; - } - - quant_idx = size * quant - 0.5; - idx_left = (int)floor(quant_idx); - idx_right = (int)ceil(quant_idx); - value = tmp[idx_left] + (quant_idx - idx_left) * (tmp[idx_right] - tmp[idx_left]) / (idx_right - idx_left); - free(tmp); - return value; -} - -void binarize(const double a[], const int size, int b[], const char how[]) -{ - double m = 0.0; - if (strcmp(how, "mean") == 0) { - m = mean(a, size); - } else if (strcmp(how, "median") == 0) { - m = median(a, size); - } - for (int i = 0; i < size; i++) { - b[i] = (a[i] > m) ? 1 : 0; - } - return; -} - -double f_entropy(const double a[], const int size) -{ - double f = 0.0; - for (int i = 0; i < size; i++) { - if (a[i] > 0) { - f += a[i] * log(a[i]); - } - } - return -1 * f; -} - -void subset(const int a[], int b[], const int start, const int end) -{ - int j = 0; - for (int i = start; i < end; i++) { - b[j++] = a[i]; - } - return; -} - -#if defined(__GNUC__) || defined(__GNUG__) - cplx _Cmulcc(const cplx x, const cplx y) { - /*double a = x._Val[0]; - double b = x._Val[1]; - - double c = y._Val[0]; - double d = y._Val[1]; - - cplx result = { (a * c - b * d), (a * d + c * b) }; - */ - return x*y; - } - - cplx _Cminuscc(const cplx x, const cplx y) { - //cplx result = { x._Val[0] - y._Val[0], x._Val[1] - y._Val[1] }; - return x - y; - } - - cplx _Caddcc(const cplx x, const cplx y) { - // cplx result = { x._Val[0] + y._Val[0], x._Val[1] + y._Val[1] }; - return x + y; - } - - cplx _Cdivcc(const cplx x, const cplx y) { - - double a = creal(x); - double b = cimag(x); - - double c = creal(y); - double d = cimag(y); - - cplx result = (a*c + b*d) / (c*c + d*d) + (b*c - a*d)/(c*c + d*d) * I; - - return result; - - // return x / y; - } - -#elif defined(_MSC_VER) - cplx _Cminuscc(const cplx x, const cplx y) { - cplx result = { x._Val[0] - y._Val[0], x._Val[1] - y._Val[1] }; - return result; - } - - cplx _Caddcc(const cplx x, const cplx y) { - cplx result = { x._Val[0] + y._Val[0], x._Val[1] + y._Val[1] }; - return result; - } - - cplx _Cdivcc(const cplx x, const cplx y) { - double a = x._Val[0]; - double b = x._Val[1]; - - double c = y._Val[0]; - double d = y._Val[1]; - - cplx result = { (a*c + b*d) / (c*c + d*d), (b*c - a*d)/(c*c + d*d)}; - - return result; - } -#endif diff --git a/src/C/helper_functions.h b/src/C/helper_functions.h deleted file mode 100644 index edf0c66..0000000 --- a/src/C/helper_functions.h +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef HELPER_FUNCTIONS_H -#define HELPER_FUNCTIONS_H -#include -#include -#include -#include "stats.h" - -#if __cplusplus -# include -typedef std::complex< double > cplx; -#else -# include -#if defined(__GNUC__) || defined(__GNUG__) -typedef double complex cplx; -#elif defined(_MSC_VER) -typedef _Dcomplex cplx; -#endif -#endif - -extern void linspace(double start, double end, int num_groups, double out[]); -extern double quantile(const double y[], const int size, const double quant); -extern void sort(double y[], int size); -extern void binarize(const double a[], const int size, int b[], const char how[]); -extern double f_entropy(const double a[], const int size); -extern void subset(const int a[], int b[], const int start, const int end); - -extern cplx _Cminuscc(const cplx x, const cplx y); -extern cplx _Caddcc(const cplx x, const cplx y); -extern cplx _Cdivcc(const cplx x, const cplx y); -#if defined(__GNUC__) || defined(__GNUG__) -extern cplx _Cmulcc(const cplx x, const cplx y); -#endif - -#endif diff --git a/src/C/histcounts.c b/src/C/histcounts.c deleted file mode 100644 index 47c4388..0000000 --- a/src/C/histcounts.c +++ /dev/null @@ -1,200 +0,0 @@ -// -// histcounts.c -// C_polished -// -// Created by Carl Henning Lubba on 19/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#include -#include - -#include "stats.h" -#include "histcounts.h" - -int num_bins_auto(const double y[], const int size){ - - double maxVal = max_(y, size); - double minVal = min_(y, size); - - if (stddev(y, size) < 0.001){ - return 0; - } - - return ceil((maxVal-minVal)/(3.5*stddev(y, size)/pow(size, 1/3.))); - -} - -int histcounts_preallocated(const double y[], const int size, int nBins, int * binCounts, double * binEdges) -{ - - int i = 0; - - // check min and max of input array - double minVal = DBL_MAX, maxVal=-DBL_MAX; - for(int i = 0; i < size; i++) - { - // printf("histcountInput %i: %1.3f\n", i, y[i]); - - if (y[i] < minVal) - { - minVal = y[i]; - } - if (y[i] > maxVal) - { - maxVal = y[i]; - } - } - - // and derive bin width from it - double binStep = (maxVal - minVal)/nBins; - - // variable to store counted occurances in - for(i = 0; i < nBins; i++) - { - binCounts[i] = 0; - } - - for(i = 0; i < size; i++) - { - - int binInd = (y[i]-minVal)/binStep; - if(binInd < 0) - binInd = 0; - if(binInd >= nBins) - binInd = nBins-1; - //printf("histcounts, i=%i, binInd=%i, nBins=%i\n", i, binInd, nBins); - binCounts[binInd] += 1; - - } - - for(i = 0; i < nBins+1; i++) - { - binEdges[i] = i * binStep + minVal; - } - - /* - // debug - for(i=0;i maxVal) - { - maxVal = y[i]; - } - } - - // if no number of bins given, choose spaces automatically - if (nBins <= 0){ - nBins = ceil((maxVal-minVal)/(3.5*stddev(y, size)/pow(size, 1/3.))); - } - - // and derive bin width from it - double binStep = (maxVal - minVal)/nBins; - - // variable to store counted occurances in - *binCounts = malloc(nBins * sizeof(int)); - for(i = 0; i < nBins; i++) - { - (*binCounts)[i] = 0; - } - - for(i = 0; i < size; i++) - { - - int binInd = (y[i]-minVal)/binStep; - if(binInd < 0) - binInd = 0; - if(binInd >= nBins) - binInd = nBins-1; - (*binCounts)[binInd] += 1; - - } - - *binEdges = malloc((nBins+1) * sizeof(double)); - for(i = 0; i < nBins+1; i++) - { - (*binEdges)[i] = i * binStep + minVal; - } - - /* - // debug - for(i=0;i 0 - binIdentity[i] = 0; - - // go through bin edges - for(int j = 0; j < nEdges; j++){ - if(y[i] < binEdges[j]){ - binIdentity[i] = j; - break; - } - } - } - - return binIdentity; - -} - -int * histcount_edges(const double y[], const int size, const double binEdges[], const int nEdges) -{ - - - int * histcounts = malloc(nEdges * sizeof(int)); - for(int i = 0; i < nEdges; i++){ - histcounts[i] = 0; - } - - for(int i = 0; i < size; i++) - { - // go through bin edges - for(int j = 0; j < nEdges; j++){ - if(y[i] <= binEdges[j]){ - histcounts[j] += 1; - break; - } - } - } - - return histcounts; - -} diff --git a/src/C/histcounts.h b/src/C/histcounts.h deleted file mode 100644 index 923f882..0000000 --- a/src/C/histcounts.h +++ /dev/null @@ -1,22 +0,0 @@ -// -// histcounts.h -// C_polished -// -// Created by Carl Henning Lubba on 19/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef histcounts_h -#define histcounts_h - -#include -#include -#include - -extern int num_bins_auto(const double y[], const int size); -extern int histcounts(const double y[], const int size, int nBins, int ** binCounts, double ** binEdges); -extern int histcounts_preallocated(const double y[], const int size, int nBins, int * binCounts, double * binEdges); -extern int * histcount_edges(const double y[], const int size, const double binEdges[], const int nEdges); -extern int * histbinassign(const double y[], const int size, const double binEdges[], const int nEdges); - -#endif /* histcounts_h */ diff --git a/src/C/main.c b/src/C/main.c deleted file mode 100644 index 7d44e57..0000000 --- a/src/C/main.c +++ /dev/null @@ -1,420 +0,0 @@ -/* Include files */ -#include "main.h" -#include -#include -#include -#include -//#include - -#include "DN_HistogramMode_5.h" -#include "DN_HistogramMode_10.h" -#include "DN_Mean.h" -#include "DN_Spread_Std.h" -#include "CO_AutoCorr.h" -#include "DN_OutlierInclude.h" -#include "FC_LocalSimple.h" -#include "IN_AutoMutualInfoStats.h" -#include "MD_hrv.h" -#include "SB_BinaryStats.h" -#include "SB_MotifThree.h" -#include "SC_FluctAnal.h" -#include "SP_Summaries.h" -#include "SB_TransitionMatrix.h" -#include "PD_PeriodicityWang.h" - -#include "stats.h" - -// check if data qualifies to be caught22 -int quality_check(const double y[], const int size) -{ - int minSize = 10; - - if(size < minSize) - { - return 1; - } - for(int i = 0; i < size; i++) - { - double val = y[i]; - if(val == INFINITY || -val == INFINITY) - { - return 2; - } - if(isnan(val)) - { - return 3; - } - } - return 0; -} - -void run_features(double y[], int size, FILE * outfile, bool catch24) -{ - int quality = quality_check(y, size); - if(quality != 0) - { - fprintf(stdout, "Time series quality test not passed (code %i).\n", quality); - return; - } - - double * y_zscored = malloc(size * sizeof * y_zscored); - - // variables to keep time - clock_t begin; - double timeTaken; - - // output - double result; - - // z-score first for all. - zscore_norm2(y, size, y_zscored); - - // GOOD - begin = clock(); - result = DN_OutlierInclude_n_001_mdrmd(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_OutlierInclude_n_001_mdrmd", timeTaken); - - // GOOD - begin = clock(); - result = DN_OutlierInclude_p_001_mdrmd(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_OutlierInclude_p_001_mdrmd", timeTaken); - - // GOOD - begin = clock(); - result = DN_HistogramMode_5(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_HistogramMode_5", timeTaken); - - // GOOD - begin = clock(); - result = DN_HistogramMode_10(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_HistogramMode_10", timeTaken); - - //GOOD - begin = clock(); - result = CO_Embed2_Dist_tau_d_expfit_meandiff(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "CO_Embed2_Dist_tau_d_expfit_meandiff", timeTaken); - - //GOOD (memory leak?) - begin = clock(); - result = CO_f1ecac(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "CO_f1ecac", timeTaken); - - //GOOD - begin = clock(); - result = CO_FirstMin_ac(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "CO_FirstMin_ac", timeTaken); - - // GOOD (memory leak?) - begin = clock(); - result = CO_HistogramAMI_even_2_5(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "CO_HistogramAMI_even_2_5", timeTaken); - - // GOOD - begin = clock(); - result = CO_trev_1_num(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "CO_trev_1_num", timeTaken); - - //GOOD - begin = clock(); - result = FC_LocalSimple_mean1_tauresrat(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "FC_LocalSimple_mean1_tauresrat", timeTaken); - - //GOOD - begin = clock(); - result = FC_LocalSimple_mean3_stderr(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "FC_LocalSimple_mean3_stderr", timeTaken); - - //GOOD (memory leak?) - begin = clock(); - result = IN_AutoMutualInfoStats_40_gaussian_fmmi(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "IN_AutoMutualInfoStats_40_gaussian_fmmi", timeTaken); - - //GOOD - begin = clock(); - result = MD_hrv_classic_pnn40(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "MD_hrv_classic_pnn40", timeTaken); - - //GOOD - begin = clock(); - result = SB_BinaryStats_diff_longstretch0(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SB_BinaryStats_diff_longstretch0", timeTaken); - - //GOOD - begin = clock(); - result = SB_BinaryStats_mean_longstretch1(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SB_BinaryStats_mean_longstretch1", timeTaken); - - //GOOD (memory leak?) - begin = clock(); - result = SB_MotifThree_quantile_hh(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SB_MotifThree_quantile_hh", timeTaken); - - //GOOD (memory leak?) - begin = clock(); - result = SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1", timeTaken); - - //GOOD - begin = clock(); - result = SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1", timeTaken); - - //GOOD - begin = clock(); - result = SP_Summaries_welch_rect_area_5_1(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SP_Summaries_welch_rect_area_5_1", timeTaken); - - //GOOD - begin = clock(); - result = SP_Summaries_welch_rect_centroid(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SP_Summaries_welch_rect_centroid", timeTaken); - - //OK, BUT filt in Butterworth sometimes diverges, now removed alltogether, let's see results. - begin = clock(); - result = SB_TransitionMatrix_3ac_sumdiagcov(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "SB_TransitionMatrix_3ac_sumdiagcov", timeTaken); - - // GOOD - begin = clock(); - result = PD_PeriodicityWang_th0_01(y_zscored, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "PD_PeriodicityWang_th0_01", timeTaken); - - if (catch24) { - - // GOOD - begin = clock(); - result = DN_Mean(y, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_Mean", timeTaken); - - // GOOD - begin = clock(); - result = DN_Spread_Std(y, size); - timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC; - fprintf(outfile, "%.14f, %s, %f\n", result, "DN_Spread_Std", timeTaken); - } else { - - } - - fprintf(outfile, "\n"); - - free(y_zscored); -} - -void print_help(char *argv[], char msg[]) -{ - if (strlen(msg) > 0) { - fprintf(stdout, "ERROR: %s\n", msg); - } - fprintf(stdout, "Usage is %s \n", argv[0]); - fprintf(stdout, "\n\tSpecifying outfile is optional, by default it is stdout\n"); - // fprintf(stdout, "\tOutput order is:\n%s\n", HEADER); - exit(1); -} - -// memory leak check; use with valgrind. -#if 0 -int main(int argc, char * argv[]) -{ - double * y = malloc(1000 * sizeof(double)); - - srand(42); - for (int i = 0; i < 1000; ++i) { - y[i] = rand() % RAND_MAX; - } - run_features(y, 1000, stdout); - free(y); -} -#endif - -#if 1 -int main(int argc, char * argv[]) -{ - FILE * infile, * outfile; - int array_size; - double * y; - int size; - double value; - // DIR *d; - struct dirent *dir; - - - switch (argc) { - case 1: - print_help(argv, ""); - break; - case 2: - if ((infile = fopen(argv[1], "r")) == NULL) { - print_help(argv, "Can't open input file\n"); - } - outfile = stdout; - break; - case 3: - if ((infile = fopen(argv[1], "r")) == NULL) { - print_help(argv, "Can't open input file\n"); - } - if ((outfile = fopen(argv[2], "w")) == NULL) { - print_help(argv, "Can't open output file\n"); - } - break; - } - - /* - // debug: fix these. - infile = fopen("/Users/carl/PycharmProjects/catch22/C/timeSeries/tsid0244.txt", "r"); - outfile = stdout; - */ - - // fprintf(outfile, "%s", HEADER); - array_size = 50; - size = 0; - y = malloc(array_size * sizeof *y); - - while (fscanf(infile, "%lf", &value) != EOF) { - if (size == array_size) { - y = realloc(y, 2 * array_size * sizeof *y); - array_size *= 2; - } - y[size++] = value; - } - fclose(infile); - y = realloc(y, size * sizeof *y); - //printf("size=%i\n", size); - - // catch24 specification - - int catch24; - printf("Do you want to run catch24? Enter 0 for catch22 or 1 for catch24."); - scanf("%d", &catch24); - - if (catch24 == 1) { - run_features(y, size, outfile, true); - } else { - run_features(y, size, outfile, false); - } - - fclose(outfile); - free(y); - - return 0; -} -#endif - -#if 0 -int main(int argc, char * argv[]) -{ - (void)argc; - (void)argv; - - /* - // generate some data - const int size = 31; // 211; - - double y[size]; - int i; - double sinIn=0; - for(i=0; i -#include -#include -#include - -/* Function Declarations */ -//extern int main(int argc, const char * const argv[]); -extern int main(int argc, char * argv[]); - -#endif - -/* End of code generation (main.h) */ diff --git a/src/C/runAllTS.sh b/src/C/runAllTS.sh deleted file mode 100755 index 2ba8ed4..0000000 --- a/src/C/runAllTS.sh +++ /dev/null @@ -1,71 +0,0 @@ -#!/bin/bash - -help() -{ - echo "" - echo "Usage: $0 -i indir -o outdir -a append_string -s" - echo -e "\t-h Show this help message" - echo -e "\t-i Path to a directory containing input time-series files (.txt with one time series value per line). Default: './timeSeries'" - echo -e "\t-o Path to a directory in which to save output feature values. Default: './featureOutput'" - echo -e "\t-a A string (minus extension) appended to the input file names to create the output file names. Default: 'output'" - echo -e "\t-s A switch to evaluate catch22 (0) or catch24 (1). Default: 0" - exit 1 -} - -while getopts "i:o:a:s:h" opt -do - case "$opt" in - i) indir="$OPTARG" ;; - o) outdir="$OPTARG" ;; - a) append="$OPTARG" ;; - s) catch24="$OPTARG" ;; - h) help ;; - esac -done - -srcdir=$(dirname "$0}") - -if [ -z "$indir" ] -then - indir="./timeSeries" -fi - -if [ -z "$outdir" ] -then - outdir="./featureOutput" -fi - -if [ -z "$append" ] -then - append="output" -fi - -if [ -z "$catch24" ] -then - catch24=0 -fi - -indir="$(dirname $indir)/$(basename $indir)" -outdir="$(dirname $outdir)/$(basename $outdir)" -mkdir -p $outdir - -# Loop through each file in indir and save the feature outputs -for entry in "${indir}"/*.txt -do - filename=$(basename "$entry") - extension="${filename##*.}" - filename="${filename%.*}" - fullfile="${outdir}/${filename}${append}.${extension}" - if [ "${filename: -${#append}}" != "${append}" ] - then - yes $catch24 | "${srcdir}/run_features" $entry $fullfile > /dev/null - - if [ -s $fullfile ] # Remove file if catch22 errors - then - echo "Output written to ${fullfile}" - else - rm $fullfile - fi - - fi -done diff --git a/src/C/splinefit.c b/src/C/splinefit.c deleted file mode 100644 index 802298c..0000000 --- a/src/C/splinefit.c +++ /dev/null @@ -1,801 +0,0 @@ -// Created by Carl Henning Lubba on 27/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// -// Based on the work of Jonas Lundgren in his Matlab Central contribution 'SPLINEFIT'. -// -#include -#include -#include - -#include "splinefit.h" -#include "stats.h" - -#define nCoeffs 3 -#define nPoints 4 - -#define pieces 2 -#define nBreaks 3 -#define deg 3 -#define nSpline 4 -#define piecesExt 8 //3 * deg - 1 - - -void matrix_multiply(const int sizeA1, const int sizeA2, const double *A, const int sizeB1, const int sizeB2, const double *B, double *C){ -//void matrix_multiply(int sizeA1, int sizeA2, double **A, int sizeB1, int sizeB2, double **B, double C[sizeA1][sizeB2]){ - - if(sizeA2 != sizeB1){ - return; - } - - /* - // show input - for(int i = 0; i < sizeA1; i++){ - for(int j = 0; j < sizeA2; j++){ - printf("A[%i][%i] = %1.3f\n", i, j, A[i*sizeA2 + j]); - } - } - */ - - for(int i = 0; i < sizeA1; i++){ - for(int j = 0; j < sizeB2; j++){ - - //C[i][j] = 0; - C[i*sizeB2 + j] = 0; - for(int k = 0; k < sizeB1; k++){ - // C[i][j] += A[i][k]*B[k][j]; - C[i*sizeB2 + j] += A[i * sizeA2 + k]*B[k * sizeB2 + j]; - //printf("C[%i][%i] (k=%i) = %1.3f\n", i, j, k, C[i * sizeB2 + j]); - } - - } - } - -} - -void matrix_times_vector(const int sizeA1, const int sizeA2, const double *A, const int sizeb, const double *b, double *c){ //c[sizeb] - - if(sizeA2 != sizeb){ - return; - } - - // row - for(int i = 0; i < sizeA1; i++){ - - // column - c[i] = 0; - for(int k = 0; k < sizeb; k++){ - c[i] += A[i * sizeA2 + k]*b[k]; - } - - } - -} - -void gauss_elimination(int size, double *A, double *b, double *x){ -// void gauss_elimination(int size, double A[size][size], double b[size], double x[size]){ - - double factor; - - // create temp matrix and vector - // double *AElim[size]; - double* AElim[nSpline + 1]; - for (int i = 0; i < size; i++) - AElim[i] = (double *)malloc(size * sizeof(double)); - double * bElim = malloc(size * sizeof(double)); - - // -- create triangular matrix - - // initialise to A and b - for(int i = 0; i < size; i++){ - for(int j = 0; j < size; j++){ - AElim[i][j] = A[i*size + j]; - } - bElim[i] = b[i]; - } - - /* - printf("AElim\n"); - for(int i = 0; i < size; i++){ - for(int j = 0; j < size; j++){ - printf("%1.3f, ", AElim[i][j]); - } - printf("\n"); - } - */ - - // go through columns in outer loop - for(int i = 0; i < size; i++){ - - // go through rows to eliminate - for(int j = i+1; j < size; j++){ - - factor = AElim[j][i]/AElim[i][i]; - - // subtract in vector - bElim[j] = bElim[j] - factor*bElim[i]; - - // go through entries of this row - for(int k = i; k < size; k++){ - AElim[j][k] = AElim[j][k] - factor*AElim[i][k]; - } - - /* - printf("AElim i=%i, j=%i\n", i, j); - for(int i = 0; i < size; i++){ - for(int j = 0; j < size; j++){ - printf("%1.3f, ", AElim[i][j]); - } - printf("\n"); - } - */ - - } - - } - - /* - for(int i = 0; i < size; i++){ - for(int j = 0; j < size; j++){ - printf("AElim[%i][%i] = %1.3f\n", i, j, AElim[i][j]); - } - } - for(int i = 0; i < size; i++){ - printf("bElim[%i] = %1.3f\n", i, bElim[i]); - } - */ - - - // -- go backwards through triangular matrix and solve for x - - // row - double bMinusATemp; - for(int i = size-1; i >= 0; i--){ - - bMinusATemp = bElim[i]; - for(int j = i+1; j < size; j++){ - bMinusATemp -= x[j]*AElim[i][j]; - } - - x[i] = bMinusATemp/AElim[i][i]; - } - /* - for(int j = 0; j < size; j++){ - printf("x[%i] = %1.3f\n", j, x[j]); - } - */ - - for (int i = 0; i < size; i++) - free(AElim[i]); - free(bElim); -} - -void lsqsolve_sub(const int sizeA1, const int sizeA2, const double *A, const int sizeb, const double *b, double *x) -//void lsqsolve_sub(int sizeA1, int sizeA2, double A[sizeA1][sizeA2], int sizeb, double b[sizeb], double x[sizeA1]) -{ - // create temp matrix and vector - /* - double *AT[sizeA1*sizeA2]; - for (int i = 0; i < sizeA2; i++) - AT[i] = (double *)malloc(sizeA1 * sizeof(double)); - double *ATA[sizeA2]; - for (int i = 0; i < sizeA2; i++) - ATA[i] = (double *)malloc(sizeA2 * sizeof(double)); - double * ATb = malloc(sizeA1 * sizeof(double)); - */ - - double * AT = malloc(sizeA2 * sizeA1 * sizeof(double)); - double * ATA = malloc(sizeA2 * sizeA2 * sizeof(double)); - double * ATb = malloc(sizeA2 * sizeof(double)); - - - for(int i = 0; i < sizeA1; i++){ - for(int j = 0; j < sizeA2; j++){ - //AT[i,j] = A[j,i] - AT[j * sizeA1 + i] = A[i * sizeA2 + j]; - } - } - - /* - printf("\n b \n"); - for(int i = 0; i < sizeA1; i++){ - printf("%i, %1.3f\n", i, b[i]); - } - */ - - /* - printf("\nA\n"); - for(int i = 0; i < sizeA2; i++){ - for(int j = 0; j < sizeA1; j++){ - printf("%1.3f, ", AT[i * sizeA1 + j]); - } - printf("\n"); - } - */ - - - matrix_multiply(sizeA2, sizeA1, AT, sizeA1, sizeA2, A, ATA); - - /* - printf("ATA\n"); - for(int i = 0; i < sizeA2; i++){ - for(int j = 0; j < sizeA2; j++){ - printf("%1.3f, ", ATA[i * sizeA2 + j]); - } - printf("\n"); - } - */ - - - - matrix_times_vector(sizeA2, sizeA1, AT, sizeA1, b, ATb); - - /* - for(int i = 0; i < sizeA2; i++){ - ATb[i] = 0; - for(int j = 0; j < sizeA1; j++){ - ATb[i] += AT[i*sizeA1 + j]*b[j]; - //printf("%i, ATb[%i]=%1.3f, AT[i*sizeA1 + j]=%1.3f, b[j]=%1.3f\n", i, i, ATb[i], AT[i*sizeA1 + j],b[j]); - } - } - */ - - /* - for(int i = 0; i < nCoeffs; i++){ - printf("b[%i] = %1.3f\n", i, b[i]); - } - */ - - /* - for(int i = 0; i < sizeA2; i++){ - printf("ATb[%i] = %1.3f\n", i, ATb[i]); - } - */ - - - gauss_elimination(sizeA2, ATA, ATb, x); - - free(AT); - free(ATA); - free(ATb); - -} - -/* -int lsqsolve() -{ - //const int nPoints = 4; - //const int nCoeffs = 3; - - //double A[nPoints][nCoeffs] = {}; - double A[4][3]; - A[0][0] = 1; - A[1][0] = 3; - A[2][0] = 6; - A[3][0] = 8; - A[0][1] = 4; - A[1][1] = 5; - A[2][1] = 3; - A[3][1] = 12; - A[0][2] = 4; - A[1][2] = 1; - A[2][2] = 0; - A[3][2] = 7; - //double b[nPoints] = {}; - double b[4]; - b[0] = 2; - b[1] = 8; - b[2] = 3; - b[3] = 1; - - double x[4]; - - double * Alin = malloc(nPoints * nCoeffs * sizeof(double)); - - for(int i = 0; i < nPoints; i++){ - for(int j = 0; j < nCoeffs; j++){ - //AT[i,j] = A[j,i] - Alin[i * nCoeffs + j] = A[i][j]; - } - } - - lsqsolve_sub(nPoints, nCoeffs, Alin, nPoints, b, x); - - free(Alin); - - return 0; - - -} -*/ - -int iLimit(int x, int lim){ - return x= breaks[breakInd] && breakInd= breaks[1]) - breakInd = 1; - A[(i%nSpline)+breakInd + (i/nSpline)*(nSpline+1)] = vB[i]; - } - - /* - printf("\nA:\n"); - for(int i = 0; i < size; i++){ - for(int j = 0; j < n+1; j++){ - printf("%1.5f, ", A[i * (n+1) + j]); - } - printf("\n"); - } - */ - - - - double * x = malloc((nSpline+1)*sizeof(double)); - // lsqsolve_sub(int sizeA1, int sizeA2, double *A, int sizeb, double *b, double *x) - lsqsolve_sub(size, nSpline+1, A, size, y, x); - - /* - printf("\nsolved x\n"); - for(int i = 0; i < n+1; i++){ - printf("%i, %1.4f\n", i, x[i]); - } - */ - - // coeffs of B-splines to combine by optimised weighting in x - double C[pieces+nSpline-1][nSpline*pieces]; - // initialise to 0 - for(int i = 0; i < nSpline+1; i++){ - for(int j = 0; j < nSpline*pieces; j++){ - C[i][j] = 0; - } - } - - int CRow, CCol, coefRow, coefCol; - for(int i = 0; i < nSpline*nSpline*pieces; i++){ - - CRow = i%nSpline + (i/nSpline)%2; - CCol = i/nSpline; - - coefRow = i%(nSpline*2); - coefCol =i/(nSpline*2); - - C[CRow][CCol] = coefsOut[coefRow][coefCol]; - - } - - /* - printf("\nC:\n"); - for(int i = 0; i < n+1; i++){ - for(int j = 0; j < n*pieces; j++){ - printf("%1.5f, ", C[i][j]); - } - printf("\n"); - } - */ - - // final coefficients - double coefsSpline[pieces][nSpline]; - for(int i = 0; i < pieces; i++){ - for(int j = 0; j < nSpline; j++){ - coefsSpline[i][j] = 0; - } - } - - //multiply with x - for(int j = 0; j < nSpline*pieces; j++){ - coefCol = j/pieces; - coefRow = j%pieces; - - for(int i = 0; i < nSpline+1; i++){ - - coefsSpline[coefRow][coefCol] += C[i][j]*x[i]; - - } - } - - /* - printf("\ncoefsSpline:\n"); - for(int i = 0; i < pieces; i++){ - for(int j = 0; j < n; j++){ - printf("%1.5f, ", coefsSpline[i][j]); - } - printf("\n"); - } - */ - - - // compute piecewise polynomial - - int secondHalf = 0; - for(int i = 0; i < size; i++){ - secondHalf = i < breaks[1] ? 0 : 1; - yOut[i] = coefsSpline[secondHalf][0]; - } - - /* - printf("\nvSpline first iter\n"); - for(int i = 0; i < size; i++){ - printf("%i, %1.5f\n", i, vSpline[i]); - } - */ - - for(int i = 1; i < nSpline; i ++){ - for(int j = 0; j < size; j++){ - secondHalf = j < breaks[1] ? 0 : 1; - yOut[j] = yOut[j]*(j - breaks[1]*secondHalf) + coefsSpline[secondHalf][i]; - } - - /* - printf("\nvSpline %i th iter\n", i); - for(int i = 0; i < size; i++){ - printf("%i, %1.4f\n", i, vSpline[i]); - } - */ - } - - /* - printf("\nvSpline\n"); - for(int i = 0; i < size; i++){ - printf("%i, %1.4f\n", i, yOut[i]); - } - */ - - free(xsB); - free(indexB); - free(vB); - free(A); - free(x); - - return 0; - -} - - diff --git a/src/C/splinefit.h b/src/C/splinefit.h deleted file mode 100644 index 25c74db..0000000 --- a/src/C/splinefit.h +++ /dev/null @@ -1,16 +0,0 @@ -// -// splinefit.h -// C_polished -// -// Created by Carl Henning Lubba on 27/09/2018. -// Copyright © 2018 Carl Henning Lubba. All rights reserved. -// - -#ifndef splinefit_h -#define splinefit_h - -#include - -extern int splinefit(const double *y, const int size, double *yOut); - -#endif /* splinefit_h */ diff --git a/src/C/stats.c b/src/C/stats.c deleted file mode 100644 index 5573ed4..0000000 --- a/src/C/stats.c +++ /dev/null @@ -1,270 +0,0 @@ -#include -#include -#include -#include -#include "helper_functions.h" - -double min_(const double a[], const int size) -{ - double m = a[0]; - for (int i = 1; i < size; i++) { - if (a[i] < m) { - m = a[i]; - } - } - return m; -} - -double max_(const double a[], const int size) -{ - double m = a[0]; - for (int i = 1; i < size; i++) { - if (a[i] > m) { - m = a[i]; - } - } - return m; -} - -double mean(const double a[], const int size) -{ - double m = 0.0; - for (int i = 0; i < size; i++) { - m += a[i]; - } - m /= size; - return m; -} - -double sum(const double a[], const int size) -{ - double m = 0.0; - for (int i = 0; i < size; i++) { - m += a[i]; - } - return m; -} - -void cumsum(const double a[], const int size, double b[]) -{ - b[0] = a[0]; - for (int i = 1; i < size; i++) { - b[i] = a[i] + b[i-1]; - //printf("b[%i]%1.3f = a[%i]%1.3f + b[%i-1]%1.3f\n", i, b[i], i, a[i], i, a[i-1]); - } - -} - -void icumsum(const int a[], const int size, int b[]) -{ - b[0] = a[0]; - for (int i = 1; i < size; i++) { - b[i] = a[i] + b[i-1]; - //printf("b[%i]%1.3f = a[%i]%1.3f + b[%i-1]%1.3f\n", i, b[i], i, a[i], i, a[i-1]); - } - -} - -double isum(const int a[], const int size) -{ - double m = 0.0; - for (int i = 0; i < size; i++) { - m += a[i]; - } - return m; -} - -double median(const double a[], const int size) -{ - double m; - double * b = malloc(size * sizeof *b); - memcpy(b, a, size * sizeof *b); - sort(b, size); - if (size % 2 == 1) { - m = b[size / 2]; - } else { - int m1 = size / 2; - int m2 = m1 - 1; - m = (b[m1] + b[m2]) / (double)2.0; - } - free(b); - return m; -} - -double stddev(const double a[], const int size) -{ - double m = mean(a, size); - double sd = 0.0; - for (int i = 0; i < size; i++) { - sd += pow(a[i] - m, 2); - } - sd = sqrt(sd / (size - 1)); - return sd; -} - -double cov(const double x[], const double y[], const int size){ - - double covariance = 0; - - double meanX = mean(x, size); - double meanY = mean(y, size); - - for(int i = 0; i < size; i++){ - // double xi =x[i]; - // double yi =y[i]; - covariance += (x[i] - meanX) * (y[i] - meanY); - - } - - return covariance/(size-1); - -} - -double cov_mean(const double x[], const double y[], const int size){ - - double covariance = 0; - - for(int i = 0; i < size; i++){ - // double xi =x[i]; - // double yi =y[i]; - covariance += x[i] * y[i]; - - } - - return covariance/size; - -} - -double corr(const double x[], const double y[], const int size){ - - double nom = 0; - double denomX = 0; - double denomY = 0; - - double meanX = mean(x, size); - double meanY = mean(y, size); - - for(int i = 0; i < size; i++){ - nom += (x[i] - meanX) * (y[i] - meanY); - denomX += (x[i] - meanX) * (x[i] - meanX); - denomY += (y[i] - meanY) * (y[i] - meanY); - - //printf("x[%i]=%1.3f, y[%i]=%1.3f, nom[%i]=%1.3f, denomX[%i]=%1.3f, denomY[%i]=%1.3f\n", i, x[i], i, y[i], i, nom, i, denomX, i, denomY); - } - - return nom/sqrt(denomX * denomY); - -} - -double autocorr_lag(const double x[], const int size, const int lag){ - - return corr(x, &(x[lag]), size-lag); - -} - -double autocov_lag(const double x[], const int size, const int lag){ - - return cov_mean(x, &(x[lag]), size-lag); - -} - -void zscore_norm(double a[], int size) -{ - double m = mean(a, size); - double sd = stddev(a, size); - for (int i = 0; i < size; i++) { - a[i] = (a[i] - m) / sd; - } - return; -} - -void zscore_norm2(const double a[], const int size, double b[]) -{ - double m = mean(a, size); - double sd = stddev(a, size); - for (int i = 0; i < size; i++) { - b[i] = (a[i] - m) / sd; - } - return; -} - -double moment(const double a[], const int size, const int start, const int end, const int r) -{ - int win_size = end - start + 1; - a += start; - double m = mean(a, win_size); - double mr = 0.0; - for (int i = 0; i < win_size; i++) { - mr += pow(a[i] - m, r); - } - mr /= win_size; - mr /= stddev(a, win_size); //normalize - return mr; -} - -void diff(const double a[], const int size, double b[]) -{ - for (int i = 1; i < size; i++) { - b[i - 1] = a[i] - a[i - 1]; - } -} - -int linreg(const int n, const double x[], const double y[], double* m, double* b) //, double* r) -{ - double sumx = 0.0; /* sum of x */ - double sumx2 = 0.0; /* sum of x**2 */ - double sumxy = 0.0; /* sum of x * y */ - double sumy = 0.0; /* sum of y */ - double sumy2 = 0.0; /* sum of y**2 */ - - /* - for (int i = 0; i < n; i++) - { - fprintf(stdout, "x[%i] = %f, y[%i] = %f\n", i, x[i], i, y[i]); - } - */ - - for (int i=0;i -#include -#include - -extern double max_(const double a[], const int size); -extern double min_(const double a[], const int size); -extern double mean(const double a[], const int size); -extern double sum(const double a[], const int size); -extern void cumsum(const double a[], const int size, double b[]); -extern void icumsum(const int a[], const int size, int b[]); -extern double isum(const int a[], const int size); -extern double median(const double a[], const int size); -extern double stddev(const double a[], const int size); -extern double corr(const double x[], const double y[], const int size); -extern double cov(const double x[], const double y[], const int size); -extern double cov_mean(const double x[], const double y[], const int size); -extern double autocorr_lag(const double x[], const int size, const int lag); -extern double autocov_lag(const double x[], const int size, const int lag); -extern void zscore_norm(double a[], int size); -extern void zscore_norm2(const double a[], const int size, double b[]); -extern double moment(const double a[], const int size, const int start, const int end, const int r); -extern void diff(const double a[], const int size, double b[]); -extern int linreg(const int n, const double x[], const double y[], double* m, double* b); //, double* r); -extern double norm_(const double a[], const int size); - -#endif diff --git a/src/catch22 b/src/catch22 new file mode 160000 index 0000000..2e1a271 --- /dev/null +++ b/src/catch22 @@ -0,0 +1 @@ +Subproject commit 2e1a271c6a7437b6a4a754e1adc7e34d7a224c01 diff --git a/src/C/catch22_wrap.c b/src/wrapper/catch22_wrap.c similarity index 100% rename from src/C/catch22_wrap.c rename to src/wrapper/catch22_wrap.c