function [yUMF,uUMF,yLMF,uLMF]=lwaPR(X, W, nAlpha)

%
% [yUMF,uUMF,yLMF,uLMF]=lwaPR(X, W, nAlpha)
%
% LWA computation for PR
%
% Dongrui Wu (dongruiw@usc.edu), 4/25/2009
%
% X: M*9 matrix storing the FOUs of the rule consequents
% W: Firing levels of the rules computed by the Jaccard similarity measure
% nAlpha: numer of alpha-cuts. nAlpha>1. Default 2.
%
% yUMF: y-coordinates of the UMF of Y_PR
% uUMF: u-coordinates of the UMF of Y_PR
% yLMF: y-coordinates of the LMF of Y_PR
% uLMF: u-coordinates of the LMF of Y_PR

%% select only the fired rules
index=find(W(:,2)>0);
X=X(index,:);
W=W(index,:);

M=size(X,1);
UMFXx=X(:,1:4);
LMFXx=X(:,5:8);
UMFXmu=zeros(M,4);
LMFXmu=zeros(M,4);

for i=1:M
    UMFXmu(i,:) =[0 1 1 0];
    LMFXmu(i,:) =[0 X(i,9) X(i,9) 0];
end

a=zeros(M,1);
b=zeros(M,1);
c=W(:,1);
d=W(:,2);

mu=[0:1/(nAlpha-1):1 1:-1/(nAlpha-1):0]; % Descretize Y domain
uUMF=mu;  %% u-coordinates of the UMF of the LWA

l1=zeros(1,nAlpha);
l2=zeros(1,nAlpha);
r1=zeros(1,nAlpha);
r2=zeros(1,nAlpha);

%% %%%%%%%%%%%%%%%%%%%% UMF   %%%%%%%%%%%%%%%%%%%%
%%% alpha==0
for j=1:M  %% for each input, compute the alpha-cut
    %% a,b: alpha-cut on X
    temp=find(UMFXmu(j,:)>0);
    indexa=min(temp)-1;  %% locate the interval where mu(i) lies in
    a(j)=UMFXx(j,indexa);
    indexb=max(temp)+1;
    b(j)=UMFXx(j,indexb);
end
l1(1)=EKM(a,c,d,-1); %% EKM algorithm to compute the alpha-cut on the UMF of the LWA
r2(1)=EKM(b,c,d,1);

%% alpha>0
for i=2:nAlpha %% for each alpha-cut
    for j=1:M  %% for each input, compute the alpha-cut
        %% a,b: alpha-cut on X
        temp=find(UMFXmu(j,:)>=mu(i));
        indexa=min(temp);  %% locate the interval where mu(i) lies in
        a(j)=UMFXx(j,indexa-1)+(UMFXx(j,indexa)-UMFXx(j,indexa-1))...
            *(mu(i)-UMFXmu(j,indexa-1))/(UMFXmu(j,indexa)-UMFXmu(j,indexa-1));
        indexb=max(temp);
        b(j)=UMFXx(j,indexb)+(UMFXx(j,indexb+1)-UMFXx(j,indexb))...
            *(mu(i)-UMFXmu(j,indexb))/(UMFXmu(j,indexb+1)-UMFXmu(j,indexb));
    end
    l1(i)=EKM(a,c,d,-1); %% EKM algorithm to compute the alpha-cut on the UMF of the LWA
    r2(i)=EKM(b,c,d,1);
end
yUMF=[l1 r2(nAlpha:-1:1)];

%% %%%%%%%%%%%%%%%%%%%%% LMF %%%%%%%%%%%%%%%%%%%%%%%
mu=min(X(:,9))*[0:1/(nAlpha-1):1 1:-1/(nAlpha-1):0]; % Descretize Y domain
uLMF=mu;  %% u-coordinates of the UMF of the LWA
%% alpha==0
for j=1:M  %% for each input, compute the alpha-cut
    %% a,b: alpha-cut on X
    temp=find(LMFXmu(j,:)>0);
    indexa=min(temp)-1;  %% locate the interval where mu(i) lies in
    a(j)=LMFXx(j,indexa);
    indexb=max(temp)+1;
    b(j)=LMFXx(j,indexb);
end
l2(1)=EKM(a,c,d,-1);
r1(1)=EKM(b,c,d,1);

%% alpha>0
for i=2:nAlpha
    %% a,b: alpha-cut on X
    for j=1:M  %% for each input, compute the alpha-cut
        %% a,b: alpha-cut on X
        temp=find(LMFXmu(j,:)>=mu(i));
        indexa=min(temp);  %% locate the interval where mu(i) lies in
        a(j)=LMFXx(j,indexa-1)+(LMFXx(j,indexa)-LMFXx(j,indexa-1))...
            *(mu(i)-LMFXmu(j,indexa-1))/(LMFXmu(j,indexa)-LMFXmu(j,indexa-1));
        indexb=max(temp);
        b(j)=LMFXx(j,indexb)+(LMFXx(j,indexb+1)-LMFXx(j,indexb))...
            *(mu(i)-LMFXmu(j,indexb))/(LMFXmu(j,indexb+1)-LMFXmu(j,indexb));
    end
    l2(i)=EKM(a,c,d,-1);
    r1(i)=EKM(b,c,d,1);
end
yLMF=[l2 r1(nAlpha:-1:1)];
