Merge branch 'master' of github.com:chlubba/catch22

C/CO_AutoCorr.c

@@ -1,4 +1,15 @@
-#include <complex.h>
+#if __cplusplus
+#   include <complex>
+typedef std::complex< double > cplx;
+#else
+#   include <complex.h>
+#if defined(__GNUC__) || defined(__GNUG__)
+typedef double complex cplx;
+#elif defined(_MSC_VER)
+typedef _Dcomplex cplx;
+#endif
+#endif
+
 #include <math.h>
 #include <stdio.h>
 #include <string.h>
@@ -8,11 +19,12 @@
 #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)
-typedef double complex cplx;
 
 int nextpow2(int n)
 {
@@ -47,7 +59,7 @@ static void apply_conj(cplx a[], int size, int normalize)
 void dot_multiply(cplx a[], cplx b[], int size)
 {
     for (int i = 0; i < size; i++) {
-        a[i] = a[i] * conj(b[i]);
+        a[i] = _Cmulcc(a[i], conj(b[i]));
     }
 }
 
@@ -60,10 +72,23 @@ double * CO_AutoCorr(const double y[], const int size, const int tau[], const in
     cplx * F = malloc(nFFT * sizeof *F);
     cplx * tw = malloc(nFFT * sizeof *tw);
     for (int i = 0; i < size; i++) {
-        F[i] = CMPLX(y[i] - m, 0.0);
+
+        #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++) {
-        F[i] = CMPLX(0.0, 0.0);
+        #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;
 
@@ -73,7 +98,8 @@ double * CO_AutoCorr(const double y[], const int size, const int tau[], const in
     fft(F, nFFT, tw);
     cplx divisor = F[0];
     for (int i = 0; i < nFFT; i++) {
-        F[i] = F[i] / divisor;
+        //F[i] = F[i] / divisor;
+        F[i] = _Cdivcc(F[i], divisor);
     }
 
     double * out = malloc(tau_size * sizeof(out));
@@ -91,13 +117,25 @@ double * co_autocorrs(const double y[], const int size)
     m = mean(y, size);
     nFFT = nextpow2(size) << 1;
 
-    cplx * F = malloc(nFFT * sizeof *F);
-    cplx * tw = malloc(nFFT * sizeof *tw);
+    cplx * F = malloc(nFFT * 2 * sizeof *F);
+    cplx * tw = malloc(nFFT * 2 * sizeof *tw);
     for (int i = 0; i < size; i++) {
-        F[i] = CMPLX(y[i] - m, 0.0);
+
+        #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++) {
-        F[i] = CMPLX(0.0, 0.0);
+
+        #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;
 
@@ -107,10 +145,10 @@ double * co_autocorrs(const double y[], const int size)
     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); // F[i] / divisor;
     }
 
-    double * out = malloc(nFFT * sizeof(out));
+    double * out = malloc(nFFT * 2 * sizeof(out));
     for (int i = 0; i < nFFT; i++) {
         out[i] = creal(F[i]);
     }
@@ -367,6 +405,9 @@ double CO_trev_1_num(const double y[], const int size)
     return out;
 }
 
+#define tau 2
+#define numBins 5
+
 double CO_HistogramAMI_even_2_5(const double y[], const int size)
 {
 
@@ -379,8 +420,8 @@ double CO_HistogramAMI_even_2_5(const double y[], const int size)
         }
     }
 
-    const int tau = 2;
-    const int numBins = 5;
+    //const int tau = 2;
+    //const int numBins = 5;
 
     double * y1 = malloc((size-tau) * sizeof(double));
     double * y2 = malloc((size-tau) * sizeof(double));
@@ -391,8 +432,8 @@ double CO_HistogramAMI_even_2_5(const double y[], const int size)
     }
 
     // set bin edges
-    double maxValue = max(y, size);
-    double minValue = min(y, size);
+    const double maxValue = max_(y, size);
+    const double minValue = min_(y, size);
 
     double binStep = (maxValue - minValue + 0.2)/5;
     //double binEdges[numBins+1] = {0};

C/CO_AutoCorr.h

@@ -1,6 +1,18 @@
 #ifndef CO_AUTOCORR_H
 #define CO_AUTOCORR_H
-#include <complex.h>
+
+#if __cplusplus
+#   include <complex>
+typedef std::complex< double > cplx;
+#else
+#   include <complex.h>
+#if defined(__GNUC__) || defined(__GNUG__)
+typedef double complex cplx;
+#elif defined(_MSC_VER)
+typedef _Dcomplex cplx;
+#endif
+#endif
+
 #include <math.h>
 #include <stdio.h>
 #include <string.h>

