% fitfrap
% For Joao Firmino @ Knust lab
%
% Jean-Yves Tinevez - January 2009


%% Parameters

% Draw figures?
show_figure = true;

% Fit by single expoenential?
fit_single_exp = true;
% Starting value for tau
tau_start = 30;

% Fit by two exponentials?
fit_two_exp = true;
% Starting values for taus
tau1_start = 3;
tau2_start = 300;

% Tolerance over beginning of recoevery, for fitting
tolerance = 0.1;

% Column names
time_column_name = 'time';
bleach_column_name = 'bleach';
nonbleach_column_name = 'nonbleach';

% File location
data_dir    = '/Users/tinevez/Projects/Data/J. Firmino - FRAP analysis';
xls_file    = '13Jan_embryo1_anterior_withalignstacks.xls';
file_path   = fullfile(data_dir,xls_file);


%% Prepare fit model
exp1_type = fittype( ' A * (1-exp( - (t-t0)/tau ) )' ,...
    'coefficients' , {'tau', 'A','t0' }, ...
    'independent', 't' );

exp2_type = fittype( ' A1 * (1-exp( - (t-t0)/tau1 ) ) + A2 * (1-exp( - (t-t0)/tau2 ) )' ,...
    'coefficients' , {'tau1', 'tau2' ,'A1', 'A2', 't0' }, ...
    'independent', 't' );


%% Operate

% Import xls
warning('off','MATLAB:xlsread:ActiveX')
[ data headers ] = xlsread(file_path);

% Collect relevant data
time_cols       = strmatch(time_column_name,headers,'exact');
bleach_cols     = strmatch(bleach_column_name,headers,'exact');
nonbleach_cols  = strmatch(nonbleach_column_name,headers,'exact');

time_col        = time_cols(1);
bleach_col      = bleach_cols(1);
nonbleach_col   = nonbleach_cols(1);

X  = data(2:end, time_col)/60;
Y1 = data(2:end, bleach_col);
Y2 = data(2:end, nonbleach_col);

% Look for minima in intensity
Yn = Y1 ./ Y2;

[junk post_frap_time] = min(diff(Yn));
Ymin  = min(Yn);
mean_prebleach      = nanmean( Yn(1:post_frap_time) - Ymin);


Y = (Yn-Ymin) ./ mean_prebleach;

% Find beginning of recovery
i_begin = find ( Y == 0, 1, 'first');
if isempty(i_begin)
    error('MATLAB:fitfrap:BoundNotFound',...
        'Could not find the beginning of recovery (In=0).')
end

% Find end of interesting part, before the slow decay at the end of the
% experiment.
% TODO: so far, nothing special, we take until the end of the movie
i_end = find(X > 10,1);% numel(Y);

% Prepare massaged-data
Xf = X(i_begin : i_end);
Yf = Y(i_begin : i_end);
Xb = X(1:i_begin-1);
Yb = Y(1:i_begin-1);
Xa = X(i_end+1:end);
Ya = Y(i_end+1:end);

% Prepare figure
if show_figure
    hf = figure;
    title(xls_file,'Interpreter','none')
    ha = cla;
    hold(ha,'on');
    plot(ha,Xb,Yb,'.','Color',[0.3 0.3 0.3])
    plot(ha,Xa,Ya,'.','Color',[0.3 0.3 0.3])
    plot(ha,Xf,Yf,'k.')
end

% Do single fit
if fit_single_exp
    
    % Prepare fit options
    exp1_options = fitoptions(exp1_type);
    exp1_options.Lower      = [ 0   0 (1-tolerance)*X(i_begin) ];
    exp1_options.Upper      = [ Inf 1 (1+tolerance)*X(i_begin) ];
    exp1_options.StartPoint = [ 30  0.5 X(i_begin) ];
    
    % Do fit
    [exp1_fit gof1] = fit(Xf,Yf,exp1_type,exp1_options);
    
    % Draw figure
    if show_figure
        plot(Xf,feval(exp1_fit,Xf),'r','LineWidth',1.5)
        text_handle1 = text(0.3*X(end),1,{
            [ '\tau = ' sprintf('%.1f s',exp1_fit.tau) ]
            sprintf('A  = %.2f',exp1_fit.A)
            sprintf('R^2  = %.2f',gof1.rsquare)
            });
        set(text_handle1,'FontName','FixedWidth','VerticalAlignmen','top','Color','r')

    end
    
end

% Do 2 exp fit
if fit_two_exp
    
    % Prepare fit options
    exp2_options = fitoptions(exp2_type);
    exp2_options.Lower      = [ 0   0   0 0 (1-tolerance)*X(i_begin) ];
    exp2_options.Upper      = [ Inf Inf 1 1 (1+tolerance)*X(i_begin) ];
    exp2_options.StartPoint = [ tau1_start tau2_start 0.5 0.5 X(i_begin) ];
    
    % Do fit
    [exp2_fit gof2] = fit(Xf,Yf,exp2_type,exp2_options);
    
    % Draw figure
    if show_figure
        plot(Xf,feval(exp2_fit,Xf),'b','LineWidth',1.5)
        text_handle2 = text(0.6*X(end),1,{
            [ '\tau_1 = ' sprintf('%.1f s',exp2_fit.tau1) ]
            sprintf('A_1  = %.2f',exp2_fit.A1) 
            [ '\tau_2 = ' sprintf('%.1f s',exp2_fit.tau2) ]
            sprintf('A_2  = %.2f',exp2_fit.A2) 
            sprintf('R^2  = %.2f',gof2.rsquare) 
            });
        set(text_handle2,'FontName','FixedWidth','VerticalAlignmen','top','Color','b')
    end
    
end