% AutoAssociationDemo.m 
% Autoassociation or attractor network demonstration program
%
% Please do not redistribute this code or modifications of it without the following notice:
% (c) E T Rolls www.oxcns.org 2015
% The background to these networks is provided in
% Rolls ET (2016) Cerebral Cortex: Principles of Operation (Oxford University Press) Appendix 2 (B)
% which is available at http://www.oxcns.org by following the link to the book.
% The instructions and exercises for this code are in 
% Rolls ET (2016) Cerebral Cortex: Principles of Operation (Oxford University Press) Appendix 4 (D) 
% which is available at http://www.oxcns.org by following the link to the book.
% The background to these networks is also available in
% Rolls ET (2008) Memory, % Attention and Decision-Making (Oxford University Press) Appendix 2, which is available at
% http://www.oxcns.org/papers/RollsMemoryAttentionAndDecisionMakingContents+Appendices1+2.pdf

% This Matlab program is intended to demonstrate some of the principles of how autoassociation or attractor networks operate,
% and for that reason Matlab code optimizations are frequently eschewed in favor of illustrating principles of neuronal operations 
% cd c:\np\papers\visn\compcap
% The synaptic weight matrix (SynMat) has column vectors to represent the synapses on a single neuronn.
% There is one column vector therefore for each of the N neurons.
% Thus SynMat(Row, Column) implements SynMat(synapse, neuron)

clear all; close all hidden; format compact; % format bank;
fig = 1;
% rng(1);           % comment this in only to set the random number generator for code development
N = 100;            % number of neurons in the fully connected autoassociation net
nSyn = N;           % number of synapses on each neuron
SynMat = zeros(nSyn, N); % (Synapses = Rows, neurons = Columns)
nPatts = 10;        % the number of training and testng patterns. Suggest 10 for sparseness = 0.5, and 30 for sparseness = 0.1
Learnrate = 1 / nPatts; % Not of especial significance in an autoassociation network, but this keeps the weights within a range 
Sparseness = 0.5;   % the sparseness of the trained representation, i.e. the proportion of neurons that have high firing rates of 1
                    % Investigate values of 0.5 and 0.1
display = 1;        % 0 no display; 1 display the network
nFlipBits = 14;     % The number of bits that are flipped to produce distorted recall cues.
                    % It is suggested that this be set to close to 0.2 * N * Sparseness
nepochs = 9;       % the number of times that the network is allowed to update during testng

%
% generate random binary training patterns
TrainPatts = zeros(N, nPatts); % This matrix stores the training patterns. Each pattern vector has N elements
for patt = 1 : nPatts
    TrainPatts(1 : N * Sparseness, patt) = 1; % the number of bits set to 1 in each pattern
    p = randperm(N); % rearrange the elements of this pattern vector in random order
    TrainPatts(:, patt) = TrainPatts(p, patt);
end

if display == 1
    fig = fig + 1;
    figure(fig)
    imagesc(TrainPatts);
    colormap(gray)
    title('Random Binary Training Patterns. Each column is a pattern');
    xlabel('Pattern number'); ylabel('Training pattern');
    pause 
end

% Now make some testing patterns with some bits flipped from the training patterns
TrainPattsFlipped = TrainPatts;
for patt = 1 : nPatts
    synarray = randperm(nSyn);
    el = 1;
    for bit = 1 : nFlipBits
        while TrainPatts(synarray(el),patt) ~= 1
            el = el + 1;
            if el > nSyn
                disp('Error: too many bits being flipped');
                el = 1;
            end
        end
        TrainPattsFlipped(synarray(el),patt) = 0;
        el = el + 1;
    end
    synarray = randperm(nSyn);
    el = 1;
    for bit = 1 : nFlipBits
        while TrainPatts(synarray(el),patt) ~= 0
            el = el + 1;
            if el > nSyn
                disp('Error: too many bits being flipped');
                el = 1;
            end
        end
        TrainPattsFlipped(synarray(el),patt) = 1;
        el = el + 1;
    end
end
if display == 1
    fig = fig + 1;
    figure(fig)
    imagesc(TrainPattsFlipped);
    colormap(gray)
    title(['Training patterns with ', num2str(nFlipBits), ' flipped bits for testing']);
    xlabel('Pattern number'); ylabel('Distorted Test pattern');
    pause 
end

