require(lattice)
require(plot3D)
phi = seq(.1,2.1,.1)
gama = seq(0,2,.1)

#########################
# EQM(a)
mEQMa = read.table("mEQMa.txt",header=FALSE)
EQMa = as.matrix(mEQMa[,2:22])
head(EQMa)

hist3D (gama, phi, EQMa, phi = 40, theta = -50,
        col = "grey", NAcol = "white", breaks = NULL,
        border = 1, facets = TRUE, colkey = NULL,
        image = FALSE, contour = FALSE, ticktype = "detailed",
        xlab="",ylab="",zlab="",bty = "g",shade=.9,zlim=c(0,.1))
text3D(1,-0.2,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.4,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.7,2.4,0.07,labels=expression(paste("EQM(",hat(alpha),paste(")"))),add=TRUE,label="x",phi = 40, theta = 40)


#########################
# EQM(b)

mEQMb = read.table("mEQMb.txt",header=FALSE)
EQMb = as.matrix(mEQMb[,2:22])
head(EQMb)

hist3D(gama, phi, EQMb, phi = 40, theta = -50,
        col = "grey", NAcol = "white", breaks = NULL,
        border = 1, facets = TRUE, colkey = NULL,
        image = FALSE, contour = FALSE, ticktype = "detailed",
        xlab="",ylab="",zlab="",bty = "g",shade=.9,zlim=c(0,.1))
text3D(1,-0.2,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.4,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.7,2.4,0.07,labels=expression(paste("EQM(",hat(beta),paste(")"))),add=TRUE,label="x",phi = 40, theta = 40)


#########################
# EQM(g)

mEQMg = read.table("mEQMg.txt",header=FALSE)
EQMg = as.matrix(mEQMg[,2:22])
head(EQMg)

hist3D(gama, phi, EQMg, phi = 40, theta = -50,
       col = "grey", NAcol = "white", breaks = NULL,
       border = 1, facets = TRUE, colkey = NULL,
       image = FALSE, contour = FALSE, ticktype = "detailed",
       xlab="",ylab="",zlab="",bty = "g",shade=.9)
text3D(1,-0.2,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.4,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.7,2.4,3,labels=expression(paste("EQM(",hat(gamma),paste(")"))),add=TRUE,label="x",phi = 40, theta = 40)

#########################
# Número de não convergências

mNF = read.table("mNF.txt",header=FALSE)
NF = as.matrix(mNF[,2:22])
head(NF)

hist3D(gama, phi, NF, phi = 40, theta = -50,
       col = "grey", NAcol = "white", breaks = NULL,
       border = 1, facets = TRUE, colkey = NULL,
       image = FALSE, contour = FALSE, ticktype = "detailed",
       xlab="",ylab="",zlab="",bty = "g",shade=.9)
text3D(1,-0.2,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.4,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.5,2.2,600,labels=expression(paste("nf")),add=TRUE,label="x",phi = 40, theta = 40)

#########################
# gráficos do EQM de g e nf para certos valores de phi

par(mfrow=c(1,3))
plot(gama,NF[,5],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab=NA,xlab=expression(gamma))
text(.3,1100,label=expression(paste(phi,"=0.5")))
par(new=TRUE)
plot(gama,EQMg[,5],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2)
axis(side=4)

plot(gama,NF[,10],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab=NA,xlab=expression(gamma))
text(.3,1100,label=expression(paste(phi,"=1.0")))
par(new=TRUE)
plot(gama,EQMg[,10],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2)
axis(side=4)

plot(gama,NF[,15],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab=NA,xlab=expression(gamma))
text(.3,1100,label=expression(paste(phi,"=1.5")))
par(new=TRUE)
plot(gama,EQMg[,15],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2)
axis(side=4)

##################################################
##################################################

#########################
# Preparando os gŕaficos de EQM e nf
postscript(file="EQMs_nf_phi.eps",horizontal = FALSE, width = 10.0, height = 10.0,
           paper="special")
par(mfrow=c(2,2))
# EQM(a)
hist3D (gama, phi, EQMa, phi = 40, theta = -50,
        col = "grey", NAcol = "white", breaks = NULL,
        border = 1, facets = TRUE, colkey = NULL,
        image = FALSE, contour = FALSE, ticktype = "detailed",
        xlab="",ylab="",zlab="",bty = "g",shade=.9,zlim=c(0,.1),
        sub="(a)")
