Using sleep portion of rrgen

Author: ascakmak

.gitignore

@@ -2,3 +2,4 @@
 artificialDataAnalysis/rrGenerator/.DS_Store
 .DS_Store
 *.c~
+*.m~

artificialDataAnalysis/pipelineArtificialData.m

@@ -1,4 +1,4 @@
-function pipelineArtificialData(pathToRepo)
+%function pipelineArtificialData(pathToRepo)
 %   
 %   Overview:
 %   	Loads artificially generated data and applies change point
@@ -31,8 +31,8 @@ function pipelineArtificialData(pathToRepo)
 settings.tolerance = 5;     % tolerance in seconds for determining if change point is detected;
                             % if estimated change point (ecp) is within true change point (tcp),
                             % the change point is detected.
-settings.act = 1;           % 1: work with actigraphy data 
-settings.figs = 0;          % 1: plot tcp & ecp on time series
+settings.act = 0;           % 1: work with actigraphy data 
+settings.figs = 1;          % 1: plot tcp & ecp on time series
 settings.iterations = 1000;
 
 % Choose parameter set according to each method
@@ -55,16 +55,34 @@ function pipelineArtificialData(pathToRepo)
 
     % Create random seed
     seed = rand(1)*1e3;
+    
+    % Generate artificial data
+    if settings.act ~= 1
+        % 24 hour RR time series
+        vector_length = 86400;
 
-    % Generate a random vector length between 1e3 and 9e3
-    vector_length = 1000;
-	
-    % Create RR time series with 'rrgen_sys' exe (compiled from c file)
-    [data, tcp] = rrgenV3_wrapper(seed, vector_length, 0, 0, pathToRepo);
+        % Create RR time series with 'rrgen_sys' exe (compiled from c file)
+        [data, tcp, sleepStart, sleepEnd] = rrgenV3_wrapper(seed, vector_length, 0, 0, pathToRepo);
+        
+        % Calculate sleep length
+        sleepLength = (sleepEnd - sleepStart) / 3600;
+        
+        hold off;
+        plot(data); hold on; plot(sleepStart,data(sleepStart),'ro'); hold on; plot(sleepEnd,data(sleepEnd),'ro');
 
-    % Can also specify probability of ectopy and noise: rrgen_2003(seed, vector_length, prob_ectopy, prob_noise);
-    
-    time = cumsum(data); % Time of RR time series is found by cumulative summation
+        % Extract sleep portion of data 
+        data = data(sleepStart:sleepEnd);
+                
+        % Adjust changepoints for sleep portion
+        tcp_idx = find((tcp > sleepStart) & (tcp < sleepEnd));
+        tcp = tcp(tcp_idx) - sleepStart;
+        
+        time = cumsum(data); % Time of RR time series is found by cumulative summation
+    else
+        vector_length = 24000;
+        [data, tcp, sleepStart, sleepEnd] = rrgenV3_wrapper(seed, vector_length, 0, 0, pathToRepo);
+        time = cumsum(data); 
+    end
 
     % Converting rrgen data to actigraphy data
     if settings.act == 1
@@ -122,11 +140,7 @@ function pipelineArtificialData(pathToRepo)
 
     %% Modified Bayesian Online Changepoint Detection
 
-    rl = mbocd(data, lambda, 'gamma');
-    
-    % Locate estimated changepoints
-    [~, ecp] = findpeaks(rl(:,1));
-    ecp(end) = []; % last index is errenous
+    ecp = msegWin(data, lambda, 30000);
 
     % Plot true and estimated change points
     if settings.figs == 1    
@@ -141,9 +155,7 @@ function pipelineArtificialData(pathToRepo)
 
     %% BOCD Original Code
 
-    [~, maxes] = bocd(data, lambda2); 
-    
-    [~, ecp] = findpeaks(maxes(:,1));
+    ecp = segWin(data,lambda2, 30000);
 
     % Plot true and estimated change points
     if settings.figs == 1    
@@ -367,4 +379,4 @@ function pipelineArtificialData(pathToRepo)
 
 tightfig;
 
-end
+%end

artificialDataAnalysis/rrGenerator/rrgenV3_wrapper.m

