%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% 	  MODULE: bpWMissingBackProp.m
%%
%%      FUNCTION: code for problem 5 of homework assignment 3
%%		  missing backprop code
%%
%% last modified: 10/16/1997
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% constants
eta = 0.01;		% learning rate
maxIter=10000;		% number of iterations
numDataPoints=200;	% number of data points
errPlotInt=50;		% plot averaged error every errPlotInt iterations
numInp=5;		% # of input units
numHid=2;		% # of hidden units
numOut=5;		% # of output units

% initialize weights with random values around 0
WHI=rand(numHid,numInp+1)*0.1; %weights from input to hidden layer (incl. bias)
WOH=rand(numOut,numHid+1)*0.1; %weights from hidden to output layer (incl. bias)

% create data set
d=createGaussData([0 0 0 0 0]',[1 1 0.25 .1 .1]',numDataPoints);
oos=createGaussData([0 0 0 0 0]',[1 1 0.25 .1 .1]',numDataPoints);
% uncomment next line for rotation of data set
% d=[.707 -.707 0; .707 .707 0; 0 0 1]*d;

figure(1);clf;hold on;
title(sprintf('error averaged over last %d epochs',errPlotInt))

% train using online learning
avgErr=0;
avgoosErr=0;
for iter=1:maxIter
  % pick stimulus at random & do forward pass
  t=d(:,floor(rand*numDataPoints)+1);
  x=[1 t']';	% including bias
  z=[1 g(WHI*x)']';	% hidden activation
  y=WOH*z;

  oost=d(:,floor(rand*numDataPoints)+1);
  oosx=[1 oost']';	% including bias
  oosz=[1 g(WHI*oosx)']';	% hidden activation
  oosy=WOH*oosz;


  % plot weights, error every errPlotInt iterations
  avgErr=avgErr+sum((y-t).^2);
  avgoosErr = avgoosErr+sum((oosy-oost).^2);
  if(rem(iter,errPlotInt)==0)
    figure(1);plot(iter,avgErr/errPlotInt,'+');
    plot(iter,avgoosErr/errPlotInt,'blue*');drawnow;
    avgErr=0;
    avgoosErr = 0;
    figure(2);clf;
    axes('position',[0.05,0.75,0.9,0.2]);hinton(WHI);title('WHI');
    axes('position',[0.05,0.35,0.9,0.2]);hinton(WOH);title('WOH');
  end

  % now backward pass: calculate deltas and adapt
  % Oops, the evil TA took it out...
  deltaO = y-t;
  z2 = z(2:size(z));
  tmp = WOH'*deltaO;
  deltaH = (z2.*(1-z2)).*tmp(2:size(tmp));
  WHI = WHI - eta*deltaH*x';
  WOH = WOH - eta*deltaO*z';
end
figure(1);hold off;figure(2);hold off;

