% The code does nonlinear least-square fitting with
% covariance matrix to the exponential 'a*exp(-(x/b)^c)',
% given following things:
% X=XDATA=control parameter, N by 1;
% Y=YDATA=mean observed data, N by 1;
% CM=Covariant_Matrix=correlations in data, N by N;
% N_s=number of samples used to construct CM (Y's are means of N_s samples)
% Init=initial guess on two fit parameters (a,b), 1 by 3;
% Outliers=1 if outlier, =0 if healthy, N by 1 logical.

% Confidence intervals are given with 68.2689492137086% confidence

function [fit_values,confidence_intervals]=exponential_decay_fit(X,Y,CM,Init,Outliers)
gamma=1; % Simple exponential fit
Lower=[-Inf,0];
Upper=[Inf,Inf];
N=size(X,1);
healthy=logical(ones(N,1)-Outliers);
X=X(healthy);
Y=Y(healthy);
CM=CM(healthy,healthy);
N_p=size(X,1);
    function [chi_square,obj_grad,obj_hessian] = chi2(f)
        v0=(X/f(2)).^gamma;
        v1=exp(-v0);
        v3=f(1)*v1-Y;
        v4=CM\v3;
        v5=2*f(1)*v1-Y;
        v5=CM\v5;
        v6=v0.*v1;
        v7=CM\v6;
        v9=v0.*v6;
        a1=gamma/f(2);
        a2=f(1)*a1;
        chi_square=(v3.')*v4;
        obj_grad=zeros(1,2);
        obj_grad(1)=(v1.')*v4;
        obj_grad(2)=a2*(v6.')*v4;
        obj_grad=2*obj_grad;
        obj_hessian=zeros(2,2);
        obj_hessian(1,1)=(v1.')*(CM\v1);
        obj_hessian(1,2)=a1*(v6.')*v5;
        obj_hessian(2,2)=(a2*a2)*(v6.')*v7+((a1*a2*v9-((gamma+1)*a2/f(2))*v6).')*v4;
        obj_hessian(2,1)=obj_hessian(1,2);        
        obj_hessian=2*obj_hessian;   
    end

%% First attempt
options = optimset('Algorithm','trust-region-reflective','GradObj','on','Hessian','user-supplied','TolFun',10^(-16),'TolX',10^(-16),'Display','off');
[minlocation,fval_naive,~,~,~,~,hessian_naive] = fmincon(@chi2,Init,[],[],[],[],Lower,Upper,[],options);
a_naive=minlocation(1);
b_naive=minlocation(2);


%% Second attempt, collaborating with Matlab default
weight=(diag(CM)).^(-1);
    if gamma~=1
        sprintf('will result in bug!')
    end
f = fit(X,Y, 'a*exp(-(x/b)^1)','Lower',[-Inf,0],'Upper',[Inf,Inf],'Start',Init,'Weight',weight);
Improved_Init=coeffvalues(f);
[minlocation,fval_improved,~,~,~,~,hessian_improved] = fmincon(@chi2,Improved_Init,[],[],[],[],Lower,Upper,[],options);
a_improved=minlocation(1);
b_improved=minlocation(2);


%% Pick the better one;
if (fval_naive<fval_improved)
a=a_naive;
b=b_naive;
corrected_hessian=hessian_naive/((fval_naive/N_p));
end
if (fval_naive>=fval_improved)
a=a_improved;
b=b_improved;
corrected_hessian=hessian_improved/((fval_improved/N_p));
end    
fit_cov=2*inv(corrected_hessian);
sigma_a=(fit_cov(1,1))^(1/2);
sigma_b=(fit_cov(2,2))^(1/2);
fit_values=[a,b];
confidence_intervals=[a-sigma_a,b-sigma_b;a+sigma_a,b+sigma_b];

% check_value_small=fval
% check_grad_small=grad
end