text3D(1,-0.4,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.6,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40,theta = 40)
text3D(-0.5,1.9,.05,labels=expression(paste("EQM(",hat(alpha),")")),add=TRUE,label="x",srt=110)
# EQM(b)
hist3D(gama, phi, EQMb, phi = 40, theta = -50,
       col = "grey", NAcol = "white", breaks = NULL,
       border = 1, facets = TRUE, colkey = NULL,
       image = FALSE, contour = FALSE, ticktype = "detailed",
       xlab="",ylab="",zlab="",bty = "g",shade=.9,zlim=c(0,.1),sub="(b)")
text3D(1,-0.4,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.6,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40,theta = 40)
text3D(-0.5,1.9,.05,labels=expression(paste("EQM(",hat(beta),")")),add=TRUE,label="x",srt=110)
# EQM(g)
hist3D(gama, phi, EQMg, phi = 40, theta = -50,
       col = "grey", NAcol = "white", breaks = NULL,
       border = 1, facets = TRUE, colkey = NULL,
       image = FALSE, contour = FALSE, ticktype = "detailed",
       xlab="",ylab="",zlab="",bty = "g",shade=.9,sub="(c)")
text3D(1,-0.4,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.6,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40,theta = 40)
text3D(-0.5,1.9,2.5,labels=expression(paste("EQM(",hat(gamma),")")),add=TRUE,label="x",srt=110)
# Número de não convergências
hist3D(gama, phi, NF, phi = 40, theta = -50,
       col = "grey", NAcol = "white", breaks = NULL,
       border = 1, facets = TRUE, colkey = NULL,
       image = FALSE, contour = FALSE, ticktype = "detailed",
       xlab="",ylab="",zlab="",bty = "g",shade=.9,sub="(d)")
text3D(1,-0.4,0,labels=expression(gamma),add=TRUE,label="x",phi = 40, theta = 40)
text3D(-0.6,1.3,0,labels=expression(phi),add=TRUE,label="x",phi = 40,theta = 40)
text3D(-0.6,2.0,800,labels=expression(paste("nf")),add=TRUE,label="x",phi = 40, theta = 40)
dev.off()

#########################
# Preparando os gŕaficos do EQM de gama e nf para alguns valores de phi

postscript(file="EQMg_nf_phi.eps",horizontal = FALSE, width = 10.0, height = 6.0,
           paper="special")
par(mfrow=c(1,3))
plot(gama,NF[,5],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab="",xlab=expression(gamma),sub="(a)",lwd=2)
text(.3,1200,label=expression(paste(phi,"=0.5")),cex=1.5)
legend(0,1150,lty=c(1,2),c("nf",expression(paste("EQM(",hat(gamma),")"))),bty='n',cex=1.5)
par(new=TRUE)
plot(gama,EQMg[,5],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2,lwd=2)
text(2.5,2.5,labels=expression(paste("EQM(",hat(gamma),")")),srt=90)
axis(side=4)

plot(gama,NF[,10],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab=NA,xlab=expression(gamma),sub="(b)",lwd=2)
text(.3,1200,label=expression(paste(phi,"=1.0")),cex=1.5)
legend(0,1150,lty=c(1,2),c("nf",expression(paste("EQM(",hat(gamma),")"))),bty='n',cex=1.5)
par(new=TRUE)
plot(gama,EQMg[,10],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2,lwd=2)
axis(side=4)

plot(gama,NF[,15],type='l',ylim=c(0,1300),panel.first=grid(lty=3,lwd=1.5),
     ylab=NA,xlab=expression(gamma),sub="(c)",lwd=2)
text(.3,1200,label=expression(paste(phi,"=1.5")),cex=1.5)
legend(0,1150,lty=c(1,2),c("nf",expression(paste("EQM(",hat(gamma),")"))),bty='n',cex=1.5)
par(new=TRUE)
plot(gama,EQMg[,15],type='l',ylim=c(0,5),axes=FALSE,xlab=NA,ylab=NA,lty=2,lwd=2)
axis(side=4)
dev.off()
