%'因子數n1,預報對象數 n2,試報年數n,學習樣本數 m,，mm為總樣本數，學習因子a1，動量因子 a2,總體誤差ee1
%個體誤差ee2,迭代次數 num,隱層節點數 yj，
close all
clear  all
clc
[file,path]=uigetfile('*.*','打開數據文件');
prompt={'包括Y的列數','樣本總個數','建模樣本數','預報因子數','預報對象數','學習因子','動量因子','訓練次數','總體誤差','隱節點數'};
def={'4','45','40','3','1','0.9','0.7','5','0.001','3'};
dlgTitle='參數輸入';
lineNo=1;
answer=inputdlg(prompt,dlgTitle,lineNo,def);
answer=char(answer);
h=str2num(answer(1,:));
mm=str2num(answer(2,:));
m=str2num(answer(3,:));
n1=str2num(answer(4,:));
n2=str2num(answer(5,:));
a1=str2num(answer(6,:));
a2=str2num(answer(7,:));
num=str2num(answer(8,:));
ee1=str2num(answer(9,:));
yj=str2num(answer(10,:));
filename1=strcat(path,file);
fid1=fopen(filename1,'r');
a=fscanf(fid1,'%f',h*mm);
fclose(fid1);
a=reshape(a,h,mm);
a=a';
p=a;
[file,path]=uiputfile('*.*','保存');
filename2=strcat(path,file);
fid2=fopen(filename2,'w');
n=mm-m;
hdata=p;
hda(3,(h))=0;
fdata(mm,(h))=0;                                                                                                                                                                                                                                                                                                                                                                           
s1=min(hdata(1:m,:));
s2=max(hdata(1:m,:));
hda(1,1:h)=mean(hdata(1:m,1:(h)));
mean(hdata(1:m,1:(h)));
hda(2,1:h)=(9*s2-s1)./8;
hda(3,1:h)=(9*s1-s2)./8;
for i=1:mm
    for j=1:(n1+n2)
        if hdata(i,j)~=-999
        fdata(i,j)=(hdata(i,j)-hda(3,j))/(hda(2,j)-hda(3,j));
        else
            fdata(i,j)=(hda(1,j)-hda(3,j))/(hda(2,j)-hda(3,j));
        end
    end
end
%-----------------------------------------fdata 中放的是標準化的資料，gdata 中放的是預報量資料
 gdata=fdata(1:mm,(n1+1):h)  
 fdata=fdata(:,1:n1) 
 ggdata=p(1:mm,(n1+1):h)
 %----------------------------------------------------------  
v=randn(n1,yj);
vc=randn(yj);
w=randn(yj,n2);
rc=randn(n2);
f=fdata;
g=gdata;
m_data(m+n,n2)=0;
n_dd=0;
ee=0.1;
eee=0
%計算
n_dd=1
bar=waitbar(0,'please wait....');
for n_dd=1:num
 waitbar(n_dd/num,bar);
        %給連接權付初值0
        v1=zeros(n1,yj);
        vc1=zeros(yj);
        w1=zeros(yj,n2);
        rc1=zeros(n2);
                                              
 for i=1:m          
    % '計算隱含層的激活值  
    for j=1:yj                             
       xx=0 ;                               
        for k=1:n1
          xx=xx+f(i,k)*v(k,j);  
        end
        xx=xx+vc(j);
       b(j)=1/(1+exp(-xx));
    end                                      
    % 計算輸出層單元的激活值                                       
       for j=1:n2
       xx=0;
       for k=1:yj
       xx=xx+b(k)*w(k,j);
       end
       xx=xx+rc(j);
       c(j)=1/(1+exp(-xx));
       mdata(i,j)=c(j);
       end
    % 計算輸出層單元的一般化誤差
    for j=1:n2
        d(j)=c(j)*(1-c(j))*(g(i,j)-c(j));
    end
     % 計算隱含層對于每個dj誤差
     for j=1:yj
         xx=0;
         for k=1:n2
             xx=xx+w(j,k)*d(k);
             
         end
         e(j)=b(j)*(1-b(j))*xx;
     end
      %調整隱含層單元 到輸出層單元的連接權
      for j=1:yj
          for k=1:n2
              w1(j,k)=w1(j,k)+a1*b(j)*d(k);
          end
      end
       %調整 輸出層單元的閾值
       for j=1:n2
         rc1(j)=rc1(j)+a1*d(j);
       end
      % 調整輸入層單元 到隱含層單元的連接權
      for j=1:n1
          for k=1:yj
              v1(j,k)=v1(j,k)+a2*f(i,j)*e(k);
          end
      end
       % 計算隱含層單元的閾值
       for j=1:yj
           vc1(j)=vc1(j)+a2*e(j);
       end
 end
 %-------------------------------------  
   labuda=0.15;
   %調整連接權
   for j=1:n1
       for k=1:yj
           v(j,k)=v(j,k)+labuda*v1(j,k);
       end
   end
   for j=1:yj
       vc(j)=vc(j)+labuda*vc1(j);
   end
   for j=1:yj
       for k=1:n2
           w(j,k)=w(j,k)+labuda*w1(j,k);
       end
   end
   for j=1:n2
       rc(j)=rc(j)+labuda*rc1(j);
   end
   %計算總體誤差
   
   ee=0;
   for i=1:m
       for j=1:n2
           ee=ee+(mdata(i,j)-g(i,j))^2/2;
       end
   end
   ee=ee/m;
   eee(n_dd)=ee;
   
