wilcox_R: What the package does (short line)

cbmhdMC <-
function(x,y,alpha=.05,q=.25,plotit=FALSE,pop=0,fr=.8,rval=15,xlab="",ylab="",nboot=600,SEED=TRUE){
#
#  Compute a confidence interval for the sum of the qth and (1-q)th quantiles
#  of the distribution of D=X-Y, where X and Y are two independent random variables.
#  The Harrell-Davis estimator is used
#  If the distribution of X and Y are identical, then in particular the
#  distribution of D=X-Y is symmetric about zero.
#
#  plotit=TRUE causes a plot of the difference scores to be created
#  pop=0 adaptive kernel density estimate
#  pop=1 results in the expected frequency curve.
#  pop=2 kernel density estimate (Rosenblatt's shifted histogram)
#  pop=3 boxplot
#  pop=4 stem-and-leaf
#  pop=5 histogram
#
library(parallel)
if(SEED)set.seed(2)
if(q>=.5)stop("q should be less than .5")
if(q<=0)stop("q should be greater than 0")
x<-x[!is.na(x)]
y<-y[!is.na(y)]
n1=length(x)
n2=length(y)
m<-outer(x,y,FUN="-")
q2=1-q
est1=hd(m,q)
est2=hd(m,q2)
data1<-matrix(sample(n1,size=n1*nboot,replace=TRUE),nrow=nboot)
data2<-matrix(sample(n2,size=n2*nboot,replace=TRUE),nrow=nboot)
data=cbind(data1,data2)
data=listm(t(data))
bvec=NA
bvec<-mclapply(data,cbmhd_subMC,x=x,y=y,q=q,q2=q2,n1=n1,n2=n2,mc.preschedule=TRUE)
bvec=list2vec(bvec)
p=mean(bvec>0)+.5*mean(bvec==0)
p=2*min(c(p,1-p))
sbv=sort(bvec)
ilow<-round((alpha/2) * nboot)
ihi<-nboot - ilow
ilow<-ilow+1
ci=sbv[ilow]
ci[2]=sbv[ihi]
if(plotit){
if(pop==1 || pop==0){
if(length(x)*length(y)>2500){
print("Product of sample sizes exceeds 2500.")
print("Execution time might be high when using pop=0 or 1")
print("If this is case, might consider changing the argument pop")
print("pop=2 might be better")
}}
MM=as.vector(m)
if(pop==0)akerd(MM,xlab=xlab,ylab=ylab)
if(pop==1)rdplot(MM,fr=fr,xlab=xlab,ylab=ylab)
if(pop==2)kdplot(MM,rval=rval,xlab=xlab,ylab=ylab)
if(pop==3)boxplot(MM)
if(pop==4)stem(MM)
if(pop==5)hist(MM,xlab=xlab)
if(pop==6)skerd(MM)
}
list(q=q,Est1=est1,Est2=est2,sum=est1+est2,ci=ci,p.value=p)
}