library(ergm)
# source("change.stat.r")
stat<-3 # 1,2,3
N<-6 #number nodes

#edge 1
i1<-1
i2<-2

#edge 2
j1<-3
j2<-4

#edge 3
k1<-5
k2<-6


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))
         form1.stat<-as.formula(net1~degree(0:N))
       }
       
)

net1<-network(N,density=0.1,directed=T)
net1<-simulate(form1,nsim=1, coef=c(0.1,0.5))

summary(form1.stat)




net1[i1,i2]
net1[i2,i1]

net1[j1,j2]
net1[j2,j1]

#net1.mat<-as.sociomatrix(net1)
#net1.mat[i1,i2]
#net1.mat[i2,i1]

net2<-net1

delete.edges(net1,intersect(get.edgeIDs(net1,i1),get.edgeIDs(net1,i2)))
delete.edges(net1,intersect(get.edgeIDs(net1,j1),get.edgeIDs(net1,j2)))

delete.edges(net2,intersect(get.edgeIDs(net2,i1),get.edgeIDs(net2,i2)))
delete.edges(net2,intersect(get.edgeIDs(net2,i1),get.edgeIDs(net2,j2)))

net1.11<-net1.10<-net1.01<-net1.00<-net1
net2.11<-net2.10<-net2.01<-net2.00<-net2
#now add 2 edges, calculate summary
#add edges for net1
add.edges(net1.11,i1,i2)
add.edges(net1.11,j1,j2)

add.edges(net1.10,i1,i2)

add.edges(net1.01,j1,j2)
#add edges for net2
add.edges(net2.11,i1,i2)
add.edges(net2.11,i1,j2)

add.edges(net2.10,i1,i2)

add.edges(net2.01,i1,j2)

switch(stat,
       {
         stat1.11<-summary(net1.11~esp(0:(N-2)))
         stat1.10<-summary(net1.10~esp(0:(N-2)))
         stat1.01<-summary(net1.01~esp(0:(N-2)))
         stat1.00<-summary(net1.00~esp(0:(N-2)))
         
         stat2.11<-summary(net2.11~esp(0:(N-2)))
         stat2.10<-summary(net2.10~esp(0:(N-2)))
         stat2.01<-summary(net2.01~esp(0:(N-2)))
         stat2.00<-summary(net2.00~esp(0:(N-2)))
       }
       ,
       {
         stat1.11<-summary(net1.11~triangle)
         stat1.10<-summary(net1.10~triangle)
         stat1.01<-summary(net1.01~triangle)
         stat1.00<-summary(net1.00~triangle)
         
         stat2.11<-summary(net2.11~triangle)
         stat2.10<-summary(net2.10~triangle)
         stat2.01<-summary(net2.01~triangle)
         stat2.00<-summary(net2.00~triangle)
       }
       ,
       {
         stat1.11<-summary(net1.11~degree(0:N))
         stat1.10<-summary(net1.10~degree(0:N))
         stat1.01<-summary(net1.01~degree(0:N))
         stat1.00<-summary(net1.00~degree(0:N))
         
         stat2.11<-summary(net2.11~degree(0:N))
         stat2.10<-summary(net2.10~degree(0:N))
         stat2.01<-summary(net2.01~degree(0:N))
         stat2.00<-summary(net2.00~degree(0:N))
       }
)


delta1<-stat1.11-stat1.10-stat1.01+stat1.00
delta2<-stat2.11-stat2.10-stat2.01+stat2.00
cat("Delta statistics for dyads having no nodes in common",delta1,"\n")
cat("Delta statistics for dyads having 1 node in common",delta2,"\n")