###############################################################################
#
#               BEGIN FUNCTION get.para
#
###############################################################################
get.para<-function(formula,theta0="MPLE",MPLEonly=FALSE,MLestimate=TRUE,seed=NULL,
burnin=10000,MCMCsamplesize=10000,interval=100,maxit=3,constraints=~.,meanstats=NULL,
control=control.ergm(),verbose=FALSE){
###############################################################################

current.warn <- options()$warn
options(warn=0)
  
nw <- ergm.getnetwork(formula)
  
 
  
  

  #control.drop=TRUE
   model.initial <- ergm.getmodel(formula, nw, drop=FALSE, initialfit=TRUE)
   model.initial.drop <- ergm.getmodel(formula, nw, drop=TRUE, initialfit=TRUE)
   namesmatch <- match(model.initial$coef.names, model.initial.drop$coef.names)
   droppedterms <- rep(FALSE, length=length(model.initial$etamap$offsettheta))
   droppedterms[is.na(namesmatch)] <- TRUE
   model.initial$etamap$offsettheta[is.na(namesmatch)] <- TRUE
  
  
  MHproposal <- MHproposal(constraints, weights=control$prop.weights, control$prop.args, nw, model.initial)
  MHproposal.miss <- MHproposal("randomtoggleNonObserved", control$prop.args, nw, model.initial)

  Clist.initial <- ergm.Cprepare(nw, model.initial)
  Clist.miss.initial <- ergm.design(nw, model.initial, verbose=verbose)
  
  Clist.initial$meanstats=meanstats
  theta0copy <- theta0
  initialfit <- ergm.initialfit(theta0copy, MLestimate, Clist.initial,
                                Clist.miss.initial, model.initial,
                                MPLEtype=control$MPLEtype, 
                                initial.loglik=control$initial.loglik,
                                verbose=verbose, compressflag = control$compress, 
                                maxNumDyadTypes=control$maxNumDyadTypes,
                                force.MPLE=(ergm.independencemodel(model.initial)
                                            && constraints==(~.)))
  MCMCflag <- ((MLestimate && (!ergm.independencemodel(model.initial)
                               || !is.null(meanstats)
                               || constraints!=(~.)))
                || control$force.mcmc)
  #MCMCflag is true              
                
  
    theta0 <- initialfit$coef
    names(theta0) <- model.initial$coef.names
    theta0[is.na(theta0)] <- 0
  
  
  
  
   model <- ergm.getmodel(formula, nw, drop=FALSE, expanded=TRUE, silent=TRUE)
   # revise theta0 to reflect additional parameters
   theta0 <- ergm.revisetheta0(model, theta0)
   model.drop <- ergm.getmodel(formula, nw, drop=TRUE, expanded=TRUE, silent=TRUE)
#            silent="MPLE" %in% theta0copy)
   namesdrop <- model$coef.names[is.na(match(model$coef.names, model.drop$coef.names))]
   names(model$etamap$offsettheta) <- names(theta0)
   droppedterms <- rep(FALSE, length=length(model$etamap$offsettheta))
   droppedterms[is.na(namesmatch)] <- TRUE
   theta0[droppedterms] <- -Inf
   model$etamap$offsettheta[names(model$etamap$offsettheta) %in% namesdrop] <- TRUE
  
  
  

  Clist <- ergm.Cprepare(nw, model)
  Clist.miss <- ergm.design(nw, model, verbose=verbose)
  Clist$obs <- summary(model$formula, drop=FALSE)


  Clist$meanstats <- Clist$obs
  

  MCMCparams=c(control, list(samplesize=MCMCsamplesize, burnin=burnin,
  interval=interval, stepMCMCsize=control$stepMCMCsize, gridsize=control$gridsize,
  maxit=maxit,Clist.miss=Clist.miss, mcmc.precision=control$mcmc.precision))

 return(list(theta0=theta0,nw=nw,model=model,Clist=Clist,initialfit=initialfit,MCMCparams=MCMCparams,
  MHproposal=MHproposal,MHproposal.miss=MHproposal,verbose=verbose))
 
}

 #other functions to transform parameters

para1<-function(Clist,MCMCparams){
statshift <- Clist$obs - Clist$meanstats
MCMCparams$meanstats <- Clist$meanstats
MCMCparams$nmatrixentries <- MCMCparams$samplesize * Clist$nstats
return(list(MCMCparams=MCMCparams,statshift=statshift))
}


#obtain weights, probabilities for network given parameters
obtain.weights<-function(Z.eta,eta.new,eta0,Z.Psi.new,Z.Psi0){
#statsmatrix and eta have same length  
N<-dim(Z.eta)[1]
p1<-length(eta.new)
p2<-dim(Z.eta)[2]
if(p1!=p2){cat("eta and statsmatrix have different length - error!\n");return}
w<-exp(Z.eta%*%as.matrix(c(eta.new-eta0)) +  as.matrix(Z.Psi.new-Z.Psi0))
#for(i in 1:N){
#w[i]<-exp(t(eta-eta0)%*%Z[i,] + Z.Psi.new - Z.Psi.old)
#end  for(i in 1:N){
den<-sum(w)
U<-w
w<-w/sum(w)
return(list(w=c(w),U=c(U),den=den))
}

ERGM.ratio<-function(U,Z.obs,eta.new,eta.old,Y.obs,Psi.new,Psi.old,Z){
ratio<-(eta.new-eta.old)%*%Z.obs + sum((Psi.new-Psi.old)*Y.obs) - log(mean(U))
return(ratio)
}



