clc
clear all
close all
%%
 I1 = importdata('WEAKV1.txt');
 V = I1(:,2)  ;
 time = I1(:,1) ;
 
 
%%
 load ('V')
 load ('N')
%%
% I1 = importdata('V22.txt');
%  V = I1(:,2)  ;
%  time = I1(:,1) ;
% 
  I2 = importdata('N22.txt');
  N = I2(:,2)  ;


W=zeros(length(N),1)  ;

MEAN = mean(N) ;

N = N - MEAN ;

%%
ms = 10 ;
T = 160*10 ;

c=1 ;
p=1 ;

w=1 ;

% figure
% findpeaks(V,'MinPeakDistance',1600, 'MinPeakHeight',20)  ;
% hold on
%%
[pks , locs] = findpeaks(V,'MinPeakDistance',1600, 'MinPeakHeight',20)  ;
%

K = find( pks>20) ;
length(K)

%
for i=3:length(pks)
     
  if pks(i,1) > 20
        
  A(1,:) = V(locs(i,1)-100:locs(i,1) , 1) ;
  AA = A(A>=-20) ;
  [L]=knnsearch(AA',-20) ;
       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 

%%

stc_Matrix = S ;
m_s = mean(stc_Matrix,2);


time = -time(1:160*10+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(:,1601);%
v2 = U(:,1600);%
%v3 = U(:,2);%


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

legend('STA','feature 1', 'feature 2')
xlabel(' t before spike (ms)')

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

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

plot(length(d):-1:1,d(1:end)/norm(d,2),'og','LineWidth',1)


hold on

plot(1 , d(1601)/norm(d,2),'oc','LineWidth',1)
plot(2 , d(1600)/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.25 0.25])
ylim([-0.25 0.25])

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


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