% General Entropy based Possibilistic Fuzzy C-Means (GEPFCM)-Relaxed.
 
% Note: GEPFCM is a generalized entropy based form of PFCM.
% In GEPFCM, the distance ||Xj-Vi||A is replaced with arbitrary functions
% fFCM,f1PCM f1E, f2E and f2PCM of ||Xj-Vi||A.
% The function fFCM can be chosen such that it be nonzero for all x.
% This choice will remove the problem of FCM and PFCM when ||Xj-Vi||A=0 
% and membership grades uij become undetermined. Note that by changing
% values of cFCM, cPCM, c1E and c2E the GEPFCM can be changed into
% FCM, PCM and PFCM. More over, by changing this coefficients,
% importance of memberships, typicalities and entropy can be changed.

clear all
clc
load Data1

N=length(yd(1,:));
D=length(yd(:,1));
e=0.00001;
cFCM=1;
cPCM=1;
c1E=1;
c2E=1;
eta=2;
m=2;
% C=2;
KF=1;
% ruopt=10;
Ti=1000;

% FCM algorithm
% Using results of the FCM, initial cluster centers (prototypes),
% intial partition matrix and values of gama for PFCM are obtained.
% Applying PFCM, new U,V and gama matrices are obtained by which
% GPFCM algorithm is initialized.

% Calculating covariance norm matrix, A.

sumYD=0;
for j=1:N
    sumYD=sumYD+yd(:,j);
end
YDave=sumYD/N;
Cov=0;
for j=1:N
    Cov=Cov+(yd(:,j)-YDave)*(yd(:,j)-YDave)';
end
Cov=Cov/N;
A=inv(Cov);
% A=eye(D,D);

% FCM.

U=rand(C,N);
% load U
Norm(1)=0;
for j=1:N
    for i=1:C
        Norm(1)=Norm(1)+U(i,j)^2;
    end
end
Norm(1)=sqrt(Norm(1));

for t=2:Ti
    
%     Cluster centers

    for i=1:C
        sumV=0;
        sumU=0;
        for j=1:N
            sumV=sumV+((U(i,j))^m)*yd(:,j);
            sumU=sumU+((U(i,j))^m);
        end
        V0=sumV/sumU;
        for j0=1:D
            V(j0,i)=V0(j0);
        end
    end
    
%     Calculating partition matrix.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            for k=1:C
                djk=(yd(:,j)-V(:,k))'*A*(yd(:,j)-V(:,k));
                mf(k)=(dji/djk)^(1/(m-1));
            end
            sumU=0;
            for k=1:C
                sumU=sumU+mf(k);
            end
            U(i,j)=1/sumU;
        end
    end

%     Checking convergence
    
    Norm(t)=0;
    for j=1:N
        for i=1:C
            Norm(t)=Norm(t)+U(i,j)^2;
        end
    end
    Norm(t)=sqrt(Norm(t));
    E0(t-1)=abs(Norm(t)-Norm(t-1));
    Eo=E0(t-1);
    t
    Eo
    if E0(t-1)<e
        break
    else
    end
end

% PFCM algorithm

% Calculating gama

K=1;
for i=1:C
    sumU=0;
    for j=1:N
        sumU=sumU+((U(i,j))^m);
    end
    sumG=0;
    for j=1:N
        sumG=sumG+(U(i,j)^m)*((yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i)));
    end
    gama(i)=K*sumG/sumU;
end

Norm(1)=0;
for j=1:N
    for i=1:C
        Norm(1)=Norm(1)+U(i,j)^2;
    end
end
Norm(1)=sqrt(Norm(1));

for t=2:Ti
    
    %     Calculating typicality matrix,T.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            T(i,j)=(1+(((cPCM/gama(i))*dji)^(1/(eta-1))))^(-1);
        end
    end
    
%     Cluster centers

    for i=1:C
        sumV=0;
        sumUT=0;
        for j=1:N
            sumV=sumV+((cFCM*(U(i,j))^m)+(cPCM*(T(i,j))^eta))*yd(:,j);
            sumUT=sumUT+((cFCM*(U(i,j))^m)+(cPCM*(T(i,j))^eta));
        end
        V0=sumV/sumUT;
        for j0=1:D
            V(j0,i)=V0(j0);
        end
    end

%     Calculating partition matrix,U.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            for k=1:C
                djk=(yd(:,j)-V(:,k))'*A*(yd(:,j)-V(:,k));
                mf(k)=(dji/djk)^(1/(m-1));
            end
            sumU=0;
            for k=1:C
                sumU=sumU+mf(k);
            end
            U(i,j)=1/sumU;
        end
    end

%     Checking convergence
    
    Norm(t)=0;
    for j=1:N
        for i=1:C
            Norm(t)=Norm(t)+U(i,j)^2;
        end
    end
    Norm(t)=sqrt(Norm(t));
    E0(t-1)=abs(Norm(t)-Norm(t-1));
    Eo=E0(t-1);
    t
    Eo
    if E0(t-1)<e
        break
    else
    end
end

K=1;
for i=1:C
    sumU=0;
    for j=1:N
        sumU=sumU+((U(i,j))^m);
    end
    sumG=0;
    for j=1:N
        sumG=sumG+(U(i,j)^m)*((yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i)));
    end
    gama(i)=K*sumG/sumU;
end

% GEPFCM algorithm

for iR=1:KF
    
