crqa.m

  set(h8, 'HandleVis','CallBack')
  tags={'crqa_Fig';'axes_logo';'text_logo';'crqa_theiler';'frame';'text';'crqa_axes_Data';'crqa_axes_Var';'crqa_axes_CoVar';'crqa_axes_RR';'crqa_axes_DET';'crqa_axes_L';'crqa_axes_ENTR';'crqa_axes_LAM';'crqa_axes_TT';'crqa_axes_RTE';'crqa_axes_T2';'crqa_m';'crqa_maxLag';'crqa_method';'crqa_eps';'crqa_lmin';'crqa_vmin';'crqa_w';'crqa_ws';'crqa_button_store';'crqa_button_print';'crqa_button_close';'crqa_button_apply'};
  h=[];
  for i=1:length(tags); h=[h; findobj('Tag',tags{i})]; end
  set(h,'Units','Norm')
  if nargout; xout=h8; end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close windows 

case 'close'
  errcode=101;
  set(0,props.root)
  h=findobj('Tag','crqa_Fig');
  if ~isempty(h), close(h(1)), end 
  %clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% store

case 'store'
  errcode=0;
  if ~isempty(findobj('Tag','crqa_button_store'))
    h=findobj('Tag','crqa_button_store');
    h1=findobj('Tag','crqa_button_close');
    if ~isempty(h1), vname_old=get(h1(1),'UserData'); else vname_old=''; end
    if isempty(vname_old), vname_old=''; end
    vname=char(inputdlg('Choose a variable name.','Store output',1,{vname_old}));
    if isempty(vname)
      return
    else
      crqa_values=get(h(1),'UserData');
      assignin('base',vname, [crqa_values])
      warndlg(['CRQA measures have been assigned to the workspace variable ''',vname,'''.'],'Store output');
           waitforbuttonpress
      set(h1(1),'UserData',vname)
    end
  end
  set(0,props.root)

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print
  case 'print'

    errcode=91;
    h=findobj('Tag','crqa_axes_Data','Parent',gcf); h_axes.h(1)=h(1);
    h=findobj('Tag','crqa_axes_Var','Parent',gcf); h_axes.h(2)=h(1);
    h=findobj('Tag','crqa_axes_CoVar','Parent',gcf); if ~isempty(h), h_axes.h(11)=h(1); end
    h=findobj('Tag','crqa_axes_RR','Parent',gcf); h_axes.h(3)=h(1);
    h=findobj('Tag','crqa_axes_DET','Parent',gcf); h_axes.h(4)=h(1);
    h=findobj('Tag','crqa_axes_L','Parent',gcf); h_axes.h(5)=h(1);
    h=findobj('Tag','crqa_axes_ENTR','Parent',gcf); h_axes.h(6)=h(1);
    h=findobj('Tag','crqa_axes_LAM','Parent',gcf); h_axes.h(7)=h(1);
    h=findobj('Tag','crqa_axes_TT','Parent',gcf); h_axes.h(8)=h(1);
    h=findobj('Tag','crqa_axes_RTE','Parent',gcf); h_axes.h(9)=h(1);
    h=findobj('Tag','crqa_axes_T2','Parent',gcf); h_axes.h(10)=h(1);
    h=findobj('Tag','uniLogo');
    tags={'text_logo';'frame';'text';'crqa_m';'crqa_maxLag';'crqa_method';'crqa_eps';'crqa_lmin';'crqa_vmin';'crqa_theiler';'crqa_w';'crqa_ws';'crqa_button_store';'crqa_button_print';'crqa_button_close';'crqa_button_apply';};
    for i=1:length(tags); h=[h; findobj('Tag',tags{i},'Parent',gcf)]; end
    set(h,'Visible','Off')
    
    set(h_axes.h,'Units','Character');
    h_axes.old_pos=get(h_axes.h,'Position');

    axes_height = .13;
    axes_base = 0.065;
    axes_hoffset = .06;
    
    for i=2:2:10
        set(h_axes.h(i-1),  'Units','normalize','Position',[0.1300    axes_base+(5-i/2)*(axes_height+axes_hoffset)    0.3270    axes_height])
        set(h_axes.h(i),'Units','normalize','Position',[0.5780    axes_base+(5-i/2)*(axes_height+axes_hoffset)    0.3270    axes_height])
    end
    if length(h_axes.h) > 10
        set(h_axes.h(11),  'Units','normalize','Position',[0.5780    axes_base+(5-2/2)*(axes_height+axes_hoffset)    0.3270    axes_height])
    end
    h_dlg=printdlg;
    waitfor(h_dlg)

    for i=1:10, set(h_axes.h(i),'Units','Character','Position',h_axes.old_pos{i}), set(h_axes.h(i),'Units','Norm'),end
    if length(h_axes.h) > 10
        set(h_axes.h(11),'Units','Character','Position',h_axes.old_pos{11}), set(h_axes.h(11),'Units','Norm')
    end
    set(h,'Visible','On')
    set(0,props.root)







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute
case 'compute'

  errcode=11;
  if length(method)>1 & strcmpi(method(1:2),'di')
      disp('Warning: RQA from distance plot not possible!')
      return
  end
  
  if ~nogui
    h_fig=findobj('tag','crqa_Fig');
    setptr(gcf,'watch'), 
    obj=({'text';'crqa_m';'crqa_maxLag';'crqa_method';'crqa_eps';'crqa_lmin';'crqa_vmin';'crqa_theiler';'crqa_w';'crqa_ws';'crqa_button_store';'crqa_button_print';'crqa_button_close'});
h_button_apply=findobj('tag','crqa_button_apply','Parent',h_fig(1));

    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','Off')
        end
    end
    set(h_button_apply(1),'ToolTip','Stops the computation.','String','Stop','Callback','set(0,''ShowHidden'',''on'');h=findobj(''tag'',''crqa_button_apply'');set(h(1),''String'',''Stopped'');set(0,''ShowHidden'',''off'')')
  end

if Nx==w & wstep<2, wstep=1; Nx=w; end
if Nx==w, Nx=w; end
if nogui~=2 
    hw=waitbar(0,['0/',num2str(Nx-w)]);
    set(hw,'Name','Please Wait!');
    drawnow 
end
if strcmpi(method,'Order Pattern') method = 'op'; end
if strcmpi(method,'Order Matrx') method = 'om'; end
if strcmpi(method,'Maximum Norm, fixed RR') method = 'rr'; end


% check if plugin exist and is executable
[plugin_exist, plugin_name, plugin_path] = is_crp_plugin;
if nogui == 1 & plugin_exist & ( method_n < 4 | method_n == 8 ) & length(x) == length(y)
    disp('(plugin used)')
end

errcode=20;

% general histograms of line structures for significance test
hist_l = []; hist_v = []; hist_w = [];


% main loop through all windows
Y = zeros(length(xscale)-w+1,13);
if undocumented
   Y = zeros(length(xscale)-w+1,17);
end

if ~nogui
    set(0,'ShowHidden','on')
end
   
for i=1:wstep:Nx-w+1; 
     if ~nogui
         if strcmpi(get(h_button_apply(1),'string'),'stopped')
              Y(i:Nx-w,1:13)=NaN;
              drawnow
              break
         end
     end


     if time_scale_flag
        x_var=var(x(i:i+w-1,2)); y_var=var(y(i:i+w-1,2)); 
        temp=cov(x(i:i+w-1,2),y(i:i+w-1,2));
     else
        x_var=var(x(i:i+w-1,1)); y_var=var(y(i:i+w-1,1)); 
        temp=cov(x(i:i+w-1,1),y(i:i+w-1,1));
     end
     xy_var=temp(1,2);

          
     do_norm = {'non';'nor'};

     % if plugin exist and method is MAX, MIN or EUC
     if plugin_exist & ( method_n < 4 | method_n == 9 ) & length(x) == length(y) 

         errcode=21;
         warning off
         while 1
             tmp_xdatafile = tempname;
             if ~exist(tmp_xdatafile), break, end
         end
         while 1
             tmp_ydatafile = tempname;
             if ~exist(tmp_ydatafile), break, end
         end
         while 1
             tmp_rqadatafile = tempname;
             if ~exist(tmp_rqadatafile), break, end
         end
         
         
         x_tmp = x(i:i+w-1,:);
         y_tmp = y(i:i+w-1,:);
         
         % save data in temporary file
         save(tmp_xdatafile,'x_tmp','-ascii','-tabs');
         save(tmp_ydatafile,'y_tmp','-ascii','-tabs');

         % call extern rp programme
         m_str = {'MAX', 'EUC', 'MIN', 'NR', 'RR', 'FAN', 'IN', 'OM', 'OP', 'EUC'};

         [status ] = system([plugin_path,filesep,plugin_name,' -m ',num2str(m), ...
                                       ' -t ',num2str(t), ...
                                       ' -e ',num2str(e), ...
                                       ' -n ',m_str{method_n}, ...
                                       ' -w ',num2str(theiler_window), ...
                                       ' -l ',num2str(lmin), ...
                                       ' -v ',num2str(vmin), ...
                                       ' -i ',tmp_xdatafile, ...
                                       ' -j ',tmp_ydatafile, ...
                                       ' -o ',tmp_rqadatafile, ...
                                       ' -s']);
         errcode=22;
         % import RQA
         rqa_in = [];
         try
             fid = fopen(tmp_rqadatafile,'r'); % open RQA data file
             while 1
                 dataStr = fgetl(fid); % read RQA data
                 if ~ischar(dataStr), break, end % leave the loop if end of file
                 if isempty(findstr(dataStr, '#')) % neglect comments line (e.g. header)
                     rqa_in = str2num(dataStr); % import RQA measures into local variable rqa
                 end
             end
             fclose(fid); % close RQA data file

             RR = rqa_in(1);
             DET = rqa_in(2);
             LAM = rqa_in(4);
             Lmax = rqa_in(6);
             L = rqa_in(7);
             ENTR = rqa_in(8);
             Vmax = rqa_in(10);
             TT = rqa_in(11);
             t1 = rqa_in(13);
             t2 = rqa_in(14);
             RT = rqa_in(16);
             RTmax = rqa_in(15);
             RF = rqa_in(19);
             ENTW = rqa_in(17);
             warning off
             delete(tmp_rqadatafile);
             delete(tmp_xdatafile);
             delete(tmp_ydatafile);
             warning on
         catch 
             warning off
             delete(tmp_rqadatafile);
             delete(tmp_xdatafile);
             delete(tmp_ydatafile);
             warning on
         end

         if nogui~=2
              waitbar(i/(Nx-w),hw,[num2str(i),'/',num2str(Nx-w)]); 
              drawnow 
         end



     % use builtin implementation
     else
         errcode=25;
         try
             if ~flag_pdist
                 if time_scale_flag
                   if length(x(i:i+w-1,:)) > 2000
                      X=crp_big(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
                   else
                      X=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
                   end
                 else
                   X=crp2(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
                   if isnan(X), break; end
                 end
           %  X=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,varargin{i_char},'silent');
           else
               %%%%%%%%%%%%
               % alternative using pdist

               xcor = @(x,y) sqrt(sum((repmat(x,size(y,1),1)-y).^2,2));

               x_dist = pdist(x(i:i+w-1,:),xcor);
               X = squareform(x_dist);
           end
           
           %%%%%%%%%%%%
           
           warning off 
           if nogui~=2 
               waitbar(i/(Nx-w),hw,[num2str(i),'/',num2str(Nx-w)]); 
           end

         %if 0
         catch
         lasterr
             error(lasterr)
             if nogui~=2 & ishandle(hw), close(hw), end
         end

         N=size(X);
         if theiler_window > 0
             X_theiler=double(triu(X,theiler_window))+double(tril(X,-theiler_window));
         else
             X_theiler=X;
         end
         errcode=26;

         % compute recurrence times of 1st and 2nd type
         if size(X_theiler,2) > 1000
             t1=[];t2=[];
             rps2=find(diff(double(X_theiler(:)))==1);
             rps=find(X_theiler(:));
             t1=diff(rps);
             t2=diff(rps2);
         else
             t1 = []; t2 = [];
             for i2=1:size(X_theiler,2)
                 if nogui~=2 & (Nx-w < 2) & ~rem(i2,25), waitbar(i2/size(X_theiler,2)), end
                 rps2=find(diff(double(X_theiler(:,i2)))==1);
                 rps=find(X_theiler(:,i2));
                 t1=[t1;diff(rps)];
                 t2=[t2;diff(rps2)];
             end
         end
         t1=mean(t1);
         t2=mean(t2);

         errcode=27;
         [a b]=dl(X_theiler);
         hist_l = [hist_l; b];
         N_hist_l(i) = length(b);

         warning off 
         errcode=271;
         b(b<lmin)=[];
         [c d dw]=tt(X_theiler);
         hist_v = [hist_v; d];
         hist_w = [hist_w; dw];         
         N_hist_v(i) = length(d);
         N_hist_w(i) = length(dw);
         
         warning off 
         errcode=272;
         d(find(d<vmin))=[];

         errcode=273;
         N_all = (N(1)*N(2));
         % reduce the number of possible states by the Theiler window;
         % because the Theiler window is applied symmetrically (on the LOI)
         % we only use N(1) and not N(2)
         if theiler_window >= 1
             N_all = N_all - N(1) - 2*((theiler_window-1)*N(1) - sum(1:(theiler_window-1)));
         end
         
         RR=sum(X_theiler(:))/N_all;
         %b(find(b>=max(N)-lmin))=[]; if isempty(b), b=0; end
         if isempty(b), b=0; end
         errcode=274;
         if sum(X_theiler(:)) > 0
           DET=sum(b)/sum(X_theiler(:));
         else
           DET=NaN;
         end
         errcode=275;
         L=mean(b);
         histL=hist(b(:),[1:min(N)]);
         ENTR=entropy(histL(:));
         errcode=276;
         if sum(X_theiler(:))>0
           LAM=sum(d)/sum(sum(X_theiler));
         else
           LAM=NaN;
         end
       
         % recurrence times
         RT = mean(dw);
         bins = [.5:max(dw)+.5];
         if ~isnan(dw) % check for empty dw
             dwh = histc(dw,bins); dwh = dwh(:);
             %[dwh dwi] = hist(dw,[1:max(dw)]');
             if dwh
                 [dws dwsi] = sort(dwh);
                 %RTp = dwi(dwsi(end)); % most probable recurrence time
                 RTp = bins(dwsi(end))+.5; % most probable recurrence time
             else
                 RTp = 0;
             end
             RTmax = max(dw); % maximal recurrence time
             RF = 1/RTmax; % minimal recurrence frequency
             ENTW = entropy(dwh(:));
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             %% new measures
             % RPDE
             errcode=278;
             if max(dw) > 0 && numel(bins) > 3
                RTE = entropy(dwh) / log(max(dw));
             else
                RTE = NaN;
             end
         else
             RTmax = NaN;
             RF = NaN;
             ENTW = NaN;
             RTE = NaN;
         end

         errcode=277;
         TT=mean(d);
         b=[b;0]; Lmax=max(b);
         d=[d;0]; Vmax=max(d);
         
         
        
        
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         %% network measures
         % global clustering coefficient
         errcode=279;
         kv = sum(X_theiler,1); % degree of nodes
         Clust = mean(diag(double(X_theiler)*double(X_theiler)*double(X_theiler))' ./ (kv .* (kv-1)));
         % transitivity
         denom = sum(sum(double(X_theiler) * double(X_theiler)));
         Trans = trace(double(X_theiler)*double(X_theiler)*double(X_theiler))/denom;
         
     end % end plugin
  
     warning on

     errcode=28;
     Y(i,1)=RR; 
     Y(i,2)=DET;
     Y(i,3)=L;
     Y(i,4)=Lmax;
     Y(i,5)=ENTR;
     Y(i,6)=LAM;
     Y(i,7)=TT;
     Y(i,8)=Vmax;
     Y(i,9)=t1;
     Y(i,10)=t2;
     Y(i,11)=RTE;
     Y(i,12)=Clust;
     Y(i,13)=Trans;
     if undocumented
         Y(i,14)=RT;
         Y(i,15)=RTmax;
         Y(i,16)=RF;
         Y(i,17)=ENTW;
     else
         Y(i,15)=x_var;
         Y(i,16)=y_var;
         Y(i,17)=xy_var;
     end
 end % end window loop
 
% significance by bootstrap
 
% nboot = 100; 
% nD = round(mean(N_hist_l(N_hist_l>0)));
% nV = round(mean(N_hist_v(N_hist_v>0)));
% nW = round(mean(N_hist_w(N_hist_w>0)));
% for i = 1:nboot
%       onesample = ceil(length(hist_l)*rand(length(hist_l),1));
%       onesample = onesample(ceil(nD*rand(nD,1)));
%       tmp = sum(hist_l(onesample,:) >= lmin )/ sum(hist_l(onesample,:));
%       bootstatDET(i,:) = (tmp(:))';
%       tmp = mean(hist_l(onesample,:) >= lmin);
%       bootstatL(i,:) = (tmp(:))';
% 
%       onesample = ceil(length(hist_v)*rand(length(hist_v),1));
%       onesample = onesample(ceil(nV*rand(nV,1)));
%       tmp = sum(hist_v(onesample,:) >= vmin )/ sum(hist_v(onesample,:));
%       bootstatLAM(i,:) = (tmp(:))';
%       tmp = mean(hist_v(onesample,:) >= vmin);
%       bootstatTT(i,:) = (tmp(:))';
% 
%       onesample = ceil(length(hist_w)*rand(length(hist_w),1));
%       onesample = onesample(ceil(nW*rand(nW,1)));
%       tmp = mean(hist_w(onesample,:));
%       bootstatT2(i,:) = (tmp(:))';
% end

 
if ishandle(hw), waitbar(1); drawnow; close(hw), end % close waitbar

  if ~nogui
    set(h_button_apply(1),'ToolTip','Starts the computation.','String','Apply','Callback','crqa compute')
    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','On')
        end
    end
  end

  if ~nogui, 
    h=findobj('Tag','crqa_Fig'); 
    if ~isempty(h), set(0,'CurrentFigure',h(1)); end
    setptr(gcf,'arrow'), 
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot
if ~nogui
    set(0,'showhidden','on')
    errcode=30;
    h=findobj('Tag','crqa_Fig'); if ~isempty(h), set(0,'CurrentFigure',h(1)); end
    tx={'RR';'DET';'L';'ENTR';'LAM';'TT';'RTE';'T_2';'Variance';'Covariance'};
    index=[1,2,3,5,6,7,11,10,15,16,17];
    tags={'crqa_axes_RR','crqa_axes_DET','crqa_axes_L','crqa_axes_ENTR','crqa_axes_LAM','crqa_axes_TT','crqa_axes_RTE','crqa_axes_T2','crqa_axes_Var','crqa_axes_CoVar'};
    h=findobj('Tag','crqa_axes_RR','Parent',gcf); h_axes.h(1)=h(1);
    h=findobj('Tag','crqa_axes_DET','Parent',gcf); h_axes.h(2)=h(1);
    h=findobj('Tag','crqa_axes_L','Parent',gcf); h_axes.h(3)=h(1);
    h=findobj('Tag','crqa_axes_ENTR','Parent',gcf); h_axes.h(4)=h(1);
    h=findobj('Tag','crqa_axes_LAM','Parent',gcf); h_axes.h(5)=h(1);
    h=findobj('Tag','crqa_axes_TT','Parent',gcf); h_axes.h(6)=h(1);
    h=findobj('Tag','crqa_axes_RTE','Parent',gcf); h_axes.h(7)=h(1);
    h=findobj('Tag','crqa_axes_T2','Parent',gcf); h_axes.h(8)=h(1);
    h=findobj('Tag','crqa_axes_Var','Parent',gcf); h_axes.h(9)=h(1);
    if ~all(x(:)==y(:)), h=findobj('Tag','crqa_axes_CoVar','Parent',gcf); h_axes.h(10)=h(1); end
   for i=1:9,
     set(gcf,'CurrentAxes',h_axes.h(i))
     if size(Y,1)==1
       cla
       text(0.5,0.5,sprintf('%6.4f',Y(index(i))),'FontWeight','bold','HorizontalAlign','Center')
%        if ~all(x(:)==y(:)) & i==9
%          cla 
%              text(0.5,0.5,sprintf('%6.5f, %6.5f',Y(index(i)),Y(index(i+1))),'FontWeight','bold','HorizontalAlign','Center')
%              set(gcf,'CurrentAxes',h_axes.h(i+1));cla
%          set(h_axes.h(i+1),'visible','off');
%        end
     else
       if i==9
         cla
             h2=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(i)));
             set(h2,'color',props.line.Color)
         if ~all(x(:)==y(:))
           set(gca,'Color','none')
           hold on; h2=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(10)),'r');
           set(gcf,'CurrentAxes',h_axes.h(i+1));cla
           set(h_axes.h(10),'visible','on');
           
           h3=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(11))); set(h3,'color',[0 .4 0]);
           ylabel(tx(i+1));
           set(h_axes.h(i+1),'Tag',tags{i+1},'Units','Norm','Color','none','YAxisLocation','right','YColor',[0 .4 0])
         end
       else
         plot(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(i)),'color',props.line.Color)
       end
     end
     set(gcf,'CurrentAxes',h_axes.h(i));
     ylabel(tx(i));
     set(gca,'Tag',tags{i},'color',props.axes.Color,'Units','Norm',...
        'xlim',[xscale(1) xscale(end)])
     if length(Y(:,1)) == 1, axis([0 1 0 1]), end
   end
   h=findobj('Tag','crqa_button_store');
   set(h(1),'UserData',Y,'Enable','On')
   xlim = get(h_axes.h(9),'xlim');
   if length(Y(:,1)) > 1, set(h_crqa_axes_Data,'xlim',xlim), end
   if ~all(x(:)==y(:)) & length(Y(:,1)) == 1, set(h_axes.h(10),'visible','off'); delete(get(h_axes.h(10),'children')), end
else
  if nargout, xout=Y(:,argout_index); end
end

if nargout, xout=Y(:,argout_index); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end 

end




%if 0
catch
  if nogui~=2, if ishandle(hw), close(hw), end , end
  if nargout, xout = NaN; end
      
  z=whos;x=lasterr;y=lastwarn;in=varargin{1};
  print_error('crqa',z,x,y,in,method,action)
  try, if ~nogui
    h=findobj('tag','crqa_button_apply');
    set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','crqa compute')
    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','On')
        end
    end, end
  end
  set(0,'showhidden','on')
    h=findobj('Tag','crqa_Fig'); if ~isempty(h), setptr(h(1),'arrow'); end
  set(0,'showhidden','off')
end


try, set(0,props.root), end



try
  set(0,'ShowHidden','off')
end