%%%%%%%%%%%%%%%%%%%%%%%%%
                              %          程序運(yùn)行過程中調(diào)用:dd函數(shù),這個(gè)函數(shù)計(jì)算在本聚類模式下的兩個(gè)模糊數(shù)的距離       
                              %%%%%%%%%%%%%%%%%%%%%%%%%

function FKMCCFV(k,mm,x,e)            %其中，e表示預(yù)設(shè)的閾值；k類，mm表示模糊關(guān)系數(shù)，x為待聚類的模糊數(shù)據(jù)矩陣，
                                      %每一行表示一個(gè)模糊數(shù).若每一個(gè)模糊數(shù)包含m個(gè)不同的數(shù)據(jù)單元,則要求x為一個(gè)n*(2m)的矩陣.其中n表示
                                      %模糊數(shù)的個(gè)數(shù),每一行中前m個(gè)表示該模糊數(shù)據(jù)的m個(gè)中心位置,后m個(gè)元素分別
                                      %對(duì)應(yīng)每一個(gè)模糊數(shù)據(jù)的左右區(qū)間長(zhǎng),所以這種方法只能考慮左右對(duì)稱的模糊數(shù)據(jù)向量類型.
n=size(x);
m=n(2)/2;
n=n(1);  %求出x的模糊數(shù)的個(gè)數(shù)n
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %       判斷類別數(shù)是不是大于數(shù)據(jù)個(gè)數(shù)      %
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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：為初始關(guān)系矩陣
                              %    c：為數(shù)據(jù)的中心矩陣（這是一個(gè)矩陣，因?yàn)閤為一個(gè)向量模糊數(shù)）                
                              %    s：為數(shù)據(jù)的區(qū)間長(zhǎng)矩陣
                              %    C：為k類初始中心矩陣
                              %    S：為k類初始區(qū)間長(zhǎng)矩陣
                              %    CC：為k類經(jīng)優(yōu)化后的中心矩陣
                              %    uu：為經(jīng)優(yōu)化后的關(guān)系矩陣
                              %    a：為三角關(guān)系運(yùn)算所對(duì)應(yīng)的對(duì)角的矩陣，只是在matlab中無法聲明a(i)，所以用一個(gè)三維數(shù)組表示
                              %    AA：為經(jīng)優(yōu)化后的三角關(guān)系運(yùn)算所對(duì)應(yīng)的對(duì)角的矩陣。
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%

u=ones(k,n)/n;
c=ones(n,m);
s=ones(n,m);
C=ones(k,m);
S=ones(k,m);
CC=ones(k,m);
uu=ones(k,n);
a=zeros(m,m,n);
AA=zeros(m,m,k);



                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %                         初始賦值
                              %   第一部分，將x分為中心和區(qū)間長(zhǎng)，分別賦值給c和u
                              %   第二部分，選取k個(gè)聚類中心，和區(qū)間長(zhǎng)。不妨取為前k個(gè)元素，其中為了計(jì)算需要對(duì)中心做一個(gè)平移變換
                              %   第三部分，根據(jù)定義，賦值三角關(guān)系運(yùn)算矩陣a
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n                               %第一部分
    c(i,1:m)=x(i,1:m);
    s(i,1:m)=x(i,m+1:2*m);
end
for i=1:k                               %第二部分
    C(i,:)=c(i,:)-0.1;
    S(i,:)=0.1;
end
for i=1:n                               %第三部分
    for j=1:m
        a(j,j,i)=1/s(i,j)^2;
    end
end

                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %                  計(jì)算部分
                              %     第四部分，計(jì)算優(yōu)化后的關(guān)系矩陣。其中用到的算法是根據(jù)論文中提供的方法。
                              %     第五部分和第六部分分別為計(jì)算優(yōu)化后的聚類中心和區(qū)間長(zhǎng)度矩陣。其中用到的算法是根據(jù)論文中提供的方法。
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%
d=ones(1,k);
for i=1:n                                %第四部分                  
    for j=1:k
        d(j)=1/(dd(x(i,:),[C(j,:),S(j,:)])^(2/(mm-1)));
    end
    for j=1:k
        uu(j,i)=1/(dd(x(i,:),[C(j,:),S(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*c(j,:)+CC(i,:);
    end
    CC(i,:)=CC(i,:)/(sum(uu(i,:).^mm));
end
for i=1:k                               %第六部分
    AA(:,:,i)=0;
    for j=1:n
        AA(:,:,i)=uu(i,j).^mm*a(:,:,j)+AA(:,:,i);
    end
    AA(:,:,i)=AA(:,:,i)/(sum(uu(i,:).^mm));
end



                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                              %             進(jìn)入循環(huán)項(xiàng)：若u和uu相應(yīng)元素的差的最大值大于預(yù)設(shè)閾值e，則進(jìn)入循環(huán)
                              %  將新算出的結(jié)果賦給就變量，重新計(jì)算相應(yīng)的關(guān)系矩陣uu和聚類中心CC及三角關(guān)系矩陣AA直到循環(huán)條件不被滿足，跳出循環(huán)。
                              %%%%%%%%%%%%%%%%%%%%%%%%%%%%
while max(max(abs(u-uu)))>e 
    C=CC;
    u=uu;
    d=ones(1,k);
    for i=1:n                   
    for j=1:k
        d(j)=1/(dd(x(i,:),[C(j,:),S(j,:)])^(2/(mm-1)));
    end
    for j=1:k
        uu(j,i)=1/(dd(x(i,:),[C(j,:),S(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*c(j,:)+CC(i,:);
    end
    CC(i,:)=CC(i,:)/(sum(uu(i,:).^mm));
    end
for i=1:k
    AA(:,:,i)=0;
    for j=1:n
        AA(:,:,i)=uu(i,j).^mm*a(:,:,j)+AA(:,:,i);
    end
    AA(:,:,i)=AA(:,:,i)/(sum(uu(i,:).^mm));
end
end    


                              %%%%%%%%%%%%%%%%%%%%%%%%%%
                              %    輸出計(jì)算的結(jié)果，分別為聚類中心和關(guān)系矩陣
                              %%%%%%%%%%%%%%%%%%%%%%%%%%
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')
fprintf('Attention: Because of the axis in the two graphs is not constrained!.\n')

                              %%%%%%%%%%%%%%%%%%%%%%%%%
                              %     畫圖
                              %%%%%%%%%%%%%%%%%%%%%%%%%
hold on
if m==2
for i=1:n
    plot([c(i,1)],[c(i,2)],'*')
    plot([c(i,1)-s(i,1),c(i,1)+s(i,1),c(i,1)+s(i,1),c(i,1)-s(i,1),c(i,1)-s(i,1)],[c(i,2)-s(i,2),c(i,2)-s(i,2),c(i,2)+s(i,2),c(i,2)+s(i,2),c(i,2)-s(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