%--------------------------------------------------------------------------
%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 fitleaking_3
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%Data structure of 'FittedParameters'
%Column 1: Experimental identifier
%Column 2: TL number
%Column 3: Vesicle ID
%Column 4: pH0
%Column 5: Time0 (just the PotentialOff value!)
%Column 6: amplitude
%Column 7: lambda
%Column 8: resnorm form lsqnonlin
%          fit function: pH0+amplitude*exp(-(Time-Time0)/lambda)
%          half-life = -lambda*ln(0.5)
%Column 9: Standard Deviation in the pre-potential window
%Column 10:Standard Deviation between the data points and the fitted
%          exponential line in the post-potential window
%--------------------------------------------------------------------------

function FittedParameters=open_data_fit_pp(Exp_ID,TL_No,Traces,Column,Potential_Window)
%This function fits exponential function to the section of traces after the
%potential was turned off.

%Refer to the paper mentioned in the license for details.

%Parameters:
%This number of points at the beginning of the post-potential section will
%be averaged for the parameter initialisation
No_Mean=3;

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

Potential_On=Potential_Window(1);
Potential_Off=Potential_Window(2);
Traces_Sorted=sortrows(Traces,[1,2,3]);

%Identify how many vesicles are in the data
Row_No=size(Traces_Sorted,1);
Vesicle_No=1;
for i=2:Row_No
    if Traces_Sorted(i,2)~=Traces_Sorted(i-1,2)
        Vesicle_No=Vesicle_No+1;
    end
end

%Identify the location of each vesilce's traces in the data (the length of
%each vesicle's traces may be different because if the HPTS ratio is
%outside the boundary of the HPTS-to-pH conversion calibration function,
%the data point is discarded entirely.
S_E=zeros(Vesicle_No,5);
S_E(1,1)=Traces_Sorted(1,2);
S_E(1,2)=1;
i=2;
l=1;
while i<Row_No+1
    j=i;
    while j<Row_No+1 && Traces_Sorted(j,2)==Traces_Sorted(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_Sorted(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_Sorted(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_Sorted(j+1,2);
        S_E(l,2)=j+1;
    end
    i=j+2;
end

%Prepare temporary variables for the parfor loop
Trace_Length_Max=max(S_E(:,4)-S_E(:,3));
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);

Time_Par=zeros(Trace_Length_Max,1,Vesicle_No);
pH_Par=Time_Par;
Trace_Length_Par=zeros(Vesicle_No,1);

Para_Par_1(:,1)=S_E(:,1);
Para_Par_3(:,1)=Potential_Off;
Options0=optimset('lsqnonlin');
Options=optimset(Options0,'Display','off');
Para0=zeros(1,3,Vesicle_No);
Para1=[-Inf,1e-20,-Inf];
Para2=[Inf,Inf,Inf];

%Prepare data for the parfor loop, again, the variable trace length has to
%be dealt with to enable the use of the parfor loop
for i=1:Vesicle_No
    pH0=mean(Traces_Sorted(S_E(i,2):S_E(i,3),Column));
    Para_Par_7(i,1)=std(Traces_Sorted(S_E(i,2):S_E(i,3),Column));
    
    Time=Traces_Sorted((S_E(i,4)+1):S_E(i,5),3);
    pH=Traces_Sorted((S_E(i,4)+1):S_E(i,5),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;
    
    Para0(1,1,i)=mean(pH(1:No_Mean))-pH0;
    pH_End=mean(pH((end-No_Mean):end));
    Para0(1,2,i)=(Time(end)-Potential_Off)/log(abs(Para0(1,1,i)/(pH_End-pH0)));
    
    if Para0(1,2,i)<=0
            Para0(1,2,i)=(Time(end)-Potential_Off);
    end
    
    Para0(1,3,i)=pH0;
end

%fitting the post-potential section to an expotential function, the 
%variable trace length is dealt with by the exp_leaking_obj_4
parfor i=1:Vesicle_No
    [Temp_Para,Para_Par_6(i,1)]=lsqnonlin(@(Para)...
        exp_leaking_obj_4(Para,Para_Par_3(i,1),...
        Time_Par(:,:,i),pH_Par(:,:,i),Trace_Length_Par(i,1)),...
        Para0(:,:,i),Para1,Para2,Options);
    Para_Par_4(i,1)=Temp_Para(1);
    Para_Par_5(i,1)=Temp_Para(2);
    Para_Par_2(i,1)=Temp_Para(3);
end

Para_Par_8(:,1)=sqrt(Para_Par_6(:,1)./(Trace_Length_Par(:,1)-1));

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];

%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