%BP算法作業
%題目一：訓練函數y=sinx，x∈（0，2π）
%%%%%%%%%輸入輸出變換采用sigmoid函數，采用9個訓練樣本，361個測試樣本
%%%%%%%%%性能指標精度ε取0.01
%%%%%%%%%BP采用單隱層
function MLP_BP1
global N Nt IN HN ON miu eita ep;
[w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold]=ini_BP;
[P,T]=get_Trainsample();
judge=0;%判斷標志
epoch=1;%訓練周期數
iN=1;%標志訓練樣本
delta_O=zeros(ON,N);%訓練誤差
delta_H=zeros(HN,N);%訓練誤差
while judge==0  
    for iN=1:N
        [input_H,output_H,input_O,output_O]=get_Trainoutput(w_H,P,theta_H,w_O,theta_O,iN);
        [delta_O,delta_H]=train_Error(T,output_O,w_O,output_H,iN);
        [w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold]=modify_Para(delta_H,P,delta_O,output_H,w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold,iN);
    end   %完成一次訓練（一個訓練周期）
    [judge,error]=E_Judge(T,output_O);
    error_Store(epoch)=error;
    epoch=epoch+1;
end
test_Output;


%初始化BP網絡
function [w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold]=ini_BP
N=9;%訓練樣本個數
Nt=361;%測試樣本個數

IN=1;%輸入層神經元數目
HN=4;%隱層神經元個數
ON=1;%輸出層神經元個數

%alpha_O=0.01;%輸出層至隱層的學習效率
%alpha_H=0.01;%隱層至輸入層學習效率

ep=0.01;%設定精度

miu=0.1;%μ
eita=0.05;%η

%以下初始化權值和閾值（Nguyen和Widrow的權值初始化算法）
lamda=0.7^IN*sqrt(HN);
w_H=rand(HN,IN)-0.5*ones(HN,IN);%在-0.5至0.5范圍內
w_Hold=zeros(HN,IN);
w_O=rand(ON,HN)-0.5*ones(HN,IN);%在-0.5至0.5范圍內
w_Oold=zeros(ON,HN);
%輸入層到隱層的權值初始化
for i=1:HN    
    for j=1:IN
        sum=0;
        for k=1:HN
            sum=w_H(k,j)^2+sum;
        end
        w_H(i,j)=lamda.*w_H(i,j)./sqrt(sum);
    end
end
%隱層到輸出層的權值初始化
for i=1:ON    
    for j=1:HN
        sum=0;
        for k=1:ON
            sum=w_O(k,j)^2+sum;
        end
        w_O(i,j)=lamda.*w_O(i,j)./sqrt(sum);
    end
end
%隱層閾值初始化
theta_H=zeros(HN,1);
theta_Hold=theta_H;
for i=1:HN
    j=ceil(IN*rand);
    theta_H(i,1)=w_H(i,j)*(2*rand-1);
end
%輸出層閾值初始化
theta_O=zeros(ON,1);
theta_Oold=theta_O;
for i=1:ON
    j=ceil(HN*rand);
    theta_O(i,1)=w_O(i,j)*(2*rand-1);
end
%訓練樣本初始化
function [P,T]=get_Trainsample()
global IN ON N;
P=zeros(IN,N);%單個樣本輸入數據
T=zeros(ON,N);%單個樣本教師數據,期望輸出
for i=1:IN    
    P(i,:)=0:(2*pi/(N-1)):(2*pi);  
end
for i=1:ON
    T(i,:)=sin(P(i,:));
end
plor(P,T,'+');
%計算各層輸出
function  [input_H,output_H,input_O,output_O]=get_Trainoutput(w_H,P,theta_H,w_O,theta_O,iN)
global IN HN ON N;
input_H=zeros(HN,N);%隱層的輸入
output_H=zeros(HN,N);%隱層的輸出
input_O=zeros(ON,N);%輸出層的輸入
output_O=zeros(ON,N);%輸出層的輸出
%計算隱層輸入輸出
input_H(1)=0;
for i=1:HN
    for j=1:IN
        input_H(i,iN)=w_H(i,j)*(j,iN)+input_H(i,iN);
    end
    output_H(i,iN)=input_H(i,iN)-theta_H(i);
    output_H(i,iN)=1/(1+exp(-output_H(i,iN)));
end
%計算輸出層輸入輸出
input_O(1)=0;
for i=1:ON
    for j=1:HN
        input_O(i,iN)=w_O(i,j)*input_H(j,iN)+input_O(i,iN);
    end
    output_O(i,iN)=input_O(i,iN)-theta_O(i);
    output_O(i,iN)=1/(1+exp(-output_O(i,iN)));
end
return;
%計算各層訓練誤差
function [delta_O,delta_H]=train_Error(T,output_O,w_O,output_H,iN)
global ON HN;
%計算輸出層訓練誤差
for i=1:ON    
    delta_O(i,iN)=(T(i,iN)-output_O(i,iN))*output_O(i,iN)*(1-output_O(i,iN));
end
%計算隱層訓練誤差
for i=1:HN
    for j=1:ON
    delta_H(i,iN)=delta_H(i,iN)+delta_O(j,iN)*w_O(j,i);
    end
    delta_H(i,iN)=output_H(i,iN)*(1-output_H(i,iN))*delta_H(i,iN);
end

%修正權值和閾值
function [w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold]=modify_Para(delta_H,P,delta_O,output_H,w_H,w_Hold,w_O,w_Oold,theta_H,theta_Hold,theta_O,theta_Oold,iN)
global IN HN ON miu eita;
%修正權值
%輸入層到隱層
w_Hold=w_H;
for i=1:HN
    for j=1:IN
        w_H(i,j)=w_H(i,j)+miu*delta_H(i,iN)*P(j,iN)+eita*(w_H(i,j)-w_Hold(i,j));
    end
end
%隱層到輸出層
w_Oold=w_O;
for i=1:ON
    for j=1:HN
        w_O(i,j)=w_O(i,j)+miu*delta_O(i,iN)*output_H(j,iN)+eita*(w_O(i,j)-w_Oold(i,j));
    end
end
%修正閾值
%隱層
theta_Hold=theta_H;
for i=1:HN
    theta_H(i)=theta_H(i)+miu*delta_H(i,iN)+eita*(theta_H(i)-theta_Hold(i));
end
%輸出層
theta_Oold=theta_O;
for i=1:ON
    theta_O(i)=theta_O(i)+miu*delta_O(i,iN)+eita*(theta_O(i)-theta_Oold(i));
end
%判斷精度有無滿足要求
function [judge,error]=E_Judge(T,output_O)
global N ON ep;
error=0;%樣本誤差(代價函數的值)
for j=1:N
    for i=1:ON
        error=0.5*(T(i,j)-output_O(i,j))^2+error;
    end
end
if error <= ep
    judge=1;
else
    judge=0;
end
return;

function test_Output(epoch,error_Store,)
global IN ON Nt;
Pt=zeros(IN,Nt);%測試單個樣本輸入數據
Tt=zeros(ON,Nt);
Tt=zeros(ON,Nt);%測試單個樣本輸出數據
for i=1:IN    
    P(i,:)=0:(2*pi/(Nt-1)):(2*pi);  
end
for i=1:ON
    T(i,:)=sin(P(i,:));
end
plor(P,T,'+');