%--------------------------------------------------------------------------
%Copyright 2016, Mengqiu Li

%The following code is licensed under a Creative Commons 
%Attribution-NonCommercial-ShareAlike 4.0 International Public License.

%Under the following terms:

%1) Attribution - You must give appropriate credit, provide a link to the 
%license, and indicate if changes were made. You may do so in any 
%reasonable manner, but not in any way that suggests the licensor 
%endorses you or your use. 

%Preferred method of attribution is citing or referring to the published 
%paper associated with this work: 
%Journal of the American Chemical Society, 2015, 137 (51), pp 16055–16063.

%2) NonCommercial - You may not use the material for commercial purposes.

%3) ShareAlike - If you remix, transform, or build upon the material, 
%you must distribute your contributions under the same license as the 
%original.

%4) No additional restrictions - You may not apply legal terms or 
%technological measures that legally restrict others from doing anything 
%the license permits.

%Under the above terms, you are free to:

%1) Share - copy and redistribute the material in any medium or format
%2) Adapt - remix, transform, and build upon the material

%The license in full: 
%   http://creativecommons.org/licenses/by-nc-sa/4.0/legalcode
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%All the code are written for MATLAB and tested in the following
%environment:
%MATLAB Version:                            8.5.0.197613  (R2015a)
%Bioinformatics Toolbox                     Version 4.5.1 (R2015a)
%Curve Fitting Toolbox                      Version 3.5.1 (R2015a)
%Image Processing Toolbox                   Version 9.2   (R2015a)
%Optimization Toolbox                       Version 7.2   (R2015a)
%Parallel Computing Toolbox                 Version 6.6   (R2015a)
%Statistics and Machine Learning Toolbox    Version 10.0  (R2015a)
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%Contact information:
%   Mengqiu Li (author): mengqiu.li@outlook.com
%   Lars Jeuken (group leader of the project): l.j.c.jeuken@leeds.ac.uk
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%Note for author and internal use:
%Based on fitpumping_9
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%Output data structure:
%
%           Column      1
%           Data        Identifier for the experimental dataset
%
%    		Column		2        	3	        4
%           Data		TL_number	Vesicle_ID	Type (how many lines)
%																		
%One_line	Column		5	6	7	8												
%           Data		a	b	ARS	SD												
%																		
%Two_line	Column		9	10	11	12	13	14	15									
%           Data		a	b	a	b*	TP	ARS	SD									
%																		
%Three_line	Column		16	17	18	19	20	21	22	23	24	25						
%           Data		a	b	a	b*	a	b*	TP1	TP2	ARS	SD						
%																		
%Four_line	Column		26	27	28	29	30	31	32	33	34	35	36	37	38			
%       	Data		a	b	a	b*	a	b*	a	b*	TP1	TP2	TP3	ARS	SD			
%																		
%Five_line	Column		39	40	41	42	43	44	45	46	47	48	49	50
%   		Data		a	b	a	b*	a	b*	a	b*	a	b*	TP1	TP2
%           Column		51	52	53	54
%   		Data		TP3	TP4	ARS	SD

%a, b: are the parameters in the simple line function y=a*x+b
%b*: Not fitted, because the lines are joined
%ARS: adjusted R-square
%SD: standard deviation
%TP: turning points (the time where turning ocurred)			
%--------------------------------------------------------------------------

function FittedParameters=open_data_fit_wp(Exp_ID,TL_No,Traces,Column,Potential_Window,Fitted_PP)
%This function fits multiple joined lines to the data and finds the best
%fitted model. 

%The results from all the fitting are stored in the output, even if the 
%model is eventually abandoned because another model fits better. The best
%model (how many lines) is indicated in Column 4.

%Refer to the paper mentioned in the license for detailed design rationale
%of this fitting procedure.

%Parameters:

%Time resolution for finding turning points (in second)
Resolution=3;

%Time gap between turning points or turning points and the boundary of the
%potential window (in second).
TimeGap=15;

%This criterion determines when to stop fitting data to higher number of
%lines, this should be at least 1
SDCriterion=1.1;

%Prepares data for the the parfor loops
WB=waitbar(0,'Preparing data...');

Traces=Traces(Traces(:,1)==TL_No,:);

