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

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

extern long int start;
extern long int end;
extern int sample_rate;
extern float ac_threshold;
extern float fnn_threshold;
extern int min_embedding;
extern int max_embedding;

/*
    Extracts maximum supine period during the 24ph recording.
    This feature is used when we want to substract the effect of
    patient's motility to our system.

    Caution: all operations are referred to the original ph file
    and not to the resampled one!
*/
void
get_supine_period(char *fname) {
    int cnt, flag;
    long int iRefPoint, iStart, iEnd, diff, diffmax, iStartMax;
    char dummy[MAX_LEN], dummy2[MAX_LEN];
    char sStart[MAX_LEN], sEnd[MAX_LEN];
    FILE *fp, *tmp;

    /* Open file descriptor */
    fp = fopen(fname, "r");
    if (!fp) {
        fprintf(stderr, "An error occured while trying to extract supine period from file: %s\n", fname);
        exit(0);
    }

    cnt = 0;
    diffmax = -1;
    flag = -1;
    while (fscanf(fp, "%s", dummy) != EOF) {
    /* Get the beginning time of recording. This will be used as a reference point. */
      if (!strcmp(dummy, "Start:")) {
            fscanf(fp, "%s %s", dummy, dummy);
            iRefPoint = time_to_index(dummy);
      }

        if (!strcmp(dummy, "Supine")) {
            flag = 1;
            fscanf(fp, "%s %s", sStart, sEnd);

            iStart = time_to_index(sStart);
            iEnd = time_to_index(sEnd);
    	    diff = iEnd - iStart;

            /* Have we passed midnight? If yes, add 24h ours to iEnd
                i.e. 23:00 - 08:00
                diff = (24 + 8) - 23
            */
	        if (diff < 0) {
	            diff += 24*60*60;   /* diff is in seconds */
	        }

            /* Have we found a bigger duration? If yes, capture if */
            if (diff > diffmax) {
                diffmax = diff;
                iStartMax = iStart -iRefPoint + 1;
            }
            printf("%s\t%s\t%ld\t%ld\t%ld\t%ld\n", sStart, sEnd, iStart, iEnd, diff, iRefPoint); /* Good for debugging purposes */
        }

        /* Break if done */
        if (!strcmp(dummy, "Findings:")) {
            break;
        }
    }

    iStart = iStartMax;
    iEnd = iStart + diffmax;

    printf("%ld\t%ld\n", iStart, iEnd);

    /* Close file descriptor */
    fclose(fp);

    if (flag == -1) {
        /* Bad diary - run complete test
            Should we assume the sleep period (i.e. 00:00-06:00) as the max supine period?
        */
        fprintf(stderr, "Bad diary\n");
        iStart = 0;
        iEnd = 100000;
    }

    /* Update ph.prot file accordingly*/
    /* Open file descriptors */
    fp = fopen("ph.prot", "r");
    tmp = fopen("tempfile", "w");
    if (!fp || !tmp) {
        fprintf(stderr, "An error occured while trying to update ph.prot\n");
        exit(0); // Is it too much?
    }

    while(fscanf(fp, "%s %s", dummy, dummy2) != EOF) {
        if (!strcmp(dummy, "start")) {
            fprintf(tmp, "%s %ld\n", dummy, iStart);
        }
        else if (!strcmp(dummy, "end")) {
            fprintf(tmp, "%s %ld\n", dummy, iEnd);
        }
        else {
            fprintf(tmp, "%s %s\n", dummy, dummy2);
        }
    }

    /* Close file descriptors */
    fclose(fp);
    fclose(tmp);

    /* Swap old and new file */
    remove("ph.prot");
    rename("tempfile", "ph.prot");
}

/*
.    Extracts time delay by finding where the autocorrelation function
    reaches ac_threshold.
*/
int
get_delay(char *fname) {
    int i;
    float c;
    FILE *fp;

    /* Open file descriptor */
    fp = fopen(fname, "r");
    if (!fp) {
        fprintf(stderr, "An error occured while trying to extract time delay from file: %s\n", fname);
        exit(0);
    }

    while (fscanf(fp, "%d%f", &i, &c) != EOF) {
        if (c <= ac_threshold)
            return i;
    }

    /* Close file descriptor */
    fclose(fp);
    return ERROR_OCCURED;
}

