function prob2a
% Problem Two

% Graph the error vector length vs. iterations count for a fixed value of the
% learning rate.
iterations = 50:200;
learningRate = 10:10:50;
for i=1:length(iterations);
for j=1:length(learningRate);
errorLength(i,j) = backpropMLP(iterations(i),learningRate(j));
end;
end;
surf(learningRate, iterations, errorLength)
shading interp
view(0,90)
hold on
errorLength=errorLength+.01;
contour3(learningRate, iterations, errorLength)
hold off
return

function errorLength = backpropMLP(iterations, learningRate)

% This program uses back-propagation to train a neural network to respond to a
% given (single) input vector with a given (single) output vector.
%
% iterations is the number of times that the weights will be updated.
% learningRate is the "strength" of each weight update.
%
% this program returns the length of the final error vector

% Set up some constants
nodeLayers = 2; % number of node layers
nodeCount = [2, 2]; % number of nodes for each layer (comma-separated)
%learningRate = 10.0; % learning rate
logisticWeight = 1; % weight of the logistic activation functions
minError = 1e-8; % iterate 'till error vector length is less than this
input{1} = [1]; % elements of training data for input one (comma-separated)
input{2} = [3]; % elements of training data for input two (comma-separated)
desired{1} = [0.95]; % desired response for output one (comma-separated)
desired{2} = [0.05]; % desired response for output two (comma-separated)

% Make sure our inputs and desired outputs are the same size
trainingLength = length(input{1});
for i=1:length(input);
if length(input{i}) ~= trainingLength;
error('The input vectors must all be the same size.');
end;
end;
for i=1:length(desired);
if length(desired{i}) ~= trainingLength;
error('The output vectors must all be the same size.');
end;
end;

% Set things up so that the inputs act like a first layer of nodes - since
% outputs are denoted by 'y,' these are the outputs of layer 1, e.g. y{1}.
nodeLayers = nodeLayers + 1;
nodeCount = [length(input) nodeCount];

% Make room for a bias input in each layer (except the output layer)
nodeCount = nodeCount + 1;
nodeCount(length(nodeCount)) = nodeCount(length(nodeCount)) - 1;

% Set up the weights
% W{} contains a set of matrices, one for each layer. The matrix for
% weight-layer k contains elements W{k}(j,i) connecting the output of
% neuron i to the input of neuron j. The weights are initially random.
for layer=2:nodeLayers;
W{layer-1} = 2*rand(nodeCount(layer), nodeCount(layer-1))-1;
end;

% Set the weights manually
W{1}(1,1) = 1;
W{1}(1,2) = -2;
W{1}(1,3) = -3;
W{1}(2,1) = 2;
W{1}(2,2) = 0;
W{1}(2,3) = 1;
W{2}(1,1) = 1;
W{2}(1,2) = 0;
W{2}(1,3) = 2;
W{2}(2,1) = 1;
W{2}(2,2) = -2;
W{2}(2,3) = -3;

% Iterate over the set of training data, updating the weights each time
for iteration = 1:iterations;

% Get this input vector and desired output vector
for i=1:length(input);
x(i) = input{i}(1);
d(i) = desired{i}(1);
end;

% Insert the input as the 1st layer
y{1} = x;

% Process each layer of weights
for layer=1:nodeLayers-1;

% For this layer, find the output of each destination node
for outNode = 1:nodeCount(layer+1);

% Find the local field for each destination node
v{layer}(outNode) = 0;
for inNode = 1:nodeCount(layer);

% If this input is the bias, make sure its value is one
if inNode == nodeCount(layer);
y{layer}(inNode) = 1;
end;

% Update the local field
v{layer}(outNode) = v{layer}(outNode) + ...
y{layer}(inNode)*W{layer}(outNode, inNode);
end;

% Apply a logistic activation function to the local field
y{layer+1}(outNode) = activation(v{layer}(outNode), logisticWeight);

end;
end;

% Back propagate
for layer=nodeLayers:-1:2;

for jNode = 1:nodeCount(layer);

% Find the 'error' of the next layer and store it so we can use it
% to determine whether we're done iterating.
if layer == nodeLayers;
e = d(jNode) - y{layer}(jNode);
error(jNode) = e;
else
e = 0;
for kNode = 1:nodeCount(layer+1);
e = e + delta{layer}(kNode) * W{layer}(kNode, jNode);
end;
end;

% Find the derivative of the activation function
phiPrime = logisticWeight * y{layer}(jNode) * ( 1 - y{layer}(jNode) );

% Compute the local gradient for this node
delta{layer-1}(jNode) = e*phiPrime;
end;
end;

% Update the weights
for layer=1:nodeLayers-1;
for inNode = 1:nodeCount(layer);
for outNode = 1:nodeCount(layer+1);
W{layer}(outNode, inNode) = W{layer}(outNode, inNode) ...
+ learningRate ...
* delta{layer}(outNode) ...
* y{layer}(inNode);
end;
end;
end;

% We're done iterating if the length of the error vector drops below minError
errorLength = sqrt(dot(error, error));
% if errorLength < minError;
% break;
% end;
end; % the iteration over training data
return