R/simulate.R

In ComputationalIntelligenceGroup/hbp-dendrite-arborization-simulation: basal arborization reconstruction and simulation

#' Given a terminal node, simulate the grow of the dendrite one step more
#'
#' Given a terminal node, simulate the grow of the dendrite one step more
#'
#' @param model is a bn model generated by bnlearn
#' @param node_info is a data.frame containing the terminal nodes that must grow
#'
#' @return a data.frame with the spherical coordinates of the new node with respect to the terminal node introduced to the function
#'
#' @example
#' path<-"/home/universidad/Documents/neuron_data/datos/All"
#' data<-get_features_files(path,60)
#' data<-compute_clusters(data)
#' simulation_model<-compute_simulation_models(data)
#' simulation<-simulate_continuation(simulation_model$continue,data[data$node_num_descendant==1])
simulate_continuation<-function(model,node_info)
{
  #Get dataset of nodes to grow and get topological order of the variables according to the BN structure
  node_info<-data.frame(node_info)
  node_order_model<-node.ordering(model$structure)
  first_desc<-c("desc_azimuth_angle","desc_elevation_angle","desc_length")
  first_desc<-node_order_model[node_order_model%in%first_desc]

  #Save position of simulated nodes
  simulated_nodes<-matrix(0,ncol=length(first_desc),nrow=nrow(node_info))

  #For each feature to simulate
  for(i in 1:length(first_desc))
  {
    #Initiate mean and sd matrices to sample
    mean_matrix<-matrix(0,ncol=1,nrow=nrow(node_info))
    sd_matrix<-matrix(0,ncol=1,nrow=nrow(node_info))

    #Get the information of the parents of the node (are discrete and their names)
    node<-first_desc[i]
    parent_names<-parents(model$structure,node)
    is_discrete_p<-sapply(parent_names,function(x){class(model$params[[x]])=="bn.fit.dnode"})
    discrete_p_names<-parent_names[is_discrete_p]
    continuous_p_names<-parent_names[!is_discrete_p]

    #If the node is Gaussian, that is, all its parents are Gaussian or dont have any parent
    if(class(model$params[[node]])=="bn.fit.gnode"){
      mean_matrix[,1]<-coefficients(model$params[[node]])%*%t(cbind(1,node_info[,colnames(node_info)%in%continuous_p_names]))
      sd_matrix[,1]<-model$params[[node]]$sd
    }

    #If the node have discrete parents
    if(class(model$params[[node]])=="bn.fit.cgnode"){

      #Get the configuration of the parents to select the correct coefficients for the normal distribution
      pos_parents<-lapply(discrete_p_names,function(x){match(node_info[,x],as.numeric(model$params[[node]]$dlevels[[x]]))})
      num_poss_parents<-sapply(model$params[[node]]$dlevels,length)#Numbre of possible values for each parent

      #The position in dlevels is a+num_poss_parents[1]*(b-1)+num_poss_parents[2]*num_poss_parents[1]*(c-1)+...
      match_idx<-pos_parents[[1]]#If there is only 1 parent then the position is a
      if(length(num_poss_parents)>1)#If there is more than one parent the position is a combination
      {
        prev_parent_num_combinations<-1
        for(j in 2:length(num_poss_parents))
        {
          prev_parent_num_combinations<-prev_parent_num_combinations*num_poss_parents[j-1]
          match_idx<-match_idx+prev_parent_num_combinations*(pos_parents[[j]]-1)
        }
      }

      #If all the parent nodes are discrete we only have the mean of the node
      if(length(parent_names)==length(discrete_p_names))
      {
        mean_matrix<-model$params[[node]]$coefficients[match_idx]
      }else{#If there is at least one gaussian parent the mean is a linear combination
        mean_matrix<-colSums(model$params[[node]]$coefficients[,match_idx]*t(cbind(1,node_info[,colnames(node_info)%in%continuous_p_names])))
      }
        sd_matrix<-model$params[[node]]$sd[match_idx]
    }

    #Simulate values for the node
      #simulated_nodes[,i]<-rnorm(mean_matrix,mean=mean_matrix,sd=sd_matrix)
      if(node=="desc_length")
      {
        simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,upper=5.75)
      }else if(node=="desc_azimuth_angle")
      {
        simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,lower=-pi,upper=pi)
      }else{
        simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,lower=0,upper=pi)
      }

    #Save the simulation to use them for the next simulations
      node_info[[first_desc[i]]]<-simulated_nodes[,i]
      node_info<-node_info[,order(names(node_info))]
  }
  simulated_nodes<-data.frame(simulated_nodes)
  colnames(simulated_nodes)<-first_desc
  return(simulated_nodes)
}


