function quahogs_driver
% January 12, 2024

    format short
    format compact
    
    [t,y] = quahogs_data;
    
    % fitnlm
    
        fun = @(param,xval) param(1) + param(2) .* exp(param(3) .* xval);
        b1_init = prctile(y(:),50);
        beta0 = [b1_init 0 0];
        mdl = fitnlm(t,y,fun,beta0) 
    
    % varpro
     
        fun = @(lambda) model(t,lambda);
        w = ones(size(y));
        lambda0 = -20;
        n = 2;  
        [~,~,~,~,~,Regression] = varpro(y, w, lambda0, n, fun);
        t_ratio = Regression.t_ratio
        std_param = Regression.std_param
        fprintf(newline)
        fprintf(newline)
    
    % quahogs fit, modified expfitdemo from NCM
     
        fun = @(lambda) expfitfun(lambda,t,y);
        lambda = fminbnd(fun,-40,0,[])
        X = [ones(size(t)) exp(lambda*t)];
        condX = cond(X)
        beta = X\y
        normres = norm(X*beta-y)
        
    % quahogs plot
    
        quahogs_plot(lambda,t,y)
    
    % ------,~,------------------------------------

    function [t,y] = quahogs_data
    
       AT = [1833 1833 1833 1833 1833 1833 1833 1862 1862 1862 1862 1862 1856 ...
             1856 1856 1856 1856 1856 1856 2063 2063 2063 2063 2063 2063]';
    
       DIC = [1598 1598 1598 1598 1598 1598 1598 1684 1684 1684 1684 1684 1749 ...
              1749 1749 1749 1749 1749 1749 2071 2071 2071 2071 2071 2071]';
         
       Calc_rate = ...
             10000\[1.08236  0.33303  1.33213  0.66607  0.41629  0.54118 ...
                    0.24978  0.58281  0.49955  0.37466  0.74933  0.29140 ...
                    0.58281 -0.16652  0.29140  0.62444  0.49955  0.37466 ...
                    0.04163 -0.95747 -1.37376 -1.24888 -0.66607 -0.79095 ...
                   -0.70770]';
    
       t = AT ./ DIC;
       t = t - mean(t);      % Center
       y = 10000*Calc_rate;  % Scale 

    end

    function [phi,dphi,ind] = model(t,lambda)
       phi = [ones(size(t)) exp(lambda*t)];
       dphi = t.*exp(lambda*t);
       ind = [2;1];
    end
    
    function res = expfitfun(lambda,t,y)
        X = [ones(size(t)) exp(lambda*t)];
        beta = X\y;
        res = norm(X*beta-y); 
    end
    
    function errest = errorestimates(lambda,t,y,s)
        X = [ones(size(t)) exp(lambda*t)];
        [Q,R] = qr(X);
        beta = R\(Q'*y);
        r = y - X*beta;
        sig = 0.2*norm(r);
        Xs = [ones(size(s)) exp(lambda*s)];
        E = Xs/R(1:2,1:2);
        e = sig*sqrt(1+diag(E*E'));
        f = Xs*beta;
        errest = [f-e f+e];
    end

    function quahogs_plot(lambda,t,y)
        s = (min(t):(max(t)-min(t))/64:max(t))';
        z = beta(1) + beta(2)*exp(lambda*s);
        err = errorestimates(lambda,t,y,s);
    
        shg
        set(gcf,'menu','none','numbertitle','off') 
        h = plot(t,y,'.',s,z,'-');
        xlabel('t')
        ylabel('y')
        axis([-0.1 0.1 -2.0 2.0])
        h(1).MarkerSize = 24;
        h(2).LineWidth = 2;
        title(['\lambda = ' num2str(lambda)], ...
            'fontsize',16, ...
            'fontweight','bold')
        
        pcolor = get(h(2),'color'); 
        patch([s;flip(s)], [err(:,1);flip(err(:,2))], pcolor, ...
            'facealpha',0.25, ...
            'edgecolor','w');
    end

end