@@ -1,4 +1,5 @@
- function [rr, tcp] = rrgenV3_wrapper(seed, vector_length, prob_ectopy, prob_noise, pathToRepo)
+ function [rr, tcp, sleepStart, sleepEnd] = rrgenV3_wrapper(seed, vector_length, ...
+                                        prob_ectopy, prob_noise, pathToRepo)
 %   
 %   OVERVIEW:
 %       Matlab wrapper for calling compiled rrgen executable
@@ -17,6 +18,8 @@
 %   OUTPUT
 %       rr - RR interval time series
 %       tcp - True change points in time series
+%       sleepStart - Index of sleep start in artificial RR interval code
+%       sleepEnd - Index of sleep end in artificial RR interval code
 %
 %   REPO
 %       https://github.com/cliffordlab/rrgen
@@ -47,24 +50,37 @@
     vector_length = 200;
 end
 
-command = [pathToRepo filesep 'artificialDataAnalysis' filesep 'rrGenerator' filesep 'rrgenV3' ' ' num2str(seed) ' ' num2str(vector_length) ' ' num2str(prob_ectopy) ' ' num2str(prob_noise)];
+generate = 1;
+while generate == 1
+    command = [pathToRepo filesep 'artificialDataAnalysis' filesep 'rrGenerator' filesep 'rrgenV3' ' ' num2str(seed) ' ' num2str(vector_length) ' ' num2str(prob_ectopy) ' ' num2str(prob_noise)];
 
-% Call rrgen2 via system and save output
-[~, systemout] = system(command);
+    % Call rrgen2 via system and save output
+    [~, systemout] = system(command);
 
-% Read the system output into cells
-systemout_cell = textscan(systemout,'%f %f %f','Delimiter',',');
+    % Read the system output into cells
+    systemout_cell = textscan(systemout,'%f %f %f','Delimiter',',');
 
-% Convert cells of doubles into matrix
-rrgen_mat = cell2mat(systemout_cell);
+    % Convert cells of doubles into matrix
+    rrgen_mat = cell2mat(systemout_cell);
 
-% Isolate columns of matrix into output vectors
-trpeaks = rrgen_mat(rrgen_mat(:,3)==7,1);
-rr = rrgen_mat(rrgen_mat(:,3)~=7,1);
+    % Isolate columns of matrix into output vectors
+    trpeaks = rrgen_mat(rrgen_mat(:,3)==7,1);
+    rr = rrgen_mat(rrgen_mat(:,3)~=7,1);
 
-% Output true change points depending on noise probability
-tcp_idx = find(rrgen_mat(rrgen_mat(:,3)~=7,2)==1);
+    % Output true change points depending on noise probability
+    tcp_idx = find(rrgen_mat(rrgen_mat(:,3)~=7,2)==1);
 
+    % Find sleep start and end indexes
+    sleep_idx = find(rrgen_mat(rrgen_mat(:,3)~=7,3)==0);
+
+    if isempty(sleep_idx)
+        disp('Regenerating RR intervals...');
+    else 
+        generate = 0; % Artificial data generation complete;
+    end
+end
+
+% Adjust for noise if necessary
 tcp = zeros(length(tcp_idx),1);
 if prob_noise==0
     tcp = tcp_idx;
@@ -77,4 +93,19 @@
     end
 end
 
+sleep = zeros(length(sleep_idx),1);
+if prob_noise==0
+    sleep = sleep_idx;
+else
+    sleep_time = trpeaks(sleep_idx);
+    t = cumsum(rr);
+    for i=1:length(sleep)
+        [~,I] = min(abs(t-sleep_time(i)));
+        sleep(i) = I(1);
+    end
+end
+
+sleepStart = sleep(1);
+sleepEnd = sleep(end-1);
+
 end 

readme.md

@@ -68,7 +68,7 @@ Add where R is located to your MATLAB environment <br />
 setenv('PATH', [getenv('PATH'),':','/usr/local/bin']); 
 
 + In artificial data analysis, no data is generated
-Make sure your rrgen_2003_v2 executable is suitable for your computer type. 
+Make sure your rrgenV3 executable is suitable for your computer type. 
 If not, delete this and re compile as decribed above.
 
 ## Useful Links