% NEWTC new time-course neural network model
%    [net,P,T] = newtc(p,t,D,S,TF)
%       p : input array
%       t : target array (time-course data)
%       D : layer delays
%       S : layer sizes
%      TF : transfer functions
%     net : feedforward network
%       P : sequential inputs for net
%       T : sequential outputs for net
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [net,P,T] = newmltc(p,t,D,S,TF)

l = size(p,1); % number of inputs
n = size(p,2); % number of observations
q = size(t,1); % number of time points
d = length(D); % number of delay lines

% Create extended input array
P = zeros(l+d,n*q);
for i = 1:n
   % duplicate inputs for all time points
   ps = p(:,i)*ones(1,q);
   P(1:l,(1:q)+((i-1)*q)) = ps;
end

% Sequentialize time-course data
T = reshape(t,1,n*q);

% Create additional inputs that represent time-delays
P((l+1):(l+d),:) = repmat(reshape(t,1,n*q),d,1);
for i = 1:d
   P(i+l,:) = circshift(P(i+l,:),[0 D(i)]);
   for j = 0:(n-1)
      P(i+l,(j*q+1):(j*q+D(i))) = 0;
   end
end


% Create feedforward network. The first two layers
% split the inputs
TF = cat(2,{'purelin','purelin'},TF);
net = newff(minmax(P),[d l S],TF);
net.inputConnect = zeros(size(net.inputConnect));
net.inputConnect(1) = 1;
net.inputConnect(2) = 1;
net.biasConnect = zeros(size(net.biasConnect));
%net.biasConnect(1) = 0;
%net.biasConnect(2) = 0;
%net.biasConnect(net.numLayers) = 0;

% Connect delay inputs to first layer and turn off
% learning
net.IW{1} = cat(2,zeros(d,l),eye(d,d));
net.inputWeights{1}.learn = 0;

% Connect real inputs to second layer and turn off
% learning
net.IW{2} = cat(2,eye(l),zeros(l,d));
net.inputWeights{2}.learn = 0;

% Disconnect all layers
net.layerConnect = zeros(size(net.layerConnect));

dl = net.numLayers;

% Connect feedforward layers
for i = 3:dl
   net.layerConnect(i,i-1) = 1;
   net.layerWeights{i,i-1}.learnFcn = 'learngdm';
   net.LW{i,i-1} = rands(size(net.LW{i,i-1},1),...
                         size(net.LW{i,i-1},2));
end

% Connect delay inputs to all layers
for i = 3:dl-1
   net.layerConnect(i,1) = 1;
   net.layerWeights{i,1}.learnFcn = 'learngdm';  
   net.LW{i,1} = rands(size(net.LW{i,1},1),...
                       size(net.LW{i,1},2));   
end