% SOMdemo.m
% Self-Organizing Map demonstration derived from CompetitiveNetDemo.m 
% 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 self-organizing map competitive 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)

% cd c:\np\papers\visn\compcap
% SynMat(Row, Col) is SynMat(synapse, neuron)
% i.e. the synapses on a neuron form a column vector.
% missing synapses are indicated by a value of NoSyn (e.g. -1), for transparency
% compared to Matlab's sparse matrices

clear all; close all hidden; format compact; % format bank;
fig = 0;
% rng(1);         % set random number generator for code development
N = 100;        % number of neurons in the competitive net
nSyn = N;       % number of synapses on each neuron that are not deleted. Set the dilution by the divisor
nepochs = 20;   % number of training epochs. default 20
Learnrate = 0.1;   % ~ 0.1 works. Aim for weightchange = 0.02 - 0.2
if Learnrate == 0
    nepochs = 0; 
end;
OutputSparseness = 0.3; % the sparseness of the output firing
display = 1; % 0 no display; 1 display during learning epoch 1
Actvn = zeros(1, N);  % activation
Rate = zeros(1, N);  % output firing rates

FilterRange = 11; % choose an odd number to set the size of the spatial smoothing filter for the firing rates              
SpatialFilter = (1/FilterRange) * ones(1,FilterRange); % the filter kernel that will be used by conv to smooth the rate array during training

% generate training patterns overlapping by shift
leng = 20;       % length of each training pattern
InputSparseness = leng / nSyn;
shift = 3;      % shift between training patterns
nPatts = 30;
TrainPatts = zeros(N, nPatts);

for patt = 1 : nPatts
    for x = 0 : leng -1
        syn = (patt -1) * shift + x;
        syn = mod(syn, nSyn);
        TrainPatts(syn + 1 , patt) = 1;
    end    
    if (patt - 1) * shift + leng >= nSyn
        break
    end
end
nPatts = patt; % this will set nPatts to be the number that can be fitted into the array without much wrap-around
disp(['Number of patterns generated= ', num2str(nPatts), ' with shift of ', num2str(shift), ' and length of ', num2str(leng)]);
for patt = 1 : nPatts
    TrainPatts(:, patt) = TrainPatts(:, patt) / norm(TrainPatts(:, patt));     % normalize the length of each training pattern    
end

if display == 1
    fig = fig + 1;
    figure(fig)
    imagesc(TrainPatts);
    colormap(gray);
    title(['Training patterns']); xlabel('Pattern number'); ylabel('Training pattern');
    pause
end

SynMat = rand(nSyn, N); % (Synapses = Rows, neurons = Columns)
for neuron = 1 : N
    SynMat(:, neuron) = NormVeclen(SynMat(:, neuron));  % normalize the synatic weight vector on each neuron
end

% Train the synaptic matrix
MeanPreSynapticRate = sum(TrainPatts(:,1)) / nSyn;
for epoch = 1 : nepochs
    order = randperm(nPatts); % present the patterns in random order during training
    for i = 1 : nPatts
        patt = order(i);
        patt = i; % or in sequence if this is not commented out
        % calculate activation
        for neuron = 1 : N
            PreSynInput = TrainPatts(:, patt); % the retrieval cue is the conditioned stimulus
            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 or linear activation function
                % Rate(neuron) = 1; % comment out this line for a threshold linear activation function
            else
                Rate(neuron) = 0;
            end
        end
        
        Rate = conv(Rate, SpatialFilter, 'same'); % simulate positive lateral interactions between nearby neurons in the Rate array
                                                  % Older versions of Matlab may not support the 3rd argument, and may fail here.
        scale = 1.0 / (max(Rate) - min(Rate));
        Rate = (Rate - min(Rate)) * scale;            % scale Rate to 0-1.
        
        % now change the weights in the synaptic matrix
        for neuron = 1 : N
            for syn = 1 : nSyn
                preSynRate = TrainPatts(syn, patt);
                weight_change = Learnrate * Rate(neuron) * preSynRate;
                % weight_change = Learnrate * Rate(neuron) * (preSynRate - MeanPreSynapticRate); % OR try an associative rule with heterosynaptic LTD
                if weight_change > 0
                    % disp(['wtchange= ',num2str(weight_change), ' preSynrate=', num2str(preSynRate), ' Rate=', num2str(Rate(neuron)), ' neuron=', num2str(neuron), ' synapse=', num2str(syn)]); 
                end
                SynMat(syn, neuron) = SynMat(syn, neuron) + weight_change;
            end
            SynMat(:, neuron) = NormVeclen(SynMat(:, neuron)); % normalize the neuron's synaptic weight vector    
        end
        
        if display == 1
            figure(100);
            %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
    end
end
%SynMat

