/*
  Subject: Non linear time series analysis of ph recording at LES
  Programmer: Stathis Kamperis
  Email: stathisk@med.auth.gr

  Added from fork
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "definitions.h"
#include "etc.h"
#include "kernel.h"
#include "parameters.h"
#include "patient.h"

/* Global variables */
long int start = START;
long int end = END;
float sample_rate = SAMPLE_RATE;
float corr_threshold = CORR_THRESHOLD;
float fnn_threshold = FNN_THRESHOLD;
int min_embedding = MIN_EMBEDDING;
int max_embedding = MAX_EMBEDDING;
int study_id = -1;

int iDelay;
int iEmbedDim;

/* Function prototypes */
void reconstruct_phase_space();
void run_surrogates_test();
void calc_max_lyap_exp();
void full_output_results();

/*******************************************************************************
*                                                                              *
*    Phase space reconstruction                                                *
*                                                                              *
*    1) Calculate auto-correlation function                                    *
*    2) Extract delay                                                          *
*    3) Run false nearest neighbors algorithm                                  *
*    4) Extract embedding dimension                                            *
*                                                                              *
*******************************************************************************/
void
reconstruct_phase_space() {
    task my_task;
    
    /* Calculate auto-correlation function */
    strcpy(my_task.info, "Calculating auto-correlation function...");
    sprintf(my_task.cmd, "./autocor -V0 ph.txt.ex.rs.stat -o");
    exec_task(my_task);
    
    /* Extract delay */
    iDelay = get_delay("ph.txt.ex.rs.stat_co");
    if (iDelay == ERROR_OCCURED) {
        fprintf(stderr, "Could not determine a valid time delay.\n");

        /* Arrange directories */
        strcpy(my_task.info, "Arranging directories...");
        sprintf(my_task.cmd, "rm -rf %d; mkdir %d/", study_id, study_id);
        exec_task(my_task);

        /* Produce graph plots */
        strcpy(my_task.info, "Producing graphing plots...");
        sprintf(my_task.cmd, \
        "gnuplot plot.1 > %d/1.png; \
         gnuplot plot.2 > %d/2.png; \
         gnuplot plot.3 > %d/3.png; \
         gnuplot plot.his > %d/his.png", \
         study_id, study_id, study_id, study_id);
        exec_task(my_task);
        
        /* Exit program */
        exit(0);
    }
    printf("Delay: %d\n\n", iDelay);
   
    /* Run false nearest neighbors algorithm
        Note: 2>&1 redirects stderr to stdout */
    strcpy(my_task.info, "Running false nearest neighbors algorithm...");    
    sprintf(my_task.cmd, "./false_nearest -V0 -d%d -m%d -M%d ph.txt.ex.rs.stat -o", iDelay, min_embedding, max_embedding);
    exec_task(my_task);
    
    /* Extract embedding */
    iEmbedDim = get_embedding("ph.txt.ex.rs.stat.fnn");
    if (iEmbedDim == ERROR_OCCURED) {
        fprintf(stderr, "Could not determine a valid embedding dimension. Using max_embedding(%d) instead.\n", \
        max_embedding);
        iEmbedDim = max_embedding;
    }
    printf("Embedding: %d\n\n", iEmbedDim);
}

