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")