% PatternAssociationDemo.m
% Pattern association 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 pattern association 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 reset the random number generator for code development

N = 8;         % number of neurons in the pattern associator
nSyn = 64;      % number of synapses on each neuron
% parameters for some training patterns that overlap by shift
leng = 16;       % length of each training pattern
InputSparseness = leng / nSyn; % the sparseness of the input CS patterns is the number of bits on for any one pattern
shift = 8;      % shift between training patterns
nPatts = 8;
Learnrate = 1 / nPatts; % Not of especial significance in an autoassociation network, but this keeps the weights within a range
OutputSparseness = 1 / N; % the US patterns have one neuron on in this simple simulation 
display = 1;        % 0 no display; 1 display the network
nFlipBits = 8;     % The number of bits that are flipped to produce distorted recall cues.

% Make some training patterns
TrainPatts = zeros(nSyn, nPatts); % the conditioned stimuli
for patt = 1 : nPatts
    for x = 0 : leng -1
        syn = (patt -1) * shift + x;
        syn = mod(syn, nSyn);
        TrainPatts(syn + 1 , patt) = 1;
    end    
end
% TrainPatts
if display == 1
    fig = fig + 1;
    figure(fig)
    imagesc(TrainPatts);
    colormap(gray)
    title(['Training patterns']); xlabel('CS pattern number'); ylabel('CS 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('CS pattern number'); ylabel('CS distorted testing pattern');

    pause 
end


SynMat = zeros(nSyn, N);
Actvn = zeros(1, N);  % activation
Rate = zeros(1, N);  % output firing rates
US = zeros(nPatts, N); % Unconditioned Stimulus, also the target for retrieval
for patt = 1 : nPatts
    US(patt, patt) = 1;
end
% US    
  
% Train. See Rolls 2008 B.2.7.1 for why an associative learning rule with heterosynaptic long term depression is useful 
for patt = 1 : nPatts
    for neuron = 1 : N
        for syn = 1 : nSyn
            preSynRate = TrainPatts(syn, patt); % the CS patterns used for training
            postSynRate = US(patt, neuron);
            % weight_change = postSynRate * preSynRate; % an associative, Hebbian, synaptic modification rule
            weight_change = Learnrate * (postSynRate) * (preSynRate - InputSparseness); % 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
    % pause
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)))]);

% Test
disp('Testing');
nCorrect = 0;   % the number of correctly recalled patterns
for patt = 1 : nPatts
    % calculate activation
    for neuron = 1 : N
        % PreSynInput = TrainPatts(:, patt); % the retrieval cue is the conditioned stimulus
        PreSynInput = TrainPattsFlipped(:, patt); % OR the distorted retrieval cue to test for generalization
        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 or linear threshold 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) * OutputSparseness);
    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; % comment out this line for a threshold linear activation function
        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);
        clims = [0 1];
        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);
        clims = [0 1];
        %imagesc(Rate, clims); title('Postsynaptic Rate'); % the internal recall term
        imagesc(Rate); title('Postsynaptic Rate'); % the internal recall term
        % 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);
        title('Presynaptic input');
        set (gca, 'XtickLabel', '');
        drawnow;
        pause;
    end
    
    % check performance
    R = corr(US(:, patt), Rate'); % correlation coefficient
    disp(['Patt ', num2str(patt), ' Correlation of recalled firing rate vector with training pattern = ', num2str(R)]);
    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)]);
    