%% nsfls1.m

%% Compute the output(s) 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.

%% 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;
%% R1, R2, R are the left, right, and average output

function [R1,R2,R]=nsfls_type2_2(X,D,M1,M2,sigma,c1,c2,sn1,sn2);

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


R1=[];
R2=[];
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

h=(UU+LL)/2;
delta=UU-h;

[l_out,r_out] = interval_wtdavg(c0',s',h,delta);

R1=[R1,l_out];
R2=[R2,r_out];
end

R3=(R1+R2)/2;
