%#											
%# function [w1f,w2f,tablein,tableout] = nnmodel(topo,xtrain,ytrain,xmon,ymon,xtest,
%#	   					ytest,epochs,trials,setseed)	
%#										
%# AIM:	Backpropagation neural network with one hidden layer for multivariate
%#		calibration. (Designed to model only one response y at a time)	
%#										
%# PRINCIPLE:	The neural network is trained with the Levenberg-Marquardt algorithm.
%#		The activation functions can be either linear ('L') or hyperbolic
%#		tangent ('H'). The network has one hidden layer. The network topology 
%#		is determined by the matrix 'topo', which has two rows : one for
%#		specifying the hidden layer and one for the output layer.	
%#		Example: topo = ['HHHH';'L---'] defines a network with 4 nodes in the 
%#		hidden layer (hyperbolic tangent transfer functions) and one node in
%#		the output layer (linear transfer function). A weight can be pruned by
%#		setting it to zero. The bias is automatically included as the last
%#		column of the weight matrices. The model is built with the training
%#		set, but a monitoring set is necessary to avoid overtraining.
%#									
%#		At the end of training, a menu is available to summarize results 
%#		(training, monitoring and test if available). Statistical estimations
%#		are obtained using the solution with monitoring error closest to the 
%#		median (over all replicate training runs performed).
%#										
%#		A diagnostic is performed to assess the sensitivity (influence) of each 
%#		input variable and hidden node in the model. Maps of projections of the
%# 		training samples on hidden nodes allow to visualize clusters and 
%#		outliers. Deviations from linearity are calculated for each hidden node	
%#		in order to estimate the degree of nonlinearity of a data set.	
%#		Input variables with lowest sensitivities should be removed and the 
%#		neural network re-trained. 					 
%#									
%# REFERENCE:	"Tutorial Review: Neural Networks in Multivariate Calibration"	
%#		         F.Despagne, D.L.Massart, Analyst 123 (1998) 157R-178R		
%# 										
%# INPUT:	topo (2 x nh) : Matrix containing structure of the NN (nh hidden nodes)
%#			        (H: hyperbolic tangent, L: linear, -: no node)	
%#		    xtrain (ntr x p) : Matrix of training inputs (ntr objects, p variables)
%#		    ytrain (ntr x 1) : Vector of training responses			
%#		    xmon   (nm x p) : Matrix of monitoring inputs (nm objects, p variables)
%#		    ymon   (nm x 1) : Vector of monitoring responses		
%#		    xtest  (nte x p) : Matrix of test inputs (if not available, enter [])
%#		    ytest  (nte x 1) : Vector of test responses (if not available, enter [])
%#		    epochs : Maximum number of epochs (iterations) for the training	
%#		    trials : Number of trials. If trials>1, several trials are performed
%#			   with the same topology but different initial random weights.
%#			   In this case, average results are reported.		
%#		    setseed : Seed for generating initial random weights.	
%#			  (Optional. Only valid if trials=1. Default: random seed)							
%#									
%# OUTPUT:	w1f (nh x (p+1)) : matrix of weights between input and hidden layer
%#              		   corresponding to the trial with monitoring error closest to the median.
%#				   		
%#		    w2f (1 x (nh+1)) : matrix of weights between hidden and output layer 
%#              		   corresponding to the trial with monitoring error closest to the median.
%#				   			
%#		    tablein (mxtr x 2):input range-scaling parameters.	
%#		    tableout (1 x 2):  output range-scaling parameters.	
%# 										
%# SUBROUTINES: range.m : range-scaling of training, monitoring and test data
%#		        levmarq.m : training of the network with Levenberg-Marquardt	
%#		        lmeval.m : estimation of responses with the final neural network model
%#		        invrange.m : returns range-scaled data to original scale	
%#		        rms.m : calculates root mean squared error		
%#		        corrplot.m : draws correlation plot			
%#		        drawnn.m : draws the neural network structure		
%#		        linfit.m : builds univariate linear regression model		
%#		        kolnorm.m : Kolmogorov-Smirnov test to check if replicate y-predicted 
%#			              values are normally distributed.		
%#										
%# AUTHOR:	Frederic Despagne					
%#		Copyright(c) 1997 for ChemoAC					
%#		Dienst FABI, Vrije Universiteit Brussel			
%#		Laarbeeklaan 103, 1090 Jette			
%#									
%# VERSION: 1.1 (28/01/1999)				
%#									
%# TEST:									
%#										

function [w1f,w2f,tablein,tableout] = nnmodel(topo,xtrain,ytrain,xmon,ymon,xtest,ytest,epochs,trials,setseed)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INITIALIZATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Test to check if all necessary input variables are available
if nargin < 9		
	error('Some input arguments are missing.')
end

% Test to check that only one trial is performed with the defined random number generator seed
if nargin == 10
	if trials ~= 1
		error('The "setseed" input variable is only valid if "trials" is set to 1.')
	end
end