/*
    Extracts embedding dimension by finding when the percentage of false
    nearest neighbors reaches fnn_threshold.
*/
int
get_embedding(char *fname) {
    int i;
    float f, dummy;
    FILE *fp;

    /* Open file descriptor */
    fp = fopen(fname, "r");
    if (!fp) {
        fprintf(stderr, "An error occured while trying to extract embedding dimension from file: %s\n", fname);
        exit(0);
    }

    while (fscanf(fp, "%d%f%f%f", &i, &f, &dummy, &dummy) != EOF) {
        if (f <= fnn_threshold)
            return i;
    }

    /* Close file descritpor */
    fclose(fp);
    return ERROR_OCCURED;
}

/*
    Extracts the slopes from the d2 data file and returns the number of logr
    values for which the slope has been calculated.
*/
int get_d2_data(char *fname, int m, float *logr, float *slope) {
    FILE *fp;
    char tmp[128], from[128], to[128], slogr[128], sslope[128];
    int cnt;

    fp = fopen(fname, "r");

    if (!fp) {
        fprintf(stderr, "Error while extracting d2 data from file: %s\n", fname);
        exit(0);
    }

    cnt = 0;
    sprintf(from, "#dim= %d", m);
    strcpy(to, "#dim=");
    while (fgets(tmp, 64, fp) != NULL) {
        if (strstr(tmp, from)) {
            while ((!strstr(slogr, to) && fscanf(fp, "%s%s", slogr, sslope) != EOF)) {
                logr[cnt] = atof(slogr);
                slope[cnt] = atof(sslope);
                //printf("%f\t%f\n", logr[cnt], slope[cnt]);
                cnt++;
            }
        }
    }

    /* Close file descriptor */
    fclose(fp);
    return cnt;
}

/*
    Extracts sub-sequence (offset, length) with a low end effect
*/
void
get_subsequence(char *fname, int *offset, int *length) {
    char dummy[MAX_LEN];
    FILE *fp;

    /* Open file descriptor */
    fp = fopen(fname, "r");
    if (!fp) {
        fprintf(stderr, "An error occured while trying to extract sub-sequence from file: %s\n", fname);
        exit(0);
    }

    *length = -1;
    *offset = -1;
    while (fscanf(fp, "%s", dummy) != EOF) {
        if (strstr(dummy, "length") && (*length == -1)) {
            fscanf(fp, "%d", length);
        }

        if (strstr(dummy, "offset") && (*offset == -1)) {
            fscanf(fp, "%d", offset);
        }
    }

    /* Close file descriptor */
    fclose(fp);
}

/*
    Swap function used by sort_surrogates() routine.
*/
void
swap(surrogates *x, surrogates *y) {
    surrogates temp;

    temp = *x;
    *x = *y;
    *y = temp;
}

/*
    Sorts surrogates data and returns the desired result.
    Bubble sort is used as the sorting algorithm,
    which is considered efficient given the small number of items.
*/
int
sort_surrogates(char *fname, int n) {
    int i, flag;
    char new_fname[MAX_LEN];
    surrogates *my_surrogates;
    FILE *fp1, *fp2;

    strcpy(new_fname, fname);
    strcat(new_fname, ".sort");

    /* Open file descriptors */
    fp1 = fopen(fname, "r");
    fp2 = fopen(new_fname, "w");
    if (!fp1 || !fp2) {
        fprintf(stderr, "An error occured while trying to sort surrogates (%s, %s)\n", fname, new_fname);
        exit(0);
    }

    /* Allocate memory for the surrogates data */
    my_surrogates = (surrogates *) malloc((n + 1) * sizeof(surrogates));
    if (!my_surrogates) {
        fprintf(stderr, "Error allocating memory for the surrogates data.\n");
        exit(0);
    }

    /* How many are there? */
    i = 0;
    while (fscanf(fp1, "%lg%s", &my_surrogates[i].score, my_surrogates[i].fname) != EOF)
        i++;

    /* Bubble sort */
    flag = 0;
    while (!flag) {
        flag = 1;
        for (i=0; i<n; i++) {
            if (my_surrogates[i].score < my_surrogates[i+1].score) {
                swap(&my_surrogates[i], &my_surrogates[i+1]);
                flag = 0;
            }
        }
    }

    /* Save sorted surrogates data */
    for(i=0; i<n+1; i++)
        fprintf(fp2, "%lg\t%s\n", my_surrogates[i].score, my_surrogates[i].fname);

    /* Close file descriptors */
    fclose(fp1);
    fclose(fp2);

    /* Check if the test was successful.
    Don't forget to free the allocated memory! */
    if (!strstr(my_surrogates[0].fname, "surr") || !strstr(my_surrogates[n].fname, "surr")) {
        free(my_surrogates);
        return 1;
    }
    else {
        free(my_surrogates);
        return 0;
    }
}

