clc
clear all
close all
%% Importing Data  (Voltage and Noise Data generated in NEURON)  
% 
% I1 = importdata('VAIS1.txt');
% V1 = I1(:,2)  ;
% 
% L1 = length(V1) ;
%  
%  
% I2 = importdata('NAIS1.txt');
% N1 = I2(:,2)  ;
% E1 = length(N1) ;
% 
% I3 = importdata('VAIS2.txt');
% V2 = I3(:,2)  ;
% L2 = length(V2) ;
%  
%  
% I4 = importdata('NAIS2.txt');
% N2 = I4(:,2)  ;
% E2 = length(N2) ;
% 
% I7 = importdata('VAIS4.txt');
% V3 = I7(:,2)  ;
% L3 = length(V3) ;
% time = I7(:,1) ;
%  
% I8 = importdata('NAIS4.txt');
% N3 = I8(:,2)  ;
% E3 = length(N3) ;
%  
%  
% V (1:L1 , 1) = V1  ;
% V (L1+1:L1+L2 , 1) = V2  ; 
% V (L1+L2+1:L1+L2+L3 , 1) = V3  ; 
% %  V (L1+L2+L3+1:L1+L2+L3+L4 , 1) = V4  ; 
% 
% N (1:E1 , 1) = N1  ;
% N (E1+1:E1+E2 , 1) = N2  ;
% N (E1+E2+1:E1+E2+E3 , 1) = N3  ;
% % N (E1+E2+E3+1:E1+E2+E3+E4 , 1) = N4  ;
%  
% 
% 
% 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',400, 'MinPeakHeight',0)  ;
% % hold on
% 
% %%  Finding the spikes and the time of spiking
% 
% % [pks , locs] = findpeaks(V,'MinPeakDistance',400, '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) ;
% % 
% %      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
% % 
% % 
% % save ('S.mat' , 'S')
% % save ('W.mat' , 'W')
% % 
% % 
% % [m,n]=size(S) ;
% % 
% % NumberSpike = n ;
% 
% %%
% 
% 
%  load ('S')
%  load ('W')
% 
% stc_Matrix = S ;
% m_s = mean(stc_Matrix,2);
% 
% 
% time = -time(1:40*20+1,1) ;
% time = sort(time) ;
% 
% 
% 
% %% Shuffling spikes
% 
% 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 ;
% 
%  
% %% Eigenvalues 
%  
% 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')







%%  Plotting STA, feature 1, feature 2


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])


%%  Plotting eigenvalues and eigenmodes

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','LineWidth',1)
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')


%% Plotting projections onto features

%[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')