%% with one extra data pre-processing step.

function [MF,n, Stat, NNum, Shift] = Fuzzistics_FL(L,R)

%% This program is to find the FOU based on the collected data from JPL.
%% The objective is to test our approach on finding the FOU.
%% MF contains the parameters of the FOU
%% n is the number of usable data after pre-processing


%% Step 1: Data Processing (Outlier processing)
%%
L1 = L; %% Left data after preprocessing stage 1
R1 = R; %% Stage 1

%%% delete these interval data with Lower<0 or Upper >10 and Upper<Lower.
n = size(L,1);
t =0;
for i=1:n
    if L(i)>=0 & L(i)<=10 & R(i)>=0 & R(i)<=10 & R(i)>=L(i)
        t=t+1;
        newL(t) = L(i);
        newR(t) = R(i);
        tempNum1(t)=i;
    else
        L1(i) = -1;
        R1(i) = -1;
    end
end
clear L R
L = newL(1:t);
R = newR(1:t);
n = t;
Num1 = t;
clear newL newR


%% --------- Outlier Processing ----------------------------------------
%% Outlier processing by Quantiles IQR= Q(0.75)-Q(0.25)  [Q(0.25)-1.5IQR, Q(0.75)+1.5IQR]
L2 = L1; %% Left data after preprocessing stage 2
R2 = R1; %% Stage 2
InLen = R-L;

Left = sort(L);
Right = sort(R);
Length = sort(InLen);

%%-- Compute Q(0.25) and Q(0.75)
NN1 = floor(n * 0.25 + 0.5);
NN2 = floor(n * 0.75 + 0.5);
QL25 = (0.5 - n * 0.25 + NN1) * Left(NN1) + (n * 0.25 + 0.5 - NN1) * Left(NN1+1);
QL75 = (0.5 - n * 0.75 + NN2) * Left(NN2) + (n * 0.75 + 0.5 - NN2) * Left(NN2+1);
LIQR = QL75 - QL25;

QU25 = (0.5 - n * 0.25 + NN1) * Right(NN1) + (n * 0.25 + 0.5 - NN1) * Right(NN1+1);
QU75 = (0.5 - n * 0.75 + NN2) * Right(NN2) + (n * 0.75 + 0.5 - NN2) * Right(NN2+1);
UIQR = QU75 - QU25;

QLen25 = (0.5 - n * 0.25 + NN1) * Length(NN1) + (n * 0.25 + 0.5 - NN1) * Length(NN1+1);
QLen75 = (0.5 - n * 0.75 + NN2) * Length(NN2) + (n * 0.75 + 0.5 - NN2) * Length(NN2+1);
LenIQR = QLen75 - QLen25;

% disp('QL25, QL75,LIQR; QU25, QU75, UIQR; QLen25, QLen75, LenIQR')
% [QL25, QL75,LIQR; QU25, QU75, UIQR; QLen25, QLen75, LenIQR]
% disp('QL25-1.5*LIQR, QL75+1.5*LIQR; QU25-1.5*UIQR, QU75+1.5*UIQR; QLen25 - 1.5*LenIQR, QLen75 + 1.5*LenIQR')
% [QL25-1.5*LIQR, QL75+1.5*LIQR; QU25-1.5*UIQR, QU75+1.5*UIQR; QLen25 - 1.5*LenIQR, QLen75 + 1.5*LenIQR]
%InLen

tempNum = Num1;
Bound = 0.25;
%%---- For left-end point outlier processing -------------
t =0;
if (LIQR > Bound)
    for i=1:n
        if L(i)>=QL25-1.5*LIQR & L(i)<=QL75+1.5*LIQR
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            newInLen(t) = InLen(i);
            tempNum2(t)= tempNum1(i);
        else
            L2(tempNum1(i)) = -1;
            R2(tempNum1(i)) = -1;
        end
    end
    clear L R InLen;
    L = newL;
    R = newR;
    InLen = newInLen;
    n = t;
    tempNum =t;
    clear newL newR newInLen;
else
    tempNum2 = tempNum1;
end

