clc,clear all;close all;
mu1 = [0 0]';
sigma1 = 1;
mu2 = [2 0]';
sigma2 = 2;
Num_Input=500;
Ro = 2;          %Row
Cl = Num_Input;  %Coulmn
Num_Epoch =100;  %Num of Epochs
mu1 = repmat(mu1,[1 Cl]);
mu2 = repmat(mu2,[1 Cl]);
class1 = normrnd(mu1,sigma1,Ro,Cl);
class2 = normrnd(mu2,sigma2,Ro,Cl);


% fclose(fid1);fclose(fid2);
if(1)
    figure;
    plot(class1(1,:),class1(2,:),'o');
    figure;
    plot(class2(1,:),class2(2,:),'r*');
    figure;
    plot(class1(1,:),class1(2,:),'o');hold on;
    plot(class2(1,:),class2(2,:),'r*');
    axis([-4 10 -6 6]);
end
%% ==============   Design   ================
Input1 = [ones(1,floor(Num_Input/2))*-1;class1(1:floor(Num_Input/2))];
Num_Neuron = 10;
WL1 = 0.1 * rand(3,Num_Neuron)-.05;
WL2 = 0.1 * rand(Num_Neuron,1)-.05;
D=[0,1];      % for Phi=1/1+exp(-V)
Eta = 0.6;    % Rate of Learning
Alpha=0.2;    % Momentom
Del_wl2=0.1 * rand(Num_Neuron,1)-0.05;
Del_wl1=0.1 * rand(Num_Neuron,3)-0.05;
d1=zeros(1,floor(Num_Input/2));d2=zeros(1,floor(Num_Input/2))+1;
data = [ones(1,Num_Input)*-1;class1(:,1:floor(Num_Input/2)),class2(:,1:floor(Num_Input/2));d1,d2];
Learning_Curve = zeros(Num_Epoch,1);
h = waitbar(0,'Please wait...');
for Epoch=1:Num_Epoch
    NN = round(Epoch/Num_Epoch*100);
    waitbar(NN/100);
    Index = randperm(Num_Input);
    Ave_Err=0;
    for i=1:Num_Input
        Idx = Index(i);
        % ======  Forward  =========
        V1 =  WL1' * data(1:3,Idx);
        Y1=Phi(V1);
        V2 =  WL2' * Y1;
        Y2=Phi(V2);
        e=data(4,Idx)-Y2;
        Ave_Err = Ave_Err + e^2;

        % ======= Back ward ===========
        Grad1 = e*Y2.*(1-Y2);
        Del_wl2=Alpha*Del_wl2+Eta*Grad1*Y1;
        WL2=WL2+Del_wl2;

        Grad2 = Grad1*WL2.*Y1.*(1-Y1);
        Del_wl1=Alpha*Del_wl1+Eta*Grad2*data(1:3,Idx)';
        WL1=WL1+Del_wl1';
    end
    Learning_Curve(Epoch) = Ave_Err/Num_Input;
end
close(h) 
MSE = sum(Learning_Curve)/Num_Epoch
figure(1);
plot(1:Num_Epoch,Learning_Curve)
hold on; grid on;
xlabel('Number of Epochs')
ylabel('MSE over Epoch')

%%  =========================  Test  ============================
Test_Input = [ones(1,Num_Input)*-1;class1(:,floor(Num_Input/2)+1:Num_Input),class2(:,floor(Num_Input/2)+1:Num_Input);d1,d2];

Error = 0;
for i=1:Num_Input
    % ======  Forward  =========
    V1 =  WL1' * Test_Input(1:3,i);
    Y1=Phi(V1);
    V2 =  WL2' * Y1;
    Y2=Phi(V2);
    e=abs(Test_Input(4,i)-Y2);
    if(e>0.5)
        Error = Error+1;
    end

end
Ec = (Num_Input-Error)/Num_Input*100

%% = ================  Ml Algorithm  =================
mu1 = [0 0]';
sigma1 = 1;
mu2 = [2 0]';
sigma2 = 2;
r = sqrt(sigma1^2*sigma2^2*(norm(mu1+mu2)^2/(sigma2^2-sigma1^2) + 4*log(sigma2/sigma1))/(sigma2^2-sigma1^2) );
Xc = (sigma2^2*mu1 - sigma1^2*mu2)/(sigma2-sigma1);
Err_ML=0;
for i=1:Num_Input
    if( norm(Test_Input([2 3],i) - Xc )>r ) % Select Class 2
        if(Test_Input(4,i) ==  0)
            Err_ML=Err_ML+1;
        end
    else  % Select Class 1
        if(Test_Input(4,i) ==  1)
            Err_ML=Err_ML+1;
        end
    end
end
Ec_ML = (Num_Input-Err_ML) /Num_Input*100

function [Y]= Phi(V)
Y = 1 ./ ( 1 + exp(-V) );
% a = 1.7;
% b = 2/3;
% Y = a * tanh( b * V);