show.elapsed.time<-function(start.time,string=""){
 end.time<-proc.time()[3]
 time<- end.time-start.time;
 hours<-floor(time/3600);
 minutes<-floor(time/60)-hours*60;
 seconds<-time-hours*3600-minutes*60;
 cat(string,"Time: ",hours,"h",minutes,"min",seconds,"sec","\n")
 return(time)
}#End FUNCTION: show.elapsed.time


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){

cov.w<-function(Y1,Y2,w){
N<-length(w)
m.Y1<-t(matrix(1,1,N)%*%(w*Y1))
m.Y2<-t(matrix(1,1,N)%*%(w*Y2))
m.Y1.Y2<- t(Y2)%*%(w*Y1)
cov.Y1.Y2<- m.Y1.Y2 -  m.Y2%*%t(m.Y1)
return(list(m.Y1=m.Y1,m.Y2=m.Y2,cov.Y1.Y2=cov.Y1.Y2))
}#end  cov.w<-function(Y1,Y2,w){

auto.cov<-function(Y,K){
M<-dim(Y)[1] 
N<-dim(Y)[2]
auto.cov.Y<-array(0,dim=c(K+1,N,N))
m.Y  <- t(matrix(1,1,M)%*%Y)/M
m.YY <- m.Y%*%t(m.Y)
for(k in 1:(K+1)){
Yi  <- Y[1:(M-k+1),]
Yik <- Y[k:M,]
auto.cov.Y[k,,]<- (1/(M-k+1))*(t(Yi)%*%Yik) - m.YY 
}#end k
cov.Y<- auto.cov.Y[1,,]
for(k in 1:K){
cov.Y <- cov.Y + 2*auto.cov.Y[k+1,,]
}
#cov.Y0<-(t(Y)%*%(Y))/M - m.YY))
#print(sum(abs(auto.cov.Y[1,,]-cov.Y0)))
return(list(auto.cov.Y=auto.cov.Y,cov.Y=cov.Y,m.Y=m.Y))
}#end auto.cov


get.updated.Psi.statistic<-function(Y.MCMC,n,Psi.new){
MCMCsamplesize<-dim(Y.MCMC)[1]
ind.lt<-lower.tri(diag(n))
stats.Psi<-matrix(0,MCMCsamplesize,1)
#Yobs.Psi.new<-sum(Y.obs*c(Psi.new))
for(i in 1:MCMCsamplesize){
stats.Psi[i]<- sum(Y.MCMC[i,]*c(Psi.new))  
#cat("i: ",i,"\n")
}#end  for(i in 1:MCMCsamplesize){
return(stats.Psi=stats.Psi)
}#end 







logit<-function(mu){log(mu)-log(1-mu)}

expit<-function(eta){exp(eta)/(1+exp(eta))}

norm<-function(d){d<-as.matrix(d);t(d)%*%d}



#############################################################################################################
#                 BEGIN OBTAIN MCMC sample
############################################################################################################

get.MCMC<-function(formula="work ~ gwesp(1.099,fixed=TRUE)",Psi.vector,theta,n,MCMCsamplesize=1e4,nw=work,MCMCstep=1e3){
#update Psi
     formula1<-as.formula(paste(formula," + edgecov(Psi)",sep=""))
     print(formula1)

     Psi<-matrix(0,n,n);
     Psi[ind.lt]<- Psi.vector
     Psi<-Psi+t(Psi)
     p1<-get.para(formula1)
     p1$MCMCparams$samplesize<-1 #e3
    
     hilf<-para1(p1$Clist,p1$MCMCparams)
     #statshift1<-hilf$statshift
     MCMCparams1<-hilf$MCMC
     #MCMCparams1$stats<-MCMCparams1$meanstats
     Z.obs<-MCMCparams1$meanstats[1]
     
     theta01<-c(theta,1)  #eta-parameter for Psi is fixed to 1
     model1<-p1$model
     
    
    
    #Model 1 is with gwesp and Psi
    eta01  <-  ergm.eta(theta01, model1$etamap)
    Clist1 <- ergm.Cprepare(nw, model1)
	  # Obtain MCMC sample
	  
	  cat("eta: ", eta01,"\n")
	  cat("theta: ",theta01,"\n")
	  
     #print(Clist1[1:12])
     #print(eta01)
     #print(MCMCparams1)
    
     #last parameter stands for 
	  
	  Z<-matrix(0,MCMCsamplesize,l.p+1) #plus 1 additional term for Psi
    Y.MCMC<-matrix(0,MCMCsamplesize,Ny)
	  
	  #Y.networks<-list()
	  start.time<-proc.time()[3]
	  #to make it shorter take only 5,000 burnin's otherwise 10,000
	  MCMCparams1$burnin<-10e3
    for(j in 1:MCMCsamplesize){
     z1 <- ergm.getMCMCsample(Clist1, p1$MHproposal, eta01, MCMCparams1, verbose=F)
     #z1 <- ergm.getMCMCsample1(Clist1, p1$MHproposal, eta01, MCMCparams1, verbose=F)
     MCMCparams1$burnin<-MCMCstep
     Z[j,]<-z1[[1]]
     updated.network<-network.update(nw, z1$newedgelist, "edgelist")
     #statsmatrix2[j,]<-summary(updated.network~gwesp(1.099,fixed=TRUE) + edgecov(Psi))
     #update network, restart at previous sampled network
     Clist1 <- ergm.Cprepare(updated.network, model1)
     #Y.networks[[j]]<-hilf.network
     Y.MCMC[j,]<- as.sociomatrix(updated.network)[ind.lt]
     #cat("j: ",j,"\n")
    if(j%%500==0 | j==1){end.time<-show.elapsed.time(start.time,string=paste("j: ",j,sep=""))}
    }#  for(j in 1:MCMCsamplesize){
    
    Z.hilf<-Z
    
    #compute true Z
    Z<-apply(rbind(MCMCparams1$meanstats,Z),2,cumsum)[2:(MCMCsamplesize+1),]
    
    Z.Psi.old<-Z[,l.p+1]
    Z<-as.matrix(Z[,1:l.p])
    
    #MCMCparams1$meanstats[2]
    
    eta0<-eta01[1:l.p]
    
    return(list(Z.Psi=Z.Psi.old,Y.MCMC=Y.MCMC,eta0=eta0,Z=Z,Clist=Clist1,MCMCparams=MCMCparams1,eta01=eta01,Z.obs=Z.obs))
    
}#enf function get.MCMC



#############################################################################################################
#                END  OBTAIN MCMC sample
############################################################################################################


################################################################################
get.shared.partner<-function(netw,n){
#max.esp<-max(summary(netw~esp(1:(n-2))))
ESP<-matrix(0,n-2,n*(n-1)/2)
k<-0
for(i in 1:(n-1)){
 for(j in (i+1):n){
  k<-k+1
  y1<-netw
  y0<-netw
  y1[i,j]<-1
  y0[i,j]<-0
  ESP[,k]<-summary(y1~esp(1:(n-2)))-summary(y0~esp(1:(n-2)))  
 }#  for(j in (i+1):n){
}#  for(i in 1:(n-1)){
return(ESP)
}#end function get.shared.partner
################################################################################

################################################################################
get.shared.partner.own<-function(y,n){
#network must be a matrix
if(is.network(y)){y<-as.matrix.network(y)}

Ny<-n*(n-1)/2

ind.lt<-lower.tri(y)


#max.esp<-max(summary(as.network(y,directed=F)~esp(0:(n-2))))
delta.ESP1<-delta.ESP<-matrix(0,Ny,n-2)

ind<-indexfct(n)

SP<-(y%*%y)
diag(SP)<-0

for(k in 1:Ny){
   i<-ind[1,k]
   j<-ind[2,k]
   y1<-y0<-y
   y1[i,j]<-y1[j,i]<-1
   y0[i,j]<-y0[j,i]<-0
 
   SP0<- SP1 <- SP
   SP1[i,]<-SP1[,i]<- SP[i,] + (1-y[i,j])*y[j,]
   SP1[j,]<-SP1[,j]<- SP[j,] + (1-y[i,j])*y[i,]
    
   SP0[i,]<-SP0[,i]<- SP[i,] - y[i,j]*y[j,]
   SP0[j,]<-SP0[,j]<- SP[j,] - y[i,j]*y[i,]
   
   #SP1[i,j]<-SP[i,j];SP0[i,j]<-SP[i,j]
   #SP1[i,i]<-SP1[j,j]<-SP0[i,i]<-SP0[j,j]<-0 
   diag(SP0)<-0
   diag(SP1)<-0
   
   hilf1<-(y1*SP1)[ind.lt]
   hilf0<-(y0*SP0)[ind.lt]
   ESP1<-tabulate(hilf1[hilf1>0])
   ESP0<-tabulate(hilf0[hilf0>0])
   l.ESP1<-length(ESP1)
   l.ESP0<-length(ESP0)
   
   if(l.ESP1==l.ESP0){
   hilf<-ESP1-ESP0
   }else{
   if(l.ESP1>l.ESP0){
    hilf<-ESP1-c(ESP0,rep(0,l.ESP1-l.ESP0))
   }else{
    hilf<- -ESP0+c(ESP1,rep(0,l.ESP0-l.ESP1))
   }
   }
      
   delta.ESP[k,1:length(hilf)]<-hilf
   
   #delta.ESP1[k,]<-summary(network(y1,directed=F)~esp(1:max.esp))-summary(network(y0,directed=F)~esp(1:max.esp))
     
 }#  for(i in 1:Ny){

return(delta.ESP)

}#end function get.shared.partner
################################################################################


################################################################################
get.shared.partner.own.ij<-function(y,n){
#network must be a matrix
if(is.network(y)){y<-as.matrix.network(y)}
ind.lt<-lower.tri(y)
Ny<-n*(n-1)/2

ind.lt<-lower.tri(y)


deltaESPij<-array(0,dim=c(Ny,Ny,n-2))

ind<-indexfct(n)

SP<-(y%*%y)
diag(SP)<-0

start.time<-proc.time()[3]



for(k in 1:(Ny-1)){
#cat(k,"th iteration \n")
if(k%%10==0 | k==1){end.time<-show.elapsed.time(start.time,string=paste("k: ",k, "out of ",Ny,sep=""))}
   for(l in (k+1):Ny){
   
   i<-ind[1,k]
   j<-ind[2,k]
   s<-ind[1,l]
   t<-ind[2,l]
   
   y11<-y00<-y10<-y01<-y
   y11[i,j]<-y11[j,i]<-1;y11[s,t]<-y11[t,s]<-1
   y00[i,j]<-y00[j,i]<-0;y00[s,t]<-y00[t,s]<-0
   y10[i,j]<-y10[j,i]<-1;y10[s,t]<-y10[t,s]<-0
   y01[i,j]<-y01[j,i]<-0;y01[s,t]<-y01[t,s]<-1
   
   #first index (1) for ij and 2nd index (1) for st
 
   SP10<-SP01<-SP00<- SP11 <- SP
   
   
   SP11[i,]<-SP11[,i]<- SP[i,] + (1-y[i,j])*y[j,]
   SP11[j,]<-SP11[,j]<- SP[j,] + (1-y[i,j])*y[i,]
   SP11[s,]<-SP11[,s]<- SP[s,] + (1-y[i,j])*y[t,]
   SP11[t,]<-SP11[,t]<- SP[t,] + (1-y[i,j])*y[s,]
   
   SP10[i,]<-SP10[,i]<- SP[i,] + (1-y[i,j])*y[j,]
   SP10[j,]<-SP10[,j]<- SP[j,] + (1-y[i,j])*y[i,]
   SP10[s,]<-SP10[,s]<- SP[s,]      -y[i,j]*y[t,]
   SP10[t,]<-SP10[,t]<- SP[t,]      -y[i,j]*y[s,]
   
   SP01[i,]<-SP01[,i]<- SP[i,]     -y[i,j] *y[j,]
   SP01[j,]<-SP01[,j]<- SP[j,]     -y[i,j] *y[i,]
   SP01[s,]<-SP01[,s]<- SP[s,] + (1-y[i,j])*y[t,]
   SP01[t,]<-SP01[,t]<- SP[t,] + (1-y[i,j])*y[s,]
   
   SP00[i,]<-SP00[,i]<- SP[i,]     -y[i,j] *y[j,]
   SP00[j,]<-SP00[,j]<- SP[j,]     -y[i,j] *y[i,]
   SP00[s,]<-SP00[,s]<- SP[s,]     -y[i,j] *y[t,]
   SP00[t,]<-SP00[,t]<- SP[t,]     -y[i,j] *y[s,] 
   
   #i,s
   SP11[s,i]<-SP11[i,s]<-SP[s,i]+(1-y[i,j])*y[s,j] + (1-y[s,t])*y[i,t]
   SP10[s,i]<-SP10[i,s]<-SP[s,i]+(1-y[i,j])*y[s,j]     -y[s,t] *y[i,t]
   SP01[s,i]<-SP01[i,s]<-SP[s,i]   -y[i,j] *y[s,j] + (1-y[s,t])*y[i,t]
   SP00[s,i]<-SP00[i,s]<-SP[s,i]   -y[i,j] *y[s,j]     -y[s,t] *y[i,t]
   
   #i,t
   SP11[t,i]<-SP11[i,t]<-SP[t,i]+(1-y[i,j])*y[t,j] + (1-y[s,t])*y[i,s]
   SP10[t,i]<-SP10[i,t]<-SP[t,i]+(1-y[i,j])*y[t,j]     -y[s,t] *y[i,s]
   SP01[t,i]<-SP01[i,t]<-SP[t,i]   -y[i,j] *y[t,j] + (1-y[s,t])*y[i,s]
   SP00[t,i]<-SP00[i,t]<-SP[t,i]   -y[i,j] *y[t,j]     -y[s,t] *y[i,s]
   
   #j,s
   SP11[s,j]<-SP11[j,s]<-SP[s,j]+(1-y[i,j])*y[s,i] + (1-y[s,t])*y[j,t]
   SP10[s,j]<-SP10[j,s]<-SP[s,j]+(1-y[i,j])*y[s,i]     -y[s,t] *y[j,t]
   SP01[s,j]<-SP01[j,s]<-SP[s,j]   -y[i,j] *y[s,i] + (1-y[s,t])*y[j,t]
   SP00[s,j]<-SP00[j,s]<-SP[s,j]   -y[i,j] *y[s,i]     -y[s,t] *y[j,t]
   
   
   #j,t
   SP11[t,j]<-SP11[j,t]<-SP[t,j]+(1-y[i,j])*y[t,i] + (1-y[s,t])*y[i,t]
   SP10[t,j]<-SP10[j,t]<-SP[t,j]+(1-y[i,j])*y[t,i]     -y[s,t] *y[i,t]
   SP01[t,j]<-SP01[j,t]<-SP[t,j]   -y[i,j] *y[t,i] + (1-y[s,t])*y[i,t]
   SP00[t,j]<-SP00[j,t]<-SP[t,j]   -y[i,j] *y[t,i]     -y[s,t] *y[i,t]
   
   
   diag(SP11)<-0
   diag(SP10)<-0
   diag(SP01)<-0
   diag(SP00)<-0
   
   hilf11<-(y11*SP11)[ind.lt]
   hilf10<-(y10*SP10)[ind.lt]
   hilf01<-(y01*SP01)[ind.lt]
   hilf00<-(y00*SP00)[ind.lt]
   
   ESP11<-tabulate(hilf11[hilf11>0])
   ESP10<-tabulate(hilf10[hilf10>0])
   ESP01<-tabulate(hilf01[hilf01>0])
   ESP00<-tabulate(hilf00[hilf00>0])
   
   l.ESP11<-length(ESP11);l.ESP01<-length(ESP01)
   l.ESP10<-length(ESP10);l.ESP00<-length(ESP00)
   
   ml<-max(c(l.ESP11,l.ESP10,l.ESP01,l.ESP00))
   
   hilf<-c(ESP11,rep(0,ml-l.ESP11))-c(ESP10,rep(0,ml-l.ESP10))-c(ESP01,rep(0,ml-l.ESP01))+c(ESP00,rep(0,ml-l.ESP00))
   
      
   deltaESPij[k,l,1:ml]<- deltaESPij[l,k,1:ml]<-hilf
   
     
   }#  for(l in 1:Ny){ 
 }#  for(k in 1:Ny){

return(deltaESPij)

}#end function get.shared.partner
################################################################################



################################################################################
get.degree<-function(network,n){
Deg<-matrix(0,n-1,n*(n-1)/2)
k<-0
for(i in 1:(n-1)){
 for(j in (i+1):n){
  k<-k+1
  y1<-network
  y0<-network
  y1[i,j]<-1
  y0[i,j]<-0
  Deg[,k]<-summary(y1~degree(1:(n-1)))-summary(y0~degree(1:(n-1)))  
 }#  for(j in (i+1):n){
}#  for(i in 1:(n-1)){
return(Deg)
}#end function get.shared.partner
################################################################################ 

################################################################################
get.degree.own<-function(network,n){
Deg<-matrix(0,n-1,n*(n-1)/2)
k<-0
for(i in 1:(n-1)){
 for(j in (i+1):n){
  k<-k+1
  y1<-network
  y0<-network
  y1[i,j]<-1
  y0[i,j]<-0
  Deg[,k]<-summary(y1~degree(1:(n-1)))-summary(y0~degree(1:(n-1)))  
 }#  for(j in (i+1):n){
}#  for(i in 1:(n-1)){
return(Deg)
}#end function get.shared.partner
################################################################################



# function for indexing odds
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


#####################################################################################################
#   GWESP
#######################################################################################

gwesp<-function(gwesp.tau,gwesp.alpha,n){
 eta<-gwesp.tau*exp(gwesp.alpha)*(1-(1-exp(-gwesp.alpha))^(1:(n-2)))
 return(eta)
}#end function 

gwesp.der<-function(gwesp.tau,gwesp.alpha,n){
 eta<-gwesp.tau*exp(gwesp.alpha)*(1-(1-exp(-gwesp.alpha))^(1:(n-2))) - gwesp.tau*(1:(n-2))*(1-exp(-gwesp.alpha))^(0:(n-3))
 return(eta)
}#end function 

gwesp.der.2<-function(gwesp.tau,gwesp.alpha,n){
 eta<- gwesp.tau*exp(gwesp.alpha)*(1-(1-exp(-gwesp.alpha))^(1:(n-2))) - gwesp.tau*(1:(n-2))*(1-exp(-gwesp.alpha))^(0:(n-3)) - 
 gwesp.tau*((1:(n-2)))*(0:(n-3))*(1-exp(-gwesp.alpha))^(c(0,0:(n-4)))
return(eta)
}




################################################################################
Psi.matrix<-function(Psi.vector,n){
hilf.matrix<-diag(rep(0,n))
ind.lt<-lower.tri(hilf.matrix)
hilf.matrix[ind.lt]<-Psi.vector
hilf.matrix<-hilf.matrix+t(hilf.matrix)
hilf.matrix
}#end  Psi.matrix<-function(Psi.old,n){



################################################################################
#              get solution for Psi subject to model for mu
################################################################################
solve.Psi.GWESP<-function(X,beta.hat,gwesp.tau,gwesp.alpha=0,ESP,Y.obs,nw0=work,Psi0=NULL,fit="NR",MCMCsamplesize=1e4,iter.max=20,MCMCstep=1e3){
Psi.old<-Psi0
Y.MCMC<-NULL

Ny<-length(Y.obs)
n<-round(abs(polyroot(c(-2*630,-1,1)))[1])

#obtain estimates for Psi via pseudolikelihood and change-statistic ESP (change of shared partner statistic)
eta<-gwesp(gwesp.tau,gwesp.alpha,n)
X.beta.hat<- X%*%beta.hat
Psi.pseudo <- X.beta.hat - ESP%*%eta
 

#within get.MCMC edgecov(Psi) is added
formula1<-paste("nw~gwesp(",gwesp.alpha,",fixed=TRUE)",sep="")

expit.X.beta.hat<-expit(X.beta.hat)
mu.X.beta<-expit.X.beta.hat


#get starting values for Psi
#if(is.null(Psi.old)){
if(is.null(Psi.old)){
Psi.old<-Psi.pseudo
Psi.old<-pmax(Psi.old,-7)
Psi.old<-pmin(Psi.old,+3)
}







mu.Y.Psi<-mu.X.beta+0.01

Psi0<-Psi.old


iter<-0
iter.non.MCMC<-0

diff.mu<-1
diff.mu.unscaled<-1
diff.Psi<-1
var.MC<-1

diff.mu.list<-NULL
diff.mu.unscaled.list<-NULL
diff.Psi.list<-NULL
Psi.list<-NULL

#iter.max<-20
X2<-Inf

M<-MCMCsamplesize
p<-Ny

#diff.mu is scaled difference 
while( ((X2 >  qchisq(0.95,df=Ny))  & iter<iter.max & diff.mu>0.1) | iter<3){
#if(1){
iter<-iter+1
cat("########################################################################\n")
cat("   Iteration ",iter,"\n")
cat("########################################################################\n")

cat("beta.hat:",beta.hat,"\n")
cat("gwesp.tau: ",gwesp.tau,"  gwesp.alpha: ",gwesp.alpha,"\n")
cat("norm Difference between mu.X.beta and mu.Psi:",diff.mu,"\n")
cat("max unscaled Difference between mu.X.beta and mu.Psi:",diff.mu.unscaled,"\n")
cat("100*norm unscaled Difference between mu.X.beta and mu.Psi:",100*norm(mu.X.beta - mu.Y.Psi)," compared to sum(Var(Y)): ",sum(mu.Y.Psi*(1-mu.Y.Psi)/MCMCsamplesize),"\n")
cat("mean unscaled Difference between mu.X.beta and mu.Psi:",mean(abs(mu.X.beta - mu.Y.Psi)),"\n")

cat("diff Psi: ",diff.Psi,"\n")

#we always have an old Psi and old gwesp and possibly new gwesp, but NO new Psi, since we want to update it

Psi.old<-Psi.new  

w<-c(0,1,2);var.MC<-1

if(iter>1){
#obtain Z.Psi.new for new parameters Psi.new
Z.Psi.new<-get.updated.Psi.statistic(Y.MCMC=Y.MCMC,n,Psi.new=Psi.new)
hilf<-obtain.weights(Z.eta=Z,eta.new=gwesp.tau,eta0=gwesp.tau,Z.Psi.new=Z.Psi.new,Z.Psi0=Z.Psi.old) 
if(sum(is.na(hilf$w))==0){
w<-hilf$w
U<-hilf$U
den<-hilf$den
U.bar<-mean(U)
auto.lags<-acf(U,plot=FALSE,type="covariance")$acf
k<-length(auto.lags)
#formula says sum_k=K to K (M -|k|)*gamma k, where gamma_k=gamma_ -k 
var.MC<-(1/(MCMCsamplesize^2*U.bar^2))*sum(c(2*(MCMCsamplesize-(1:(k-1)))*auto.lags[2:k],MCMCsamplesize*auto.lags[1]))

#ratio.0.eta.old<-ERGM.ratio(U,Z.obs,eta.new=gwesp.tau,eta.old=0,Y.obs,Psi.new=Psi.old,Psi.old=rep(0,Ny),Z)
#ratio.eta.new.eta.old<-ERGM.ratio(U,Z.obs,eta.new=gwesp.tau,eta.old=gwesp.tau,Y.obs,Psi.new=Psi.new,Psi.old=Psi.old,Z)
#lik<-ratio.eta.new.eta.old - ratio.0.eta.old - Ny*log(2)

cat("var - w: ", var(w), "  mean - w: ",mean(w)," max -w: ", max(w) ," var.MC: ",var.MC," \n")
}else{w<-c(0,1,2);var.MC<-1}
}
          
new.MCMC.sample<-FALSE
if(var.MC > 0.01 | max(w)>0.1| max(w)==0  | iter==1 | iter.non.MCMC>=2  | diff.mu.unscaled>0.05){
iter.non.MCMC<-0
new.MCMC.sample<-TRUE
Y.MCMC.old<-Y.MCMC
cat("\n ############    obtain MCMC sample    ############ \n\n")
w<-rep(1/MCMCsamplesize,MCMCsamplesize)
hilf<-get.MCMC(formula=formula1,Psi.vector=Psi.old,theta=gwesp.tau,n=n,MCMCsamplesize,nw=nw0,MCMCstep=MCMCstep)
Y.MCMC<-hilf$Y.MCMC
Z.Psi.old<-hilf$Z.Psi #assign Z.Psi.old only here
Z<-hilf$Z
Z.obs<-hilf$Z.obs
}else{
iter.non.MCMC<-iter.non.MCMC+1
cat("iterations without getting new MCMC sample: ",iter.non.MCMC," \n")
}

# obtain MCMC.sample

cat("similarity between MCMC sample and Y.obs:  Z.obs=",Z.obs," mean(Z)=",sum(w*Z),"\n")




###############################################################################
#separate solving for each Psi
##############################################################################
if(fit=="direct"){
#compute change statistics for all Y's in Y.MCMC
formula.ij.1<-as.formula(paste("Y.network.ij.1~gwesp(",gwesp.alpha,",fixed=TRUE)",sep=""))
formula.ij.0<-as.formula(paste("Y.network.ij.0~gwesp(",gwesp.alpha,",fixed=TRUE)",sep=""))


if(new.MCMC.sample){
cat("\n#      get delta.Z for MCMC sample          #\n\n")
#create new delta.Z.gwesp for new MCMC sample
delta.Z.gwesp<-matrix(0,MCMCsamplesize,Ny)
start.time<-proc.time()[3]
for(i in 1:MCMCsamplesize){
Y.matrix.i<-Psi.matrix(Y.MCMC[i,],n)
ESP<-get.shared.partner.own(Y.matrix.i,n)
delta.Z.gwesp[i,]<-ESP%*%eta
if(i%%200==0 | i==1){end.time<-show.elapsed.time(start.time,string=paste("i: ",i, "out of ",MCMCsamplesize,sep=""))}
}#end  for(i in 1:MCMCsamplesize){
}#end if new.MCMC.sample


obj.fct<-function(x,mu.X.beta,delta.Z,w){mu.X.beta-sum(expit(x + delta.Z)*w)}

cat("\n#      Solve for Psi          #\n\n")
Psi.hilf<-rep(0,Ny,1)
for(j in 1:Ny){
Psi.hilf[j]<-Psi.old[j] + 0.2*(uniroot(obj.fct,c(-50,+50),mu.X.beta=mu.X.beta[j],delta.Z=delta.Z.gwesp[,j],w=w)$root-Psi.old[j])
if(j%%100==0 | j==1 | j==Ny){end.time<-show.elapsed.time(start.time,string=paste(j,"jth parameter out of ",Ny,sep=""))}
}
Psi.new<-Psi.hilf
Psi.global<<-Psi.new

hilf3<-mean.var.w(Y.MCMC,w,binary=TRUE)
mu.Y.Psi<-hilf3$m.Y
C11<-cov.Y<-hilf3$cov.Y

#save(Y.MCMC,Y.MCMC.old,Psi.new,Psi.old,file="SolvePsi.RData")
#joint solving - doesn't work    
}

################################################################################
#                   simple method 
###############################################################################
if(fit=="simple"){
#modify Psi.old slightly
hilf3<-mean.var.w(Y.MCMC,w,binary=TRUE)
C11<-hilf3$cov.Y
mu.Y.Psi <- hilf3$m.Y
hilf<-  (mu.X.beta - mu.Y.Psi)/(mu.Y.Psi*(1-mu.Y.Psi))   #(mu.X.beta*(1-mu.X.beta)) 
#print(c(max(Psi.old),min(Psi.old)))
cat("stepsize for simple method: ",0.1*(iter/(iter.max))^(0.8),"\n")
Psi.new  <- Psi.old + sign(hilf)*pmin(abs(hilf),sqrt(abs(hilf)))*0.1*(iter/(4*iter.max))^(0.8)
Psi.new<-pmax(Psi.new,-7)
Psi.new<-pmin(Psi.new,+3)
cat("max(delta-Psi): ",max(hilf)," min(delta-Psi): ",min(hilf),"\n")
Psi.global<<-Psi.new
#save(Y.MCMC,Y.MCMC.old,Psi.new,Psi.old,file="SolvePsi.RData")
}


if(fit=="NR"){
#################################################################################
#iteration scheme to solve for Psi 
#################################################################################
cat("\n#####Least squares scheme for Psi ### \n\n")
#print(dim(Y.MCMC))
#print(length(w))

#compute covariance from MCMC-sample


hilf3<-mean.var.w(Y.MCMC,w,binary=TRUE)
#hilf2<-auto.cov(Y.MCMC,K=20)

#C11.auto<-hilf2$cov.Y    
#C11.auto.inv<-solve(C11.auto)
    
mu.Y.Psi<-hilf3$m.Y
C11<-cov.Y<-hilf3$cov.Y
v.Y<-hilf3$v.Y
v.Y.sqrt<-sqrt(v.Y)
v.Y.sqrt.inv<-1/v.Y.sqrt

#v.Y.sqrt.Xbeta<-sqrt(c(mu.X.beta)*(1-c(mu.X.beta)))
#v.Y.sqrt.inv.Xbeta<-1/sqrt(c(mu.X.beta)*(1-c(mu.X.beta)))   

#lambda determines how much we use  NR vs. simple
lambda <- 1 - (max(abs(mu.X.beta - mu.Y.Psi))/0.2) 
if(max(abs(mu.X.beta - mu.Y.Psi))>0.2){lambda<-0}
if(max(abs(mu.X.beta - mu.Y.Psi))<0.05){lambda<-1}
    
R11<- v.Y.sqrt.inv*t(v.Y.sqrt.inv*C11) #==diag(v.Y.sqrt.inv)%*%cov.Y%*%diag(v.Y.sqrt.inv)
R11<-lambda*R11 + (1-lambda)*diag(Ny)
R11.inv<-solve(R11)

C11    <-  v.Y.sqrt*t(v.Y.sqrt*R11) 
C11.inv<-  v.Y.sqrt.inv*t(v.Y.sqrt.inv*R11.inv)


delta.Y<-(mu.X.beta-mu.Y.Psi) 
#eigen.value<-min(eigen(C11)$values)

delta.Psi<-C11.inv%*%delta.Y
#delta.Psi<-C11.auto.inv%*%delta.Y

cat("max(delta-Psi): ",max(delta.Psi)," min(delta-Psi): ",min(delta.Psi)," mean(delta.Psi): ",mean(delta.Psi),"\n")
#delta.Psi<- 1/min(norm(delta.Psi),100) * delta.Psi

off.set<-0.01
#step.size<- 0.01 + min(2*(max(max(abs(mu.X.beta - mu.Y.Psi)),0.01)),1)
step.size<- 0.005+ 0.495*((iter.max-iter+1)/iter.max)^2
#step.size<-1
cat("Stepsize: ",step.size,"\n")
cat("NR-simple-Parameter (0-simple,1 -NR): ",lambda,"\n")


Psi.new<- Psi.old + step.size*delta.Psi 
Psi.new<-pmax(Psi.new,-7)
Psi.new<-pmin(Psi.new,+2)
Psi.global<<-Psi.new


#save(Y.MCMC,Y.MCMC.old,Psi.new,Psi.old,file="SolvePsi.RData")
#end iteration scheme
##############################################################################
}#end if 

#compute C11.inv = Cov.Y.inv
if(fit!="NR"){
v.Y<-hilf3$v.Y
v.Y.sqrt<-sqrt(v.Y)
v.Y.sqrt.inv<-1/v.Y.sqrt
R11<- v.Y.sqrt.inv*t(v.Y.sqrt.inv*C11) #==diag(v.Y.sqrt.inv)%*%cov.Y%*%diag(v.Y.sqrt.inv)
R11.inv<-solve(R11)
C11    <-  v.Y.sqrt*t(v.Y.sqrt*R11) 
C11.inv<-  v.Y.sqrt.inv*t(v.Y.sqrt.inv*R11.inv)
}


if(iter==1){
#obtain  update for Z.Psi
Z.Psi.old<-get.updated.Psi.statistic(Y.MCMC=Y.MCMC,n,Psi.new=Psi.old)
}else{
Z.Psi.old<-Z.Psi.new
} 

cat("  mu.Y.Psi - mean of ERGM specified mainly by Psi, because theta is fixed \n")
cat("  mu.X.beta - mean of GLM is fixed, because beta is fixed \n")   
cat(" Aim is to find Psi such that mu.Y.Psi=mu.X.beta\n")
cat(" max(mu.Y.Psi): ",max(mu.Y.Psi),"  min(mu.Y.Psi): ", min(mu.Y.Psi),"\n") 
cat(" max(mu.X.beta): ",max(mu.X.beta),"  min(mu.X.beta): ", min(mu.X.beta),"\n") 
 
diff.mu<-norm((mu.X.beta - mu.Y.Psi)/(mu.Y.Psi*(1-mu.Y.Psi)))    #(mu.X.beta*(1-mu.X.beta)))
diff.mu.unscaled<-max(abs(mu.X.beta - mu.Y.Psi)) 
diff.Psi<- norm(Psi.new-Psi.old)



cat("Eucledean norm of Psi.new-Psi.old: ",diff.Psi,"\n") 
cat("Max(Psi.new): :",max(Psi.new),"  Min(Psi.new): ",min(Psi.new),"\n") 
cat("Max(Psi.old): :",max(Psi.old),"  Min(Psi.old): ",min(Psi.old),"\n") 
  


diff.mu.list  <-  c(diff.mu.list,diff.mu)
diff.mu.unscaled.list<- c(diff.mu.unscaled.list,diff.mu.unscaled)
diff.Psi.list  <- c(diff.Psi.list,diff.Psi)
Psi.list<-cbind(Psi.list,Psi.new)

save(Y.MCMC,Y.MCMC.old,Psi.new,Psi.old,beta.hat,gwesp.alpha,gwesp.tau,file="SolvePsi.RData")



#X2<-M*t(mu.X.beta - mu.Y.Psi)%*%C11.inv%*%(mu.X.beta - mu.Y.Psi)                                
X2<-norm((M*mu.X.beta - M*mu.Y.Psi)/sqrt(M*mu.Y.Psi))
#cat("M*t(mu.X.beta - mu.Y.Psi)%*%C11.inv%*%(mu.X.beta - mu.Y.Psi): ",X2,"\n")
cat("test-statistic comparing  mu.X.beta with  mu.Y.Psi :(M*mu.X.beta - M*mu.Y.Psi)/sqrt(M*mu.Y.Psi)",X2,"\n")
cat("X2 -df=",Ny," alpha 95%: ",qchisq(0.95,df=Ny),"\n")



}#end while



return(list(Y.MCMC=Y.MCMC,Psi.new=Psi.old,Psi.old=Psi0,diff.mu=diff.mu,w=w,U=U,
diff.Psi=diff.Psi,Z=Z,Z.Psi=Z.Psi.new,Z.obs=Z.obs,diff.mu.list1=diff.mu.list,
diff.mu.list2=diff.mu.unscaled.list,diff.Psi.list=diff.Psi.list,Psi.list=Psi.list))



############################################################################################################
}#end function  solve.Psi.GWESP
#######################################################################################################


#logistic regression for gwesp0 and gwesp
logistic.gwesp<-function(g,ESP,Psi,Y,n){
gwesp.tau<-g[1]
gwesp.alpha<-g[2]
Psi.par<-g[3]
eta<-gwesp(gwesp.tau,gwesp.alpha,n)
ESP.eta<-ESP%*%eta
Xbeta<- g[3]*Psi + ESP.eta 

#because we want to maximize function, but nle only minimises, we need "minus" 
res<-  - sum(Y*Xbeta - log(1+exp(Xbeta)) )

eta.der.alpha<-gwesp.der(gwesp.tau,gwesp.alpha,n)
eta.der.tau<-gwesp(1,gwesp.alpha,n)

ESP.eta.alpha <- ESP%*%eta.der.alpha
ESP.eta.tau   <- ESP%*%eta.der.tau

gr<-Y-expit(Xbeta)
gr<- - c(sum(gr*ESP.eta.tau),sum(gr*ESP.eta.alpha),sum(gr*Psi))
attr(res, "gradient") <- gr

return(res)
}#end  logistic.gwesp<-function(g,ESP,Psi,Y.obs){


compute.cov.gwesp<-function(gwesp.tau,gwesp.alpha,Psi.par,ESP,Psi,Y,n){
eta<-gwesp(gwesp.tau,gwesp.alpha,n)
eta.der.alpha<-gwesp.der(gwesp.tau,gwesp.alpha,n)
eta.der2.alpha<-gwesp.der.2(gwesp.tau,gwesp.alpha,n)
eta.der.tau<-gwesp(1,gwesp.alpha,n)

ESP.eta <- ESP%*%eta
Xbeta <- Psi + ESP.eta 
mu<-expit(Xbeta)

ESP.eta.der.tau <- ESP%*%eta.der.tau 
ESP.eta.der.alpha  <-ESP%*%eta.der.alpha
ESP.eta.der2.alpha  <-ESP%*%eta.der2.alpha


XWX<-matrix(0,3,3)
XWX[1,1] <- -sum( mu*(1-mu)*ESP.eta.der.tau^2)
XWX[2,2] <- -sum( mu*(1-mu)*ESP.eta.der.alpha^2) - sum((Y-mu)*ESP.eta.der2.alpha) 
XWX[1,2] <- XWX[2,1] <- -sum(mu*(1-mu)*ESP.eta.der.tau*ESP.eta.der.alpha)
XWX[3,3] <- -sum( mu*(1-mu)*Psi^2)
XWX[1,3] <-XWX[3,1]<- -sum( mu*(1-mu)*Psi*ESP.eta.der.tau)
XWX[2,3] <-XWX[3,2]<- -sum( mu*(1-mu)*Psi*ESP.eta.der.alpha)

Cov.gwesp<-solve(-XWX)
return(list(Cov.gwesp=Cov.gwesp,XWX=XWX))
}#end




compute.delta.Zij<-function(Y,n){
Ny<-length(Y)
deltaZij<-array(0,dim=c(Ny,Ny,n-2))
for(i in 1:(Ny-1)){
cat(i,"th iteration \n")
   for(j in (i+1):Ny){
     Y11<-Y00<-Y01<-Y10<-Y
     Y11[i]<-1;Y11[j]<-1
     Y10[i]<-1;Y10[j]<-0
     Y01[i]<-0;Y01[j]<-1
     Y00[i]<-0;Y00[j]<-0
     Y11<-as.network(Psi.matrix(Y11,n),directed=F)
     Y10<-as.network(Psi.matrix(Y10,n),directed=F)
     Y01<-as.network(Psi.matrix(Y01,n),directed=F)
     Y00<-as.network(Psi.matrix(Y00,n),directed=F)
     deltaZij[j,i,]<-deltaZij[i,j,]<-summary(Y11~esp(1:(n-2)))+summary(Y00~esp(1:(n-2)))-summary(Y10~esp(1:(n-2)))-summary(Y01~esp(1:(n-2)))
     #compute E(Y_i=1,Y_j=1)
     }#end for(j in  (i+1):Ny
   }#end for(i in 1:(Ny-1)){
return(deltaZij)
}#end compute.delta.Zij


compute.approx.Cov<-function(deltaZij,gwesp.tau,gwesp.alpha,mu,n){
Ny<-n*(n-1)/2
E11<-CovY<-matrix(0,Ny,Ny)
eta<-gwesp(gwesp.tau,gwesp.alpha,n)
for(i in 1:(Ny-1)){
   for(j in (i+1):Ny){
     alpha<-exp(deltaZij[i,j,]%*%eta)
     mux<-mu[i];muy<-mu[j]
     coeff<-c(alpha-1,-(1+(alpha-1)*(mux+muy)),alpha*mux*muy);
     coeff<-coeff[3:1]
     hilf <- Re(polyroot(coeff));
     s1<-min(mux,muy);s2<-max(0,mux+muy-1);
     cond<-(hilf<=s1 & hilf>=s2);
      #print(cond);print(hilf);print(s1);print(s2)
      #if(cond[1]){
      E11[i,j]<-E11[j,i]<-hilf[1]
      #cat("first root\n")
      #}else{
      #cond[2] 
      #E11[i,j]<-E11[j,i]<-hilf[2]
      #cat("second root\n")}
      CovY[i,j]<-CovY[j,i]<-E11[i,j]-mux*muy
   }#end for(j in  (i+1):Ny
}#end for(i in 1:(Ny-1)){
diag(CovY)<-mu*(1-mu)
return(list(CovY=CovY,E11=E11))
}


output<-function(Y){
n<-length(Y)
out<-NULL
for(i in 1:n){
out<-paste(out," & ",Y[i],sep="")
}
out
}