t=0;
if (UIQR > Bound)
    for i=1:n
        if R(i)>=QU25-1.5*UIQR & R(i)<=QU75+1.5*UIQR
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            newInLen(t) = InLen(i);
            tempNum2(t)= tempNum2(i);
        else
            L2(tempNum2(i)) = -1;
            R2(tempNum2(i)) = -1;
        end
    end
    clear L R InLen;
    L = newL;
    R = newR;
    InLen = newInLen;
    n = t;
    tempNum =t;
    clear newL newR
end

t=0;
if (LenIQR > Bound)
    for i=1:n
        if  (InLen(i)>= QLen25 - 1.5*LenIQR) & (InLen(i)<= QLen75 + 1.5*LenIQR)
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            tempNum2(t)= tempNum2(i);
        else
            L2(tempNum2(i)) = -1;
            R2(tempNum2(i)) = -1;
        end
    end
    clear L R;
    L = newL;
    R = newR;
    n = t;
    tempNum =t;
    clear newL newR
end
Num2 = tempNum;
clear InLen;

%% ------------Stage 3: Confidence interval ------------
%%%------------------------------------------------------------------------
L3 = L2; %% Left data after preprocessing stage 3
R3 = R2; %% Stage 3
[n1, n2] = size(L');
NN = 2000;
rand('state', 22331);
AA = floor(n1*rand(n1, NN))+1;
Resample_L = L(AA);
Resample_R = R(AA);
%size(Resample_L)


temp_m_L = mean(Resample_L);
temp_m_R = mean(Resample_R);

%size(temp_m_L)
m_L = mean(temp_m_L);
v_L = sqrt(n1)* std(temp_m_L);
m_R = mean(temp_m_R) ;
v_R = sqrt(n1)* std(temp_m_R);

InLen = R-L;
m_Len = mean(InLen);
v_Len = std(InLen);

% m_L=mean(L);
% v_L=std(L);
% m_R=mean(R);
% v_R=std(R);

%% K(95%, 90%, n)
%% K(95%, 90%, 1000) =1.709
%% K(95%, 90%, 25) = 2.208
%% K(95%, 90%, 30) = 2.140
%k=2.208;
%k=1.709;
if n<=2
    k = 32.019;
elseif n==3
    k=8.380;
elseif n==4
    k = 5.369;
elseif n==5
    k = 4.275;
elseif n==6
    k = 3.712;
elseif n==7
    k = 3.369;
elseif n == 8
    k = 3.136;
end

if n==9
    k = 2.967;
elseif n ==10
    k = 2.839;
elseif n==11
    k =2.737;
elseif n==12
    k = 2.655;
elseif n==13
    k = 2.587;
elseif n == 14
    k = 2.529;
elseif n == 15
    k = 2.48;
elseif n==16
    k = 2.437;
end

if n==17
    k = 2.4;
elseif n ==18
    k = 2.366;
elseif n==19
    k =2.337;
elseif n==20
    k = 2.31;
end

if n>=21 & n<25
    k = 2.31;
end

if n>=25
    k = 2.208;
end

tempNum = Num2;
Bound = 0.25;
%%---- For left-end point confident interval processing -------------
t =0;
if (v_L > Bound)
    for i=1:n
        if L(i)>=m_L-k*v_L & L(i)<=m_L + k*v_L
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            newInLen(t) = InLen(i);
            tempNum3(t)= tempNum2(i);
        else
            L3(tempNum2(i)) = -1;
            R3(tempNum2(i)) = -1;
        end
    end
    clear L R InLen;
    L = newL;
    R = newR;
    InLen = newInLen;
    n = t;
    tempNum =t;
    clear newL newR newInLen;
else
    tempNum3 = tempNum2;
end

t=0;
if (v_R > Bound)
    for i=1:n
        if (R(i)>= m_R-k*v_R) & (R(i)<= m_R + k*v_R)
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            newInLen(t) = InLen(i);
            tempNum3(t)= tempNum3(i);
        else
            L3(tempNum3(i)) = -1;
            R3(tempNum3(i)) = -1;
        end
    end
    clear L R InLen;
    L = newL;
    R = newR;
    InLen = newInLen;
    n = t;
    tempNum =t;
    clear newL newR newInLen;
end

t=0;
if (v_Len > Bound)
    for i=1:n
        if  (InLen(i)>= m_Len - k*v_Len) & (InLen(i)<= m_Len + k*v_Len)
            t=t+1;
            newL(t) = L(i);
            newR(t) = R(i);
            tempNum3(t)= tempNum3(i);
        else
            L3(tempNum3(i)) = -1;
            R3(tempNum3(i)) = -1;
        end
    end
    clear L R;
    L = newL;
    R = newR;
    n = t;
    tempNum =t;
    clear newL newR
end
Num3 = tempNum;
clear InLen

%%
%%% -------Constraint for reasonable interval data (with overlap)----------------
L4 = L3; %% Left data after preprocessing stage 3
R4 = R3; %% Stage 3
[n1, n2] = size(L');
NN = 2000;
rand('state', 231);
AA = floor(n1*rand(n1, NN))+1;
Resample_L = L(AA);
Resample_R = R(AA);
%size(Resample_L)

temp_m_L = mean(Resample_L);
temp_m_R = mean(Resample_R);

%size(temp_m_L)
m_L = mean(temp_m_L);
v_L = sqrt(n1)* std(temp_m_L);
m_R = mean(temp_m_R) ;
v_R = sqrt(n1)* std(temp_m_R);
% 

%Stat=[m_L, v_L, m_R, v_R];

% m_L = mean(L);
% v_L = std(L);
% m_R = mean(R);
% v_R = std(R);

% disp('m_L, v_L, m_R, v_R')
[m_L, v_L, m_R, v_R];

%--- Method #1
if (v_L+v_R==0)
    %disp('This is For Barrier 1')
    Barrier = (m_L + m_R)/2;
elseif v_L==0
    %disp('This is For Barrier 2')
    Barrier = m_L+ 0.01;
elseif v_R==0
   % disp('This is For Barrier 3')
    Barrier = m_R - 0.01;
else
    %disp('This is For Barrier 4')
    Barrier1 =(-(m_L*v_R^2-m_R*v_L^2) + v_L*v_R*sqrt((m_L-m_R)^2+2*(v_R^2-v_L^2)*log(v_L/v_R)))/(v_L^2-v_R^2);
    Barrier2 =(-(m_L*v_R^2-m_R*v_L^2) - v_L*v_R*sqrt((m_L-m_R)^2+2*(v_R^2-v_L^2)*log(v_L/v_R)))/(v_L^2-v_R^2);
%     Barrier1 =(-(m_L*v_R^2-m_R*v_L^2) + v_L*v_R*sqrt((m_L-m_R)^2-2*(v_R^2-v_L^2)*log(v_L/v_R)))/(v_L^2-v_R^2);
%     Barrier2 =(-(m_L*v_R^2-m_R*v_L^2) - v_L*v_R*sqrt((m_L-m_R)^2-2*(v_R^2-v_L^2)*log(v_L/v_R)))/(v_L^2-v_R^2);
    if (v_L-v_R <=0)&(Barrier1>=m_L) & (Barrier1<=0.5*(m_R+m_L))
        Barrier = Barrier1;
    else
        Barrier = Barrier2;
    end
end


% %-- Method #2
% if (v_L+v_R==0)
%     disp('This is For Barrier 1')
%     Barrier = (m_L + m_R)/2;
% else
%     disp('This is For Barrier 2')
%     Barrier =m_L + (v_L/(v_L+v_R))*(m_R - m_L);
% end

Barrier;

t =0;
for i=1:n
    if (L(i)<= Barrier) & (R(i)>= Barrier)
        t=t+1;
        newL(t) = L(i);
        newR(t) = R(i);
        tempNum4(t)= tempNum3(i);
    else
        L4(tempNum3(i)) = -1;
        R4(tempNum3(i)) = -1;
    end
end
clear L R

L = newL;
R = newR;
n=t;
Num4 =t;
clear newL newR

Data1 = [L1, R1, L2, R2, L3, R3, L4, R4];

L5 = L4; %% Left data after preprocessing stage 4
R5 = R4; %% Stage 4

%% ----------- Admissible region determination ----------
%% -- Table for t_alpha t_0.1
n;
alpha = 0.05;
t_alpha = T_function(n, alpha);

%% --- COmpute the mean of left-end and right-end 
[n1, n2] = size(L');
NN = 2000;
rand('state', 1231);
AA = floor(n1*rand(n1, NN))+1;
Resample_L = L(AA);
Resample_R = R(AA);
%size(Resample_L)

temp_m_L = mean(Resample_L);
temp_m_R = mean(Resample_R);

%size(temp_m_L)
m_L = mean(temp_m_L);
v_L = sqrt(n1)* std(temp_m_L);
m_R = mean(temp_m_R) ;
v_R = sqrt(n1)* std(temp_m_R);

Stat=[m_L, v_L, m_R, v_R];

a = mean(L);
b = mean(R);
v_L = std(L);
v_R = std(R);
%[a, v_L, b, v_R]
C = R - 5.831*L;
D = R - 0.171*L - 8.29;
%[C', D']
%t_alpha = 2.3* t_alpha;
shift1 = t_alpha * std(C)/sqrt(n);
shift2 = t_alpha * std(D)/sqrt(n);

Shift = [shift1, shift2];

% disp('std(C), std(D), shift1, shift2')
[std(C), std(D), shift1, shift2];
% disp('b, a, 5.831*a-shift1, 8.29+0.171*a-shift2')
[b, a, 5.831*a-shift1, 8.29+0.171*a-shift2];
if (b<5.831*a -shift1)&(b<8.29+0.171*a-shift2)
%     disp('Interior FOU')
    clear FS_L FS_R
    j =0;
    for i=1:n
        tempL = 0.5*(L(i)+R(i))-(R(i)-L(i))*sqrt(2)/2;
        tempR = 0.5*(L(i)+R(i))+(R(i)-L(i))*sqrt(2)/2;
        L5(tempNum4(i)) = tempL;
        R5(tempNum4(i)) = tempR;
        
        if (tempL> 0)&(tempR<10)
            j=j+1;
            FS_L(j)= tempL;
            FS_R(j)= tempR;
            FS_C(j)= 0.5*(tempL + tempR);
        end
    end
    n=j;
    %[FS_L; FS_R]    
    L1 = min(FS_L);
    L2 = max(FS_L);
    R1 = min(FS_R);
    R2 = max(FS_R);
    C1 = min(FS_C);
    C2 = max(FS_C);

    temp = (R1-min(FS_C))/(max(FS_C)-L2);
    Apex = (R1+temp*L2)/(1+temp);
    Height = (R1-Apex)/(R1-min(FS_C));
    UMF =[L1, min(FS_C), max(FS_C), R2];
    LMF = [L2, Apex, Apex, R1, Height];

elseif (b>5.831*a-shift1)
%     disp('Left Shoulder')
    clear FS_L FS_R
    j=0;
    for i=1:n
        [tempL, tempR]= leftShoulder(L(i), R(i));
        L5(tempNum4(i)) = tempL;
        R5(tempNum4(i)) = tempR;
        if tempL>=0 & tempR<10
            j=j+1;
            FS_L(j)=tempL;
            FS_R(j)=tempR;
        end
    end
    n=j;
    UMF =[0,  0, max(FS_L), max(FS_R)];
    LMF = [0, 0, min(FS_L), min(FS_R) 1];

else
%     disp('Right Shoulder')
    clear FS_L FS_R
    j=0;
    for i=1:n
        [tempL, tempR]= rightShoulder(L(i), R(i));
        L5(tempNum4(i)) = tempL;
        R5(tempNum4(i)) = tempR;
        if tempL>0 & tempR<=10
            j=j+1;
            FS_L(j)=tempL;
            FS_R(j)=tempR;
        end
    end
    n=j;
    UMF =[min(FS_L), min(FS_R), 10, 10];
    LMF = [max(FS_L), max(FS_R), 10, 10 1];
end


MF=[UMF LMF];

NNum = [Num1, Num2, Num3, Num4, n];
%NNum = [Num1, Num3, Num4, n]

Data2 = [L5, R5];
Data = [Data1, Data2];

%Data = [L1, R1, L2, R2, L3, R3, L4, R4, L5, R5]

