% 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 

% Declaimer: The program has been for the elastin-like polypetides (ELPs).
% The use of the program for other molecules may require substntial
% modifications. User accepts to use this program at his or her own risk.


% Program requires the following inputs:
% a) force_matrix - is the matric that contains force valuess for the set of your curves 
% b) SepNew_matrix - is the matric that contains separation valuess for the set of your curves 
% c) fit_high_force - is the normalization force 
% d) Short_fname_matr - optional parameter that may specify the information
% about your files, for example the numbers of your files

% 
% dirname - the name of the directory where you store the force curves 
% curve_ind - indices of the force curves that were defined to be removed based on the maximum force
% rr_matr - curve indices that fall outside the average curve in the normalization filtering 
%will be filtered 
% force and SepNew - input force and Separation vectors 
% NormSep - normalized separation at the specified force 
% NormSep_low - is the valued of the normalized sepration that defines our low force region 
% force_matr_final - is the matrix that contains forces smoothed with lowess

% The program will output the 


function [NormSep,median_force,rr_matr,Short_fname_matr]= normalize_force(force_matrix,SepNew_matrix,fit_high_force,Short_fname_matr);
a_ind=1;
tt=1;
curve_ind=zeros(1,1);
NormSep_matr=zeros(size(force_matrix,1),1);
force_matr_final=zeros(NormSep_matr)
NormSep_low=0.75;

while tt<=size(force_matrix,2)
    force=force_matrix(:,tt);
    SepNew=SepNew_matrix(:,tt);  
    % Here we remove the curves with maximum force (above 100 nm) below our
    % normalization force. You may want to change both the separation and force value for your data 
    index=find(force>=fit_high_force & SepNew>=100);
    if length(index)==0;
        curve_ind=[curve_ind,tt];
        tt=tt+1;
        a_ind=a_ind+1;
        continue;
    end 
    %Here we perform nornmalization of the force curve
    index2=index(1);
    Sepone=SepNew(index2);
    NormSep=SepNew./Sepone; 
    NormSep_matr(:,tt)=NormSep;
        
    [NormSep_low_act,NormSep_low_ind]=min(abs(NormSep-NormSep_low));
    NormSep_selected=sort(NormSep(NormSep_low_ind:index2));
    force_selected=force(NormSep_low_ind:index2);
    Force_Smooth=smooth(NormSep_selected,force_selected,0.5,'lowess');
    force_matr_final=[force_matr_final, Force_Smooth];
    tt=tt+1;
end

median_force=median(force_matr_final,2);
st=1;
while st<=size(force_mean_matr_final,1);
    nn=1:size(force_mean_matr_final,2);
    var_force(st,:)=std(median_force(st)-force_matr_final(st,nn));
    st=st+1;
    nn=1;
end
% average_curve defines the interval in which the normalized curve should fall
average_curve=[median_force-2*var_force, median_force+2*var_force];

%the lines below will find either to except or regect a given curve based
%on the deviation from allowed deviation from the median curve
fg=1;
rr_matr=zeros(1,1);
while fg<=size(force_mean_matr_final,2);
    up_force_check=force_matr_final(:,fg)-average_curve(:,1);  %THE LOW BOUND ON FORCE 
    down_force_check=force_matr_final(:,fg)-average_curve(:,2); %THE HIGH BOUND ON FORCE 
    if min(up_force_check)<0 | max(down_force_check)>0;
        rr_matr=[rr_matr,fg];
    end
    fg=fg+1;
end

%rr_matr is the curve indices that fall outside the average curve and thus
%will be filtered 
rr_matr(:,1)=[];
curve_ind(:,1)=[];
Short_fname_matr(curve_ind,:)=[]
Short_fname_matr(rr_matr,:)=[];
force_mean_matr_final(:,rr_matr)=[];
NormSep_matr(:,rr_matr)=[];
force_mean_matr_filt=force_mean_matr_final;
NormSep_matr_filt=NormSep_matr;