% Lab 6, problem 2
% you must first load in apples and oranges data

% initialize data array
data=[apples apples2 oranges oranges2];
[N K]=size(data);

% initialize teacher
teacher=[ones(1,K/2) zeros(1,K/2)];
half1=1:K/2;
half2=K/2+half1;

% learning rate
eta=.05;

% number of trials - you may want to make this longer
num_trials=1000;

% initialize weights
v=randn(2,1);
v0=randn(1);
w=randn(2,2);
w0=randn(2,1);

% initialize data plot
figure(1)
clf
subplot(121)
plot(data(1,half1),data(2,half1),'r+',data(1,half2),data(2,half2), 'go')
hold on
x1=[0 5];
x2=-(w(1,1)*x1+w0(1))/w(2,1);
h1=plot(x1,x2,'k--','LineWidth',2,'EraseMode','xor');
x3=-(w(1,2)*x1+w0(2))/w(2,2);
h2=plot(x1,x3,'k-.','LineWidth',2,'EraseMode','xor');
axis image, axis([0 5 -2 3])

subplot(122)
datay=sigmoid(w'*data+w0*ones(1,K));
hy(1)=plot(datay(1,half1),datay(2,half1),'r+','EraseMode','xor');
hold on
hy(2)=plot(datay(1,half2),datay(2,half2),'go','EraseMode','xor');
y1=[-0.1 1.1];
y2=-(v(1)*y1+v0)/v(2);
h3=plot(y1,y2,'LineWidth',2,'EraseMode','xor');
axis image, axis([-0.1 1.1 -0.1 1.1])

% loop over trials
for t=1:num_trials
   
   dw=0;
   dw0=0;
   dv=0;
   dv0=0;
   e2=0;
   
   % loop over training set
   for alpha=1:K

     x=data(:,alpha);			% get this data vector
     nety=w'*x+w0;			% net input for y layer
     y=sigmoid(nety);                   % output for y layer
     netz=v'*y+v0;			% net input for z layer
     z=sigmoid(netz);			% output for z layer
     e=teacher(alpha)-z;	        % error
     e2=e2+e^2;                         % accumulate squared error
     delta_z=e*z*(1-z);                 % delta_z (eq. 9 of handout)
     delta_y=y.*(1-y).*v*delta_z;	% delta_y (eq. 16 of handout)
     dv=dv+2*eta*y*delta_z;		% dv (eq. 10)
     dv0=dv0+2*eta*delta_z;             
     dw=dw+2*eta*x*delta_y';		% dw (eq. 17)
     dw0=dw0+2*eta*delta_y;
   end
   
   % update weights
   w=w+dw;        
   w0=w0+dw0;
   v=v+dv;
   v0=v0+dv0;
   
   % save E for this trial
   E(t)=e2;
   
   % update display of separating hyperplane
   x2=-(w(1,1)*x1+w0(1))/w(2,1);
   set(h1,'YData',x2)
   x3=-(w(1,2)*x1+w0(2))/w(2,2);
   set(h2,'YData',x3)
   
   datay=sigmoid(w'*data+w0*ones(1,K));
   set(hy(1),'XData',datay(1,half1),'YData',datay(2,half1))
   set(hy(2),'XData',datay(1,half2),'YData',datay(2,half2))
   y2=-(v(1)*y1+v0)/v(2);
   set(h3,'YData',y2)
   
   % update E plot
   if t>1
     set(hE,'XData',[t-1 t],'YData',E(t-1:t))
   else
     figure(2)
     hE=plot(1,E(1),'EraseMode','none');
     axis([1 num_trials 0 E(1)])
   end
   
   drawnow
   
end

figure(1)
subplot(121)
hold off
subplot(122)
hold off
