text:
#######################################
#-------------PENETRATION-------------#
#######################################

pen.start.time <- Sys.time()

#set response variable
y_var <- "Buyer_Pos_P1"

# drop variables with only one level
finaldata_pen <- as.data.table(initial_data)
indices <- sapply(finaldata_pen,function(x) length(unique(x))==1)
finaldata_pen <- finaldata_pen[,!indices,with=F]

#finaldata_pen <- finaldata_pen[sample(nrow(finaldata_pen),0.50*nrow(finaldata_pen)),]

# set mandatory variables
pen_mandatory_covariate <- c("Buyer_Pre_P1")

# final data for the modeling process
finaldata_pen <- as.data.frame(finaldata_pen)
finaldata_pen <- finaldata_pen[, !(names(finaldata_pen) %in% c("Trps_POS_P1","Dol_per_Trip_POS_P1","Trps_PRE_P1","Dol_per_Trip_PRE_P1","Nonbuyer_Pre_P1"))]
finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% occ_dolocc_outliers_panid_to_remove,]

# set initial model formula
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_pen[,!names(finaldata_pen) %in% c(y_var,"panid",random_demos,random_campaigns,"Buyer_Pre_P1",scoring_vars)]),collapse=" + ")," "),"+ offset(log(Buyer_Pre_P1+1))")
options(warn=-1)

# run glm model
model_glm_total <- glm(Formula_fixed,data=finaldata_pen,family=binomial(link="logit"))
summary(model_glm_total)

##### Variable selection - Condition 1 #####
# Remove pvalues greater than 0.7 and rerun the glm model
pval_excluded_vars <- setdiff(c(rownames(as.data.frame(which(summary(model_glm_total)$coefficient[,4]>0.7))),setdiff(names(model_glm_total$coefficients),rownames(summary(model_glm_total)$coefficients))),pen_mandatory_covariate)
finaldata_pen <- finaldata_pen[, !names(finaldata_pen) %in% c(pval_excluded_vars)]

# re-run glm model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_pen[,!names(finaldata_pen)%in%c(y_var,"panid",random_demos,random_campaigns,"Buyer_Pre_P1",scoring_vars)]),collapse=" + ")," "),"+ offset(log(Buyer_Pre_P1+1))")
model_glm_condition1 <- glm(Formula_fixed,data=finaldata_pen,family=binomial(link="logit"))
summary(model_glm_condition1)
rm(model_glm_total)

##### Variable selection - Condition 2 #####
# Remove variables that have z value lower than 2 and VIF greater than 15

# get variables with z value lower than 2
z_excluded_vars <- rownames(as.data.frame(which(abs(summary(model_glm_condition1)$coefficient[,3])<1.96)))

# get the variables with VIF greater than 15
multicol <- alias(model_glm_condition1)
multicols_toremove <- rownames(as.data.frame(multicol$Complete))

if (length(multicols_toremove)>0){
  finaldata_pen <- finaldata_pen[ , !names(finaldata_pen) %in% multicols_toremove]
  Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_pen[,!names(finaldata_pen)%in%c(y_var,"panid",random_demos,random_campaigns,"Buyer_Pre_P1")]),collapse=" + ")," "),"+ offset(log(Buyer_Pre_P1+1))")
  model_glm_new <- glm(Formula_fixed,data=finaldata_pen,family=binomial(link="logit"))
} else {
  model_glm_new <- model_glm_condition1
}

vif_output <- as.data.frame(vif(model_glm_new))
vif_to_remove <- rownames(vif_output)[(vif(model_glm_new)) > 15]

# final vars to exclude under condition 2
vars_condition2 <- setdiff(intersect(z_excluded_vars,vif_to_remove),pen_mandatory_covariate)

# re-run glm model
finaldata_pen <- finaldata_pen[ , !names(finaldata_pen) %in% c(vars_condition2,pval_excluded_vars)]
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_pen[,!names(finaldata_pen)%in%c(y_var,"panid",random_demos,random_campaigns,"Buyer_Pre_P1",scoring_vars)]),collapse=" + ")," ")," + offset(log(Buyer_Pre_P1+1))")
model_glm <- glm(Formula_fixed,data=finaldata_pen,family=binomial(link="logit"))
summary(model_glm)
rm(model_glm_condition1)

# p-values and signs checks
pvalues <- summary(model_glm)$coefficients[,4][-1]
keep_covariate <- names(pvalues[pvalues<pvalue_lvl])

# final cleaning of first model based on signs and pvalues
neutralvars <- setdiff(names(finaldata_pen),c(positivevars,negativevars))
coeffs <- sign(summary(model_glm)$coefficients[,1])
varstokeep <- c();initialvars <- c()
for (i in negativevars) if (i %in% names(coeffs[coeffs<0])) varstokeep <- c(varstokeep,i)
for (j in positivevars) if (j %in% names(coeffs[coeffs>0])) varstokeep <- c(varstokeep,j)
for (l in neutralvars) varstokeep <- c(varstokeep,l)
for (k in varstokeep) if (k %in% keep_covariate) initialvars <- c(initialvars,k)

#set the formula for the second model after cleaning the first model
Formula_final <- paste0(paste0(y_var,"~ ",paste0(initialvars,collapse=" + ")," "),"+ group + offset(log(Buyer_Pre_P1+1))")
model_glm_2 <- glm(Formula_final,data=finaldata_pen,family=binomial(link="logit"))
summary(model_glm_2)
rm(model_glm)

# cleaning of second model based on signs and pvalues
pvalues_2 <- summary(model_glm_2)$coefficients[,4][-1]
keep_covariate_2 <- names(pvalues_2[pvalues_2<pvalue_lvl])
coeffs_2 <- sign(summary(model_glm_2)$coefficients[,1])
varstokeep_2 <- c();vars_2 <- c()
for (i in negativevars) if (i %in% names(coeffs_2[coeffs_2<0])) varstokeep_2 <- c(varstokeep_2,i)
for (j in positivevars) if (j %in% names(coeffs_2[coeffs_2>0])) varstokeep_2 <- c(varstokeep_2,j)
for (l in neutralvars) varstokeep_2 <- c(varstokeep_2,l)
for (k in varstokeep_2) if (k %in% keep_covariate_2) vars_2 <- c(vars_2,k)

#######################################
#------------ Final Model ------------#
#######################################

#set the formula for the final model
Formula_final <- paste0(paste0(y_var,"~ ",paste0(vars_2,collapse=" + ")," "),"+ group + offset(log(Buyer_Pre_P1+1))")
pen_final_model <- glm(Formula_final,data=finaldata_pen,family=binomial(link="logit"))
summary(pen_final_model)
rm(model_glm_2)

# final model
# final cleaning of model based on signs
pvalues_2 <- summary(pen_final_model)$coefficients[,4][-1]
keep_covariate_2 <- names(pvalues_2[pvalues_2<pvalue_lvl])
coeffs <- sign(summary(pen_final_model)$coefficients[,1])
finalvars1 <- c();finalvars <- c()
for (i in negativevars) if (i %in% names(coeffs[coeffs<0])) finalvars1 <- c(finalvars1,i)
for (j in positivevars) if (j %in% names(coeffs[coeffs>0])) finalvars1 <- c(finalvars1,j)
for (l in neutralvars) if (l %in% keep_covariate_2) finalvars1 <- c(finalvars1,l)
for (k in finalvars1) if (k %in% keep_covariate_2) finalvars <- c(finalvars,k)

# final correlation check (find pair vars with correlation greater then 0.8 and drop the one with bigger pvalue)
finaldata_pen_num <- as.data.table(finaldata_pen[,names(finaldata_pen) %in% finalvars])
t1 <- as.data.frame(finaldata_pen_num[,sapply(finaldata_pen_num,is.numeric)])
finaldata_pen_num <- as.data.frame(finaldata_pen_num)
finaldata_pen_num <- finaldata_pen_num[,c(which(t1==T))]
temp_cor <- as.matrix(cor(finaldata_pen_num))
temp_cor[lower.tri(temp_cor, diag = T)] <- NA
temp_corIndex <- which(temp_cor > 0.8 |temp_cor < (-0.8), arr.ind = TRUE)
cordata <- data.table(col1=colnames(temp_cor)[temp_corIndex[,1]],col2=colnames(temp_cor)[temp_corIndex[,2]])

if(nrow(cordata)>0) {
  drop_cor_vars <- c()
  for (i in 1:nrow(cordata)) {
    if(cordata[i][,1,with=F]=="group"|cordata[i][,2,with=F]=="group"){
      if(cordata[i][,1,with=F]=="group") drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,2,with=F]))
      if(cordata[i][,2,with=F]=="group") drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,1,with=F]))
    } else {
      if (summary(pen_final_model)$coefficients[which(rownames(summary(pen_final_model)$coefficients)==as.character(cordata[i][,1,with=F])),4] > summary(pen_final_model)$coefficients[which(rownames(summary(pen_final_model)$coefficients)==as.character(cordata[i][,2,with=F])),4]) {
        drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,1,with=F]))
      } else {
        drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,2,with=F]))
      }
    }
  }
  finalvars <- setdiff(finalvars,drop_cor_vars)
}

#set the final formula
Formula_final <- paste0(paste0(y_var,"~ ",paste0(finalvars,collapse=" + ")," ")," + group  + offset(log(Buyer_Pre_P1+1))")
pen_final_model <- glm(Formula_final,data=finaldata_pen,family=binomial(link="logit"))
summary(pen_final_model)

options(warn=0)

# # kruskal wallis test for pair comparison of publisher's levels only for penetration
# # Pen model
# finaldata_pen <- as.data.table(finaldata_pen)
# finaldata_pen[,Prediction:=predict(pen_final_model, type="response")]
# finaldata_pen <- droplevels(finaldata_pen)
#
# kruskal_var <- 'Prediction'
# kruskal_randoms <- c()
# for(i in random_campaigns){
#   kt <- kruskal.test(eval(parse(text=kruskal_var))~eval(parse(text=i)),data=finaldata_pen)
#   print(i);print(kt)
#   if(kt$p.value<0.10){
#     kruskal_randoms <- append(kruskal_randoms,i)
#   }
# }
#
# warning('In case of aggregation based on Kruskal-Wallis test, make sure that the other level of each sub-campaign is named as Other_XXX')
# for(i in kruskal_randoms){
#   lev <- setdiff(levels(finaldata_pen[,eval(parse(text=i))]),'none')
#   if(length(lev)>2){
#     aggr <- c()
#     for(j in 1:(length(lev)-1)){
#       for( k in (j+1):(length(lev))){
#         df <- subset(finaldata_pen, eval(parse(text=i))  %in% c(lev[j],lev[k]))
#         l1 <- list(c(lev[j],lev[k]),P_Value=try(round(kruskal.test(df[,eval(parse(text=kruskal_var))]~df[,eval(parse(text=i))])$p.value,5),silent=TRUE))
#         if(l1[2]=="NaN") next
#         vectorElements <- unlist(l1)
#         listnames <- names(vectorElements) <- c(paste0(i,':'),paste0(i,':'),'P_value:')
#         listout <- paste(listnames,vectorElements)
#         write.table(listout,file=paste0(final_Occ_modelling_output,'Pen_Kruskal.csv'), row.names = FALSE, col.names = FALSE, quote = FALSE,append=TRUE)
#         if(l1[['P_Value']]<0.05){
#           aggr <- append(aggr,c(sapply(l1[1],`[`,1),sapply(l1[1],`[`,2)))
#           aggr <- setdiff(aggr,c('none',paste0("Other_",i)))
#         }
#       }
#     }
#     if(length(unique(aggr))>0){
#       finaldata_pen[!eval(parse(text=i)) %in% c(aggr,'none'),eval(parse(text=i)):=paste0("Other_",i)]
#     } else {
#       finaldata_pen[,eval(parse(text=i)):=NULL]
#     }
#   }
# }

finaldata_pen <- as.data.table(finaldata_pen)
finaldata_pen <- droplevels(finaldata_pen)
pen_random_campaigns <- intersect(random_campaigns,names(finaldata_pen))

