% NNet_toolbox.m
%
% creates a user-defined 3-layer neural network
% with read-out
%
% Written for CoSMO 2012 by Gunnar Blohm
% www.compneurosci.com 
%
%=============================================

clear;
warning off;

%% creating training set
%----------------------
N = 200; % length of training set
input = 15*(randn(2,N)); % random set of retinal target and eye positions
output = input(1,:) - input(2,:); % spatial target position

% network variables
x = -50:2:50; % preferred directions of input units
Ni = length(x); % number input units (map)
Ne = 2; % number eye position coding input units
Nh = 91; % number hidden layer units
No = length(x); % number population output units

% network input
Pd{1,1} = exp(-(repmat(x',1,N)-repmat(input(1,:),Ni,1)).^2./10.^2./2); % activations of 1-D retinal map
Pd{2,1} = [(input(2,:)+50)/100; (input(2,:)-50)/100]; % push-pull eye position coding
Ai = [];
Q = 0;
TS = 1;
VV = [];
TV = [];

% network output
Tl{1,1} = exp(-(repmat(x',1,N)-repmat(input(1,:) - input(2,:),Ni,1)).^2./10.^2./2); % 1-D spatial output map

%% general network definition
%---------------------------
net = network; % create a network

net.numInputs = length(Pd); % number of physically different inputs
net.numLayers = 2; % number of layers (not including the input)
net.biasConnect = [0; 0]; % bias enable of layers (0: no, 1: yes)
net.inputConnect = [ones(1, net.numInputs); zeros(1, net.numInputs)]; 
% inputs are connected to layer indicated by "1"
net.layerConnect = [0 0; 1 0]; % interconnection of layers (row: origin
                               %                            line: target)
net.outputConnect = [0 1]; % outputs are layers indicated by "1" 

% different potential training algorithms
net.trainFcn = 'trainrp'; % resilient back-prop
% net.trainFcn = 'trainscg'; % scaled conjugate gradient back-prop
% net.trainFcn = 'trainoss'; % one step secant back-prop
% net.trainFcn = 'traincgb'; % Powell-Beale conjugate gradient back-prop
% net.trainFcn = 'traingdx'; % gradient descent with momentum and adaptive learning back-prop

net.adaptFcn = 'trainb'; % batch training
net.performFcn = 'mse'; % error function for training performance monitoring

% specific network parameters
scD = 1;
net.inputs{1}.size = Ni; % target direction
net.inputs{1}.range = repmat([-scD scD], net.inputs{1}.size, 1);
net.inputs{2}.size = Ne; % eye position coding
net.inputs{2}.range = repmat([-scD scD], net.inputs{2}.size, 1);
net.layers{1}.size = Nh; % number of hidden layer units
net.layers{1}.transferFcn = 'logsig';
net.layers{2}.size = No; % number of population output units
net.layers{2}.transferFcn = 'logsig';
net.trainParam.goal = 1e-3; % maximum training steps
net.trainparam.epochs = 1e4; % desired error
net.trainParam.show = 1; % shows training progress at ever step
net.trainParam.min_grad = 1e-7; % minimum allowed gradient

% initialize network (i.e. weights)
net.initFcn = 'initlay';
for j = 1:net.numInputs,
    net.inputWeights{1,j}.initFcn = 'rands'; % random initial input weights
    net.inputWeights{1,j}.delays = 0;
end
for j = 1:net.numLayers,
    net.layers{j}.initFcn = 'initwb';
%     net.biases{j}.initFcn = 'rands';
    for k = 1:net.numLayers,
        net.layerWeights{j,k}.initFcn = 'rands'; % random initial layer weights
        net.layerWeights{j,k}.delays = 0;
    end
end
net = init(net);

%% traing the net
%---------------
tic;

[net,TR] = train(net,Pd,Tl); % batch training 
% (use 'adapt' for sequential training)
disp(['Training with ' net.trainFcn ' needed '...
    num2str(toc,'%10.3f') ' s'])

%% simulate and save network
%--------------------------
Ntest = 1000; % length of test set
Tinput = 15*(randn(2,Ntest));

% test input
Pdtest{1,1} = exp(-(repmat(x',1,Ntest)-repmat(Tinput(1,:),Ni,1)).^2./10.^2./2);
Pdtest{2,1} = [(Tinput(2,:)+50)/100; (Tinput(2,:)-50)/100];

% test output
Tltest{1,1} = exp(-(repmat(x',1,Ntest)-repmat(Tinput(1,:) - Tinput(2,:),Ni,1)).^2./10.^2./2);

% test net
net.outputConnect = [1 1]; % change outputs so that we can also simulate hidden layer activity
Z = sim(net,Pdtest); % simulate trained network
figure
OUTdesired = Tinput(1,:) - Tinput(2,:);
OUTdecoded = Z{2}'*x'./sum(Z{2},1)'; % centre of mass decoding
[m, b, r] = postreg(OUTdecoded', OUTdesired);

net.userdata.epochs = TR; % save training progress in userdata structure of network
save('net')

%% analyze net
%-------------