library(ergm)
library(network)
# setwd( "D:/Thomas/UoW/Topics/ERGM/R/Marginal Fitting")
# setwd("/home/tsuesse/R/ERGM/Marginal")
# source("Cov_struc_ergm.r")
rm(list=ls())
dir<-F
use.Psi<-T  # TRUE


# source("fcts.Lazega.r")




############################################################################################################
# DON'T RUN
if(0){
  
#statistic
# stat=1: ESP (edgewise shared partners)
# stat=2: triangle
# stat=3: degree
# stat=4: GWESP 
stat<-4


N<-6 #number nodes


# BEGIN SIMULATION

switch(stat,
{
#stat=1 edgewise-shared partner
form1<-as.formula(net1 ~ edges + esp(3))
form1.stat<-as.formula(net1~esp(0:(N-2)))
}
,
# stat=2: triangle
{
form1<-as.formula(net1 ~ edges + triangle)
form1.stat<-as.formula(net1~triangle)
}
,
# stat=3: degree
{
form1<-as.formula(net1 ~ edges + degree(3))  # undirected only
form1.stat<-as.formula(net1~degree(0:N))
}
,
#g.sim1<-simulate(network(N,density=expit(edges0),directed=F) ~ edges+gwesp(fixed = T), coef=c(edges0,gwesp0))
form1<-as.formula(net1 ~ edges+gwesp(fixed = T))
)

net1<-network(N,density=0.1,directed=T)
#net1<-network(N)
#net1<-simulate.ergm(net1~ edges + mutual,nsim=1, coef=c(0.1,0.5))
#net1<-simulate(network(16) ~ edges + mutual,coef=c(0,0))
net1<-simulate(form1,nsim=1, coef=c(0.1,0.5))

summary(form1.stat)






#calculate exactly
library(sfsmisc)
#now simulate several networks and record correlation


sim<-1e5
Y<-array(0,dim=c(sim,N,N))
Y.vec<-matrix(0,sim,N*(N-1)/2)

names.Y.vec<-NULL
for(i  in 1:(N-1)){
names.Y.vec<-c(names.Y.vec,paste((i+1):N,i,sep=""))
}
colnames(Y.vec)<-names.Y.vec


ind.lt<-lower.tri(diag(N))



net2<-simulate(form1,nsim=sim, coef=c(0.1,0.5),control=control.simulate(MCMC.burnin=1e4, MCMC.interval=1e4))

for(i in 1:sim){

if(i%%100==0){cat("iteration at i=",i,"\n")}

 hilf<-as.sociomatrix(net2[[i]])
 Y[i,,]<-hilf
 Y.vec[i,]<-c(hilf[ind.lt])
}# end for loop
 
rm(net2)

cor(Y.vec)
 
 
 break
 break
 
}# end if(0)
 # END SIMULATION


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

##############################################################################################
 # exact probabilities without simulation

library(sfsmisc)


if(dir){N<-5}else{N<-6}

names.Y.vec<-NULL
if(dir){
for(i  in 1:(N-1)){
  names.Y.vec<-c(names.Y.vec,paste((i+1):N,i,sep=""))
}
}
for(i  in 1:(N-1)){
  names.Y.vec<-c(names.Y.vec,paste(i,(i+1):N,sep=""))
}


ind.lt<-lower.tri(diag(N))

theta1<- -0.1
theta2<-0.5
theta3<-2*theta2


dyads<-N*(N-1)/(2-(dir)*1)
sim<-2^dyads

Y.vec1<-matrix(0,sim,dyads)
colnames(Y.vec1)<-names.Y.vec

stats.matrix<-matrix(0,sim,6) # 6 statistics



if(use.Psi){
  
Psi.nodes<- 1:N
  
#Psi<-Psi.matrix(runif(dyads),N) 
Psi <- matrix(Psi.nodes,N,N,byrow=F) + matrix(Psi.nodes,N,N,byrow=T) 
diag(Psi)<-0
#Psi<- Psi - mean(Psi)*N^2/(N*(N-1))
#diag(Psi)<-0
#Psi<-Psi/N
Psi<- (Psi - 7)/5 - 1 
}

 # binary