%     Determination of radii of the clusters based on PCNM.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            Tf(i,j)=(1+(((cPCM/gama(i))*dji)^(1/(eta-1))))^(-1);
        end
    end
    
    for i=1:C
        sum1=0;
        sum2=0;
        for j=1:N
            sum1=sum1+(Tf(i,j)^eta)*(yd(:,j)-V(:,i))*(yd(:,j)-V(:,i))';
            sum2=sum2+(Tf(i,j)^eta);
        end
        PCNM=inv(sum1/sum2);
        sum1=0;
        sum2=0;
        for j=1:N
            sum1=sum1+(Tf(i,j)^eta)*(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            sum2=sum2+(Tf(i,j)^eta);
        end
        R(i)=sqrt(sum1/sum2);
        omega(i)=ruopt/(R(i)^2);
    end
    
    Norm(1)=0;
    for j=1:N
        for i=1:C
            Norm(1)=Norm(1)+U(i,j)^2;
        end
    end
    Norm(1)=sqrt(Norm(1));
    
%     Calculating gama
    
    for j=1:length(yd(1,:))
        Xj=yd(:,j);
        for i=1:C
            dji=(Xj-V(:,i))'*A*(Xj-V(:,i));
            for k=1:C
                djk=(Xj-V(:,k))'*A*(Xj-V(:,k));
                mf(k)=(dji/djk)^(1/(m-1));
            end
            sumU=0;
            for k=1:C
                sumU=sumU+mf(k);
            end
            Uf(i,j)=1/sumU;
        end
    end
    
    for i=1:C
        sumU=0;
        for j=1:N
            sumU=sumU+((Uf(i,j))^m);
        end
        sumG=0;
        for j=1:N
            sumG=sumG+(Uf(i,j)^m)*((yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i)));
        end
        gama(i)=K*sumG/sumU;
    end
    
%     Controlling entropy term for switching between algorithms.
    
    if c1E~=0
        for i=1:C
            c1Ei(i)=omega(i);
        end
    else
        for i=1:C
            c1Ei(i)=0;
        end
    end
    
    for t=2:Ti

    %     Calculating S matrix.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            S(i,j)=exp(-1-(c1Ei(i)/c2E)*f1E(dji));
        end
    end
    
    %     Cluster centers.
    
    for i=1:C
        sumV=0;
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            sumV=sumV+((cFCM*(U(i,j))^m)*fFCMp(omega(i),dji)+(cPCM*(T(i,j))^eta)*f1PCMp(omega(i),dji)+gama(i)*((1-T(i,j))^eta)*f2PCMp(omega(i),dji)+c1Ei(i)*S(i,j)*f1Ep(dji))*yd(:,j);
        end
        sumUTS=0;
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            sumUTS=sumUTS+(cFCM*(U(i,j))^m)*fFCMp(omega(i),dji)+(cPCM*(T(i,j))^eta)*f1PCMp(omega(i),dji)+gama(i)*((1-T(i,j))^eta)*f2PCMp(omega(i),dji)+c1Ei(i)*S(i,j)*f1Ep(dji);
        end
        V0=sumV/sumUTS;
        for j0=1:D
            V(j0,i)=V0(j0);
        end
    end
    
    %     Calculating U matrix.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            for k=1:C
                djk=(yd(:,j)-V(:,k))'*A*(yd(:,j)-V(:,k));
                mf(k)=(fFCM(omega(i),dji)/fFCM(omega(i),djk))^(1/(m-1));
            end
            sumU=0;
            for k=1:C
                sumU=sumU+mf(k);
            end
            U(i,j)=1/sumU;
        end
    end
    
%     Calculating T matrix.
    
    for i=1:C
        for j=1:N
            dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
            T(i,j)=(1+(((cPCM/gama(i))*(f1PCM(omega(i),dji)/f2PCM(omega(i),dji)))^(1/(eta-1))))^(-1);
        end
    end

%     Checking convergence
    
    Norm(t)=0;
    for j=1:N
        for i=1:C
            Norm(t)=Norm(t)+U(i,j)^2;
        end
    end
    Norm(t)=sqrt(Norm(t));
    E0(t-1)=abs(Norm(t)-Norm(t-1));
    Eo=E0(t-1);
    t
    Eo
    if E0(t-1)<e
        break
    else
    end
    end
end

% Elimination of outliers: For this purpose, we use above method of
% clustering with appropriate functions fFCM, f1PCM and f2PCM to find
% accurate cluster centers and then using GPFCM with
% fFCM=f1PCM=x and f2PCM=1, we find T and U with Euclidean norm.
% Because covariance norm considers effects of outliers it should not
% be used. Note that if the data is noisy, in calculation of
% final U and T matrices, using covariance norm matrix is not recommended
% because it brings effects of the noises into these matrices.

for j=1:length(yd(1,:))
    Xj=yd(:,j);
    for i=1:C
        dji=(Xj-V(:,i))'*A*(Xj-V(:,i));
        for k=1:C
            djk=(Xj-V(:,k))'*A*(Xj-V(:,k));
            mf(k)=(dji/djk)^(1/(m-1));
        end
        sumU=0;
        for k=1:C
            sumU=sumU+mf(k);
        end
        Uf(i,j)=1/sumU;
    end
end
K=1;
for i=1:C
    sumU=0;
    for j=1:N
        sumU=sumU+((Uf(i,j))^m);
    end
    sumG=0;
    for j=1:N
        sumG=sumG+(Uf(i,j)^m)*((yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i)));
    end
    gama(i)=K*sumG/sumU;
end
for i=1:C
    for j=1:N
        dji=(yd(:,j)-V(:,i))'*A*(yd(:,j)-V(:,i));
        Tf(i,j)=(1+(((cPCM/gama(i))*dji)^(1/(eta-1))))^(-1);
    end
end
U=Uf;
T=Tf;

plot(yd(1,:),yd(2,:),'bo')
hold on
plot(V(1,:),V(2,:),'r*')

% plot3(yd(1,:),yd(2,:),yd(3,:),'o')
% hold on
% plot3(V(1,:),V(2,:),V(3,:),'r*')