if(length(pen_random_campaigns)>0){
  pen_random_campaigns_factors <- paste0("+(0+group|",pen_random_campaigns,")")
} else {
  pen_random_campaigns_factors <- NULL
}
random_demos_factors <- paste0( "+(0+group|",random_demos,")" )

##################################################################################################################

# run glmer model
pen.mixed.start.time <- Sys.time()
Formula_final <- paste0(c(paste0(paste0(y_var,"~ ",paste0(finalvars,collapse=" + ")," "),"+ offset(log(Buyer_Pre_P1+1)) + group"),pen_random_campaigns_factors),collapse='')
pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
pen.mixed.model.execution.time <- Sys.time()-pen.mixed.start.time
out_resid <- order(residuals(pen_final_model),decreasing = T)
finaldata_pen <- as.data.table(finaldata_pen)
pen_total_outliers_panid <- finaldata_pen[out_resid,panid]
segment1_panid <- data_NB_NE_B[,panid]
#segment2_panid <- data_B_E_NB[,panid]
outliers1_panid <- intersect(pen_total_outliers_panid,segment1_panid)
outliers2_panid <- c()
#outliers2_panid <- intersect(pen_total_outliers_panid,segment2_panid)
pen_outliers_panid_segm <- c(outliers1_panid,outliers2_panid)

# create a txt file with the summary of the first glmer model
logsavefname <- "Pen_First_Run_Results.txt"
sink(paste0(initial_Pen_modelling_output,logsavefname))
print(Sys.time()-pen.start.time)
print(summary(pen_final_model)$coef)
print(ranef(pen_final_model))
sink()

# remove outliers if needed (negative estimation for exposure)
pen_outliers_panid_to_remove <- c()
pen_outliers_panid_segm_new <- c()
pen_outliers_panid_to_remove_threshold <- c()

flag <- c()
coef_flag <- c()
for(i in setdiff(levels(finaldata_pen$group),"0")) flag <- append(flag,paste0('group',i))

for(j in flag){
  coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
}

if(length(pen_random_campaigns)>0){
  min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  max_random_coeff <- max(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=max)))
} else {
  min_random_coeff <- 0
}

if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
  out_resid_positive <- order(residuals(pen_final_model),decreasing = T)
  out_resid_negative <- order(residuals(pen_final_model),decreasing = F)
  pen_outliers_panid_positive_extreme_panid <- finaldata_pen[out_resid_positive,panid]
  pen_outliers_panid_negative_extreme_panid <- finaldata_pen[out_resid_negative,panid]
  for(i in 1:length(pen_random_campaigns)){
    min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=min)))
    if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
      max_random_coeff <- max(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=max)))
      random_frame <- as.data.frame(ranef(pen_final_model)[pen_random_campaigns][i])
      random_table <- as.data.table(random_frame)
      random_table[,level:=row.names(random_frame)]
      negative_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
      positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))== max_random_coeff]$level
      negative_extreme_panid <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
      positive_extreme_panid  <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
      pen_outliers_panid_segm1 <- intersect(pen_outliers_panid_negative_extreme_panid,negative_extreme_panid)
      pen_outliers_panid_segm1 <- pen_outliers_panid_segm1[1:round(0.02*(length(negative_extreme_panid)))]
      if(max_random_coeff > 0.05) {
        pen_outliers_panid_segm2 <- intersect(pen_outliers_panid_positive_extreme_panid,positive_extreme_panid)
        pen_outliers_panid_segm2 <- pen_outliers_panid_segm2[1:round(0.02*(length(positive_extreme_panid)))]
      } else {
        pen_outliers_panid_segm2 <- c()
      }
      pen_outliers_panid_segm_new <- c(pen_outliers_panid_segm1,pen_outliers_panid_segm2)
      finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_segm_new,]
    }
  }
  pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
  coef_flag <- c()
  for(j in flag){
    coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
  }
  min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
    for(i in 1:length(pen_random_campaigns)){
      min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=min)))
      if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
        max_random_coeff <- max(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=max)))
        random_frame <- as.data.frame(ranef(pen_final_model)[pen_random_campaigns][i])
        random_table <- as.data.table(random_frame)
        random_table[,level:=row.names(random_frame)]
        negative_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
        positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))== max_random_coeff]$level
        negative_extreme_panid <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
        positive_extreme_panid  <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
        pen_outliers_panid_segm1 <- intersect(pen_outliers_panid_negative_extreme_panid,negative_extreme_panid)
        pen_outliers_panid_segm1 <- pen_outliers_panid_segm1[1:round(0.05*(length(negative_extreme_panid)))]
        if(max_random_coeff > 0.05) {
          pen_outliers_panid_segm2 <- intersect(pen_outliers_panid_positive_extreme_panid,positive_extreme_panid)
          pen_outliers_panid_segm2 <- pen_outliers_panid_segm2[1:round(0.05*(length(positive_extreme_panid)))]
        } else {
          pen_outliers_panid_segm2 <- c()
        }
        pen_outliers_panid_segm_new <- c(pen_outliers_panid_segm1,pen_outliers_panid_segm2)
        finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_segm_new,]
      }
    }
    pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  }
  if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
    for(i in 1:length(pen_random_campaigns)){
      min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=min)))
      if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
        max_random_coeff <- max(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns][i],FUN=max)))
        random_frame <- as.data.frame(ranef(pen_final_model[pen_random_campaigns])[i])
        random_table <- as.data.table(random_frame)
        random_table[,level:=row.names(random_frame)]
        negative_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
        positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i]),'.group1')))== max_random_coeff]$level
        negative_extreme_panid <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
        positive_extreme_panid  <- finaldata_pen[eval(parse(text=paste0(names(ranef(pen_final_model)[pen_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
        pen_outliers_panid_segm1 <- intersect(pen_outliers_panid_negative_extreme_panid,negative_extreme_panid)
        pen_outliers_panid_segm1 <- pen_outliers_panid_segm1[1:round(0.10*(length(negative_extreme_panid)))]
        if(max_random_coeff > 0.05) {
          pen_outliers_panid_segm2 <- intersect(pen_outliers_panid_positive_extreme_panid,positive_extreme_panid)
          pen_outliers_panid_segm2 <- pen_outliers_panid_segm2[1:round(0.10*(length(positive_extreme_panid)))]
        } else {
          pen_outliers_panid_segm2 <- c()
        }
        pen_outliers_panid_segm_new <- c(pen_outliers_panid_segm1,pen_outliers_panid_segm2)
        finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_segm_new,]
      }
    }
    pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  }
}

pvalue_flag <- c()
for(j in flag){
  pvalue_flag <- append(pvalue_flag,coef(summary(pen_final_model))[j,"Pr(>|z|)"])
}

if(length(pen_random_campaigns)>0){
  min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  max_random_coeff <- max(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=max)))
} else {
  min_random_coeff <- 0
}

# drop outliers in case of high exposure estimation (in case we want a lower group1 coef to meet business requirements)
# segment_panid <- pen_reduction[,panid]
# outliers_panid <- intersect(pen_total_outliers_panid,segment_panid)
# pen_outliers_panid_to_remove_threshold <- outliers_panid[c(1:round(0.10*(length(outliers_panid))))]
# finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_to_remove_threshold,]
# pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))

if(min(coef_flag) + min_random_coeff < 0 | max(pvalue_flag)>0.40){
  cat('\n')
  cat("removing outliers ...\nEstimated Iteration Execution ");print(pen.mixed.model.execution.time)
  cat("\nIteration 1 ... Start Time: ");print(Sys.time())
  pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.05*(nrow(data_NB_NE_B))))]
  #pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.05*(nrow(data_NB_NE_B)+nrow(data_B_E_NB))))]
  finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_to_remove,]
  pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
  coef_flag <- c()
  pvalue_flag <- c()
  for(j in flag){
    coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
    pvalue_flag <- append(pvalue_flag,coef(summary(pen_final_model))[j,"Pr(>|z|)"])
  }
  min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
  if(min(coef_flag) + min_random_coeff < 0 | max(pvalue_flag)>0.40){
    pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.10*(nrow(data_NB_NE_B))))]
    #pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.10*(nrow(data_NB_NE_B)+nrow(data_B_E_NB))))]
    finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_to_remove,]
    cat("\nIteration 2 ... Start Time: ");print(Sys.time())
    pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    pvalue_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
      pvalue_flag <- append(pvalue_flag,coef(summary(pen_final_model))[j,"Pr(>|z|)"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
    if(min(coef_flag) + min_random_coeff < 0 | max(pvalue_flag)>0.40){
      pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.15*(nrow(data_NB_NE_B))))]
      #pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.15*(nrow(data_NB_NE_B)+nrow(data_B_E_NB))))]
      finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_to_remove,]
      cat("\nIteration 3 ... Start Time: ");print(Sys.time())
      pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
      coef_flag <- c()
      pvalue_flag <- c()
      for(j in flag){
        coef_flag <- append(coef_flag,coef(summary(pen_final_model))[j,"Estimate"])
        pvalue_flag <- append(pvalue_flag,coef(summary(pen_final_model))[j,"Pr(>|z|)"])
      }
      min_random_coeff <- min(unlist(lapply(ranef(pen_final_model)[pen_random_campaigns],FUN=min)))
      if(min(coef_flag) + min_random_coeff < 0 | max(pvalue_flag)>0.40){
        pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.20*(nrow(data_NB_NE_B))))]
        #pen_outliers_panid_to_remove <- pen_outliers_panid_segm[c(1:round(0.20*(nrow(data_NB_NE_B)+nrow(data_B_E_NB))))]
        finaldata_pen <- finaldata_pen[!finaldata_pen$panid %in% pen_outliers_panid_to_remove,]
        cat("\nFinal Iteration ... Start Time: ");print(Sys.time())
        pen_final_model <- glmer(Formula_final,data=finaldata_pen,family=binomial(link="logit"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
      }
    }
  }
}

write.csv(as.data.frame(summary(pen_final_model)$coefficients),file=paste0(final_Pen_modelling_output,"Pen_fixed_effects.csv"))
write.csv(as.data.frame(summary(pen_final_model)$coefficients),file=paste0(scoring_Pen_output,"Pen_fixed_effects.csv"))

# create output for random effects
randoms <- as.data.frame(augment.ranef.mer(ranef(pen_final_model,condVar=TRUE)))
random.ef_temp.report <- randoms[,!names(randoms) %in% c("qq")]
write.csv(random.ef_temp.report,file=paste0(final_Pen_modelling_output,"Pen_random_effects_pvalue.csv"),row.names=FALSE)

