rendered paste body/*
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;
}