C/DN_OutlierInclude.c

@@ -43,7 +43,7 @@ double DN_OutlierInclude_np_001_mdrmd(const double y[], const int size, const in
     if(constantFlag) return 0; // if constant, return 0
 
     // find maximum (or minimum, depending on sign)
-    double maxVal = max(yWork, size);
+    double maxVal = max_(yWork, size);
 
     // maximum value too small? return 0
     if(maxVal < inc){

C/SB_CoarseGrain.c

@@ -19,8 +19,8 @@ void sb_coarsegrain(const double y[], const int size, const char how[], const in
     }
     */
 
-    double * th = malloc((num_groups + 1) * sizeof(th));
-    double * ls = malloc((num_groups + 1) * sizeof(th));
+    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]);

C/SB_MotifThree.c

@@ -52,23 +52,44 @@ double SB_MotifThree_quantile_hh(const double y[], const int size)
     // 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] - 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]));
+        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]));
+            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,
@@ -78,37 +99,46 @@ double SB_MotifThree_quantile_hh(const double y[], const int size)
                 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);
                 }
             }
-
-            out2[i][j] = (double)sizes_r2[i][j] / (size - 1);
+            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(sizes_r1);
+
     // free nested array
     for (int i = 0; i < alphabet_size; i++) {
         free(r1[i]);
     }
     free(r1);
-    free(sizes_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]);
         }
-    }
-    for (int i = 0; i < alphabet_size; i++) {
         free(r2[i]);
-        free(sizes_r2[i]);
-        free(out2[i]);
     }
+
     free(r2);
     free(sizes_r2);
     free(out2);

C/SC_FluctAnal.c

@@ -108,7 +108,7 @@ double SC_FluctAnal_2_50_1_logi_prop_r1(const double y[], const int size, const
             }
 
             if (strcmp(how, "rsrangefit") == 0) {
-                F[i] += pow(max(buffer, tau[i]) - min(buffer, tau[i]), 2);
+                F[i] += pow(max_(buffer, tau[i]) - min_(buffer, tau[i]), 2);
             }
             else if (strcmp(how, "dfa") == 0) {
                 for(int k = 0; k<tau[i]; k++){
@@ -162,19 +162,19 @@ double SC_FluctAnal_2_50_1_logi_prop_r1(const double y[], const int size, const
             buffer[j] = logtt[j] * m1 + b1 - logFF[j];
         }
 
-        sserr[i - minPoints] += norm(buffer, i);
+        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);
+        sserr[i - minPoints] += norm_(buffer, ntt-i+1);
 
     }
 
     double firstMinInd = 0.0;
-    double minimum = min(sserr, nsserr);
+    double minimum = min_(sserr, nsserr);
     for(int i = 0; i < nsserr; i++)
     {
         if(sserr[i] == minimum)

C/SP_Summaries.c

@@ -18,7 +18,7 @@ int welch(const double y[], const int size, const int NFFT, const double Fs, con
     int k = floor((double)size/((double)windowWidth/2.0))-1;
 
     // normalising scale factor
-    double KMU = k * pow(norm(window, windowWidth),2);
+    double KMU = k * pow(norm_(window, windowWidth),2);
 
     double * P = malloc(NFFT * sizeof(double));
     for(int i = 0; i < NFFT; i++){
@@ -40,10 +40,13 @@ int welch(const double y[], const int size, const int NFFT, const double Fs, con
 
         // initialise F (
         for (int i = 0; i < windowWidth; i++) {
-            F[i] = CMPLX(xw[i] - m, 0.0);
+            cplx tmp = { xw[i] - m, 0.0 };
+            F[i] = tmp; // CMPLX(xw[i] - m, 0.0);
         }
         for (int i = windowWidth; i < NFFT; i++) {
-            F[i] = CMPLX(0.0, 0.0);
+            // F[i] = CMPLX(0.0, 0.0);
+            cplx tmp = { 0.0, 0.0 };
+            F[i] = tmp;
         }
 
         fft(F, NFFT, tw);
@@ -128,9 +131,10 @@ double SP_Summaries_welch_rect(const double y[], const int size, const char what
     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*M_PI*f[i];
-        Sw[i] = S[i]/(2*M_PI);
+        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;