for(i in pen_random_campaigns){
  random.ef_temp <- as.data.frame(ranef(pen_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp[2],keep.rownames=TRUE)
  random.coef <- random.coef_temp[,label:=names(random.coef_temp)[2]]
  write.table(random.coef,file=paste0(final_Pen_modelling_output,'Pen_random_effects_coeffs.csv'), quote = FALSE,append    = i > 1,
              sep= ",", row.names = FALSE,col.names=FALSE)
}

for(i in pen_random_campaigns){
  random.ef_temp <- as.data.frame(se.ranef(pen_final_model)[i])
  colnames(random.ef_temp)[2] <- paste0(i,' S.E.')
  options(warn=-1);write.table(random.ef_temp[2],file=paste0(final_Pen_modelling_output,'Pen_SE_random_effects.csv'), quote = FALSE,append    = i > 1,
                               sep= ",", row.names = FALSE);options(warn=0)
}

pen_random_demos <- c()
random_groupe <- c(pen_random_demos,pen_random_campaigns)

if(length(pen_random_demos)>0){
  random_group_demos <- data.frame(random_demos)
  names.random_group_demos <- c("class")
  setnames(random_group_demos,names.random_group_demos)
  random_group_demos$exposed <- 0
} else {
  random_group_demos <- c()
}

random_group_campaigns <- data.frame(pen_random_campaigns)
names.random_group_campaigns <- c("class")
setnames(random_group_campaigns,names.random_group_campaigns)
random_group_campaigns$exposed <- 1

random_group_exposed <- rbind( random_group_campaigns,random_group_demos)

random_group_data <- data.frame()
for(i in random_groupe){
  random.ef_temp <- as.data.frame(ranef(pen_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp,keep.rownames=TRUE)
  random.coef <- random.coef_temp[,class:=i]
  names.random.coef <- c("level","group0","group1","class")
  setnames(random.coef ,names.random.coef)
  random.coef$row <- row.names(random.coef)
  random_group_data <-  rbind( random.coef,random_group_data)
}

random_group_SEdata <- data.frame()
for(i in random_groupe){
  random.sef_temp <- as.data.frame(se.ranef(pen_final_model)[i])
  random.secoef_temp <- as.data.table(random.sef_temp,keep.rownames=TRUE)
  random.secoef <- random.secoef_temp[,class:=i]
  names.random.secoef <- c("levelSE","SE_group0","SE_group1","class")
  setnames(random.secoef ,names.random.secoef)
  random.secoef$row <- row.names(random.secoef)
  random_group_SEdata <-  rbind( random.secoef,random_group_SEdata)
}

final_random <- merge(random_group_data,random_group_SEdata, by=c("class","row"))
final_random <- merge(final_random,random_group_exposed, by=c("class"))
View(final_random)
write.csv(final_random,file=paste0(scoring_Pen_output,'Pen_random_effects.csv'), row.names = FALSE)

# create a txt file with the summary of the glmer model
logsavefname <- "Pen_Final_Results.txt"
sink(paste0(final_Pen_modelling_output,logsavefname))
print(Sys.time()-pen.start.time)
print(summary(pen_final_model))
print(ranef(pen_final_model))
sink()

for(i in pen_random_campaigns){
  pdf(paste0(final_Pen_modelling_output,"Pen_Caterpillar_Plots ",i,".pdf"),width=7,height=5)
  print(ggCaterpillar(ranef(pen_final_model, condVar=TRUE)[i]))
  dev.off()
}

# create descriptive plots for the Pen breaks based on  the predicted values
finaldata_pen <- as.data.table(finaldata_pen)
finaldata_pen[,Prediction_Pen:=predict(pen_final_model, type="response")]
finaldata_pen <- droplevels(finaldata_pen)

response_var <- c("Prediction_Pen")
for(i in random_campaigns){
  data <- as.data.table(copy(finaldata_pen))
  data <- data[eval(parse(text=i))!="none",]
  data <- droplevels(data)
  data <- as.data.frame(data)
  rr <- c(response_var,i)
  df <- data[,rr]
  res <- ldply(levels(df[,i]),function(x) round(mean(df[df[,i]==x,]$Prediction_Pen),2))
  df <-cbind(res,levels(df[,i]))
  colnames(df) <- c('Buyer_Pos_Prob','Levels')
  p <- ggplot(df,aes(x=Levels,y=Buyer_Pos_Prob))+geom_bar(stat='identity')+ theme(axis.text.x = element_text(angle=90, vjust=0.5))+ylab('Buyer_Pos_Prob')+
    theme(legend.position="none")+theme(panel.background = element_rect(fill = 'white'))+ #geom_hline(yintercept=mean(finaldata_pen$Prediction_Pen),color="red",size=0.5)+
    geom_text(stat='identity',aes(label=Buyer_Pos_Prob,group=Levels), position = position_dodge(0.9), vjust = 0)
  ggsave(filename=paste0(output_checks,i,'_Buyer_Pos_Prob.pdf'), width=4, height=5,p)
}

# save Penetration final dataset and final model
Pen_Formula_final <- Formula_final
all_pen_outliers_panid_to_remove <- c(pen_outliers_panid_to_remove,pen_outliers_panid_segm_new,pen_outliers_panid_to_remove_threshold)
rm(finaldata_pen_num,temp_cor,cordata)

# save data in .Rdata file
save(finaldata_pen,Pen_Formula_final,pen_final_model,all_pen_outliers_panid_to_remove,file=paste0(output,"saved_data/","pen_final_model_data.Rdata"))

# --------------------------------------------------------------------------------------------------------------------------------------------------
# ---- END OF PENNETRATION --------
# --------------------------------------------------------------------------------------------------------------------------------------------------

# Check and Rerun Occasion and DolOcc models due to outliers removal if needed
if(length(c(pen_outliers_panid_segm_new,pen_outliers_panid_to_remove))>0){
  #load("Occasion_initial_model_data.Rdata")
  cat('\n')
  cat("Re-Running Occasion Model ...")
  finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% c(pen_outliers_panid_segm_new,pen_outliers_panid_to_remove),]

  # run glmer model
  occ_final_model <- glmer(Occ_Random_Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
}

# save data in .Rdata file
save(finaldata_occ,occ_final_model,Occ_Random_Formula_final,file=paste0(output,"saved_data/","Occasion_final_model_data.Rdata"))

write.csv(as.data.frame(summary(occ_final_model)$coefficients),file=paste0(final_Occ_modelling_output,"Occ_fixed_effects.csv"))
write.csv(as.data.frame(summary(occ_final_model)$coefficients),file=paste0(scoring_Occ_output,"Occ_fixed_effects.csv"))

# create output for random effects
randoms <- as.data.frame(augment.ranef.mer(ranef(occ_final_model,condVar=TRUE)))
random.ef_temp.report <- randoms[,!names(randoms) %in% c("qq")]
write.csv(random.ef_temp.report,file=paste0(final_Occ_modelling_output,"Occ_random_effects_pvalue.csv"),row.names=FALSE)

for(i in random_demos){
  random.ef_temp <- as.data.frame(ranef(occ_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp,keep.rownames=TRUE)
  random.coef <- random.coef_temp[,class:=i]
  names(random.coef) <- c("","group0","group1","class")
  options(warn=-1);write.table(random.coef,file=paste0(final_Occ_modelling_output,'Occ_random_demographics_effects.csv'), quote = FALSE,append    = i > 1,
                               sep= ",", row.names = FALSE,col.names=ifelse(i %in% random_demos[1], TRUE, FALSE));options(warn=0)
}

for(i in occ_random_campaigns){
  random.ef_temp <- as.data.frame(ranef(occ_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp[2],keep.rownames=TRUE)
  random.coef <- random.coef_temp[,label:=names(random.coef_temp)[2]]
  write.table(random.coef,file=paste0(final_Occ_modelling_output,'Occ_random_effects_coeffs.csv'), quote = FALSE,append    = i > 1,
              sep= ",", row.names = FALSE,col.names=FALSE)
}

for(i in random_demos){
  random.ef_temp <- as.data.frame(se.ranef(occ_final_model)[i])
  options(warn=-1);write.table(random.ef_temp,file=paste0(final_Occ_modelling_output,'Occ_SE_random_demographics_effects.csv'), quote = FALSE,append    = i > 1,
                               sep= ",", row.names = FALSE);options(warn=0)
}

for(i in occ_random_campaigns){
  random.ef_temp <- as.data.frame(se.ranef(occ_final_model)[i])
  colnames(random.ef_temp)[2] <- paste0(i,' S.E.')
  options(warn=-1);write.table(random.ef_temp[2],file=paste0(final_Occ_modelling_output,'Occ_SE_random_campaigns_effects.csv'), quote = FALSE,append    = i > 1,
                               sep= ",", row.names = FALSE);options(warn=0)
}

random_groupe <- c(random_demos,occ_random_campaigns)

if(length(random_demos)>0){
  random_group_demos <- data.frame(random_demos)
  names.random_group_demos <- c("class")
  setnames(random_group_demos,names.random_group_demos)
  random_group_demos$exposed <- 0
} else {
  random_group_demos <- c()
}

random_group_campaigns <- data.frame(occ_random_campaigns)
names.random_group_campaigns <- c("class")
setnames(random_group_campaigns,names.random_group_campaigns)
random_group_campaigns$exposed <- 1

random_group_exposed <- rbind( random_group_campaigns,random_group_demos)

random_group_data <- data.frame()
for(i in random_groupe){
  random.ef_temp <- as.data.frame(ranef(occ_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp,keep.rownames=TRUE)
  random.coef <- random.coef_temp[,class:=i]
  names.random.coef <- c("level","group0","group1","class")
  setnames(random.coef ,names.random.coef)
  random.coef$row <- row.names(random.coef)
  random_group_data <-  rbind( random.coef,random_group_data)
}

random_group_SEdata <- data.frame()
for(i in random_groupe){
  random.sef_temp <- as.data.frame(se.ranef(occ_final_model)[i])
  random.secoef_temp <- as.data.table(random.sef_temp,keep.rownames=TRUE)
  random.secoef <- random.secoef_temp[,class:=i]
  names.random.secoef <- c("levelSE","SE_group0","SE_group1","class")
  setnames(random.secoef ,names.random.secoef)
  random.secoef$row <- row.names(random.secoef)
  random_group_SEdata <-  rbind( random.secoef,random_group_SEdata)
}

final_random <- merge(random_group_data,random_group_SEdata, by=c("class","row"))
final_random <- merge(final_random,random_group_exposed, by=c("class"))
View(final_random)
write.csv(final_random,file=paste0(scoring_Occ_output,'Occ_random_effects.csv'), row.names = FALSE)

# create a txt file with the summary of the glmer model
logsavefname <- "Occasion_Final_Results.txt"
sink(paste0(final_Occ_modelling_output,logsavefname))
print(Sys.time()-occ.start.time)
print(summary(occ_final_model))
print(ranef(occ_final_model))
sink()

for(i in c(occ_random_campaigns,random_demos)){
  pdf(paste0(final_Occ_modelling_output,"Occ_Caterpillar_Plots ",i,".pdf"),width=7,height=5)
  print(ggCaterpillar(ranef(occ_final_model, condVar=TRUE)[i]))
  dev.off()
}

if(length(unique(c(occ_outliers_panid_to_remove,occ_outliers_panid_to_remove_threshold,occ_outliers_panid_segm_new,pen_outliers_panid_segm_new,pen_outliers_panid_to_remove)))>0){
  #load("DolOcc_initial_model_data.Rdata")
  cat('\n')
  cat("Re-Running DolOcc Model ...")
  finaldata_dolocc <- finaldata_dolocc[!finaldata_dolocc$panid %in% unique(c(occ_outliers_panid_to_remove,occ_outliers_panid_to_remove_threshold,occ_outliers_panid_segm_new,pen_outliers_panid_segm_new,pen_outliers_panid_to_remove)),]

  # run glmer model
  dolocc_final_model <- glmer(dolocc_Formula_final,data=finaldata_dolocc,family=Gamma(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
}

# save data in .Rdata file
save(finaldata_dolocc,dolocc_final_model,dolocc_Formula_final,final_DollOcc_modelling_output,file=paste0(output,"saved_data/","DolOcc_final_model_data.Rdata"))

write.csv(as.data.frame(summary(dolocc_final_model)$coefficients),file=paste0(final_DollOcc_modelling_output,"DolOcc_fixed_effects.csv"))
write.csv(as.data.frame(summary(dolocc_final_model)$coefficients),file=paste0(scoring_DolOcc_output,"DolOcc_fixed_effects.csv"))

# create output for random effects
randoms <- as.data.frame(augment.ranef.mer(ranef(dolocc_final_model,condVar=TRUE)))
random.ef_temp.report <- randoms[,!names(randoms) %in% c("qq")]
write.csv(random.ef_temp.report,file=paste0(final_DollOcc_modelling_output,"DolOcc_random_effects_pvalue.csv"),row.names=FALSE)

for(i in dolocc_random_campaigns){
  random.ef_temp <- as.data.frame(ranef(dolocc_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp[2],keep.rownames=TRUE)
  random.coef <- random.coef_temp[,label:=names(random.coef_temp)[2]]
  write.table(random.coef,file=paste0(final_DollOcc_modelling_output,'DolOcc_random_effects_coeffs.csv'), quote = FALSE,append    = i > 1,
              sep= ",", row.names = FALSE,col.names=FALSE)
}

for(i in dolocc_random_campaigns){
  random.ef_temp <- as.data.frame(se.ranef(dolocc_final_model)[i])
  colnames(random.ef_temp)[2] <- paste0(i,' S.E.')
  options(warn=-1);write.table(random.ef_temp[2],file=paste0(final_DollOcc_modelling_output,'DolOcc_SE_random_effects.csv'), quote = FALSE,append    = i > 1,
                               sep= ",", row.names = FALSE);options(warn=0)
}

dolocc_random_demos <- c()
random_groupe <- c(dolocc_random_demos,dolocc_random_campaigns)

if(length(dolocc_random_demos)>0){
  random_group_demos <- data.frame(random_demos)
  names.random_group_demos <- c("class")
  setnames(random_group_demos,names.random_group_demos)
  random_group_demos$exposed <- 0
} else {
  random_group_demos <- c()
}

random_group_campaigns <- data.frame(dolocc_random_campaigns)
names.random_group_campaigns <- c("class")
setnames(random_group_campaigns,names.random_group_campaigns)
random_group_campaigns$exposed <- 1

random_group_exposed <- rbind( random_group_campaigns,random_group_demos)

random_group_data <- data.frame()
for(i in random_groupe){
  random.ef_temp <- as.data.frame(ranef(dolocc_final_model)[i])
  random.coef_temp <- as.data.table(random.ef_temp,keep.rownames=TRUE)
  random.coef <- random.coef_temp[,class:=i]
  names.random.coef <- c("level","group0","group1","class")
  setnames(random.coef ,names.random.coef)
  random.coef$row <- row.names(random.coef)
  random_group_data <-  rbind( random.coef,random_group_data)
}

random_group_SEdata <- data.frame()
for(i in random_groupe){
  random.sef_temp <- as.data.frame(se.ranef(dolocc_final_model)[i])
  random.secoef_temp <- as.data.table(random.sef_temp,keep.rownames=TRUE)
  random.secoef <- random.secoef_temp[,class:=i]
  names.random.secoef <- c("levelSE","SE_group0","SE_group1","class")
  setnames(random.secoef ,names.random.secoef)
  random.secoef$row <- row.names(random.secoef)
  random_group_SEdata <-  rbind( random.secoef,random_group_SEdata)
}

final_random <- merge(random_group_data,random_group_SEdata, by=c("class","row"))
final_random <- merge(final_random,random_group_exposed, by=c("class"))
View(final_random)
write.csv(final_random,file=paste0(scoring_DolOcc_output,'DolOcc_random_effects.csv'), row.names = FALSE)

# create a txt file with the summary of the glmer model
logsavefname <- "DolOcc_Final_Results.txt"
sink(paste0(final_DollOcc_modelling_output,logsavefname))
print(Sys.time()-dolocc.start.time)
print(summary(dolocc_final_model))
print(ranef(dolocc_final_model))
sink()

for(i in dolocc_random_campaigns){
  pdf(paste0(final_DollOcc_modelling_output,"DolOcc_Caterpillar_Plots ",i,".pdf"),width=7,height=5)
  print(ggCaterpillar(ranef(dolocc_final_model, condVar=TRUE)[i]))
  dev.off()
}

total_outliers_panid_to_remove <- unique(c(panid_out_all,occ_outliers_panid_to_remove,dolocc_outliers_panid_to_remove,all_pen_outliers_panid_to_remove))

# save all modeling data in .Rdata file
save(finaldata_dolocc,finaldata_pen,finaldata_occ,total_outliers_panid_to_remove,file=paste0(output,"saved_data/","Total_Modeling_data.Rdata"))

initial_data <- initial_data[!initial_data$panid %in% total_outliers_panid_to_remove,]
write.csv(initial_data,file=paste0(scoring_output,"final_data.csv"),row.names=F)

# create a summary for the outliers
logsavefname <- "Outliers Summary.txt"
sink(paste0(output,logsavefname))
print(paste0("Modeling Execution in: ",Sys.time()-Start.Time))
print(paste0("Number of outliers HH: ",length(unique(na.omit(total_outliers_panid_to_remove)))," from ",nrow(initial_data)))
print(cat("List of removed Experian IDs: "))
print(cat(na.omit(total_outliers_panid_to_remove)))
sink()

# calculate an estimation of total lift
total_lift_estimation <- ((1+coef(summary(pen_final_model))["group1","Estimate"])*(1+coef(summary(dolocc_final_model))["group1","Estimate"])*(1+coef(summary(occ_final_model))["group1","Estimate"])-1)
cat(paste0("last lift estimation: ",100*total_lift_estimation,"%"))

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
"""
## fromSkyline

Helper function to convert a Skyline vector to a full matrix.

### Type
```julia
fromSkyline(skyline::Vector{Float64}, kdiag::Vector{Int})
```

### Arguments
```julia
* skyline::Vector{Float64}     : 1D Line(nod, nodof)
* kdiag::Vector{Int}         : 2D Triangle(nod, nodof)
```
"""
function fromSkyline(skyline::Vector{Float64}, kdiag::Vector{Int})
	neq = size(kdiag, 1)
	km = zeros(neq, neq)
	km[1, 1] = skyline[kdiag[1]]
	for i in 2:neq
		km[i, i] = skyline[kdiag[i]]
		for j in (kdiag[i-1] + 1):kdiag[i]
			km[i, i - (kdiag[i]-j)] = skyline[j]
			km[i - (kdiag[i]-j), i] = skyline[j]
    end
  end
	km
end

ext: jl


text:
####
#### this is for extra drawings that are not meant to be overlaid with anything
#### else such as boxplot, heatmap, polarplot, pieplot
####

"""
    boxplot(...)

Erase previous drawings and add a boxplot. Missing values are allowed but not Infinities or Nans.
"""
function boxplot(ys...; axes=nothing, o...)
    axes = check_axes(axes)
    reset!(axes; parent=axes.parent) # always on fresh axes

    isempty(first(ys)) && throw(ArgumentError("Cannot display empty vectors."))

    # setting an empty struct first so that we can exploit the options
    # the actual data will be provided after analysis
    nobj = sum(size(y, 2) for y ∈ ys)
    bp = Boxplot(Matrix{Float64}(undef,0,0), nobj)
    set_properties!(bp; defer_preview=true, o...)

    # analyse data
    stats = Matrix{Float64}(undef, nobj, 6) # 1:wlow, 2:q25, 3:q50, 4:q75, 5:whigh, 6:mean
    outliers = Vector{Vector{Float64}}(undef, nobj)

    boxcounter = 1
    axmin = Inf
    axmax = -Inf

    for y ∈ ys
        for k ∈ Base.axes(y, 2)
            yk = collect(skipmissing(view(y, :, k)))
            if any(isnan, yk) || any(isinf, yk)
                throw(ArgumentError("Inf or NaN values not allowed in boxplot."))
            end

            q00, q25, q50, q75, q100 = quantile(yk, [.0, .25, .5, .75, 1.0])
            iqr = q75 - q25
            mu  = mean(yk)

            wrlength = bp.boxstyles[k].wrlength
            wlow  = q25 - wrlength * iqr
            whigh = q75 + wrlength * iqr
            if isinf(wrlength)
                wlow, whigh = q00, q100 # min/max values
            end

            stats[k, :] = [wlow, q25, q50, q75, whigh, mu]

            # outliers
            outliers[k] = filter(e->(e<wlow || whigh<e), yk)

            # Keep track of extremes to set axis limits afterwards
            axmin_ = min(q00, wlow)
            axmin_ -= 0.5abs(axmin_)
            axmin_ < axmin && (axmin = axmin_)
            axmax_ = max(q100, whigh)
            axmax_ += 0.5abs(axmax_)
            axmax_ > axmax && (axmax = axmax_)

            boxcounter += 1
        end
    end
    #
    bp.stats = stats
    push!(axes.drawings, bp)

    if bp.horiz # horizontal boxplot
        ylim(0, nobj+1)
        yticks(1:nobj)
        xlim(axmin, axmax)
        for k ∈ 1:nobj
            bp.boxstyles[k].oshow || continue
            nok = length(outliers[k])
            s = bp.boxstyles[k].omstyle # style of outliers
            if nok > 0
                scatter!(outliers[k], fill(k, nok); marker=s.marker, msize=s.msize, mcol=s.color)
            end
        end
    else # vertical boxplot
        xlim(0, nobj+1)
        xticks(1:nobj)
        ylim(axmin, axmax)
        for k ∈ 1:nobj
            bp.boxstyles[k].oshow || continue
            nok = length(outliers[k])
            s = bp.boxstyles[k].omstyle # style of outliers
            if nok > 0
                scatter!(fill(k, nok), outliers[k]; marker=s.marker, msize=s.msize, mcol=s.color)
            end
        end
    end
    return DrawingHandle(bp)
end

####
#### HEATMAP
####

"""
    heatmap_ticks(ax)

Returns ticks position centered appropriately for a heatmap.

### Example

```julia
heatmap(randn(2, 2))
xticks(heatmap_ticks("x"), ["var1", "var2"]; angle=45)
```
"""
function heatmap_ticks(ax::String; axes=gca())
    isempty(axes.drawings) && throw(ArgumentError("No heatmap found."))
    # it must necessarily be the first object since heatmap resets the axes
    axes.drawings[1] isa Heatmap || throw(ArgumentError("No heatmap found."))
    nrows, ncols = size(axes.drawings[1].data)
    ax_lc = lowercase(ax)
    ax_lc ∈ ["x", "x2"] && return (collect(0:(ncols-1)) .+ 0.5) ./ ncols
    ax_lc ∈ ["y", "y2"] && return (collect((nrows-1):-1:0) .+ 0.5) ./ nrows
    throw(ArgumentError("Unrecognised ax descriptor expected "*
                                               "one of [x, x2, y, y2]"))
end

"""
    heatmap(X; opts...)

Creates a heatmap corresponding to matrix `X`. Doing this with matrices larger than 100x100
with the GLE backend can be slow. As an indication, 500x500 takes about 15-20 seconds on a standard
laptop. Time scales with the number of elements (so 100x100 takes under a second).
You can provide ticks using `xticks(heatmap_ticks("x"), ["tick1", ...]; opts...)`.
See also [`heatmap_ticks`](@ref).
Note: handles missing values but not NaN or Inf.
"""
function heatmap(data::Matrix{<:CanMiss{Real}}; axes=nothing, o...)
    # TODO:
    #  - deal with tables instead of forcing matrices?
    # ================================================

    nrows, ncols = size(data)
    min(nrows, ncols) <= 1000 || throw(ArgumentError("The matrix is too large to be displayed " *
                                                     "as a heatmap."))

    axes = check_axes(axes)
    reset!(axes; parent=axes.parent) # always on fresh axes

    zmin, zmax = extrema(skipmissing(data))

    h = Heatmap(data=Matrix{Int}(undef, 0, 0), zmin=zmin, zmax=zmax)
    h.transpose = (nrows < ncols)
    set_properties!(h; defer_preview=true, o...)
    push!(axes.drawings, h)

    ncolors = length(h.cmap)

    minv, maxv = extrema(skipmissing(data))
    incr = (maxv - minv) / ncolors

    data2 = Matrix{Int}(undef, nrows, ncols)
    if iszero(incr)
        data2 .= ceil(ncols/2)
    else
        @inbounds for i=1:nrows, j=1:ncols
            dij = data[i,j]
            if ismissing(dij)
                data2[i,j] = ncolors + 1 # offset doesn't matter, anything ≥ 1 is fine
            else
                data2[i,j] = ceil(Integer, (dij - minv) / incr)
            end
        end
    end
    # NOTE: there may be a few 0 (min values), will be treated as 1 see add_sub_heatmap!

    # store the color assignments
    h.data = data2

    # set axis
    xlim(0,1); ylim(0,1)

    # add labels if provided to the heatmap object
    xplaces = heatmap_ticks("x")
    yplaces = heatmap_ticks("y")

    # used for sparse labelling in case ncols > 20
    stepx    = (ncols<=50)*10 + (50<ncols<=200)*50 + (200<ncols)*100
    lastx    = ncols - mod(ncols, stepx) + 1
    maskx    = collect(1:stepx:lastx) .- 1
    maskx[1] = 1
    stepy    = (nrows<=50)*10 + (50<nrows<=200)*50 + (200<nrows)*100
    lasty    = nrows - mod(nrows, stepy) + 1
    masky    = collect(1:stepy:lasty) .- 1
    masky[1] = 1

    if ncols <= 20
        xticks(xplaces, ["$i" for i ∈ 1:ncols]; axes=axes, fontsize=12)
    else
        xticks(xplaces[maskx], ["$i" for i ∈ maskx]; axes=axes, fontsize=12)
    end
    if nrows <= 20
        yticks(yplaces, ["$i" for i ∈ 1:nrows]; axes=axes, fontsize=12)
    else
        yticks(yplaces[masky], ["$i" for i ∈ masky]; axes=axes, fontsize=12)
    end
    x2ticks("off")
    y2ticks("off")
    return DrawingHandle(h)
end


text:
#changed generate_default_data_dictionary
"""
TODO: Fill me in with some stuff ...
"""
function optimize_specific_growth_rate(data_dictionary::Dict{String,Any})

    # get some stuff from dd -
    number_of_reactions = data_dictionary["number_of_reactions"]

    # update the c-vector -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]
    objective_coefficient_array[13] = -1

    # return -
    return data_dictionary
end

"""
TODO: Fill me in with some stuff ...
"""
function optimize_flux_at_index(index_array::Array{Int64,1}, data_dictionary::Dict{String,Any})

    # update the c-vector -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]

    # iterate through the index array, and flip the flag on the coefficient array -
    for index_value in index_array
        objective_coefficient_array[index_value] = -1
    end

    # return -
    return data_dictionary
end

function constrain_measured_metabolites(data_dictionary::Dict{String,Any}, path_to_measurements_file::String)

    # from the data dictionary, get the list of metabolite symbols -
    list_of_metabolite_symbols = data_dictionary["list_of_metabolite_symbols"]

    # load the measurements file -
    measurements_dictionary = JSON.parsefile(path_to_measurements_file)

    # get the stoichiometric_matrix -
    stoichiometric_matrix = copy(data_dictionary["stoichiometric_matrix"])

    # get the metabolite measurement information -
    absolute_metabolite_measurement_dictionary = measurements_dictionary["absolute_metabolite_measurements"]
    for (metabolite_key, local_measurement_dict) in absolute_metabolite_measurement_dictionary

        # ok, so we have a metabolite_key - find the index of this metabolite in the metabolite list -
        idx_metabolite_match = findall(list_of_metabolite_symbols .== metabolite_key)

        # if we have this metabolite, update the STM -
        if (isempty(idx_metabolite_match) == false)

            # get experimental values -
            mean_ratio_value = parse(Float64,local_measurement_dict["mean_value"])
            lower_bound = parse(Float64,local_measurement_dict["lower_bound"])
            upper_bound = parse(Float64,local_measurement_dict["upper_bound"])

            # compute a value -
            r_value = rand()
            value = (1 - r_value)*lower_bound + r_value*(upper_bound)

            # update the stm -
            metabolite_index = (getindex(idx_metabolite_match))[1]
            old_value = stoichiometric_matrix[metabolite_index,growth_rate_reaction_index]
            stoichiometric_matrix[metabolite_index,growth_rate_reaction_index] = (old_value - value)
        end
    end

    data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
    return data_dictionary
end

function constrain_measured_fluxes(data_dictionary::Dict{String,Any}, path_to_measurements_file::String)

    # TODO: is the path_to_measurements_file legit?
    measurements_dictionary = JSON.parsefile(path_to_measurements_file)

    # get some stuff from the dd -
    rxn_name_list = data_dictionary["list_of_reaction_name_strings"]
    flux_bounds_array = data_dictionary["flux_bounds_array"]
    (number_of_reactions,nc) = size(flux_bounds_array)

    # do we have any flux ratios?
    flux_ratio_dictionary_array = measurements_dictionary["flux_ratio_measurements"]

    # initialize -
    number_of_additional_constraints = length(flux_ratio_dictionary_array)
    constraint_counter = 1
    additional_constraint_array = zeros(number_of_additional_constraints,number_of_reactions)
    for (reaction_key, local_measurement_dict) in flux_ratio_dictionary_array

        # ok, so we have a reaction key - find the index of this reaction in the reaction list -
        idx_reaction_match = findall(rxn_name_list .== reaction_key)

        # if we have this reaction, then update the bounds array -
        if (isempty(idx_reaction_match) == false)

            # all fluxes are relative to glucose -
            glucose_uptake_index = 50   # GLCpts -

            # get ratio value -
            mean_ratio_value = parse(Float64,local_measurement_dict["mean_value"])

            # get the measurement CV -
            measurement_cv = parse(Float64,local_measurement_dict["coefficient_of_variation"])

            # build a distribution, and draw a sample -
            distrubution = Distributions.Normal(mean_ratio_value, measurement_cv*mean_ratio_value)
            sampled_ratio_value = rand(distrubution)

            # get the actual index -
            idx_reaction = (getindex(idx_reaction_match))[1]

            # update the constraint -
            additional_constraint_array[constraint_counter,idx_reaction] = 1.0
            additional_constraint_array[constraint_counter,glucose_uptake_index] = -1*sampled_ratio_value

            # update -
            constraint_counter = constraint_counter + 1
        end
    end

    # cache the additional constraints -
    data_dictionary["additional_constraint_array"] = additional_constraint_array

    # we need to add some additional "species bounds" -
    tmp_species_bounds_array = data_dictionary["species_bounds_array"]
    additional_constraint_block = zeros(number_of_additional_constraints,2)
    species_bounds_array = [tmp_species_bounds_array ; additional_constraint_block]
    data_dictionary["species_bounds_array"] = species_bounds_array

    # iterate through the measurements_dictionary and set values -
    individual_measuerment_dictionaries = measurements_dictionary["exchange_flux_measurements"]
    for (reaction_key, local_measurement_dict) in individual_measuerment_dictionaries

        # ok, so we have a reaction key - find the index of this reaction in the reaction list -
        idx_reaction_match = findall(rxn_name_list .== reaction_key)

        # if we have this reaction, then update the bounds array -
        if (isempty(idx_reaction_match) == false)

            # get the measured value -
            mean_measured_value = parse(Float64,local_measurement_dict["mean_value"])

            # get the CV value -
            coefficient_of_variation = parse(Float64, local_measurement_dict["coefficient_of_variation"])

            # what is the measured value?
            distrubution = Distributions.Normal(mean_measured_value, coefficient_of_variation*mean_measured_value)
            measured_value = abs(rand(distrubution))

            # get the actual index -
            idx_reaction = (getindex(idx_reaction_match))[1]

            # what is the directionality of this measurement?
            directionality = Symbol(local_measurement_dict["directionality"])
            if directionality == :input

                flux_bounds_array[idx_reaction,1] = -1*measured_value
                flux_bounds_array[idx_reaction,2] = 0.0

            elseif directionality == :output

                flux_bounds_array[idx_reaction,1] = 0
                flux_bounds_array[idx_reaction,2] = measured_value

            elseif directionality == :free

                flux_bounds_array[idx_reaction,1] = -1*measured_value
                flux_bounds_array[idx_reaction,2] = measured_value

            elseif directionality == :input_fixed

                flux_bounds_array[idx_reaction,1] = -1*measured_value
                flux_bounds_array[idx_reaction,2] = -1*measured_value

            elseif directionality == :output_fixed

                flux_bounds_array[idx_reaction,1] = measured_value
                flux_bounds_array[idx_reaction,2] = measured_value

            elseif directionality == :nil

                flux_bounds_array[idx_reaction,1] = -1e-6
                flux_bounds_array[idx_reaction,2] = 1e-6

            else
                # TODO: issue a warning ...
            end
        end
    end

    # update -
    data_dictionary["flux_bounds_array"] = flux_bounds_array

    # return modified dictionary -
    return data_dictionary
end

"""
TODO: Fill me in with some stuff ...
"""
function generate_exchange_flux_index_array(reaction_name_array::Array{String,1},pattern::String)

    # initialize -
    exchange_flux_index_array = Int64[]

    # we are looking for a name starting w/EX_
    number_of_reactions = length(reaction_name_array)
    for reaction_index = 1:number_of_reactions

        # grab the name -
        reaction_name = reaction_name_array[reaction_index]

        # check -
        if (startswith(reaction_name,pattern) == true)
            push!(exchange_flux_index_array,reaction_index)
        end
    end

    # return -
    return exchange_flux_index_array
end

"""
TODO: Fill me in with some stuff ...
"""
function generate_default_data_dictionary() #Updted to include rules

    # hardcode paths to stuff (paths internal to package)
    path_to_cobra_mat_file = "$(path_to_package)/cobra/config/matlab_cobra_files/modelReg.mat"
    model_name = "modelReg"
    organism_id = :eco    # we use the KEGG organism symbols

    # Wher are we?
    top_level_path = pwd()

    # load the biophysical_constants dictionary -
    default_biophysical_dictionary = load_default_biophysical_dictionary(organism_id)

    # TODO: check if string is a legit path to the cobra file, and the model name is ok
    # load the *.mat code from the cobra code folder -
    file = matopen(path_to_cobra_mat_file)
    cobra_dictionary = read(file,model_name)
    close(file)

    # Setup: the stoichiometric matrix -
    stoichiometric_matrix = Matrix(cobra_dictionary["S"])
    (number_of_species,number_of_reactions) = size(stoichiometric_matrix)

    # get the species symbol list -
    list_of_metabolite_symbols = cobra_dictionary["mets"]

    # get list of gene symbols -
    list_of_gene_symbols = cobra_dictionary["genes"]

    # setup the objective_coefficient_array -
    objective_coefficient_array = cobra_dictionary["c"]

    #generate list of regulatoryRules and rxn names
    regulatoryRules = cobra_dictionary["regulatoryRules"]
    rxns = cobra_dictionary["rxns"]

    #generate grRules and rules and rxnGeneMat
    grRules = cobra_dictionary["grRules"]
    irules = cobra_dictionary["rules"]
    rxnGeneMat = cobra_dictionary["rxnGeneMat"]

    #Generate regulatoryGenes
    regulatoryGenes = cobra_dictionary["regulatoryGenes"]

    #generate regulatory inputs 1 and 2
    regulatoryInputs1 = cobra_dictionary["regulatoryInputs1"]
    regulatoryInputs2 = cobra_dictionary["regulatoryInputs2"]
    # Setup: default flux bounds array, including the "exhange fluxes"
    default_flux_bounds_array = zeros(number_of_reactions,2)
    lb = cobra_dictionary["lb"] # lower bound -
    ub = cobra_dictionary["ub"] # upper bound -
    for reaction_index = 1:number_of_reactions
        default_flux_bounds_array[reaction_index,1] = lb[reaction_index]
        default_flux_bounds_array[reaction_index,2] = ub[reaction_index]
    end

    # calculate the default vamx (based upon biophysical_constants/literature) -
    (default_vmax, default_enzyme_concentration) = calculate_default_vmax(default_biophysical_dictionary)

    # load the ec_number file -
    path_to_ec_file::String = "$(path_to_package)/cobra/config/data/ec_numbers.dat"
    ec_number_dictionary =  load_ec_mapping_file(path_to_ec_file)

    # load the gene order array -
    path_to_gene_file = "$(path_to_package)/cobra/config/data/gene_list.dat"
    gene_order_array = load_gene_order_file(path_to_gene_file)

    # calculate the vamx array -
    path_to_brenda_data = "$(path_to_package)/cobra/config/data/ECN-EcoliCoreModel-Palsson-EcoSal-Plus-2013.json"
    total_vmax_array = generate_vmax_array(path_to_brenda_data, default_biophysical_dictionary, gene_order_array, ec_number_dictionary)

    # correct the total_vmax_array for this model, using the rules -
    local_data_dictionary = Dict{String,Any}()
    local_data_dictionary["default_vmax_value"] = default_vmax
    model_vmax_array = calculate_rules_vector(local_data_dictionary, total_vmax_array)

    # update the default bounds array w/our "default" biophysical_constants -
    flux_bounds_array = update_default_flux_bounds_array(default_flux_bounds_array, model_vmax_array, cobra_dictionary)
    flux_bounds_array = default_flux_bounds_array #remove this LATER
    # species bounds array - default, everything is bounded 0,0
    # species in the [e] (extracellular) compartment are unbounded
    # we have *updated* the STM w/exchange reactions [] -> a and a -> []
    # all species are bounded to 0
    species_bounds_array = zeros(number_of_species,2)

    # create list of reaction strings -
    list_of_reaction_name_strings = cobra_dictionary["rxns"]

    # What sense do we do? (by default we min)
    is_minimum_flag = true

    # construct the reaction string list -
    list_of_chemical_reaction_strings = reconstruct_reaction_string_list(cobra_dictionary)

    # grab the reversible reaction list -
    reversible_reaction_flag_array = cobra_dictionary["rev"]

    # calculate the reaction name -
    reaction_name_array_tmp = cobra_dictionary["rxns"]
    reaction_name_array = String[]
    for rxn_name in reaction_name_array_tmp
        push!(reaction_name_array,rxn_name)
    end
    exchange_flux_index_array = generate_exchange_flux_index_array(reaction_name_array,"EX_")
    exchange_flux_index_array = [exchange_flux_index_array ; 13]    # add the growth rate -

    # =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{String,Any}()
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
    data_dictionary["number_of_species"] = number_of_species
	data_dictionary["number_of_reactions"] = number_of_reactions
    data_dictionary["flux_bounds_array"] = flux_bounds_array
    data_dictionary["species_bounds_array"] = species_bounds_array
    data_dictionary["objective_coefficient_array"] = objective_coefficient_array
	data_dictionary["list_of_metabolite_symbols"] = list_of_metabolite_symbols
    data_dictionary["list_of_gene_symbols"] = list_of_gene_symbols
    data_dictionary["list_of_reaction_name_strings"] = list_of_reaction_name_strings
    data_dictionary["list_of_chemical_reaction_strings"] = list_of_chemical_reaction_strings
	data_dictionary["is_minimum_flag"] = is_minimum_flag
    data_dictionary["reversible_reaction_flag_array"] = reversible_reaction_flag_array
    data_dictionary["exchange_flux_index_array"] = exchange_flux_index_array

    # stuff for rules -
    data_dictionary["default_vmax"] = default_vmax
    data_dictionary["regulatoryRules"] = regulatoryRules
    data_dictionary["grRules"] = grRules
    data_dictionary["irules"] = irules
    data_dictionary["rxnGeneMat"] = rxnGeneMat
    data_dictionary["regulatoryGenes"] = regulatoryGenes
    data_dictionary["regulatoryInputs1"] = regulatoryInputs1
    data_dictionary["regulatoryInputs2"] = regulatoryInputs2
    data_dictionary["rxns"] = rxns

    # in case we need something later -
    data_dictionary["cobra_dictionary"] = cobra_dictionary
    # ========================================================================================= #
    return data_dictionary
end

# - PRIVATE HELPER METHODS -------------------------------------------------------------------- #
function load_gene_order_file(path_to_gene_file)

    # gene array -
    gene_symbol_array = String[]

    # open the gene file -
    # TODO: We need to check to see if this is a legit path -
    open(path_to_gene_file) do f

        # load the file into the buffer -
        buffer = read(f, String)

        # split along new line -
        list_of_records = split(buffer,'\n')

        # how many records do we have?
        number_of_records = length(list_of_records)
        for record_index = 1:number_of_records

            # convert the record to a string -
            record = string(list_of_records[record_index])

            # is the record empty?
            if (isempty(record) == false)

                # split into fields -
                field_array = split(record,',')

                # grab the ec number -
                order_index = field_array[1]
                gene_symbol = field_array[2]

                # add a record -
                push!(gene_symbol_array, gene_symbol)
            end
        end
    end

    # return -
    return gene_symbol_array
end

function load_ec_mapping_file(path_to_ec_file::String)

    # initialize the ec number array -
    ec_number_dictionary = Dict{String,String}()

    # open the ec file -
    open(path_to_ec_file) do f

        # load the file into the buffer -
        buffer = read(f, String)

        # split along new line -
        list_of_records = split(buffer,'\n')

        # how many records do we have?
        number_of_records = length(list_of_records)
        for record_index = 1:number_of_records

            # convert the record to a string -
            record = string(list_of_records[record_index])

            # is the record empty?
            if (isempty(record) == false)

                # split into fields -
                field_array = split(record,'\t')

                # grab the ec number -
                bg_number = field_array[1]
                ec_number = field_array[2]

                # add a record -
                ec_number_dictionary[bg_number] = ec_number

            end
        end
    end

    # return -
    return ec_number_dictionary
end

function generate_vmax_array(path_to_brenda_data::String, biophysical_dictionary::Dict{String,Any}, gene_symbol_array::Array{String,1}, ec_number_dictionary::Dict{String,String})

    # TODO: is the path_to_brenda_data legit?
    kcat_measurements_dictionary = JSON.parsefile(path_to_brenda_data)

    # initialize -
    vmax_array = Float64[]

    # calculate a default value -
    (default_vmax, default_enzyme_concentration) = calculate_default_vmax(biophysical_dictionary)

    # for each gene symbol, lookup an ec number, then get a range of measured kcats -
    for gene_symbol in gene_symbol_array

        # do we have this gene key?
        if (haskey(ec_number_dictionary, gene_symbol) == true)

            # what is the ec number for this key?
            ec_number_key = ec_number_dictionary[gene_symbol]

            # ok, so we have an ec number. Do we have a kcat measurement for this ec number?
            if (haskey(kcat_measurements_dictionary, ec_number_key) == true)

                # ok, we have an ec_number w/data - get the data
                # this will come back as an array of dictionaries
                local_vmax_array = Float64[]
                array_data_dictionaries = kcat_measurements_dictionary[ec_number_key]
                for local_data_dictionary::Dict{String,Any} in array_data_dictionaries

                    # get the turnoverNumber -
                    turnoverNumber = parse(Float64,local_data_dictionary["turnoverNumber"])*(3600) # hr^-1

                    # there is a "strange" data issue w/the Brenda data. They use -999 to mean no value
                    if (turnoverNumber>0)
                        # cache -
                        push!(local_vmax_array, turnoverNumber*default_enzyme_concentration)
                    end
                end

                if isempty(local_vmax_array) == true
                    push!(local_vmax_array, default_vmax)
                end

                # what is the min/max -
                min_value = minimum(local_vmax_array)
                max_value = maximum(local_vmax_array)

                # calculate a value (uniform random between the measurements)
                sampled_vmax_value = min_value + (max_value - min_value)*rand()

                # yes - we have an experimental value -
                push!(vmax_array, sampled_vmax_value)
            else

                # nope - bummer. Use the default value -
                push!(vmax_array, default_vmax)
            end
        else

            # nope - bummer. Use the default value -
            push!(vmax_array, default_vmax)
        end
    end

    # return -
    return vmax_array
end

function calculate_default_vmax(biophysical_dictionary::Dict{String,Any})

    # ok, so we need to get some stuff from the dictionary -
    # TODO: Check these keys are contained in the dictionary
    default_turnover_number = parse(Float64,biophysical_dictionary["biophysical_constants"]["default_turnover_number"]["value"])              # convert to h^-1
    default_enzyme_concentration = parse(Float64,biophysical_dictionary["biophysical_constants"]["default_enzyme_concentation"]["value"])     # mumol/gDW

    # calculate the default VMax -
    default_vmax = (default_turnover_number)*(default_enzyme_concentration)*(3600)

    return (default_vmax, default_enzyme_concentration)
end

function update_default_flux_bounds_array(flux_bounds_array::Array{Float64,2}, vmax_array::Array{Float64,1}, cobra_dictionary::Dict{String,Any})

    # get the reverse flag array from the cobra dictionary -
    rev_flag_array = cobra_dictionary["rev"]

    # how many bounds do we have?
    (number_of_bounds, number_of_columns) = size(flux_bounds_array)

    # initialize and update the array -
    updated_flux_bounds_array = zeros(number_of_bounds,number_of_columns)
    for flux_index = 1:number_of_bounds

        # get vmax value -
        vmax_value = vmax_array[flux_index]

        # setup new bounds -
        updated_flux_bounds_array[flux_index,1] = 0.0
        updated_flux_bounds_array[flux_index,2] = vmax_value
        if (rev_flag_array[flux_index] == 1.0)
            updated_flux_bounds_array[flux_index,1] = -1*vmax_value
        end
    end

    # return -
    return updated_flux_bounds_array
end

function load_default_biophysical_dictionary(organism_id::Symbol)

    # load the organism_id mapping file -
    path_to_mapping_file = "$(path_to_package)/flux/config/Mapping.json"

    # TODO: Check do we have a mapping file?
    organism_id_map = JSON.parsefile(path_to_mapping_file)

    # TODO: Throw an error if we are missing this organism
    local_path_to_biophysical_dictionary = organism_id_map[String(organism_id)]["path_to_default_constants"]

    # TODO: check - does this path point to a file?
    global_path_to_biophysical_dictionary = "$(path_to_package)"*"$(local_path_to_biophysical_dictionary)"
    default_biophysical_dictionary = JSON.parsefile(global_path_to_biophysical_dictionary)

    # return -
    return default_biophysical_dictionary
end


text:
using ModelingToolkit
using Sundials, Test, SparseArrays

# Comparing solution obtained by defining explicit Jacobian function with solution obtained from
# symbolically generated Jacobian

function testjac(res, du, u, p, t) #System of equations
    res[1] = du[1] - 1.5 * u[1] + 1.0 * u[1] * u[2]
    res[2] = du[2] + 3 * u[2] - u[1] * u[2]
end

function testjac_jac(J, du, u, p, gamma, t) #Explicit Jacobian
    J[1, 1] = gamma - 1.5 + 1.0 * u[2]
    J[1, 2] = 1.0 * u[1]
    J[2, 1] = -1 * u[2]
    J[2, 2] = gamma + 3 - u[1]
    nothing
end

testjac_f = DAEFunction(testjac, jac = testjac_jac, jac_prototype = sparse([1, 2, 1, 2], [1, 1, 2, 2], zeros(4)))

prob1 = DAEProblem(
    testjac_f,
    [0.5, -2.0],
    ones(2),
    (0.0, 10.0),
    differential_vars = [true, true],
)
sol1 = solve(prob1, IDA(linear_solver = :KLU))

# Now MTK style solution with generated Jacobian

@variables t u1(t) u2(t)
@parameters p1 p2

D = Differential(t)

eqs = [D(u1) ~ p1*u1 - u1 * u2,
    D(u2) ~ u1*u2 - p2*u2]

@named sys = ODESystem(eqs)

u0 = [u1 => 1.0,
    u2 => 1.0]

tspan = (0.0, 10.0)

du0 = [0.5, -2.0]

p = [p1 => 1.5,
     p2 => 3.0]

prob = DAEProblem(sys, du0, u0, tspan, p, jac=true, sparse=true)
sol = solve(prob, IDA(linear_solver= :KLU))

@test maximum(sol - sol1) < 1e-12


text:
using Oceananigans.Fields: cpudata, FieldSlicer, interior_copy

"""
    correct_field_size(arch, grid, FieldType, Tx, Ty, Tz)

Test that the field initialized by the FieldType constructor on `arch` and `grid`
has size `(Tx, Ty, Tz)`.
"""
correct_field_size(a, g, FieldType, Tx, Ty, Tz) = size(parent(FieldType(a, g))) == (Tx, Ty, Tz)

function run_similar_field_tests(f)
    g = similar(f)
    @test typeof(f) == typeof(g)
    @test f.grid == g.grid
    return nothing
end

"""
     correct_field_value_was_set(N, L, ftf, val)

Test that the field initialized by the field type function `ftf` on the grid g
can be correctly filled with the value `val` using the `set!(f::AbstractField, v)`
function.
"""
function correct_field_value_was_set(arch, grid, FieldType, val::Number)
    f = FieldType(arch, grid)
    set!(f, val)
    return all(interior(f) .≈ val * arch_array(arch, ones(size(f))))
end

function run_field_reduction_tests(FT, arch)
    N = 8
    topo = (Bounded, Bounded, Bounded)
    grid = RegularRectilinearGrid(FT, topology=topo, size=(N, N, N), x=(-1, 1), y=(0, 2π), z=(-1, 1))

    u = XFaceField(arch, grid)
    v = YFaceField(arch, grid)
    w = ZFaceField(arch, grid)
    c = CenterField(arch, grid)

    f(x, y, z) = exp(x) * sin(y) * tanh(z)

    ϕs = (u, v, w, c)
    [set!(ϕ, f) for ϕ in ϕs]

    u_vals = f.(nodes(u, reshape=true)...)
    v_vals = f.(nodes(v, reshape=true)...)
    w_vals = f.(nodes(w, reshape=true)...)
    c_vals = f.(nodes(c, reshape=true)...)

    # Convert to CuArray if needed.
    u_vals = arch_array(arch, u_vals)
    v_vals = arch_array(arch, v_vals)
    w_vals = arch_array(arch, w_vals)
    c_vals = arch_array(arch, c_vals)

    ϕs_vals = (u_vals, v_vals, w_vals, c_vals)

    dims_to_test = (1, 2, 3, (1, 2), (1, 3), (2, 3))

    for (ϕ, ϕ_vals) in zip(ϕs, ϕs_vals)

        @test all(ϕ .== ϕ_vals) # if this isn't true, reduction tests can't pass

        # Important to make sure no CUDA scalar operations occur!
        CUDA.allowscalar(false)
        @test minimum(ϕ) == minimum(ϕ_vals)
        @test maximum(ϕ) == maximum(ϕ_vals)
        @test mean(ϕ) == mean(ϕ_vals)
        @test minimum(∛, ϕ) == minimum(∛, ϕ_vals)
        @test maximum(abs, ϕ) == maximum(abs, ϕ_vals)
        @test mean(abs2, ϕ) == mean(abs2, ϕ)

        for dims in dims_to_test
            @test minimum(ϕ, dims=dims) == minimum(ϕ_vals, dims=dims)
            @test maximum(ϕ, dims=dims) == maximum(ϕ_vals, dims=dims)
            @test mean(ϕ, dims=dims) == mean(ϕ_vals, dims=dims)

            @test minimum(sin, ϕ, dims=dims) == minimum(sin, ϕ_vals, dims=dims)
            @test maximum(cos, ϕ, dims=dims) == maximum(cos, ϕ_vals, dims=dims)
            @test mean(cosh, ϕ, dims=dims) == mean(cosh, ϕ, dims=dims)
        end
        CUDA.allowscalar(true)
    end

    return nothing
end

function run_field_interpolation_tests(arch, FT)

    grid = RegularRectilinearGrid(size=(4, 5, 7), x=(0, 1), y=(-π, π), z=(-5.3, 2.7))

    velocities = VelocityFields(arch, grid)
    tracers = TracerFields((:c,), arch, grid)

    (u, v, w), c = velocities, tracers.c

    # Choose a trilinear function so trilinear interpolation can return values that
    # are exactly correct.
    f(x, y, z) = exp(-1) + 3x - y/7 + z + 2x*y - 3x*z + 4y*z - 5x*y*z

    # Maximum expected rounding error is the unit in last place of the maximum value
    # of f over the domain of the grid.
    ε_max = f.(nodes((Face, Face, Face), grid, reshape=true)...) |> maximum |> eps

    set!(u, f)
    set!(v, f)
    set!(w, f)
    set!(c, f)

    # Check that interpolating to the field's own grid points returns
    # the same value as the field itself.

    ℑu = interpolate.(Ref(u), nodes(u, reshape=true)...)
    ℑv = interpolate.(Ref(v), nodes(v, reshape=true)...)
    ℑw = interpolate.(Ref(w), nodes(w, reshape=true)...)
    ℑc = interpolate.(Ref(c), nodes(c, reshape=true)...)

    @test all(isapprox.(ℑu, Array(interior(u)), atol=ε_max))
    @test all(isapprox.(ℑv, Array(interior(v)), atol=ε_max))
    @test all(isapprox.(ℑw, Array(interior(w)), atol=ε_max))
    @test all(isapprox.(ℑc, Array(interior(c)), atol=ε_max))

    # Check that interpolating between grid points works as expected.

    xs = reshape([0.3, 0.55, 0.73], (3, 1, 1))
    ys = reshape([-π/6, 0, 1+1e-7], (1, 3, 1))
    zs = reshape([-1.3, 1.23, 2.1], (1, 1, 3))

    ℑu = interpolate.(Ref(u), xs, ys, zs)
    ℑv = interpolate.(Ref(v), xs, ys, zs)
    ℑw = interpolate.(Ref(w), xs, ys, zs)
    ℑc = interpolate.(Ref(c), xs, ys, zs)

    F = f.(xs, ys, zs)

    @test all(isapprox.(ℑu, F, atol=ε_max))
    @test all(isapprox.(ℑv, F, atol=ε_max))
    @test all(isapprox.(ℑw, F, atol=ε_max))
    @test all(isapprox.(ℑc, F, atol=ε_max))

    return nothing
end

@testset "Fields" begin
    @info "Testing Fields..."

    @testset "Field initialization" begin
        @info "  Testing Field initialization..."

        N = (4, 6, 8)
        L = (2π, 3π, 5π)
        H = (1, 1, 1)

        for arch in archs, FT in float_types
            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic))
            @test correct_field_size(arch, grid, CenterField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])

            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded))
            @test correct_field_size(arch, grid, CenterField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3] + 1)

            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded))
            @test correct_field_size(arch, grid, CenterField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField,  N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3])

            grid = RegularRectilinearGrid(FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded))
            @test correct_field_size(arch, grid, CenterField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField,  N[1] + 1 + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField,  N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3])
        end
    end

    @testset "Setting fields" begin
        @info "  Testing field setting..."

        CUDA.allowscalar(true)

        FieldTypes = (CenterField, XFaceField, YFaceField, ZFaceField)

        N = (4, 6, 8)
        L = (2π, 3π, 5π)
        H = (1, 1, 1)

        int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4), Int128(-5)]
        uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10), UInt128(11)]
        float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
        rational_vals = Any[1//11, -23//7]
        other_vals = Any[π]
        vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

        for arch in archs, FT in float_types
            ArrayType = array_type(arch)
            grid = RegularRectilinearGrid(FT, size=N, extent=L, topology=(Periodic, Periodic, Bounded))

            for FieldType in FieldTypes, val in vals
                @test correct_field_value_was_set(arch, grid, FieldType, val)
            end

            for FieldType in FieldTypes
                field = FieldType(arch, grid)
                sz = size(field)
                A = rand(FT, sz...)
                set!(field, A)
                @test field.data[2, 4, 6] == A[2, 4, 6]
            end

            Nx = 8
            topo = (Bounded, Bounded, Bounded)
            grid = RegularRectilinearGrid(FT, topology=topo, size=(Nx, Nx, Nx), x=(-1, 1), y=(0, 2π), z=(-1, 1))

            u = XFaceField(arch, grid)
            v = YFaceField(arch, grid)
            w = ZFaceField(arch, grid)
            c = CenterField(arch, grid)

            f(x, y, z) = exp(x) * sin(y) * tanh(z)

            ϕs = (u, v, w, c)
            [set!(ϕ, f) for ϕ in ϕs]

            @test u[1, 2, 3] == f(grid.xF[1], grid.yC[2], grid.zC[3])
            @test v[1, 2, 3] == f(grid.xC[1], grid.yF[2], grid.zC[3])
            @test w[1, 2, 3] == f(grid.xC[1], grid.yC[2], grid.zF[3])
            @test c[1, 2, 3] == f(grid.xC[1], grid.yC[2], grid.zC[3])
        end
    end

    @testset "Field reductions" begin
        @info "  Testing field reductions..."

        for arch in archs, FT in float_types
            run_field_reduction_tests(FT, arch)
        end
    end

    @testset "Field interpolation" begin
        @info "  Testing field interpolation..."

        for arch in archs, FT in float_types
            run_field_interpolation_tests(arch, FT)
        end
    end

    @testset "Field utils" begin
        @info "  Testing field utils..."

        @test Fields.has_velocities(()) == false
        @test Fields.has_velocities((:u,)) == false
        @test Fields.has_velocities((:u, :v)) == false
        @test Fields.has_velocities((:u, :v, :w)) == true

        grid = RegularRectilinearGrid(size=(4, 6, 8), extent=(1, 1, 1))
        ϕ = CenterField(CPU(), grid)
        @test cpudata(ϕ).parent isa Array

        if CUDA.has_cuda()
            ϕ = CenterField(GPU(), grid)
            @test cpudata(ϕ).parent isa Array
        end

        @test FieldSlicer() isa FieldSlicer

        @info "    Testing similar(f) for f::Union(Field, ReducedField)..."

        grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))

        for X in (Center, Face), Y in (Center, Face), Z in (Center, Face)
            for arch in archs
                f = Field(X, Y, Z, arch, grid)
                run_similar_field_tests(f)

                for dims in (3, (1, 2), (1, 2, 3))
                    f = ReducedField(X, Y, Z, arch, grid, dims=dims)
                    run_similar_field_tests(f)
                end
            end
        end
    end
