%catch bond multistart -- use the multistart algorithm to fit catch bond
%data
    
file_path = './data_to_fit_to_two_state_catch.xlsx';
% file_path = '/Users/amywang/Desktop/analysis/trap-analysis/first_step.xlsx'
datatable = readtable(file_path, 'Sheet', 'last_step multistep ABD'); % read sheet 1, ABD data, single step
data = datatable{:,1:2};
% zeroconstraint = 100;           %CHANGE zero force constraint 0.0754 for tern, 0.0586 for vh, 100 if no using
                                %forces in column 1, corresponding dwell times in column 2
num_fits = 100;                  %CHANGE number of of times to fit the data. Most data need more than one try to find a global minimum

k10 = 11;               %parameters from Buckley et al. 
k20 = 0.14;
k12 = 3;
k21 = 20;
x10 = 0.1;
x20 = 0.4;
x12 = 0.2;
x21 = 4;

initial = [k10 k20 k12 k21 x10 x20 x12 abs(x21)];   %important, pass abs value 
% of x21 so that can apply positivity constraint. Also, please note variable order.


%%%%fit the data%%%%

%lots of possible options to set on the optimization --read the
%documentation
%options = optimoptions('fmincon','Display','iter')

%A b Aeq and beq need to be specified as empty or fmincon will be mad
A = [];
b = [];
Aeq = [];
beq = [];
f = @(x)catchbond(x, abs(data));  %anonymous function so that optimization occurs
%on x, but the catchbond function (see below) is passed  x, data, and norm. Here, x
%corresponds to the parameter vector

%Set upper and lower bounds on fit parameters
%For the lower bound, nonnegativity seems to be important in getting good 
%fits. One can argue that it is physically justified. However, when a
%parameter is very close to zero it probably means the fit want so to make
%it negative, and this needs to be reported as fixed at a bound. The same
%applies to fit parameters that hit the upper bound. Pegging out at the
%upper bound also strongly suggests that the fit is not physically
%reasonable. GA doesn't occasionally like zero as a lower bound, so better
%to make it small but finite.

% general lower bounds 
lb = initial*0.02;
ub = initial*50; 

%commented out line does a local optimization using fmincon
% [output, fval] = fmincon(f, initial, A, b, Aeq, beq, lb, ub, f2)

garesults = [];
gafvals = [];
opts = optimoptions('ga', 'Display', 'off', 'ConstraintTolerance', 1e-6, 'FunctionTolerance', 1e-9, 'UseVectorized', true); %'UseParallel', true);
tic
for i = 1:num_fits
%     [results,fval,exitflag,output] = ga(f, 8, A, b, Aeq, beq, lb, ub, f2, opts);
    [results,fval,exitflag,output] = ga(f, 8, A, b, Aeq, beq, lb, ub, [], opts);
    garesults = [garesults; results];
    gafvals = [gafvals; fval];
    disp(i)
end
toc

sorted_results = [garesults gafvals];
sorted_results = sortrows(sorted_results, 9);

function probs = catchbond(x, data)
% perform a maximum likilihood calculate for a 2-state, single-direction
% catch bond. Follows the nomenclature of the Thomas/Vogel 2006 BPJ paper

k10 = x(:,1);
k20 = x(:,2);
k12 = x(:,3);
k21 = x(:,4);
x10 = x(:,5);
x20 = x(:,6);
x12 = x(:,7);
x21 = -x(:,8);

beta = 4.1;
F= abs(data(:,1));
dwells = data(:,2);

%these rates have n elements, where n is the number of individual force points

r10 = k10.*exp(F'.*x10/beta);
r20 = k20.*exp(F'.*x20/beta);
r12 = k12.*exp(F'.*x12/beta);
r21 = k21.*exp(F'.*x21/beta);


B1 = 0.001;
B2 = 1-B1;

%this module assumes flux balance
% B1 = k21.*k10./(k21.*k10+k12.*k20);
% B2 = 1-B1;
 
%this module assume equilibration on F-actin
% B1 = r21./(r21+r12);       %fraction weak at t = 0
% B2 = 1-B1;

%this module assumes bound in strong state
% B2 = 0;
% B1 = 1;

b = r10 + r20 + r12 + r21;
c = r21.*r10 + r10.*r20 + r12.*r20;

lambda1 = (b + (b.^2 - 4*c).^0.5)./2;
lambda2 = (b - (b.^2 - 4*c).^0.5)./2;
    
C1 = (r21+r12+B1.*r20+B2.*r10-lambda1)./(lambda2-lambda1);
C2 = (r21+r12+B1.*r20+B2.*r10-lambda2)./(lambda1-lambda2);

%probs is the non-normalized probability density function.
probs = lambda1.*C1.*exp(-lambda1.*dwells') + lambda2.*C2.*exp(-lambda2.*dwells');
probs = probs + 1e-100; %deal with Matlab barfing when it takes log of a negative number

% probs = B1.*(r10+r12).*exp(-(r10+r12).*dwells/beta) + B2.*(r20+r21).*exp(-(r20+r21).*dwells/beta);

probs = -sum(log(probs'));

end