Potential_On=Potential_Window(1);
Potential_Off=Potential_Window(2);
Traces=sortrows(Traces,[1,2,3]);
Row_No=size(Traces,1);
Vesicle_No=1;
for i=2:Row_No
    if Traces(i,2)~=Traces(i-1,2)
        Vesicle_No=Vesicle_No+1;
    end
end

S_E=zeros(Vesicle_No,5);
S_E(1,1)=Traces(1,2);
S_E(1,2)=1;
i=2;
l=1;
while i<Row_No+1
    j=i;
    while j<Row_No+1 && Traces(j,2)==Traces(j-1,2)
        j=j+1;
    end
    if j~=Row_No
        j=j-1;
    end
    S_E(l,5)=j;
    k=i;
    while k<j+1
        if Traces(k,3)>=Potential_On
            S_E(l,3)=k-1;
            break;
        else
            k=k+1;
        end
    end
    k=j;
    while k>i-1
        if Traces(k,3)<=Potential_Off
            S_E(l,4)=k;
            break;
        else
            k=k-1;
        end
    end
    if l<Vesicle_No
        l=l+1;
        S_E(l,1)=Traces(j+1,2);
        S_E(l,2)=j+1;
    end
    i=j+2;
end

Options0=optimset('lsqnonlin');
Options=optimset(Options0,'Display','off');

Trace_Length_Max=max(S_E(:,4)-S_E(:,3));
Time_Par=zeros(Trace_Length_Max,1,Vesicle_No);
pH_Par=Time_Par;
Trace_Length_Par=zeros(Vesicle_No,1);

Leaking_C8=zeros(Vesicle_No,1);
Leaking_C9=zeros(Vesicle_No,1);

Para0_Par_L1=zeros(1,2,Vesicle_No);
Para1_Par_L1=[-Inf,-Inf];
Para2_Par_L1=[Inf,Inf];

Para0_Par_L2=zeros(1,3,Vesicle_No);
Para1_Par_L2=[-Inf,-Inf,-Inf];
Para2_Par_L2=[Inf,Inf,Inf];

Para0_Par_L3=zeros(1,4,Vesicle_No);
Para1_Par_L3=[-Inf,-Inf,-Inf,-Inf];
Para2_Par_L3=[Inf,Inf,Inf,Inf];

Para0_Par_L4=zeros(1,5,Vesicle_No);
Para1_Par_L4=[-Inf,-Inf,-Inf,-Inf,-Inf];
Para2_Par_L4=[Inf,Inf,Inf,Inf,Inf];

Para0_Par_L5=zeros(1,6,Vesicle_No);
Para1_Par_L5=[-Inf,-Inf,-Inf,-Inf,-Inf,-Inf];
Para2_Par_L5=[Inf,Inf,Inf,Inf,Inf,Inf];

Para_Par_1=zeros(Vesicle_No,1);Para_Par_2=zeros(Vesicle_No,1);Para_Par_3=zeros(Vesicle_No,1);Para_Par_4=zeros(Vesicle_No,1);Para_Par_5=zeros(Vesicle_No,1);
Para_Par_6=zeros(Vesicle_No,1);Para_Par_7=zeros(Vesicle_No,1);Para_Par_8=zeros(Vesicle_No,1);Para_Par_9=zeros(Vesicle_No,1);Para_Par_10=zeros(Vesicle_No,1);
Para_Par_11=zeros(Vesicle_No,1);Para_Par_12=zeros(Vesicle_No,1);Para_Par_13=zeros(Vesicle_No,1);Para_Par_14=zeros(Vesicle_No,1);Para_Par_15=zeros(Vesicle_No,1);
Para_Par_16=zeros(Vesicle_No,1);Para_Par_17=zeros(Vesicle_No,1);Para_Par_18=zeros(Vesicle_No,1);Para_Par_19=zeros(Vesicle_No,1);Para_Par_20=zeros(Vesicle_No,1);
Para_Par_21=zeros(Vesicle_No,1);Para_Par_22=zeros(Vesicle_No,1);Para_Par_23=zeros(Vesicle_No,1);Para_Par_24=zeros(Vesicle_No,1);Para_Par_25=zeros(Vesicle_No,1);
Para_Par_26=zeros(Vesicle_No,1);Para_Par_27=zeros(Vesicle_No,1);Para_Par_28=zeros(Vesicle_No,1);Para_Par_29=zeros(Vesicle_No,1);Para_Par_30=zeros(Vesicle_No,1);
Para_Par_31=zeros(Vesicle_No,1);Para_Par_32=zeros(Vesicle_No,1);Para_Par_33=zeros(Vesicle_No,1);Para_Par_34=zeros(Vesicle_No,1);Para_Par_35=zeros(Vesicle_No,1);
Para_Par_36=zeros(Vesicle_No,1);Para_Par_37=zeros(Vesicle_No,1);Para_Par_38=zeros(Vesicle_No,1);Para_Par_39=zeros(Vesicle_No,1);Para_Par_40=zeros(Vesicle_No,1);
Para_Par_41=zeros(Vesicle_No,1);Para_Par_42=zeros(Vesicle_No,1);Para_Par_43=zeros(Vesicle_No,1);Para_Par_44=zeros(Vesicle_No,1);Para_Par_45=zeros(Vesicle_No,1);
Para_Par_46=zeros(Vesicle_No,1);Para_Par_47=zeros(Vesicle_No,1);Para_Par_48=zeros(Vesicle_No,1);Para_Par_49=zeros(Vesicle_No,1);Para_Par_50=zeros(Vesicle_No,1);
Para_Par_51=zeros(Vesicle_No,1);Para_Par_52=zeros(Vesicle_No,1);