% Test to check if a test set is available
if isempty(xtest)
	testflag = 0;		
else
	testflag = 1;
end  

[nxtr,mxtr] = size(xtrain);	   % Size of the training input data
[ndef,mdef] = size(topo);	   % Topology of the network
nxm = size(xmon,1);		       % Number of monitoring samples
nxte = size(xtest,1);		   % Number of test samples
init = 1/(1*mxtr);  		   % Determination of weights initialization range
periodisp = ceil(epochs/10);   % Determination of the periodicity of display
indisp = zeros(1,trials);	% Vector containing optimum number of epochs, initialization
partyh = zeros(mdef,nxtr);	% Matrix of partial responses with one hidden node, training set
w1in = [];
w2in = [];
w1fin = [];			% The matrix where replicate sets of weights input-hidden will be stored
w2fin = [];			% The matrix where replicate sets of weights hidden-output will be stored
conctrx = [];			% Matrix of partial responses with one input variable, training set, first iteration 
concytrhat = [];		% Matrix of replicate predicted y-values with partial models, training set
setfig = [312 288 560 420];	% Vector of position indices for figure display
periodisp = [];


%%%%%%%%%%%%%%%%%%%%%%%%%% SCALING OF INPUT AND OUTPUT DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
for i = 1:mxtr

	% Range-scaling of training and monitoring input data between -1 and 1
	[xtra(:,i),tablein(i,:),xmonit(:,i)] = range(xtrain(:,i),-1,1,xmon(:,i));

	% Range-scaling of test input data between -1 and 1	
	if testflag		
		[xtra(:,i),tablein(i,:),xtes(:,i)] = range(xtrain(:,i),-1,1,xtest(:,i));
	end
end

% Range-scaling of training and monitoring output data between 0.2 and 0.8
[ytra,tableout,ymonit] = range(ytrain,0.2,0.8,ymon);

xtr = xtra'; ytr = ytra';	% Transposition of training input and output data	
xm = xmonit'; ym = ymonit';	% Transposition of monitoring input and output data


