perturbRecovery

0001 function [] = perturbRecovery(animals)
0002 
0003 if nargin < 1
0004   animals = 1:9;
0005 end
0006 
0007 % tmts = {'control'};
0008 tmts = {'control','mass','inertia'};
0009 
0010 ylims = {20*[-1,1],30*[-1,1],40*[-1,1],[-5,20],[-4,4],[],[-70,-10],40*[-1,1],[],[],1000*[-1,1],1000*[-1,1],1000*[-1,1]};
0011 
0012 states = {'pitch','roll','yaw','x','y','z','vx','vy','ax','ay','dpitch','droll','dyaw'};
0013 % plot_state = 1-[1,1,1,0,0,0,1,1,0,0,1,1,1];
0014 plot_state = [1,1,1,0,0,0,1,1,0,0,1,1,1];
0015 % plot_state = [0,0,0,0,0,0,0,0,0,0,0,1,1];
0016 % plot_state = [0,1,0,0,0,0,0,0,0,0,0,0,0];
0017 % plot_state = [0,0,0,0,0,0,0,1,0,0,0,0,0];
0018 units = {'deg','deg','deg','cm','cm','cm','cm/sec','cm/sec','g','g','deg/sec','deg/sec','deg/sec'};
0019 
0020 sig_plot = 1;
0021 save_pre = 0;
0022 save_post = 0;
0023 save_postc = 1;
0024 
0025 cols = hsv(length(tmts));
0026 
0027 for t = 1:length(tmts)
0028   load(['data/',tmts{t},'_lp.mat'],'lp');
0029   eval(['lp',tmts{t},' = lp;']);
0030   load(['data/',tmts{t},'_rec.mat'],'rec');
0031   eval(['rec',tmts{t},' = rec;']);
0032 end
0033 
0034 % Phase over one cycle
0035 np = size(rec(1).Ppre,2);
0036 pinc = (2*pi/np);
0037 p = 0:pinc:(2*pi-pinc);
0038 
0039 ani = struct('c',[],'m',[],'i',[]);
0040 
0041 % Loop through animals
0042 for a = animals
0043   
0044   n = [];
0045   ti = '';
0046   leg = {};
0047   
0048   if length(animals) == 1
0049     ti = [ti,', a = ',num2str(animals)];
0050   end
0051   
0052   % Loop through treatments
0053   for t = 1:length(tmts)
0054     
0055     % Get data for treatment
0056     eval(['lp = lp',tmts{t},';']);
0057     eval(['rec = rec',tmts{t},';']);
0058     
0059     Bpr = [];
0060     
0061     % Loop through states
0062     for s = 1:length(states)
0063   
0064       if ~plot_state(s)
0065         continue;
0066       end
0067 
0068       if save_post
0069         postfig = figure(1);
0070         clf
0071         if sig_plot
0072           s1 = subplot(3,2,[1,2,3,4]);
0073           hold on;
0074           s2 = subplot(3,2,[5,6]);
0075           hold on;
0076         else
0077           s1 = gca;
0078           hold on;
0079         end
0080       end
0081       
0082       Ppre = [];
0083       Bpre = [];
0084 
0085       % Loop through trials for animal
0086       for k = find(lp.animal == a)'
0087         
0088         if ~isempty(rec(k).Ppre)
0089           
0090           Ppre = [Ppre, rec(k).Ppre'];
0091           Bpre = [Bpre, rec(k).Bpre(:,:,s)'];
0092         end
0093         
0094       end % Loop through trials for each animal
0095       
0096       Bpr(1:length(Bpre(:)),s) = Bpre(:);
0097 
0098       ord = 5;
0099 
0100       Bo(s) = ftFit( 5, mod(Ppre(:),2*pi), Bpre(:) );
0101       Bm = [];
0102       for k = 1:size(Ppre,2)
0103         Bm(:,k) = ftVal( Bo(s), Ppre(:,k) );
0104       end
0105       Be = Bpre - Bm;
0106       L20 = sqrt(sum(Be.^2,1));
0107       
0108       m0 = mean(L20);
0109       s0 = std(L20);
0110       
0111       z0 = 1.65*s0+m0;
0112       
0113       L2 = [];
0114       
0115 %       keyboard
0116 
0117       cols = hsv(length(find(lp.animal == a)));
0118       
0119       K = find(lp.animal == a)';
0120       
0121       for k = 1:length(K);
0122         
0123         if ~isempty(rec(K(k)).P) && length(rec(K(k)).P)-rec(K(k)).post >= np
0124           
0125           Bp = ftVal( Bo(s), rec(K(k)).P(rec(K(k)).post+(0:np-1)) )';
0126           
0127           for n = 0:length(rec(K(k)).P)-np-rec(K(k)).post-1
0128             
0129             L2(k,n+1) = sqrt(sum((rec(K(k)).B(1,rec(K(k)).post+n+(0:np-1),s) - Bp(mod((0:np-1)+n,np)+1)).^2));
0130             
0131           end
0132           
0133           L2(k,find(L2(k,:)==0)) = nan;
0134       
0135           rs = find(L2(k,:) <= z0,1,'first');
0136           if isempty(rs) || rs == 1
0137             rp = nan;
0138             r = nan;
0139           else
0140             
0141             rp = rs*2*pi/np;
0142             
0143             Bp = ftVal( Bo(s), rec(K(k)).P(rec(K(k)).post + (0:rs)) )';
0144             r = sum(abs(rec(K(k)).B(1,rec(K(k)).post+(0:rs),s)-Bp));
0145           end
0146 %           eval(['ani(a).',tmts{t}(1),'(s,k) = rp;']);
0147           eval(['ani(a).',tmts{t}(1),'(s,k) = r;']);
0148           
0149           if save_post
0150 
0151             % Plot backpack variable
0152             plot(s1,rec(K(k)).P,rec(K(k)).B(:,:,s),'color',cols(k,:),'LineWidth',1);
0153             
0154             % Plot L2 norm
0155             plot(s2,rec(K(k)).P(rec(K(k)).post+(1:sum(~isnan(L2(k,:))))),L2(k,1:sum(~isnan(L2(k,:)))),'color',cols(k,:),'LineWidth',1);
0156 
0157             if ~isempty(ylims{s})
0158               ylim1 = ylims{s};
0159             else
0160               ylim1 = get(s1,'YLim');
0161             end
0162             
0163             ylim2 = get(s2,'YLim');
0164             
0165             % Plot perturbation
0166             plot(s1,rec(K(k)).P(rec(K(k)).post)*[1,1],ylim1,'--','color',cols(k,:));
0167             plot(s2,rec(K(k)).P(rec(K(k)).post)*[1,1],ylim2,'--','color',cols(k,:));
0168           end
0169           
0170         end
0171       end
0172       
0173 %       keyboard
0174         
0175       
0176       if save_post
0177       
0178         figure(postfig);
0179 
0180         if ~isempty(ylims{s})
0181           ylim = ylims{s};
0182         else
0183           ylim = get(s1,'YLim');
0184         end
0185         
0186         xlim = [-4*pi,8*pi];
0187         
0188         Bf = ftVal( Bo(s), xlim(1):pinc:xlim(2) );
0189         plot(s1,xlim(1):pinc:xlim(2),Bf,'k','LineWidth',2);
0190         
0191         plot(s2,xlim,z0*[1,1],'k','LineWidth',2);
0192 
0193         xlabel(s2,'Phase (rad)','FontName','Times','FontSize',18);
0194         ylabel(s1,[states{s},' (',units{s},')'],'FontName','Times','FontSize',18);
0195         ylabel(s2,'{L_2}','FontName','Times','FontSize',18);
0196         title(s1,[tmts{t},': a#',num2str(a),', ',states{s},' (',units{s},')'],'FontName','Times','FontSize',22);
0197         xtick = [fliplr(0:-6.28:xlim(1)),6.28:6.28:xlim(2)];
0198         set(s1,'XLim',xlim,'XTick',xtick,'YLim',ylim,'FontName','Times','FontSize',14);
0199         set(s2,'XLim',xlim,'XTick',xtick,'FontName','Times','FontSize',14);
0200         
0201 %         keyboard
0202 
0203         figprefix = [tmts{t}(1),'_a',num2str(a),'_'];
0204 
0205         dirs = {['fig/rec/'],['fig/rec/fig'],['fig/rec/eps']};
0206         for d = 1:length(dirs)
0207           if ~exist(dirs{d},'dir')
0208             mkdir(dirs{d});
0209           end
0210         end
0211         saveas(gca,['fig/rec/',figprefix,states{s},'.png'],'png');
0212 %         saveas(gca,['fig/fig/',figprefix,states{s},'.fig'],'fig');
0213 %         saveas(gca,['fig/eps/',figprefix,states{s},'.eps'],'eps');
0214       end      
0215       
0216 %       keyboard
0217 
0218       if save_pre
0219       
0220         figure(2);
0221         plot(mod(Ppre,2*pi),Bpre,'.');
0222 
0223 
0224         hold on;
0225         plot(p,Bp,'k','LineWidth',6);
0226 
0227         if ~isempty(ylims{s})
0228           ylim = ylims{s};
0229         else
0230           ylim = get(gca,'YLim');
0231         end
0232 
0233         xlabel('Phase (rad)','FontName','Times','FontSize',18);
0234         ylabel([states{s},' (',units{s},')'],'FontName','Times','FontSize',18);
0235         title([tmts{t},': a#',num2str(a),', ',states{s},' (',units{s},')'],'FontName','Times','FontSize',22);
0236         set(gca,'XLim',[0,2*pi],'YLim',ylim,'FontName','Times','FontSize',14);
0237 
0238         figprefix = [tmts{t}(1),'_a',num2str(a),'_'];
0239 
0240         dirs = {['fig/pre/'],['fig/pre/fig'],['fig/pre/eps']};
0241         for d = 1:length(dirs)
0242           if ~exist(dirs{d},'dir')
0243             mkdir(dirs{d});
0244           end
0245         end
0246         saveas(gca,['fig/',figprefix,states{s},'.png'],'png');
0247 %         saveas(gca,['fig/fig/',figprefix,states{s},'.fig'],'fig');
0248 %         saveas(gca,['fig/eps/',figprefix,states{s},'.eps'],'eps');
0249       end
0250     end % Loop through states
0251     
0252     Ppr = Ppre(:);
0253 
0254     ord = 5;
0255     
0256     % Which states to combine
0257     st = 10+[1,2,3];
0258 
0259     Bo = ftFit( 5, mod(Ppre(:),2*pi), Bpr );
0260     Bm = ftVal( Bo, mod(Ppre(:),2*pi) );
0261     Be = Bpr - Bm;
0262     L20 = sqrt(sum(reshape(sum(Be(:,st).^2,2),np,size(Be,1)/np),1));
0263 
0264     m0 = mean(L20);
0265     s0 = std(L20);
0266 
0267     z0 = 1.65*s0+m0;
0268     
0269     
0270     if save_postc
0271       
0272       figure(137);
0273       clf;
0274       hold on;
0275     end
0276 
0277 
0278     K = find(lp.animal == a)';
0279 
0280     for k = 1:length(K);
0281 
0282       if ~isempty(rec(K(k)).P) && length(rec(K(k)).P)-rec(K(k)).post >= np
0283 
0284         Bp = ftVal( Bo, rec(K(k)).P(rec(K(k)).post+(0:np-1)) );
0285 
0286         for n = 0:length(rec(K(k)).P)-np-rec(K(k)).post-1
0287 
0288           L2(k,n+1) = sqrt(sum(sum((squeeze(rec(K(k)).B(1,rec(K(k)).post+n+(0:np-1),st)) - Bp(mod((0:np-1)+n,np)+1,st)).^2,2),1));
0289 
0290         end
0291 
0292         L2(k,find(L2(k,:)==0)) = nan;
0293 
0294         rs = find(L2(k,:) <= z0,1,'first');
0295         if isempty(rs) || rs == 1
0296           rp = nan;
0297 %           r = nan;
0298         else
0299 
0300           rp = rs*2*pi/np;
0301 
0302 %           Bp = ftVal( Bo, rec(K(k)).P(rec(K(k)).post + (0:rs)) )';
0303 %           r = sum(abs(rec(K(k)).B(1,rec(K(k)).post+(0:rs),s)-Bp));
0304         end
0305         eval(['ani(a).',tmts{t}(1),'(14,k) = rp;']);
0306 %         eval(['ani(a).',tmts{t}(1),'(s,k) = r;']);
0307 
0308         if save_postc
0309           % Plot L2 norm
0310           plot(rec(K(k)).P(rec(K(k)).post+(1:sum(~isnan(L2(k,:))))),L2(k,1:sum(~isnan(L2(k,:)))),'color',cols(k,:),'LineWidth',1);
0311 
0312         end
0313 
0314       end
0315     end
0316     
0317     if save_postc
0318 
0319       xlabel('Phase (rad)','FontName','Times','FontSize',18);
0320       ylabel('{L_2}','FontName','Times','FontSize',18);
0321       title([tmts{t},': a#',num2str(a)],'FontName','Times','FontSize',22);
0322       set(gca,'FontName','Times','FontSize',14);
0323 
0324       %         keyboard
0325 
0326       figprefix = [tmts{t}(1),'_a',num2str(a),'_c'];
0327 
0328       dirs = {['fig/recc/'],['fig/recc/fig'],['fig/recc/eps']};
0329       for d = 1:length(dirs)
0330         if ~exist(dirs{d},'dir')
0331           mkdir(dirs{d});
0332         end
0333       end
0334       saveas(gca,['fig/rec/',figprefix,'.png'],'png');
0335     end
0336     
0337     keyboard
0338 
0339   end % Loop through treatments
0340 %
0341 %   if sig_plot
0342 %     % Compare Mass and Inertia to Control using Mann-Whitney U statistic
0343 %     um = [];
0344 %     pm = [];
0345 %     uj = [];
0346 %     pj = [];
0347 %     for k = 1:length(P)
0348 %       cm = [ones(size(X{1},1),1);2*ones(size(X{2},1),1)];
0349 %       [um(k),pm(k)] = MWU(cm,[X{1}(:,k);X{2}(:,k)]);
0350 %       cj = [ones(size(X{1},1),1);2*ones(size(X{3},1),1)];
0351 %       [uj(k),pj(k)] = MWU(cj,[X{1}(:,k);X{3}(:,k)]);
0352 %     end
0353 %
0354 %     pmf = filtfilt(B,A,pm);
0355 %     pjf = filtfilt(B,A,pj);
0356 %     plot(s2,P,pmf,'LineWidth',2,'color',cols(2,:));
0357 %     plot(s2,P,pjf,'LineWidth',2,'color',cols(3,:));
0358 %
0359 %     alpha = 0.05;
0360 %     plot(s2,P([1,end]),alpha*[1,1],'k','LineWidth',2);
0361 %
0362 %
0363 %     % Print recovery times
0364 % %     disp([states{s},' recovery times (msec)']);
0365 % %     pt = find(TT(sa)>0,1,'first');
0366 % %     disp(['mass : ',num2str((find(diff(pmf(pt:end) > alpha)==1))/0.5)]);
0367 % %     disp(['inertia : ',num2str((find(diff(pif(pt:end) > alpha)==1))/0.5)]);
0368 %
0369 % %     keyboard
0370 %   end
0371 %
0372 %
0373 %   grid on;
0374 %   ylabel(s1,[states{s},' (',units{s},')'],'FontSize',18,'FontName','Times');
0375 %   title(s1,[states{s},' (',units{s},')',ti],'FontSize',18,'FontName','Times');
0376 %   legend(s1,leg);
0377 %   set(s1,'XLim',rng,'FontSize',14,'FontName','Times');
0378 %
0379 %   if sig_plot
0380 %     xlabel(s2,'Time (sec)','FontSize',18,'FontName','Times');
0381 %     ylabel(s2,'p-value','FontSize',18,'FontName','Times');
0382 %     legend(s2,{'mass','inertia'});
0383 %     set(s2,'XLim',rng,'YLim',[0,0.55],'FontSize',14,'FontName','Times');
0384 %   end
0385 %
0386 %
0387 %   if save_plots
0388 %     if length(animals) == 1
0389 %       figprefix = ['a',num2str(animals),'_'];
0390 %     else
0391 %       figprefix = '';
0392 %     end
0393 %
0394 %     figprefix = [figprefix,'tmtvp_'];
0395 %
0396 %     dirs = {['fig/'],['fig/fig'],['fig/eps']};
0397 %     for d = 1:length(dirs)
0398 %       if ~exist(dirs{d},'dir')
0399 %         mkdir(dirs{d});
0400 %       end
0401 %     end
0402 %     saveas(gca,['fig/fig/',figprefix,lp.states{s},'.fig'],'fig');
0403 %     saveas(gca,['fig/',figprefix,lp.states{s},'.png'],'png');
0404 %     saveas(gca,['fig/eps/',figprefix,lp.states{s},'.eps'],'eps');
0405 %   end
0406 end % Loop through animals
0407 
0408 save(['data/ani.mat'],'ani');