%solution_sheet12

load features_lecture_adaptiveBCI

%Ex 1
icl1=find(fv_train.y(1,:));
icl2=find(fv_train.y(2,:));
icl1_test=find(fv_test.y(1,:));
icl2_test=find(fv_test.y(2,:));

%plot distributions of train and test datasets.
plot(fv_train.x(1,icl1),fv_train.x(2,icl1),'rx')
hold on
plot(fv_train.x(1,icl2),fv_train.x(2,icl2),'go')
plot(fv_test.x(1,icl1_test),fv_test.x(2,icl1_test),'kx')
plot(fv_test.x(1,icl2_test),fv_test.x(2,icl2_test),'bo')

%calculate LDA classifier of train data. new version
mean_train(:,1)=mean(fv_train.x(:,icl1)')';
mean_train(:,2)=mean(fv_train.x(:,icl2)')';

w=pinv(cov(fv_train.x'));
w=-w*(diff(mean_train,1,2));
b=-mean(fv_train.x')*w;

%calculate LDA classifier old version.
w_old=pinv((cov(fv_train.x(:,icl1)')+cov(fv_train.x(:,icl2)'))/2)*(diff(mean_train,1,2));
b_old=-mean(mean_train')*w_old;

%both are equivalent
w./w_old
b/b_old

%plot classifier boundary
yl= (-w(1)*xlim-b)/w(2);
plot(xlim, yl, 'k');
yl= (-w_old(1)*xlim-b_old)/w_old(2);
plot(xlim, yl, 'm');

%calculate LDA error with train classifier.
out=w'*fv_test.x+b*ones(1,size(fv_test.x,2));
out_old=w_old'*fv_test.x+b_old*ones(1,size(fv_test.x,2));
out_label=sign(out);
err_orig=mean(out_label~=[-1 1]*fv_test.y);
out_label=sign(out_old);
err_old=mean(out_label~=[-1 1]*fv_test.y);

%adapt trial by trial
UC=0.01;
%type one, using extended covariance matrix.
[nf,N]=size(fv_train.x);
fvold.x=[ones(1,N) ;fv_train.x];
Cold=fvold.x*fvold.x'/N;
invCold=pinv(Cold);
invCnew=invCold;
factor=1/(1-UC);
out_test=zeros(1,size(fv_test.y,2));
cadap.w=w;
cadap.b=b;
 for itrial=1:size(fv_test.x,2)
     out_test(itrial)=cadap.w'*fv_test.x(:,itrial)+cadap.b;
     v=invCnew*[1 ;fv_test.x(:,itrial)];
     denominator=((1-UC)/UC)+[1;fv_test.x(:,itrial)]'*v;
     invCnew=factor*(invCnew-(v*v'/denominator));
     invCnew=(invCnew+invCnew')/2;
     cadap.w=invCnew(2:end,2:end)*diff(mean_train,1,2);
     cadap.b=invCnew(1,2:end)*diff(mean_train,1,2);
 end;

 yl= (-cadap.w(1)*xlim-cadap.b)/cadap.w(2);
plot(xlim, yl, 'k');

%calculate new errors
out_label=sign(out_test);
err= mean(out_label~=[-1 1]*fv_test.y);

 
 %second way: using normal covariance and updating the pooled mean
 
icovnew=inv(cov(fv_train.x'));
pmold=mean(fv_train.x');
pmean=pmold';
out_test=zeros(1,size(fv_test.y,2));
cadap.w=w;
cadap.b=b;
for itrial=1:size(fv_test.x,2)
     out_test(itrial)=cadap.w'*fv_test.x(:,itrial)+cadap.b;
     pmean=(1-UC)*pmean+UC*fv_test.x(:,itrial);
     sample=fv_test.x(:,itrial)-pmean;
     v=icovnew*sample;
     denominator=((1-UC)/UC)+sample'*v;
     icovnew=factor*(icovnew-(v*v'/denominator));
     icovnew=(icovnew+icovnew')/2;
     cadap.w=icovnew*diff(mean_train,1,2);
     cadap.b=-cadap.w'*pmean;
 end

 yl= (-cadap.w(1)*xlim-cadap.b)/cadap.w(2);
plot(xlim, yl, 'b');

%calculte new errors
out_label=sign(out_test);
err= mean(out_label~=[-1 1]*fv_test.y);
