##set random seed
set.seed(1)

##some functions
source("dist.r")           

##library mva provides cdmscale, used to get starting values
library(mva)         

##read in Padgett and Ansell's (1993) Florentine family marriage data
Y_as.matrix(read.table("ff.m.dat"))

##take out unconnected actors
zro_(1:dim(Y)[1])[ Y%*%rep(1,dim(Y)[1])==0 ] 
if(length(zro)>0) {Y_Y[-zro,-zro]}


n_dim(Y)[1]    #number of nodes
k_2            #dimension of latent space


## starting values of Z and alpha
D_dst(Y)
Z_cmdscale(D, 2)
Z[,1]_Z[,1]-mean(Z[,1])
Z[,2]_Z[,2]-mean(Z[,2])

alpha_2

## now find a good starting point for the MCMC
avZ_c(alpha,c(Z))
#simulated annealing 
avZ_optim(avZ,mlpY,Y=Y,method="SANN")$par
#BFGS 
avZ_optim(avZ,mlpY,Y=Y,method="BFGS")$par

avZ_avZ*2/(avZ[1])  #MLE too extreme, keep shape but rescale
alpha_avZ[1]
Z_Z.mle_matrix(avZ[-1],nrow=n,ncol=k)


##see if this Z represents Y
(-1)^(1-Y) * (alpha+lpz.dist(Z) ) >=0

## MCMC params

adelta_.5          #params for
zdelta_.2          #proposal distribution

nscan_10^6         #number of scans   
odens_10^3         #save output every odens step

aca_acz_0          #keep track of acceptance rates
Alpha_alpha        #keep track of alpha
Lik_lpY(Y,Z,alpha) #keep track of alpha and likelihood 

Z.post_list()    #keep track of positions
for(i in 1:k){ Z.post[[i]]_t(Z[,i]) }


## MCMC

for(ns in 1:nscan){

tmp_Z.up(Y,Z,alpha,zdelta)                     #update z's
if(tmp$Z!=Z){ 
              acz_acz+1/odens
              Z_proc.crr(tmp$Z,Z.mle)  
             }

tmp_alpha.up(Y,Z,alpha,adelta,a.a=2,a.b=1)  #update alpha
if(tmp$alpha!=alpha) { 
                       aca_aca+1/odens
                       alpha_tmp$alpha 
                       lik_tmp$lik
                      }

if( ns%%odens==0 ){                            #output
                     Alpha_c(Alpha,alpha)
                     cat(ns,aca,acz,alpha,lik,"\n")
                     Lik_c(Lik,lik)
                     acz_aca_0 
                     for(i in 1:k){ Z.post[[i]]_rbind(Z.post[[i]],t(Z[,i])) } 
                   }

                  }

Post_list(Z=Z.post,Alpha=Alpha,Lik=Lik)
dput(Post,"ff.m.post") #output to a file

##plot results, for k=2
Zp_Post$Z
Lik_Post$Lik
Alpha_Post$Alpha


##posterior mean 
Z.pm_Z.mle
for(i in 1:k){  Z.pm[,i]_rep(1/dim(Zp[[i]])[1],dim(Zp[[i]])[1])%*%Zp[[i]] }
Z.mle_Z.mle*sum(diag((Z.pm%*%t(Z.mle))))/sum(diag(( t(Z.mle)%*%Z.mle)))
Z.mle_proc.crr(Z.mle,Z.pm)

##plot likelihood and alpha
par(mfrow=c(1,2))
plot(Lik,type="l")
plot(Alpha,type="l")

##plot positions and confidence sets
par(mfrow=c(1,1))
par(mar=c(3,3,1,1))
par(mgp=c(2,1,0))

##set up color scheme based on a scaled Z.mle


r_atan2(Z.mle[,2],Z.mle[,1])
r_r+abs(min(r))
r_r/max(r)
g_1-r
b_(Z.mle[,2]^2+Z.mle[,1]^2)
b_b/max(b)

plot(Z.mle[,1],Z.mle[,2],xaxt="n",yaxt="n",xlab="",ylab="",type="n", 
                          xlim=range(Zp[[1]]),ylim=range(Zp[[2]]))

for(i in 1:nrow(Y)){
for(j in 1:nrow(Y)){ lines(  Z.mle[c(i,j),1] , Z.mle[c(i,j),2] ,lty=Y[i,j]) }}

text(Z.mle[,1]*1.1,Z.mle[,2]*1.1,labels(Y)[[1]])

for(i in 1:n) {
      points( Zp[[1]][,i],Zp[[2]][,i],pch=46,cex=.75,
      col=rgb(r[i],g[i],b[i]) ) }