% clc
% clear all
% close all
%%

% 
%  I7 = importdata('VData.txt');
%  V = I7(:,2)  ;
%  L4 = length(V) ;
% 
%  I8 = importdata('NData.txt');
%  N = I8(:,2)  ;
% E4 = length(N) ;
%  
% 
% 
% W=zeros(length(N),1)  ;
% 
% MEAN = mean(N) ;
% 
% N = N - MEAN ;
% 
% %%
% ms = 20 ;
% T = 40*20 ;
% 
% c=1 ;
% p=1 ;
% 
% w=1 ;
% 
% % figure
% % findpeaks(V,'MinPeakDistance',1600, 'MinPeakHeight',20)  ;
% % hold on
% %%
% % [pks , locs] = findpeaks(V,'MinPeakDistance',2600, 'MinPeakHeight',0)  ;
% % %
% % 
% % K = find( pks>0) ;
% % length(K)
% % 
% % %
% % for i=3:length(pks)
% %      
% %   if pks(i,1) > 0
% %         
% %   A(1,:) = V(locs(i,1)-100:locs(i,1) , 1) ;
% %   AA = A(A>=0) ;
% %   [L]=knnsearch(AA',0) ;
% %        B = AA(1,L) ;
% % 
% %        
% %   [row1]=find( V==B) ;
% %   
% %  
% %   for jjj=1:size(row1,1)
% %   if (row1(jjj,1) < locs(i,1)-5) && (row1(jjj,1) > locs(i,1)-100)
% %   row2 = row1(jjj,1) ;
% %  if V(row2)<5
% %      S(:,c)= N(row2-T:row2,1) ;
% %        c=c+1 ;
% %        
% %        
% %      W(row2,1) = 1  ; 
% %      
% %      
% % %   plot(row2, V(row2) , '*g')
% % %   hold on  
% %      
% %   end 
% %   end      
% %   end
% % end
% % end
% 
% %%
%   %save ('S.mat' , 'S')
%   %save ('W.mat' , 'W')
% %%
% load ('S')
% load ('W')
% 
% %%
% 
% [m,n]=size(S) ;
% 
% NumberSpike = n 
% 
% %%
% 
% stc_Matrix = S ;
% m_s = mean(stc_Matrix,2);
% 
% %%
% 
%  time = I7(:,1) ;
% 
% time = -time(1:40*20+1,1) ;
% time = sort(time) ;
% 
% 
% 
% %%
% 
% cc=1 ;
% spp = W ;
% 
% sp_binary_shuf(:,1) = Shuffle_(spp);
% 
% spike_indx = find( sp_binary_shuf >0);
% 
% maxindx = spike_indx(end);
% indx_int = find(spike_indx>T & spike_indx<maxindx-T);
% 
% for i=1:length( indx_int   )
% 
% SS(:,cc)= N( spike_indx(indx_int(i)) -T: spike_indx(indx_int(i))   ,1) ;
%        cc=cc+1 ;
%     
% end
% 
%  stc_Matrix_shuf = SS ;
% 
% %%
% 
% ccc=1 ;
% spp = W ;
% 
% sp_binary_shuf2(:,1) = Shuffle_(spp);
% 
% spike_indx2 = find( sp_binary_shuf2 >0);
% 
% maxindx2 = spike_indx2(end);
% indx_int2 = find(spike_indx2>T & spike_indx2<maxindx2-T);
% 
% for i=1:length( indx_int2   )
% 
% SSS(:,ccc)= N( spike_indx2(indx_int2(i)) -T: spike_indx2(indx_int2(i))   ,1) ;
%        ccc=ccc+1 ;
%     
% end
% 
%  stc_Matrix_shuf2 = SSS ;
% 
%  
% %% 
%  
% inp_mean = mean(stc_Matrix,2); INP_mean = repmat(inp_mean,1,size(stc_Matrix,2)); 
% inp_mean_shuf = mean(stc_Matrix_shuf,2); INP_mean_shuf = repmat(inp_mean_shuf,1,size(stc_Matrix_shuf,2)); 
% inp_mean_shuf2 = mean(stc_Matrix_shuf2,2); INP_mean_shuf2 = repmat(inp_mean_shuf2,1,size(stc_Matrix_shuf2,2));
% C_shuf2 = cov((stc_Matrix_shuf2 - INP_mean_shuf2)');
% C_shuf = cov((stc_Matrix_shuf - INP_mean_shuf)') - C_shuf2;
% C = cov((stc_Matrix - INP_mean)') - C_shuf2;
% [U_shuf,D_shuf] = eig(C_shuf);%[u_shuf,s_shuf,v] = svd(stc_new_shuf);
% [U,D] = eig(C);%-C_shuf);
% d = diag(D)/norm(diag(D));
% v1 = U(:,1);%
% v2 = U(:,2);%
% %v3 = U(:,2);%


%%
% figure
% plot(time,m_s/norm(m_s),'k','LineWidth',2)
% xlabel(' t before spike (ms)')
% legend ('STA')
%%
figure


 plot(time,v2/norm(v2),'r','LineWidth',2)
 hold on
plot(time,m_s/norm(m_s),'k','LineWidth',2)
hold on, plot(time,v1/norm(v1),'c','LineWidth',2)
legend('STA','feature 1', 'feature 2')
xlabel(' t before spike (ms)')


%ylim([-0.04 0.12])


%%
d_shuf = diag(D_shuf)/norm(diag(D_shuf));

figure


plot(length(d):-1:1,d(1:end)/norm(d,2),'Color',[0.75 0.75 0.75], 'linestyle','none' , 'Marker', 'o')

hold on
plot(length(d_shuf):-1:1,d_shuf(1:end)/norm(d_shuf,2),'-k','LineWidth',2)


hold on

plot(801 , d(1)/norm(d,2),'oc','LineWidth',1)
plot(800 , d(2)/norm(d,2),'or','LineWidth',1)

xlabel('eigenmode index')
ylabel('eigenvalue')


%%
%[s1,s2] = proj1(v1 , v2, S+0.06 ,NumberSpike );
% 
% [s1,s2] = Projection_Vector(W, v1 , v2 , N );
% 
% [s1_shuf,s2_shuf] = Projection_Vector(sp_binary_shuf,v1,v2, N);%
% [s1_shuf2,s2_shuf2] = Projection_Vector(sp_binary_shuf2,v1,v2, N );%

figure



plot(s1_shuf,s2_shuf,'Color',[0.75 0.75 0.75], 'linestyle','none' , 'Marker', '.')
hold on
plot(s1_shuf2,s2_shuf2,'Color',[0.75 0.75 0.75], 'linestyle','none' , 'Marker', '.')

 xlim([-0.2 0.2])
 ylim([-0.2 0.2])

hold on

plot(s1,s2,'.k')

xlabel('Projection onto feature 1')
ylabel('Projection onto feature 2')




% figure
% plot(diag(D_shuf),'og')
% hold on
% plot(diag(D),'or')


% %%
% 
% figure
% 
% plot(time,U(:,1598)/norm(U(:,1599)),'b','LineWidth',2)


% ssp_binary_shuf = find(sp_binary_shuf2>0)  ;
% [ss1 ss2]=Projection_ (ssp_binary_shuf, C_shuf2  , N+0.06) ;
% 
% WW=find(W>0) ;
% [s1,s2] = Projection_(WW,C,N+0.06) ;
% figure
% 
% 
% plot(ss1,ss2,'.k')
% hold on
% plot(s1,s2,'.r')


%%

% %%
% figure
% plot(v1)
% %%
% figure
% plot(S(:,50))