N_Combination=ceil((Potential_Off-Potential_On-2*TimeGap)/Resolution);
LoopNo_Max_L2=int8(factorial(N_Combination)/(factorial(N_Combination-1)*factorial(1)));
LoopNo_Max_L3=int8(factorial(N_Combination)/(factorial(N_Combination-2)*factorial(2)));
LoopNo_Max_L4=int8(factorial(N_Combination)/(factorial(N_Combination-3)*factorial(3)));
LoopNo_Max_L5=int8(factorial(N_Combination)/(factorial(N_Combination-4)*factorial(4)));

ResNorm1=zeros(1,Vesicle_No);

ResNorm2_1=zeros(LoopNo_Max_L2,Vesicle_No);
ResNorm2_2=zeros(LoopNo_Max_L2,Vesicle_No);
ResNorm2_Final=zeros(1,Vesicle_No);
ResNorm2_1(1,:)=Potential_On+TimeGap;
LoopNo_L2=zeros(Vesicle_No,1);
LoopNo_L2(:,:)=1;


ResNorm3_1=zeros(LoopNo_Max_L3,Vesicle_No);
ResNorm3_2=zeros(LoopNo_Max_L3,Vesicle_No);
ResNorm3_3=zeros(LoopNo_Max_L3,Vesicle_No);
ResNorm3_Final=zeros(1,Vesicle_No);
ResNorm3_1(1,:)=Potential_On+TimeGap;
ResNorm3_2(1,:)=ResNorm3_1(1,:)+TimeGap;
LoopNo_L3=zeros(Vesicle_No,1);
LoopNo_L3(:,:)=1;

ResNorm4_1=zeros(LoopNo_Max_L4,Vesicle_No);
ResNorm4_2=zeros(LoopNo_Max_L4,Vesicle_No);
ResNorm4_3=zeros(LoopNo_Max_L4,Vesicle_No);
ResNorm4_4=zeros(LoopNo_Max_L4,Vesicle_No);
ResNorm4_Final=zeros(1,Vesicle_No);
ResNorm4_1(1,:)=Potential_On+TimeGap;
ResNorm4_2(1,:)=ResNorm4_1(1,:)+TimeGap;
ResNorm4_3(1,:)=ResNorm4_2(1,:)+TimeGap;
LoopNo_L4=zeros(Vesicle_No,1);
LoopNo_L4(:,:)=1;

ResNorm5_1=zeros(LoopNo_Max_L5,Vesicle_No);
ResNorm5_2=zeros(LoopNo_Max_L5,Vesicle_No);
ResNorm5_3=zeros(LoopNo_Max_L5,Vesicle_No);
ResNorm5_4=zeros(LoopNo_Max_L5,Vesicle_No);
ResNorm5_5=zeros(LoopNo_Max_L5,Vesicle_No);
ResNorm5_Final=zeros(1,Vesicle_No);
ResNorm5_1(1,:)=Potential_On+TimeGap;
ResNorm5_2(1,:)=ResNorm5_1(1,:)+TimeGap;
ResNorm5_3(1,:)=ResNorm5_2(1,:)+TimeGap;
ResNorm5_4(1,:)=ResNorm5_3(1,:)+TimeGap;
LoopNo_L5=zeros(Vesicle_No,1);
LoopNo_L5(:,:)=1;

