## Process the results in parallel

library(snowfall)

run.once<-function(i){
newbind<-function(x,y){
	nx<-nrow(x)
	ny<-nrow(y)
	rmax<-max(nx,ny)
	cmax<-ncol(x)+ncol(y)
	temp<-array(NA, c(rmax, cmax))
	temp[1:nx, 1:ncol(x)]<-x
	temp[1:ny, (ncol(x)+1):cmax]<-y
	temp
}
	library(coda)
	temp.bind<-NULL
	nfile<-paste('histnorm',i,'.txt', sep='')
	tfile<-paste('histt',i,'.txt', sep='')
	if (file.exists(nfile) & file.exists(tfile)){
	ncoda<-read.table(nfile)
	tcoda<-read.table(tfile)
	ncoda<-ncoda[2:ncol(ncoda)]
	tcoda<-tcoda[2:ncol(tcoda)]
	ncoda<-as.mcmc(ncoda)
	tcoda<-as.mcmc(tcoda)
	norm.sum<-summary(ncoda)
	t.sum<-summary(tcoda)
	temp.norm<-cbind(norm.sum$statistics[,1:2], norm.sum$quantiles[,c(1,5)], geweke.diag(ncoda)$z,effectiveSize(ncoda) )
	temp.t<-cbind(t.sum$statistics[,1:2], t.sum$quantiles[,c(1,5)], geweke.diag(tcoda)$z,effectiveSize(tcoda) )
	temp.bind<-newbind(temp.norm, temp.t)

	}
	temp.bind
}

m<-20

sfInit(parallel=TRUE, cpus=m)
library(coda)

allres<-sfLapply(1:100, run.once)

dput(allres, 'nast-N500T5.txt')

sfStop()

take1<-function(x, id){
	n<-length(x)
	m<-nrow(x[[1]])
	res<-array(NA, c(n, m))
	for (i in 1:n){
		res[i, ]<-x[[i]][, id]
	}
	res
}

nsd<-take1(allres, 2)
tsd<-take1(allres, 8)

cbind(apply(nsd, 2, mean), apply(tsd, 2, mean))

(apply(tsd, 2, mean)/apply(nsd, 2, mean)-1)*100



library(snowfall)

run.once<-function(i){
newbind<-function(x,y){
	nx<-nrow(x)
	ny<-nrow(y)
	rmax<-max(nx,ny)
	cmax<-ncol(x)+ncol(y)
	temp<-array(NA, c(rmax, cmax))
	temp[1:nx, 1:ncol(x)]<-x
	temp[1:ny, (ncol(x)+1):cmax]<-y
	temp
}
	library(coda)
	temp.bind<-NULL
	nfile<-paste('R',i,'/histnorm.txt', sep='')
	tfile<-paste('R',i,'/histt.txt', sep='')
	if (file.exists(nfile) & file.exists(tfile)){
	ncoda<-read.table(nfile)
	tcoda<-read.table(tfile)
	ncoda<-ncoda[2:ncol(ncoda)]
	tcoda<-tcoda[2:ncol(tcoda)]
	ncoda<-as.mcmc(ncoda)
	tcoda<-as.mcmc(tcoda)
	norm.sum<-summary(ncoda)
	t.sum<-summary(tcoda)
	temp.norm<-cbind(norm.sum$statistics[,1:2], norm.sum$quantiles[,c(1,5)], geweke.diag(ncoda)$z,effectiveSize(ncoda) )
	temp.t<-cbind(t.sum$statistics[,1:2], t.sum$quantiles[,c(1,5)], geweke.diag(tcoda)$z,effectiveSize(tcoda) )
	temp.bind<-newbind(temp.norm, temp.t)

	}
	temp.bind
}

m<-20

sfInit(parallel=TRUE, cpus=m)
library(coda)

allres<-sfLapply(1:100, run.once)

dput(allres, 'nast-N100T4.txt')

sfStop()

take1<-function(x, id){
	n<-length(x)
	m<-nrow(x[[1]])
	res<-array(NA, c(n, m))
	for (i in 1:n){
		res[i, ]<-x[[i]][, id]
	}
	res
}

nsd<-take1(allres, 2)
tsd<-take1(allres, 8)

cbind(apply(nsd, 2, mean), apply(tsd, 2, mean))

(apply(tsd, 2, mean)/apply(nsd, 2, mean)-1)*100