end


text:
function OutlookClient(email, id, secret)
    if isfile("credentials.json")

        open("credentials.json", "r") do file
            creds = JSON.parse(readstring(file))
            auth = OutlookAuthentication(creds["access_token"], creds["refresh_token"])
            close(file)
            client = OutlookClient(email, id, secret, auth, Array{OutlookFolder, 1}())
            refresh_tokens(client)
            load_folders(client)
            client
        end

    else

        println("Please open $(authurl(id)) in your browser to authenticate this application.")
        client = OutlookClient(email, id, secret)
        start_auth_server(client)
        load_folders(client)
        client
    end
end

function finish_authentication(params, client::OutlookClient)
    code = replace(params, "/?code=", "")
    get_tokens(client, code)
end

clientsecret(client::OutlookClient) = client.clientsecret
clientid(client::OutlookClient)     = client.clientid
refreshtoken(client::OutlookClient) = client.authentication.refreshtoken
accesstoken(client::OutlookClient)  = client.authentication.accesstoken
emailaddress(client::OutlookClient) = client.email
folders(client::OutlookClient)      = client.folders


text:
module Openresty

using Openresty_jll

export OpenrestyCtx
export setup, start, stop, restart, isrunning, reopen, reload

const nginxbindir = joinpath(nginx_dir, "sbin")
const htmltemplatedir = joinpath(nginx_dir, "html")
const conftemplatedir = joinpath(nginx_dir, "conf")
const luapath = joinpath(lualib_dir, "?.lua")
const luacpath = joinpath(lualib_dir, "?.so")

