Nothing
.one_cycle <- function(cv_data ,
r ,
s ,
estimated_fitness_start = NULL ,
estimated_PA_start = NULL ,
alpha_start = NULL ,
stop_cond = 10^-6 ,
print_out = FALSE ,
cv_deg_thresh = c(1,10) ,
normal_start_f = TRUE ,
weight_f = 0 ,
...) {
oopts <- options(scipen = 999)
on.exit(options(oopts))
FitMultinomial <- function(true,dat){
true[true == 0] <- 1
return(sum(dat*log10(true)))
}
rate_PAFit <- s
if (is.null(estimated_PA_start)) {
alpha_optimal <- rep(-1,length(cv_deg_thresh))
}
alpha_each <- matrix(0,nrow = length(cv_deg_thresh), ncol = length(rate_PAFit))
colnames(alpha_each) <- rate_PAFit
estimated_fitness <- vector(mode = "list", length = length(cv_deg_thresh))
#estimated_PA <- NULL
max_val <- rep(-Inf,length(cv_deg_thresh))
s_optimal <- rep(0,length(cv_deg_thresh))
for (j in 1:length(rate_PAFit)) {
#print(paste0("Processing case ",count, " of ",total))
if (is.null(estimated_PA_start)) {
if (!is.null(estimated_fitness_start)) {
result_PAFit <- PAFit(cv_data$stats , mode_f = "Log_linear" , s = rate_PAFit[j],
auto_stop = TRUE , stop_cond = stop_cond , normalized_f = FALSE ,
normal_start_f = normal_start_f,
weight_power = weight_f, start_f = estimated_fitness_start,
...);
}
else
result_PAFit <- PAFit(cv_data$stats , mode_f = "Log_linear" , s = rate_PAFit[j],
auto_stop = TRUE , stop_cond = stop_cond , normalized_f = FALSE ,
normal_start_f = normal_start_f,
weight_power = weight_f,
...);
}
else {
#print("inside here");
if (!is.null(estimated_fitness)) {
result_PAFit <- PAFit(cv_data$stats , true_A = estimated_PA_start,
s = rate_PAFit[j] , only_f = TRUE,
auto_stop = TRUE ,
stop_cond = stop_cond , normalized_f = FALSE ,
normal_start_f = normal_start_f,
weight_power = weight_f, start_f = estimated_fitness_start,
...);
#plot(result_PAFit$center_k,result_PAFit$theta,log = "xy", main = "Inside loop");
#print(paste0("Max of theta inside one cycle: ",max(estimated_PA_start)))
}
else
result_PAFit <- PAFit(cv_data$stats , true_A = estimated_PA_start,
s = rate_PAFit[j] , only_f = TRUE,
auto_stop = TRUE ,
stop_cond = stop_cond , normalized_f = FALSE ,
normal_start_f = normal_start_f,
weight_power = weight_f,
...);
}
if (is.null(estimated_PA_start)) {
alpha_temp <- result_PAFit$alpha
#print(alpha_temp)
if (TRUE == result_PAFit$diverge_zero) {
alpha_each[,j] <- -Inf
}
else {
for (dd in 1:length(cv_deg_thresh)) {
alpha_each[dd,j] <- 0
chosen_node_big <- names(cv_data$stats$z_j[cv_data$stats$z_j >= cv_deg_thresh[dd]])
for (k in 1:length(cv_data$m_each))
if (cv_data$m_each[k] != 0) {
#PA <- result_PAFit$A[cv_data$deg_each[k,] + 1]
#PA[is.na(PA)] <- result_PAFit$A[length(result_PAFit$A)]
chosen_node <- chosen_node_big[cv_data$deg_each[k,chosen_node_big] != 0]
PA <- cv_data$deg_each[k,chosen_node]^alpha_temp
#print(PA)
if (sum(PA == 0) > 0)
print("This should not happen")
#PA[PA == 0] <- 1
fitness <- rep(1,dim(cv_data$deg_each[,chosen_node, drop = FALSE])[2])
names(fitness) <- colnames(cv_data$deg_each[,chosen_node, drop = FALSE])
fitness[chosen_node] <- result_PAFit$f[chosen_node]
prob_PAFit <- PA * fitness
prob_PAFit <- prob_PAFit / sum(prob_PAFit,na.rm = TRUE)
#temp_sapply <- sapply(prob_PAFit,is.na)
#print(prob_PAFit)
#print(temp_sapply)
if (length(prob_PAFit) > 0) {
prob_PAFit[sapply(prob_PAFit,is.na)] <- 0
alpha_each[dd,j] <- alpha_each[dd, j] +
FitMultinomial(true = as.vector(prob_PAFit),
dat = as.vector(cv_data$prob_em_each[k,chosen_node] *
cv_data$m_each[k]))
}
}
}
}
for (dd in 1:length(cv_deg_thresh))
if (alpha_each[dd,j] > max_val[dd]) {
estimated_fitness[[dd]] <- result_PAFit$f[as.character(as.numeric(cv_data$stats$f_position))]
if (sum(is.na(estimated_fitness[[dd]])) > 0) {
if (sum(is.na(result_PAFit$f)) > 0) {print("both result_PAFit$f and estimated_fitness[dd] have NA")}
else {print("only estimated_fitness[dd] have NA"); f_pos <- cv_data$stats$f_position;
save(result_PAFit,cv_data,f_pos,estimated_fitness, file = "net_stat_bug_detailed.Rdata")}
}
s_optimal[dd] <- rate_PAFit[j]
max_val[dd] <- alpha_each[dd,j]
alpha_optimal[dd] <- result_PAFit$alpha
}
} else {
### estimated_PA is provided
if (TRUE == result_PAFit$diverge_zero) {
alpha_each[,j] <- -Inf
} else {
for (dd in 1:length(cv_deg_thresh)) {
alpha_each[dd,j] <- 0
chosen_node_big <- names(cv_data$stats$z_j[cv_data$stats$z_j >= cv_deg_thresh[dd]])
for (k in 1:length(cv_data$m_each))
if (cv_data$m_each[k] != 0) {
chosen_node <- chosen_node_big[!is.na(result_PAFit$A[cv_data$deg_each[k,chosen_node_big] + 1])]
chosen_node <- chosen_node[result_PAFit$A[cv_data$deg_each[k,chosen_node] + 1] != 0]
#PA <- result_PAFit$A[cv_data$deg_each[k,chosen_node] + 1]
#PA <- result_PAFit$A[cv_data$deg_each[k,chosen_node] + 1]
PA <- (cv_data$deg_each[k,chosen_node])^alpha_start
if (sum(is.na(PA)) > 0) {
#PA[is.na(PA)] <- max(estimated_PA_start)
#PA[is.na(PA)] <- cv_data$deg_each[k,chosen_node][is.na(PA)]^alpha_start
print("It should not happen")
}
if (sum(PA == 0) > 0)
PA[PA == 0] <- 1
fitness <- rep(1,dim(cv_data$deg_each[,chosen_node, drop = FALSE])[2])
names(fitness) <- colnames(cv_data$deg_each[,chosen_node, drop = FALSE])
fitness[chosen_node] <- result_PAFit$f[chosen_node]
prob_PAFit <- PA * fitness
prob_PAFit <- prob_PAFit / sum(prob_PAFit,na.rm = TRUE)
#temp_sapply <- sapply(prob_PAFit,is.na)
#print(prob_PAFit)
#print(temp_sapply)
if (length(prob_PAFit) > 0) {
prob_PAFit[sapply(prob_PAFit,is.na)] <- 0
alpha_each[dd,j] <- alpha_each[dd, j] +
FitMultinomial(true = as.vector(prob_PAFit),
dat = as.vector(cv_data$prob_em_each[k,chosen_node] *
cv_data$m_each[k]))
}
}
}
}
for (dd in 1:length(cv_deg_thresh))
if (alpha_each[dd,j] > max_val[dd]) {
estimated_fitness[[dd]] <- result_PAFit$f[as.character(as.numeric(cv_data$stats$f_position))]
if (sum(is.na(estimated_fitness[[dd]])) > 0) {
if (sum(is.na(result_PAFit$f)) > 0) {print("both result_PAFit$f and estimated_fitness[dd] have NA")}
else {print("only estimated_fitness[dd] have NA")}
}
s_optimal[dd] <- rate_PAFit[j]
max_val[dd] <- alpha_each[dd,j]
#alpha_optimal[dd] <- result_PAFit$alpha
}
}
}
ok_index <- 1:length(cv_deg_thresh)
s_optimal_final <- mean(s_optimal[ok_index])
#s_optimal_final <- 10^(1/length(s_optimal[ok_index])* log10(prod(s_optimal[ok_index])))
#print("NA in estimated_fitness[[ok_index[1]]]:")
#print(sum(is.na(estimated_fitness[[ok_index[1]]])))
if (is.null(estimated_PA_start))
alpha_optimal_final <- mean(alpha_optimal[ok_index])
estimated_fitness_final <- estimated_fitness[[ok_index[1]]]
#print("NA in estimated_fitness_final before sum:")
#print(sum(is.na(estimated_fitness_final)))
if (length(ok_index) > 1)
for (temp_jj in 2:length(ok_index))
estimated_fitness_final <- estimated_fitness_final + estimated_fitness[[ok_index[temp_jj]]]
estimated_fitness_final <- estimated_fitness_final / length(ok_index)
# print("NA in estimated_fitness_final after sum:")
#print(sum(is.na(estimated_fitness_final)))
### One pass to find optimal r #######
chosen_node_big <- names(cv_data$stats$z_j[cv_data$stats$z_j >= cv_deg_thresh[1]])
ratio_vec_PAFit <- sort(r,decreasing = TRUE)
PA_each <- rep(0,length(ratio_vec_PAFit))
names(PA_each) <- ratio_vec_PAFit
max_val <- -Inf
estimated_PA <- NULL
switch_flag <- 0
for (i in 1:length(ratio_vec_PAFit)) {
#print(paste0("Processing case ",count, " of a maximum of ",total))
if (!is.null(estimated_PA))
result_PAFit <- PAFit(cv_data$stats,
s = s_optimal_final,
r = ratio_vec_PAFit[i],
only_PA = TRUE,
true_f = estimated_fitness_final,
auto_stop = TRUE,
start_A = estimated_PA,
stop_cond = stop_cond,
normal_start_f = normal_start_f,
weight_power = weight_f ,
normalized_f = FALSE, ...)
else {
if (is.null(estimated_PA_start)) {
result_PAFit <- PAFit(cv_data$stats,
s = s_optimal_final,
r = ratio_vec_PAFit[i],
only_PA = TRUE,
true_f = estimated_fitness_final,
auto_stop = TRUE,
alpha_start = alpha_optimal_final,
stop_cond = stop_cond,
normal_start_f = normal_start_f,
normalized_f = FALSE, ...)
}
else
result_PAFit <- PAFit(cv_data$stats,
s = s_optimal_final,
r = ratio_vec_PAFit[i],
only_PA = TRUE,
true_f = estimated_fitness_final,
auto_stop = TRUE,
start_A = estimated_PA_start,
stop_cond = stop_cond,
normal_start_f = normal_start_f,
normalized_f = FALSE, ...)
}
#alpha_temp <- result_PAFit$alpha
if (TRUE == result_PAFit$diverge_zero)
PA_each[i] <- -Inf
else for (k in 1:length(cv_data$m_each))
if (cv_data$m_each[k] != 0) {
chosen_node <- chosen_node_big[!is.na(result_PAFit$A[cv_data$deg_each[k,chosen_node_big] + 1])]
chosen_node <- chosen_node[result_PAFit$A[cv_data$deg_each[k,chosen_node] + 1] != 0]
#PA <- cv_data$deg_each[k,chosen_node]^alpha_temp
PA <- result_PAFit$A[cv_data$deg_each[k,chosen_node] + 1]
if (sum(is.na(PA)) > 0) {
print("it should not happen here")
PA[is.na(PA)] <- cv_data$deg_each[k,chosen_node][is.na(PA)]^alpha_temp
}
#print(PA)
#PA[PA == 0] <- mean(result_PAFit$A)
if (sum(PA == 0) > 0) {
print("This should not happen")
PA[PA == 0] <- cv_data$deg_each[k,chosen_node][PA == 0]^alpha_temp
}
fitness <- rep(1,dim(cv_data$deg_each[,chosen_node, drop = FALSE])[2])
names(fitness) <- colnames(cv_data$deg_each[,chosen_node, drop = FALSE])
fitness[chosen_node] <- result_PAFit$f[chosen_node]
prob_PAFit <- PA * fitness
prob_PAFit <- prob_PAFit / sum(prob_PAFit,na.rm = TRUE)
#temp_sapply <- sapply(prob_PAFit,is.na)
#print(prob_PAFit)
#print(temp_sapply)
if (length(prob_PAFit) > 0) {
prob_PAFit[sapply(prob_PAFit,is.na)] <- 0
PA_each[i] <- PA_each[i] +
FitMultinomial(true = as.vector(prob_PAFit), dat = as.vector(cv_data$prob_em_each[k,chosen_node] *
cv_data$m_each[k]))
}
}
if (i == 1) {
lambda_optimal <- result_PAFit$lambda
r_optimal <- ratio_vec_PAFit[i]
max_val <- PA_each[i]
estimated_PA <- result_PAFit$theta
alpha_op_final <- result_PAFit$alpha
} else if (PA_each[i] > max_val) {
lambda_optimal <- result_PAFit$lambda
r_optimal <- ratio_vec_PAFit[i]
max_val <- PA_each[i]
estimated_PA <- result_PAFit$theta
alpha_op_final <- result_PAFit$alpha
}
}
#if (print_out == TRUE)
#print(paste0("s_optimal = ",s_optimal_final, " ,r_optimal = ", r_optimal))
result <- list(r_optimal = r_optimal,
lambda_optimal = lambda_optimal,
s_optimal = s_optimal_final,
estimated_fitness = estimated_fitness_final,
estimated_PA = estimated_PA,
cv_deg_thresh = cv_deg_thresh,
alpha_optimal = alpha_op_final,
center_k = result_PAFit$center_k)
class(result) <- "CV_Result"
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.