% Train the synaptic matrix
for patt = 1 : nPatts
    for neuron = 1 : N
        postSynRate = TrainPatts(neuron, patt); % postsynaptic firing rate. The external input to the neurons.
        for syn = 1 : nSyn
            if syn ~= neuron % avoid self connections of a recurrent collateral axon onto its sending neuron
                preSynRate = TrainPatts(syn, patt); % the presynaptic rate is the same as the postsynaptic rate because of the recurrent collaterals 
                weight_change = Learnrate * (postSynRate - Sparseness) * (preSynRate - Sparseness); % use a covariance rule. 
                % The sparseness is the average firing rate
                % weight_change = Learnrate * (postSynRate) * (preSynRate); % OR use a Hebb rule
                % weight_change = Learnrate * (postSynRate) * (preSynRate - Sparseness); % OR use a Hebb rule but with also heterosynaptic LTD see Rolls (2008) B.3.3.6. 
                SynMat(syn, neuron) = SynMat(syn, neuron) + weight_change;
            end
        end
    end
    if display == 1
        figure(100);
        positionVector1 = [0.2, 0.2, 0.7, 0.7]; % left, bottom, width, height
        subplot('Position', positionVector1); % Top right
        imagesc(SynMat);
        colormap(gray);
        title('synaptic matrix during learning');
        positionVector2 = [0.2, 0.05, 0.7, 0.05]; % left, bottom, width, height
        subplot('Position', positionVector2);
        imagesc(TrainPatts(:, patt)');
        title('Postsynaptic Rate');
        set (gca, 'YtickLabel', '')
        positionVector3 = [0.1, 0.2, 0.03, 0.7]; % left, bottom, width, height
        subplot('Position', positionVector3);
        imagesc(TrainPatts(:, patt));
        title('Presynaptic input');
        set (gca, 'XtickLabel', '');
        drawnow;
    end
    % pause % comment this out if you wish to see the training pattern by pattern
end
%SynMat
disp(['Synaptic matrix minimum= ', num2str(min(min(SynMat))), ' maximum= ', num2str(max(max(SynMat)))]);
% comment in the following two lines to test whether the network operates with positive-only synaptic weights
% SynMat = SynMat - min(min(SynMat)); 
% disp(['Synaptic matrix minimum= ', num2str(min(min(SynMat))), ' maximum= ', num2str(max(max(SynMat)))]);
if display == 1
    pause
end

% Now test the network to show whether it completes incomplete patterns and
% corrects distortions in the applied recall cue
disp('Testing');
nCorrect = 0;   % the number of correctly recalled patterns

for patt = 1 : nPatts
    Rate = zeros(1, N);  % output firing rates
    for epoch = 1 : nepochs
        if epoch <=1
            clamp = 1; % the recall cue or pattern is applied only at the beginning of retrieval
        else
            clamp = 0;
        end
        % calculate activation
        Actvn = zeros(1, N); % the activations of the N output neurons
        for neuron = 1 : N
            PreSynInput = Rate' + (TrainPattsFlipped(:, patt) * clamp); % the firing on the recurrent collaterals is the firing rate of the neurons
            % plus at the beginning of retrieval the cue recall pattern
            Actvn(neuron) = dot(PreSynInput,SynMat(:, neuron)); % dot product of input pattern and synaptic weight vector for one neuron
        end
        scale = 1.0 / (max(Actvn) - min(Actvn));
        Actvn = (Actvn - min(Actvn)) * scale;            % scale Activation to 0-1.
        % Now convert the activation into a firing rate using a binary thereshold activation functon
        % The threshold is selected in this case artificially by sorting
        % the array to determine where the threshold should be set to achieve the required sparseness
        % In the brain, inhibitory neurons are important in setting the proportion of excitatory neurons that are active
        tmp = sort(Actvn, 'descend');
        el = floor(length(Rate) * Sparseness);
        if el < 1
            el = 1;
        end
        Threshold = tmp(el);
        Rate = Actvn;
        % Threshold
        % Actvn
        for neuron = 1 : length(Rate)
            if Rate(neuron) >= Threshold % threshold binary activation function
                Rate(neuron) = 1;
            else
                Rate(neuron) = 0;
            end
        end
        if display == 1
            figure(101);
            %subplot(2,2,2);
            positionVector1 = [0.2, 0.25, 0.7, 0.7]; % left, bottom, width, height. SynMat
            subplot('Position', positionVector1);
            imagesc(SynMat); title('synaptic matrix');
            colormap(gray);
            %subplot(2,2,4);
            positionVector2 = [0.2, 0.125, 0.7, 0.03]; % left, bottom, width, height. Activation
            subplot('Position', positionVector2);
            imagesc(Actvn); title('Activation'); % the internal recall term
            % disp(['Actvn max = ', num2str(max(Actvn)), ' min = ', num2str(min(Actvn))]);
            set (gca, 'YtickLabel', '');
            positionVector3 = [0.2, 0.0, 0.7, 0.03]; % left, bottom, width, height. Rate
            subplot('Position', positionVector3);
            imagesc(Rate); title('Postsynaptic Rate'); % the firing rate
            % disp(['Rate max = ', num2str(max(Rate)), ' min = ', num2str(min(Rate))]);
            set (gca, 'YtickLabel', '');
            %subplot(2,2,1);
            positionVector4 = [0.1, 0.25, 0.02, 0.7]; % left, bottom, width, height. PreSynaptic Input
            subplot('Position', positionVector4);
            imagesc(PreSynInput);   % on epoch 1 this will be the distorted recall cue
                                    % on subsequent epochs it is the Rate from the preceding epoch
            title('Presynaptic input');
            set (gca, 'XtickLabel', '');
            drawnow;
        end
        % sparse =  spars(Rate) % calc sparseness of the firing Rate representation using the function in spars.m
        %Rate
        R = corr(TrainPatts(:, patt), Rate'); % correlation coefficient
        disp(['Patt ', num2str(patt), ' epoch ', num2str(epoch), ' Correlation of recalled firing rate vector with training pattern = ', num2str(R)]);
        pause
    end
    if R > 0.98 % the criterion of correct is a correlation of 0.98 between a recalled pattern and the stored pattern
        nCorrect = nCorrect + 1; % accumulate the number of correct recalls across patterns
    end
end

PercentCorrect = nCorrect / nPatts * 100;
disp(['Percentage correct recall = ', num2str(PercentCorrect)]);

