%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %          程序運行過程中調用:ddd函數,這個函數計算在本聚類模式下的兩個模糊
                              %          數的距離          
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FKMCMD(k,mm,x,e)      %其中，e表示預設的閾值；k類，mm表示模糊關系數，x為待聚類的模糊數據矩陣，
                               %每一行表示一個模糊數.若每一個模糊數包含m個不同的數據單元,則要求x為一個n*(4m)的矩陣.其中n表示
                               %模糊數的個數,每一行中前m個表示該模糊數據的m個中心位置,隨后的m個數據表示內徑，在后
                               %面的m個元素表示左區間長，最后面的m個元素表右區間長。
n=size(x);
m=n(2)/4;                      %求出模糊數據向量的規模，即其包含幾個部分
n=n(1);                        %求出x的模糊數的個數n

                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %       判斷類別數是不是大于數據個數      %
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k>n
    fprintf('ERROR! The number of clusters is greater than the number of the data, namely %d>%d. Please revise your cluster number!\n',k,n)
else
    
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %                             聲明變量空間                                                 
                              %    u：為初始關系矩陣
                              %    uu：為經優化后的關系矩陣
                              %    c：為數據的中心矩陣（這是一個矩陣，因為x為一個向量模糊數）                
                              %    g：為數據的內徑矩陣
                              %    l：為數據的左區間長矩陣
                              %    r：為數據的右區間長矩陣
                              %    C：為k類初始中心矩陣
                              %    G：為k類初始內徑矩陣
                              %    L：為k類初始左區間長矩陣
                              %    R：為k類初始右區間長矩陣
                              %    CC：為k類經優化后的中心矩陣
                              %    GG：為k類經優化后的內徑矩陣
                              %    LL：k類經優化后的左區間長矩陣
                              %    RR：k類經優化后的右區間長矩陣
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=ones(k,n)/n;
uu=ones(k,n);
c=ones(n,m);
g=ones(n,m);
l=ones(n,m);
r=ones(n,m);
C=ones(k,m);
G=ones(n,m);
L=ones(n,m);
R=ones(n,m);
CC=ones(k,m);
GG=ones(n,m);
LL=ones(n,m);
RR=ones(n,m);

                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %                         初始賦值
                              %   第一部分，將x分為中心，內徑，左區間長和右區間長，分別賦值給c，g，l和r
                              %   第二部分，選取k個聚類中心，內徑，左區間長和右區間長。不妨取為前k個元素，其中為了計算需要對中心做一個平移變換
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              
for i=1:n                                %第一部分                         
    c(i,:)=x(i,1:m);
    g(i,:)=x(i,(m+1):2*m);
    l(i,:)=x(i,(2*m+1):3*m);
    r(i,:)=x(i,(3*m+1):4*m);
end
for j=1:k                                %第二部分  
    C(j,:)=c(j,:)-0.1;
    G(j,:)=g(j,:);
    L(j,:)=l(j,:);
    R(j,:)=r(j,:);
end


                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %                  計算部分
                              %     第三部分，計算優化后的關系矩陣。其中用到的算法是根據論文中提供的方法。
                              %     第四部分，第五部分，第六部，和第七部分分分別為計算優化后的聚類中心，內徑，
                              %                    左區間和右區間長度矩陣。其中用到的算法是根據論文中提供的方法。
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d=ones(1,k);
for i=1:n                                %第三部分                  
    for j=1:k
        d(j)=1/(ddd(x(i,:),[C(j,:), G(j,:),L(j,:),R(j,:)])^(2/(mm-1)));
    end
    for j=1:k
        uu(j,i)=1/(ddd(x(i,:),[C(j,:), G(j,:),L(j,:),R(j,:)])^(2/(mm-1))*sum(d));
    end
    d=ones(1,k);
end

for i=1:k                               %第四部分
    CC(i,:)=0;
    for j=1:n
        CC(i,:)=uu(i,j).^mm*(8*c(j,:)-l(j,:)+r(j,:)+L(i,:)-R(i,:))+CC(i,:);
    end
    CC(i,:)=CC(i,:)/(8*sum(uu(i,:).^mm));
end
for i=1:k                               %第五部分
    GG(i,:)=0;
    for j=1:n
        GG(i,:)=uu(i,j).^mm*(4*g(j,:)+l(j,:)+r(j,:)-L(i,:)-R(i,:))+GG(i,:);
    end
    GG(i,:)=GG(i,:)/(8*sum(uu(i,:).^mm));
