function [est_param, est_mixing, h_fun, ls_err] = deconv_fluor(F,G,xt, n_fluoro, init_param_guess, method_condlinear, xt_weight)
%[est_param, est_mixing, h_fun, ls_err] = deconv_fluor(F,G,xt, n_fluoro, init_param_guess, method_condlinear, xt_weight)
% INPUT:
%   F                   Convoluted fluorescence (measured)
%   G                   Instrument impule response function (measured)
%   t                   Time
%   n_fluoro            Number of fluorophores
%   init_param_guess:   Initial guess for the non-linear parameters (tau)
%   method_condlinear:  Method for solving the conditional linear
%                       subproblem (e.g. mixing coefficients). Options are:
%                       'nnls'  non-negative least squares
%                       'ls' least squares
%                       'wnnls'  weighted non-negative least squares
%                       'wls' weighted least squares
%                       Default is 'nnls'.
%   xt_weights          Weights (Tx1) given in connection with 'wnnls' or
%                       'wls'. Default is [], not weighted.
%
%
% OUTPUT:
%   h   fluorescence impulse response function (estimated)
%
%
%  F(t_i) = int[0,t_i](G(t_i-s)h(s) ds

% Least squares approach (one emission and one excitation)
% y = A c + epsilon, (y is Tx1, A is TxP, c is Px1)
% solve for c=A\y
% Note, y  is equally spaced F (observed)
%       A  is h(t,theta) calculated with the current theta-parameter values
%       (theta is either known a priori or are what we are optimizing
%       non-linearly
%       c is the "mixing" coeffients and are solved for.

% Copyright, 2023 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. 
%
% Jesper L?ve Hinrich
% Department of Food Science
% Faculty of Sciences
% University of Copenhagen
% Rolighedsvej 26, 1958 Frederiksberg C, Denmark

if nargin < 6 || isempty(method_condlinear)
    method_condlinear = 'nnls';
    xt_weight = [];
end
if nargin < 7 || isempty(xt_weight)
    xt_weight = [];
    if any(strcmpi(method_condlinear,{'wls','wnnls'})) && isempty(xt_weight)
        error('Weights must be specified when using ''wls'' or ''wnnls''.' )
    end
end


%% Initialization
[h_fun, h_init_param] = getHfunction('multiexponential');
if nargin < 5 || isempty(init_param_guess)
    param = initialize_tau_by_slope(F,xt,n_fluoro);
else
    param = init_param_guess(:)';
    n_fluoro=length(param);
end


%% Search for param based on the following and the error
opt_fun = 'fminsearch';

if strcmpi(opt_fun, 'fminsearch')
    options=optimset('Display','off','TolX',1e-4,'TolFun',1e-4,'MaxFunEvals',100,'MaxIter',100); %
elseif strcmpi(opt_fun, 'fmincond')
    options =  optimoptions('fmincon');
else
    error('Unknown optimization function "%s"!',opt_fun)
end

fit_param_simultaneously = true;
est_param = param;
maxiter=50;
conv_crit = 1e-8;
delta_rel_cost = inf;
old_cost = sum(F(:).^2);

iter=0;

while delta_rel_cost >= conv_crit && iter <= maxiter
    iter = iter+1;
    
    if fit_param_simultaneously
        if length(unique(est_param)) ~= n_fluoro
            est_param = est_param.*(randn(size(est_param))*0.1+1);
        end
        
        jfitfun = @(p) model_fitting(xt,F,G,h_fun,p, method_condlinear, xt_weight);
        if strcmpi(opt_fun, 'fminsearch')
            
            [est_param,err,ExitFlag,OutPut] = fminsearch(jfitfun,est_param,options);
        elseif strcmpi(opt_fun, 'fmincond')
            [opttxt,est_param,err,ExitFlag,OutPut] = evalc('fmincon(jfitfun,est_param,[],[],[],[],zeros(1,length(param)),inf(1,length(param)),[], options);');
        end
        [ls_err,est_mixing] = jfitfun(est_param);
        
        
    else
        param = param(randperm(n_fluoro)); 
        for i = 1:length(param)
            sub_param = param(1:i);
            fixed_param = -inf*ones(1,length(param)-i);
            sub_jfit = @(p) model_fitting(xt,F,G,h_fun,[p,fixed_param], method_condlinear, xt_weight);
            if strcmpi(opt_fun, 'fminsearch')
                [est_param_sub,err,ExitFlag,OutPut] = fminsearch(sub_jfit,sub_param,options);
                
            elseif strcmpi(opt_fun, 'fmincond')
                
                [opttxt,est_param_sub,err,ExitFlag,OutPut] = evalc('fmincon(sub_jfit,sub_param,[],[],[],[],zeros(1,length(sub_param)),Inf(1,length(param)-length(sub_param)),[], options);');
            end
            param(1:i) = est_param_sub;
        end
        est_param = param;
        [ls_err,est_mixing] = model_fitting(xt,F,G,h_fun,est_param);
        
    end
    cost = ls_err;
    delta_rel_cost = abs(cost-old_cost)/abs(cost); 
    old_cost = cost;    
      
end

%% Output results
[est_param, idx_sort] = sort(est_param, 'ascend');
est_mixing = est_mixing(idx_sort,:);

end

function [err_ls, est_mixing] = model_fitting(xt,Ft,Gt,ht,g_nl_param, opt_method, w)
% Solves for the conditional linear parameters using either
% non-negative least squares or least squares or their weighted
% versions.
% Input:
%   xt:     Time
%   Ft:     Observed ?fluorescence?
%   Gt:     Instrument Impulse Response Function
%   ht:     A function taking (xt,init_param) that gives a matrix
%   g_nl_param: Guess for the non-linear parameters (tau in the
%               multi-exponential.
%
% Output:
%   err_ls:     Error after solving the least squares problem
%   est_mixing  The estimated mixing coefficients

if nargin < 6 || isempty(opt_method)
    opt_method = 'nnls';
    w = [];
end
y = Ft;

A = ht(xt,g_nl_param);
% Convolute the signal
A2 = zeros(size(A));
for i = 1:size(A,2)
    tmp = conv(Gt,A(:,i));
    A2(:,i) = round(tmp(1:size(A2,1)));
end
A=A2;

% Find mixing coefficients
if strcmpi(opt_method,'nnls')
    c = fcnnls(A,y);
elseif strcmpi(opt_method,'ls')
    c = A\y;
elseif strcmpi(opt_method,'wnnls')
    c = fcnnls([],[],(A'*diag(w)*A),A'*diag(w)*y);
    
elseif strcmpi(opt_method,'wls')
    c = (A'*diag(w)*A)\A'*diag(w)*y;
else
    error('Method "%s" is unknown.',opt_method)
end
est_mixing=c;

% Error
if contains(opt_method,'w')
    err_ls = sqrt(sum(w.*sum((y-A*c).^2,2))); %sqrt of Eq. 2 in Grinvald & Steinberg
    %   weight is predefined, by expert knowledge or experimentally measured  (e.g. statistical error)
    
else
    err_ls = norm(y-A*c,'fro');
end
end




