%% train_nsfls.m

%% Tune the parameters 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 inputÐoutput training data.

%% M1,M2, sigma are mxn matrix denotes the mean and std of
%% antecedent Gaussian MFs (N rules, with n antecedent in each rule);
%% c1,c2 are Nx1 vectors, which denote the COS of consequents;
%% sn is the input std;
%% X is input matrix, Lxant matrix, each row is onw input.
%% D is Lx1 vector which denotes the desired output
%% alpha is the step size for tuning M, sigma, c1,c2;
%% alpha4  is the step size for tuning sn;


function [M1,M2,c1,c2,sigma,sn]=train_nsfls_type2_1(X,D,M1,M2,sigma,c1,c2,sn,alpha,alpha4);

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

%%%%%% Following are the training for type-2 FLS %%%%%%%

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)^2*(M2(j,m)-M1(j,m))/(2*sigma(j,m)^2) &...
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=(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

%% compute the output of FLS

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

f=(l_out+r_out)/2;
e=D(i)-f;


fa1=wl/sum(wl);
fa2=wr/sum(wr);

%%% Select MFs contributed to the right-most 

ME10=M1;
ME20=M2;
sn0=sn;

LE=length(I2u);

for t=1:LE
for k=1:ant 
if X(i,k)<M1(I2u(t),k)
l=I2u(t);
M1(l,k)=M1(l,k)+alpha*e*0.5*((X(i,k)-M1(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c2(l)-r_out)*wr(l)/sum(wr);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));
elseif  X(i,k)>M2(I2u(t),k)
l=I2u(t);
M2(l,k)=M2(l,k)+alpha*e*0.5*((X(i,k)-M2(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c2(l)-r_out)*wr(l)/sum(wr);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));
end
end
end


LE=length(I2l);

for t=1:LE
for k=1:ant
l=I2l(t);
if X(i,k)<=(M1(l,k)+M2(l,k))/2-sn(k)^2*(M2(l,k)-M1(l,k))/(2*sigma(l,k)^2)
M2(l,k)=M2(l,k)+alpha*e*0.5*((X(i,k)-M2(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c2(l)-r_out)*wr(l)/sum(wr);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));
elseif X(i,k)>=(M1(l,k)+M2(l,k))/2+sn(k)^2*(M2(l,k)-M1(l,k))/(2*sigma(l,k)^2)
M1(l,k)=M1(l,k)+alpha*e*0.5*((X(i,k)-M1(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c2(l)-r_out)*wr(l)/sum(wr);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c2(l)-r_out)*wr(l)/sum(wr));
end
end
end   


%%% Select MFs contributed to the left-most 

LE=length(I1l);

for t=1:LE
for k=1:ant 
l=I1l(t);
if X(i,k)<=(M1(l,k)+M2(l,k))/2-sn(k)^2*(M2(l,k)-M1(l,k))/(2*sigma(l,k)^2)
M2(l,k)=M2(l,k)+alpha*e*0.5*((X(i,k)-M2(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c1(l)-l_out)*wl(l)/sum(wl);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));
elseif X(i,k)>=(M1(l,k)+M2(l,k))/2+sn(k)^2*(M2(l,k)-M1(l,k))/(2*sigma(l,k)^2)
M1(l,k)=M1(l,k)+alpha*e*0.5*((X(i,k)-M1(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c1(l)-l_out)*wl(l)/sum(wl);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));
end
end
end



LE=length(I1u);

for t=1:LE
for k=1:ant
if  X(i,k)<M1(I1u(t),k)
l=I1u(t);
M1(l,k)=M1(l,k)+alpha*e*0.5*((X(i,k)-M1(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c1(l)-l_out)*wl(l)/sum(wl);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME10(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));
elseif X(i,k)>M2(I1u(t),k) 
l=I1u(t);
M2(l,k)=M2(l,k)+alpha*e*0.5*((X(i,k)-M2(l,k))/(sigma(l,k)^2+sn(k)^2))...
*(c1(l)-l_out)*wl(l)/sum(wl);
sigma(l,k)=sigma(l,k)+alpha*e*0.5*(((sigma(l,k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));

sn0(k)=sn0(k)+alpha4*e*0.5*(((sn(k)*(X(i,k)-ME20(l,k))^2)/(sigma(l,k)^2+sn(k)^2)^2)...
*(c1(l)-l_out)*wl(l)/sum(wl));
end
end
end   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sn=sn0;

fa1=wr'/sum(wr);
fa2=wl'/sum(wl);
c1=c1+alpha*e*fa1/2;
c2=c2+alpha*e*fa2/2;

%%%% After Training, The type-2 MFs should be reconstructed


for l=1:N
for k=1:ant

if sigma(l,k) < 0
sigma(l,k)=abs(sigma(l,k));
end

if sn(k) < 0
sn(k)=abs(sn(k));
end

if sn(k) > sigma(l,k)
sn(k)=sigma(l,k);
end

end
end

mc2=max([c2,c1]');
mc1=min([c2,c1]');
c2=mc2';
c1=mc1';

ME1=[];
ME2=[];

for t=1:N
P=[M2(t,:)',M1(t,:)']';
m2=max(P);
m1=min(P);
ME1=[ME1,m1'];
ME2=[ME2,m2'];
end

M1=ME1';
M2=ME2';

end