if testflag

	% Range-scaling of test output data between 0.2 and 0.8
	[ytra,tableout,ytes] = range(ytrain,0.2,0.8,ytest);

	xte = xtes'; yte = ytes';	% Transposition of test input and output data
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NETWORK INITIALIZATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for nit = 1:trials			% Loop on the number of trials
	close
	clc
	disp(' ')

	% Set the random number generator seed
	if nargin == 10
		randseed = setseed;
	else
		t0 = clock;			% The clock parameters vector
		randseed = 210610+100*nit;	% sum(100*t0(:,1:6)) A number that changes every second, as a seed to generate random weights
	end
	rseed(nit) = randseed;
	rand('seed',randseed);			% Random generator seed set on the clock

	w10 = rand(mdef,mxtr+1);		% Generation of the first matrix of random weights
	w20 = rand(1,mdef+1);			% Generation of the second matrix of random weights
	[w1,w1tab] = range(w10,-init,+init);	% Weights scaled to the initialization range
	[w2,w2tab] = range(w20,-init,+init);	% Weights scaled to the initialization range


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRAINING AND PREDICTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	% Training of the network. The monitoring set is used to stop the training
	[w1f,w2f,indisp(nit)] = levmarq(topo,w1,w2,xtr,ytr,xm,ym,epochs,periodisp);

	w1in = [w1in;w1];	% The replicate initial w1 matrices are accumulated
	w2in = [w2in;w2];	% The replicate initial w2 matrices are accumulated
	w1fin = [w1fin;w1f];	% The replicate final w1 matrices are accumulated
	w2fin = [w2fin;w2f];	% The replicate final w2 matrices are accumulated

	% Estimation of training responses
	[ytr1,ytr2,intr1,intr2] = lmeval(topo,w1f,w2f,xtr);

	% Inverse-scaling of estimated training responses
	ytrhat(nit,:) = invrange(ytr2,0.2,0.8,tableout);

	% Determination of root mean squared error of calibration
	rmstra(nit) = rms(ytrain,ytrhat(nit,:)');

	% Estimation of monitoring responses
	[ymon1,ymon2,inmon1,inmon2] = lmeval(topo,w1f,w2f,xm);

	% Inverse-scaling of estimated monitoring responses
	ymonhat(nit,:) = invrange(ymon2,0.2,0.8,tableout);

	% Determination of root mean squared error of prediction (monitoring set)
	rmsmon(nit) = rms(ymon,ymonhat(nit,:)');

	if testflag

		% Estimation of the test responses
		[ytes1,ytes2,intes1,intes2] = lmeval(topo,w1f,w2f,xte);

		% Inverse-scaling of estimated test responses
		yteshat(nit,:) = invrange(ytes2,0.2,0.8,tableout);

		% Determination of root mean squared error of prediction (test set)
		rmstest(nit) = rms(ytest,yteshat(nit,:)');
	end

end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%% DETERMINATION OF THE REFERENCE TRIAL (MEDIAN RMSEM) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

medrmsmon = median(rmsmon);			% Determination of the RMSEM median
diffmedian = abs(rmsmon-medrmsmon);		% All RMSEM values are compared to the median
[indiffref,sortind] = sort(diffmedian);		% The RMSEM values are sorted according to their distance to the median
indref = sortind(1);				% The index of the reference trial (closest to the median) is found		
retainw1 = 1+(indref-1)*mdef;			% Indices of lines corresponding to the reference trial in w1fin
w1inref = w1in(retainw1:retainw1+mdef-1,:);	% The initial w1 weight matrix for the reference trial
w2inref = w2in(indref,:);			% The initial w2 weight matrix for the reference trial
w1ref = w1fin(retainw1:retainw1+mdef-1,:);	% The final w1 weight matrix for the reference trial
w2ref = w2fin(indref,:);			% The final w2 weight matrix for the reference trial
clear w1fin w2fin w1in w2in

% Estimation of training responses
[ytr1,ytr2,intr1,intr2] = lmeval(topo,w1ref,w2ref,xtr);

% Inverse-scaling of estimated training responses
ytrref = invrange(ytr2,0.2,0.8,tableout);

% Determination of root mean squared error of calibration
rmstraref = rms(ytrain,ytrref');

% Estimation of monitoring responses
[ymon1,ymon2,inmon1,inmon2] = lmeval(topo,w1ref,w2ref,xm);

% Inverse-scaling of estimated monitoring responses
ymonref = invrange(ymon2,0.2,0.8,tableout);

% Determination of root mean squared error of prediction (monitoring set)
rmsmonref = rms(ymon,ymonref');

if testflag

	% Estimation of the test responses
	[ytes1,ytes2,intes1,intes2] = lmeval(topo,w1ref,w2ref,xte);

	% Inverse-scaling of estimated test responses
	ytesref = invrange(ytes2,0.2,0.8,tableout);

	% Determination of root mean squared error of prediction (test set)
	rmstestref = rms(ytest,ytesref');
end

residtrain = ytrain-ytrref';		% training set residuals for the reference trial
residmon = ymon-ymonref';		% Monitoring set residuals for the reference trial
mrmstra = mean(rmstra); 		% Determination of average RMSEC on retained trials	
stdtra = std(rmstra);			% Standard deviation of RMSEC on retained trials
mrmsmon = mean(rmsmon);			% Determination of average RMSEM on retained trials
stdmon = std(rmsmon);			% Standard deviation of RMSEM on retained trials

if testflag
	residtest = ytest-ytesref';	% Test set residuals for the reference trial	
	mrmstest = mean(rmstest);	% Determination of average RMSEP on retained trials
	stdtest = std(rmstest);		% Standard deviation of RMSEP on retained trials
end


%%%%%%% ESTIMATION OF INPUT VARIABLES SENSITIVITIES BY PARTIAL MODELING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:mxtr
	xtrunc = 0*xtr;		% All training input variables but one
	xtrunc(i,:) = xtr(i,:);	% are set to 0

	% Training of the network. The monitoring set is used to stop the training
	[w1n,w2n] = levmarq(topo,w1inref,w2inref,xtr,ytr,xm,ym,epochs,[]);

	% Estimation of training responses using only one input in the model
	[ytrunc1,ytrunc2] = lmeval(topo,w1n,w2n,xtrunc);
	
	% Inverse-scaling of training responses estimated using one input only
	ytrunchat = invrange(ytrunc2,0.2,0.8,tableout);
	concyhtr(i,:) = ytrunchat;
		
	
	xtruncm = 0*xm;		% All monitoring input variables but one
	xtruncm(i,:) =xm(i,:);	% are set to 0

	% Estimation of monitoring responses using only one input in the model
	[ytruncm1,ytruncm2] = lmeval(topo,w1n,w2n,xtruncm);

	% Inverse-scaling of monitoring responses estimated using one input only
	ytruncmhat = invrange(ytruncm2,0.2,0.8,tableout);
	concymhat(i,:) = ytruncmhat;

	% Determination of variance of the predicted training responses without noise
	varytrunc0(1,i) = cov(ytrunchat);

	if testflag
		xtrunct = 0*xte;	% All test input variables but one
		xtrunct(i,:) = xte(i,:);% are set to 0	

		% Estimation of test responses using only one input in the model
		[ytrunct1,ytrunct2] = lmeval(topo,w1n,w2n,xtrunct);
		
		% Inverse-scaling of test responses estimated using one input only
		ytruncthat = invrange(ytrunct2,0.2,0.8,tableout);
		concyteshat(i,:) = ytruncthat;
	end
end

% The sensitivity of each input variable is the scaled variance of the training
% responses estimated using this variable only
varytrunc = 100*(varytrunc0./sum(varytrunc0));			


%%%%%%% ESTIMATION OF HIDDEN NODES SENSITIVITIES BY PARTIAL MODELING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 			
for i = 1:mdef
	w2trunc = 0*w2ref;				% All hidden nodes but one