%徑向基網絡（先聚類，后訓練）
%先以K－均值聚類算法聚類，確定聚類中心，后訓練網絡，高斯函數的寬度參數取定值sigma=max([sum((z1-z2).^2),sum((z1-z3).^2),sum((z3-z2).^2)])/sqrt(2*3)
% 運行時間較長，查全率、查準率也比較差，原因可能是某些不同樣本之間特征值差別不大，導致聚類不準確，聚類中心與實際中心不一致，（結果與聚類訓練同時進行
% 的差別不大）
clear Re Pr
c=load('ccc.txt');
q1=3 ;%q1為常向量，表示神經元的個數
Ir=0.0121; max_epoch=1500;
EE=zeros(length(q1),max_epoch);
res=zeros(size(q1));%res將要存儲誤差向量，這里先置零 
 aa=1;
clear oi ok wwji wwkj b1 b2 deltak deltao yy result cout1 count2 E
tic  
ceshiyangben=c;err_goal=0.05;
  for i=1:59
  target(i,:)=[0 0 1];
  end
  for i=60:131
  target(i,:)=[0 1 0];
  end
  for i=132:178
  target(i,:)=[1 0 0];
  end
  E=[];T=target';%mubiaoxiangliang
  q=q1(aa);
  [m,n]=size(ceshiyangben);[s2,s1]=size(T); 
  for i=2:n;
   for j=1:m
   ceshiyangben(j,i)=ceshiyangben(j,i)/max(ceshiyangben(:,i));
   end
  end %歸一化
 w=ceshiyangben';
  z1=w(:,1);z2=w(:,2);z3=w(:,3); I=1;   %1
for j=1:15
w1=[];w2=[];w3=[];
for i=1:178
if (sqrt(sum(w(:,i)-z1(:,j))^2)<sqrt(sum(w(:,i)-z2(:,j))^2))&(sqrt(sum(w(:,i)-z1(:,j))^2)<sqrt(sum(w(:,i)-z3(:,j))^2))
    w1=[w1,w(:,i)];[h,n1]=size(w1);
     else  if (sqrt(sum(w(:,i)-z2(:,j))^2)<sqrt(sum(w(:,i)-z1(:,j))^2))&(sqrt(sum(w(:,i)-z2(:,j))^2)<sqrt(sum(w(:,i)-z3(:,j))^2))
              w2=[w2,w(:,i)];[h,n2]=size(w2);
          else w3=[w3,w(:,i)];[h,n3]=size(w3);
          end
end
end
z1=[z1,mean(w1')'];z2=[z2,mean(w2')'];z3=[z3,mean(w3')'];
w11=w1;w22=w2;w33=w3;
[a,i]=size(z1); [a,k]=size(z2); [a,m]=size(z3);
if (z1(:,i)==z1(:,i-1))&(z2(:,k)==z2(:,k-1))&(z3(:,m)==z3(:,m-1))
    break
end
end
   Wki=0.2*rands(3,q); %
   b2=0.1*rands(3,1);%—隨機給定隱含層、輸出層偏值 
   wwkj=[b2 Wki];
   P=ceshiyangben';
   count1=zeros(1,3);
   count2=zeros(1,3);
   z1=z1(:,j);z2=z2(:,j);z3=z3(:,j);   %聚合中心
   r1=sqrt(sum((P-z1*ones(1,178)).^2));r2=sqrt(sum((P-z2*ones(1,178)).^2));r3=sqrt(sum((P-z3*ones(1,178)).^2));  % ||X-W||
   for kk=1:max_epoch
     sigma=max([sum((z1-z2).^2),sum((z1-z3).^2),sum((z3-z2).^2)])/sqrt(2*3);
     ojj1=exp(-r1.^2/(2*sigma^2));               %高斯函數
     ojj2=exp(-r2.^2/(2*sigma^2));
     ojj3=exp(-r3.^2/(2*sigma^2));
     ojj=[ojj1;ojj2;ojj3];
     oj1=[-ones(178,1) ojj'];
     netk=wwkj*oj1';
     oo=logsig(netk);%3-178
     errr=T-oo;   
     sse=mse(errr);
     deltao=errr.*oo.*(ones(3,178)-oo);%3-178
     if(sse<err_goal)
     break;
     end
     wwkj=wwkj+Ir*(deltao*oj1);
     E=[E,sse];
 end
  res(aa)=norm(errr);
  EE(aa,:)=E;
  toc
   %yangbenceshi 
    r1=sqrt(sum((P-z1*ones(1,178)).^2));r2=sqrt(sum((P-z2*ones(1,178)).^2));r3=sqrt(sum((P-z3*ones(1,178)).^2)); 
     ojjj1=exp(-r1.^2/(2*sigma^2));  %11 178
     ojjj2=exp(-r2.^2/(2*sigma^2));
     ojjj3=exp(-r3.^2/(2*sigma^2));
     ojjj=[ojj1;ojj2;ojj3];
     ojj1=[-ones(178,1) ojj'];
     netkk=wwkj*oj1';
     ooo=logsig(netk);%3-178
     for j=1:s2
        for i=1:178
       yy(j,i)=round(oo(j,i));
        end
     end
  for h=1:178
   if yy(:,h)==T(:,h)
      if h<=59
      result(1,h)=1;count1(1,1)=count1(1,1)+1;
      elseif h<=131&h>59
      result(1,h)=2;count1(1,2)=count1(1,2)+1;
      elseif h<=178&h>131
      result(1,h)=3;count1(1,3)=count1(1,3)+1;
      end
    else
    result(1,h)=0;   
   end
  end
    Rresult(aa,:)=result;

  for hh=1:178
    if yy(:,hh)==[0 0 1]'
        count2(1,1)=count2(1,1)+1;
    elseif  yy(:,hh)==[0 1 0]' 
        count2(1,2)=count2(1,2)+1;
    elseif  yy(:,hh)==[1 0 0]' 
       count2(1,3)=count2(1,3)+1;  
    end
  end
  Re(aa,:)=count1./[59 72 47];
  Pr(aa,:)=count1./count2;
  for dd=1:length(q1)
    plot(EE(dd,:))
    hold on
  end
 number=find(res==min(res));
Re3=Re
Pr3=Pr
figure
plot(Rresult(1,:)')