% Single Curve Model Fitting Fitting 
% By Alexei Valiaev  October 2006
% Please send comments or bug reports to Alexei Valiaev at av12 at duke dot edu

% dirname - the name of the directory where you store the force curves 
% par_i - the initiatial values for the parameters of the FJC model
% T - is the value of the room temperature
% HF - the maximum force value estimate used in the FJC fit 
% LF- the minimum force value estimate used in the FJC fit 
% Since the points on the force curves are discrete values for HF and LF will be adjusted for each force curve. See AHF and ALF description
% BN - number of the bootstrap samples used in the fitting
% i - is just an index which shows the number of the current file
% kuhn_BN - vector which shows the value of the kuhn length for each bootstrap sample for an individual curve
% kuhn_BN_A - is a matrix which shows the value of the kuhn length for each bootstrap sample and for does this for each force curve in the set
% AHF - an actual value of the maximum force used in the fitting (it willbe close to value of HF)
% ALF- - an actual value of the minimum force used in the fitting (it will be close to value of LF)
% F - is the force vector from your force curve
% S - is the distance vector from your force curve
% FS- force in the selected region of the force curve i.e., the segment between the minimum and maximum force
% SS - separation in the selected region of the force curve
% fwdrag - is the force without with no drag contribution subtracted.
% NP - number of points in your force curve 
% ukb - uncertainty in the estimated Kuhn length calculated from the bootstrap samples
% ucb  - uncertainty in the estimated Contour length calculated from the bootstrap samples
% mkb - mean value of the Kuhn length from the bootstrap samples 
% ukcb - is a matrix with values  ukb and ucb
% pFJC - the parameters obtained from FJC fit to the selected region of the force curve
% ApFJC - the parameters obtained from FJC fit to the selected region for the actual data (i.e, before bootstrapping 
% f_fit - linearly spaced vector for the force 
% Z_fit - FJC model fit data
% f_fitM and Z_fitM - collects the fit data for all curves into a matrix
% fname - is the number of your file i.e  1,2 or 12 etc. The program assumes that you number the files at the end of the filename 
% points - the matrix of the point indexes correponding for minimum and maximum forces for the set of your force curves 
% nb - number of the boostrap samples 
% ww_matr - Matrix of the separation residuals which are obtained by iterative reighting procedure
% sep_perms - matrix of the bootstrap samples created from the separation data
% force_perms - matrix of the bootstrap samples created from the force data
% delta_pFJC - a difference between the Kuhn lengthts of the two adjacent  iteractions when the interactive process will be terminated 
% sigma_square and sigma_square_hat - scale parameters to be estimated 
% ww - weigths calculated from sepration residuals 
% Vd_f - is the variance function 

% FUNCTIONS
% find_points - uses the moving window approach to find the minimum and maximum force  defined as AHF and ALF
% find_drag - the program subtracts the drag force contribution
% bootstrap_calc - performs data bootstrapping and estimates the FJC parameters from the bootstrap samples 
% we fit the data to the FJC model using interative procedure.
% At each iteration we use the new weigths for the separation residuals
% find_file_name - this function mataches each set of the FJC parameters with the corresponding file name that you have in your dirname. 
% bootstrp - performs the boostrapping of the data. See paper references for details 
% fjcmin_noweight - performs FJC parameter estimation. There are no weights used in this function. This is our first step 
% fjcmin_weight - performs FJC parameter estimation. No weights are used in the minimization
% fjc - calcualtes the FJC distance based on the FJC fit parameters

% The program saves:
%a) the separation and force vectors of the FJC fits for each force curve
%b) The Kuhn lengths obtained from the bootstrap samples for each curve
%c) The uncertainty in the calculated Kuhn length and Contour length using the bootstrap approach  

function built_force_curve
close all
clear all

dirname='D:\LAPTOP\DISK_C\DIWORK\2003\15March_2003\test';
curdir = pwd;
files = dir(dirname);
numfiles = length(files) - 2;