end
for i=1:k                               %第六部分
    LL(i,:)=0;
    for j=1:n
        LL(i,:)=uu(i,j).^mm*(-2*c(j,:)+g(j,:)+l(j,:)+2*C(i,:)-G(i,:))+LL(i,:);
    end
    LL(i,:)=LL(i,:)/(sum(uu(i,:).^mm));
end
for i=1:k                               %第七部分
    RR(i,:)=0;
    for j=1:n
        RR(i,:)=uu(i,j).^mm*(2*c(j,:)+g(j,:)+r(j,:)-2*C(i,:)-G(i,:))+RR(i,:);
    end
    RR(i,:)=RR(i,:)/(sum(uu(i,:).^mm));
end

                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %             進入循環項：若u和uu相應元素的差的最大值大于預設閾值e，則進入循環
                              %  將新算出的結果賦給就變量，重新計算相應的關系矩陣uu和聚類中心，內徑，左區間和右區間矩陣，直到循環條件不被滿足，跳出循環。
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while max(max(abs(u-uu)))>e 
    L=LL;
    R=RR;
    G=GG;
    C=CC;
    u=uu;
    d=ones(1,k);
for i=1:n                                           
    for j=1:k
        d(j)=1/(ddd(x(i,:),[C(j,:), G(j,:),L(j,:),R(j,:)])^(2/(mm-1)));
    end
    for j=1:k
        uu(j,i)=1/(ddd(x(i,:),[C(j,:), G(j,:),L(j,:),R(j,:)])^(2/(mm-1))*sum(d));
    end
    d=ones(1,k);
end

for i=1:k                               
    CC(i,:)=0;
    for j=1:n
        CC(i,:)=uu(i,j).^mm*(8*c(j,:)-l(j,:)+r(j,:)+L(i,:)-R(i,:))+CC(i,:);
    end
    CC(i,:)=CC(i,:)/(8*sum(uu(i,:).^mm));
end
for i=1:k                              
    GG(i,:)=0;
    for j=1:n
        GG(i,:)=uu(i,j).^mm*(4*g(j,:)+l(j,:)+r(j,:)-L(i,:)-R(i,:))+GG(i,:);
    end
    GG(i,:)=GG(i,:)/(8*sum(uu(i,:).^mm));
end
for i=1:k                            
    LL(i,:)=0;
    for j=1:n
        LL(i,:)=uu(i,j).^mm*(-2*c(j,:)+g(j,:)+l(j,:)+2*C(i,:)-G(i,:))+LL(i,:);
    end
    LL(i,:)=LL(i,:)/(sum(uu(i,:).^mm));
end
for i=1:k                               
    RR(i,:)=0;
    for j=1:n
        RR(i,:)=uu(i,j).^mm*(2*c(j,:)+g(j,:)+r(j,:)-2*C(i,:)-G(i,:))+RR(i,:);
    end
    RR(i,:)=RR(i,:)/(sum(uu(i,:).^mm));
end
end

                              %%%%%%%%%%%%%%%%%%%%%%%%%%
                              %    輸出計算的結果，分別為聚類中心和關系矩陣
                              %%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('The centers of the %d clustering is:\n',k)
CC
fprintf('The membership of the elements relating to the %d clusters is:\n',k)
uu=uu'
fprintf('The graph has shown the result of the clustering, in which the red point is the center of the clusters.\n')
                              %%%%%%%%%%%%%%%%%%%%%%%%%
                              %     畫圖
                              %%%%%%%%%%%%%%%%%%%%%%%%%
hold on
if m==2
for i=1:n
    plot([c(i,1)],[c(i,2)],'*')
    plot([c(i,1)-l(i,1),c(i,1)+l(i,1),c(i,1)+l(i,1),c(i,1)-l(i,1),c(i,1)-l(i,1)],[c(i,2)-l(i,2),c(i,2)-l(i,2),c(i,2)+l(i,2),c(i,2)+l(i,2),c(i,2)-l(i,2)])
end
else 
    for i=1:n
    plot([c(i,1)],[c(i,2)],'*')

    end
end
for i=1:k
plot([CC(i,1)],[CC(i,2)],'Marker','o','LineStyle','none','Color',[1 0 0])
end
hold off
end