Merge pull request #23 from DynamicsAndNeuralSystems/catch24-interpol…

…ation-dev

`catch24` and `CO_f1ecac` feature interpolation

.vscode/settings.json

@@ -0,0 +1,6 @@
+{
+    "files.associations": {
+        "*.rmd": "markdown",
+        "dn_outlierinclude.h": "c"
+    }
+}

C/CO_AutoCorr.c

@@ -176,7 +176,7 @@ int co_firstzero(const double y[], const int size, const int maxtau)
 
 }
 
-int CO_f1ecac(const double y[], const int size)
+double CO_f1ecac(const double y[], const int size)
 {
 
     // NaN check
@@ -190,17 +190,19 @@ int CO_f1ecac(const double y[], const int size)
 
     // compute autocorrelations
     double * autocorrs = co_autocorrs(y, size);
-    
+
     // threshold to cross
     double thresh = 1.0/exp(1);
 
-    int out = size;
-    for(int i = 0; i < size-1; i++){
-
-        // printf("autocorrs_i: %1.3f autocorrs_i+1: %1.3f\n", autocorrs[i], autocorrs[i+1]);
-
-        if ((autocorrs[i] - thresh)*(autocorrs[i+1] - thresh) < 0){
-            out = i + 1;
+    double out = (double)size;
+    for(int i = 1; i < size-1; i++){
+        // printf("i=%d autocorrs_i=%1.3f\n", i, autocorrs[i]);
+        if ( autocorrs[i] < thresh ){
+            double m = autocorrs[i] - autocorrs[i-1];
+            double dy = thresh - autocorrs[i-1];
+            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;
         }

C/CO_AutoCorr.h

@@ -29,7 +29,7 @@ extern double CO_Embed2_Basic_tau_incircle(const double y[], const int size, con
 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 int CO_f1ecac(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

C/DN_Mean.c

@@ -0,0 +1,11 @@
+#include <stdio.h>
+
+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;
+}

C/DN_Mean.h

@@ -0,0 +1,12 @@
+//
+//  Created by Trent Henderson 27 September 2021
+//
+
+#ifndef DN_MEAN
+#define DN_MEAN
+
+#include <stdio.h>
+
+extern double DN_Mean(const double a[], const int size);
+
+#endif /* DN_MEAN */

C/DN_Spread_Std.c

@@ -0,0 +1,13 @@
+#include <stdio.h>
+#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;
+}

C/DN_Spread_Std.h

@@ -0,0 +1,12 @@
+//
+//  Created by Trent Henderson 27 September 2021
+//
+
+#ifndef DN_SPREADSTD
+#define DN_SPREADSTD
+
+#include <stdio.h>
+
+extern double DN_Spread_Std(const double a[], const int size);
+
+#endif /* DN_SPREADSTD */

C/main.c

@@ -3,10 +3,13 @@
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <stdbool.h>
 //#include <dirent.h>
 
 #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"
@@ -45,7 +48,7 @@ int quality_check(const double y[], const int size)
     return 0;
 }
 
-void run_features(double y[], int size, FILE * outfile)
+void run_features(double y[], int size, FILE * outfile, bool catch24)
 {
     int quality = quality_check(y, size);
     if(quality != 0)
@@ -66,6 +69,18 @@ void run_features(double y[], int size, FILE * outfile)
     // 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);
@@ -108,18 +123,6 @@ void run_features(double y[], int size, FILE * outfile)
     timeTaken = (double)(clock()-begin)*1000/CLOCKS_PER_SEC;
     fprintf(outfile, "%.14f, %s, %f\n", result, "CO_trev_1_num", 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_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 = FC_LocalSimple_mean1_tauresrat(y_zscored, size);
@@ -198,6 +201,23 @@ void run_features(double y[], int size, FILE * outfile)
     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);
@@ -282,7 +302,19 @@ int main(int argc, char * argv[])
     fclose(infile);
     y = realloc(y, size * sizeof *y);
     //printf("size=%i\n", size);
-    run_features(y, size, outfile);
+
+    // 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);
 
@@ -325,52 +357,63 @@ int main(int argc, char * argv[])
     zscore_norm(y, size);
 
     double result;
-
+
+    result = DN_OutlierInclude_n_001_mdrmd(y, size);
+    printf("DN_OutlierInclude_n_001_mdrmd: %1.5f\n", result);
+    result = DN_OutlierInclude_p_001_mdrmd(y, size);
+    printf("DN_OutlierInclude_p_001_mdrmd: %1.5f\n", result);
     result = DN_HistogramMode_5(y, size);
     printf("DN_HistogramMode_5: %1.3f\n", result);
     result = DN_HistogramMode_10(y, size);
     printf("DN_HistogramMode_10: %1.3f\n", result);
-    result = CO_Embed2_Dist_tau_d_expfit_meandiff(y, size);
-    printf("CO_Embed2_Dist_tau_d_expfit_meandiff: %1.3f\n", result);
-    result = CO_f1ecac(y, size);
-    printf("CO_f1ecac: %1.f\n", result);
-    result = CO_FirstMin_ac(y, size);
-    printf("CO_FirstMin_ac: %1.f\n", result);
-    result = CO_HistogramAMI_even_2_5(y, size);
-    printf("CO_HistogramAMI_even_2_5: %1.3f\n", result);
-    result = CO_trev_1_num(y, size);
-    printf("CO_trev_1_num: %1.5f\n", result);
-    result = DN_OutlierInclude_p_001_mdrmd(y, size);
-    printf("DN_OutlierInclude_p_001_mdrmd: %1.5f\n", result);
-    result = DN_OutlierInclude_n_001_mdrmd(y, size);
-    printf("DN_OutlierInclude_n_001_mdrmd: %1.5f\n", result);
+    result = SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1(y, size);
+    printf("SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1: %1.5f\n", result);
+    result = SB_TransitionMatrix_3ac_sumdiagcov(y, size);
+    printf("SB_TransitionMatrix_3ac_sumdiagcov: %1.5f\n", result);
     result = FC_LocalSimple_mean1_tauresrat(y, size);
     printf("FC_LocalSimple_mean1_tauresrat: %1.5f\n", result);
-    result = FC_LocalSimple_mean3_stderr(y, size);
-    printf("FC_LocalSimple_mean3_stderr: %1.5f\n", result);
-    result = IN_AutoMutualInfoStats_40_gaussian_fmmi(y, size);
-    printf("IN_AutoMutualInfoStats_40_gaussian_fmmi: %1.5f\n", result);
-    result = MD_hrv_classic_pnn40(y, size);
-    printf("MD_hrv_classic_pnn40: %1.5f\n", result);
+    result = SB_MotifThree_quantile_hh(y, size);
+    printf("SB_MotifThree_quantile_hh: %1.5f\n", result);
+    result = CO_HistogramAMI_even_2_5(y, size);
+    printf("CO_HistogramAMI_even_2_5: %1.3f\n", result);
+    result = CO_Embed2_Dist_tau_d_expfit_meandiff(y, size);
+    printf("CO_Embed2_Dist_tau_d_expfit_meandiff: %1.3f\n", result);
     result = SB_BinaryStats_diff_longstretch0(y, size);
     printf("SB_BinaryStats_diff_longstretch0: %1.5f\n", result);
+    result = MD_hrv_classic_pnn40(y, size);
+    printf("MD_hrv_classic_pnn40: %1.5f\n", result);
     result = SB_BinaryStats_mean_longstretch1(y, size);
     printf("SB_BinaryStats_mean_longstretch1: %1.5f\n", result);
-    result = SB_MotifThree_quantile_hh(y, size);
-    printf("SB_MotifThree_quantile_hh: %1.5f\n", result);
-    result = SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1(y, size);
-    printf("SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1: %1.5f\n", result);
-    result = SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1(y, size);
-    printf("SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1: %1.5f\n", result);
+    result = FC_LocalSimple_mean3_stderr(y, size);
+    printf("FC_LocalSimple_mean3_stderr: %1.5f\n", result);
     result = SP_Summaries_welch_rect_area_5_1(y, size);
     printf("SP_Summaries_welch_rect_area_5_1: %1.5f\n", result);
     result = SP_Summaries_welch_rect_centroid(y, size);
     printf("SP_Summaries_welch_rect_centroid: %1.5f\n", result);
-    result = SB_TransitionMatrix_3ac_sumdiagcov(y, size);
-    printf("SB_TransitionMatrix_3ac_sumdiagcov: %1.5f\n", result);
+    result = CO_f1ecac(y, size);
+    printf("CO_f1ecac: %1.f\n", result);
+    result = CO_FirstMin_ac(y, size);
+    printf("CO_FirstMin_ac: %1.f\n", result);
+    result = IN_AutoMutualInfoStats_40_gaussian_fmmi(y, size);
+    printf("IN_AutoMutualInfoStats_40_gaussian_fmmi: %1.5f\n", result);
     result = PD_PeriodicityWang_th0_01(y, size);
     printf("PD_PeriodicityWang_th0_01: %1.f\n", result);
+    result = SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1(y, size);
+    printf("SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1: %1.5f\n", result);
+    result = CO_trev_1_num(y, size);
+    printf("CO_trev_1_num: %1.5f\n", result);
 
+    int catch24;
+    printf("Do you want to run catch24? Enter 0 for catch22 or 1 for catch24.");
+    scanf("%d", &catch24);
+
+    if (catch24 == 1) {
+        result = DN_Mean(y, size);
+        printf("DN_Mean: %1.5f\n", result);
+        result = DN_Spread_Std(y, size);
+        printf("DN_Spread_Std: %1.5f\n", result);
+    } else {
+    }
 
   return 0;
 }

Matlab/DN_Mean.m

@@ -0,0 +1,6 @@
+function out = DN_Mean(y)
+
+    % no combination of single functions
+    coder.inline('never');
+
+    out = DN_Mean(y, 5);

Matlab/DN_Spread_Std.m

@@ -0,0 +1,6 @@
+function out = DN_Spread_Std(y)
+
+    % no combination of single functions
+    coder.inline('never');
+
+    out = DN_Spread_Std(y, 5);

README.md

@@ -33,5 +33,5 @@ And Julia users can use [this Julia package](https://github.com/brendanjohnharri
 
 - See language-specific usage information in the [wiki](https://github.com/chlubba/catch22/wiki/Installation-and-Testing).
 - __Important Note:__ _catch22_ features only evaluate _dynamical_ properties of time series and do not respond to basic differences in the location (e.g., mean) or spread (e.g., variance).
-  - If you think features of the raw distribution may be important for your application, we suggest you add them (in the simplest case, two additional features: the mean and standard deviation) to this feature set.
+  - If you think features of the raw distribution may be important for your application, we suggest you apply the function argument `catch24 = true` (`TRUE` in R, `True` in Python) to your call to the _catch22_ function in the language of your choice.
 - Note that time series are _z_-scored internally which means e.g., constant time series will lead to `NaN` outputs.

featureList.txt

@@ -1,5 +1,7 @@
 DN_HistogramMode_5
 DN_HistogramMode_10
+DN_Mean
+DN_Spread_Std
 CO_f1ecac
 CO_FirstMin_ac
 CO_HistogramAMI_even_2_5

wrap_Matlab/M_wrapper.c

@@ -3,6 +3,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"

wrap_Matlab/catch22_CO_Embed2_Dist_tau_d_expfit_meandiff.c

@@ -10,6 +10,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"

wrap_Matlab/catch22_CO_FirstMin_ac.c

@@ -10,6 +10,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"

wrap_Matlab/catch22_CO_HistogramAMI_even_2_5.c

@@ -10,6 +10,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"

wrap_Matlab/catch22_CO_f1ecac.c

@@ -10,6 +10,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"

wrap_Matlab/catch22_CO_trev_1_num.c

@@ -10,6 +10,8 @@
 #include "CO_AutoCorr.h"
 #include "DN_HistogramMode_10.h"
 #include "DN_HistogramMode_5.h"
+#include "DN_Mean.h"
+#include "DN_Spread_Std.h"
 #include "DN_OutlierInclude.h"
 #include "FC_LocalSimple.h"
 #include "IN_AutoMutualInfoStats.h"