% Modified version of the MATLAB template file FitNumAlgMR.m that was
% modified for the solution of Example 4 of Supplemental Unit S4 of "A
% First Course on Kinetics and Reaction Engineering."
%
function S4_Example_4(p_guess)
    % Set variables
    x = [0.4809	36.0
    0.4809	143.0
    0.4809	296.0
    0.4809	561.0
    0.4725	25.0
    0.4725	185.0
    0.4725	335.5
    0.4725	462.0];
    
    % Experimental response variables
    y_hat = [0.4188	0.0063	0.0298
    0.2882	0.0046	0.1200
    0.2188	0.0027	0.1682
    0.0868	0.0011	0.2578
    0.4163	0.0067	0.0279
    0.2629	0.0043	0.1334
    0.1598	0.0021	0.2056
    0.1074	0.0014	0.2326];
    
    % Function that evaluates the model equations
    function g = nlaeqns(u)
        % Current parameter values are available in column vector p
        % Current set variables are available in column vector x_set
        x1 = x_set(1);
        x2 = x_set(2);
        p1 = 10^p(1);
        p2 = 10^p(2);
        p3 = 10^p(3);
        s1 = p2 + p3;
        s2 = p1-p2-p3;
		g = [
			u(1) - x1*exp(-p1*x2);
            u(2) - p1*x1*((exp(-s1*x2)-exp(-p1*x2))/s2);
            u(3) - p2*x1*(1/s1+exp(-p1*x2)/s2-p1*exp(-s1*x2)/s1/s2);
		];
    end % of internal function nlaeqns

    % Function that calculates the responses using the model solutions
    function y = mrmodel(p_current)
        % Make the current parameters available to the model equations
        p = p_current;
        y = zeros(n_data,n_resp);
        
        % loop through the data points
        for i = 1:n_data
            % Get the set variables
            for j = 1:n_set
                x_set(j) = x(i,j);
            end
            
            % provide guesses for the solution to the model equations
            u_guess = [
                y_hat(i,1);
                y_hat(i,2);
                y_hat(i,3);
            ];
            
            % Solve the set of nonlinear model equations
            % u = fsolve(@nlaeqns, u_guess);
            
            % The next line suppresses the fsolve message
            options = optimoptions('fsolve','Display','off');
            u = fsolve(@nlaeqns, u_guess, options);

            % Use the solution to calculate the responses
            y(i,:) = [
                u(1);
                u(2);
                u(3);
            ];
        end
    end % of internal function mrmodel

    % Function that evaluates the objective function
    function phi = objfcn(p_curr)
        % calculate the model responses using the current parameter values
        y = mrmodel(p_curr);
        
        % calculate the errors
        err = y_hat - y;
        
        % form the error matrix
        Emat = zeros(n_resp);
        for i = 1:n_resp
            % diagonal element
            for k = 1:n_data
                Emat(i,i) = Emat(i,i) + err(k,i)^2;
            end
            % below diagonal elements
            for j = i+1:n_resp
                for k = 1:n_data
                    Emat(j,i) = Emat(j,i) + err(k,j)*err(k,i);
                end
            end
            % above diagonal elements
            for j = i+1:n_resp
                Emat(i,j) = Emat(j,i);
            end
        end
        
        % evaluate the objective function
        phi = det(Emat);
    end % of internal function objfcn

    % All the necessary internal functions are now defined

    % get the number of data, parameters response and set variables
    n_data = size(y_hat,1);
    n_param = length(p_guess);
    n_resp = size(y_hat,2);
    n_set = size(x,2);
    
    % create global variables to make parameters and set variables
    % available to the internal function that evaluates the model equations
    p = zeros(n_param,1);
    x_set = zeros(n_set,1);

    % minimize the objective function and compute the responses
    y = zeros(n_data,n_resp);
    pf = fminsearch(@objfcn,p_guess);
    y = mrmodel(pf);
    resid = y - y_hat;
    
    % create a parity plot for each response variable
    for i = 1:n_resp
        figure;
        plot(y_hat(:,i),y_hat(:,i),'r',y_hat(:,i),y(:,i),'o')
        xlabel('Measured Response')
        ylabel('Predicted Response')
        title(['Parity Plot for Response variable ',num2str(i)])
        
        % create a residuals plot for this response variable versus each
        % set variable
        for n = 1:n_set
            figure;
            plot(x(:,n),zeros(size(y,1),1),'r',x(:,n),resid(:,i),'o')
            xlabel(['x',num2str(n)])
            ylabel('Residual')
             title(['y',num2str(i),' Residuals Plot vs. x',num2str(n)])
       end
    end

    % Display the parameters 
    display('Best Values for the Parameters:')
    pf
end % of S4_Example_4.m