for(i in 1:sim){
 
 Y.vec1[i,]<-empty.dimnames(digitsBase(i-1,2,dyads))
 Ynet<-matrix(0,N,N)
 
 if(!dir){
 Ynet[ind.lt]<-Y.vec1[i,]
 Ynet<-Ynet+t(Ynet)
 }else{
 Ynet[ind.lt]<-Y.vec1[i,1:(dyads/2)]
 Ynet<-t(Ynet)
 Ynet[ind.lt]<-Y.vec1[i,((dyads/2)+1):dyads]
 
 }
 
 
 var(Psi[ind.lt])  # 0.2
 mean(Psi[ind.lt]) # -1
 
 
 if(i%%100==0){cat("iteration at i=",i,"\n");print(Ynet)}

 net1<-as.network(Ynet,directed=dir)
 
 stats.matrix[i,1]<-exp(summary(net1~edges)*theta1+ summary(net1~esp(3))*theta2)  # edges + esp(3)
 stats.matrix[i,2]<-exp(summary(net1~edges)*theta2)                              # just edges
 stats.matrix[i,3]<-exp(summary(net1~gwesp(fixed=T))*theta2)                        # just gwesp
 stats.matrix[i,4]<-exp(summary(net1~edges)*theta1+ summary(net1~triangle)*theta2)  # edges + triangle

  if(!dir){
  stats.matrix[i,5]<-exp(summary(net1~edges)*theta1+ summary(net1~degree(3))*theta2)   # edges+ degree (undirected)
  }else{
  stats.matrix[i,5]<-exp(summary(net1~edges)*theta1+ summary(net1~idegree(3))*theta2 + summary(net1~odegree(3))*theta3)   # edges+ idegree (directed)
  }
  stats.matrix[i,6]<- exp(summary(net1~edges)*theta1 + summary(net1~gwesp(fixed = T))*theta2) # edges + gwesp
 
  if(use.Psi){
    # adjust the weighting caharacterised by network statistics with exp(Psi*Y)
  stats.matrix[i,]<-stats.matrix[i,]*exp(sum(Psi*Ynet))
  }
  
}#end for


break
break

# stats<-stats.degree3   #stats.degree3<-stats.triangle<-stats.gwesp<-stats.edges.esp3<-stats.edges
# stats<-stats.triangle
# stats<-stats.gwesp
# stats<-stats.esp3
# stats<-stats.edges
# stats<-stats.gwesp
library(tensor)

probs<-stats.matrix/matrix(colSums(stats.matrix),sim,6,byrow=TRUE)

#compute correlations
EY<-Corr<-list()
for(j in 1:6){
EY[[j]]<-t(Y.vec1)%*%probs[,j]
EYYt<-t(Y.vec1)%*%(probs[,j]*Y.vec1)
Cov <-EYYt- EY[[j]]%*%t(EY[[j]])
Corr[[j]]<- diag(sqrt(1/diag(Cov)))%*%Cov%*%diag(sqrt(1/diag(Cov)))
cat("Var equal to Bernoulli???\n")
print(all.equal(diag(Cov),c(EY[[j]]*(1-EY[[j]]))))
}
Corr



save(Corr,EY,file=paste("Corr_dir",dir,"_Psi",use.Psi,"new.RData",sep=""))
##############################################################################################################
#                          END SIMULATION EXACT
##############################################################################################################

##############################################################################################################
#                               BEGIN PLOTS
##############################################################################################################
use.Psi<-F
dir<-T
load(paste("Corr_dir",dir,"_Psi",use.Psi,".RData",sep=""))

source("fcts.Lazega.r")
library(graphics)
#library(scatterplot3d)
#library(lattice)
library(fields)

