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

dst_function(Y,infd=F){
#calculates distances between nodes based on path length
#returns d=g if nodes are not connected

g_dim(Y)[1]
Dst_Yr_Y
Dst_Y*(Y==1) + g*(Y==0)
for(r in 2:(g-1)) {
Yr_Yr%*%Y
Dst_Dst+(r-g)*( Yr>0 & Dst==g )
                  }

if(infd==T){for(i in 1:g){ for(j in 1:g) { if( Dst[i,j]==g ){ Dst[i,j]_Inf} }
  }}
diag(Dst)_0

Dst  }


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


proc.crr_function(Z,Z0,k=2){
#Procrustes transform; gives rotation,reflection,trranslation 
#of Z closest to Z0

for( i in 1:k) { Z[,i]_Z[,i]-mean(Z[,i]) + mean(Z0[,i]) } #translation

A_t(Z)%*%(  Z0%*%t(Z0)  )%*%Z
eA_eigen(A,symmetric=T)
Ahalf_eA$vec[,1:k]%*%diag(sqrt(eA$val[1:k]))%*%t(eA$vec[,1:k])

t(t(Z0)%*%Z%*%solve(Ahalf)%*%t(Z)) }



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

lpY_function(Y,Z=0,alpha){
#log probability of the graph 

#Y is sociomatrix
#Z are latent positions
#alpha is intercept


lpZ_lpz.dist(Z)    #the function of Z that is in the linear predictor

lpg_alpha+lpZ

diag(lpg)_0
sum( Y*lpg - log( 1+exp(lpg) ) ) + dim(Y)[1]*log(2) 
}

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

mlpY_function(avZ,Y){

#minus log prob of graph, with Z in vector form, 
#to be used in by optim or nlm

alpha_avZ[1]
Z_matrix( avZ[-1],nrow=n,ncol=k)
lpZ_lpz.dist(Z)    #the function of Z that is in the linear predictor
lpg_alpha+lpZ
diag(lpg)_0
-(  sum( Y*lpg - log( 1+exp(lpg) ) ) + dim(Y)[1]*log(2)  )
}


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

lpz.dist_function(Z){
##gives the negative distance between nodes
ZtZ_Z%*%t(Z)

mg_as.matrix(diag(ZtZ))%*%rep(1,length(Z[,1])) #distances
mg_mg+t(mg)
d_sqrt((mg-2*ZtZ))
-d             }

########
Z.up_function(Y,Z,alpha,zdelta,mu.z=0,sd.z=10){
##update Z
Znew_Z+matrix(rnorm(k*n,0,zdelta),nrow=n,ncol=k)

lnew_lpY(Y,Znew,alpha)
lold_lpY(Y,Z,alpha)

hr_lnew-lold+sum( dnorm(Znew,mu.z,sd.z,log=T) )-sum( dnorm(Z,mu.z,sd.z,log=T) )
if( runif(1)>exp(hr) ){Znew_Z
                       lnew_lold }

list(Z=Znew,lik=lnew)}

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

alpha.up_function(Y,Z,alpha,adelta,a.a=1,a.b=1){
##update alpha
alphanew_abs(alpha+runif(1,-adelta,adelta) )

lnew_lpY(Y,Z,alphanew)
lold_lpY(Y,Z,alpha)

hr_exp( lnew-lold )*
      ( dgamma(alphanew,a.a,scale=a.b)/dgamma(alpha,a.a,scale=a.b) ) 

if(runif(1)>hr){ alphanew_alpha
                 lnew_lold }

list(alpha=alphanew,lik=lnew)  }

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