% Now test it
disp('Touch Enter to start Testing');
pause
RecalledRates = zeros(N, nPatts); % the firing rate output vectors have N elements
sparsenessSum = 0;
for patt = 1 : nPatts
    % calculate activation
    Rate = zeros(1, N);  % output firing rates
    for neuron = 1 : N
        PreSynInput = TrainPatts(:, patt); % the retrieval cue
        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 : N
        if Rate(neuron) >= Threshold % threshold binary or linear activation function
            % Rate(neuron) = 1; % comment out this line for a threshold linear acivation function
        else
            Rate(neuron) = 0;
        end
    end
    RecalledRates(:, patt) = Rate; % save the recalled rate for later analysis
    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;
        sparse = spars(Rate);
        sparsenessSum = sparsenessSum + sparse;
        pause;
    end

    
    
%     
%     figure(101);
%     %subplot(2,2,2);
%     positionVector1 = [0.55, 0.52, 0.4, 0.4]; % left, bottom, width, height
%     subplot('Position', positionVector1);
%     imagesc(SynMat); title('synaptic matrix');
%     colormap(gray);
%     %subplot(2,2,4);
%     positionVector2 = [0.55, 0.37, 0.4, 0.05]; % left, bottom, width, height
%     subplot('Position', positionVector2);
%     imagesc(Rate); title('Postsynaptic Rate');
%     set (gca, 'YtickLabel', '');
%     %subplot(2,2,1);
%     positionVector3 = [0.45, 0.52, 0.03, 0.4]; % left, bottom, width, height
%     subplot('Position', positionVector3);
%     imagesc(TrainPatts(:, patt));
%     title('Presynaptic input');
%     set (gca, 'XtickLabel', '');
%     subplot(2, 2, 3); 
%     [pixelCount] = hist(Rate); % firing rate histogram. 10 bins in range rate 0-1
%     bar(pixelCount, 'FaceColor', 'none');
%     grid on;
%     title('Firing Rate Distribution'); xlabel('Firing Rate x10'); ylabel('Frequency');
%     subplot(2, 2, 3);
%     hold on;
%     xdata = 1:length(pixelCount);
%     plot(xdata, pixelCount, '*');
%     %[estimates, model] = fitcurve(xdata, pixelCount);
%     %[sse, FittedCurve] = model(estimates);
%     %plot(xdata, FittedCurve, 'r');
%     title('Exponential fit to Rate Distribution'); xlabel('Firing Rate x10'); ylabel('Frequency');
%     hold off;
%     drawnow;
%     sparse = spars(Rate);
%     sparsenessSum = sparsenessSum + sparse;
%     % pause
end

MeanSparseness = sparsenessSum / nPatts;
pause

% show cross correlation matrix between the different training patterns for comparison
PattCorMat = zeros(nPatts);
for i = 1 : nPatts
    for j = 1 : nPatts
        R = corr(TrainPatts(:, i), TrainPatts(:, j));
        PattCorMat(j,i) = R;
    end
end
fig = fig + 1;
figure(fig)
x = linspace(-1, 1, 21);
imagesc(PattCorMat);
colormap(gray(256));
title([' Correlation between patterns ', '  min=', num2str(min(min(PattCorMat))), '  max=', num2str(max(max(PattCorMat)))]);
xlabel('Pattern number'); ylabel('Pattern number');

% and get the mean abs value of the off-diagonal terms, which are the
% correlations between the different patterns
numvals = 0;
PattsumCor = 0;
for col = 1 : nPatts
    for row = 1 : nPatts
        if row ~= col
            PattsumCor = PattsumCor + abs(PattCorMat(row, col));
            numvals = numvals + 1;
        end
    end
end
disp(['Average of absolute values of correlation between different input patterns=', num2str(PattsumCor / numvals)]); 

% show cross correlation matrix between the output rates for the different training patterns
CorMat = zeros(nPatts);
for i = 1 : nPatts
    for j = 1 : nPatts
        R = corr(RecalledRates(:, i), RecalledRates(:, j));
        CorMat(j,i) = R;
    end
end
fig = fig + 1;
figure(fig)
x = linspace(-1, 1, 21);
imagesc(CorMat);
colormap(gray(256));
title([' Correlation between rates ', '  min=', num2str(min(min(CorMat))), '  max=', num2str(max(max(CorMat)))]);
xlabel('Pattern number'); ylabel('Pattern number');

% and get the mean abs value of the off-diagonal terms, which are the
% correlations between the rates output for the different patterns
numvals = 0;
sumCor = 0;
for col = 1 : nPatts
    for row = 1 : nPatts
        if row ~= col
            sumCor = sumCor + abs(CorMat(row, col));
            numvals = numvals + 1;
        end
    end
end
disp(['Average of absolute values of correlation between output Rates to different patterns=', num2str(sumCor / numvals)]); 