par_i =[0.25 200]; 
T=25;
HF=198;
LF=60;
BN=3;
indL=7;
i=1;
while i<=numfiles
    if i==1;
        kuhn_BN=zeros(1,1);
    end
    %%%%%%%%% LOAD FILES %%%%%%%%%
    % You will need to change the two lines with eval function based on the
    % names for your force curves. For example my files are called something like forSepav150303_000, forSepav150303_001 etc... 
    filename = files(i+2).name;
    eval(['cd ',dirname]); 
    load(filename);
    eval(['cd ',curdir]);
    endmat=length(filename)-4; 
    eval(['S=SepNew',filename(indL:endmat),';']);
    eval(['F=force',filename(indL:endmat),';']);
    %%%%%%%%% END LOAD FILES %%%%%%%%%
    NP=length(S);
    
    [x,AHF,ALF]=find_points(F,HF,LF,S);
    [F,fwdrag,FS,dF,S_Fit] = find_drag(x,S,F,NP);
    qq=1;
    Y=zeros(qq,3);
    U=zeros(qq,2);    
    G=zeros(qq,2);
    
    k=1;
    FS = F(x(1):x(2));
    SS = S(x(1):x(2));
    [ukb,ucb,sigma_square,Vd_f,pFJC,kuhn_BN]=bootstrap_calc(SS,FS,T,par_i,BN,i,kuhn_BN);
    
    ApFJC(1,:)=pFJC; 
    %%%%%%%%%%%%%%Fit the data%%%%%%%%%%%%%%%%
    f_fit=linspace(0,-FS(end),length(FS))';
    Z_fit = fjc(f_fit,pFJC,T);
    
    if i==1;
        ukcb=zeros(1,2);
        ukcb=[ukcb; ukb ucb];
        Z_fitM=zeros(length(F),1);  
        f_fitM=zeros(length(F),1);  
        f_fitM(:,i)=[zeros(x(1)-1,1); f_fit; zeros(length(S)-x(2),1)];
        Z_fitM(:,i)=[zeros(x(1)-1,1); Z_fit; zeros(length(S)-x(2),1)];
    else  
        ukcb=[ukcb; ukb ucb];
        Z_fitM(:,i)=[zeros(x(1)-1,1); Z_fit; zeros(length(S)-x(2),1)];
        f_fitM(:,i)=[zeros(x(1)-1,1); f_fit; zeros(length(S)-x(2),1)];
    end    
    fname=find_file_name(k,filename);
    if i==1
        ApFJC_M=zeros(1,3);
        points=zeros(1,2);
        ApFJC_M=[ApFJC_M; ApFJC fname];
        points=[points;x'];
    else 
        ApFJC_M=[ApFJC_M; ApFJC fname];
        points=[points;x'];
    end   
    i=i+1;    
end 
end

%%  Build the box plot with the Kuhn lengths obtained from the bootstrap samples and 95% confidence intervals %%%%
ukcb(1,:)=[];
ApFJC_M(1,:)=[];
points(1,:)=[];
[nb numfiles] = size(kuhn_BN);
Grid=repmat(1:numfiles,[nb,1]);
ukb=ukcb(:,1);

figure(1)
hold on
a=1;
kuhn_BN_A=zeros(size(kuhn_BN,1),1);
for i=1:numfiles
    conf_int_max=(ApFJC_M(a,1)+2*ukb(a)); 
    conf_int_min=(ApFJC_M(a,1)-2*ukb(a));
    % now we will draw the short red lines indicating the conf. intrevals 
    line_conf_int=[a-0.2,a+0.2];
    plot(Grid(:,a),kuhn_BN(:,i),'b*',line_conf_int,conf_int_min*ones(1,2),'r',line_conf_int,conf_int_max*ones(1,2),'r');
    kuhn_BN_A=[kuhn_BN_A kuhn_BN(:,i)];
    a=a+1;
end
kuhn_BN_A(:,1)=[];
hold off
set(gca,'FontSize',24);
xlabel('Curve','FontSize',20,'Color',[1 0.1 0.3])
ylabel('Kuhn Length [nm]','FontSize',20,'Color',[1 0.1 0.3])   

%%%%%%%%%%%%%%Save the data in mat files to the directory par%%%%%%%%%%%
convertdir = 'par'; 
mkdir(dirname,convertdir)
savename=[dirname,'\',convertdir,'\','par',filename(1:(end-8))]; 
pFJCname=['pFJC',filename(9:(endmat-4))];
sep_r_name=['sep_r',filename(9:(endmat-4))];   
Pname=['points',filename(9:(endmat-4))];
Z_fitM_name=['FJC_fit_dM',filename(9:(endmat-4))];
f_fitM_name=['FJC_fit_fM',filename(9:(endmat-4))];
ukcb_name=['ukcb',filename(9:(endmat-4))];
kuhn_BN_name=['kuhn_BN',filename(9:(endmat-4))];

eval([pFJCname,' = ApFJC_M;']);
eval([Pname,' = points;']);
eval([Z_fitM_name,' = Z_fitM;']);
eval([f_fitM_name,' = f_fitM;']);
eval([kuhn_BN_name,' = kuhn_BN;']);
eval([ukcb_name,' = ukcb;']);
eval(['save(savename,''',pFJCname,''',''',Pname,''',''',Z_fitM_name,''',''',f_fitM_name,''',''',ukcb_name,''',''',kuhn_BN_name,''')'])
disp('Finished')


function [ukb,ucb,sigma_square,Vd_f,pFJC,kuhn_BN,sep_r]=bootstrap_calc(SS,FS,T,par_i,BN,i,kuhn_BN);
[bootstat,bootsam] = bootstrp(BN,'mean',SS); 
sep_perms=SS(bootsam);
force_perms=FS(bootsam);
sep_perms_sort=sort(sep_perms);
force_perms_sort=sort(force_perms);

pFJC_M=zeros(1,2);
pb=zeros(1,2);
sigma_square_matr=zeros(1,1);
ww_matr=zeros(length(FS),1);

% Here we perform wieghted least squares minimization to find the FJC fit
% parameters for both original data and the bootstrap data 
rr=1;
while rr<=BN
    %%%%%%%%%%
    ss=1;
    delta_pFJC=1;
    while ss<=50 
        if  abs(delta_pFJC)<=0.01
            break;
        else 
            if ss==1;
                [p,fval,exitflag,output] = fminsearch(@fjcmin_noweight,par_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),sep_perms(:,rr),-force_perms(:,rr),T);
            else
                [p,fval,exitflag,output] = fminsearch(@fjcmin_weight,par_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),sep_perms(:,rr),-force_perms(:,rr),T,ww,sigma_square);
                delta_pFJC=pFJC_M(ss)-pFJC_M(ss-1);    
            end
            pFJC_M=[pFJC_M; p];            
            [sep_r,ww,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(p,force_perms(:,rr),sep_perms(:,rr));
        end 
        ss=ss+1;
    end
    %%%%%%%%%%%%%pb will keep the values of the parameters only on the last iteration
    pb=[pb; p];
    kuhn_BN(rr,i)=p(1);
    rr=rr+1;
end            

ss=1;
delta_pFJC=1;
while ss<=50 
    if  abs(delta_pFJC)<=0.01
        break;
    else 
        if ss==1;
            [pFJC,fval,exitflag,output] = fminsearch(@fjcmin_noweight,par_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,T);
        else
            [pFJC,fval,exitflag,output] = fminsearch(@fjcmin_weight,par_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,T,ww,sigma_square);
            delta_pFJC=pFJC_M(ss)-pFJC_M(ss-1);    
        end
        [sep_r,ww,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(pFJC,FS,SS);
        pFJC_M=[pFJC_M; pFJC];     
        ww_matr=[ww_matr,ww];
        sigma_square_matr=[sigma_square_matr,sigma_square]; 
    end 
    ss=ss+1;
end
pFJC_M(1,:)=[];
pb(1,:)=[];
%%%%%%Estimate the uncertaintity using bootstrap samples%%%
mkb=sum(pb(:,1))/BN;
contour_mean_bootstrap=sum(pb(:,2))/BN;
for pp=1:BN
    BB_kuhn(pp)=(pb(pp,1)-mkb)^2;
    BB_contour(pp)=(pb(pp,2)-contour_mean_bootstrap)^2;
end 
ukb=((1/(BN-1))*sum(BB_kuhn))^0.5;
ucb=((1/(BN-1))*sum(BB_contour))^0.5; 

function [sep_r,ww,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(pFJC,FS,SS,ww);
T=25;
KT = 4.1*((T+273)/298);
f=-FS;
z = pFJC(2)*(coth(f.*pFJC(1)/KT)-KT./(f*pFJC(1)));
sep_r=SS-z; 

%%Weight sep residuals
%   # beta
b=pFJC(1)*(-FS)./KT;
%  # Langevin(beta)
Lng_b=(coth(b)-1./b);
Vd_f=pFJC(2)*pFJC(1)*(1-Lng_b.^2-2*Lng_b./b);
sigma_square_hat=(1/(length(FS)-2))*sum(sep_r.^2);
sigma_square=(1/(length(FS)-2))*sum(sep_r.^2./Vd_f);
ww=(1./Vd_f);

function [x,AHF,ALF]=find_points(F,HF,LF,S)
% find_points - uses the moving window approach to find the minimum and maximum force defined as AHF and ALF
% ws - the window size for the moving average filter. The size of this window should be adjusted based on the noise in your data
% wi - window increment 
% fm- the mean value of force in the window
% x - is the vector that contains the index of the points corresponding to minimum and maximum force
x=zeros(2,1);
F=-F;
%Here we will find the high force point 
[qq,ww]=find(F>=HF & S>=40);
find_force=F(qq);
x(2)=qq(1);
AHF=find_force(1);
%Here we will find the low force point
ws=10;
wi=1;
ii=x(2);
while ii>0
    fm=mean(F((ii-ws):ii));
    if fm>LF
        ii=ii-wi;
        continue;
    else
        ALF=F(ii-wi);
        x(1)=ii-wi;
        break;
    end
end 

function [F,fwdrag,FS,dF,S_Fit] = find_drag(x,S,F,NP);
% find_drag - the program subtracts the drag force contribution
% time - the number of seconds it takes to go 1 micron distance. This value will depend on your pulling speed 
% fm- the mean value of force in the window
% x - is the vector that contains the index of the points corresponding to minimum and maximum force
% dF - is the drag force contribution. Note that multiplier 0.0069 is for a
% specific contilever. You will need to find this coefficient for the
% cantilever you will use.

% The cantilever retraction velocity (speed) relative to the surface will
% change depending on force exerted on the ecantilever. We thus defined a velocity at each point of
%the force separation data by differentiating separation with respect
%to time. One can estimate this derivative based on the actual
%separation data. However, sensor noise influences the estimate from the actual data.
%Instead we fit the data selected 5{th} order polynomial. We found that 
%this polynomial described the data reasonably well. Then we differentiated
%that polynomial to find relative velosity.
time=1;
time_step=1/NP;
t_vec=time_step:time_step:time;
delta=t_vec(2)-t_vec(1);
Tn=(t_vec(x(1):x(2)))';

nn=5;
warning off MATLAB:polyfit:RepeatedPointsOrRescale
[pp,s]=polyfit(Tn,S(x(1):x(2)),nn);
S_Fit=polyval(pp,Tn);
Sd=S_Fit-S(x(1):x(2));
DSep_Fit=gradient(S_Fit,delta);  %%%%%%%%%%%FOR THE VELOcITY OF 1 MICRON/SEC

dF=0.0069*DSep_Fit; %%%%%%for the TM in water 
FS=F(x(1):x(2))+dF/2;
fwdrag=F;
F=[F(1:x(1)-1);FS;F(x(2)+1:end)];

function err = fjcmin_weight(p1,x,f,T,ww,sigma_square,sigma_square_hat)
KT = 4.1*((T+273.15)/298); % pNnm (4.114E-21 J at 298 K J=Nm=1E21 pNnm 298 - 273.15 = 24.85 C = 76.73 F)
z = zeros(length(f),1);
z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
err = sum(ww.*(z-x).^2-sigma_square*log(ww));

function err = fjcmin_noweight(p1,x,f,T)
KT = 4.1*((T+273.15)/298); % pNnm (4.114E-21 J at 298 K J=Nm=1E21 pNnm 298 - 273.15 = 24.85 C = 76.73 F)
z = zeros(length(f),1);
z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
err = norm(z-x);

function z = fjc(f,p1,T)
KT = 4.1*((T+273.15)/298); 
z = zeros(length(f),1);
warning off MATLAB:divideByZero
z =p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));

function fname=find_file_name(k,filename);
if k==1;
    Short_fname=filename(end-6:end-4);
    fname = str2num(Short_fname);
elseif k==2
    Short_fname=filename(end-6:end-4);
    fname = str2num(Short_fname);
    Short_fname_num_k2=str2num(strcat(Short_fname,'.2'));
    fname = [fname; Short_fname_num_k2];
elseif k==3
    Short_fname=filename(end-6:end-4);
    fname = str2num(Short_fname);
    Short_fname_num_k3=str2num(strcat(Short_fname,'.3'));
    fname = [fname; Short_fname_num_k3];
end 

function [bootstat, bootsam] = bootstrp(nboot,bootfun,varargin)
%BOOTSTRP Bootstrap statistics.
%   Reference:
%      Efron, Bradley, & Tibshirani, Robert, J.
%      "An Introduction to the Bootstrap", 
%      Chapman and Hall, New York. 1993.

%   B.A. Jones 9-27-95, ZP You 8-13-98
%   Copyright 1993-2002 The MathWorks, Inc. 
%   $Revision: 2.11 $  $Date: 2002/01/17 21:30:04 $

% Initialize matrix to identify scalar arguments to bootfun.
la = length(varargin);
scalard = zeros(la,1); db = cell(la,1);

% find out the size information in varargin.
n = 1;
for k = 1:la
    [row,col] = size(varargin{k});
    if max(row,col) == 1
        scalard(k) = 1;
    end
    if row == 1 & col ~= 1
        row = col;
        varargin{k} = varargin{k}(:);
    end
    n = max(n,row);
end
% Create index matrix of bootstrap samples.
bootsam = unidrnd(n,n,nboot);
% Get result of bootfun on actual data and find its size. 
thetafit = feval(bootfun,varargin{:});
[ntheta ptheta] = size(thetafit);
% Initialize a matrix to contain the results of all the bootstrap calculations.
bootstat = zeros(nboot,ntheta*ptheta);
% Do bootfun - nboot times.
for bootiter = 1:nboot
    for k = 1:la
        if scalard(k) == 0
            db{k} = varargin{k}(bootsam(:,bootiter),:);
        else
            db{k} = varargin{k};
        end
    end
    tmp = feval(bootfun,db{:});
    bootstat(bootiter,:) = (tmp(:))';
end