mutable struct OpenrestyCtx
    workdir::String
    sudo::Bool
    pid::Union{Int,Nothing}

    function OpenrestyCtx(workdir; sudo::Bool=false)
        ctx = new(workdir, sudo, nothing)
        readpid(ctx)
        ctx
    end
end

htmldir(ctx::OpenrestyCtx) = joinpath(ctx.workdir, "html")
confdir(ctx::OpenrestyCtx) = joinpath(ctx.workdir, "conf")
conffile(ctx::OpenrestyCtx) = joinpath(confdir(ctx), "nginx.conf")
logsdir(ctx::OpenrestyCtx) = joinpath(ctx.workdir, "logs")
pidfile(ctx::OpenrestyCtx) = joinpath(logsdir(ctx), "nginx.pid")
accesslogfile(ctx::OpenrestyCtx) = joinpath(logsdir(ctx), "access.log")
errorlogfile(ctx::OpenrestyCtx) = joinpath(logsdir(ctx), "error.log")

function setup(ctx::OpenrestyCtx, configfile::Union{String,Nothing}=nothing; force::Bool=false, reset_templates::Bool=force, lua_package_path::Union{String,Vector{String},Nothing}=nothing, lua_package_cpath::Union{String,Vector{String},Nothing}=nothing)
    existing_setup = isdir(confdir(ctx)) && isdir(htmldir(ctx)) && isdir(logsdir(ctx))

    existing_setup && !force && error("setup already exists, specify force=true to overwrite")

    # make the workdir
    for path in (htmldir(ctx), confdir(ctx), logsdir(ctx))
        isdir(path) || mkpath(path)
    end

    # place configuration file
    if force || !isfile(conffile(ctx))
        # copy over bundled base configurations
        cp(conftemplatedir, confdir(ctx); force=true)
        # copy over provided configuration
        (configfile !== nothing) && cp(configfile, conffile(ctx); force=true)
        set_lua_package_path(ctx, lua_package_path, lua_package_cpath)
    end

    # copy over bundled templates
    (reset_templates || !existing_setup) && cp(htmltemplatedir, htmldir(ctx); force=true)

    nothing
