%----------	uses MATLAB Neural Network Toolbox
%----------	Version 2.0.4 02-Jun-1997

trial		= 1;
lay		= [3 3];		% # hidden, output layers neurons
eta		= 0.001;		% learning rate
winsize		= 10;			% epoch window size for stopping
maxepochs	= 200;			% maximum epochs
tol		= 0.005;		% stop if change in error < .05%

f1		= 'tansig';		% hidden layer activation function
f2		= 'logsig';		% output layer activation function

trainfraction	= .5;			% fraction of total as training set
msegoal		= .01;			% MSE goal

filename1	= sprintf('trainbp%d_%.4f_%da.ps',lay(1),eta,trial);
filename2	= sprintf('trainbp%d_%.4f_%db.ps',lay(1),eta,trial);
%----------	load data
if ~exist('build');
	load build.dat;
end
if ~exist('bo');
	load bo.dat;
end
d	= 14;				% # dimensions
dy	= 3;				% # outputs

%----------	load ordered data (cheating)
io	= bo(:,1:d)';
oo	= bo(:,d+1:d+dy)';

%----------	all training data
inpt	= build(:,1:d)';
oupt	= build(:,d+1:d+dy)';
[d N]	= size(inpt);

trange	= 1:round(trainfraction*N);	% training set index
tsize	= length(trange);

vrange	= (max(trange)+1):N;		% validation set index
vsize	= length(vrange);

tin	= inpt(:,trange);		% training i/o
tout	= oupt(:,trange);

vin	= inpt(:,vrange);		% validation i/o
vout	= oupt(:,vrange);

range	= [min(build(:,1:d));max(build(:,1:d))]';

%----------	init network
%rand('state',1);

[w1,b1,w2,b2]	= initff(range,lay(1),f1,lay(2),f2);

%----------	train network
update	= 2;
maxepoch= 1;
ssegoal	= msegoal * tsize;
tp	= [update maxepoch ssegoal eta];
%tp	= [update maxepoch ssegoal];
%tp(5)= eta;

clear vmse;
clear tmse;
%-----	Run 1st winsize epochs without termination
for epochs=1:winsize
	epochs
	[w1,b1,w2,b2,te,tr]	= trainbp(w1,b1,f1,w2,b2,f2,tin,tout,tp);
%	[w1,b1,w2,b2,te,tr]	= trainlm(w1,b1,f1,w2,b2,f2,tin,tout,tp);
	tsse			= tr(2);
	tmse(epochs)		= tsse/tsize;
	
	vouthat			= simuff(vin,w1,b1,f1,w2,b2,f2);
	verr			= vouthat-vout;
	vserr			= verr.*verr;
	vsse			= sum(sum(vserr));
	vmse(epochs)		= vsse/vsize;		
end

oldavg_vmse	= mean(vmse((epochs-winsize+1):epochs));


%-----	Terminate if winsize average increases or vmse change < tol
for epochs=winsize+1:maxepochs
	epochs
	[w1,b1,w2,b2,te,tr]	= trainbp(w1,b1,f1,w2,b2,f2,tin,tout,tp);
%	[w1,b1,w2,b2,te,tr]	= trainlm(w1,b1,f1,w2,b2,f2,tin,tout,tp);
	tsse			= tr(2);
	tmse(epochs)		= tsse/tsize;
	
	vouthat			= simuff(vin,w1,b1,f1,w2,b2,f2);
	verr			= vouthat-vout;
	vserr			= verr.*verr;
	vsse			= sum(sum(vserr));
	vmse(epochs)		= vsse/vsize;		

	avg_vmse		= mean(vmse((epochs-winsize+1):epochs));
	if avg_vmse>oldavg_vmse
		disp('Average VMSE increased, exiting');
		break
	end;
	oldavg_vmse		= avg_vmse;

	delta_vmse		= diff(vmse(epochs:-1:epochs-1))/vmse(epochs-1)
	if abs(delta_vmse)<tol
		disp('Change in VMSE < tol, exiting');
		break
	end
end

figure(1);
clf;
subplot(2,2,1);
plot(1:epochs,tmse(1:epochs),'r-',1:epochs,vmse(1:epochs),'b--');
legend('Train MSE','Validation MSE');
xlabel('Epochs');
ylabel('MSE');
ts	= sprintf('trainbp 14-%d-3 : lr=%0.4f MSE_{T} = %0.4f MSE_{V} = %0.4f',lay(1),eta,tmse(epochs),vmse(epochs));
title(ts);

fstr	= ['print -dpsc ' filename1];
eval(fstr);
%-----------	plot true vs predicted
yhat	= simuff(inpt,w1,b1,f1,w2,b2,f2);
figure(1);
for y	= 1:dy
	subplot(2,2,y+1);
	plot(oupt(y,trange),yhat(y,trange),'rx',oupt(y,vrange),yhat(y,vrange),'b.',[0 1],[0,1],'k--');
	axis([0 1 0 1]);
	ylabel('True');
	xlabel('NN');
	axis square;
	tstring	= sprintf('(%d) True vs NN output',y);
	title(tstring);
end

fstr	= ['print -dpsc ' filename1];
eval(fstr);
%-----	plot some of the outputs w.r.t. ordered outputs
yo	= simuff(io,w1,b1,f1,w2,b2,f2);
prange	= 1:500;
figure(2);
clf;
for y	= 1:dy
	subplot(dy,1,y);
	plot(oo(y,prange),'r-');
	hold on;
	plot(yo(y,prange),'b--');
	xlabel('Sample');
	ystring	= sprintf('Output %d',y);
	ylabel(ystring);
	%legend('Original','NN');
end
subplot(dy,1,1);
title('- True -- NN');
fstr	= ['print -dpsc ' filename2];
eval(fstr);