% Simple test script that creates probabilistic Gaussian receptive fields
%
% 
% LR 19/10/2012

clear all;

% Design goals:
% 1. Layer size independence.
% 2. Torroidal layer topology.
% 3. Gaussian probabilistic connectivity.
% 4. Sparse matrix utilisation.
% 5. Layer connectivity scaling - differing sized layers should be handled
%    transparently by a scaling factor so the topology is preserved
% 5*. Pixel/Neuron independent co-ordinate frame for computing the topology
%    - this is a possible implementation of design goal 5.
% 6. Flexibility and extensibility for recurrent conenctions
% 

% Variables
numConnections = 100;
receptiveFieldWidth = 3;     % odd values give us symmetric envelopes
layerA = [28 28];
layerB = [128 128 ];

scaleFactor = layerA ./ layerB;

connectivityMatrix = spalloc(layerA(1)*layerA(2), layerB(1)*layerB(2), (numConnections+2)*layerB(1)*layerB(2)); 

fieldWidth = floor(layerA / 2);

tic;
idx = 1;
for col = 1:layerB(1)
    fprintf('.');        
    for row = 1:layerB(2)
        % fprintf('Creating (%i, %i)\n', col, row);
        % build up the connectivity row by row - i.e each receptive field
        
        tmpRow = zeros( [ layerA(1) layerA(2) ], 'double');
        connections = 0;
        randWeights = rand(1, numConnections); % define the initial weights
        
        while connections < numConnections
            randPos = mod(floor(randn(1,2) * receptiveFieldWidth) + fieldWidth, layerA) + 1; 
            if (tmpRow(randPos(1), randPos(2)) == 0)
                tmpRow(randPos(1), randPos(2)) = randWeights(connections + 1);
                connections = connections + 1;
            end
        end
        tmpRow = tmpRow ./ norm(tmpRow);
        tmpRow = circshift(tmpRow, [ (ceil(row*scaleFactor(1)) - fieldWidth(1)-1) (ceil(col*scaleFactor(2)) - fieldWidth(2)-1) ]);
        tmpRow = sparse(reshape(tmpRow, 1, layerA(1)*layerA(2)));

        connectivityMatrix(:, idx) = tmpRow;
        idx = idx + 1;
    end
end
fprintf('\nConnectivity creation time: %f sec\n', toc);
save('connectivityTest.mat', 'connectivityMatrix');


figure(1);
imagesc(reshape(connectivityMatrix(:,1), layerA(1), layerA(2)));
title('Connectivity of the 1st Neuron');

load('mnist_double01_50.mat');

testing = train_x(10,:);  % get a character from MNIST

figure(2);
imagesc(reshape(testing, layerA(1), layerA(2))');
title('Input data');
tmp = testing * connectivityMatrix;
figure(3);
imagesc(reshape(tmp, layerB(1), layerB(2))');
title('Layer Activations')

% Some code for pulling out columns, etc and stuffing them back
%for i = 1:128
%imagesc(reshape(connectivityMatrix(:, i), 28, 28));
%pause
%end

% todo: insert activation function code bsxfun with custom inline function?
%{
figure(3);
imagesc(reshape(connectivityMatrix(16*32 - 16,:), 32, 32));
%}
%{
tic
for i = 1:size(connectivityMatrix, 2)
synapseIndex = find(connectivityMatrix(:,i));
preRates = testing(synapseIndex);
connectivityMatrix(synapseIndex,i) = repmat( 50, 1, numel(synapseIndex));
end
toc
%}
%%


figure(1);
%imagesc(reshape(connectivityMatrix(,(16*16)), layerA(1), layerA(2)));

%%
clear all;

load('connectivityTest.mat');
%load('mnist_double01.mat');
load('mnist_double01_50.mat');


%%

testing = (test_x(1,:) * connectivityMatrix) ./ 100;

% hebbian update
%dWij = 0.01 * x_ij * w_ij * y_ij





figure(1);
imagesc(reshape(testing(1,:), 56, 56));
figure(2);
imagesc(reshape(test_x(1,:), 28, 28));

%% 
clear all;
hitImage = zeros(128,128);
conns = 150;

counter = 0;
while counter < conns
x = floor(randn*5) + 64;
y = floor(randn*5) + 64;
                
if (hitImage(x, y) == 0)
    hitImage(x, y) = 1;
    counter = counter + 1;
end


end

imagesc(hitImage);

%% Simple Hebbian 2D
clear all;

% Generate 3 clusters
MU1 = [-1 -1]; SIGMA1 = [.5 0; 0 .5];
MU2 = [2 2];   SIGMA2 = [.7 0; 0 .7];
MU3 = [-3 3];  SIGMA3 = [0.1 0; 0 0.1];
input = [mvnrnd(MU1,SIGMA1,100);mvnrnd(MU2,SIGMA2,100);mvnrnd(MU3,SIGMA3,100)];
clear MU1 MU2 MU3 SIGMA1 SIGMA2 SIGMA3


input = input + abs(min(min(input)));
input = input ./ abs(max(max(input)));

for i = 1:size(input, 1)
    %input(i,:) = input(i,:) ./ norm(input(i,:));
end

connectivity = rand(2,3);

for i = 1:size(connectivity,2)
    %connectivity(:,i) = connectivity(:,i) ./ norm(connectivity(:,i));
end

figure(1);
plot(connectivity(1,:),connectivity(2,:),'LineStyle','o', 'color', 'g');



%input = connectivity';

hold on;
plot(input(:,1),input(:,2),'LineStyle','.');

epochs = 30;

deltaSizes = [];
for kk = 1:epochs

    sum = 0;
    for j = 1%1:size(input,1)
        
        y = input(j,:) * connectivity;
        %tmp = zeros(size(y));
        %tmp(y>median(y)) = 1;
        %y = tmp;
        
        %y = y ./ norm(y);
        
        dW = 0.01 * (connectivity .* (input(j,:)' * y));
        sum = sum + dW;
        
        connectivity = connectivity + dW;
        
        for i = 1:size(connectivity,2)
         %   connectivity(:,i) = connectivity(:,i) ./ norm(connectivity(:,i));
        end
        
        deltaSizes(epochs) = norm(sum);
        
    end
    hold on;
    plot(connectivity(1,:),connectivity(2,:),'LineStyle','+', 'color', 'r');
    
end
    
    hold on;
plot(connectivity(1,:),connectivity(2,:),'LineStyle','+', 'color', 'r');
hold off;

figure(2);
plot(deltaSizes);

%%



clear all;

% Generate 3 clusters
MU1 = [-1 -1]; SIGMA1 = [.5 0; 0 .5];
MU2 = [2 2];   SIGMA2 = [.7 0; 0 .7];
MU3 = [-3 3];  SIGMA3 = [0.1 0; 0 0.1];
input = [mvnrnd(MU1,SIGMA1,100);mvnrnd(MU2,SIGMA2,100);mvnrnd(MU3,SIGMA3,100)];
clear MU1 MU2 MU3 SIGMA1 SIGMA2 SIGMA3





net = competlayer(8, 0.01, 0.01);
done = train(net, input');


figure(1);
plot(input(:,1),input(:,2),'LineStyle','.');
W = cell2mat(done.IW);

hold on;
plot(W(:,1),W(:,2),'LineStyle','+', 'color', 'r');
hold off;





%% 
train_x = train_x(1:50,:);
save('mnist_double01_50.mat', 'train_x');

