% Nois Resistant Fuzzy Inference System without Membership Function (NRFIS-NOMF)

clear all
clc
load TAIEX1999
r=length(XI(:,1));
N=length(XI(1,:));
e=0.00001;
C=5;
m=2;
Ti=1000;

% FCM algorithm

% Calculating covariance norm matrix, A.

sumXI=0;
for j=1:N
    sumXI=sumXI+XI(:,j);
end
XIave=sumXI/N;
Cov=0;
for j=1:N
    Cov=Cov+(XI(:,j)-XIave)*(XI(:,j)-XIave)';
end
Cov=Cov/N;
A=inv(Cov);
% A=eye(r,r);

% Starting FCM.

U=rand(C,N);
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)*XI(:,j);
            sumU=sumU+((U(i,j))^m);
        end
        V0=sumV/sumU;
        for j0=1:r
            V(j0,i)=V0(j0);
        end
    end
    
%     Calculating partition matrix.
    
    for i=1:C
        for j=1:N
            dji=(XI(:,j)-V(:,i))'*A*(XI(:,j)-V(:,i));
            for k=1:C
                djk=(XI(:,j)-V(:,k))'*A*(XI(:,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

% Calculating partition matrix,U.

clear U
for q=1:r
    for i=1:C
        for j=1:N
            dji=(XI(q,j)-V(q,i))'*(XI(q,j)-V(q,i));
            for k=1:C
                djk(k)=(XI(q,j)-V(q,k))'*(XI(q,j)-V(q,k));
                mf(k)=(dji/djk(k))^(1/(m-1));
            end
            sum=0;
            for k=1:C
                sum=sum+mf(k);
            end
            U(q,i,j)=1/sum;
        end
    end
end
    
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************
% ***********************************************************************

% Tunning THEN part of the system.

% Generating XIm matrix.

for j=1:N
    XIm(1,j)=1;
end
for q=2:r+1
    for j=1:N
        XIm(q,j)=XI(q-1,j);
    end
end

% Calculating X matrix.

for i=1:C
    for j=1:N
        tuU(i,j)=1;
        for q=1:r
            tuU(i,j)=tuU(i,j)*U(q,i,j);
        end
    end
end

for j=1:N
    sumtuU(j)=0;
    for i=1:C
        sumtuU(j)=sumtuU(j)+tuU(i,j);
    end
end

for j=1:N
    for i=1:C
        tuUw(i,j)=tuU(i,j)/sumtuU(j);
    end
end

for i=1:C
    for j=1:N
        tuw(i,j)=tuUw(i,j);
    end
end

for j=1:N
    cc=0;
    for i=1:C
        for q=1:r+1
            cc=cc+1;
            X(j,cc)=tuw(i,j)*XIm(q,j);
        end
    end
end

P=pinv(X'*X)*X'*yd';

% Extracting THEN part parameters.

cc=0;
for i=1:C
    for j=1:r+1
        cc=cc+1;
        PC(i,j)=P(cc);
    end
end

% Calculating output of the system for a given input.

 for j=1:N
    Xj=XIm(:,j);
    for q=1:r
        for i=1:C
            sumU0=0;
            for k=1:C
                dji=(XI(q,j)-V(q,i))'*(XI(q,j)-V(q,i));
                djk=(XI(q,j)-V(q,k))'*(XI(q,j)-V(q,k));
                sumU0=sumU0+((dji/djk)^(1/(m-1)));
            end
            U0(q,i)=1/sumU0;
        end
    end
    for i=1:C
        tuU0(i)=1;
        for q=1:r
            tuU0(i)=tuU0(i)*U0(q,i);
        end
    end
    sumtuU0=0;
    for i=1:C
        sumtuU0=sumtuU0+tuU0(i);
    end
    for i=1:C
        tuU0(i)=tuU0(i)/sumtuU0;
    end
    for i=1:C
        tu0(i)=tuU0(i);
    end
    for i=1:C
        fo(i)=PC(i,:)*Xj;
    end
    ym(j)=0;
    for i=1:C
        ym(j)=ym(j)+tu0(i)*fo(i);
    end
end

% RMSEtrain

t0=1:N;
plot(t0,yd,'b',t0,ym,'r')

RMSEtrain=(yd-ym)*(yd-ym)';
RMSEtrain=sqrt(RMSEtrain/N);
MSEtrain=RMSEtrain^2;

% Testing data evaluation.

% Generating XIm matrix.

for j=1:length(XIorig)
    XImorig(1,j)=1;
end
for q=2:r+1
    for j=1:length(XIorig)
        XImorig(q,j)=XIorig(q-1,j);
    end
end

for j=1:length(XIorig)
    Xj=XImorig(:,j);
    for q=1:r
        for i=1:C
            sumU0=0;
            for k=1:C
                dji=(XIorig(q,j)-V(q,i))'*(XIorig(q,j)-V(q,i));
                djk=(XIorig(q,j)-V(q,k))'*(XIorig(q,j)-V(q,k));
                sumU0=sumU0+((dji/djk)^(1/(m-1)));
            end
            U0(q,i)=1/sumU0;
        end
    end
    for i=1:C
        tuU0(i)=1;
        for q=1:r
            tuU0(i)=tuU0(i)*U0(q,i);
        end
    end
    sumtuU0=0;
    for i=1:C
        sumtuU0=sumtuU0+tuU0(i);
    end
    for i=1:C
        tuU0(i)=tuU0(i)/sumtuU0;
    end
    for i=1:C
        tu0(i)=tuU0(i);
    end    
    for i=1:C
        fo(i)=PC(i,:)*Xj;
    end
    yo(j)=0;
    for i=1:C
        yo(j)=yo(j)+tu0(i)*fo(i);
    end
end

% Test evaluataion.

for j=N+1:length(XIorig)
    yts(j-N)=yo(j);
    ydts(j-N)=ydorig(j);
end
RMSEtest=(ydts-yts)*(ydts-yts)';
RMSEtest=sqrt(RMSEtest/length(yts));
MSEtest=RMSEtest^2;
RMSEtrain
RMSEtest

t1=1:N;
t2=N+1:length(ydorig);
L=min(ydorig):0.001*max(ydorig):max(ydorig);
for i=1:length(L)
    t3(i)=N;
end
plot(t1,yd,'b',t1,ym,'r',t2,ydts,'b',t2,yts,'r',t3,L,'s')