#' Given a terminal node, simulate the bifurcation of the dendrite
#'
#' Given a terminal node, simulate the bifurcation of the dendrite
#'
#' @param model is a bn model generated by bnlearn
#' @param node_info is a data.frame containing the terminal nodes that must grow
#'
#' @return a data.frame with the spherical coordinates of the two new nodes with respect to the terminal node introduced to the function
#'
#' @example
#' path<-"/home/universidad/Documents/neuron_data/datos/All"
#' data<-get_features_files(path,60)
#' data<-compute_clusters(data)
#' simulation_model<-compute_simulation_models(data)
#' simulation<-simulate_bifurcation(simulation_model$bifurcation,data[data$node_num_descendant==2,])
simulate_bifurcation<-function(model,node_info)
{
  node_info<-data.frame(node_info)
  node_order_model<-node.ordering(model$structure)
  first_desc<-c("desc_azimuth_angle","desc_elevation_angle","desc_length","desc_azimuth_angle2","desc_elevation_angle2","desc_length2")
  first_desc<-node_order_model[node_order_model%in%first_desc]

  #Simulate the discrete nodes
  disc_sim<-c("desc_longer","dendrite_diameter")
  disc_sim<-node_order_model[node_order_model%in%disc_sim]
  for(i in 1:length(disc_sim))
  {
  node<-disc_sim[i]
  parents_longer<-parents(model$structure,node)
  if(length(parents_longer)==0)
  {
    node_info[[node]]<-c(rMultinom(matrix(model$params[[node]]$prob,nrow=1),nrow(node_info)))
  }else{
    #Get the configuration of the parents to select the correct probability
    pos_parents<-lapply(parents_longer,function(x){match(node_info[,x],as.numeric(attributes(model$params[[node]]$prob)$dimnames[[x]]))})
    pos_parents<-c(list(rep(1,length(pos_parents[[1]]))),pos_parents)
    num_poss_parents<-sapply(attributes(model$params[[node]]$prob)$dimnames,length)#Number of possible values for each parent

    #The position in dlevels is a+num_poss_parents[1]*(b-1)+num_poss_parents[2]*num_poss_parents[1]*(c-1)+...
    match_idx<-pos_parents[[1]]#If there is only 1 parent then the position is a
    if(length(num_poss_parents)>1)#If there is more than one parent the position is a combination
    {
      prev_parent_num_combinations<-1
      for(j in 2:length(num_poss_parents))
      {
        prev_parent_num_combinations<-prev_parent_num_combinations*num_poss_parents[j-1]
        match_idx<-match_idx+prev_parent_num_combinations*(pos_parents[[j]]-1)
      }
    }
    match_matrix<-sapply(match_idx,function(x){x:(x+num_poss_parents[1]-1)})
    node_info[[node]]<-c(rMultinom(probs=matrix(model$params[[node]]$prob[c(match_matrix)],ncol=num_poss_parents[1],byrow=T),m=1))
  }
  node_info<-node_info[,order(names(node_info))]
  }
  #Save position of simulated nodes
  simulated_nodes<-matrix(0,ncol=length(first_desc),nrow=nrow(node_info))
  node_info<-node_info[,order(names(node_info))]

  #For each feature to simulate
  for(i in 1:length(first_desc))
  {
    #Initiate mean and sd matrices to save
    mean_matrix<-matrix(0,ncol=1,nrow=nrow(node_info))
    sd_matrix<-matrix(0,ncol=1,nrow=nrow(node_info))

    #Get the information of the parents of the node (are discrete and their names)
    node<-first_desc[i]
    parent_names<-parents(model$structure,node)
    is_discrete_p<-sapply(parent_names,function(x){class(model$params[[x]])=="bn.fit.dnode"})
    discrete_p_names<-parent_names[is_discrete_p]
    continuous_p_names<-parent_names[!is_discrete_p]

    #If the node is Gaussian, that is, all its parents are Gaussian or dont have any parent
    if(class(model$params[[node]])=="bn.fit.gnode"){
      mean_matrix[,1]<-coefficients(model$params[[node]])%*%t(cbind(1,node_info[,colnames(node_info)%in%continuous_p_names]))
      sd_matrix[,1]<-model$params[[node]]$sd
    }

    #If the node have discrete parents
    if(class(model$params[[node]])=="bn.fit.cgnode"){

      #Get the configuration of the parents to select the correct coefficients for the normal distribution
      pos_parents<-lapply(discrete_p_names,function(x){match(node_info[,x],as.numeric(model$params[[node]]$dlevels[[x]]))})
      num_poss_parents<-sapply(model$params[[node]]$dlevels,length)#Numbre of possible values for each parent

      #The position in dlevels is a+num_poss_parents[1]*(b-1)+num_poss_parents[2]*num_poss_parents[1]*(c-1)+...
      match_idx<-pos_parents[[1]]#If there is only 1 parent then the position is a
      if(length(num_poss_parents)>1)#If there is more than one parent the position is a combination
      {
        prev_parent_num_combinations<-1
        for(j in 2:length(num_poss_parents))
        {
          prev_parent_num_combinations<-prev_parent_num_combinations*num_poss_parents[j-1]
          match_idx<-match_idx+prev_parent_num_combinations*(pos_parents[[j]]-1)
        }
      }

      if(length(parent_names)==length(discrete_p_names))#If all the parent nodes are discrete we only have the mean of the node
      {
        mean_matrix<-model$params[[node]]$coefficients[match_idx]
      }else{#If there is at least one gaussian parent the mean is a linear combination
        mean_matrix<-colSums(model$params[[node]]$coefficients[,match_idx]*t(cbind(1,node_info[,colnames(node_info)%in%continuous_p_names])))
      }
      sd_matrix<-model$params[[node]]$sd[match_idx]
    }

    #simulated_nodes[,i]<-rnorm(mean_matrix,mean=mean_matrix,sd=sd_matrix)
    if(grepl("desc_length*",node)){
      simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,upper=5.75)
    }else if(grepl("desc_azimuth_angle*",node)){
      simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,lower=-pi,upper=pi)
    }else{
      simulated_nodes[,i]<-rtnorm(length(mean_matrix),mean=mean_matrix,sd=sd_matrix,lower=0,upper=pi)
    }
    node_info[[first_desc[i]]]<-simulated_nodes[,i]
    node_info<-node_info[,order(names(node_info))]
  }
  simulated_nodes<-data.frame(simulated_nodes)
  colnames(simulated_nodes)<-first_desc
  simulated_nodes<-simulated_nodes[,order(names(simulated_nodes))]
  return(simulated_nodes)
}