/*******************************************************************************
*                                                                              *
*    Surrogates data test                                                      *
*                                                                              *
*    1) Determine end-to-end mismatch                                          *
*    2) Extract sub-sequence (offset, length) with a low end effect            *
*    3) Create surrogates data                                                 *
*    4) Run actual surrogates data test                                        *
*    5) Read the result                                                        *
*                                                                              *
*******************************************************************************/
void
run_surrogates_test() {
    int offset, length;
    task my_task;
    
    /* Determine end-to-end mismatch before making surrogating data */
    strcpy(my_task.info, "Determining end-to-end mismatch...");
    sprintf(my_task.cmd, "./endtoend -V0 ph.txt.ex.rs.stat -o");
    exec_task(my_task);
    
    /* Extract sub-sequence (offset, length) with a low end effect */
    get_subsequence("ph.txt.ex.rs.stat_end", &offset, &length);
    
    strcpy(my_task.info, "Extracting sub-sequence with a low end effect");
    sprintf(my_task.cmd, "./choose -V0 -o ph.txt.ex.rs.stat.sel -x%d -l%d ph.txt.ex.rs.stat", offset, length);
    exec_task(my_task);

    /* Create surrogates data */
    strcpy(my_task.info, "Creating surrogates data...");
    sprintf(my_task.cmd, "./surrogates -V0 -n39 ph.txt.ex.rs.stat.sel -o");
    exec_task(my_task);
    
    /* Run actual surrogates data test */
    strcpy(my_task.info, "Running actual surrogates data test...");
    sprintf(my_task.cmd, "./timerev -V0 -d%d ph.txt.ex.rs.stat.sel ph.txt.ex.rs.stat.sel_surr_??? &> surrogates.dat", iDelay);
    exec_task(my_task);
    
    /* Read the result */
    if (sort_surrogates("surrogates.dat", 39) == SURR_SUCCESS)
        printf("PASSED!\n");
    else {
        fprintf(stderr, "FAILED!\n\n");
        
        /* Arrange directories */
        strcpy(my_task.info, "Arranging directories...");
        sprintf(my_task.cmd, "rm -rf %d; mkdir %d; cp surrogates.dat* %d/", study_id, study_id, study_id);
        exec_task(my_task);

        /* Produce graph plots */
        strcpy(my_task.info, "Producing graphing plots...");
        sprintf(my_task.cmd,             \
        "gnuplot plot.1 > %d/1.png;      \
         gnuplot plot.2 > %d/2.png;      \
         gnuplot plot.3 > %d/3.png;      \
         gnuplot plot.his > %d/his.png;  \
         gnuplot plot.fnn > %d/fnn.png", \
         study_id, study_id, study_id, study_id, study_id);
        exec_task(my_task);
        
        /* Exit program */
        exit(0);
    }
}

/*******************************************************************************
*                                                                              *
*    Calculate maximum lyapunov exponent                                       *
*                                                                              *
*    1) Calculate the whole lyapunov spectrum                                  *
*    2) Extract the maximum lyapunov exponent                                  *
*                                                                              *
*******************************************************************************/
void
calc_max_lyap_exp() {
    float fMaxLyapunov;
    task my_task;

    /* Calculate the whole lyapunov spectrum */
    strcpy(my_task.info, "Calculating lyapunov spectrum...");
    sprintf(my_task.cmd, "./lyap_spec -V0 -d%d -m1,%d ph.txt.ex.rs.stat -o", iDelay, iEmbedDim);    
    exec_task(my_task);
    
    /* Extract the maximum lyapunov exponent */   
    fMaxLyapunov = get_max_lyapunov("ph.txt.ex.rs.stat.lyaps");
    printf("Max lyapunov exponent: %f\n\n", fMaxLyapunov);
}

/*******************************************************************************
*                                                                              *
*    Arrange directories and produce graph plots                               *
*                                                                              *
*   1) Arrange directories                                                     *
*   2) Produce graph plots                                                     *
*                                                                              *
*******************************************************************************/
void
dirs_and_plots() {
    task my_task;
    
    /* Arrange directories */
    strcpy(my_task.info, "Arranging directories...");
    sprintf(my_task.cmd, "rm -rf %d; mkdir %d; cp rms.* surrogates.dat* %d/", study_id, study_id, study_id);
    exec_task(my_task);

    /* Produce graph plots */
    strcpy(my_task.info, "Producing graphing plots...");
    sprintf(my_task.cmd, \
    "gnuplot plot.1 > %d/1.png; \
     gnuplot plot.2 > %d/2.png; \
     gnuplot plot.3 > %d/3.png; \
     gnuplot plot.his > %d/his.png; \
     gnuplot plot.fnn > %d/fnn.png; \
     gnuplot plot.lyap > %d/lyap.png; \
     gnuplot plot.d2 > %d/d2.png; \
     gnuplot plot.avd2 > %d/avd2.png; \
     gnuplot plot.rec1 > %d/rec1.png; \
     gnuplot plot.rec2 > %d/rec2.png", \
     study_id, study_id, study_id, study_id, study_id, study_id, study_id, study_id, study_id, study_id);
    exec_task(my_task);   
}