end

"""
Replace all occurrences of OPENRESTY_LUA_PACKAGE_PATH in the configuration file by
actual lua package path
"""
function set_lua_package_path(ctx::OpenrestyCtx, lua_package_path::Union{String,Vector{String},Nothing}=nothing, lua_package_cpath::Union{String,Vector{String},Nothing}=nothing)
    config = read(conffile(ctx), String)
    all_lua_paths = [luapath]
    all_lua_cpaths = [luacpath]
    if isa(lua_package_path, String)
        push!(all_lua_paths, lua_package_path)
    elseif isa(lua_package_path, Vector{String})
        append!(all_lua_paths, lua_package_path)
    end
    if isa(lua_package_cpath, String)
        push!(all_lua_cpaths, lua_package_cpath)
    elseif isa(lua_package_cpath, Vector{String})
        append!(all_lua_cpaths, lua_package_cpath)
    end
    pathstr = join(all_lua_paths, ';') * ";;"
    cpathstr = join(all_lua_cpaths, ';') * ";;"

    config = replace(config, "OPENRESTY_LUA_PACKAGE_PATH" => pathstr)
    config = replace(config, "OPENRESTY_LUA_PACKAGE_CPATH" => cpathstr)
    open(conffile(ctx), "w") do f
        write(f, config)
    end
    nothing
