######################################################
## A sample code to run CARw in Lu et al. (2005)   ###
## By Haolan Lu, 8-3-2005                          ###
######################################################

##################################################
### read a dataset that is generated based    ####
### on Figure 1 template in Lu et al. (2005)  ####
### set up adjacent information , Z           ####
##################################################



true.geo.num.ngb<-c(2,3,3,2,3,4,4,3,3,4,4,3,2,3,3,2)
true.geo.adj<-c(2,5,1,3,6,2,4,7,3,8,1,6,9,2,5,7,10,3,6,8,11,4,7,12,5,10,13,
6,9,11,14,7,10,12,15,8,11,16,9,14,10,13,15,11,14,16,12,15)

geo.ind1<-geo.ind2<-rep(0, 16)
for (i in 1:16){
    if(i==1){geo.ind1[i]<-1; geo.ind2[i]<-true.geo.num.ngb[i]}
    else{
        geo.ind2[i]<-sum(true.geo.num.ngb[1:i])
        geo.ind1[i]<-geo.ind2[i-1]+1
        }
}

y<-c(19,13,97,87,10,11,103,91,11,690,103,94,759,775,695,709)

Z<-matrix(0, ncol=16, nrow=16)

for( i in 1:numcty){
    for(j in geo.ind1[i]:geo.ind2[i]){
        Z[i,true.geo.adj[j]]<-5
        }
}

Z[2,3]<-Z[3,2]<-Z[6,7]<-Z[7,6]<-Z[6,10]<-Z[10,6]<-Z[9,10]<-Z[10,9]<-10
Z[9,13]<-Z[13,9]<-Z[10,11]<-Z[11,10]<-Z[11,15]<-Z[15,11]<-Z[12,16]<-Z[16,12]<-10

byrow.matrix<-function(x,nr, nc){
    new<-matrix(0, nrow=nr, ncol=nc)
    for(i in 1:nr){
        a<-(i-1)*nc+1
        b<-i*nc
        new[i,]<-x[a:b]
        }
    return(new)
}

##############################################################
#### set up the initial values for calling method 3     ######
##############################################################

numcty<-16
Zarray<-byrow.matrix(Z, nr=1, nc=16*16)
E<-rep(100,16)


X1<-rep(0,16)
area1<-c(1,2,5,6,9)
area2<-c(3,4,7,8,11,12)
area3<-c(10,13,14,15,16)
X1[area1]<--2
X1[area2]<-0
X1[area3]<-2

Xarray<-X1

epsilon<-0.5
tau<-10
tauPriorAlpha<-1
tauPriorBeta<-1
Warray<-matrix(true.geo.W,nrow=1, byrow=T)
Dwarray.eps<-matrix(diag(true.geo.num.ngb+epsilon), nrow=1, byrow=T)
Dwarray<-matrix(diag(true.geo.num.ngb), nrow=1, byrow=T)

phiAccept<-rep(0,numcty)
gammaAccept<-0

memory.limit(size=9999999999)

nupdate<-11000
burnin<-1000

dimBeta<-1
betaPriorVar<-rep(1,dimBeta)
betaPriorMean<-c(1,dimBeta)

betaUpdated<-rep(0, dimBeta*nupdate)
betaAccept<-rep(0, dimBeta)
phiUpdated<-rep(0,numcty*nupdate)

tauUpdated<-rep(0, nupdate)
tempphi<-rep(0,nupdate)
phiVect<-true.theta
avgPhiVect<-true.theta%*%true.geo.W/true.geo.num.ngb

phiScaleVect[area3]<-0.1
phiScaleVect[area1]<-0.6
phiScaleVect[area2]<-0.3

beta<-1
betaScale<-rep(0.05, dimBeta)
phiScaleVect[4]<-0.35
gammaPriormean<-c(6,-1)
gammaPriorVar<-c(1,1)
gscaleVect<-c(0.7, 0.7)
beta0<-0
gamma<-c(1,-4)
lengthwij<-24
What<-rep(0, lengthwij)

method3<-function(){
r<-.C("carlin", as.double(beta),as.integer(dimBeta), as.double(betaScale), as.double(betaPriorVar),as.double(betaPriorMean),as.double(y), as.double(E),
            as.double(Xarray), as.double(newZarray), as.integer(geo.ind1), as.integer(geo.ind2),as.integer(true.geo.adj),
            as.integer(true.geo.num.ngb), as.double(phiVect), as.double(avgPhiVect), as.double(phiScaleVect),
            as.double(tau),as.double(tauPriorAlpha), as.double(tauPriorBeta), as.double(Warray), as.double(Dwarray.eps),
            as.double(gamma), as.double(gscaleVect),  as.double(gammaPriorVar), as.double(gammaPriormean),as.integer(nupdate),
            phiUpdated=as.double(phiUpdated),
            phiAccept=as.double(phiAccept), tauUpdated=as.double(tauUpdated),
            gamma0=as.double(tempphi),gamma1=as.double(tempphi),betaUpdated=as.double(betaUpdated), betaAccept=as.double(betaAccept),
            gammaAccept=as.double(gammaAccept),  epsilon=as.double(epsilon),
            numcty=as.integer(numcty), What=as.double(What), lengthwij=as.integer(lengthwij), burnin=as.integer(burnin)
            )
    phi<-matrix(r$phiUpdated, nrow=nupdate, ncol=numcty, byrow=T)
    par(mfrow=c(3,2),oma=c(0.5,0.5,3,0.5))
    ts.plot(r$tauUpdated[(burnin+1):nupdate])
    ts.plot(r$betaUpdated[(burnin+1):nupdate])
    ts.plot(r$gamma0[(burnin+1):nupdate])
    ts.plot(r$gamma1[(burnin+1):nupdate])
    ts.plot(phi[((burnin+1):nupdate),4]+r$betaUpdated[(burnin+1):nupdate]*X1[4])
    ts.plot(phi[((burnin+1):nupdate),14]+r$betaUpdated[(burnin+1):nupdate]*X1[14])
     mtext(i, side=3, outer=T)
}

###########################################################
### Note: you need download carlin-beta-wombling.dll    ###
### and save it under "c:/" directory.                  ###
### in order to call .C("carlin") function from R.      ###
### Please visit r-project.org for the details about    ###
### the set-up for calling a C function in R.           ###
###########################################################

dyn.load("c:\\carlin-beta-wombling.dll")
method3()
dyn.load("c:\\carlin-beta-wombling.dll")