devtools::install_github("amvallone/estdaR")
library(estdaR)
library(rgdal)

# FOR CITIES ####

load("b.rdata")
load("Wd.rdata")

spatial.markov <- sp.mkv(b,Wd,fixed=FALSE)

p1 <- sp.mkv(b,Wd,fixed=FALSE)[[1]]
St <- do.call(rbind, lapply(1:5,function(x)st.st(spatial.markov[[1]][,,x])))

m<-rbind(p1[,,1],St[1,],p1[,,2],St[2,],p1[,,3],St[3,],p1[,,4],St[4,],p1[,,5],St[5,])
l.l <- c(rep("l(1) ",6),rep("l(2) ",6),rep("l(3) ",6),rep("l(4) ",6),rep("l(5) ",6))
rownames(m) <- paste(l.l, rep(c(paste("Q ",1:5,sep=""), "\u03C0"),5),sep=" ")
colnames(m) <- paste("Q",1:5,sep=" ")
col<-c("white","red")
m<-round(m,3)
data <- reshape2::melt(m)
names(data) <- c("y","x","Probabilities")
s <- rep(.2,150)
s[seq(6,150,6)]<-.7
g <- ggplot(data,aes(x,y,Probabilities,label=sprintf("%0.3f", round(Probabilities, digits = 3)))) +
  geom_tile(aes(fill=Probabilities),color="black",size=s)+theme_bw() +
  scale_fill_gradient(low=col[1], high=col[2]) +
  labs(x=NULL, y=NULL,title=NULL) +
  theme(axis.title=element_text(size=16,face="bold"),
        panel.border = element_rect(colour="white"),
        axis.text = element_text(size=14),
        axis.line = element_line(),
        legend.text =element_text(size=14),
        legend.title=element_text(size=.8*16,face="bold"),
        plot.title = element_text(size=2*16,hjust=0.5))+
  geom_text(size=6)+scale_x_discrete(position = "top")+
  scale_y_discrete(limits = rev(levels(data$y)))
g

# FOR LMAs ####

b2 <-read.delim("b2.txt") # import population data panel
FLMAs <- readOGR(".", "FLMAs") # Import FLMA map
Wd2_nb <- read.gal("Wd2.gal") # Import contiguity W matrix
# Wd2 <- nb2listw(Wd2_nb, glist = NULL, style = "W")

# Create an inverse distance matrix
coor <- coordinates(FLMAs)
dist <- nbdists(Wd2_nb, coor)
idw <- lapply(dist, function(x) 1/x)
Wd2 <- nb2listw(Wd2_nb, glist = idw, style = "W")

# Computation of SMC
spatial.markov <- sp.mkv(b2,Wd2,fixed=FALSE)

p1 <- sp.mkv(b2,Wd2,fixed=FALSE)[[1]]
St <- do.call(rbind, lapply(1:5,function(x)st.st(spatial.markov[[1]][,,x])))

m<-rbind(p1[,,1],St[1,],p1[,,2],St[2,],p1[,,3],St[3,],p1[,,4],St[4,],p1[,,5],St[5,])
l.l <- c(rep("l(1) ",6),rep("l(2) ",6),rep("l(3) ",6),rep("l(4) ",6),rep("l(5) ",6))
rownames(m) <- paste(l.l, rep(c(paste("Q ",1:5,sep=""), "\u03C0"),5),sep=" ")
colnames(m) <- paste("Q",1:5,sep=" ")
col<-c("white","red")
m<-round(m,3)
data <- reshape2::melt(m)
names(data) <- c("y","x","Probabilities")
s <- rep(.2,150)
s[seq(6,150,6)]<-.7

# Plot of SMC
g <- ggplot(data,aes(x,y,Probabilities,label=sprintf("%0.3f", round(Probabilities, digits = 3)))) +
  geom_tile(aes(fill=Probabilities),color="black",size=s)+theme_bw() +
  scale_fill_gradient(low=col[1], high=col[2]) +
  labs(x=NULL, y=NULL,title=NULL) +
  theme(axis.title=element_text(size=16,face="bold"),
        panel.border = element_rect(colour="white"),
        axis.text = element_text(size=14),
        axis.line = element_line(),
        legend.text =element_text(size=14),
        legend.title=element_text(size=.8*16,face="bold"),
        plot.title = element_text(size=2*16,hjust=0.5))+
  geom_text(size=6)+scale_x_discrete(position = "top")+
  scale_y_discrete(limits = rev(levels(data$y)))
g