end

function start(ctx::OpenrestyCtx; accesslog=nothing, errorlog=nothing, append::Bool=(isa(accesslog,AbstractString)||isa(errorlog,AbstractString)))
    proc = openresty() do openresty_path
        @debug("starting", openresty_path, workdir=ctx.workdir, nginxbindir, sudo=ctx.sudo)
        command = Cmd(ctx.sudo ? `sudo $openresty_path -p $(ctx.workdir)` : `$openresty_path -p $(ctx.workdir)`; detach=true, dir=nginxbindir)
        if (accesslog === nothing) && (errorlog === nothing)
            run(command; wait=false)
        else
            redirected_command = pipeline(command, stdout=accesslog, stderr=errorlog, append=append)
            run(redirected_command; wait=false)
        end
    end

    # Wait until nginx comes up on a slow system.
    # Does not throw error though if a running nginx process is not detected,
    # leaves it to the caller to do the final check
    timedwait(15.0; pollint=1.0) do
        readpid(ctx)
        !((ctx.pid === nothing) && process_running(proc))
    end
    nothing
end

isrunning(ctx::OpenrestyCtx) = (ctx.pid !== nothing) ? isrunning(ctx, ctx.pid) : false
function isrunning(ctx::OpenrestyCtx, pid::Int)
    # we do have read permission on cmdline even if process was started with sudo
    cmdlinefile = "/proc/$pid/cmdline"
    openresty_path = openresty(p->p)
    if isfile(cmdlinefile)
        cmdline = read(cmdlinefile, String)
        @debug("found command line", cmdlinefile, cmdline , openresty_path, ctx.workdir)
        if occursin(openresty_path, cmdline) && occursin(ctx.workdir, cmdline)
            # process still running
            return true
        end
    end

    ctx.pid = nothing
    rm(pidfile(ctx); force=true)
    false
