% Multiple Curve Model  Fitting 
% By Alexei Valiaev  October 2006
% Please send comments or bug reports to Alexei Valiaev at av12 at duke dot edu
% Program requires loading of the matrix with point indices that would define the
% region in which you fit the data. The output of tbe Single Curve Model
% provides you the matrix of point indices

% dirname - the name of the directory where you store the force curves 
% kuhn_i - the initiatial values for the Kuhn length parameter of the FJC model (Value of kuhn_i does not change 
% nd- specified the number of steps in the stagewise minimization process. Increase in "nd"
%increase the simulation time. 
%execution time 
% kuhn - the Kuhn value that we increment during our stagewise minimization.
% dk - is the Kuhn length increment that we use in the minimization. The
% incremenet can be large or smaller than given depending on accuracy you
% want and computer power. Smaller increment will give better accuracy
% but takes it takes longer to simulate 
% err_MA - the matrix that contains squared separation residuals
% T - is the value of the room temperature
% NP - number of points in your force curve 
% N_iter - number of iterations used in the stagewise minimization 
% iterN - is the current iteration number 
%  ii -  index of the forces in the SF_M matrix
% 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
% SF_M- is a large matric that contains all your data in the folder "dirname"
% L_i - initial value for the contour length. 
% L_kuhn - is a 3D matrix it keeps values from all the steps of the  minimization procedure. Referring to this matric maybe helpful for
% debugging if something goes wrong with your data 
% f_fit - linearly spaced vector for the force 
% Z_fit - FJC model fit data
% sep_r - residuals in separation (distance)
% perm_vec - the vector of odd numbers start from 1 and goes to numfiles*2
% kuhn_min - Value of the Kuhn length at which the sepraration residuals are minimized at each value of the bootsrap sample
% rand_perms - matrix that contains the indices for the force curves obtained using bootstrap sampling 
% sigma_square and sigma_square_hat - scale parameters to be estimated 
% Vd_f - is the variance function 


% FUNCTIONS

% find_drag - the program subtracts the drag force contribution
%sep_r_noweight - calculates residuals in separation at the first iteration of the simulation
% fit_FJC_weight - this finction fits the FJC model to the force curve data
% for each value of the Kuhn length. Then the separation residuals are
% calcualated from which we estimate sum of squared residuals over all curves in the data set. Then the value of the  Kuhn
%length corresponding to the smallest sum of squared residuals
% bootstrp - performs the boostrapping of the data. See paper references for details 
% fjcmin_fix_kuhn_no_weight - - performs FJC parameter estimation. No weights are used in the minimization
% fjcmin_fix_kuhn_no_weight - - performs FJC parameter estimation. No weights are used in the minimization

%Additional commnents 
% Try running the simulation witb a course step first. i.e make dk larger
% this will help to zoom in at the right kuhn length range in the next run.

% The program saves:
%a) The values of the Kuhn lengths at which the sum squared of separation residuals are minimized 
%b) Uncertainty in the Kuhn length obtained from the bootstrap samples 
%c) mean value of the kuhn length calculated from the bootstrap samples

function Fit_FJC_fic_kuhn
clear all
close all
dirname='D:\LAPTOP\DISK_C\DIWORK\2003\15March_2003\test'; % put the folder name with your force curves

curdir = pwd;
files = dir(dirname);
numfiles = length(files) - 2;

load points150303_test_program
points=points150303;

T=25;
indL=7;
kuhn_i=0.35;
kuhn=kuhn_i;
dk=0.002; 
NP=4096;

ii=1;
nd=5;
N_iter=2;
iterN=1;

L_matr=zeros(1,1);
SF_M=zeros(NP,1);

while ii<=numfiles
    filename = files(ii+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),';']);
        % Here we put the force and separation from the your force curves into one big matrix
    S=sort(S);
    SF=[S,F];
    SF_M=[SF_M, SF];
    ii=ii+1;
end
SF_M(:,1)=[];

while  iterN <= N_iter; 
    ii=1;
    rr=1;
    %  ii is the index of the force curve
    while ii < size(SF_M,2)
        close all
        S=SF_M(:,ii);
        F=SF_M(:,ii+1);
        b=[points(rr,1) points(rr,2)]; 
        [F,fwdrag,FS,dF,S_Fit] = find_drag(b,S,F,NP);
        FS = F(b(1):b(2));
        SS = S(b(1):b(2));
        SS= sort(SS);
        L_i =200; 
        if ii==1;
            sigma_square_matr=zeros(1,1);
            err_M=zeros(1,1);
        end
        [L,L_matr,sigma_square_matr,err_M]=fit_FJC_weight(SS,FS,ii,kuhn,L_matr,sigma_square_matr,err_M);
        L_kuhn_matr=[kuhn,L];      
        f_fit = linspace(0,-FS(end),length(SS))';
        Z_fit = fjc_fix_kuhn(f_fit,L,kuhn);
        warning off MATLAB:divideByZero
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%
        sep_r=sep_r_noweight(L_kuhn_matr,FS,SS);        
        if rr==1 & iterN==1
            par=zeros(1,2,1);  
        else   
            par(rr,:,iterN)=L_kuhn_matr;
        end 
        ii=ii+2;    
        rr=rr+1;
    end  
    if iterN==1;
        err_MA=zeros(1,1);
    end
    %Here is the error function which is being minimized
    err_MA=[err_MA;sum(sigma_square_matr)];
    %Sum of squared residuals are saved at each iteration (coolumn) and at each boostrap (row)
    iterN=iterN+1;
    kuhn=kuhn+dk;
    kuhn_matrix(:,iterN)= kuhn;
    disp('OK')
