%simulate a backprop network on an XOR task

%set up trial sequence
n = 50000; %number of trials
stimVals = [0 1 1 0 1 1 1 1 1 1; ... %top horizontal line (values for numerals 1,2,3,4,5,6,7,8,9,0)
            1 1 1 1 0 0 1 1 1 1; ... %top-right vertical line
            1 0 1 1 1 1 1 1 1 1; ... %bottom-right vertical line
            0 1 1 0 1 1 0 1 1 1; ... %bottom horizontal line
            0 1 0 0 0 1 0 1 0 1; ... %bottom-left vertical line
            0 0 0 1 1 1 0 1 1 1; ... %top-left vertical line
            0 1 1 1 1 1 0 1 1 0]; %middle horizontal line
teachingVals = eye(10); %for each stimulus, correct response vector has 1 corresponding to entry for that stimulus and zero elsewhere
nstim = size(stimVals,2); %number of stimulus types
stimSeq = randi(nstim,n,1); %random sequence of n stimuli

%set up network
nhidden = [200]; %vector of numbers of hidden units per layer, for layers 1 through M-1
units = [size(stimVals,1) nhidden size(teachingVals,1)]; %number of units for layers 0 through M; layers 0 and M determined by sizes of stimVals and teachingVals, respectively
M = length(units)-1; %number of weight layers
W = cell(M,1); %weight matrices; W is a cell array, each cell containing the weight matrix for one layer; bias term (b) is appended as final column of each matrix
for m=1:M
    W{m} = randn(units(m+1),units(m)+1); %create weight matrix for layer m; (size of postsynaptic layer) x (size of presynaptic layer plus bias term)
end
e = .01; %learning rate

%simulate trial sequence
perf = zeros(n,nstim); %will track model's test performance on all stimulus types after every learning trial
for i=1:n %loop over trials
    stim = stimVals(:,stimSeq(i)); %input to network for this trial<br>
    T = teachingVals(:,stimSeq(i)); %teaching signal for this trial<br>
    [v,a] = feedforward(stim,W); %get network activation by forward computation<br>
    Delta = backprop(stim,v,a,W,T); %get weight updates by backpropagation
    for m=1:M
        W{m} = W{m} + e*Delta{m}; %update weights
    end
    for j=1:nstim %loop through test trials; will test model's performance on every stimulus type
        [v,a] = feedforward(stimVals(:,j),W); %get network activation for stimulus j by forward computation
        perf(i,j) = sum((v{M}-teachingVals(:,j)).^2); %sum squared error: activation at output layer vs. correct response
    end
end

plot(perf) %shows one learning curve (error as a function of trials) for each stimulus; perfect learning is indicated by all curves reaching zero
axis([0 n 0 2]) %fix vertical scaling so that convergence to zero is visible