%Setting up the time points table for brute-forcing the fitting of 
%multiple lines
Para_Par_1(:,1)=S_E(:,1);
for i=1:Vesicle_No
    Time=Traces((S_E(i,3)+1):S_E(i,4),3);
    MaxTime=max(Time);
    pH=Traces((S_E(i,3)+1):S_E(i,4),Column);
    Trace_Length_Par(i,1)=size(Time,1);
    Time_Par(1:Trace_Length_Par(i,1),1,i)=Time;
    pH_Par(1:Trace_Length_Par(i,1),1,i)=pH;
    MeanpH=mean(pH);
    Para0_Par_L1(1,2,i)=MeanpH;
    Para0_Par_L2(1,2,i)=MeanpH;
    Para0_Par_L3(1,2,i)=MeanpH;
    Para0_Par_L4(1,2,i)=MeanpH;
    Para0_Par_L5(1,2,i)=MeanpH;
    Row_Leaking=find(Fitted_PP(:,2)==S_E(i,1),1);
    Leaking_C8(i,1)=Fitted_PP(Row_Leaking,8);
    Leaking_C9(i,1)=Fitted_PP(Row_Leaking,9);
    
    while ResNorm2_1(LoopNo_L2(i,1),i)<(MaxTime-TimeGap)
        LoopNo_L2(i,1)=LoopNo_L2(i,1)+1;
        ResNorm2_1(LoopNo_L2(i,1),i)=ResNorm2_1(LoopNo_L2(i,1)-1,i)+Resolution;
    end
    
    while ResNorm3_2(LoopNo_L3(i,1),i)<(MaxTime-TimeGap)
        while ResNorm3_2(LoopNo_L3(i,1),i)<(MaxTime-TimeGap)
            LoopNo_L3(i,1)=LoopNo_L3(i,1)+1;
            ResNorm3_1(LoopNo_L3(i,1),i)=ResNorm3_1(LoopNo_L3(i,1)-1,i);
            ResNorm3_2(LoopNo_L3(i,1),i)=ResNorm3_2(LoopNo_L3(i,1)-1,i)+Resolution;
        end
        ResNorm3_1(LoopNo_L3(i,1),i)=ResNorm3_1(LoopNo_L3(i,1)-1,i)+Resolution;
        ResNorm3_2(LoopNo_L3(i,1),i)=ResNorm3_1(LoopNo_L3(i,1),i)+TimeGap;
    end
    
    while ResNorm4_3(LoopNo_L4(i,1),i)<(MaxTime-TimeGap)
        while ResNorm4_3(LoopNo_L4(i,1),i)<(MaxTime-TimeGap)
            while ResNorm4_3(LoopNo_L4(i,1),i)<(MaxTime-TimeGap)
                LoopNo_L4(i,1)=LoopNo_L4(i,1)+1;
                ResNorm4_1(LoopNo_L4(i,1),i)=ResNorm4_1(LoopNo_L4(i,1)-1,i);
                ResNorm4_2(LoopNo_L4(i,1),i)=ResNorm4_2(LoopNo_L4(i,1)-1,i);
                ResNorm4_3(LoopNo_L4(i,1),i)=ResNorm4_3(LoopNo_L4(i,1)-1,i)+Resolution;
            end
            ResNorm4_2(LoopNo_L4(i,1),i)=ResNorm4_2(LoopNo_L4(i,1)-1,i)+Resolution;
            ResNorm4_3(LoopNo_L4(i,1),i)=ResNorm4_2(LoopNo_L4(i,1),i)+TimeGap;
        end
        ResNorm4_1(LoopNo_L4(i,1),i)=ResNorm4_1(LoopNo_L4(i,1)-1,i)+Resolution;
        ResNorm4_2(LoopNo_L4(i,1),i)=ResNorm4_1(LoopNo_L4(i,1),i)+TimeGap;
        ResNorm4_3(LoopNo_L4(i,1),i)=ResNorm4_2(LoopNo_L4(i,1),i)+TimeGap;
    end
    
    while ResNorm5_4(LoopNo_L5(i,1),i)<(MaxTime-TimeGap)
        while ResNorm5_4(LoopNo_L5(i,1),i)<(MaxTime-TimeGap)
            while ResNorm5_4(LoopNo_L5(i,1),i)<(MaxTime-TimeGap)
                while ResNorm5_4(LoopNo_L5(i,1),i)<(MaxTime-TimeGap)
                    LoopNo_L5(i,1)=LoopNo_L5(i,1)+1;
                    ResNorm5_1(LoopNo_L5(i,1),i)=ResNorm5_1(LoopNo_L5(i,1)-1,i);
                    ResNorm5_2(LoopNo_L5(i,1),i)=ResNorm5_2(LoopNo_L5(i,1)-1,i);
                    ResNorm5_3(LoopNo_L5(i,1),i)=ResNorm5_3(LoopNo_L5(i,1)-1,i);
                    ResNorm5_4(LoopNo_L5(i,1),i)=ResNorm5_4(LoopNo_L5(i,1)-1,i)+Resolution;
                end
                ResNorm5_3(LoopNo_L5(i,1),i)=ResNorm5_3(LoopNo_L5(i,1)-1,i)+Resolution;
                ResNorm5_4(LoopNo_L5(i,1),i)=ResNorm5_3(LoopNo_L5(i,1),i)+TimeGap;
            end
            ResNorm5_2(LoopNo_L5(i,1),i)=ResNorm5_2(LoopNo_L5(i,1)-1,i)+Resolution;
            ResNorm5_3(LoopNo_L5(i,1),i)=ResNorm5_2(LoopNo_L5(i,1),i)+TimeGap;
            ResNorm5_4(LoopNo_L5(i,1),i)=ResNorm5_3(LoopNo_L5(i,1),i)+TimeGap;
        end
        ResNorm5_1(LoopNo_L5(i,1),i)=ResNorm5_1(LoopNo_L5(i,1)-1,1)+Resolution;
        ResNorm5_2(LoopNo_L5(i,1),i)=ResNorm5_1(LoopNo_L5(i,1),i)+TimeGap;
        ResNorm5_3(LoopNo_L5(i,1),i)=ResNorm5_2(LoopNo_L5(i,1),i)+TimeGap;
        ResNorm5_4(LoopNo_L5(i,1),i)=ResNorm5_3(LoopNo_L5(i,1),i)+TimeGap;
    end
