function R = testTS(B, dataX, dataY, pnorm,alpha)
%Peforms the two-sample test. Set pnorm = inf for the
%sup-norm statistic
%INPUTS:
% *B - number of Monte-Carlo samples to use to estimate the null distribution.
% *dataX - r1 x d matrix. Data for sample one. Each observation is dx1 vector and there are r1
% replicates.
% *dataY - r2 x d matrix. Data for sample two. Each observation is dx1 vector and there are r2
% replicates.
% *pnorm - norm of the test (e.g. pnorm = inf)
% *alpha - size of the test (e.g. alpha = .05)
%OUTPUTS:
% *R - indicator for rejection of the null hypothesis, equality of mean
% vectors

%%%%%%%%%%%%%%%%%%%%
%Step 1. Estimate the covariance matrices
sigmaX = covshrinkKPM(dataX, 1);

sigmaY = covshrinkKPM(dataY, 1);

%form the linear combo to work with
FullSigma = size(dataX,1)^(-1) * sigmaX + size(dataY,1)^(-1) * sigmaY; 
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%Step 2. Compute test stat
T = statGeneral(dataX, dataY, pnorm);
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%Step 3. Approximate the null distribution
nullv = null_distAC(B,FullSigma,pnorm);
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%Step 4. Decide to reject or not
tail = 0;

[CV R] = BSHTcv(T, nullv, alpha, tail);
%%%%%%%%%%%%%%%%%%%%