% MATLAB template file for fitting a model to multiple response 
% experimental data of the form
% (y_hat(1),...,y_hat(n_resp),x(1),...,x(n_set)). Explicit expressions, 
% y(1) = f1(x(1),...,x(n_set),u(1),...,u(n)),
% y(2) = f2(x(1),...,x(n_set),u(1),...,u(n)),etc. are used to calculate the 
% model-predicted responses from the experimental set variables (x) and the
% results (u) of solving a set of n model equations of the form
%   du(1)/dv = g1(v,u(1),...,u(n);p(1),...,p(n_par)); u(1)(v_0) = u1_0
%   du(2)/dv = g2(v,u(1),...,u(n);p(1),...,p(n_par)); u(2)(v_0) = u2_0
%              .
%              .
%              .
%   du(n)/dv = gn(v,u(1),...,u(n);p(1),...,p(n_par)); u(n)(v_0) = un_0
% where p represents the parameters in the model.
%
% See Supplemental Unit S4 of "A First Course on Kinetics and Reaction
% Engineering" and/or the MATLAB documentation for additional information
% or if questions arise.
%
% Carl R. F. Lund, 4/22/13
%   last rev. 5/28/14
%
% Arguments
%   p_guess (vector) - guesses for the model parameters
%
% This file requires six modifications before each use; the locations
% where editing is required are indicated by the comment "% EDIT HERE"
%
function FitNumDifMR(p_guess)
    % Known quantities and constants (in consistent units)
% EDIT HERE (Required modification 1 of 6):
    % define universal and experimental constants here

% EDIT HERE (Required modification 2 of 6) enter the set variable data
% as a matrix named x with one column for each set variable and one row per
% data point

% EDIT HERE (Required modification 3 of 6) enter the measured response
% variable data as a matrix named y_hat with one column for each 
% response variable and one row per data point
    
	% Function that evaluates the ODEs
    function dudv = odeqns(v,u)
        % Current parameter values are available in column vector p
        % Current set variable values are available in column vector x_set
% EDIT HERE (Required modification 4 of 6):
        dudv = [
            % Evaluate du1/dv = g1(v,u(1),...,u(n);p(1),...,p(n_par)) here
			% Evaluate du2/dv = g2(v,u(1),...,u(n);p(1),...,p(n_par)) here
            % and so on through gn, one per line
        ];
    end % of internal function odeqns

    % Function that solves the model equations and calculates the model-
    % predicted responses
    function y = mrmodel(p_current)
        % Declare y
        y = zeros(n_data,n_resp);
        % Make the current parameters available to the model equations
        p = p_current;
        
        for i = 1:n_data
            % Make the set variables available to the model equations
            for j = 1:n_set
                x_set(j) = x(i,j);
            end
            
            % Initial and final values
% EDIT HERE (Required modification 5 of 6):
            v0 = ; % insert the independent variable initial value here
            u0 = [
                % insert the initial value of dependent variable u1 here
                % insert initial values for u2, u3, ..., one per line
            ];
            vf = ; % insert the independent variable final value here
            
            % Solve the set of ordinary differential model equations
            [v,uu] = ode45(@odeqns,[v0 vf],u0);
% NOTE: For stiff equations replace ode45 with ode15s
            last_value = length(v);
            u = uu(last_value,:);
            
            % Use the solution to calculate the responses, y(1), y(2), ...
            % by evaluating f1(u(1),...,u(n),x_set(1),...,x_set(n_set)),
            % f2(u(1),...,u(n),x_set(1),...,x_set(n_set)), etc.
% EDIT HERE (Required modification 6 of 6):
            y(i,:) = [
                % evaluate response variable y(1) here
                % evaluate response variable y(2) here
                % and so on, one response variable per line
            ];
        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 FitNumDifMR.m