Home > polypedal > cartperturb > original_shai.m

original_shai

PURPOSE ^

% Only need to do this once per analysis session

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

% Only need to do this once per analysis session

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 
0002 %% Only need to do this once per analysis session
0003 load original/most.mat
0004 sam = 0;
0005 do_xr
0006 do_xrg
0007 
0008 %% Identify animals
0009 f = fopen('roachruns-fix.csv','rt');
0010 ani=zeros(length(xr),0);
0011 while ~feof(f);
0012     s = fscanf(f,'%s,');
0013     trial = s(2:end-3);
0014     animal = s(end);
0015     for k=1:length(xr)
0016         if strcmp(xr{k}.filePath,trial)
0017             xr{k}.animal=animal;
0018             ani(k,str2num(animal))=1;
0019             fprintf('[%d] %s is %s\n',k,trial,animal);
0020         end
0021     end
0022 end
0023 fclose(f)
0024 
0025 %% Phase separation using hybrid Shai/Sam code
0026 
0027 if 0
0028   close all;
0029   sam = 0;
0030 %  for sam = linspace(0,pi,6)
0031     do_xrg
0032     %  end
0033 end
0034 
0035 if 1
0036   close all;
0037 
0038   % t = 950:1425; % Samples
0039   % t0 = 1060;    % Perturbation
0040 
0041   % Problem size
0042   N = size(phi,2);
0043   Nt = size(phi,1);
0044 
0045   % Detrend (project constant angular frequency hypothesis forward in time)
0046   % using robust linear regression
0047   samples = t(1):t0+ofs;
0048   for f = 1:Nt
0049     detrend(f,:) = robustfit(samples,phi(f,samples)).';
0050 %     detrend(f,2) = diff(phi(f,samples([1,end])))/length(samples);
0051 %     detrend(f,1) = phi(f,samples(1))-detrend(f,2)*samples(1);
0052   end
0053 
0054   % Detrend phase
0055   phires = phi - (detrend(:,2)*(1:N) + detrend(:,1)*ones(1,N));
0056   phires0 = phires;
0057   save('phires0.mat','phires0');  
0058   
0059   % Remove certain datasets
0060   phires([1,35,39],:) = nan;
0061 
0062   % Drop data before perturbation
0063   phires = phires(:,t0:t(end));
0064   N = size(phires,2);
0065 
0066 %   figure;
0067 %   plot(phires');
0068 %   return;
0069 
0070   % Phase at perturbation
0071   p0 = mod(mean(dp(:,t0+[-13:13]+ofs),2),2*pi);
0072   save('p0.mat','p0','dp');
0073     
0074 %   trials = ceil(rand(1,25)*size(phires,1));
0075   trials = 1:length(p0);
0076 
0077   % Create matrix of random phase offsets
0078   Nrands = 100;
0079   dp0 = 2*pi*rand(length(p0),Nrands+1);
0080   dp0(:,1) = 0;
0081   
0082   % Trials to include in analysis
0083   good = ~isnan(phires(:,end));
0084   
0085   % Initialize matrix of classification qualities
0086   q = [];
0087   
0088   pr = phires(trials,50:125);
0089   save('pr.mat','pr','good');
0090   
0091   tic
0092   % Loop through perturbation phase randomizations
0093   for n = 1:Nrands+1
0094     
0095     if n == 1
0096       figure;
0097       plot(phires(trials,50:125)')
0098     end
0099     
0100     [c,q(n,:)] = phaseClass(p0(trials)+dp0(trials,n), phires(trials,50:125), good(trials));
0101     
0102     % Store classification for real data
0103     if n == 1
0104       c0 = c;
0105     end
0106     if mod(n,10)==0
0107         dt = toc;
0108         fprintf('.');
0109         if mod(n,100)==0
0110             fprintf(' %d\n',n);
0111         end
0112         if dt>0.2
0113             pause(dt/20)
0114         end
0115         tic
0116     end
0117   end
0118   fprintf('done\n');
0119   dt = toc;
0120   
0121 
0122   figure;
0123   hold on;
0124   plot((1:size(q,2))*180/size(q,2),q(1,:),'k','LineWidth',4);
0125   xlabel('Partition angle (degrees)', 'FontSize', 18,'FontName','Times');
0126   ylabel('Quality of Classification','FontSize',18,'FontName','Times');
0127 
0128   plot((1:size(q,2))*180/size(q,2),q(2:end,:),'color',0.75*[1,1,1]);
0129 
0130   ylim = get(gca,'YLim');
0131   set(gca,'XLim',[0,180],'FontSize',14,'FontName','Times');
0132 %   set(gca,'XLim',[0,180],'YLim',[0, ylim(2)], 'FontSize',14,'FontName','Times');
0133   title(['<<foo>>',' ',num2str(Nrands),' random offsets'],'FontSize',18,'FontName','Times');
0134     
0135 %  legend({'Coherent at Perturbation','Incoherent at Perturbation'}, 'Location','SouthWest')
0136   % Compute empirical distribution function
0137   [F,X] = ecdf(max(q,[],2));
0138 
0139   pvalue = (length(find(max(q,[],2) >= max(q(1,:))))-1)/Nrands;
0140   disp(['P-value = ',num2str(pvalue)]);
0141 
0142 %   figure;
0143 %   hold on;
0144 %   plot(X,F,'LineWidth',4)
0145 %   ylim = get(gca,'YLim');
0146 %   plot(max(q(1,:))*[1,1],ylim,'k--','LineWidth',4);
0147 %   xlabel('Maximum classification error', 'FontSize', 18,'FontName','Times');
0148 %   ylabel('Empirical Probability of Minimum','FontSize',18,'FontName','Times');
0149 %
0150 %   set(gca,'FontSize',14,'FontName','Times');
0151 %   title([dataset,' EDF of classification error, N = ',num2str(Nrands), ', P = ',num2str(pvalue)],'FontSize',18,'FontName','Times');
0152 
0153   
0154   
0155 %   close all;
0156 %   figure;
0157 %   hold on;
0158 %
0159 %   q = [];
0160 %
0161 %   for k = 1:300
0162 %
0163 %     trials = ceil(rand(1,35)*size(phires,1));
0164 %
0165 %     % Trials to include
0166 %     inc = ~isnan(phires(trials,end));
0167 % %     inc = inc(trials);
0168 %
0169 %     tic
0170 %     [c,q(k,:),Q] = phaseClass(p0(trials), phires(trials,50:125), inc);
0171 %     q(k,:) = q(k,:) - mean(q(k,:));
0172 % %     disp(['Phase classification took ',num2str(toc),' seconds']);
0173 %
0174 %   end
0175 %
0176 % %   plot(q')
0177 %
0178 %   plot(mean(q));
0179 %   plot(mean(q)+std(q),'--');
0180 %   plot(mean(q)-std(q),'--');
0181 
0182   
0183 %   figure;
0184 %   boxplot(Q);
0185 end

Generated on Mon 02-Aug-2010 16:44:38 by m2html © 2003