#' Given the coordinates of a root node, simulate the coordinates of the first branch
#'
#' Given the coordinates of a root node, simulate the coordinates of the first branch
#'
#' @param model is a bn model generated by bnlearn
#' @param evidence is a list containing the coordinates of the root nodes that must grow
#'
#' @return a list with the cartesian coordinates of the new nodes with respect to the coordinates of the root node introduced to the function
#'
#' @example
#' path<-"/home/universidad/Documents/neuron_data/datos/All"
#' data <-extract_first_node_coordinates(path,eps)
#' model <- train_first_node_coordinates(data)
#' path<-"/home/universidad/Documents/neuron_data/datos/All/neuron.DAT"
#' neuron<-neuro_converter(path,eps=eps)
#' root_nodes_position[[i]][1] = neuron$neurons$neurites[[1]][i,]$tree$children[[1]]$root$x
#' root_nodes_position[[i]][2] = neuron$neurons$neurites[[1]][i,]$tree$children[[1]]$root$y
#' root_nodes_position[[i]][3] = neuron$neurons$neurites[[1]][i,]$tree$children[[1]]$root$z
#' first_bifurcation <- simulate_first_bifurcation(model,list(x2.x = root_nodes_position[[i]][1],y2.y = root_nodes_position[[i]][2],z2.z = root_nodes_position[[i]][3]))
simulate_first_bifurcation <- function(model,evidence){

  model_fit <-model$params
  model_2_rbmn <- bnfit2nbn(model_fit)
  model_2_rbmn_mn <- nbn2mn(model_2_rbmn)

  # This is the names vector for my nodes:
  names_variables <- names(model$structure$nodes)
  # Names of given variables
  obs_names <- c("x2.x","y2.y","z2.z")
  # Names of variables to be predicted
  pred_names <- setdiff(names_variables, obs_names)
  evidence <- as.data.frame(evidence)
  list_prev_models <- list()
  for(i in 1:length(evidence[,1])){
    obs_val <- evidence[i,]
    predicted_values <- condi4joint(model_2_rbmn_mn, par = pred_names, pour = obs_names, x2 = obs_val)
    generated_data<-predicted_values$mu
    list_prev_models[[i]] <- generated_data
  }

  return(list_prev_models)




}


