%
% Application Programming in MatLab-environment
% Autumn 2000
% 
% Exercise 9 & 10
%
% Clear the variables and close all windows
clear
close all

% ----- Problem 2 -----
n0 = 2; N = 6;
x1 = randn(n0,N); x2 = x1 + 1;
%
% output-vectors of the training data: class 1: [1 ; 0]
%                                      class 2: [0 ; 1]
y1 = [ones(1,N) ; zeros(1,N)];
y2 = [zeros(1,N) ; ones(1,N)];

%
% Thus, the final training data is the following:
%
x = [ x1 x2 ];
y = [ y1 y2 ];
n2 = size(y,1);

%
% Illustration of the data
%
p = plot(x1(1,:),x1(2,:),'bx',x2(1,:),x2(2,:),'rx'); 
set(p,'MarkerSize',12); hold on
%
clear x1 x2 y1 y2

%
% Scaling of the training data
% (unnecessary for y since its values are already in [0 1])
%
[sx,ax,bx] = scaledata(x,[ 0 1 ]);
[sy,ay,by] = scaledata(y,[ 0 1 ]);

%
% Perform the training of the network
% 
% First optimization is done without using the gradient,
% and after that with the gradient. Please compare the 
% minimum values of both optimizations.
%
n1 = 2;
k  = 1:n1; k = k';
w1_dim = n1*(n0+1);
w2_dim = n2*(n1+1);
v0 = rand(w1_dim+w2_dim,1);

% 1. without gradient
options = optimset('Diagnostics','on','Display','iter');
tic;
v = fminunc('mlp_cgd',v0,options,k,n1,sx,sy);
toc

% 2. with gradient
options = optimset('Diagnostics','on','Display','iter',...
                   'GradObj','on');
tic;
v = fminunc('mlp_cgd',v0,options,k,n1,sx,sy);
toc

% reshape the weight matrices from the vector v
w1 = reshape(v(1:w1_dim),n1,n0+1);
w2 = reshape(v(w1_dim+1:w1_dim+w2_dim),n2,n1+1);

% Let us illustrate the classification of the MLP-network.
% The system resembles the illustration of the quadratic form 
% in Exercise 2.

% We form a 100*100-grid in the current axis and scale the 
% coordinate values into the interval [0 1].
lv = get(gca,'XLim');
gridstep = (lv(2)-lv(1))/99;
xp = lv(1):gridstep:lv(2);
lv = get(gca,'YLim');
gridstep = (lv(2)-lv(1))/99;
yp = lv(1):gridstep:lv(2);
[sxp,ax,bx] = scaledata(xp,[ 0 1 ]);
[syp,ay,by] = scaledata(yp,[ 0 1 ]);
[X,Y] = meshgrid(sxp,syp);
[m1,m2] = size(X);

% We compute the output of the MLP-network in each grid point
% and then calculate the difference between the two components of
% the output (stored in Z).
Z  = zeros(m1,m2);
vp = zeros(n0,1);
o  = zeros(n2,1);
for i=1:m1
    for j=1:m2
        vp(1) = X(i,j);
        vp(2) = Y(i,j);
        o = mlp_out(vp,k,w1,w2);
        Z(i,j) = o(1) - o(2);
    end
end

% We add in the previous picture the two class boundaries given by
% the two optimization algorithms. This is achieved by searching for
% the coordinates (X,Y) for which the difference o(1)-o(2) between
% the output-components changes sign.
i1 = 2:(m1-1);
i2 = 2:(m2-1);
[I,J] = find(Z(i1-1,i2).*Z(i1+1,i2)<0 | Z(i1,i2-1).*Z(i1,i2+1)<0);
I = I + 1;
J = J + 1;
ix = I + m1*(J-1);
xc = invscaledata(X(ix),ax,bx);
yc = invscaledata(Y(ix),ay,by);
plot(xc,yc,'ro')
hold off

figure
[X,Y] = meshgrid(xp,yp);
mesh(X,Y,Z)
rotate3d


