install.packages("foreach")
install.packages("doParallel")

######################################################################################prepare
rm(list = ls())

if(TRUE){
library(foreach)#并行计算
library(doParallel)#实现多核并行
no_cores<- detectCores(logical=FALSE)+15
cl<- makeCluster(no_cores)
registerDoParallel(cl)
}

###function:compute q2
q2_estimation<- function(xi_0,lam,m0,p,nn,k,xm,ym,nm){

epsilon1=0.001
sigma<- matrix(0,nrow=p+1,ncol=p+1)
xi<- rep(0,p+1)
z<- rep(0,m0*nn)
zeta<- rep(0,m0*nn)
w<- rep(0,m0*nn)
ww<- matrix(0,nrow=m0*nn,ncol=m0*nn)

l=0
flag=0
while(flag==0){
l<- l+1
xi0<- xi
for(i in 1:m0){
for(j in 1:nn){
s<- (i-1)*nn+j
zeta[s]<- sum(xm[s,]*xi)
bp<- 1.0/(exp(-zeta[s])+1.0)
z[s]<- zeta[s]+(ym[i,j]-nm[i,j]*bp)/(nm[i,j]*bp*(1.0-bp))
w[s]<- lam*((1.0-lam)^(m0-i))*nm[i,j]*bp*(1.0-bp)
ww[s,s]<- w[s]
}
}
sigma<- t(xm)%*%ww%*%xm
xi<- solve(sigma)%*%t(xm)%*%ww%*%z
print(xi)
###
absxi<- xi-xi0
s1=0.0
s2=0.0
for(i in 1:(p+1)){
s1<- s1+(absxi[i])^2
s2<- s2+(xi0[i])^2
}
eps<- (sqrt(s1))/(sqrt(s2))
###
if(eps<epsilon1 || l==30) flag=1
}
result<- (2-lam)/lam*t(xi-xi_0)%*%sigma%*%(xi-xi_0)
return(result)
}


###function: compute run length
rlfun<- function(mu0,cov0,alpha,belta,lam,m0,p,a1,b1,nn,climit,tau,shift,case){
library(MASS)
kmax=3000
xi_0<- rep(0,p+1)
xi_0[1]<- alpha
xi_0[2:(p+1)]<- belta[1:p]

alpha1<- alpha+shift
belta1<- belta
belta1[1]<- belta[1]+shift
mu1<- mu0
mu1[1]<- mu0[1]+shift

x<- matrix(0,nrow=nn,ncol=p)
y<- rep(0,nn)
n<- rep(0,nn)
ew_x<- rep(0,p)

xs<- matrix(0,nrow=kmax*nn,ncol=p)
ys<- matrix(0,nrow=kmax,ncol=nn)
ns<- matrix(0,nrow=kmax,ncol=nn)

k=0
flag=0
rl=0
while(k<kmax && flag==0){
k<- k+1
print(k)
x<- mvrnorm(nn,rep(0,p),cov0)

for(i in 1:nn){
if(k<=(m0+tau)){
x[i,]<- x[i,]+mu0
ab<- alpha+sum(x[i,]*belta)
}else{
switch (case,
{
x[i,]<- x[i,]+mu0
ab<- alpha1+sum(x[i,]*belta)
},
{
x[i,]<- x[i,]+mu0
ab<- alpha+sum(x[i,]*belta1)
},
{
x[i,]<- x[i,]+mu1
ab<- alpha+sum(x[i,]*belta)
}
)
}
bp<- 1/(exp(-ab)+1.0)
n[i]=5
y[i]<- rbinom(1,n[i],bp)
xs[(k-1)*nn+i,]<- x[i,]
ys[k,i]<- y[i]
ns[k,i]<- n[i]
}

if(k>m0){
m<- k-m0
ew_x<- mu0
xm<- matrix(1.0,nrow=m0*nn,ncol=p+1)
ym<- matrix(0,nrow=m0,ncol=nn)
nm<- matrix(0,nrow=m0,ncol=nn)
xbar<- rep(0,p)
for(i in 1:m0){
xm[((i-1)*nn+1):(i*nn),2:(p+1)]=xs[((m+i-1)*nn+1):((m+i)*nn),1:p]
ym[i,]=ys[m+i,]
nm[i,]=ns[m+i,]
for(j in 1:p){
xbar[j]<- 1.0*sum(xs[((m+i-1)*nn+1):((m+i)*nn),j])/nn
}
ew_x<- (1.0-lam)*ew_x+lam*xbar
}
q1<- nn*((2-lam)/lam)*(t(ew_x-mu0))%*%(solve(cov0))%*%(ew_x-mu0)
q2<- q2_estimation(xi_0,lam,m0,p,nn,k,xm,ym,nm)
q<- q1+q2
if(q>climit) flag=1
}
}
rl<- k
#print(rl)
return(rl)
}
###################################################################################main project
#######################################################generate data