simulate_neo_tortosity <- function(tortosity_model,eps60pred,eps60datapoint,num_of_descendant,simorcont){

    model <- tortosity_model
    prev_model<-eps60pred
    node_info<-eps60datapoint

    #prev_model$desc_length <- exp(prev_model$desc_length)

    model_fit <- model$params
    model_2_rbmn <- bnfit2nbn(model_fit)
    model_2_rbmn_mn <- nbn2mn(model_2_rbmn)

    # This is the names vector for my nodes:
    names_variables <- names(model$structure$nodes)

    if(simorcont==1){
      #names_next_branch <- names(eps60pred)
      next_branch <- eps60pred
    }
    else{
      eps60pred <- eps60pred[order(eps60pred[,4]),]
      eps60pred$desc_azimuth_angle2 <- 0
      eps60pred$desc_elevation_angle2<-0
      eps60pred$desc_length2 <- 0
      mergedrows<-c()
      for(i in 1:(length(eps60pred[,4])/2)){
        i1 <- (i-1)*2+1
        i2 <- i1+1
        next_branch <- eps60pred[i1,]
        next_branch2 <- eps60pred[i2,]
        eps60pred[i1,c(5,6,7)]<- next_branch2[c(1,2,3)]
        mergedrows<- c(mergedrows,i2)
      }

      eps60pred <- eps60pred[-mergedrows,]
      next_branch <- eps60pred
    }

    datapointnames <- names(eps60datapoint)

    # Names of given variables
    obs_names <-  c( datapointnames[1:2],names(next_branch[,-4]), datapointnames[3:9])#specify observation names here

    # Names of variables to be predicted
    pred_names <- setdiff(names_variables, obs_names)
    prev_model_final <- prev_model[1,]
    list_prev_models <- list()
    for(i in 1:length(eps60datapoint[[1]])){
        obs_val <- unlist(c(eps60datapoint[i,1:2],next_branch[i,-4],eps60datapoint[i,3:9]))

      # Then predict all the unknown variables
      predicted_values <- condi4joint(model_2_rbmn_mn, par = pred_names, pour = obs_names, x2 = obs_val)
      generated_data<-predicted_values$mu

      #PROPERLY SIMULATE

      if(simorcont == 1){
      N_nodes_branch <-generated_data[13]
      }
      else{
      N_nodes_branch <-generated_data[14]
      generated_data<- generated_data[-1]
      }

      if(N_nodes_branch <0){
        N_nodes_branch <- abs(N_nodes_branch)
      }


      if(N_nodes_branch >1){

        N<-ceiling(N_nodes_branch/4+1)
        #N <- 2*ceiling(log(N_nodes_branch)+1)
        #N <- ceiling(log(N_nodes_branch)+1)
        angle_distr <- rnorm(N,mean = eps60datapoint$azimuth_angle[i], abs(generated_data[17]))*0
        comp_length_distr <- rnorm(N,mean = eps60datapoint$elevation_angle[i], abs(generated_data[19]))*0#because negative values sometimes appear


#variables for knowledge extraction of the problem. TESTING
      testing_env = 1
      if(testing_env){
        length_distr <- rnorm(N,mean = generated_data[14] ,sd = abs(generated_data[15]))
        #angle_distr <- rnorm(N,mean = 0, sd = abs(generated_data[17]/5))
        comp_length_distr <- rnorm(N,mean = pi/2, sd = abs(generated_data[19]/3))#
        length_distr <- abs(length_distr)
        #angle_distr <- angle_distr*0
        angle_distr[1]<-eps60datapoint$azimuth_angle[i]
        comp_length_distr[1] <- pi/2
        for(k in 2:N){
          shift <- sum(angle_distr[2:N])
          angle_distr[k] <- rnorm(1,mean = -shift,sd = abs(generated_data[17])/1)
        }

      }

#END OF TESTING ENVIRONMENT

        ids = prev_model[i,4]
        datapoints<-list()
        prev_model_2 <- matrix(0,nrow = (N),ncol = 4)
        prev_model_2 <- as.data.frame(prev_model_2)
        setnames(prev_model_2,names(prev_model))

        for(j in 1:N){
          final_length <- length_distr[j]
          final_azimuth <- angle_distr[j]
          final_elevation <- comp_length_distr[j]
          final_id <- ids[j]
          prev_model_2[j,] <- c(final_azimuth,final_elevation,final_length,final_id)
        }


        prev_model_final <- rbind(prev_model_final,prev_model_2)
      }

      else{
        prev_model_2 <- prev_model[i,]
        prev_model_final <- rbind(prev_model_final,prev_model[i,])

      }

      list_prev_models <- append(list_prev_models,list(prev_model_2))

    }

    prev_model_final <- prev_model_final[-1,]

    return(list_prev_models)

}