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 ** ********************************************************************************/voidreconstruct_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 ** ********************************************************************************/voidrun_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 ** ********************************************************************************/voidcalc_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 ** ********************************************************************************/voiddirs_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 ** ********************************************************************************/intmain(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;}