end  
err_MA(1,:)=[];
[min_err,ind_err]=min(err_MA);

%BELOW CACULATE THE KUHN LENGTH THAT CORRESSPONDS TO THE MINIMUM SUM OF SQUARED RESIDUALS
[min_err,ind_err]=min(err_MA);

%BELOW ARE THE BOOTSTRAPS CALCULATIONS BOOTSTRPING AND ESTIMATING
%CONFIDENCE INTERVALS% % Now we bootstrap the curves 
perm_vec=[1:2:numfiles*2]';
%%Bootsam containes the indexes randomly drawn nd times from the perm_vec
[bootstat,bootsam] = bootstrp(nd,'mean',perm_vec); 
rand_perms=perm_vec(bootsam); %

p1_kernel_matrix=zeros(N_iter,3);
L_matr=zeros(1,1);
L_kuhn=zeros(numfiles,N_iter,nd);
kuhn_matrix=zeros(1,1,1);
hh=1;
while hh<=nd
    ss=1;
    SF_perm=zeros(length(F),1);
    points_perm=zeros(1,2);
    while ss<=numfiles
        SF_perm=[SF_perm, SF_M(:,rand_perms(ss,hh)),SF_M(:,rand_perms(ss,hh)+1)];
        points_perm=[points_perm; points(0.5*(rand_perms(ss,hh)+1),1) points(0.5*(rand_perms(ss,hh)+1),2)];
        ss=ss+1;
    end 
    SF_perm(:,1)=[];
    points_perm(1,:)=[];
    iterN=1;
    kuhn=kuhn_i;
    while  iterN <= N_iter; 
        ii=1;
        rr=1;
        while ii <= size(SF_perm,2)
            S=SF_perm(:,ii);
            F=SF_perm(:,ii+1);
            b=[points_perm(rr,1) points_perm(rr,2)]; 
            SS = S(b(1):b(2));
            FS=F(b(1):b(2));
              if ii==1;
                sigma_square_matr=zeros(1,1);
            end
            [L,L_matr,sigma_square_matr,err_M]=fit_FJC_weight(SS,FS,ii,kuhn,L_matr,sigma_square_matr,err_M);
            L_kuhn_matr=[kuhn,L];      
            f_fit=linspace(0,-FS(end),length(SS))';
            Z_fit = fjc_fix_kuhn(f_fit,L,kuhn);
            warning off MATLAB:divideByZero
            sep_r=sep_r_noweight(L_kuhn_matr,FS,SS);        
            if rr==1 
                par=zeros(1,2);  
                par=[par; L_kuhn_matr];
            else   
                par=[par; L_kuhn_matr];  
            end 
            L_kuhn(rr,iterN,hh)= par(rr+1,2);  
            ii=ii+2;    
            rr=rr+1;
        end  
        err_bootstrap(iterN,hh)=sum(sigma_square_matr);
        iterN=iterN+1;
        kuhn=kuhn+dk;
        kuhn_matrix(:,iterN,hh)= kuhn;
        disp('OK')
    end  
    hh=hh+1; 
    iter=num2str(hh); disp(iter)
end 
CL=clock; digits(4); Hour=vpa(CL(4)); Minute=vpa(CL(5));

% Calculate bootstrap uncertinity%
[min_chi2_sep,ind_min_chi2_sep]=min(err_bootstrap);
%We find the value of the kuhn lenght that minimize the squared residuals at each iteration%%%
for ff=1:length(ind_min_chi2_sep)
    kuhn_min(ff)=kuhn_i+dk*(ind_min_chi2_sep(ff)-1);
end
mean_kuhn=sum(kuhn_min)/length(kuhn_min);
for pp=1:nd
    BB_kuhn(pp)=(kuhn_min(1,pp)-mean_kuhn)^2;
end 
ukb=((1/(nd-1))*sum(BB_kuhn))^0.5;

