.Random.seed<-c(9,20,10,33,13,3,61,59,10,2,25,10)
library(yags)
source("myags.ssc")
source("myags.pmat.fit.mancl.ssc")
source("saws2final.ssc")
### modify the LZ.exchalp function by V. Carey to new one
source("LZ.exchalp.ssc")

### Set initial values

nsim<-1000
K<-20
MAXITER<-25
#CORR<-"exchangeable"
CORR<-"independence"
#FIXEDTREAT<-T
FIXEDTREAT<-F
SCALEFUN<-LZ.scalefun
#SCALEFUN<-one.scalefun


## Start main Program

one.scalefun<-function(R,ID,cor.met,p){ 1 }


now<-proc.time()
make.pois.data<-function(K,n,theta,tau,rcluster=T,
              fixedtreat=FIXEDTREAT){
if (rcluster){
	varn<- 2*n
   ni<-ceiling(rgamma(2*K,n^2/varn,n/varn))
 if (FIXEDTREAT){ 
   ## set ni[c(2,4,6,...,K)] to zero
   ## i.e., first K/2 clusters treat=-1 only
    ni[2*(1:(K/2))]<-0
   ## then set ni[c(K+1,K+3,...2K-1)] to zero
   ## i.e., last K/2 clusters treat=1 only
    ni[K-1 + 2*(1:(K/2))]<-0
}
	y<-strata<-treat<-rep(NA,sum(ni))
	for (i in 1:K){
	   NI<-sum(ni[(2*i-1):(2*i)])
		if (i==1) index<-1:NI
		else index<-sum(ni[1:(2*i-2)]) + 1:NI
		strata[index]<-rep(i,NI)
		treat[index]<-c(rep(-1,ni[2*i-1]),rep(1,ni[2*i]))
    if (tau==0)		y[index]<-rpois(NI,theta)
    else 		y[index]<-rpois(NI,exp(log(theta)+treat*rnorm(1,mean=0,sd=tau) ))
	}
}
else{
#	if (tau==0) thetai<-rep(theta,K)
#	else 	thetai<- theta*tau*rgamma(K,1/tau )
		#print(thetai)
		treat<-c(rep( c(rep(-1,n/2),rep(1,n/2)), K ) )
		#,rep(0,n*(K-5)))
	y<-strata<-rep(NA,K*n)
	for (i in 1:K){
		index<-(i-1)*n + 1:n
		strata[index]<-rep(i,n)
if (tau==0)		y[index]<-rpois(n,theta)
else 		y[index]<-rpois(n,exp(log(theta)+treat*rnorm(1,mean=0,sd=tau) ))
	}
}
	data.frame(Y=y,strata=strata,treat=treat)
}

#temp<-make.pois.data(6,4,10,.5,rcluster=T)