/*******************************************************************************
*                                                                              *
*    Main function                                                             *
*                                                                              *
*******************************************************************************/
int
main(int argc, char **argv) {
    int i, j;
    float fMaxLyapunov, results[4], logr[10000], slopes[10000];
    task my_task;
    
    /* Clear screen */
    system("clear");

    /* Extract supine repiod */
    get_supine_period("ph.txt");
    
    /* Get parameters */
    if (read_protocol("ph.prot") == PROTOCOL_NOT_FOUND) {
        get_parameters(argc, argv);
    }
        
    /* Print patient data */
    print_patient_data("ph.txt");
    
    /* Print parameters */
    print_parameters();    
    
    /* Clean old files from previous executions. */
    strcpy(my_task.info, "Cleaning old files...");
    sprintf(my_task.cmd, "rm -rf ph.txt.* rms.* surrogates.*");
    exec_task(my_task);
    
    /* Extract data from file */
    printf("Extracting data from pH file...\n\n");
    extract_data("ph.txt");
        
    /* Cut out and resample
    Note: Given an erroneous filename, `resample` reads from the stdin
    and, thus, "pauses" the execution. */
    strcpy(my_task.info, "Cutting out and resampling...");
    sprintf(my_task.cmd, "./resample -V0 -s%f -x%ld -l%ld %s -o", sample_rate, start, end - start, "ph.txt.ex");
    exec_task(my_task);
    
    /* Calculate histogram */
    strcpy(my_task.info, "Calculating histogram...");
    sprintf(my_task.cmd, "./histogram -V0 ph.txt.ex.rs -o");
    exec_task(my_task);
    
    /* Convert time series to stationary... */
    printf("Converting time series to stationary...\n\n");
    make_stationary("ph.txt.ex.rs");    
    
    /* Calculate recurrent plot */
    strcpy(my_task.info, "Calculating recurrent plot of non-stationary time series...");
    sprintf(my_task.cmd, "./recurr -V0 ph.txt.ex.rs -o");
    exec_task(my_task);
    strcpy(my_task.info, "Calculating recurrent plot of stationary time series...");
    sprintf(my_task.cmd, "./recurr -V0 ph.txt.ex.rs.stat -o");
    exec_task(my_task);
    
    /* Reconstruct phase space */
    reconstruct_phase_space();
    
    /* Run surrogates data test */
    run_surrogates_test();
    
    /* Calculate maximux lyapunov exponent */
    calc_max_lyap_exp();

    /* Remove noise */
    strcpy(my_task.info, "Removing noise...");
    sprintf(my_task.cmd, "./project -V0 -m%d -q%d -r0.1 -k40 -i1 ph.txt.ex.rs.stat -o", iEmbedDim, iEmbedDim/2);
    exec_task(my_task);

    /* Calculate rms errors */
    strcpy(my_task.info, "Calculating rms error... (1/2)");
    sprintf(my_task.cmd, "./rms -V0  ph.txt.ex.rs.stat > rms.1");
    exec_task(my_task);

    strcpy(my_task.info, "Calcuclating rms error... (2/2)");
    sprintf(my_task.cmd, "./rms -V0 ph.txt.ex.rs.stat_c > rms.2");
    exec_task(my_task);
    
    /* Estimate d2 correlation dimension */
    strcpy(my_task.info, "Estimating d2 correlation dimension...");
    sprintf(my_task.cmd, "./d2 -V0 -d%d -M1,%d -r0.001 -#100 ph.txt.ex.rs.stat -o", iDelay, max_embedding); // embedding?
    exec_task(my_task);

    sprintf(my_task.cmd, "./d2 -V0 -d%d -M1,%d -r0.001 -#100 ph.txt.ex.rs.stat_c -o", iDelay, max_embedding); // embedding?
    exec_task(my_task);

    /* Smooth d2 correlation dimension */
    strcpy(my_task.info, "Smoothing d2 correlation dimension...");
    sprintf(my_task.cmd, "./av-d2 -V0 ph.txt.ex.rs.stat.d2 -o");
    exec_task(my_task);
    
    sprintf(my_task.cmd, "./av-d2 -V0 ph.txt.ex.rs.stat_c.d2 -o");
    exec_task(my_task);
    
    /* Estimate scaling region */
    printf("Estimating scaling region...\n");
    
    j = get_d2_data("ph.txt.ex.rs.stat_c.d2", iEmbedDim, logr, slopes);
    
    /*
    The first column of .d2 data contains the r-values or the log(r)-values?
    for (i = 0; i < j; i++)
        logr[i] = log(logr[i]); */

    estimate_scaling_region(logr, slopes, j, 7, 0.8, results);
    
    printf("Mean: %f\tStd: %f\tlogr1: %f\t logr2: %f\n\n", results[0], results[1], results[2], results[3]);
    
    /* Arrange dirs and produce graph plots */
    dirs_and_plots();
    
    return 0;
}