R/deprecated/lolog_change_stats.R

In duncan-clark/Blolog: Bayesian Inference for LOLOG Network Models

#Given a network and a permuation this function calculates the LOLOG change stats
#after each edge has been considered in the edge permutation scheme

#' @export

lolog_change_stats <- function(net,         #the network of interest
                               s,           #The permutations of interest
                               formula_rhs #rhs formula of model
                               ){
  
  #define network size
  n <- network.size(net)
  e_obs <- length(net$mel)
  if(get.network.attribute(net,"directed")){e<-n*(n-1)
  }else{e <- (n)*(n-1)/2}
  
  if(length(s)!= e){stop("permutation is the wrong length")}
  
  #make edge list that we can permute
  edges <- combn(seq(1,n),2)
  if(get.network.attribute(net,"directed")){edges <- as.list(cbind(as.data.frame(edges),as.data.frame(edges[c(2,1),])))
  }else{edges <- as.list(as.data.frame(edges))}
  
  #returns on and off change stats for two thetas
  #empty <- network(x = matrix(rep(0,n**2),nrow = n),directed =y$isDirected())
  #This way makes a model that retains vertex covariates.
  empty <- net
  net <- lolog::as.BinaryNet(net)
  delete.edges(empty,seq(1,e_obs))    
    
  model <- createCppModel(as.formula(paste("empty ~",formula_rhs,sep = "")))
  order <- unique(unlist(edges[s]))
  edges <- edges[s]
    
  change_on <- list(rep(NULL,length(edges)))
  change_off <- lapply(1:length(edges),function(x){rep(0,length(model$statistics))})
    
  for(t in 1:length(edges)){
    tail <- edges[[t]][1]
    head <- edges[[t]][2]
      
    old_stats <- model$statistics()
    #need -1 here to account for different R and C++ indexing
    model$dyadUpdate((tail-1),(head-1),order,head)
    change_on[[t]] <- model$statistics() - old_stats
    
    tmp <- model$getNetwork()
      
    #don't need -1 here to account for different R and C++ indexing since the function is exposed to R properly
    tmp$setDyads(tail,head,net$getDyads(tail,head)*1)
      
    #plot(as.network(tmp),label = seq(1,n),pad =2,label.pos = 1)
      
    model$setNetwork(tmp)
    #model$calculate - don't need to calculate since already updated the dyads
    #print(t/length(edges))
  }

  return(change_on)
}