function BT = DoBootstraps(A, B, C, X, U, E, N)

% Perform bootstrapping for the MAR-1 model
%
% BT = DoBootstraps(A, B, C, X, U, E, N)
% When supplied with the A, B, C, X, U, E matrices, and N, the number of
% bootstraps one requires, this function performs bootstrapping of the
% mar-1 model. The output is a structure, BT, that contains the following
% elements:
% BT.A = the N x p matrix of growth rates A
% BT.B = the p x p x N matrix of variate interaction terms B
% BT.C = the p x q x N matrix of covariate interaction terms C
% BT.Sigma = the p x p xN matrix of Sigma values
% BT.EIGS = the N x p matrixs of eigenvalues
% BT.DetB = the N x 1 vector of determinants of the B matrix
% BT.DetB2p = the N x 1 vector of determinants of B raised to the 2/p
% BT.MaxLambda = the N x 1 vector of maximum lambda values
% BT.MaxLambdaBXB = the N x 1  vector of maximum values of lambda for B kron B
% BT.tr = the N x 1 vector of trace values
% BT.wcReact = the N x 1 vector of 'worst case' reactivities.


p = size(B,1);
q = size(U,2);

% E = randn(size(X,1),p)

boot.A = zeros(N,p);
boot.B = zeros(p,p,N);
boot.C = zeros(p,q,N);
boot.sigma = zeros(p,p,N);
boot.EIGS = zeros(N, p);
boot.DetB = zeros(N,1);
boot.DetB2p = zeros(N,1);
boot.MaxLambda = zeros(N,1);
boot.MaxLambdaBXB = zeros(N,1);
boot.tr = zeros(N,1);
boot.wcReact = zeros(N,1);

% If U is empty, the deployed version doesn't like it. So we generate a new
% one.
if ( ( isdeployed == 1 ) && ( isempty(U) == 1 ) )
    % The only thing we have to do is make it the same size as X and Y, so
    % we make a vector of zeros and then turn it empty.
    U = zeros(size(X,1));
    U(:,1) = [];
end

% 0.9.0 - make sure singular matrix warning is turned off. It should be
% from the LAMBDA function but this doesn't seem to "stick" in the compiled
% version.
warning('off', 'MATLAB:nearlySingularMatrix');
warning('off', 'MATLAB:SingularMatrix');

% 0.9.3 - We need to respect the included variates if requested by users
% (was not doing so before this).

% priori_list = 1 if there are a priori variate/covariate inclusions, 0
% otherwise. 0 is the default.... re-set in the if/then statment if
% includefile exists.
priori_list = 0;

% Is the includefile present? If so, load it (INCLUDED_VARIATES and _COVARS
% will be loaded).
loadwork;

if ( exist(includefile, 'file') == 2 )
    load(includefile);
    priori_list = 1;
end

% Set up a waitbar
wb = waitbar(0, 'Performing bootstrap resampling...');
for n = 1:N
    % Compute Xstar by essentially "simulating" the data with randomly
    % re-sampled process errors (everything else stays the same). This
    % gives us a new X matrix, Xboot.
    % NOTE - X Should be RAW X (not the lagged version) and U should be RAW
    % U (again, not the lagged version).
    Xboot = bootstrapdata(A, B, X, C, U, E);

    % Now we need to make Xstar, which is the bootstrap version of X, i.e.,
    % data from time 1 to T-1, and Ystar, i.e. data from 2 to T.
    Xstar = Xboot(1:end-1,:);
    Ystar = Xboot(2:end, :);

    % We need to have U be from 2 to T, because that is the length of Ystar
    % and Xstar now (they're 1 shorter than before).


    Ustar = U(2:end,:);

    % Now that we have Xstar (a new dataset), we put that into the CLSfit
    % function to obtain new values of A (A1), B (B1) and so on. If priors
    % exist, we are going to have to include those. This checks that.
    
    if (priori_list == 0 )
        % 0.9.3
        % No list of pre-determined include/exclude for variates/covars...
        % test them all.
        [A1, B1, C1, E1, Yhat1, varY1, varDY1] = CLSfit(Xstar, Ystar, Ustar);
    else
        % 0.9.3
        % A list of included vars/covars exists... use it to "seed" the
        % search.
        [A1, B1, C1, E1, Yhat1, varY1, varDY1] = CLSfit(Xstar, Ystar, Ustar, ...
            INCLUDED_VARIATES, INCLUDED_COVARS);
    end

    % Record these in the bootstrap data matrices.
    boot.A(n,:) = A1;
    boot.B(:,:,n) = B1;
    
    % 0.9.0
    % The compiler does not like empty arrays of large dimension so we have
    % to play with this a bit for the deployed version. Since it's already
    % zeros we just leave it... so if it's not empty (isempty is 0, means
    % it is NOT empty), we convert boot.C into C1 for this rep... otherwise
    % we just leave the zeros we started with.
    try
        % We try to do it. If it generates error then we catch it and do
        % nothing.
        boot.C(:,:,n) = C1;
    catch
        boot.C(:,:,N) = boot.C(:,:,N);
    end

    % Find sigma with assessCLSfits
    [sigma1, L, TR2, CR2] = assessCLSfits(E1, size(Xstar,1), size(Xstar,1), varY1, varDY1);

    % The only thing we really need here is sigma1
    boot.sigma(1:p,1:p,n) = sigma1;

    % Stability/Resilience
    [vEIGS1, nDetB1, nDetB2p1] = calcStability(B1);
    boot.EIGS(n,:) = vEIGS1;
    boot.DetB(n,1) = nDetB1;
    boot.DetB2p(n,1) = nDetB2p1;

    % Return rate
    [ maxLambdaB1, BXB1, maxLambdaBXB1 ] = calcReturnRate(B1, vEIGS1);
    boot.MaxLambda(n,1) = maxLambdaB1;
    boot.MaxLambdaBXB(n,1) = maxLambdaBXB1;

    % Reactivity
    [trSV1, wcReact1 ] = calcReactivity(B1, p, sigma1);
    boot.tr(n,1) = trSV1;
    boot.wcReact(n,1) = wcReact1;

    % Update waitbar
    waitbar(n/N, wb);

end  % End the n loop

% Close the wait bar
close(wb);

% and return the value
BT = boot;