/*
    Extracts the maximum lyapunov exponent from lyapunov spectrum file.
*/
float
get_max_lyapunov(char *fname) {
    long int i;
    float max_lyap;
    char dummy[MAX_LEN];
    FILE *fp;

    /* Open file descriptor */
    fp = fopen(fname, "r");
    if (!fp) {
        fprintf(stderr, "An error occured while trying to extract maximum lyapunov exponent from file: %s\n", fname);
        exit(0);
    }

    while (fscanf(fp, "%ld%f", &i, &max_lyap) != EOF) {
        fgets(dummy, MAX_LEN, fp);
        if (strstr(dummy, "Average"))
            break;
    }

    /* Close file descriptor */
    fclose(fp);
    return max_lyap;
}

/*
    Extracts data from pH file by skipping the unnecessary headers.
*/
void
extract_data(char *fname) {
    int flag;
    char str[MAX_LEN], new_fname[MAX_LEN];
    FILE *fp1, *fp2;

    strcpy(new_fname, fname);
    strcat(new_fname, ".ex");

    /* Open file decriptors */
    fp1 = fopen(fname, "r");
    fp2 = fopen(new_fname, "w");
    if (!fp1 || !fp2) {
        fprintf(stderr, "An error occured while trying to extract data from pH file.\n");
        exit(0);
    }

    flag = -1;
    while (fscanf(fp1, "%s", str) != EOF) {
        if (flag == 1 && !strchr(str, ':'))
            fprintf(fp2, "%s\n", str);

        if (!strcmp(str, "Value"))
            flag = 1;
    }

    /* Close file descriptors */
    fclose(fp1);
    fclose(fp2);
}

/*
    Converts a non-stationary time series to stationary, by implementing
    the method of differences. Keep in mind that there is a TISEAN routine
    that does pretty much the same thing and probably even better.
*/
void
make_stationary(char *fname) {
    float prev, curr;
    char new_fname[1024];
    FILE *fp1, *fp2;

    strcpy(new_fname, fname);
    strcat(new_fname, ".stat");

    fp1 = fopen(fname, "r");
    fp2 = fopen(new_fname, "w");

    if (!fp1 || !fp2) {
        fprintf(stderr, "An error occured while trying to convert time series to stationary (%s, %s).\n", fname, new_fname);
        exit(0);
    }

    prev = -1;
    while (fscanf(fp1, "%f", &curr) != EOF) {
        if (prev == -1)
            prev = curr;
        else {
            fprintf(fp2, "%f\n", curr - prev);
            prev = curr;
        }
    }

    fclose(fp1);
    fclose(fp2);
}

/*
    Executes a shell command and checks if something went terribly wrong,
    by establishing a one-directional pipe with the TISEAN routines.
*/
void
exec_task(task my_task) {
    char output[MAX_LEN]="";
    FILE *pipe;

    printf("%s\n", my_task.info);
    printf("[%s]\n\n", my_task.cmd);    /* Good for debugging purposes */

    pipe = popen(my_task.cmd, "r");
    fgets(output, MAX_LEN, pipe);
    pclose(pipe);

    printf("%s", output);

    /*
        [Ss]egmentation fault    : Virtually every TISEAN routine
        Wrong type of parameter  : Virtually every TISEAN routine
        Not enough points found! : False nearest neighbors algorithm
        Singular matrix          : Lyapunov spectrum
    */
    if (strstr(output, "fault") || strstr(output, "Wrong") || strstr(output, "Not") || \
	strstr(output, "Singular") || strstr(output, "File")) exit(0);

}