end

function stop(ctx::OpenrestyCtx; grace_seconds::Int=2)
    if isrunning(ctx)
        signalquit(ctx)
        sleep(grace_seconds)
        isrunning(ctx) && signalstop(ctx)
        while isrunning(ctx)
            @debug("waiting for shutdown...")
            sleep(grace_seconds+1)
        end
    end
    nothing
end

"""
Do a full stop and start. Disrupts the service while nginx restarts and also for the duration specified in `delay_seconds` seconds.
Logs are redirected afresh.
"""
function restart(ctx::OpenrestyCtx; delay_seconds::Int=0, accesslog=nothing, errorlog=nothing, append::Bool=(isa(accesslog,AbstractString)||isa(errorlog,AbstractString)))
    stop(ctx)
    (delay_seconds > 0) && sleep(delay_seconds)
    start(ctx; accesslog=accesslog, errorlog=errorlog, append=append)
end

"""
Reopen log files without disruption.
"""
reopen(ctx::OpenrestyCtx) = signalreopen(ctx)

"""
Reload any configuration changes without disruption.
"""
reload(ctx::OpenrestyCtx) = signalreload(ctx)

signalstop(ctx::OpenrestyCtx) = signal(ctx, "stop")
signalquit(ctx::OpenrestyCtx) = signal(ctx, "quit")
signalreopen(ctx::OpenrestyCtx) = signal(ctx, "reopen")
signalreload(ctx::OpenrestyCtx) = signal(ctx, "reload")
function signal(ctx::OpenrestyCtx, signal::String)
    if isrunning(ctx)
        @debug("sending signal $signal")
        openresty() do openresty_path
            command = Cmd(ctx.sudo ? `sudo $openresty_path -p $(ctx.workdir) -s $signal` : `$openresty_path -p $(ctx.workdir) -s $signal`; dir=nginxbindir)
            run(command)
        end
    end
    nothing
end

function readpid(ctx::OpenrestyCtx)
    pfile = pidfile(ctx)
    if isfile(pfile)
        ctx.pid = parse(Int, read(pidfile(ctx), String))
        isrunning(ctx)
    else
        ctx.pid = nothing
    end
    nothing
end

end


text:
using MetNets
using Test
# using Random
using Serialization
# using ProgressMeter

@testset "MetNets.jl" begin
    # Write your tests here.
    include("Utils_tests/Utils_tests.jl")
end


text:
function knn!{T}(metric::Function, neighbors::NeighborList{T}, node::Node{T}, point_index::Int, τ::T)
    node.index == 0 && return τ

    d = metric(node.index, point_index)::T

    if d <= τ && d > 0.0001#&& node.index != point_index
        τ = insert_neighbor!(neighbors, node.index, d)
    end

    node.isleaf && return τ

    μ = node.radius
    if τ + d <= μ
        τ = knn!(metric, neighbors, node.inside, point_index, τ)
    else
        τ = knn!(metric, neighbors, node.outside, point_index, τ)
        if d <= τ + μ
            τ = knn!(metric, neighbors, node.inside, point_index, τ)
        end
    end

    return τ
end

function knn!{T}(neighbors::NeighborList{T}, tree::VPTree, point_index::Int)
    clear!(neighbors)
    knn!(tree.metric, neighbors, tree.root, point_index, T(Inf))
    return neighbors
end

function knn!(tree::VPTree, point_index::Int, k::Int)
    T = typeof(tree.metric(1, 1))
    neighbors = NeighborList{T}(k)
    knn!(neighbors, tree, point_index)
end


text:
using Test
using Statistics

baseline_length_s = 10
fsamp_hz = 10
fknee_hz = [0.01, 0.01]
σ_k = [0.1, 100]
slope = [-1,-1]

chunks = [[datachunk(1, 1, 1000, 1000), datachunk(2, 1, 1000, 1000)]]
baselines_per_process = 2000
rank = 0
comm = nothing

noise =  generate_noise_mpi(
    chunks,
    baselines_per_process,
    baseline_length_s,
    total_time,
    fsamp_hz,
    σ_k,
    fknee_hz,
    slope;
    rank = rank,
    comm = comm,
    input_seed = 1234,
)

@test length(noise) == fsamp_hz * baseline_length_s*baselines_per_process
@test std(noise[1:100000]) ≈ σ_k[1] rtol=0.01
@test std(noise[end-100000: end]) ≈ σ_k[2] rtol=0.01


text:
# a tutorial adapted from the fipy diffusion 1D example
# see: http://www.ctcms.nist.gov/fipy/examples/diffusion/index.html
# to run this file, first load JFVM and GR packages
using JFVM, GR, SpecialFunctions

function diffusion_tutorial()
    # define the domain
    L = 5  # domain length
    Nx = 100 # number of cells
    meshstruct = createMesh1D(Nx, L)
    BC = createBC(meshstruct) # all Neumann boundary condition structure
    BC.left.a[:] .= 0.0
    BC.left.b[:] .=1.0
    BC.left.c[:] .=1.0 # left boundary
    BC.right.a[:] .= 0.0
    BC.right.b[:] .=1.0
    BC.right.c[:] .=0.0 # right boundary
    x = meshstruct.cellcenters.x
    ## define the transfer coeffs
    D_val = 1.0
    alfa = createCellVariable(meshstruct, 1)
    Dave = createFaceVariable(meshstruct, D_val)
    ## define initial values
    c_old = createCellVariable(meshstruct, 0, BC) # initial values
    c = c_old
    ## loop
    dt = 0.001 # time step
    final_t = 1.5
    c_analytical=zeros(size(x))
    for t=dt:dt:final_t
        # step 1: calculate divergence term
        RHS = divergenceTerm(Dave*gradientTerm(c_old))
        c = solveExplicitPDE(c_old, dt, RHS, BC, alfa)
        # step 2: calculate the new value for internal cells
        # c_old_int=internalCells(c_old)
        # c_val_int = dt*excludeGhostRHS(meshstruct, RHS+constantSourceTerm(c_old))+c_old_int[:]
        c_analytical .= 1.0 .- erf.(x/(2*sqrt(D_val*t)))
        # GR.plot(x, [c.value[2:end-1] c_analytical])
        c_old = c
        println(t)
    end
    GR.plot(x, [c.value[2:end-1] c_analytical])
end

diffusion_tutorial()

max_stars_repo_name: simulkade/JFVM

max_forks_repo_name: simulkade/JFVM


text: print(if true 1 else 2)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct ShaderCompilationError <: Exception
    msg
end

Base.showerror(io::IO, err::ShaderCompilationError) = print(io, "ShaderCompilationError:\n\n$(err.msg)")

"""
    compile(shader)

Compile a shader file in text format to SPIR-V.
"""
function compile(source::ShaderSource; extra_flags=[], validate_spirv=true)::ShaderSource
    flags = ["-V"; "--stdin"; extra_flags]
    source.language == SPIR_V && error("Shader $source already compiled")
    source.language == HLSL && push!(flags, "-D")
    if validate_spirv && "--spirv-val" ∉ flags
        push!(flags, "--spirv-val")
    end

    if isempty(source.entry_points)
        entry_point = :main
    elseif length(source.entry_points) == 1
        entry_point = first(source.entry_points)
    else
        error("Multiple entry points not supported in uncompiled source")
    end

    dst = string(tempname(), ".spv")
    input = IOBuffer()
    write(input, source.code)
    seekstart(input)
    err = IOBuffer()

    append!(flags, [
        "-e",
        string(entry_point),
        "-S",
        string(file_ext(source.language, source.stage)),
        "-o",
        dst,
    ])
    try
        run(pipeline(`$glslangValidator $flags`, stdin=input, stdout=err))
    catch e
        if e isa ProcessFailedException
            err_str = String(take!(err))
            throw(ShaderCompilationError(err_str))
        else
            rethrow(e)
        end
    end

    code = read(dst)
    rm(dst)
    ShaderSource(code, SPIR_V, source.stage, [entry_point])
end

function Vk.ShaderModule(device, source::ShaderSource)
    length(source.code) % 4 == 0 || pad_shader_code!(source.code)
    Vk.ShaderModule(device, length(source.code), reinterpret(UInt32, source.code))
end

function Vk.PipelineShaderStageCreateInfo(shader::Shader)
    specialization_info = isempty(shader.specialization_constants) ? C_NULL : shader.specialization_constants
    Vk.PipelineShaderStageCreateInfo(shader.source.stage, shader.shader_module, string(shader.entry_point); specialization_info)
end