crqa.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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'});
    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','Off')
        end
    end
    h=findobj('tag','crqa_button_apply');
    set(h(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+1; end
if Nx==w, Nx=w+1; end
if nogui~=2, hw=waitbar(0,['0/',num2str(Nx-w)]);set(hw,'Name','Please Wait!');h1=get(hw,'chil');h1=get(h1,'title'); 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 = [];

for i=1:wstep:Nx-w; 
     if ~nogui
         set(0,'ShowHidden','on')
         h=findobj('tag','crqa_button_apply','Parent',h_fig(1));
         if strcmpi(get(h(1),'string'),'stopped')
              Y(i:Nx-w,1:6)=NaN;
              break
         end
     end


     if time_scale_flag
        x_var=var(x(i:i+w-1,2:end)); y_var=var(y(i:i+w-1,2:end)); 
        temp=cov(x(i:i+w-1,2:end),y(i:i+w-1,2:end));
     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 & ishandle(h1), set(h1,'str',[num2str(i),'/',num2str(Nx-w)]); waitbar(i/(Nx-w)); drawnow, end



     % use builtin implementation
     else
     
         errcode=25;
         try
         %  X=crp_big(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,'fan','silent');
             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');
             end
           %  X=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,varargin{i_char},'silent');
           warning off 
           if nogui~=2 & ishandle(h1), set(h1,'str',[num2str(i),'/',num2str(Nx-w)]); waitbar(i/(Nx-w)); drawnow, end

         %if 0
         catch
             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);

         [dwh dwi] = hist(dw,[1:max(dw)]);
         if dwh
             [dws dwsi] = sort(dwh);
             RTp = dwi(dwsi(end)); % most probable recurrence time
         else
             RTp = 0;
         end
         RTmax = max(dw); % maximal recurrence time
         RF = 1/RTmax; % minimal recurrence frequency
         ENTW = entropy(dwh(:));

         errcode=277;
         TT=mean(d);
         b=[b;0]; Lmax=max(b);
         d=[d;0]; Vmax=max(d);
         
     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;
     if undocumented
         Y(i,11)=RT;
         Y(i,12)=RTmax;
         Y(i,13)=RF;
         Y(i,14)=ENTW;
     end
     Y(i,15)=x_var;
     Y(i,16)=y_var;
     Y(i,17)=xy_var;
 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
    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

  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';'T_1';'T_2';'Variance';'Covariance'};
    index=[1,2,3,5,6,7,9,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_T1','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_T1','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')
   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
      
  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