%#										
%# function [y1,y2,h1,h2] = lmeval(topo,w1,w2,x)			
%#										
%# AIM:		Computes the output of a backpropagation neural network.		
%#										
%# PRINCIPLE:	The topology of the network is contained in the parameter matrices topo,
%#		w1 and w2. The training data contained in the matrix x are projected
%#		on the transfer functions which can be nonlinear (hyperbolic tangent)
%#		or linear.					
%#									
%# INPUT:	topo (2 x nh) : matrix containing structure of the NN (nh hidden nodes)
%#			      (H: hyperbolic tangent, L: linear, -: no function)
%#		w1 (nh x (p+1)) : matrix of random weights between input-hidden layer
%#		w2 (1 x (nh+1)) : matrix of random weights between hidden-output layer
%#		x (p x n) 	: matrix of inputs (n objects, p variables)								
%#										
%# OUTPUT:	y1 (nh x n) : outputs produced by the hidden nodes		
%#		y2 (1 x n) : outputs produced by the output node		
%#		h1 (nh x n) : inputs to hidden nodes				
%#		h2 (1 x n) : inputs to output node				
%# 										
%# SUBROUTINES: pmntanh.m : fast hyperbolic tangent function			
%#										
%# AUTHOR:	Programmed by : Magnus Norgaard, IAU/EI/IMM (1994)	
%#			        Neural Network-Based System Identification Toolbox
%#				Web site : http://www.iau.dk/Projects/proj/nnsysid.html
%#											
%#		Modified by Frederic Despagne				
%#		Copyright(c) 1997 for ChemoAC				
%#		Dienst FABI, Vrije Universiteit Brussel			
%#		Laarbeeklaan 103, 1090 Jette				
%#									
%# VERSION: 1.1 (28/02/1998)					
%#										
%# TEST:	Krzysztof Szczubialka 			
%#										

function [y1,y2,h1,h2] = lmeval(topo,w1,w2,x)

%%%%%%%%%%%%%%%%%% IDENTIFICATION OF THE NETWORK TOPOLOGY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NetDef=topo;
[n1,N] = size(x);		     	% Size of input matrix
[layers,dummy] = size(topo);        	% Number of hidden layers
L_hidden = find(NetDef(1,:)=='L')';   	% Location of linear hidden units
H_hidden = find(NetDef(1,:)=='H')';   	% Location of tanh hidden units
L_output = find(NetDef(2,:)=='L')';   	% Location of linear output units
H_output = find(NetDef(2,:)=='H')';   	% Location of tanh output units
[hidden,inputs] = size(w1);		% Size of the first weight matrix
inputs = inputs-1;			% The bias is not taken into account
y1 = zeros(hidden,N);			% Initialization of hidden layer outputs
y2 = zeros(1,N);			% Initialization of output layer outputs


%%%%%%%%%%%%%%%%%%%% COMPUTATION OF HIDDEN LAYER OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h1 = w1*[x;ones(1,N)];  			% Inputs are weighted
y1(H_hidden,:) = pmntanh(h1(H_hidden,:));	% Nonlinear nodes outputs
y1(L_hidden,:) = h1(L_hidden,:);		% Linear nodes outputs


%%%%%%%%%%%%%%%%%%%% COMPUTATION OF OUTPUT LAYER OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   
h2 = w2*[y1;ones(1,N)];				% Hidden layer outputs are weighted
y2(H_output,:) = pmntanh(h2(H_output,:));	% Nonlinear nodes outputs
y2(L_output,:) = h2(L_output,:);		% Linear nodes outputs


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%