function eegtutorial

datapath = './data/';
coord = -load([datapath 'coordinates.txt']);

D=64;
fs = 250; %Hz
T = 176;

%fid = fopen([datapath '4_04_errors(81).bin'],'r','b');N1 = 81;
fid = fopen([datapath '4_02_errors(46).bin'],'r','b');N1=46;

error = fread(fid,[128,T*N1],'int16');
error = error(1:D,:);
error = reshape(error,[D,T,N1]);
fclose(fid);

fid = fopen([datapath '4_04_corrects(547).bin'],'r','b');N2 = 547;
%fid = fopen([datapath '4_02_corrects(550).bin'],'r','b');N2=550;
correct = zeros(D,T,N2);
for j=1:N2, for i=1:T,
        correct(:,i,j) = fread(fid,D,'int16');
        fread(fid,D,'int16');
    end; end;
fclose(fid);

t=(1:T)/fs-0.20;
subplot(2,1,1); plot(t,mean(error,3)');
title('ERP condition 1'); xlim([t(1) t(end)])
subplot(2,1,2); plot(t,mean(correct,3)');
title('ERP condition 2'); xlim([t(1) t(end)])


stuffall{1} = 'slice';
stuffall{end+1} = 'trial-averaged-erp';
stuffall{end+1} = 'space-averaged-erp';
stuffall{end+1} = 'difference-erp';
stuffall{end+1} = 'maximum-t2';
stuffall{end+1} = 'matched-filter';
stuffall{end+1} = 'overtraining';
stuffall{end+1} = 'matched-filter-difference';
stuffall{end+1} = 'linear-discrimination';
stuffall{end+1} = 'backward-forward-model';
stuffall{end+1} = 'DSS';
stuffall{end+1} = 'PCA';
stuffall{end+1} = 'fisher';
stuffall{end+1} = 'logist';


for istuff=1:length(stuffall)
    stuff = stuffall{istuff};
    clf
    
    switch stuff
        
        case 'slice'
            subplot(2,2,2)
            imagesc(squeeze(error(12,:,:))')
            xlabel('samples'); ylabel('trials')
            title('electrode #12');
            caxis([-250 250])
            subplot(2,2,3)
            imagesc(squeeze(error(:,75,:)))
            xlabel('trials'); ylabel('electrodes')
            title('sample #75');
            h=colorbar('westoutside'); title(h,'\muV')
            caxis([-250 250])
            subplot(2,2,4)
            imagesc(squeeze(error(:,:,25)))
            xlabel('samples'); ylabel('electrodes')
            title('trial #25');
            caxis([-250 250])
            
        case 'difference-erp'
            subplot(2,2,1); plot(mean(error,3)');
            title('ERP condition 1'); xlim([1 T])
            subplot(2,2,3); plot(mean(correct,3)');
            title('ERP condition 2'); xlim([1 T])
            xlabel('samples');
            d = mean(error,3)-mean(correct,3);
            subplot(2,2,2); plot(d');
            title('ERP difference');  xlim([1 T])
            subplot(2,2,4); imagesc(d);
            title('ERP difference');
            xlabel('samples'); ylabel('electrodes');
            
        case 'trial-averaged-erp'
            subplot(2,2,1)
            plot(mean(error,3)');
            axis tight
            xlabel('samples'); ylabel('\muV');
            title('mean across trials')
            subplot(2,2,3)
            imagesc(mean(error,3));
            xlabel('samples'); ylabel('electrodes');
            title('mean across trials')
            h=colorbar('west'); title(h,'\muV')
            caxis([-1 1]*max(abs(caxis)))
            
        case 'space-averaged-erp'
            w = ones(D,1);
            y = reshape(w'*error(:,:),[T N1]);
            showerp(y,0,w,coord);
            
        case 'matched-filter'
            x = error(:,64+(-5:5),:);
            w = mean(x(:,:),2);
            w = w/norm(w);
            y = reshape(w'*error(:,:),[T N1]);
            showerp(y,0,w,coord);

        case 'fisher'
            win = 50:75;
            [D T N1]=size(error);
            [D T N2]=size(correct);
            mu1=mean(mean(error(:,win,:),3),2);
            mu2=mean(mean(correct(:,win,:),3),2);
            tmp = error(:,win,:);   R1 = cov(tmp(:,:)');
            tmp = correct(:,win,:); R2 = cov(tmp(:,:)');
            Rpool =0.5 *(R1+R2);
            w = inv(Rpool)*(mu1-mu2);
            y2 = reshape(w'*error(:,:),[T N1]);
            y1 = reshape(w'*correct(:,:),[T N2]);
            showerp2(y1,y2,0,w,coord);
            
          case 'logist'
            win = 50:75; L=length(win);
            [D T N1]=size(error);
            [D T N2]=size(correct);
            x1 = error(:,win,:);   x1 = x1(:,:)';
            x2 = correct(:,win,:); x2 = x2(:,:)';
            w = logist([x1;x2],[ones(N1*L,1); zeros(N2*L,1)]);
            y2 = reshape(w(1:D)'*error(:,:),[T N1]);
            y1 = reshape(w(1:D)'*correct(:,:),[T N2]);
            showerp2(y1,y2,0,w(1:D),coord);
            
        case 'overtraining'
            x = error(:,64+(-5:5),:); x=x(:,:);
            Nhalf=round(N1/2);
            w = mean(x(:,1:Nhalf*11),2);
            w = w/norm(w);
            ytrain = reshape(w'*error(:,      1:Nhalf*T),[T     Nhalf ]);
            ytest  = reshape(w'*error(:,Nhalf*T+1:end  ),[T (N1-Nhalf)]);
            showerp2(ytrain,ytest,0,w,coord);
            
        case 'maximum-t2'
            error = error(:,:,1:end-1); N1=N1-1;% last trial is bad.
            x = error(:,64+(-5:5),:);
            w = inv(cov(x(:,:)'))*mean(x(:,:),2);
            y = reshape(w'*error(:,:),[T N1]);
            showerp(y,0,w,coord);
            
        case 'backward-forward-model'
            error = error(:,:,1:end-1); N1=N1-1;% last trial is bad.
            x = error(:,64+(-5:5),:);
            w = inv(cov(x(:,:)'))*mean(x(:,:),2);
            y = w'*x(:,:);
            a=x(:,:)/y;
            subplot(2,2,1); scalp(coord,w);title('backward model')
            subplot(2,2,3); scalp(coord,a);title('forward model')
            
        case 'matched-filter-difference'
            indx=64+(-5:5);
            m1 = mean(  error(:,indx,:),3);
            m2 = mean(correct(:,indx,:),3);
            w = mean(m1-m2,2); w = w/norm(w);
            y1 = reshape(w'*  error(:,:),[T N1]);
            y2 = reshape(w'*correct(:,:),[T N2]);
            showerp2(y1,y2,0,w,coord);
            
        case 'linear-discrimination'
            N=15;
            x1 = randn(N,2);
            x2 = randn(N,2)+ones(N,1)*[2 2];
            x2(N,:)=[2 10]*1.5;
            
            subplot(2,4,2)
            x=[x1;x2];
            y=[ones(N,1); zeros(N,1)];
            w=logist(x,y);
            scatterplot(x,y,w)
            title('Logistic Regression')
            
            subplot(2,4,1)
            w = inv(cov(x1)+cov(x2))*(mean(x1)-mean(x2))';
            w(end+1) = -mean(x*w);
            scatterplot(x,y,w)
            title('Fisher Linear Discriminant')
            
        case 'PCA'
            [W,D] = eig(cov(error(:,:)'));
            [tmp,indx]=sort(-diag(D)); W = W(:,indx);
            Y = reshape(W'*error(:,:),size(error));
            A = inv(W'); % obviously A = W
            K=4;
            for i=1:K
                subplot(K,2,2*i-1); ploterp(squeeze(Y(i,:,:))');
                title(['y_' num2str(i)' '(t)'])
                subplot(K,4,4*i-1); scalp(coord,A(:,i));
                title(['{\bf a}_' num2str(i)'])
            end
            
        case 'DSS'
            [W,D] = eig(cov(mean(error,3)'),cov(error(:,:)'));
            [tmp,indx]=sort(-diag(D)); W = W(:,indx);
            Y = reshape(W'*error(:,:),size(error));
            A = inv(W');
            K=4;
            for i=1:K
                subplot(K,2,2*i-1); ploterp(squeeze(Y(i,:,:))');
                title(['y_' num2str(i)' '(t)'])
                subplot(K,4,4*i-1); scalp(coord,A(:,i));
                title(['{\bf a}_' num2str(i)'])
            end
            
    end
    drawnow
    saveas(gcf,stuff,'epsc');
    system(['convert -density 300x300 ' stuff '.eps ' stuff '.jpg']);
end
end

function    scatterplot(x,y,v)
x = [x ones(size(x,1),1)];
ax=[min(x(:,1)), max(x(:,1)), min(x(:,2)), max(x(:,2))];
plot(x(y==0,1),x(y==0,2),'*b',x(y==1,1),x(y==1,2),'+r');
hold on;
tmean=mean(x);
tmp = tmean; tmp(1)=0; t1=tmp; t1(2)=ax(3); t2=tmp; t2(2)=ax(4);
xmin=median([ax(1), -(t1*v)/v(1), -(t2*v)/v(1)]);
xmax=median([ax(2), -(t1*v)/v(1), -(t2*v)/v(1)]);
tmp = tmean; tmp(2)=0; t1=tmp; t1(1)=ax(1); t2=tmp; t2(1)=ax(2);
ymin=median([ax(3), -(t1*v)/v(2), -(t2*v)/v(2)]);
ymax=median([ax(4), -(t1*v)/v(2), -(t2*v)/v(2)]);
if v(1)*v(2)>0, tmp=xmax;xmax=xmin;xmin=tmp;end;
plot([xmin xmax],[ymin ymax],'k'); hold off
axis tight
ax = axis; axis([ax(1)-1 ax(2)+1 ax(3)-1 ax(4)+1]);
end



function showerp(y,showbar,w,coord)
subplot(3,2,1)
imagesc(y');
xlabel('samples'); ylabel('trials');
if showbar, h=colorbar('east'); title(h,'\muV'); end
caxis([-1 1]*max(abs(caxis)))
title('{\bf y}(t)')
subplot(3,2,3)
ploterp(y',[],0,1:size(y,1))
xlabel('samples');
legend('mean','SEM')
subplot(3,2,5)
scalp(coord,w);
caxis([-1 1]*max(abs(caxis)));
colorbar
title('{\bf w}')
end

function showerp2(y1,y2,showbar,w,coord)
subplot(3,2,2)
imagesc(y1');
xlabel('samples'); ylabel('trials');
if showbar, h=colorbar('east'); title(h,'\muV'); end
caxis([-1 1]*max(abs(caxis)))
title('{\bf y1}(t)')
subplot(3,2,4)
imagesc(y2');
xlabel('samples'); ylabel('trials');
if showbar, h=colorbar('east'); title(h,'\muV'); end
caxis([-1 1]*max(abs(caxis)))
title('{\bf y2}(t)')
subplot(3,2,6)
ploterp(y1',[],0,1:size(y1,1),[0 0 1]); hold on
ploterp(y2',[],0,1:size(y2,1),[1 0 0]); hold off
xlabel('samples');
ylabel('{\bf y1}(t) and {\bf y2}(t)')
subplot(3,2,5);
scalp(coord,w);
caxis([-1 1]*max(abs(caxis)));
colorbar
title('{\bf w}')
end

function ploterp(X,X2,offset,times,color)
% ploterp(X,X2,offset,time,color) X and X2 should be trials by
% samples. offset draws a line at that time. times is a vector that
% tells the times of each sample for plotting. color is a 3-vector
% indicating rgb values for the line and SEM area.

if nargin<4 | isempty(times), times=1:size(X,2); end
if nargin<5 | isempty(color), 
  linecolor = [0 0 1]; backcolor=[0.8 0.8 0.95]; 
else
  linecolor = color; backcolor=[0.8 0.8 0.8]+color/5; 
end

[N,L] = size(X);

if nargin<2 | isempty(X2)
  m = mean(X);
  se = std(X)/sqrt(size(X,1));
else
  m = mean(X)-mean(X2);
  se = sqrt((std(X)/sqrt(size(X,1))).^2+(std(X2)/sqrt(size(X,1))).^2);
end

plot(times([1 L]), [0 0],'k'); hold on;

fill([times(1:L)  times(L:-1:1)],[m+2*se m(end:-1:1)-2*se(end:-1:1)], ...
    backcolor,'LineStyle','n')
axis tight; 
ax=axis;
plot(times,m,'color',linecolor);    
hold off

if nargin>2 & ~isempty(offset) 
  hold on;
  plot([offset offset],ax(3:4),'k')
  axis(ax);  
  hold off
end

end

function scalp(e,a,plote)

% function scalp(electrodes,a) 
% plots scalp plot
% based on project a and electrode positions in e
% if plote is 1, then it plots the location of the electrodes

  if nargin < 3
    plote = 0;
  end
  
  a = a(:);
  if size(e,1) ~= size(a,1)
    error(['Number of colums in e must be the same as the length' ...
	  ' of a']);
  end
  
  [Xi,Yi] = meshgrid(linspace(min(e(:,1)), max(e(:,1)),50), ...
	  linspace(min(e(:,2)), max(e(:,2)),50));
  [Zi] =  griddata(e(:,1), e(:,2),a,Xi,Yi);
  maxz = max(max(Zi));
  
  surf(Xi,Yi,zeros(size(Zi)),Zi,'EdgeColor','none');
  view(2)
  %shading('interp')

  axis off;
  xlabel('x');
  
  if plote
    hold on
    
    plot3(e(:,1),e(:,2),maxz*ones(size(e(:,1))),'*');
    
    hold off
    text(4.75,10.5,'back');
    text(4.85,0.5,'front');
    text(0.45,5.5,'left');
    text(9.2,5.5,'right');
  
  end
end