Nt<-8

 pa1<-pa4<-pd1<-pd4<-matrix(NA,nsim,Nt)
 cia1<-cia4<-cid1<-cid4<-matrix(NA,nsim,Nt)
 nrejecta1<-nrejecta4<-nrejectd1<-nrejectd4<-rep(0,Nt)
 va1<-va4<-vd1<-vd4<-matrix(NA,nsim,4,dimnames=list(NULL,c("Vm","Vs","VsMD","Va")))
 wa1<-wa4<-wd1<-wd4<-rep(NA,nsim)
 ba1<-ba4<-bd1<-bd4<-rep(NA,nsim)

 converge.countd1<-converge.countd4<-0
 errors.d1<-errors.d4<-erri.d1<-erri.d4<-NULL

 for (i in 1:nsim){


print(paste("i=",i))
	

### Do random cluster size
dat<-make.pois.data(K,10,10,0,rcluster=T)
			Y<-dat[,"Y"]
		strata<-dat[,"strata"]
		treat<-dat[,"treat"]		
	mout<-	myags(Y~treat,id=strata,family=poisson,
	scalefun=SCALEFUN,corstr=CORR,maxiter=MAXITER)	
	bd1[i]<-mout$coefficients[2]
	vd1[i,]<-c(mout$naive.parmvar[2,2],mout$robust.parmvar[2,2],mout$robust.parmvar.md[2,2],NA)	
	if (is.null(mout$errorcode)){	
 sout<-saws2(mout$coefficients,mout$u,mout$omega,test=2)
   vd1[i,4]<-sout$arse1^2
   wd1[i]<-sout$ssew
   pd1[i,]<-sout$twosided.pvalues
   cid1[i,]<-sout$ci.length
   add<-rep(0,Nt)
   add[pd1[i,]<=0.05]<-1
   nrejectd1<-nrejectd1+ add
   converge.countd1<-converge.countd1+1
    }
else{ errors.d1<-c(errors.d1,mout$errorcode) 
	   erri.d1<-c(erri.d1,i)
	}
print(round(nrejectd1/converge.countd1,digits=4))


## with random effects on treatments
dat<-make.pois.data(K,10,10,.5,rcluster=T)

		Y<-dat[,"Y"]
		strata<-dat[,"strata"]
		treat<-dat[,"treat"]		
	mout<-	myags(Y~treat,id=strata,family=poisson,
		scalefun=SCALEFUN,corstr=CORR,maxiter=MAXITER)	
	bd4[i]<-mout$coefficients[2]

	vd4[i,]<-c(mout$naive.parmvar[2,2],mout$robust.parmvar[2,2],mout$robust.parmvar.md[2,2],NA)
			
	if (is.null(mout$errorcode)){	

 sout<-saws2(mout$coefficients,mout$u,mout$omega,test=2)
   vd4[i,4]<-sout$arse1^2
   wd4[i]<-sout$ssew
		
   pd4[i,]<-sout$twosided.pvalues
   cid4[i,]<-sout$ci.length
   add<-rep(0,Nt)
   add[pd4[i,]<=0.05]<-1
   nrejectd4<-nrejectd4+ add
   converge.countd4<-converge.countd4+1
}
else{ errors.d4<-c(errors.d4,mout$errorcode) 
	   erri.d4<-c(erri.d4,i)
	}

print(round(nrejectd4/converge.countd4,digits=4))




}

out<-list(pd1=pd1,cid1=cid1,vd1=vd1,wd1=wd1,bd1=bd1,
converge.countd1=converge.countd1,errors.d1=errors.d1,
erri.d1=erri.d1,
pd4=pd4,cid4=cid4,vd4=vd4,wd4=wd4,bd4=bd4,
converge.countd4=converge.countd4,errors.d4=errors.d4,
erri.d4=erri.d4)

#dump("out","sim.Pois.exchange.out")



## Now print out results

print.sim<-function(outp,outci,outerri){

test.names<-c("Vm,Chi2","Vs,Chi2","Vs,F,K-p","Vs,F,d.tilde.no.H","Vs,F,d.hat.no.H",
"Va,Chi2","Va,F,d.tilde.with.H","Va,F,d.hat.with.H")

nsim<-dim(outp)[[1]]
index<-rep(T,nsim)
if (!is.null(outerri)) index[outerri]<-F

print(paste("Number Converged=",length((1:nsim)[index]) ) )

prop.le.alpha<-function(x,alpha=.05) 
 length(x[x<=alpha])/length(x)

print("P-Values")
pval<-apply(outp[index,],2,prop.le.alpha)
names(pval)<-test.names
print(pval)
print("Mean ci")
meanci<-apply(outci[index,],2,mean)
names(meanci)<-test.names
print(meanci)
invisible()
}



quantiles<-function(x){
	out<-quantile(x)
	out<-c(out,mean(x))
	names(out)[6]<-"mean"
	out
}


quantiles(out$wd1)
print.sim(out$pd1,out$cid1,out$erri.d1)

mult<-10^0
quantiles(mult*out$vd1[,1])
quantiles(mult*out$vd1[,2])
quantiles(mult*out$vd1[,3])
quantiles(mult*out$vd1[,4])
mult*var(out$bd1)

quantiles(out$wd4)
print.sim(out$pd4,out$cid4,out$erri.d4)

quantiles(mult*out$vd4[,1])
quantiles(mult*out$vd4[,2])
quantiles(mult*out$vd4[,3])
quantiles(mult*out$vd4[,4])
mult*var(out$bd4)


out.sim.pois.indep.nofix<-out
dump("out.sim.pois.indep.nofix","out.sim.pois.indep.nofix")


proc.time()-now