[M1,M2,sigma,c1,c2]=svd_qr_nsfls_type2_1(X,M1,M2,sigma,c1,c2,sn);

%% svd_qr_nsfls1.m

%% Rule-reduction of an interval type-1 non-singleton type-2 FLS 
%% when the antecedent membership functions are Gaussian primary 
%% membership functions with uncertain means and the input sets are type-1 
%% Gaussian, using some inputoutput training data.

%% M1,M2,sigma are mxn matrix denotes the mean and std of
%% antecedent Gaussian MFs (m rules, with n antecedent in each rule)
%% c1,c2 are mx1 vectors, which denotes the COS of consequents
%% X is input matrix, Lxn matrix, each row is an input.
%% sn is the input std
%% call the function gausstype2.m
%% th is a threshold for choosing the principal components

function [M1,M2,sigma,c1,c2]=svd_qr_nsfls_type2_1(X,M1,M2,sigma,c1,c2,sn);

FBFu=[];
FBFl=[];

[L,ant]=size(X);
[N,ant]=size(M);
rule=N;

c0=(c1+c2)/2;
s=c2-c0;


for i=1:L
U=[];
MU1=[];

UU=[];
LL=[];
for j=1:N
Uu=1;
Ll=1;
for m=1:ant

if X(i,m) <= M1(j,m)
xupp=((sn(m)^2)*M1(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);
xlow=((sn(m)^2)*M2(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);

elseif X(i,m)>M1(j,m) & X(i,m)<= (M1(j,m)+M2(j,m))/2-sn(m)^2*(M2(j,m)-M1(j,m))/(2*sigma(j,m)^2)
xupp=X(i,m);
xlow=((sn(m)^2)*M2(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);

elseif X(i,m)>(M1(j,m)+M2(j,m))/2-sn(m)*(M2(j,m)-M1(j,m))/(2*sigma(j,m)) &...
X(i,m)<(M1(j,m)+M2(j,m))/2+sn(m)*(M2(j,m)-M1(j,m))/(2*sigma(j,m))
xupp=X(i,m);
xlow=(M1(j,m)+M2(j,m))/2;

elseif X(i,m)>=(M1(j,m)+M2(j,m))/2+sn(m)^2*(M2(j,m)-M1(j,m))/(2*sigma(j,m)^2) & X(i,m) < M2(j,m)
xupp=X(i,m);
xlow=((sn(m)^2)*M1(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);

elseif X(i,m) >= M2(j,m)
xupp=((sn(m)^2)*M2(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);
xlow=((sn(m)^2)*M1(j,m)+(sigma(j,m)^2)*X(i,m))/(sn(m)^2+sigma(j,m)^2);
end

P=[sigma(j,m),M1(j,m),M2(j,m)];
[uu,ll]=gausstype2(xupp,P);

uu=uu*gaussmf(xupp,[sn(m), X(i,m)]);

[tmp,ll]=gausstype2(xlow,P);
ll=ll*gaussmf(xlow,[sn(m), X(i,m)]);


Uu=Uu*uu;
Ll=Ll*ll;
end
UU=[UU,Uu];
LL=[LL,Ll];
end

[r_out,I2l,I2u,wr]= rightpoint(c2',LL,UU);
[l_out,I1u,I1l,wl]= leftpoint(c1',LL,UU);


FBFu=[FBFu,wr'/sum(wr)];
FBFl=[FBFl,wl'/sum(wl)];
end

FBFu=FBFu';
FBFl=FBFl';

%%%%%%%%%%
[U,SS,V]=svd(FBFu);

a=diag(SS);
for i=1:rule
if a(i)>th
r=i;
end
end

A=[V(1:r,1:r)', V((r+1):rule,1:r)'];
[Q,R,E]=qr(A);
[e,I]=sort(E);
I1=I(rule,1:r);
%%%%%%%%%%%%%
[U,SS,V]=svd(FBFl);

a=diag(SS);
for i=1:N
if a(i)>th
r=i;
end
end

A=[V(1:r,1:r)', V((r+1):rule,1:r)'];
[Q,R,E]=qr(A);
[e,I]=sort(E);
I2=I(rule,1:r);
%%%%%%%%%%%%
I=union(I1,I2);
r=length(I);

%%%%%%% Keep the most important parameters

M1=M1(I,:);
M2=M2(I,:);
c1=c1(I);
c2=c2(I);
sigma=sigma(I,:);