%%%% 2 classes BP algorithm in non-linear-classfing


% Produce sample set and original classfying line
% Produce samples: R0,R1
close all
clear all

% neuro number in each layers
N=2;        % Input vector dimension                      
K=2;        % First layer neuros
J=1;        % Second layer neuros


M=5000;                         % Maxmum allowably iterative times
s=30;                              
S=3*s;                          % Nomber of total samples. Learning set, with total 2*S samples, is composed of both S samples class  
alf=0.1;                       % Step length
u0=0.25;                        % Sigmoid function f=1./(1+exp(-u/u0))
clr=['r  ','g  ','y  ','b  ','c  ','m  ','r: ','g: ','y: ','b: ','c: ','m: ','r-.','g-.','y-.','b-.','c-.','m-.'];

            figure(1)
            hold on
            whitebg(1,'k')
            figure(2)
            hold on
            whitebg(2,'k')
                        
        % Sigmoid function
            figure(3)
            u=-10:0.1:10;
            plot(u,1./(1+exp(-u./u0)));
            clear u
            

        % Produce and plot sample set X0 and X1, learning samiple set R=R0+R1
            X0=randn(S,N)+1.2*ones(S,N);                        % S X0 input vector, in line vector with N dimension
            X0=[X0,ones(S,1)];                                  % Expand the pattern vector X0 with an 1

            X11=randn(s,N)-[2*ones(s,1),-1.5*ones(s,1)];        %
            X12=randn(s,N)-[-1.5*ones(s,1),2*ones(s,1)];        %  
            X13=randn(s,N)-[1.5*ones(s,1),1.5*ones(s,1)];       % X1 composed of 3 part, each in a scale of s and 3*s=S
            X1=[X11;X12;X13];                                   %
            X1=[X1,ones(S,1)];                                  % Expand the pattern vector X1 with an 1

            D0=0.9*ones(S,1);               % Expected output for X0; 
            D1=0.1*ones(S,1);               % Expected output for X1;
            
            mX=[X0;X1];             % Sample matrix with dimension S*(N+1)
            D=[D0;D1];              % Expected output matrix with dimension of S*I
            
            figure(1)
            plot(X0(:,1),X0(:,2),'b*')
            plot(X1(:,1),X1(:,2),'gx')

        % Plot axies
                        Xmax=max(mX(:,1));
                        Xmin=min(mX(:,1));
                        Ymax=max(mX(:,2));
                        Ymin=min(mX(:,2));
                        Xmax=ceil(Xmax);
                        Xmin=floor(Xmin);
                        Ymax=ceil(Ymax);
                        Ymin=floor(Ymin);
            plot([Xmin,Xmax],[0,0],'m-.')
            plot([0,0],[Ymin,Ymax],'m-.')
            text(Xmin+0.3,Ymax-0.5,'X0: *');
            text(Xmin+0.3,Ymax-0.9,'X1: x');
            
                        
            clear X11 X12 X13 nn X0 X1 D0 D1 s
            
   aw=1;
   Vn=0;    % Times
   
   while aw~=0,
   aw=input('Continue(1) or stop (any key)?')
    if aw==0,
        break;
    end
            
        %Produce original matrixes W3,W2,W1
            W1=rand(K,N+1);
            W2=rand(J,K+1);
            p=1;
            
            while p<M ,
                Ep(p)=0;                                            % Square error
                dW1=0;dW2=0;                                        % Weights correctors
                
                    for s=1:2*S
 
                        X=mX(s,:);
                        
                        % out of first layer
                        Net1=X*W1';                         % Pure input, K dimension vector;
                        u=exp(-Net1./u0);                   % K'd
                        O1=1./(1+u);                        % Output vectors of layer 1
                        O1=[O1,1];                          % Expand the output vectors
                        F1=diag(u./(u0*(1+u).^2));          % Differential coefficient, K by K dimension matrix
                
                        % out of second layer
                        Net2=O1*W2';                        % Pure input, J dimension vector;
                        u=exp(-Net2./u0);                   % J'd
                        O2=1./(1+u);                        % Output vectors of layer 2
                        F2=diag(u./(u0*(1+u).^2));          % Differential coefficient, J by J dimension matrix
                        
                        % calculate errors and amend weights
                        Es=D(s)-O2;                         % Errors
                        Ep(p)=Ep(p)+Es*Es';                 % Summation of square error
                        
                        d2=Es*F2;                           % Error component for layer 2
                        d1=d2*W2(:,1:K)*F1;                 % Error component for layer 1
                      
                        dW2=dW2+alf*d2'*O1;  dW1=dW1+alf*d1'*X;                    %
                        
                    end % of for
                        W1=W1+dW1/(2*S);W2=W2+dW2/(2*S);    % Amending with mean correctors
                        p=p+1;
                        if mod(p,100)==0,
                            disp(p)
                        end
                
             end % while
             
                figure(2)                   % Plot error evaluation line
                cln=3*(mod(Vn,18))+1;       % 
                clrnum=clr(cln:cln+2);      % Color number 'r  ' or 'y: ' or 'g-.' etc. each with 3 characters
                plot(Ep,clrnum);
                clear p s u dW1 dW2 cln
                
                % Plot the borderline of 2 classes
                figure(1)
                n=1; Lx=[];Ly=[];
                m=1; Rx=[];Ry=[];
                for y=Ymin:0.1:Ymax
                    NL='New Line';               % A new line start
                    for x=Xmin:0.1:Xmax
                        X=[x,y,1];
                        O1=[1./(1+exp(-X*W1'./u0)),1];
                        O2=1/(1+exp(-O1*W2'/u0));
                        
                        if NL~='New Line'    
                            if (pr<=0.5)&(O2>=0.5),
                                Lx(n)=x;
                                Ly(n)=y;
                                n=n+1;
                            end
                            if (pr>=0.5)&(O2<=0.5),
                                Rx(m)=x;
                                Ry(m)=y;
                                m=m+1;
                            end
                        end
                        pr=O2; NL='0';          % Storage last value
                    end
                end
                if size(Lx)>0,
                    plot(Lx,Ly,clrnum);             % Plot borderline
                    text(Lx(1),Ly(1),num2str(Vn));                          % 
                    text(Lx(fix(length(Lx)/2)),Ly(fix(length(Lx)/2)),num2str(Vn));    % Write bordline number
                    text(Lx(length(Lx)),Ly(length(Lx)),num2str(Vn));        % Repeat
                end
                if size(Rx)>0,    
                    plot(Rx,Ry,clrnum);             % Plot borderline
                    text(Rx(1),Ry(1),num2str(Vn));                          % Write bordline number
                    text(Rx(fix(length(Rx)/2)),Ry(fix(length(Rx)/2)),num2str(Vn));    % Write bordline number
                    text(Rx(length(Rx)),Ry(length(Rx)),num2str(Vn));        % Repeat 
                end
                
                Vn=Vn+1;
                clear  Lx Ly Rx Ry Ep NL pr
     end % while aw