%%% Saving the data %%%%%%%%%
convertdir = 'par_fix_Kuhn'; 
mkdir(dirname,convertdir)
savename=[dirname,'\',convertdir,'\','pFJC',filename(7:(end-8))]; 
Kuhn_Min=['kuhn_min',filename(9:(endmat-4))];
Uncert_Kuhn=['uncert_kuhn',filename(9:(endmat-4))];
Mean_Kuhn=['mean_kuhn',filename(9:(endmat-4))];
Points=['points',filename(9:(endmat-4))];
Delta_Kuhn=['dk',filename(9:(endmat-4))];
eval([Kuhn_Min,' = kuhn_min;']);
eval([Mean_Kuhn,' = mean_kuhn;']);
eval([Uncert_Kuhn,' = ukb;']);
eval([Points,' = points;']);
eval([Delta_Kuhn,' = dk;']);
eval(['save(savename,''',Kuhn_Min,''',''',Uncert_Kuhn,''',''',Mean_Kuhn,''',''',Points,''',''',Delta_Kuhn,''')']);

%%%%%%%%%%Ploting the distribution of Kuhn segment lengths %%%%%%%%%%%
kuhn_sort=sort(kuhn_min);
interval_start=min(kuhn_min)-0.1;
interval_end=max(kuhn_min)+0.1;

interval=interval_start:dk:interval_end;
kuhn_group=histc(kuhn_sort,interval);
figure(1)
bar(interval,kuhn_group,'histc'),
title('Kuhn Length','FontSize',12,'Color',[0.5 0.1 0.1])
xlabel('Kuhn Length [nm]','FontSize',20,'Color',[1 0.1 0.3])
ylabel('Ratio of occurrences','FontSize',20,'Color',[1 0.1 0.3])
set(gca,'FontSize',24)

function [L,L_matr,sigma_square_matr,err_M]=fit_FJC_weight(SS,FS,ii,kuhn,L_matr,sigma_square_matr,err_M);
L_i=200;
ss=1;
delta_kuhn=1;
T=22; 
while ss<=50 
    if  abs(delta_kuhn)<=0.01
        break;
    else 
        if ss==1;
            [L,fval,exitflag,output] = fminsearch(@fjcmin_fix_kuhn_no_weight,L_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,kuhn);
        else
            [L,fval,exitflag,output] = fminsearch(@fjcmin_fix_kuhn_weight,L_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,sep_resid_weight,kuhn,sep_resid_weight,sigma_square);
            delta_kuhn=L_matr(ss)-L_matr(ss-1);    
        end
        L_matr=[L_matr;L];
        f_fit=linspace(0,-FS(end),length(SS))';
        Z_fit = fjc_fix_kuhn(f_fit,L,kuhn);
        warning off MATLAB:divideByZero         
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%                
        [sep_r,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f,err_value]=calc_sep_residual_weight(L,FS,SS,kuhn);
        %Lets add in some control parameters 
    end 
    ss=ss+1;
end
err_M=[err_M;err_value];
sigma_square_matr=[sigma_square_matr;sigma_square];


function [sep_r,sep_resid_weight,sigma_square_mod,sigma_square_hat,Vd_f,err_value]=calc_sep_residual_weight(L,FS,SS,kuhn);
T=25;
KT = 4.1*((T+273)/298);
f=-FS;
z = L*(coth(f.*kuhn/KT)-KT./(f*kuhn));
sep_r=SS-z; 

b=kuhn*(-FS)./KT;
%  # Langevin(beta)
Lng_b=(coth(b)-1./b);
Vd_f=L*kuhn*(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);
sep_resid_weight=(1./Vd_f);
sigma_square_mod=sum((sep_r.^2.*sep_resid_weight)-sigma_square_hat*log(sep_resid_weight));
err_value=(1/sigma_square)*sum(sep_resid_weight.*sep_r.^2-sigma_square*(log(sep_resid_weight)-log(sigma_square)));

function err = fjcmin_fix_kuhn_weight(L,x,f,sep_resid_weight,kuhn,sep_r,sigma_square)
T=25;
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 = L*(coth(f.*kuhn/KT)-KT./(f*kuhn));
err = (1/sigma_square)*sum(sep_resid_weight.*(z-x).^2-sigma_square*(log(sep_resid_weight)-log(sigma_square)));

function err = fjcmin_fix_kuhn_no_weight(L,x,f,kuhn)
T=25;
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 = L*(coth(f.*kuhn/KT)-KT./(f*kuhn));
err = norm(z-x);

function z = fjc_fix_kuhn(f,L,kuhn);
T=25;
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 = L*(coth(f.*kuhn/KT)-KT./(f*kuhn));

function sep_r=sep_r_noweight(L_kuhn_matr,FS,SS);
T=25;
KT = 4.1*((T+273.15)/298);
num_point=1;
sep_r=zeros(1,1);
Index_matrix=zeros(1,1);
L=L_kuhn_matr(2);
kuhn=L_kuhn_matr(1);
while num_point < length(FS);
    f=-FS(num_point);
    z = L*(coth(f.*kuhn/KT)-KT./(f*kuhn));
    sep_ri=SS(num_point)-z; 
    sep_r=[sep_r;sep_ri]; 
    num_point=num_point+1;
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 [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