if ee<ee1
   break
end
end
close(bar);
ee;
plot(eee);
fprintf(fid2,'包括預報量的列數：%1d,樣本總數：%2d,建模數：%1d,預報因子數：%1d,預報量數：%1d,隱節點數：%1d\n',h,mm,m,n1,n2,yj);
fprintf(fid2,'訓練次數：%d:,動量因子：%4.1f,學習因子:%4.1f,總體誤差：%f\n',num,a1,a2,ee1);
fprintf(fid2,'收斂誤差為：%8.6f\n\n',ee) ;
sprintf('收斂誤差為：%12.6f',ee) 
n_dd;
%預報
for i=1:mm
    %計算隱含層新的激活值
   for j=1:yj
       xx=0;
       for k=1:n1
           xx=xx+f(i,k)*v(k,j);
       end
       xx=xx+vc(j);
       b(j)=1/(1+exp(-xx));
   end
   %
     for j=1:n2
       xx=0;
       for k=1:yj
           xx=xx+b(k)*w(k,j);
       end
       xx=xx+rc(j);
       c(j)=1/(1+exp(-xx));
      mdata(i,j)=c(j);
     end
 end

mdata;
  % transfering the original value
for i=1:mm
    for j=1:n2
        if mdata(i,j)~=-999
          mdata(i,j)=mdata(i,j)*(hda(2,h)-hda(3,h))+hda(3,h)  ;
    
      elseif mdata(i,j)==-999
            
            mdata(i,j)=(hda(1,j)-hda(3,h))/(hda(2,h)-hda(3,h));
        end
    end
end
mdata
fprintf(fid2,'predictand  fitting\n\n');
  fprintf(fid2,'              實況        預測       誤差       相對誤差\n');
for i=1:m
    for j=1:n2
        xd(i,j)=abs(ggdata(i,j)-mdata(i,j))/ggdata(i,j)*100;
        fprintf(fid2,'        %12.4f%12.4f%12.4f%12.4f%%',ggdata(i,j),mdata(i,j),(ggdata(i,j)-mdata(i,j)),abs(ggdata(i,j)-mdata(i,j))/ggdata(i,j)*100);
    end
    fprintf(fid2,'\n');
end
jd1=mean(abs(ggdata(1:m,:)-mdata(1:m,:)));
fprintf(fid2,'average %12.4f%12.4f%12.4f%12.4f%%\n\n',mean(ggdata(1:m,:)),mean(mdata(1:m,:)),jd1,mean(xd(1:m,:)));
    fprintf(fid2,'\npredictand  prediction\n\n');
for i=m+1:mm
    for j=1:n2
        xd(i,j)=abs(ggdata(i,j)-mdata(i,j))/ggdata(i,j)*100;
        fprintf(fid2,'        %12.4f%12.4f%12.4f%12.4f%%',ggdata(i,j),mdata(i,j),(ggdata(i,j)-mdata(i,j)),abs(ggdata(i,j)-mdata(i,j))/ggdata(i,j)*100);
    end
    fprintf(fid2,'\n');
end
    fprintf(fid2,'\n\n');    

jd2=mean(abs(ggdata((m+1):mm,:)-mdata((m+1):mm,:)));
       
       fprintf(fid2,'\n\n  average%12.4f%12.4f%12.4f%12.4f%%\n\n',mean(ggdata(1+m:mm,:)),mean(mdata(1+m:mm,:)),jd2,mean(xd(1+m:mm,1)));
       fprintf(fid2,'\n');
    
 fprintf(fid2,'\n\n');
 fprintf(fid2,'誤差變化：\n');
 ne=length(eee)
for i=1:ne
fprintf(fid2,'%4d %12.6f\n',i,eee(i));
end
  fclose(fid2);