if(!dir){
####################################################



# undirected


for(k in c(1,2,3,4,5,6)){
  # k=1: edges +esp(3)
  # k=2 : edges
  # k=3: gwesp
  # k=4: edges+ triangle
  # k=5: degree
  # k=6: edges+gwesp
  
  
  titles<-c("edges + esp(3)","edges","gwesp","edges + triangle","edges + degree","edges + gwesp")
  
  
  Corr1<-Corr[[k]]
  Corr1<-abs(round(Corr1,3))
  
  
  ind<-indexfct(N)
  names<-rep(NA,dyads)
  for(i in 1:dyads){
    names[i]<-paste("(",ind[1,i],",",ind[2,i],")",sep="")
  }
  
  
  breaks.ini<-as.numeric(names(table(c(Corr1))))
  #breaks<-c(breaks[breaks<0],-1e-4,0,1e-4,breaks[breaks>0])
  #breaks<-c(1e-4,breaks)-1e-5
  if(breaks.ini[1]==0){
    breaks<-c(0,breaks.ini[-1]-1e-5,1)
    breaks1<-c(0,breaks.ini[-1],1)
  }else{
    breaks<- c(0,breaks.ini-1e-5,1)
    breaks1<- c(0,breaks.ini,1)
  }
  
    
  
  
  
  l.breaks<-length(breaks)
  
  cols<-grey.colors(l.breaks,start=0,end=1)
  cols<-cols[-(l.breaks-1)]
  
  
  #(100-l.breaks+2):100
  jpeg(filename = paste("Corr",k,"_undir.jpeg",sep=""),width = 3000, height =3000,  res = 300)
  nf<-layout(matrix(c(1,2),1,2,byrow=T),widths =c(6,1) ,heights=c(4))
  par(mar=c(3,3,3,2))
  #par(mar=c(2,2,2,1))
  image(x=1:dyads,y=1:dyads,Corr1,lwd=1,col=cols,breaks=breaks,axes=FALSE,xlab="dyads",ylab="dyads",main=titles[k],cex.main=2)
  axis(1,at=1:dyads,labels=names)
  axis(2,at=1:dyads,labels=names)
  
  if(use.Psi){
  list.arg<-list( at=0.5+c(1,ceiling(l.breaks/2),l.breaks-3,l.breaks-1),
                 labels=breaks1[c(1,ceiling(l.breaks/2),l.breaks-2,l.breaks-1)] )
  }else{
  list.arg<-list( at=seq(1.5,l.breaks,by=1) , labels=breaks1[-(l.breaks-1)] )
  }
  
  
  par(mar=c(0,0,0,0))
  plot.new()
  image.plot(zlim=range(0,Corr1),legend.only=TRUE,col=cols,smallplot=c(.1,.4, .1,.9),
             cex.main=0.8,cex.axis=0.8,cex.lab=0.8,cex.sub=0.8,breaks=1:l.breaks
             #,axis.args=list(at=c(0.02,0.1,1), labels=c("0.02","0.1","1"))  
             #,axis.args=, labels=breaks1[-(l.breaks-1)] ) 
             ,axis.args=list.arg 
             ) 
  dev.off()
  
  
}#end for k



}else{


# directed
ind<-indexfct(N)



names<-rep(NA,dyads)
for(i in 1:(dyads/2)){
  names[i]<-paste("",ind[1,i],",",ind[2,i],"",sep="")
  names[i+(dyads/2)]<-paste("",ind[2,i],",",ind[1,i],"",sep="")
}









k<-1 # any number between 1 and 6
for(k in c(1,2,4,5,6)){
# k=1: edges +esp(3)
# k=2 : edges
# k=3: gwesp
# k=4: edges+ triangle
# k=5: (i)degree
# k=6: edges+gwesp


titles<-c("edges + esp(3)","edges","gwesp","edges + triangle","edges + idegree + odegree","edges + gwesp")


Corr1<-Corr[[k]]
Corr1<-abs(round(Corr1,3))

breaks.ini<-as.numeric(names(table(c(Corr1))))
#breaks<-c(breaks[breaks<0],-1e-4,0,1e-4,breaks[breaks>0])
#breaks<-c(1e-4,breaks)-1e-5
if(breaks.ini[1]==0){
  breaks<-c(0,breaks.ini[-1]-1e-5,1)
  breaks1<-c(0,breaks.ini[-1],1)
}else{
  breaks<- c(0,breaks.ini-1e-5,1)
  breaks1<- c(0,breaks.ini,1)
  }



l.breaks<-length(breaks)

cols<-grey.colors(l.breaks,start=0,end=1)
cols<-cols[-(l.breaks-1)]


#(100-l.breaks+2):100
jpeg(filename = paste("Corr",k,"_dir.jpeg",sep=""),width = 3000, height =3000,  res = 300)
nf<-layout(matrix(c(1,2),1,2,byrow=T),widths =c(6,1) ,heights=c(4))
par(mar=c(3,3,3,2))
#par(mar=c(2,2,2,1))
image(x=1:dyads,y=1:dyads,Corr1,lwd=1,col=cols,breaks=breaks,axes=FALSE,xlab="dyads",ylab="dyads",main=titles[k],cex.main=2)
axis(1,at=1:dyads,labels=names)
axis(2,at=1:dyads,labels=names)




par(mar=c(0,0,0,0))
plot.new()
image.plot(zlim=range(0,Corr1),legend.only=TRUE,col=cols,smallplot=c(.1,.4, .1,.9),
           cex.main=0.8,cex.axis=0.8,cex.lab=0.8,cex.sub=0.8,breaks=1:l.breaks,
           #axis.args=list(at=c(0.02,0.1,1), labels=c("0.02","0.1","1")) ) 
           axis.args=list(at=seq(1.5,l.breaks,by=1), labels=breaks1[-(l.breaks-1)] ) ) 
dev.off()


}#end for k

}#end if dir

################################################################################
mean.var.w<-function(Y,w,binary=FALSE){
MCMCsamplesize<-length(w)
eps<-10*(1/MCMCsamplesize);
N<-length(w)
m.Y<- t(matrix(1,1,N)%*%(w*Y))
m.Y11<-t(Y)%*%(w*Y)
if(binary){m.Y[m.Y<eps]<-eps;m.Y[m.Y>(1-eps)]<-1-eps;diag(m.Y11)<-m.Y
if(sum(m.Y<eps | m.Y>(1-eps))>0){cat("binary mean modified!!!")}
}
cov.Y <- m.Y11 - m.Y%*%t(m.Y)
return(list(m.Y=m.Y,m.Y11=m.Y11,cov.Y=cov.Y,v.Y=diag(cov.Y)))
}#end  mean.var.w<-function(Y,w,binary=FALSE){
################################################################################

hilf<- mean.var.w(Y.vec,probs,binary=TRUE)

indexfct<-function(c){
  LODDS<-c*(c-1)/2;
  hilf<-NULL;
  for(i in 1:(c-1)){
    for(j in (i+1):c){
      hilf<-cbind(hilf,matrix(c(i,j),2,1))
    }
  }
  hilf
}#end indexfct