#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.