arl0=200
smax=500
lam=0.2
delta<- c(0.0,0.05,0.06,0.07,0.1,0.2,0.3,0.5)
#delta<- c(0.0,0.15,0.2,0.25,0.3,0.6,1.0,1.5)
#delta<- c(0.0,0.025,0.03,0.035,0.07,0.12,0.15,0.2)
p=2
nn=20
m0=40
tau=40
case=1

a1=1
b1=30
h0=15.95
#alpha=-2.2
alpha=-2.94
belta=rep(0,p)
belta[1]=1.0
belta[2]=2.0
#belta[2]=1.0

mu0<- rep(0,p)
cov0<- matrix(0,nrow=p,ncol=p)
for(i in 1:p) {
for(j in 1:p){
#cov0[i,j]<- ifelse(i==j,0.1,0.0)
cov0[i,j]<- ifelse(i==j,1.0,0.0)
#cov0[i,j]<- ifelse(i==j,1,0.5^abs(i-j))
}
}
#print(cov0)

###############################################find control limit (arl0=370, lambda=0.1)
epsilon=1.0
hl=h0-0.5
hu=h0+0.5

arl=0
while(abs(arl-arl0)>epsilon){
print("***********************************")
climit=(hl+hu)/2
print(climit)
count=0
arl=0
rlvec=rep(0,smax)
rlvec<- foreach(s=1:smax, .combine="c") %dopar%
rlfun(mu0,cov0,alpha,belta,lam,m0,p,a1,b1,nn,climit,tau=0,shift=0,case=1)

count=0
rlvecc=rep(0,smax)
for(s in 1:smax){
if(rlvec[s]>m0){
count=count+1
rlvecc[count]=rlvec[s]-m0
#print(rlvecc[count])
}
}
arl=mean(rlvecc[1:count])
print(arl)
ic_out<-cbind(arl,climit)
if(arl>=arl0){
hu=climit
}else{
hl=climit
}
}
write.table(ic_out,file="C:/Users/Administrator/Desktop/n=5_ic_climits.txt",row.names=F,quote=F,sep="\t")

###################################################compute arl1 (arl0=370, lambda=0.1)
arl1<- rep(0,length(delta))
sdrl1<- rep(0,length(delta))

for(i in 1:length(delta)){
shift=delta[i]
rlvec=rep(0,smax)
rlvec<- foreach(s=1:smax, .combine="c") %dopar%
rlfun(mu0,cov0,alpha,belta,lam,m0,p,a1,b1,nn,climit,tau,shift,case)

count=0
rlvecc=rep(0,smax)
for(s in 1:smax){
if(rlvec[s]>(m0+tau)){
count=count+1
rlvecc[count]=rlvec[s]-(m0+tau)
#print(rlvecc[count])
}
}
#write.table(rlvecc[1:count],file="C:/Users/Administrator/Desktop/n=5_oc_rls.txt",row.names=F,quote=F)
#arl[i]=1.0*sum(rlvecc[1:count])/count
#sdrl[i]=sqrt(1.0*sum((rlvecc[1:count]-arl)^2)/count)
arl1[i]=mean(rlvecc[1:count])
sdrl1[i]=sd(rlvecc[1:count])
print("***********************************")
print(arl1[i])
print(sdrl1[i])
}

oc_out<- cbind(arl1,sdrl1)
write.table(oc_out,file="C:/Users/Administrator/Desktop/n=5_oc_arl&sdrl.txt",row.names=F,quote=F,sep="\t")