### Code for the correlation plot in Zheng et al (2006)
### The "How many prisoners do you know" paper. 
### [IMPORTANT] You may need to re-size your graphics window 
### after the plot is done to get better effects. 
###############
### Step 0: get plot data
attach(sims.all, name='sims.obj')

  n.keep<-n.iter/n.thin

   result.alpha<-array(NA, c(m.iter, data.n, n.keep/2))
   result.beta<-array(NA, c(m.iter, data.j, n.keep/2))
   result.omega<-array(NA, c(m.iter, data.j, n.keep/2))
   for (i in 1:data.n) result.alpha[,i,] <- t(sims[-(1:(n.keep/2)),,i])
   for (i in 1:data.j) {result.beta[,i,] <- t(sims[-(1:(n.keep/2)),,data.n+i])
                result.omega[,i,] <- t(sims[-(1:(n.keep/2)),,data.n+data.j+i])}

y.pred<-exp(outer(result.alpha[1,,105], result.beta[1,,105], '+'))

y.diff<-sqrt(y)-sqrt(y.pred)
y.std<-y.diff
for(j in 1:data.j){
      y.std[,j]<-(y.diff[,j]-mean(y.diff[,j], na.rm=T))/sqrt(var(y.diff[,j], use='complete'))      
}

### Step 1: prepare correlation matrix
z.plot<-cor(sqrt(y)-sqrt(y.pred), use='pairwise.complete')
### Step 2: calculate order of the rows/columns
###         it can be specified manually or computed using hierarchical clustering.
z.order<-c(11,9,7,5,3,1,12,10,8,6,4,2) ## first names
#z.other<-z.other[hclust(dist(t(y.std[,z.other])))$order]
z.other<-c(32,31,28,27,30,29,25,26,13,20,24,18,16,21,17,19,23,22,14,15) ## other category
z.order<-c(z.other, z.order)
### Step 3: order the correlation matrix according to z.order
z.plot<-z.plot[z.order, z.order]

### Step 4: mask out the upper half
for(i in 1:32)
 for(j in i:32)
       z.plot[i,j]<-0 ## in this project, we do not pay attention to negative correlation

### Step 5: get the range of the correlation values
max.cor<-0.7
min.cor<-min(z.plot)

### Step 6: graphics layout
layout(matrix(c(1,2), 1, 2, byrow=TRUE), c(10,1))
par(mar=c(0.1, 0.1, 2, 0.1), pty='s')
###         break values for image function
###         here the number of levels is set to 8
z.breaks<-c(min(z.plot), seq(0, 0.7, 0.1))
z.colors<-gray(8:1/8)
image(1:32, 1:32, z.plot, xaxt='n', yaxt='n', col=z.colors,
             breaks=z.breaks, 
             xlim=c(-6, 33.5), ylim=c(-6, 33.5), xaxt='n', yaxt='n', bty='n', xlab='', ylab='')

### Step 7: plot legends
###         category names
text(1:32, 1:32, category.short[z.order], adj=1, cex=0.9)
###         boxes
rect(26.5, 26.5, 32.5, 32.5, lty=2)
text(26, 32.5, 'female \n names', cex=1, adj=1, font=2)
rect(20.5, 20.5, 26.5, 26.5, lty=2)
text(20, 26.5, 'male \n names', cex=1, adj=1, font=2)
rect(8.5, 8.5, 0.5, 0.5, lty=2)
text(0, 8.5, 'negative \n experiences', adj=1, cex=1, font=2)

###         the color legend
par(mar=c(4, 0.1, 2 ,0.1), pty='m')
plot(c(0.2,1), c(min(z.breaks), max(z.breaks)), type='n', bty='n', xlab='', ylab='', xaxt='n', yaxt='n')
options(digits=1)
for(i in 2:(length(z.breaks))){
      rect(0.5, z.breaks[i-1], 1, z.breaks[i], col=z.colors[i-1])
      text(0.45, z.breaks[i-1], format(round(z.breaks[i-1],1)), cex=1, adj=1)
}
rect(0.5, z.breaks[length(z.breaks)], 1, z.breaks[length(z.breaks)], col=z.colors[length(z.colors)])
detach('sims.obj')
#######################################