end
LoopNo_L2=LoopNo_L2-1;
LoopNo_L3=LoopNo_L3-1;
LoopNo_L4=LoopNo_L4-1;
LoopNo_L5=LoopNo_L5-1;

%The parfor loops for fitting

waitbar(0/5,WB,'Fitting WP sections: One-line Model');
%fit to the 1-line model
parfor i=1:Vesicle_No
    [Temp_Para_L1,ResNorm1(1,i)]=lsqnonlin(@(Para) oneline_fitting(Para,Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L1(:,:,i),Para1_Par_L1,Para2_Par_L1,Options);
    Para_Par_3(i,1)=Temp_Para_L1(1);
    Para_Par_4(i,1)=Temp_Para_L1(2);
    Para_Par_5(i,1)=adjusted_r_squared_2(pH_Par(:,:,i),Trace_Length_Par(i,1),2,ResNorm1(1,i));
    Para_Par_6(i,1)=sqrt(ResNorm1(1,i)/(Trace_Length_Par(i,1)-1));
    if Para_Par_6(i,1)<SDCriterion*Leaking_C8(i,1) || Para_Par_6(i,1)<SDCriterion*Leaking_C9(i,1)
        Para_Par_2(i,1)=1;
    end
end

waitbar(1/5,WB,'Fitting WP sections: Two-line Model');
%fit to the 2-line model
for i=1:Vesicle_No    
    if Para_Par_2(i,1)==0 && LoopNo_L2(i,1)>=1
        parfor n2=1:LoopNo_L2(i,1)
            [~,ResNorm2_2(n2,i)]=lsqnonlin(@(Para) twoline_fitting(Para,ResNorm2_1(n2,i),Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L2(:,:,i),Para1_Par_L2,Para2_Par_L2,Options);
        end
        Temp_Row=find(ResNorm2_2(1:LoopNo_L2(i,1),i)==min(ResNorm2_2(1:LoopNo_L2(i,1),i)),1);
        ResNorm2_Final(1,i)=ResNorm2_2(Temp_Row,i);
        Para_Par_11(i,1)=ResNorm2_1(Temp_Row,i);
        [Temp_Para_L2,~]=lsqnonlin(@(Para) twoline_fitting(Para,Para_Par_11(i,1),Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L2(:,:,i),Para1_Par_L2,Para2_Par_L2,Options);
        Para_Par_7(i,1)=Temp_Para_L2(1);
        Para_Par_8(i,1)=Temp_Para_L2(2);
        Para_Par_9(i,1)=Temp_Para_L2(3);
        Para_Par_10(i,1)=(Para_Par_7(i,1)-Para_Par_9(i,1))*Para_Par_11(i,1)+Para_Par_8(i,1);
        Para_Par_12(i,1)=adjusted_r_squared_2(pH_Par(:,:,i),Trace_Length_Par(i,1),4,ResNorm2_Final(1,i));
        Para_Par_13(i,1)=sqrt(ResNorm2_Final(1,i)/(Trace_Length_Par(i,1)-1));
        if Para_Par_13(i,1)<SDCriterion*Leaking_C8(i,1) || Para_Par_13(i,1)<SDCriterion*Leaking_C9(i,1)
            ARSes=[Para_Par_5(i,1),Para_Par_12(i,1)];
            BestFit=find(ARSes==max(ARSes),1);
            Para_Par_2(i,1)=BestFit;
        end
    end
    if Para_Par_2(i,1)==0 && LoopNo_L2(i,1)<1
        Para_Par_2(i,1)=1;
    end
end

waitbar(2/5,WB,'Fitting WP sections: Three-line Model');
%fit to the 3-line model
for i=1:Vesicle_No    
    if Para_Par_2(i,1)==0 && LoopNo_L3(i,1)>=1
        parfor n3=1:LoopNo_L3(i,1)
            [~,ResNorm3_3(n3,i)]=lsqnonlin(@(Para) threeline_fitting(Para,[ResNorm3_1(n3,i),ResNorm3_2(n3,i)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L3(:,:,i),Para1_Par_L3,Para2_Par_L3,Options);
        end
        Temp_Row=find(ResNorm3_3(1:LoopNo_L3(i,1),i)==min(ResNorm3_3(1:LoopNo_L3(i,1),i)),1);
        ResNorm3_Final(1,i)=ResNorm3_3(Temp_Row,i);
        Para_Par_20(i,1)=ResNorm3_1(Temp_Row,i);
        Para_Par_21(i,1)=ResNorm3_2(Temp_Row,i);
        [Temp_Para_L3,~]=lsqnonlin(@(Para) threeline_fitting(Para,[Para_Par_20(i,1),Para_Par_21(i,1)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L3(:,:,i),Para1_Par_L3,Para2_Par_L3,Options);
        Para_Par_14(i,1)=Temp_Para_L3(1);
        Para_Par_15(i,1)=Temp_Para_L3(2);
        Para_Par_16(i,1)=Temp_Para_L3(3);
        Para_Par_18(i,1)=Temp_Para_L3(4);
        Para_Par_17(i,1)=(Para_Par_14(i,1)-Para_Par_16(i,1))*Para_Par_20(i,1)+Para_Par_15(i,1);
        Para_Par_19(i,1)=(Para_Par_16(i,1)-Para_Par_18(i,1))*Para_Par_21(i,1)+Para_Par_17(i,1);
        Para_Par_22(i,1)=adjusted_r_squared_2(pH_Par(:,:,i),Trace_Length_Par(i,1),6,ResNorm3_Final(1,i));
        Para_Par_23(i,1)=sqrt(ResNorm3_Final(1,i)/(Trace_Length_Par(i,1)-1));
        if Para_Par_23(i,1)<SDCriterion*Leaking_C8(i,1) || Para_Par_23(i,1)<SDCriterion*Leaking_C9(i,1)
            ARSes=[Para_Par_5(i,1),Para_Par_12(i,1),Para_Par_22(i,1)];
            BestFit=find(ARSes==max(ARSes),1);
            Para_Par_2(i,1)=BestFit;
        end
    end
    if Para_Par_2(i,1)==0 && LoopNo_L3(i,1)<1
        ARSes=[Para_Par_5(i,1),Para_Par_12(i,1)];
        BestFit=find(ARSes==max(ARSes),1);
        Para_Par_2(i,1)=BestFit;
    end
end

waitbar(3/5,WB,'Fitting WP sections: Four-line Model');
%fit to the 4-line model
for i=1:Vesicle_No    
    if Para_Par_2(i,1)==0 && LoopNo_L4(i,1)>=1
        parfor n4=1:LoopNo_L4(i,1)
            [~,ResNorm4_4(n4,i)]=lsqnonlin(@(Para) fourline_fitting(Para,[ResNorm4_1(n4,i),ResNorm4_2(n4,i),ResNorm4_3(n4,i)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L4(:,:,i),Para1_Par_L4,Para2_Par_L4,Options);
        end
        Temp_Row=find(ResNorm4_4(1:LoopNo_L4(i,1),i)==min(ResNorm4_4(1:LoopNo_L4(i,1),i)),1);
        ResNorm4_Final(1,i)=ResNorm4_4(Temp_Row,i);
        Para_Par_32(i,1)=ResNorm4_1(Temp_Row,i);
        Para_Par_33(i,1)=ResNorm4_2(Temp_Row,i);
        Para_Par_34(i,1)=ResNorm4_3(Temp_Row,i);
        [Temp_Para_L4,~]=lsqnonlin(@(Para) fourline_fitting(Para,[Para_Par_32(i,1),Para_Par_33(i,1),Para_Par_34(i,1)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L4(:,:,i),Para1_Par_L4,Para2_Par_L4,Options);
        Para_Par_24(i,1)=Temp_Para_L4(1);
        Para_Par_25(i,1)=Temp_Para_L4(2);
        Para_Par_26(i,1)=Temp_Para_L4(3);
        Para_Par_28(i,1)=Temp_Para_L4(4);
        Para_Par_30(i,1)=Temp_Para_L4(5);
        Para_Par_27(i,1)=(Para_Par_24(i,1)-Para_Par_26(i,1))*Para_Par_32(i,1)+Para_Par_25(i,1);
        Para_Par_29(i,1)=(Para_Par_26(i,1)-Para_Par_28(i,1))*Para_Par_33(i,1)+Para_Par_27(i,1);
        Para_Par_31(i,1)=(Para_Par_28(i,1)-Para_Par_30(i,1))*Para_Par_34(i,1)+Para_Par_29(i,1);
        Para_Par_35(i,1)=adjusted_r_squared_2(pH_Par(:,:,i),Trace_Length_Par(i,1),8,ResNorm4_Final(1,i));
        Para_Par_36(i,1)=sqrt(ResNorm4_Final(1,i)/(Trace_Length_Par(i,1)-1));
        if Para_Par_36(i,1)<SDCriterion*Leaking_C8(i,1) || Para_Par_36(i,1)<SDCriterion*Leaking_C9(i,1)
            ARSes=[Para_Par_5(i,1),Para_Par_12(i,1),Para_Par_22(i,1),Para_Par_35(i,1)];
            BestFit=find(ARSes==max(ARSes),1);
            Para_Par_2(i,1)=BestFit;
        end
    end
    if Para_Par_2(i,1)==0 && LoopNo_L4(i,1)<1
        ARSes=[Para_Par_5(i,1),Para_Par_12(i,1),Para_Par_22(i,1)];
        BestFit=find(ARSes==max(ARSes),1);
        Para_Par_2(i,1)=BestFit;
    end
end

waitbar(4/5,WB,'Fitting WP sections: Five-line Model');
%fit to the 5-line model
for i=1:Vesicle_No   
    if Para_Par_2(i,1)==0 && LoopNo_L5(i,1)>=1
        parfor n5=1:LoopNo_L5(i,1)
            [~,ResNorm5_5(n5,i)]=lsqnonlin(@(Para) fiveline_fitting(Para,[ResNorm5_1(n5,i),ResNorm5_2(n5,i),ResNorm5_3(n5,i),ResNorm5_4(n5,i)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L5(:,:,i),Para1_Par_L5,Para2_Par_L5,Options);
        end
        Temp_Row=find(ResNorm5_5(1:LoopNo_L5(i,1),i)==min(ResNorm5_5(1:LoopNo_L5(i,1),i)),1);
        ResNorm5_Final(1,i)=ResNorm5_5(Temp_Row,i);
        Para_Par_47(i,1)=ResNorm5_1(Temp_Row,i);
        Para_Par_48(i,1)=ResNorm5_2(Temp_Row,i);
        Para_Par_49(i,1)=ResNorm5_3(Temp_Row,i);
        Para_Par_50(i,1)=ResNorm5_4(Temp_Row,i);
        [Temp_Para_L5,~]=lsqnonlin(@(Para) fiveline_fitting(Para,[Para_Par_47(i,1),Para_Par_48(i,1),Para_Par_49(i,1),Para_Par_50(i,1)],Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),Para0_Par_L5(:,:,i),Para1_Par_L5,Para2_Par_L5,Options);
        Para_Par_37(i,1)=Temp_Para_L5(1);
        Para_Par_38(i,1)=Temp_Para_L5(2);
        Para_Par_39(i,1)=Temp_Para_L5(3);
        Para_Par_41(i,1)=Temp_Para_L5(4);
        Para_Par_43(i,1)=Temp_Para_L5(5);
        Para_Par_45(i,1)=Temp_Para_L5(6);
        Para_Par_40(i,1)=(Para_Par_37(i,1)-Para_Par_39(i,1))*Para_Par_47(i,1)+Para_Par_38(i,1);
        Para_Par_42(i,1)=(Para_Par_39(i,1)-Para_Par_41(i,1))*Para_Par_48(i,1)+Para_Par_40(i,1);
        Para_Par_44(i,1)=(Para_Par_41(i,1)-Para_Par_43(i,1))*Para_Par_49(i,1)+Para_Par_42(i,1);
        Para_Par_46(i,1)=(Para_Par_43(i,1)-Para_Par_45(i,1))*Para_Par_50(i,1)+Para_Par_44(i,1);
        Para_Par_51(i,1)=adjusted_r_squared_2(pH_Par(:,:,i),Trace_Length_Par(i,1),10,ResNorm5_Final(1,i));
        Para_Par_52(i,1)=sqrt(ResNorm5_Final(1,i)/(Trace_Length_Par(i,1)-1));
        %there is not IF function here to check the SD because five-line is
        %the last model, if further models are added, add the IF check here
        %if Para_Par_52(i,1)<SDCriterion*Leaking_C7(i,1) || Para_Par_52(i,1)<SDCriterion*Leaking_C8(i,1)
            ARSes=[Para_Par_5(i,1),Para_Par_12(i,1),Para_Par_22(i,1),Para_Par_35(i,1),Para_Par_51(i,1)];
            BestFit=find(ARSes==max(ARSes),1);
            Para_Par_2(i,1)=BestFit;
        %end
    end
    if Para_Par_2(i,1)==0 && LoopNo_L5(i,1)<1
        ARSes=[Para_Par_5(i,1),Para_Par_12(i,1),Para_Par_22(i,1),Para_Par_35(i,1)];
        BestFit=find(ARSes==max(ARSes),1);
        Para_Par_2(i,1)=BestFit;
    end
end

FittedParameters=[Para_Par_1,Para_Par_2,Para_Par_3,Para_Par_4,Para_Par_5,Para_Par_6,Para_Par_7,Para_Par_8,Para_Par_9,Para_Par_10,...
    Para_Par_11,Para_Par_12,Para_Par_13,Para_Par_14,Para_Par_15,Para_Par_16,Para_Par_17,Para_Par_18,Para_Par_19,Para_Par_20,...
    Para_Par_21,Para_Par_22,Para_Par_23,Para_Par_24,Para_Par_25,Para_Par_26,Para_Par_27,Para_Par_28,Para_Par_29,Para_Par_30,...
    Para_Par_31,Para_Par_32,Para_Par_33,Para_Par_34,Para_Par_35,Para_Par_36,Para_Par_37,Para_Par_38,Para_Par_39,Para_Par_40,...
    Para_Par_41,Para_Par_42,Para_Par_43,Para_Par_44,Para_Par_45,Para_Par_46,Para_Par_47,Para_Par_48,Para_Par_49,Para_Par_50,...
    Para_Par_51,Para_Par_52];

close(WB);

%add TL_No and Experimental identifier to the results, this is added in the
%end because this function is built from the legacy versions.
Temp=zeros(size(FittedParameters,1),(2+size(FittedParameters,2)));
Temp(:,3:end)=FittedParameters;
Temp(:,1)=Exp_ID;
Temp(:,2)=TL_No;
FittedParameters=Temp;
end