% ch13_nonlinear_fit.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% Fit exponential decay y(t) = A*exp(-t/tau) using MATLAB's fminsearch
% (Nelder-Mead simplex minimization). No optimization toolbox is required.
% Replace the example data and starting estimates below.
% Fit the original y measurements directly, without taking logarithms of y.
% Assumptions: independent errors in y and negligible errors in t.
% The routine minimizes in A and log(tau), ensuring a positive decay time.
% Solver stopping tests do not establish uniqueness or physical validity.
% Try different initial estimates and inspect the residuals.

clear; close all; clc;

%% Example data and settings
seed = 1;
rng(seed, 'twister');
% Alternative: comment out the line above and uncomment the following.
% rng('shuffle');
t = linspace(0,5,30)';
y = 2.0*exp(-t/1.5) + 0.1*randn(size(t));
sigma_y = 0.1*ones(size(t)); % Known standard deviations; [] if unknown
A0 = 1.0;
tau0 = 1.0;                % Must be positive
max_iter = 2000;
max_evaluations = 4000;
tol_parameters = 1e-8;     % Tolerance in A and log(tau), not in tau itself
tol_objective = 1e-8;     % Absolute tolerance in the minimized sum
make_plot = true;

% For measured data, replace the example above, for instance:
% data = readmatrix('measurements.txt');
% t = data(:,1); y = data(:,2); sigma_y = data(:,3);
% Use sigma_y = [] for a two-column file without known uncertainties.

%% Fit
options = optimset('Display','off','MaxIter',max_iter, ...
    'MaxFunEvals',max_evaluations,'TolX',tol_parameters,'TolFun',tol_objective);
[A,tau,info] = exponential_fit(t,y,sigma_y,A0,tau0,options);
fprintf('\nNonlinear exponential fit: y = A*exp(-t/tau)\n');
fprintf('Amplitude A:        %.12g\n',A);
fprintf('Decay time tau:     %.12g\n',tau);
fprintf('Residual sum of squares: %.10g\n',info.sse);
fprintf('Residual standard deviation: %.6g\n',info.residual_std);
fprintf('Degrees of freedom: %d\n',info.dof);
if ~isempty(sigma_y)
    fprintf('Chi-squared:         %.6g\n',info.objective);
    fprintf('Reduced chi-squared: %.6g\n',info.objective/info.dof);
end
fprintf('Iterations:         %d\n',info.iterations);
fprintf('Objective calls:    %d\n',info.evaluations);
fprintf('Solver exit flag:   %d\n',info.exitflag);
fprintf('Status: %s\n',info.message);
% This program reports point estimates, not parameter uncertainties.

%% Optional data, fitted curve, and residual plots
if make_plot
    tt = linspace(min(t),max(t),300);
    figure('Color','w');
    subplot(2,1,1);
    if isempty(sigma_y)
        hd = plot(t,y,'ko','MarkerFaceColor','w','MarkerSize',6);
    else
        hd = errorbar(t,y,sigma_y,'ko','MarkerFaceColor','w','MarkerSize',6);
    end
    hold on;
    hf = plot(tt,A*exp(-tt/tau),'k-','LineWidth',1.5);
    xlabel('t','FontSize',14); ylabel('y','FontSize',14);
    title('Nonlinear exponential fit','FontSize',14,'FontWeight','normal');
    legend([hd,hf],{'Data','Fitted curve'},'Location','best');
    set(gca,'FontSize',14); grid on; box on;
    subplot(2,1,2);
    if isempty(sigma_y)
        plot(t,info.residuals,'ko','MarkerFaceColor','w','MarkerSize',6);
    else
        errorbar(t,info.residuals,sigma_y,'ko','MarkerFaceColor','w','MarkerSize',6);
    end
    hold on; plot([min(t) max(t)],[0 0],'k--');
    xlabel('t','FontSize',14); ylabel('y - y_{fit}','FontSize',14);
    set(gca,'FontSize',14); grid on; box on;
end

%% Reusable fitting routine
function [A,tau,info] = exponential_fit(t,y,sigma_y,A0,tau0,options)
% Copy this routine and its helper to exponential_fit.m for reuse.
% sigma_y = [] gives unweighted least squares. Otherwise minimize the sum
% of squared residuals divided by the known measurement variances.
    validateattributes(t,{'double'},{'real','vector','finite','nonnegative'});
    validateattributes(y,{'double'},{'real','vector','finite'});
    validateattributes(A0,{'double'},{'real','scalar','finite'});
    validateattributes(tau0,{'double'},{'real','scalar','finite','positive'});
    t=t(:); y=y(:); n=numel(t);
    if n<3 || numel(y)~=n || max(t)==min(t)
        error('Supply at least three measurements with distinct time values.');
    end
    if isempty(sigma_y)
        scale=ones(n,1);
    else
        validateattributes(sigma_y,{'double'},{'real','vector','finite','positive'});
        scale=sigma_y(:);
        if numel(scale)~=n
            error('Supply one uncertainty for each measurement.');
        end
    end
    cost=@(q) decay_cost(q,t,y,scale);
    [q,value,exitflag,output]=fminsearch(cost,[A0;log(tau0)],options);
    A=q(1); tau=exp(q(2));
    if ~isfinite(value) || ~isfinite(tau) || tau<=0
        error('The solver did not return a finite usable fit.');
    end
    residuals=y-A*exp(-t/tau);
    sse=sum(residuals.^2);
    info=struct('residuals',residuals,'sse',sse,'dof',n-2, ...
        'residual_std',sqrt(sse/(n-2)),'objective',value, ...
        'iterations',output.iterations,'evaluations',output.funcCount, ...
        'exitflag',exitflag,'message',output.message);
end

function value=decay_cost(q,t,y,scale)
    tau=exp(q(2));
    if any(~isfinite(q)) || ~isfinite(tau) || tau<=0
        value=Inf;
        return;
    end
    residuals=(y-q(1)*exp(-t/tau))./scale;
    value=sum(residuals.^2);
    if ~isfinite(value)
        value=Inf;
    end
end