function y=prob2b(iterations, learningRate)

input{1} = [0 0 1 1];   % elements of training data for input one (comma-separated)
input{2} = [0 1 0 1];   % elements of training data for input two (comma-separated)
desired{1} = [1 0 0 1]; % desired response for output one (comma-separated)
desired{2} = [0 1 1 0]; % desired response for output two (comma-separated)
nodeCount = [2, 2];        % number of nodes for each non-input layer (comma-separated)
logisticWeight = 1;        % weight of the logistic activation functions

% 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}.
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;

% Find the weights that properly satisfy the XOR function
W = computeWeights( iterations, learningRate, nodeCount, input, desired,     ...
                    logisticWeight);
                
% Compute the output for every possible input
x1 = 0:.05:1;
x2 = 0:.05:1;
for i=1:length(x1);
    for j=1:length(x2);
        y = computeOutputs(nodeCount, [x1(i) x2(j)], W, logisticWeight);
        output(i,j,1) = y{3}(1);
        output(i,j,2) = y{3}(2);
    end;
end;
figure(1);
surf(x1, x2, output(:,:,1));
shading interp;
figure(2);
surf(x1, x2, output(:,:,2));
shading interp;

return;





function W = computeWeights(iterations, learningRate, nodeCount, input,      ...
                            desired, logisticWeight)
% This program uses back-propagation to train a neural network to perform the 
% XOR operation
%
% iterations is the number of times that the weights will be updated.
% learningRate is the "strength" of each weight update.
% nodeCount is an array containing the number of nodes per layer (including inputs)
% input is a cell array of training data input components
% desired is a cell array of training data desired components
% logisticWeight is the weight of the logistic activation functions
%
% this program returns the weight of every connection

% Set up some constants
weightRange = 2;           % weights will be initialized with equal probability in 
weightCenter = 0;          %    the range weightRange centered at weightCenter
rand('state',0);



% 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 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.
for layer=2:length(nodeCount);
    W{layer-1} = weightRange*rand(nodeCount(layer), nodeCount(layer-1))     ...
                 + (weightCenter - 0.5);
end;

% Iterate over the set of training data, updating the weights each time
for iteration = 1:iterations;
    
    % Randomize the order of our input data and our desired response.  This works 
    % by randomly reordering rows of the 'concatenation' matrix.  The vector 
    % 'order' contains the new row order.  The 'transform' matrix is an identity 
    % matrix in which each diagonal is shifted to the new position of the 
    % corresponding row.  Premultiplying this reorders 'concatenation.'  The 
    % X and d matrices can then be pulled back out, in a re-ordered form.
    order = randperm(length(input{1}));
    transform = zeros(length(input{1}));
    for i=1:length(input{1});
        transform(order(i), i) = 1;
    end;
    for i=1:length(input)
        input{i} = input{i} * transform;
        desired{i} = desired{i} * transform;
    end;
    
    % Present each element in the training set once
    for trainingElement = 1:length(input{1});
        
        % Get this input vector and desired output vector
		for i=1:length(input);
            x(i) = input{i}(trainingElement);
            d(i) = desired{i}(trainingElement);
		end;
        
        % Find the output for this input
        y = computeOutputs(nodeCount, x, W, logisticWeight);

		% Back propagate
		for layer=length(nodeCount):-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 == length(nodeCount);
                    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:length(nodeCount)-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;
    
        % Finally, find the length of the error vector
        errorLength = sqrt(dot(error, error));
        
    end;    % the iteration over training elements
end;        % all iterations
return






function y = computeOutputs(nodeCount, x, W, logisticWeight)

% This function will compute the output at every node for a system with 
% the given geometry, input, weights, and logistic weighting (in the activation
% function).

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

% Process each layer of weights
for layer=1:length(nodeCount)-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;

return
