function Fit_FJC_fix_kuhn_stiff

% Iterative minimization of the Extended FJC Chain model. Multiple Ccurve Model
% By Alexei Valiaev  October 2006
% Please send comments or bug reports to Alexei Valiaev at av12 at duke dot edu
% Program is not yetthoroughly tested, so please any bug reports 

% Instead of the stagewise minimization now the program performs iterative
% minimization. The Extended FJC chain model (i.e three parameters is used.
% First we fix the kuhn length (kuhn) and segmental stiffness (stiff) and
% minimize over the contour length to get the min sum of squared residuals.
% Then we fix the contour length and minimize over the kuhn length and
% segmental stiffness. We do this kind of minimization several times until
% we get convergences.  The convergence is defined by comparing the
% differences in the Kuhn length and segmental stiffness at the adjacent iterations
%From my tests about 5 iterations is usually sufficient. Same type of iterative minimization can be used for FJC model in that case the process will be even simplier


%Variables
% dirname - the name of the directory where you store the force curves 

%Functions

%fit_data_ext_FJC_kuhn_stiff_weight - here we do the minimization with respect to the contour length, so we fix kuhn and segmental stiffness
%This function includes subfunctions: fjcextmin_fix_stiff_kuhn and
%fjcextmin_fix_stiff_kuhn_weight. Remember that we also need to use the interatively reqighted minimization over the weights.  So it is kind of two iterative minimization processes happening in this  program
% fjcextmin_fix_contour - here we do minimization over the kuhn length and segmental stiffness so we dix the contour length 
% calc_sep_residual_FJC_new_appr_noplot - here we calculate the separation residuals

clear all
close all

kuhn_init=0.42;
kuhn=kuhn_init;
delta_kuhn=0.002; 
sep_force_matrix=zeros(4096,1);
%%Ok Here I load all the files%%%%%%%

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

load points150303_test_program
points=points150303;

sep_force_matrix=zeros(4096,1);
ii=1;

% load all of the force curves and store them in one big matrix 
while ii<=numfiles
    close all
    filename = files(ii+2).name;
    eval(['cd ',dirname]); 
    load(filename);
    eval(['cd ',curdir]);
    endmat=length(filename)-4; 
    Short_fname=filename(end-6:end-4);
    eval(['SepNew=SepNew',filename(7:endmat),';']);
    eval(['force=force',filename(7:endmat),';']);
    [curve_max_force(ii),ind_max_force(ii)]=max(-force);
    L_max_force(ii)=SepNew(ind_max_force(ii));  
    sep_force=[SepNew,force];
    sep_force_matrix=[sep_force_matrix, sep_force];
    ii=ii+1;
end
sep_force_matrix(:,1)=[];

T=22;
p1_init=[L_max_force]; 
contour_matrix=zeros(length(L_max_force),1);
kuhn_stiff_matrix=zeros(1,2);
kuhn_stiff=[0.3 4000];
kk=1;
while kk<=10
    par=fit_data_ext_FJC_kuhn_stiff_weight(sep_force_matrix,p1_init,kuhn_stiff,points);
    contour=par;
    clear par
    if kk==1
        p1_init=[0.3 4000];
    else 
        p1_init=kuhn_stiff;
    end
    [par,fval,exitflag,output] = fminsearch(@fjcextmin_fix_contour,p1_init,optimset('Tolx',1e-7,'MaxFunEvals',3000,'MaxIter',7000),T,sep_force_matrix,points,contour);
    kuhn_stiff=par;
    kuhn_stiff_matrix=[kuhn_stiff_matrix;kuhn_stiff];
    ii=1;
    rr=1;
    % Here we calculate the sum of squared residuals for all of the curves
    % in the data set 
    sep_residual_matrix=zeros(length(force),1);
    while ii < size(sep_force_matrix,2)
        SepNew=sep_force_matrix(:,ii);
        force=sep_force_matrix(:,ii+1);
        b=[points(rr,1) points(rr,2)]; 
        Force_new = -force(b(1):b(2));
        Sep_selected = SepNew(b(1):b(2));
        z = zeros(length(Sep_selected),1);
        ff1_mod=linspace(0,Force_new(end),length(Sep_selected))';
        KT = 4.1*((T+273.15)/298);
        z = contour(rr)*(coth(ff1_mod.*kuhn_stiff(1)/KT)-KT./(ff1_mod*kuhn_stiff(1))) + (contour(rr)*ff1_mod/(kuhn_stiff(1)*kuhn_stiff(2))); 
        sep_residual=calc_sep_residual_FJC_new_appr_noplot(contour,kuhn_stiff,Force_new,Sep_selected,rr);
        %store the martix of residuals at every iteration
        sep_residual_matrix(:,rr)=[zeros(points(rr,1)-1,1);sep_residual; zeros((4096-points(rr,2)),1)];
        ii=ii+2;
        rr=rr+1;
    end
    p1_init=contour;
    contour_matrix=[contour_matrix,contour'];
    if kk~=1;
        delta_kuhn_stiff=[(kuhn_stiff_matrix(kk,1)-kuhn_stiff_matrix(kk-1,1)),(kuhn_stiff_matrix(kk,2)-kuhn_stiff_matrix(kk-1,2))];
    else 
        delta_kuhn_stiff=[1 1000]; % give some intial values
    end
    %Here we put the condition so that if the differnce between the kuhn ans
    %segmental stiffness is small at the two consequtive iteration then we
    %finish the interative procedure
    if delta_kuhn_stiff(1)<0.01 & delta_kuhn_stiff(2)<200
        break;
    end 
    kk=(kk+1);
end
kuhn_stiff_matrix(1,:)=[];
contour_matrix(:,1)=[];
kuhn_stiff_matrix_original=kuhn_stiff_matrix;
contour_matrix_original=contour_matrix;
disp('Finished the iterative Procedure')

% Saving the data
convertdir = 'par_fix_Kuhn_stiff'; 
mkdir(dirname,convertdir)
savename=[dirname,'\',convertdir,'\','pFJC_fix_Kuhn_stiff',filename(7:(end-8))]; 
kuhn_stiff_matrix_name=['kuhn_stiff_matrix',filename(9:(endmat-4))];
contour_matrix_name=['contour_matrix',filename(9:(endmat-4))];
sep_residual_matrix_name=['sep_residual_matrix',filename(9:(endmat-4))];

eval([kuhn_stiff_matrix_name,' = kuhn_stiff_matrix;']);
eval([contour_matrix_name,' = contour_matrix;']);
eval([sep_residual_matrix_name,' = sep_residual_matrix;']);
eval(['save(savename,''',kuhn_stiff_matrix_name,''',''',contour_matrix_name,''',''',sep_residual_matrix_name,''')']);

%Here we optimize with respect to contour length
function par=fit_data_ext_FJC_kuhn_stiff_weight(sep_force_matrix,p1_init,kuhn_stiff,points);
ss=1;
delta_p1_bootstr=1;
T=22; 
while ss<=50 
    if  abs(delta_p1_bootstr)<=0.01
        break;
    else 
        if ss==1;
            [par,fval,exitflag,output] = fminsearch(@fjcextmin_fix_stiff_kuhn,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),T,sep_force_matrix,points,kuhn_stiff);
        else
            [par,fval,exitflag,output] = fminsearch(@fjcextmin_fix_stiff_kuhn_weight,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),T,sep_force_matrix,points,kuhn_stiff,sep_resid_weight);
            delta_p1_bootstr=p1_bootstr_matrix(ss)-p1_bootstr_matrix(ss-1);    
        end
        %%%Ok now I am dealing with all curves at the same time so all the
        %%%so I will need to remake the weighting
        [sep_residual,sep_resid_weight_matr,sigma_square,ww,sep_resid_weight]=calc_sep_residual_weight(par,kuhn_stiff,sep_force_matrix,points);
        warning off MATLAB:divideByZero         
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%                
    end 
    ss=ss+1;
end


function err = fjcextmin_fix_stiff_kuhn(p1_init,T,sep_force_matrix,points,kuhn_stiff)
p1=p1_init;
kuhn=kuhn_stiff(1);
stiff=kuhn_stiff(2);
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)
ii=1;
rr=1;
kuhn=kuhn_stiff(1);
stiff=kuhn_stiff(2);
while ii < size(sep_force_matrix,2)
    SepNew=sep_force_matrix(:,ii);
    force=sep_force_matrix(:,ii+1);
    b=[points(rr,1) points(rr,2)]; 
    Force_new = -force(b(1):b(2));
    Sep_selected = SepNew(b(1):b(2));
    z = zeros(length(Sep_selected),1);
    %%%Here we fixed the kuhn length and the segmental stiffness if in order to load the contour length approprietely  
    z = p1(rr)*(coth(Force_new.*kuhn/KT)-KT./(Force_new*kuhn)) + (p1(rr)*Force_new/(kuhn*stiff)); 
    ff1_mod=linspace(0,Force_new(end),length(Sep_selected))';
    warning off MATLAB:divideByZero
    ZZ1_kernel_mod = fjcext_kuhn_stiff(ff1_mod,p1,T,kuhn,stiff,rr);
    err_curve(:,rr) = norm(z-Sep_selected);   
    ii=ii+2;    
    rr=rr+1;
    clear zz
end 
err=sum(err_curve);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function err = fjcextmin_fix_stiff_kuhn_weight(p1_init,T,sep_force_matrix,points,kuhn_stiff,sep_resid_weight)
p1=p1_init;
kuhn=kuhn_stiff(1);
stiff=kuhn_stiff(2);
KT = 4.1*((T+273.15)/298); 
ii=1;
rr=1;
kuhn=kuhn_stiff(1);
stiff=kuhn_stiff(2);
while ii < size(sep_force_matrix,2)
    SepNew=sep_force_matrix(:,ii);
    force=sep_force_matrix(:,ii+1);
    b=[points(rr,1) points(rr,2)]; 
    Force_new = -force(b(1):b(2));
    Sep_selected = SepNew(b(1):b(2));
    z = zeros(length(Sep_selected),1);
    
    z = p1(rr)*(coth(Force_new.*kuhn/KT)-KT./(Force_new*kuhn)) + (p1(rr)*Force_new/(kuhn*stiff)); 
    ff1_mod=linspace(0,Force_new(end),length(Sep_selected))';
    warning off MATLAB:divideByZero
    ZZ1_kernel_mod = fjcext_kuhn_stiff(ff1_mod,p1,T,kuhn,stiff,rr);
    err_curve(:,rr) = norm(sep_resid_weight.*z-sep_resid_weight.*Sep_selected);   
    ii=ii+2;    
    rr=rr+1;
    clear zz
end 
err=sum(err_curve);

%%%%%%%%%%Here we optimize with respect to kuhn and stiffness
function err = fjcextmin_fix_contour(p1,T,sep_force_matrix,points,contour)
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)
ii=1;
rr=1;
while ii < size(sep_force_matrix,2)
    SepNew=sep_force_matrix(:,ii);
    force=sep_force_matrix(:,ii+1);
    b=[points(rr,1) points(rr,2)]; 
    Force_new = -force(b(1):b(2));
    Sep_selected = SepNew(b(1):b(2));
    z = zeros(length(Sep_selected),1);
    %%%Here we fixed the Contour lengths
    z = contour(rr)*(coth(Force_new.*p1(1)/KT)-KT./(Force_new*p1(1))) + (contour(rr)*Force_new/(p1(1)*p1(2)));
    ff1_mod=linspace(0,Force_new(end),length(Sep_selected))';
    warning off MATLAB:divideByZero
    ZZ1_kernel_mod = fjcext_contour(ff1_mod,p1,T,contour,rr);
    err_curve(:,rr) = norm(z-Sep_selected);   
    ii=ii+2;    
    rr=rr+1;
    clear zz
end  
err=sum(err_curve);

function err = fjcextmin_fix_contour_weight(p1,T,sep_force_matrix,points,contour)
%p1_matrix =p1;
%weighted_force = f./weight_factor_kernel;
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(4096,1);
ii=1;
rr=1;
%p1(3:length(p1))=L_contour;
while ii < size(sep_force_matrix,2)
    SepNew=sep_force_matrix(:,ii);
    force=sep_force_matrix(:,ii+1);
    b=[points(rr,1) points(rr,2)]; 
    Force_new = -force(b(1):b(2));
    Sep_selected = SepNew(b(1):b(2));
    z = zeros(length(Sep_selected),1);
    %%%Here we fixed the Contour lengths
    z = contour(rr)*(coth(Force_new.*p1(1)/KT)-KT./(Force_new*p1(1))) + (contour(rr)*Force_new/(p1(1)*p1(2)));
    ff1_mod=linspace(0,Force_new(end),length(Sep_selected))';
    warning off MATLAB:divideByZero
    ZZ1_kernel_mod = fjcext_contour(ff1_mod,p1,T,contour,rr);
    err_curve(:,rr) = norm(sep_resid_weight.*z-sep_resid_weight.*Sep_selected);   
    ii=ii+2;    
    rr=rr+1;
    clear zz
end  
err=sum(err_curve);

function z = fjcext_contour(f,p1,T,contour,rr)
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 = contour(rr)*(coth(f.*p1(1)/KT)-KT./(f*p1(1))) + (contour(rr)*f/(p1(1)*p1(2)));
warning off MATLAB:divideByZero
%z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1))) + (p1(2)*f/(p1(1)*p1(3)));

function z = fjcext_kuhn_stiff(f,p1,T,kuhn,stiff,rr)
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(rr)*(coth(f.*kuhn/KT)-KT./(f*kuhn)) + (p1(rr)*f/(kuhn*stiff)); 
warning off MATLAB:divideByZero


function z = fjcext_kuhn_stiff_weight(f,p1,T,kuhn,stiff,rr)
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(rr)*(coth(f.*kuhn/KT)-KT./(f*kuhn)) + (p1(rr)*f/(kuhn*stiff)); 
warning off MATLAB:divideByZero

%OK I FINISHED AT THE END OF THIS FILE  
function  [sep_residual,sep_resid_weight_matr,sigma_square,ww,sep_resid_weight]=calc_sep_residual_weight(contour,kuhn_stiff,sep_force_matrix,points);
%function sep_residual=calc_sep_residual_FJC_new_appr_noplot(contour,kuhn_stiff,Force_new,Sep_selected,rr);
ii=1;
rr=1;
sep_resid_weight_matr=zeros(size(sep_force_matrix,1),1);
while ii < size(sep_force_matrix,2)
    SepNew=sep_force_matrix(:,ii);
    force=sep_force_matrix(:,ii+1);
    b=[points(rr,1) points(rr,2)]; 
    Force_new = -force(b(1):b(2));
    Sep_selected = SepNew(b(1):b(2));
    %Force_new=-sep_force_matrix(:,rr+1);
    T=22;
    KT = 4.1*((T+273.15)/298);
    num_point=1;
    sep_residual=zeros(1,1);
    Index_matrix=zeros(1,1);
    while num_point < length(Force_new);
        f=Force_new(num_point);
        z = contour(rr)*(coth(f.*kuhn_stiff(1)/KT)-KT./(f*kuhn_stiff(1))) + (contour(rr)*f/(kuhn_stiff(1)*kuhn_stiff(2)));
        sep_residual_element=Sep_selected(num_point)-z; 
        sep_residual=[sep_residual;sep_residual_element]; 
        num_point=num_point+1;
    end
    sigma_square=sum((1/(length(Force_new)-2))*sep_residual.^2);
    L_B=(coth(Force_new.*kuhn_stiff(1)./KT)-KT./(Force_new.*kuhn_stiff(1)));
    aa=KT./(Force_new.*kuhn_stiff(1));
    coef=2.*aa.*L_B;
    ww=(contour(rr).*kuhn_stiff(1).*(1-L_B.^2-coef));
    sep_resid_weight=(1./(sigma_square*ww));
    sep_resid_weight_matr(:,rr)=[sep_resid_weight;zeros(size(sep_resid_weight_matr,1)-length(sep_resid_weight),1)];
    ii=ii+2;    
    rr=rr+1;
end  