text:
######################################
#--------------OCCASION--------------#
######################################

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

# set response variable
y_var <- "Trps_POS_P1"

# keep only the buyers
finaldata_occ <- initial_data[Buyer_Pos_P1=="1",]

# drop variables with only one level
indices <- sapply(finaldata_occ,function(x) length(unique(x))==1)
finaldata_occ <- finaldata_occ[,!indices,with=F]

#######################################
#-------- Variable Selection ---------#
#######################################

# remove dolocc outliers
finaldata_occ <- as.data.frame(finaldata_occ)
if(exists('dolocc_outliers_panid_to_remove', envir = environment())){
  finaldata_occ <- finaldata_occ[, !(names(finaldata_occ) %in% c("Dol_per_Trip_POS_P1","Nonbuyer_Pre_P1","Buyer_Pre_P1","Dol_per_Trip_PRE_P1"))]
  finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% dolocc_outliers_panid_to_remove,]
}else{
  finaldata_occ <- finaldata_occ[, !(names(finaldata_occ) %in% c("Dol_per_Trip_POS_P1","Nonbuyer_Pre_P1","Buyer_Pre_P1","Dol_per_Trip_PRE_P1"))]
}

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

# run glm model
model_glm_total <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm_total)

##### Variable selection - Condition 1 #####
# Remove pvalues greater than 0.7 and rerun the glm model
pval_excluded_vars <- 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)))
finaldata_occ <- finaldata_occ[, !names(finaldata_occ) %in% c(pval_excluded_vars)]

# re-run glm model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"Trps_PRE_P1","panid",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1))")
model_glm_condition1 <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
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_occ <- finaldata_occ[ , !names(finaldata_occ) %in% multicols_toremove]
  Formula_fixed <-  paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"panid",random_demos,random_campaigns)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1)) ")
  model_glm_new <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
} else {
  model_glm_new <- model_glm_condition1
}

vif_output <- (vif(model_glm_new))
vif_to_remove <- names(vif_output[vif_output > 15])

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

# remove the variables under the above conditions
finaldata_occ <- finaldata_occ[ , !names(finaldata_occ) %in% c(vars_condition2,pval_excluded_vars)]
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"panid","Trps_PRE_P1",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1))")

# run model
model_glm <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm)
rm(model_glm_condition1)

# p-values and signs checks
neutralvars <- setdiff(names(finaldata_occ),c(positivevars,negativevars))
pvalues <- summary(model_glm)$coefficients[,4][-1]
keep_covariate <- names(pvalues[pvalues<pvalue_lvl])
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(Trps_PRE_P1+1))")

model_glm_2 <-  glm(Formula_final,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm_2)
rm(model_glm)

# final 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 <- c();vars_2 <- c()
for (i in negativevars) if (i %in% names(coeffs_2[coeffs_2 < 0])) varstokeep <- c(varstokeep,i)
for (j in positivevars) if (j %in% names(coeffs_2[coeffs_2 > 0])) varstokeep <- c(varstokeep,j)
for (l in neutralvars) varstokeep <- c(varstokeep,l)
for (k in varstokeep) if (k %in% keep_covariate_2) vars_2 <- c(vars_2,k)

# final correlation check (find pair vars with correlation greater then 0.8 and drop the one with bigger pvalue) #### Maybe 0.6!
if(length(vars_2)>1){
  finaldata_occ_num <- as.data.table(finaldata_occ[,names(finaldata_occ) %in% vars_2])
  t1 <- as.data.frame(finaldata_occ_num[,sapply(finaldata_occ_num,is.numeric)])
  finaldata_occ_num <- as.data.frame(finaldata_occ_num)
  finaldata_occ_num <- finaldata_occ_num[,c(which(t1==T))]
  temp_cor <- as.matrix(cor(finaldata_occ_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(model_glm_2)$coefficients[which(rownames(summary(model_glm_2)$coefficients)==as.character(cordata[i][,1,with=F])),4] > summary(model_glm_2)$coefficients[which(rownames(summary(model_glm_2)$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]))
        }
      }
    }
    vars_2 <- setdiff(vars_2,drop_cor_vars)
  }
}

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

#set the formula for the final model
Formula_final <-paste0( paste0(y_var,"~ ",paste0(vars_2,collapse=" + ")," "),"+ group + offset(log(Trps_PRE_P1+1))")
occ_final_model <-  glm(Formula_final,data=finaldata_occ,family=poisson(link="log"))
summary(occ_final_model)
rm(model_glm_2)

# final model
# final cleaning of model based on signs
pvalues_2 <- summary(occ_final_model)$coefficients[,4][-1]
keep_covariate_2 <- names(pvalues_2[pvalues_2 < pvalue_lvl])
coeffs <- sign(summary(occ_final_model)$coefficients[,1])
finalvars <- c();finalvars1 <- 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) finalvars1 <- c(finalvars1,l)
for (k in finalvars1) if (k %in% keep_covariate_2) finalvars <- c(finalvars,k)

#set final model
Formula_final <- paste0(paste0(y_var,"~ ",paste0(finalvars,collapse=" + ")," "),"+ group + offset(log(Trps_PRE_P1+1))")
occ_final_model <-  glm(Formula_final,data=finaldata_occ,family=poisson(link="log"))
summary(occ_final_model)

options(warn=0)

finaldata_occ <- droplevels(finaldata_occ)
occ_random_campaigns <- intersect(random_campaigns,names(finaldata_occ))

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

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

# run glmer model
occ.mixed.start.time <- Sys.time()
Formula_final <- paste0(c(paste0(paste0(y_var,"~ ",paste0(finalvars,collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1)) + group"),random_demos_factors,occ_random_campaigns_factors),collapse='')
occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
occ.mixed.model.execution.time <- Sys.time()-occ.mixed.start.time
out_resid <- order(residuals(occ_final_model),decreasing = T)
finaldata_occ <- as.data.table(finaldata_occ)
occ_total_outliers_panid <- finaldata_occ[out_resid,panid]
segment_panid <- data_NB_NE_B[,panid]
occ_outliers_panid_segm <- intersect(occ_total_outliers_panid,segment_panid)

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

# remove outliers if needed (negative estimation for exposure)
occ_outliers_panid_to_remove <- c()
occ_outliers_panid_segm_new <- c()
occ_outliers_panid_to_remove_threshold <- c()

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

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

if(length(occ_random_campaigns)>0){
  min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  max_random_coeff <- max(unlist(lapply(ranef(occ_final_model)[occ_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(occ_final_model),decreasing = T)
  out_resid_negative <- order(residuals(occ_final_model),decreasing = F)
  occ_outliers_panid_positive_extreme_panid <- finaldata_occ[out_resid_positive,panid]
  occ_outliers_panid_negative_extreme_panid <- finaldata_occ[out_resid_negative,panid]
  for(i in 1:length(occ_random_campaigns)){
    min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i],FUN=max)))
      random_frame <- as.data.frame(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
      positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i]),'.group1')))== max_random_coeff]$level
      negative_extreme_panid <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
      positive_extreme_panid  <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
      occ_outliers_panid_segm1 <- intersect(occ_outliers_panid_negative_extreme_panid,negative_extreme_panid)
      occ_outliers_panid_segm1 <- occ_outliers_panid_segm1[1:round(0.02*(length(negative_extreme_panid)))]
      if(max_random_coeff > 0.05) {
        occ_outliers_panid_segm2 <- intersect(occ_outliers_panid_positive_extreme_panid,positive_extreme_panid)
        occ_outliers_panid_segm2 <- occ_outliers_panid_segm2[1:round(0.02*(length(positive_extreme_panid)))]
      } else {
        occ_outliers_panid_segm2 <- c()
      }
      occ_outliers_panid_segm_new <- c(occ_outliers_panid_segm1,occ_outliers_panid_segm2)
      finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_segm_new,]
    }
  }
  occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
  coef_flag <- c()
  for(j in flag){
    coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
  }
  min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
    for(i in 1:length(occ_random_campaigns)){
      min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i],FUN=max)))
        random_frame <- as.data.frame(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
        positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i]),'.group1')))== max_random_coeff]$level
        negative_extreme_panid <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
        positive_extreme_panid  <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
        occ_outliers_panid_segm1 <- intersect(occ_outliers_panid_negative_extreme_panid,negative_extreme_panid)
        occ_outliers_panid_segm1 <- occ_outliers_panid_segm1[1:round(0.05*(length(negative_extreme_panid)))]
        if(max_random_coeff > 0.05) {
          occ_outliers_panid_segm2 <- intersect(occ_outliers_panid_positive_extreme_panid,positive_extreme_panid)
          occ_outliers_panid_segm2 <- occ_outliers_panid_segm2[1:round(0.05*(length(positive_extreme_panid)))]
        } else {
          occ_outliers_panid_segm2 <- c()
        }
        occ_outliers_panid_segm_new <- c(occ_outliers_panid_segm1,occ_outliers_panid_segm2)
        finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_segm_new,]
      }
    }
    occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  }
  if(min_random_coeff < -0.05 & min(coef_flag) + min_random_coeff < 0){
    for(i in 1:length(occ_random_campaigns)){
      min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i],FUN=max)))
        random_frame <- as.data.frame(ranef(occ_final_model)[occ_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(occ_final_model)[occ_random_campaigns][i]),'.group1')))==min_random_coeff,]$level
        positive_extreme_levels <- random_table[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i]),'.group1')))== max_random_coeff]$level
        negative_extreme_panid <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% negative_extreme_levels,panid] # -13%
        positive_extreme_panid  <- finaldata_occ[eval(parse(text=paste0(names(ranef(occ_final_model)[occ_random_campaigns][i])))) %in% positive_extreme_levels,panid] # +12%
        occ_outliers_panid_segm1 <- intersect(occ_outliers_panid_negative_extreme_panid,negative_extreme_panid)
        occ_outliers_panid_segm1 <- occ_outliers_panid_segm1[1:round(0.10*(length(negative_extreme_panid)))]
        if(max_random_coeff > 0.05) {
          occ_outliers_panid_segm2 <- intersect(occ_outliers_panid_positive_extreme_panid,positive_extreme_panid)
          occ_outliers_panid_segm2 <- occ_outliers_panid_segm2[1:round(0.10*(length(positive_extreme_panid)))]
        } else {
          occ_outliers_panid_segm2 <- c()
        }
        occ_outliers_panid_segm_new <- c(occ_outliers_panid_segm1,occ_outliers_panid_segm2)
        finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_segm_new,]
      }
    }
    occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=Gamma(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  }
}

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

if(length(occ_random_campaigns)>0){
  min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  max_random_coeff <- max(unlist(lapply(ranef(occ_final_model)[occ_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 <- occ_reduction[,panid]
# outliers_panid <- intersect(occ_total_outliers_panid,segment_panid)
# occ_outliers_panid_to_remove_threshold <- outliers_panid[c(1:round(0.10*(length(outliers_panid))))]
# finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_to_remove_threshold,]
# occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))

if(min(coef_flag) + min_random_coeff < 0){
  cat("removing outliers ...\nEstimated Iteration Execution ");print(occ.mixed.model.execution.time)
  cat("\nIteration 1 ... Start Time: ");print(Sys.time())
  occ_outliers_panid_to_remove <- occ_outliers_panid_segm[1:round(0.025*(nrow(data_NB_NE_B)))]
  finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_to_remove,]
  occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
  coef_flag <- c()
  for(j in flag){
    coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
  }
  min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
  if(min(coef_flag) + min_random_coeff < 0){
    occ_outliers_panid_to_remove <- occ_outliers_panid_segm[1:round(0.05*(nrow(data_NB_NE_B)))]
    finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_to_remove,]
    cat("\nIteration 2 ... Start Time: ");print(Sys.time())
    occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
    coef_flag <- c()
    for(j in flag){
      coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
    }
    min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
    if(min(coef_flag) + min_random_coeff < 0){
      occ_outliers_panid_to_remove <- occ_outliers_panid_segm[1:round(0.10*(nrow(data_NB_NE_B)))]
      finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_to_remove,]
      cat("\nIteration 3 ... Start Time: ");print(Sys.time())
      occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
      coef_flag <- c()
      for(j in flag){
        coef_flag <- append(coef_flag,coef(summary(occ_final_model))[j,"Estimate"])
      }
      min_random_coeff <- min(unlist(lapply(ranef(occ_final_model)[occ_random_campaigns],FUN=min)))
      if(min(coef_flag) + min_random_coeff < 0){
        occ_outliers_panid_to_remove <- occ_outliers_panid_segm[1:round(0.15*(nrow(data_NB_NE_B)))]
        finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% occ_outliers_panid_to_remove,]
        cat("\nFinal Iteration ... Start Time: ");print(Sys.time())
        occ_final_model <- glmer(Formula_final,data=finaldata_occ,family=poisson(link="log"),control = glmerControl(calc.derivs = FALSE, optimizer="nloptwrap2"))
      }
    }
  }
}

occ_dolocc_outliers_panid_to_remove <- unique(c(dolocc_outliers_panid_to_remove,occ_outliers_panid_to_remove,occ_outliers_panid_segm_new,occ_outliers_panid_to_remove_threshold))

write.csv(as.data.frame(summary(occ_final_model)$coefficients),file=paste0(initial_Occ_modelling_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(initial_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(initial_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(initial_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(initial_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(initial_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(initial_Occ_modelling_output,'Occ_random_effects.csv'), row.names = FALSE)

# create a txt file with the summary of the glmer model
logsavefname <- "Occasion_Initial_Results.txt"
sink(paste0(initial_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(initial_Occ_modelling_output,"Occ_Caterpillar_Plots ",i,".pdf"),width=7,height=5)
  print(ggCaterpillar(ranef(occ_final_model, condVar=TRUE)[i]))
  dev.off()
}

# save Occasion final dataset and final model
Occ_Random_Formula_final <- Formula_final
rm(finaldata_occ_num,temp_cor,cordata)

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

Occasion.Execution.Time <- Sys.time()-occ.start.time
cat(paste0("Occasion.Execution.Time: ",Occasion.Execution.Time))


text:
# Standard composite operators

struct ProductOperator{B <: BasisSet, O1 <: LinearOperator{B}, O2 <: LinearOperator{B}} <: LinearOperator{B}
    op1 :: O1
    op2 :: O2

    function ProductOperator(o1 :: LinearOperator, o2 :: LinearOperator)
        @assert basis(o1) == basis(o2)
        B = typeof(basis(o1))
        new{B,typeof(o1),typeof(o2)}(o1, o2)
    end
end
basis(op::ProductOperator) = basis(op.op1)
apply(i::Integer, op::ProductOperator, j::Integer) =
    sum(apply(i, op.op1, k)*apply(k, op.op2, j) for k = 1:length(basis(op)))

struct SumOperator{B <: BasisSet, O1 <: LinearOperator{B}, O2 <: LinearOperator{B}} <: LinearOperator{B}
    op1 :: O1
    op2 :: O2

    function SumOperator(o1 :: LinearOperator, o2 :: LinearOperator)
        @assert basis(o1) == basis(o2)
        B = typeof(basis(o1))
        new{B,typeof(o1),typeof(o2)}(o1, o2)
    end
end
basis(op::SumOperator) = basis(op.op1)
apply(i::Integer, op::SumOperator, j::Integer) = apply(i, op.op1, j) + apply(i, op.op2, j)

struct ScalarMultiplicationOperator{S <: Number, B <: BasisSet, O <: LinearOperator{B}} <: LinearOperator{B}
    α :: S
    op :: O

    function ScalarMultiplicationOperator(α :: S, op :: LinearOperator) where {S <: Number}
        B = typeof(basis(op))
        new{S,B,typeof(op)}(α, op)
    end
end
basis(op::ScalarMultiplicationOperator) = basis(op.op)
apply(i::Integer, op::ScalarMultiplicationOperator, j::Integer) = op.α * apply(i, op.op, j)

# Extending addition and multiplication
Base.:(+)(op1::LinearOperator{B}, op2::LinearOperator{B}) where {B <: BasisSet} = SumOperator(op1, op2)
Base.:(-)(op::LinearOperator) = ScalarMultiplicationOperator(-1, op)
Base.:(-)(op1::LinearOperator{B}, op2::LinearOperator{B}) where {B <: BasisSet} = SumOperator(op1, -op2)
Base.:(*)(op1::LinearOperator{B}, op2::LinearOperator{B}) where {B <: BasisSet} = ProductOperator(op1, op2)
Base.:(*)(α::Number, op::LinearOperator) = ScalarMultiplicationOperator(α, op)
Base.:(*)(op::LinearOperator, α::Number) = ScalarMultiplicationOperator(α, op)
Base.:(/)(op::LinearOperator, α::Number) = ScalarMultiplicationOperator(1/α, op)
function Base.:(^)(op::LinearOperator, n::Integer)
    @assert n > 0
    n == 1 && return op
    ProductOperator(op, op^(n-1))
end

ext: jl


text:
EXAMPLES_DIR = "../examples/"

function run_example(replace_strings, example_name, module_suffix="")
    example_filepath = joinpath(EXAMPLES_DIR, example_name * ".jl")
    txt = read(example_filepath, String)

    for strs in replace_strings
        txt = replace(txt, strs[1] => strs[2])
    end

    test_script_filepath = example_name * "_example_test.jl"

    open(test_script_filepath, "w") do f
        write(f, "module Test_$example_name" * "_$module_suffix\n")
        write(f, txt)
        write(f, "\nend # module")
    end

    try
        include(test_script_filepath)
    catch err
        @error sprint(showerror, err)
        rm(test_script_filepath)
        return false
    end

    rm(test_script_filepath)
    return true
end

@testset "Examples" begin
    @info "Testing examples..."

    @testset "One-dimensional diffusion" begin
        @info "  Testing one-dimensional diffusion example"

        replace_strings = [
            ("size = (1, 1, 128)", "size = (1, 1, 16)"),
            ("stop_iteration = 1000", "stop_iteration = 1"),
            ("stop_iteration += 100", "stop_iteration += 1"),
            ("for i=1:100", "for i=1:1"),
            ("mp4(", "# mp4(")
        ]

        @test_skip run_example(replace_strings, "one_dimensional_diffusion")
    end

    @testset "Two-dimensional turbulence example" begin
        @info "  Testing two-dimensional turbulence example"

        replace_strings = [
            ("N=(128, 128, 1)", "N=(16, 16, 1)"),
            ("for i=1:100", "for i=1:1"),
            ("stop_iteration += 10", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test_skip run_example(replace_strings, "two_dimensional_turbulence")
    end

    for arch in archs
        @testset "Wind and convection mixing example [$(typeof(arch))]" begin
            @info "  Testing wind and convection-driving mixing example [$(typeof(arch))]"

            replace_strings = [
                ("Nz = 48", "Nz = 16"),
                ("progress_frequency=10", "progress_frequency=1"),
                ("for i in 1:100", "for i in 1:1"),
                ("stop_iteration += 10", "stop_iteration += 1"),
                ("mp4(", "# mp4(")
            ]

            if arch == GPU()
                push!(replace_strings, ("architecture = CPU()", "architecture = GPU()"))
            end

            @test_skip run_example(replace_strings, "ocean_wind_mixing_and_convection", string(typeof(arch)))

            rm("ocean_wind_mixing_and_convection.jld2", force=true)
        end
    end

    @testset "Ocean convection with plankton example" begin
        @info "  Testing ocean convection with plankton example"

        replace_strings = [
            ("Nz = 128", "Nz = 16"),
            ("progress_frequency=100", "progress_frequency=1"),
            ("for i = 1:100", "for i = 1:1"),
            ("stop_iteration += 100", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test_skip run_example(replace_strings, "ocean_convection_with_plankton")
    end

    @testset "Internal wave example" begin
        @info "  Testing internal wave example"

        replace_strings = [
            ("Nx = 128", "Nx = 16"),
            ("progress_frequency = 20", "progress_frequency = 1"),
            ("for i=1:100", "for i=1:1"),
            ("stop_iteration += 20", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test_skip run_example(replace_strings, "internal_wave")
    end

    @testset "Eady turbulence" begin
        @info "  Testing Eady turbulence example"

        replace_strings = [
            ("Nh = 64", "Nh = 16"),
            ("Nz = 32", "Nz = 16"),
            ("end_time = 3day", "end_time = 1")
        ]

        @test_skip run_example(replace_strings, "eady_turbulence")
    end
end


text:
function PVBuilder(arg0::CorrelatedRandomVectorGenerator, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::ObservableSatellite)
    return PVBuilder((CorrelatedRandomVectorGenerator, jdouble, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4)
end

function build(obj::PVBuilder, arg0::Vector{SpacecraftState})
    return jcall(obj, "build", PV, (Vector{SpacecraftState},), arg0)
end


text:
module test_julia_symbols

using Test

a,b = 1,2
@test a:b === 1:2
@test (:) isa Colon
@test 1:2 isa UnitRange{Int}
@test  :b ===  :b

@test Symbol("F") === :F
@test Symbol(".") === :.
@test Symbol(":") === :(:)
@test Symbol("::") === :(::)

@test :(:)   isa Symbol
@test :(::)  isa Symbol
@test :(:::) isa Expr
@test :(:::) == Expr(:(::), :(:))

end # module test_julia_symbols


text:
# Consensus and fine splitting functions

"""
   fine_clustering(kmer_vecs; rough_radius::Float64 = 0.01, fine_radius::Float64 = 1.0, clust_multiplier = 1.0,
                        min_devs = 6, initial_devs = 20, user_min_clust = 10, distfunc = euclidean, center = mean,
                        verbose::Int64 = 1, cycle_lim = 20, triangle::Bool = false, p_value::Float64 = 0.01)

Preforms fine clustering on the already rough clustered sequences.

"""

function fine_clustering(kmer_vecs; rough_radius::Float64=0.01, fine_radius::Float64=1.0, clust_multiplier=1.0,
                        min_devs=6, initial_devs=20, user_min_clust=10, distfunc=euclidean, center=mean, verbose::Int64=1,
                        cycle_lim=20, triangle::Bool=false, p_value::Float64=0.01)

    if verbose > 0
        println("Your inputs have begun rough clustering...")
    end
    if verbose == 1 ||verbose == 2
        v = 2
    elseif verbose == 0
        v = 0
    end
    #"rough" clustering step to get the massively different sequences into different clusters
    μs, sizes, indices, centroid_inds = cluster(kmer_vecs, rough_radius, distfunc=corrected_kmer_dist,
                                        center=center,verbose=v, cycle_lim=cycle_lim, triangle=triangle)
    #only sequences greater than indicated size may form a cluster
    keep_inds = find(user_min_clust.<sizes)
    indices = indices[keep_inds]

    fine_indices = Array{Array{Int64, 1}, 1}()
    #get array which converts kmer indices into kmer sequences
    nuc_inds = gen_arrangements(Int(log(4, length(kmer_vecs[1]))))

    if verbose > 0
        println("Your inputs have begun fine clustering...")
    end
    if verbose == 1 || verbose == 0
        in_verbose = 0
    elseif verbose == 2
        in_verbose = 1
    end

    #call cluster_split on every cluster created by rough clustering
    for ind in 1:length(indices)
        if verbose > 0
            println("Splitting cluster number $(ind)")
        end
        temp =cluster_split(kmer_vecs[indices[ind]], 1, nuc_inds, rough_radius=rough_radius, fine_radius=fine_radius,
            clust_multiplier=clust_multiplier, min_devs=min_devs, initial_devs=initial_devs, user_min_clust=user_min_clust,
            distfunc=euclidean,center=center, verbose=in_verbose, cycle_lim=cycle_lim, triangle=triangle, p_value=p_value)
        for arr in temp
            push!(fine_indices, indices[ind][arr])
        end
    end
    if verbose > 0
        println("Returning $(length(fine_indices)) templates")
    end
    return fine_indices
end

"""
     gen_arrangements(k)

Returns an array containing all possible kmer sequences of length k

"""
function gen_arrangements(k)
    alph = [["A"],["C"],["G"],["T"]]
    current = copy(alph)
    for i in 1:k-1
        current = vcat([[vcat(j,letter) for letter in alph] for j in current]...)
    end
    return [join(i) for i in current]
end
"""
	cluster_split(kmer_vecs, count, nuc_inds; rough_radius::Float64=0.01, fine_radius::Float64=1.0, clust_multiplier=1.0,
                        min_devs=6, initial_devs=20, user_min_clust=10, distfunc=euclidean,
                        center=mean, verbose::Int64=1, cycle_lim=20, triangle::Bool=false, p_value::Float64=0.01)

Determines if passed cluster is a 'real' cluster and recursively splits non-error clusters
"""
function cluster_split(kmer_vecs, count, nuc_inds; rough_radius::Float64=0.01, fine_radius::Float64=1.0, clust_multiplier=1.0,
                        min_devs=6, initial_devs=20, user_min_clust=10, distfunc=euclidean,
                        center=mean, verbose::Int64=1, cycle_lim=20, triangle::Bool=false, p_value::Float64=0.01)
    if verbose == 2
        println("______________________________________________")
        println("Current level of cluster_split recursion is $(count)")
        count+=1
    end

    k = Int(log(4, length(kmer_vecs[1])))
    num_vecs = length(kmer_vecs)

    #arr of each kmer's stdev throughout the kmer vectors
    stdevs = [std([kmer_vecs[seq_ind][kmer_ind] for seq_ind in 1:length(kmer_vecs)]) for kmer_ind in 1:4^k]
    #arr of indices of high->low stdev kmer indices
    dev_inds = reverse(sortperm(stdevs))[1:initial_devs]#setting the number of kmers high results in a lack of splitting

    #eliminate kmer indices that represent homopolymer differences
    dev_inds = get_non_homopolymers(dev_inds, nuc_inds, min_devs=min_devs)

    if length(dev_inds) == 0
        return [[i for i in 1:length(kmer_vecs)]]
    end
    #kmer counts of high stdev kmer indices
    top_kmers =[kmer[dev_inds] for kmer in kmer_vecs]
    #cluster kmers by high stdev indices
    us,sizes,indices,C_inds =cluster(top_kmers, fine_radius,distfunc=distfunc,center=center,
                                            verbose=verbose,cycle_lim=cycle_lim, triangle=triangle)
    #replace 2nd largest with each size to determine whether or not to split on it
    #p_val_arr = [(1-cdf(Poisson(maximum(sizes)*(error_rate/4)), sizes[i])) *sum(kmer_vecs[1]) for i in 1:length(sizes)]

    #minimum non-error cluster size. Estimate of the sum of poisson distributions
    #min_clust = clust_multiplier*.25*2.5*sqrt(rough_radius*sum(kmer_vecs[1])*2*log(maximum(sizes)))
    #minimum cluster size is the largest number between user input and the minimum non-error size
    keep_inds = find(user_min_clust.<sizes)
    #keep_inds = Array{Int64, 1}()
    #for ind in first_keep_inds

    #keep_inds = [ind for ind in first_keep_inds if (1-cdf(Poisson(maximum(sizes)*(rough_radius/4)), sizes[ind])) *sum(kmer_vecs[1]) < p_value]

    #if there are no clusters to keep, add everything to the largest one
    if length(keep_inds) == 0 || length(keep_inds) == 1
        return [[i for i in 1:length(kmer_vecs)]]
    else
        final_indices = Array{Array{Int64, 1}, 1}()
        keep_centers = [center(kmer_vecs[indices[keep]]) for keep in keep_inds]

        for keep in keep_inds
            push!(final_indices, indices[keep])
        end

        for i in 1:length(indices) #go through every cluster index
            if !(i in keep_inds) #if that cluster is not one of the kept clusters
                C_arr = indices[i] #get the array of seq indices corresponding to the cluster
                for index in C_arr #go through every seq indix in the non-kept cluster
                    push!(final_indices[indmin([distfunc(kmer_vecs[index], center) for center in keep_centers])],
                            index)
                end
            end
        end
    end

    #if the input has gotten this far, it must be split. Call cluster_split on each current cluster
    for clust_ind in 1:length(final_indices)
        outer_clust = copy(final_indices[clust_ind])
        new_clusts = cluster_split(kmer_vecs[outer_clust], count, nuc_inds, rough_radius=rough_radius, fine_radius=fine_radius,
                                    min_devs=min_devs, distfunc=distfunc,center=center,verbose=verbose,
                                    cycle_lim=cycle_lim, triangle=triangle, p_value=p_value)
        for new_clust in new_clusts
            push!(final_indices, outer_clust[new_clust])
        end
        final_indices = deleteat!(final_indices, clust_ind)

    end
    return final_indices
end

"""
    get_non_homopolymers(dev_inds::Array{Int64, 1}, nuc_inds; min_devs = 6)

Takes in array of kmer indices. Returns the highest standard deviation non-homopolymer indices

"""
function get_non_homopolymers(dev_inds::Array{Int64, 1}, nuc_inds; min_devs=6)
    k=log(4, length(nuc_inds))
    len = length(dev_inds)
    dev_seqs = nuc_inds[dev_inds]
    to_remove = Array{Int64, 1}()
    #go through every high deviation kmer index
    for outer_ind in 1:len
        #make sure to git every combination only once
        for inner_ind in outer_ind + 1:len
            #align corresponding sequences and determine if they represent homopolymer differences
            if are_homopolymers(nw_align(dev_seqs[outer_ind], dev_seqs[inner_ind], mismatch_cost=-100.0), k=k)
                push!(to_remove, outer_ind, inner_ind)
            end
        end
    end
    #difference between all indices and the homopolymer indices is the non homopolymer indices
    new_dev_inds = setdiff(dev_inds, dev_inds[unique(to_remove)])

    if length(new_dev_inds) < min_devs #if there aren't enough indices, empty array says "don't split!"
        return []
    else
        return new_dev_inds[1:min_devs] #pick the top non-homopolymer kmer indices
    end
end

"""
    are_homopolymers(alignment::Tuple{String,String}; k = 6, polylen = 3)

Determines if the aligned kmers represent a homopolymer difference. If the
kmers represent a homopolymer difference, the function returns true.

"""
function are_homopolymers(alignment::Tuple{String,String};k=6,polylen::Int64=3)

    al1 = alignment[1]
    al2 = alignment[2]
    len = length(alignment[1])
    if len != k + 1 #homopolymers only differ by one
        return false
    end
    if polylen > k - 1
        return false
    end

    ind1 = findfirst(alignment[1], '-')#findfirst(isequal('-'), alignment[1])
    ind2 = findfirst(alignment[2], '-')#findfirst(isequal('-'), alignment[2])
    diff = ind1-ind2
    #determine if inds are only on edges
    #        XXXXX- (ind1)
    # (ind2) -XXXXX
    if abs(diff) == k #both gaps are on edges
        str = al1[1:(1+polylen)] * al2[1:(1+polylen)]
        if (length(unique(al1[1:(1+polylen)] * al2[1:(1+polylen)])) == 2) || (length(unique(al1[(len-polylen):len] * al2[(len-polylen):len])) == 2)
            return true
        else
            return false
        end
    elseif ind1 == 1 || ind2 == 1 #one gap on left edge, one on non-edge
        if ind1 == 1
            temp = ind2
        else
            temp = ind1
        end
        if length(unique(al1[1:(1+polylen)] * al2[1:(1+polylen)])) == 2 #edge box is true. THIS MEANS NO HOMOPOLY!!!
            return false

        else #edge box is false. May be a homopoly.
            if temp + polylen <= len && length(unique(al1[temp:(temp+polylen)]*al2[temp:(temp+polylen)])) == 2
                return true
            elseif temp - polylen >= 1 && length(unique(al1[(temp-polylen):temp]*al2[(temp-polylen):temp])) == 2
                return true
            end
        end
        if ind1 == 1
            al = al2
        else
            al=al1
        end
        #if you're here, it means that the non-edge gap is inside, not bordering, the homopoly diff
        if polylen + 2 < len
            if abs(diff) == 1 && length(unique(al[1:polylen+2])) == 2
                return true
            elseif abs(diff) == 2 && length(unique(al[2:polylen+3])) == 2 && length(unique(al1[[ind1, ind2]]*al2[[ind1, ind2]])) == 2
                return true
            end
        end
    else #one of the indices is k+1
        if ind1 == len
            temp = ind2
        else
            temp = ind1
        end
        if length(unique(al1[(len-polylen):len] * al2[(len-polylen):len])) == 2 #edge box is true. THIS MEANS NO HOMOPOLY!!!
            return false
        else #edge box is false. May be a homopoly

            if temp + polylen <= len && length(unique(al1[temp:(temp+polylen)]*al2[temp:(temp+polylen)])) == 2

                return true
            elseif temp - polylen >= 1 && length(unique(al1[(temp-polylen):temp]*al2[(temp-polylen):temp])) == 2
                return true
            end
        end
        if ind1 == len
            al=al2
        else
            al=al1
        end
        if polylen + 2 < len
            if abs(diff) == 1 && length(unique(al[len-polylen-1:len])) == 2
                return true
            elseif abs(diff) == 2 && length(unique(al[len-polylen-2:len-1])) == 2 && length(unique(al1[[ind1, ind2]]*al2[[ind1, ind2]])) == 2
                return true
            end
        end
    end
    #if the alignment did not meet any of these conditions, it must be a nonhomopolymer diff
    return false
    #=
    arr = [alignment[1]... , alignment[2]... ]
    #max homopolymer length difference is 2
    if length(findin(arr, '-')) > 2
        return false
    end

    len = length(alignment[1])

    #check left->right
    for i in 1:len
        if i+polylen > len #no chance of sufficiently large polymer, break
            break;
        end
        if (alignment[1][i] == '-') #test if this is start of homopolymer
            @inbounds temp = alignment[2][(i+1):(i+polylen)]
            temp2 = alignment[1][i:i+polylen]

            if length(unique(temp)) == 1 && length(unique(temp2))== 2
                return true
            end
        elseif (alignment[2][i] == '-')
            @inbounds temp = alignment[1][(i+1):(i+polylen)]
            temp2 = alignment[2][i:i+polylen]
            if length(unique(temp)) == 1 && length(unique(temp2))== 2
                return true
            end
        end
    end
    #check right->left
    for i in len:-1:1
        if i - polylen < 1 #no chance of sufficiently large polymer, return
            return false
        end
        if (alignment[1][i] == '-')
            @inbounds temp = alignment[2][(i-polylen):(i-1)]
            temp2 = alignment[1][i-polylen:i]
            if length(unique(temp)) == 1 && length(unique(temp2)) == 2
                return true
            end
        elseif (alignment[2][i] == '-')
            @inbounds temp = alignment[1][(i-polylen):(i-1)]
            temp2 = alignment[2][i-polylen:i]
            if length(unique(temp)) == 1 && length(unique(temp2)) == 2
                return true
            end
        end
    end
    return false
    =#
end

#==
      OLD VERSION
==#

"""
    alignment_consensus(seqs)

Computes a consensus sequence for array of given sequences that have already been aligned.
Consensus is computed by site-wise mode among given sequences.
"""
function alignment_consensus(seqs)
    string_lengths = [length(i) for i in seqs]
    ncols = length(seqs[1])
    if !all([x == ncols for x in string_lengths])
        error("Aligned sequences must all be the same length.")
    end
    result = string([mode([seq[i] for seq in seqs]) for i in 1:ncols]...)
    return degap(result)
end

"""
    coords(ref, read)

Returns an array of degapped coordinates, such that
coords(ref, read)[i] gives you the position the aligned read/ref
that matches the i'th ungapped position in ref.
"""
function coords(ref, read)
    if length(ref) != length(read)
        error("Aligned strings are meant to be the same length.")
    end
    degappedRef = degap(ref)
    coordMap = Array{Int64}(length(degappedRef))
    count = 1
    for i in 1:length(degappedRef)
        while ref[count] == '-'
            count += 1
        end
        coordMap[i] = count
        count += 1
    end
    return coordMap
end

"""
    get_centroid(reads, k, distfunc = corrected_kmer_dist)

Computes a rough cluster centroid by returning the sequence in `reads`
that has the nearest distance in `k`-mer space to the mean `k`-mer vector.
"""
function get_centroid(reads, k, distfunc = corrected_kmer_dist)
    kmerVecs = [kmer_count(read, k) for read in reads]
    meanVec = mean(kmerVecs)
    dists = [distfunc(kmerVec, meanVec, k=k) for kmerVec in kmerVecs]
    return reads[indmin(dists)]
end

"""
    consensus_seq(reads; thresh = 0.7, shift = 1, k = 6,
                       distfunc = corrected_kmer_dist)

Computes a cluster centroid by taking the nearest read to the mean `k`-mer
vector (see `get_centroid(reads, k, distfunc)`) and refining it (see `refine_ref(ref, reads)`).
`shift` determines the window size of comparison between sequences when refining rough centroid
locally (actual window size is `shift+1`).
"""
function consensus_seq(reads; degap_param = true, thresh = 0.7, shift = 1, k = 6,
                       distfunc = corrected_kmer_dist)
    ref = get_centroid(reads, k, distfunc)
    return refine_ref(ref, reads, degap_param = degap_param, thresh = thresh, shift = shift)
end

#=function get_matches(candidate_ref, reads, shift; kmer_align = true)
    alignments = []
    if kmer_align
        alignments = [kmer_seeded_align(candidate_ref, i) for i in reads]
    else
        alignments = [nw_align(candidate_ref, i) for i in reads]
    end
    maps = [coords(i...) for i in alignments]
    matchContent = [[degap(alignments[i][2][maps[i][k]:maps[i][k+shift]]) for i in 1:length(maps)] for k in 1:length(candidate_ref)-shift]
    matches = [freq(matchContent[k], degap(candidate_ref[k:k+shift])) for k in 1:length(matchContent)]
    return alignments, maps, matches, matchContent
end=#
#warn("Multithreaded version active as sole get_matches")

"""
    get_matches(candidate_ref, reads, shift; kmer_align = true)

Determine alignments between each sequence of `reads` and `candidate_ref`, along with coordinate mappings
and local matches. Used internally in `refine_ref` and `peak_split`.
"""
function get_matches(candidate_ref, reads, shift; degap_param = true, kmer_align = true)
    alignments = []
    if kmer_align
        alignments = map(i -> kmer_seeded_align(candidate_ref, i), reads)
    else
        alignments = map(i -> nw_align(candidate_ref, i), reads)
    end
    maps = [coords(i...) for i in alignments]

    if (degap_param)
        matchContent = [[degap(alignments[i][2][maps[i][k]:maps[i][k+shift]]) for i in 1:length(maps)] for k in 1:length(candidate_ref)-shift]

        matches = [freq(matchContent[k], degap(candidate_ref[k:k+shift])) for k in 1:length(matchContent)]

    else
       matchContent = [[(alignments[i][2][maps[i][k]:maps[i][k+shift]]) for i in 1:length(maps)] for k in 1:length(candidate_ref)-shift]

       matches = [freq(matchContent[k], candidate_ref[k:k+shift]) for k in 1:length(matchContent)]

    end

    return alignments, maps, matches, matchContent
end

"""
    refine_ref(candidate_ref, reads;  thresh = 0.7, shift = 1)

Takes a candidate consensus sequence, `candidate_ref`, and corresponding array of `reads` and refines the candidate via majority
votes from reads at local windows of size `shift+1`. If after alignment the frequency of a local region of the candidate
is less than `thresh`, this part of the candidate is refined.

"""
function refine_ref(candidate_ref, reads; degap_param = true, thresh = 0.7, shift = 1)
    # need to add in code to handle deletions at the start and end of
    # the candidate!

    alignments, maps, matches, _ = get_matches(candidate_ref, reads, shift, degap_param = degap_param)

    boolVec = BitArray{1}(1+zeros(length(candidate_ref)))
    for i in 1:length(matches)
        if matches[i] < thresh
            boolVec[i:i+shift] = false
        end
    end

    # need to explicitly handle the case where the reference is shorter than the reads.
    # refine front
    frontStrCol = [degap((alignments[j][2])[1:maps[j][1]-1]) for j in 1:length(maps)]
    frontNonEmpties = [length(s) != 0 for s in frontStrCol]
    frontConsensus = ""
    majority = length(alignments) / 2
    while countnz(frontNonEmpties) > majority
        frontStrCol = frontStrCol[frontNonEmpties]
        frontConsensus = mode(frontStrCol)
        frontAli, frontMaps, _, _ = get_matches(frontConsensus, frontStrCol, 1, degap_param = degap_param, kmer_align=false)
        frontStrCol = [degap(frontAli[j][2][1:frontMaps[j][1]-1]) for j in 1:length(frontMaps)]
        frontNonEmpties = [length(s) != 0 for s in frontStrCol]
        majority = length(frontAli) / 2
    end
    frontStr = frontConsensus

    # refine end
    endStrCol = [degap((alignments[j][2])[maps[j][end]+1:end]) for j in 1:length(maps)]
    endNonEmpties = [length(s) != 0 for s in endStrCol]
    endConsensus = ""
    majority = length(alignments) / 2
    while countnz(endNonEmpties) > majority
        endStrCol = endStrCol[endNonEmpties]
        endConsensus = mode(endStrCol)
        endAli, endMaps, _, _ = get_matches(endConsensus, endStrCol, 1, degap_param = degap_param, kmer_align=false)
        endStrCol = [degap(endAli[j][2][endMaps[j][end]+1:end]) for j in 1:length(endMaps)]
        endNonEmpties = [length(s) != 0 for s in endStrCol]
        majority = length(endAli) / 2
    end
    endStr = endConsensus

    # # if something breaks, CHECK THIS FIRST!
    if sum(1-boolVec) < 0.5
        return frontStr * candidate_ref * endStr
    end
    ranges = []
    i=1
    while i < length(boolVec)
        if !boolVec[i]
            start = i
            i += 1
            while !boolVec[i] && (i < length(boolVec))
                i += 1
            end
            push!(ranges,[start, i-1])
        end
        i += 1
    end

    newRef = candidate_ref[1:ranges[1][1]-1]
    for i in 1:length(ranges)-1
        low, high = ranges[i][1], ranges[i][2]
        insStr = mode([degap((alignments[j][2])[maps[j][low]:maps[j][high]])
                       for j in 1:length(maps)])
        newRef = newRef*insStr
        newRef = newRef*candidate_ref[ranges[i][2]+1:ranges[i+1][1]-1]
    end

    low, high = ranges[length(ranges)][1], ranges[length(ranges)][2]
    insStr = mode([degap((alignments[j][2])[maps[j][low]:maps[j][high]]) for j in 1:length(maps)])
    newRef = newRef*insStr
    newRef = newRef*candidate_ref[ranges[length(ranges)][2]+1:length(candidate_ref)]
    return frontStr*newRef*endStr
end

"""
    recursive_split(cluster, glom; centroid = "", shift = 7, minreads = 5, polylen = 3, k = 6, reads = nothing, phreds = nothing, names = nothing)

Recursively splits input `cluster` using `peak_split` and populates `glom` with the separated sub-clusters.
Splits clusters according to `centroid` (unless not given, in which case it is inferred).
`polylen` is lower bound for run-on subsequence to by considered a homopolymer region and not be split on.
Does not split on clusters of size less than `minReads`.
If `reads` is given as a string array, this will populate `reads` with tuples of `(cluster, phreds, names)` for
determining the actual clusters that are split, rather than just corresponding consensus sequences, where `phreds`
and `names` may be arrays of corresponding values, or `nothing`.
"""
function recursive_split(cluster::Array{String}, glom; centroid = "", shift = 7, minreads = 5, polylen = 3, k = 6, reads = nothing, phreds = nothing, names = nothing)
    consensus = ""
    if length(centroid) == 0
        consensus = consensus_seq(cluster, k=k)
    else
        consensus = refine_ref(centroid, cluster)
    end
    splitClus = peak_split(consensus, cluster, shift=shift, minReads=minreads, polylen=polylen)

    if length(splitClus) > 1
        for clust in splitClus
            split_phreds = phreds==nothing ? nothing : phreds[clust]
            split_names = names==nothing  ? nothing : names[clust]
            recursive_split(cluster[clust], glom, shift=shift, minreads=minreads, polylen=polylen, k=k, reads=reads, phreds=split_phreds, names=split_names)
        end
    else
        push!(glom, [consensus, length(cluster)])
        if reads != nothing
            push!(reads, (cluster, phreds, names))
        end
    end
end

"""
    peak_split(candidate_ref, reads; thresh = 0.7, shift = 3, minReads = 5, polylen = 3)

Creates and returns sub-clusters of `reads` split by variation from `candidateRef` if, locally, the most common read is less frequent than given `thresh`, default 0.7.
If the size of the input cluster `reads` is large, then `thresh` may be changed: thresh = max(1-10/length(reads), thresh).
Size of local window/site for comparison = `shift+1`.
Return values are arrays of boolean arrays, where each boolean array can be used to index into the input reads data to give the split clusters.
Determines  treshold for a cluster split, where if the candidate has local region with frequency less than `thresh` the cluster is split,
 to be `thresh = max(1-10/length(reads), thresh)` to account for larger cluster sizes.
`polylen` is lower bound for run-on subsequence to by considered a homopolymer region and not be split on.
Does not split on clusters of size less than `minReads`.
"""
function peak_split(candidate_ref, reads; thresh = 0.7, shift = 3, minReads = 5, polylen = 3)
    candidateRef = candidate_ref
    threshProp = thresh
    if length(reads) < minReads
        return [BitArray([true for r in reads])]
    end
    # Adjusting thresh is meant to account for the fact that 30% of a very large cluster is still very large,
    #  so we assume that, for large clusters, if we have 10 seqs, that is enough to split into a new cluster.
    # However, this will fail for very very large clusters that have 10 sequences that have the same mismatch base error.
    # Need to build in the expected number of mismatched base errors at any given base to control the false positive rate here.
    thresh = max(1-10/length(reads), thresh)
    numSteps = length(candidateRef) - shift
    #Need to add in code to handle deletions at the start and end of the candidate!
    alignments, maps, matches, matchContent = get_matches(candidateRef, reads, shift)
    boolVec = BitArray{1}(1 + zeros(length(candidateRef)))
    for i in 1:length(matches)
        if (matches[i] < thresh)
            boolVec[i:i+shift] = false
        end
    end
    # TODO: make ranges statically sized to avoid push!()
    ranges = []
    i = 1
    while i <= numSteps
        if !boolVec[i]
            start = i
            while !boolVec[i] && i < numSteps
                i += 1
            end
            # Julia subranges are inclusive, => window size = (shift+1)
            if shift+1 < polylen
                start = max(1, start - (polylen - (shift+1)))
                i = min(numSteps, i + (polylen - (shift+1)))
            end
            push!(ranges, [start, i-1])
        end
        i += 1
    end

    # add 2 to lengths to account for front and end as well
    scoreList = 1 + zeros(length(ranges) + 2)
    scoreInds = zeros(UInt32, length(ranges) + 2)
    for i in 1:length(ranges)
        low, high = ranges[i][1], ranges[i][2]
        sortedFreqs = sorted_freqs([degap(alignments[j][2][maps[j][low]:maps[j][high]]) for j in 1:length(maps)])
        # The following might happen with discrepancies at ends of reads; this is taken care of later in this function
        if length(sortedFreqs) == 1
            continue
        end
        # This will not handle the case where the third most frequent variant is a true mutation and the second is noise, but this case only occurs or is important when:
        #  1) the second most frequent, noisy read is more common than reads of a mutation (ie in a very long homopolymer)
        #  2) the third variant, the true mutation, doesn't cause a split ANYWHERE else, so that the mutation is basically a single base change at a homopolymer region (as opposed to a noisy indel)
        # So it is very uncommon for this case to happen but it might need to be handled.
        ali = nw_align(sortedFreqs[1][2], sortedFreqs[2][2])
        if !differ_by_just_one_gap(ali, polylen=polylen)
            scoreList[i] = minimum(matches[ranges[i][1]:ranges[i][2]])
            scoreInds[i] = indmin(matches[ranges[i][1]:ranges[i][2]]) + ranges[i][1] - 1
        end
    end

    # take care of splits at overhanging beginnings of reads, where consensus is empty ("")
    # if any overhang is more common than threshold, potentially split on most common overhang
    frontStrCol = [degap(alignments[j][2][1:maps[j][1]-1]) for j in 1:length(maps)]
    frontOverhangs = [length(s) != 0 for s in frontStrCol]
    if freq(frontStrCol, "") < thresh
        sortedFreqs = sorted_freqs([degap(alignments[j][2][1:maps[j][polylen]])
                    for j in 1:length(maps) if frontOverhangs[j]])
        if length(sortedFreqs) == 1
            scoreList[end-1] = freq(frontStrCol, mode(frontStrCol[frontOverhangs]))
            scoreInds[end-1] = 1
        else
            ali = nw_align(sortedFreqs[1][2][1:min(polylen, end)], sortedFreqs[2][2][1:min(polylen, end)])
            if !differ_by_just_one_gap(ali, polylen=polylen)
                scoreList[end-1] = freq(frontStrCol, mode(frontStrCol[frontOverhangs]))
                scoreInds[end-1] = 1
            end
        end
    end
    # take care of splits at overhanging ends
    endStrCol = [degap(alignments[j][2][maps[j][end]+1:end]) for j in 1:length(maps)]
    endOverhangs = [length(s) != 0 for s in endStrCol]
    if freq(endStrCol, "") < thresh
        sortedFreqs = sorted_freqs([degap(alignments[j][2][maps[j][end-polylen]:end])
                    for j in 1:length(maps) if endOverhangs[j]])
        if length(sortedFreqs) == 1
            scoreList[end] = freq(endStrCol, mode(endStrCol[endOverhangs]))
            scoreInds[end] = 1
        else
            ali = nw_align(sortedFreqs[1][2][max(end-polylen, 0):end], sortedFreqs[2][2][max(end-polylen, 0):end])
            if !differ_by_just_one_gap(ali, polylen=polylen)
                scoreList[end] = freq(endStrCol, mode(endStrCol[endOverhangs]))
                scoreInds[end] = 1
            end
        end
    end
    # TODO: maybe add singleton clusters back to cluster with closest distance
    if minimum(scoreList) < thresh
        ind = indmin(scoreList)
        splitpos = scoreInds[ind]
        # if in middle of sequence
        if ind < length(scoreInds) - 1
            unio = union(matchContent[splitpos])
            splitCluster = [matchContent[splitpos] .== s for s in unio]
            return splitCluster
        # else if on front
        elseif ind == length(scoreInds) - 1
            unio = union(frontStrCol)
            splitCluster = [frontStrCol .== s for s in unio]
            return splitCluster
        # else if on end
        else
            unio = union(endStrCol)
            splitCluster = [endStrCol .== s for s in unio]
            return splitCluster
        end
    else
        return [BitArray([true for r in reads])]
    end
end

"""
    differ_by_just_one_gap(alignment; polylen = 3)

`alignment` is an array containing two already aligned strings.
Returns true if first and second sequences in given array differ by
just one gap where that gap is in a homopolymer region.
A homopolymer region is determined to be a region of a single repeated
nucleotide of length at least `polylen`.
"""
function differ_by_just_one_gap(alignment; polylen = 3)
    gapDiffCount = 0
    otherDiffCount = 0
    gapInd = 0
    noGapSeq = 0
    for i in 1:length(alignment[1])
        if (alignment[1][i] == '-')
            gapDiffCount += 1
            gapInd = i
            noGapSeq = 2
        elseif (alignment[2][i] == '-')
            gapDiffCount += 1
            gapInd = i
            noGapSeq = 1
        elseif (alignment[1][i] != alignment[2][i])
            otherDiffCount += 1
        end
    end
    if otherDiffCount == 0 && gapDiffCount == 1
        # check if homopolymer region by matching regex to local string
        lo, hi = max(gapInd - (polylen-1), 1), min(gapInd + (polylen-1), length(alignment[1]))
        regex = Regex(string(alignment[noGapSeq][gapInd]) * "{" * string(polylen) * "}")
        return ismatch(regex, alignment[noGapSeq][lo:hi])
    else
        return false
    end
end

"""
    consensus_viz(candidate_ref, reads; thresh = 0.7, shift = 3,
                  intitle = "Consensus agreement.")

Creates a plot to visualize the agreement of a consensus sequence, `candidate_ref`,
with its cluster of `reads`.
Size of local window/site for comparison = `shift+1`.
I'm pretty sure `thresh` does nothing here.
"""
function consensus_viz(candidate_ref, reads;
                       thresh = 0.7, shift = 3, intitle = "Consensus agreement.")
    #Need to add in code to handle deletions at the start and end of the candidate!
    alignments, maps, matches, _ = get_matches(candidate_ref, reads, shift)
    frontStrCol = freq([degap((alignments[j][2])[1:maps[j][1]-1]) for j in 1:length(maps)], "")
    endStrCol = freq([degap((alignments[j][2])[maps[j][end]+1:end]) for j in 1:length(maps)], "")
    w, h = plt[:figaspect](0.1)
    figure(figsize=(w, h))
    plot(1:length(matches)+2, vcat([frontStrCol], matches,[endStrCol]))
    ax = gca()
    ax[:set_xlim]([-10, length(matches)+12])
    ax[:set_ylim]([0, 1])
    title(intitle)
    xlabel("Consensus Position")
    ylabel("Agreement")
end

"""
    disagreements(candidate_ref, reads; thresh = 0.7, shift = 3)

Prints local disagreements between `candidate_ref` and each sequence of `reads` after aligning.
Size of local window/site for comparison = `shift+1`.
A disagreement is a region where the candidate has a local region with frequency less than `thresh`.
"""
function disagreements(candidate_ref, reads; thresh = 0.7, shift = 3)
    # need to add in code to handle deletions at the start and end of
    # the candidate!
    alignments = [kmer_seeded_align(candidate_ref, i) for i in reads]
    maps = [coords(i...) for i in alignments]

    matches = [freq([(alignments[i][2])[maps[i][k]:maps[i][k+shift]]
                     for i in 1:length(maps)], candidate_ref[k:k+shift])
               for k in 1:length(candidate_ref)-shift]

    boolVec = BitArray{1}(1+zeros(length(candidate_ref)))
    for i in 1:length(matches)
        if matches[i] < thresh
            boolVec[i:i+shift] = boolVec[i:i+shift] & [false for k in 1:1+shift]
        end
    end

    ranges = []
    i = 1
    while i < length(boolVec)
        if !boolVec[i]
            start = i
            i += 1
            while !boolVec[i] && i < length(boolVec)
                i += 1
            end
            push!(ranges,[start, i-1])
        end
        i += 1
    end

    frontStrCol = [degap((alignments[j][2])[1:maps[j][1]-1]) for j in 1:length(maps)]
    if freq(frontStrCol, "") < thresh
        println("Before ref starts:")
        freq_dict_print(countmap(frontStrCol), thresh=1.1)
    end

    for i in 1:length(ranges)
        low, high = ranges[i][1], ranges[i][2]
        insStrCol = [degap((alignments[j][2])[maps[j][low]:maps[j][high]]) for j in 1:length(maps)]
        if freq(insStrCol, candidate_ref[low:high]) < thresh
            println("Ref pos ", low, ":", candidate_ref[low:high])
            freq_dict_print(countmap(insStrCol), thresh=1.1)
        end
    end

    endStrCol = [degap((alignments[j][2])[maps[j][end]+1:end]) for j in 1:length(maps)]
    if freq(endStrCol, "") < thresh
        println("After ref ends:")
        freq_dict_print(countmap(endStrCol), thresh=1.1)
    end
end

"""
    diff_in_homopolymer_region(alignment::Array{String, 1}; polylen=3)

Returns true if two aligned sequences differ only by single gaps in homopolyer regions
(ie one gap per region). `alignment` is an array of two strings that have already been aligned.
A homopolymer region is determined to be a region of a single repeated
nucleotide of length at least `polylen`.
"""
function diff_in_homopolymer_region(alignment::Array{String, 1}; polylen=3)
    for i in 1:length(alignment[1])
        if alignment[1][i] == '-' || alignment[2][i] == '-'
            lo = max(1, i - polylen + 1)
            hi = min(length(alignment[1]), i + polylen - 1)
            if !differ_by_just_one_gap([alignment[1][lo:hi], alignment[2][lo:hi]], polylen=polylen)
                return false
            end
        end
    end
    return true
end

"""
    merge_consensuses(seqs)

Takes in an array of 2-array [sequence, size] pairs and returns a list with
reach sequence unique with summed sizes among all previous copies.
"""
function merge_consensuses(seqs)
    tmp = sort(seqs, by=s->s[1]);
    i = 1
    while i < length(tmp)
        if tmp[i][1] == tmp[i+1][1]
            # consolidate cluster sizes
            tmp[i][2] += tmp[i+1][2]
            deleteat!(tmp, i+1)
        else
            i += 1
        end
    end
    return tmp
end

"""
    merge_diff_homopolymers(seqs; thresh = 0.002, k = 6, polylen = 3)

Takes in an array of 2-array [sequence, size] pairs and combines sequences that are identical
or differ by single gaps in homopolymer regions, summing respective frequencies. Slower than
naive `merge_consensuses(seqs)`.
A homopolymer region is determined to be a region of a single repeated
nucleotide of length at least `polylen`.
`thresh` gives an upper bound on distance in `k`-mer space for sequences to be compared for
homopolymer differences.
"""
function merge_diff_homopolymers(seqs; thresh = 0.002, k = 6, polylen = 3)
    # sort by decreasing size
    tmp = reverse(sort(seqs, by=s->s[2]))
    kmer_vecs = [kmer_count(s[1], k) for s in tmp]
    i = 1
    while i < length(tmp)
        j = i + 1
        while j < length(tmp)
            # if identical
            if tmp[i][1] == tmp[j][1]
                # consolidate cluster sizes
                tmp[i][2] += tmp[j][2]
                deleteat!(tmp, j)
                deleteat!(kmer_vecs, j)
            # if differ in homopolyer region
            elseif corrected_kmer_dist(kmer_vecs[i], kmer_vecs[j], k=k) < thresh &&
                    diff_in_homopolymer_region(kmer_seeded_align(tmp[i][1], tmp[j][1]), polylen=polylen)
                # consolidate cluster sizes
                tmp[i][2] += tmp[j][2]
                deleteat!(tmp, j)
                deleteat!(kmer_vecs, j)
            else
                j += 1
            end
        end
        i += 1
    end
    return tmp
end

"""
    get_coarse_centroid(seqs::Array{String, 1}; subsample = 1000, k = 4)

Returns a 'master consensus' representing largest coarse cluster of sequences,
computed among a subsample of `reads`. The consensus is the closest sequence of
`reads` to the mean `k`-mer vector in the largest computed cluster.
"""
function get_coarse_centroid(seqs::Array{String, 1}; subsample = 1000, k = 4)
    if subsample >= length(seqs)
        sampled_seqs = seqs
    else
        inds = shuffle([1:length(seqs)...])[1:subsample]
        sampled_seqs = seqs[inds]
    end
    kmer_vecs = [kmer_count(seq, k) for seq in sampled_seqs]
    @time _, sizes, _, centroid_inds = dp_centers(kmer_vecs, 0.15)
    return sampled_seqs[centroid_inds[indmax(sizes)]]
end


text:
module DashJbrowse
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "1.0.4"

include("jl/circulargenomeview.jl")
include("jl/lineargenomeview.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_jbrowse",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "async-LinearGenomeView.js",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/async-LinearGenomeView.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-CircularGenomeView.js",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/async-CircularGenomeView.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-LinearGenomeView.js.map",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/async-LinearGenomeView.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-CircularGenomeView.js.map",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/async-CircularGenomeView.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_jbrowse.min.js",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/dash_jbrowse.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_jbrowse.min.js.map",
    external_url = "https://unpkg.com/dash_jbrowse@1.0.4/dash_jbrowse/dash_jbrowse.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end


text:
module StackedHourglass

using CUDA, Knet
const KnetMoment = Knet.Ops20.BNMoments

using Knet.CuArrays

#standard Library
using Distributed, Random

using Images, MAT, FFTW

#exported types
export HG2

#exported methods
export subpixel, set_testing, save_hourglass, load_hourglass, features

abstract type NN end;
const HGType = Union{KnetArray{Float32,4},AutoGrad.Result{KnetArray{Float32,4}},Array{Float32,4},AutoGrad.Result{Array{Float32,4}}}

const PType1 = Union{Param{KnetArray{Float32,1}},Param{Array{Float32,1}}}
const PType4 = Union{Param{KnetArray{Float32,4}},Param{Array{Float32,4}}}

include("residual.jl")
include("hourglass.jl")
include("cuda_resizing.jl")

include("gaussian_pyramids.jl")
include("cuda_files.jl")
include("helper.jl")
include("load.jl")
include("subpixel.jl")
include("image_preprocessing.jl")

end


text:
# Normal distribution

import StatsFuns
export Normal

import Base: eltype

# Note: `@measure Normal Lebesgue` expands to (cleaned up)

# quote
#     struct Normal{P, X} <: AbstractMeasure{X}
#         par::P
#     end

#     function Normal(nt::NamedTuple)
#         P = typeof(nt)
#         return Normal{P, eltype(Normal{P})}
#     end

#     Normal(; kwargs...) = Normal((; kwargs...))

#     (baseMeasure(μ::Normal{P, X}) where {P, X}) = Lebesgue{X}

#     (:≪(::Normal{P, X}, ::Lebesgue{X}) where {P, X}) = true

# end
# @measure Normal Lebesgue

@measure Normal(μ,σ) ≃ (1/sqrt2π) * Lebesgue(X)

function logdensity(d::Normal{P} , x::X) where {P <: NamedTuple{(:μ, :σ)}, X}
    return - (log(2) + log(π)) / 2 - log(d.par.σ)  - (x - d.par.μ)^2 / (2 * d.par.σ^2)
end

# Standard normal

function logdensity(d::Normal{EmptyNamedTuple,X} , ::Lebesgue{X}) where {X <: Real}
    return - (log(2) + log(π)) / 2  - x^2 / 2
end

Normal() = Normal(NamedTuple())

# @implement HasDensity{Normal{P,X},X} where {X, P <: NamedTuple{(:μ, :σ)}} begin
#     baseMeasure(d) = Lebesgue(X)
#     logdensity(d, x) = - (log(2) + log(π)) / 2 - log(d.par.σ)  - (x - d.par.μ)^2 / (2 * d.par.σ^2)
# end

# # μ, τ

# eltype(::Type{Normal}, ::Type{NamedTuple{(:μ, :τ), Tuple{A, B}}}) where {A,B} = promote_type(A,B)

# @implement HasDensity{Normal{P,X},X} where {X, P <: NamedTuple{(:μ, :τ)}} begin
#     baseMeasure(d) = Lebesgue(X)
#     logdensity(d, x) = - (log(2) + log(π) - log(d.par.τ)  + d.par.τ * (x - d.par.μ)^2) / 2
# end

# # σ

# eltype(::Type{Normal}, ::Type{NamedTuple{(:σ,), Tuple{B}}}) where {B} = B

# @implement HasDensity{Normal{P,X},X} where {X, P <: NamedTuple{(:σ,)}} begin
#     baseMeasure(d) = Lebesgue(X)
#     logdensity(d, x) = - (log(2) + log(π)) / 2 - log(d.par.σ)  - x^2 / (2 * d.par.σ^2)
# end

# # τ

# eltype(::Type{Normal}, ::Type{NamedTuple{(:τ,), Tuple{B}}}) where {B} = B

# @implement HasDensity{Normal{P,X},X} where {X, P <: NamedTuple{(:τ,)}} begin
#     baseMeasure(d) = Lebesgue(X)
#     logdensity(d, x) = - (log(2) + log(π) - log(d.par.τ)  + d.par.τ * x^2) / 2
# end

# @implement HasRand{Normal{P,X},X} where {X, P <: NamedTuple{(:μ, :σ)}} begin
#     rand(d) = Distributions.rand(Normal(d.par.μ,d.par.σ))
# end

# @implement IsMeasurable{Normal{P,X},S,X} where {X <: Real, S <: Interval{X}, P <: NamedTuple{(:μ, :σ)}} begin
#     measure(d, i) = StatsFuns.normcdf(d.par.μ, d.par.σ, i.last) - StatsFuns.normcdf(d.par.μ, d.par.σ, i.first)
# end

# logdensity(Normal(σ=0.5), 4.0)
#  == logdensity(Normal(μ=0.0, σ=0.5), 4.0)
#  == logdensity(Normal(μ=0.0, τ=4.0), 4.0)
#  == logdensity(Normal(τ=4.0), 4.0)


text:
module Ariadne

using Luxor

export Grid
export Maze
export PolarGrid
export binarytree!
export colidx
export dijkstra
export dijkstra!
export distances
export distances!
export east
export grid
export hasdistances
export haspath
export islinked
export link!
export links
export neighbors
export north
export path
export path!
export pathto
export rowidx
export sidewinder!
export south
export totxt
export topng
export unlink!
export west

mutable struct Cell
    rowidx::Int
    colidx::Int
    links::Set{Cell}
end

abstract type AbstractDistances <: AbstractDict{Cell,Int} end
abstract type AbstractGrid <: AbstractMatrix{Cell} end
abstract type AbstractMaze <: AbstractGrid end

Cell(rowidx, colidx) = Cell(rowidx, colidx, Set{Cell}())

Base.show(io::IO, cell::Cell) = print(io, "($(rowidx(cell)), $(colidx(cell)))")

rowidx(cell) = cell.rowidx
colidx(cell) = cell.colidx
links(cell) = cell.links

islinked(::Nothing, _) = false
islinked(_, ::Nothing) = false
islinked(a, b) = b in links(a)

function link!(a, b, reverse = true)
    push!(links(a), b)
    reverse && link!(b, a, false)
    return a
end

function unlink!(a, b, reverse = true)
    pop!(links(a), b)
    reverse && unlink!(b, a, false)
    return a
end

cellmatrix(nrows, ncols) = [Cell(r, c) for r in 1:nrows, c in 1:ncols]

struct Distances <: AbstractDistances
    root::Cell
    cells::Dict{Cell,Int}
end

function Distances(root)
    cells = Dict{Cell,Int}()
    cells[root] = 0
    frontier = [root]
    while !isempty(frontier)
        newfrontier = empty([], Cell)
        for cell in frontier
            for link in links(cell)
                if link in keys(cells)
                    continue
                end
                cells[link] = cells[cell] + 1
                push!(newfrontier, link)
            end
        end
        frontier = newfrontier
    end
    return Distances(root, cells)
end

Base.length(d::Distances) = length(cells(d))

Base.iterate(d::Distances) = iterate(cells(d))
Base.iterate(d::Distances, s) = iterate(cells(d), s)

Base.get(d::Distances, key, default) = get(cells(d), key, default)

Base.minimum(d::Distances) = findmin(cells(d))
Base.maximum(d::Distances) = findmax(cells(d))

root(d::Distances) = d.root
cells(d::Distances) = d.cells

dijkstra(root) = root |> Distances

function pathto(d, goal)
    current = goal
    path = Distances(root(d), Dict{Cell,Int}())
    cells(path)[current] = d[current]
    while current != root(d)
        for link in links(current)
            if d[link] < d[current]
                cells(path)[link] = d[link]
                current = link
                break
            end
        end
    end
    return path
end

mutable struct Maze <: AbstractMaze
    grid::AbstractGrid
    distances::Union{AbstractDistances,Nothing}
    path::Union{AbstractDistances,Nothing}
end

Maze(g) = Maze(g, nothing, nothing)

grid(m::Maze) = m.grid
distances(m::Maze) = m.distances
path(m::Maze) = m.path

Base.size(m::Maze) = size(grid(m))
Base.getindex(m::Maze, r, c) = getindex(grid(m), r, c)

function distances!(m::Maze, d)
    m.distances = d
    return m
end

function path!(m::Maze, p)
    m.path = p
    return m
end

function dijkstra!(m::Maze)
    distances!(m, dijkstra(first(m)))
    return m
end

hasdistances(m::Maze) = !isnothing(distances(m))
haspath(m::Maze) = !isnothing(path(m))

include("gamma.jl")
include("theta.jl")
include("algorithms.jl")

end # module

max_stars_repo_name: basp/ariadne

max_issues_repo_name: basp/ariadne

max_forks_repo_name: basp/ariadne


text:
# function to perform tests of the rotation functions in the Rotations module

#=
########################
# Define helper methods
########################

import Rotations: numel
numel{T <: Rotations.RotationTypes}(X::T) = Rotations.numel(Rotations.strip_eltype(T))
numel{T}(X::T) = (T <: FixedSizeArrays.FixedArray) || (T <: FixedSizeArrays.AbstractArray)  ? length(X) : error("numel undefined for input of type $(T)")

# a macro to test if the contents of two containers are approximatley equal
macro contents_approx_eq(a, b)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test typeof(ai) == typeof(bi)
            @test_approx_eq ai bi
        end
    end
end
macro contents_approx_eq_eps(a, b, eps)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test typeof(ai) == typeof(bi)
            @test_approx_eq_eps ai bi $(esc(eps))
        end
    end
end

# a macro to test if tow types are approximatey equal
macro types_approx_eq(a, b)
    quote
        @test typeof($(esc(a))) == typeof($(esc(b)))
        @contents_approx_eq($(esc(a)), $(esc(b)))
    end
end
macro types_approx_eq_eps(a, b, eps)
    quote
        @test typeof($(esc(a))) == typeof($(esc(b)))
        @contents_approx_eq_eps($(esc(a)), $(esc(b)), $(esc(eps)))
    end
end

# a macro to test if the contents of two containers are approximatley equal, without examining the element types
macro contents_approx_eq_notype(a, b)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test_approx_eq ai bi
        end
    end
end
=#

#####################################################################################
# build a full list of rotation types including the different angle ordering schemas
#####################################################################################

rot_types = (RotMatrix, Quat, SPQuat, AngleAxis, RodriguesVec,
             RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
             RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ)

one_types = (RotX, RotY, RotZ)
two_types = (RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)
taitbyran_types = (RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX)
all_types = (RotMatrix, Quat, SPQuat, AngleAxis, RodriguesVec,
             RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
             RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
             RotX, RotY, RotZ,
             RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)

###############################
# Start testing
###############################

@testset "Rotations Tests" begin

    ###############################
    # Check fixed relationships
    ###############################

    @testset "Identity rotation checks" begin
        I = eye(SMatrix{3,3,Float64})
        I32 = eye(SMatrix{3,3,Float32})
        @testset "$(R)" for R in all_types
            @test eye(R)::R == I
            @test eye(R{Float32})::R{Float32} == I32
        end
    end

    ################################
    # check on the inverse function
    ################################

    @testset "Testing inverse()" begin
        repeats = 100
        I = eye(RotMatrix{Float64})
        @testset "$(R)" for R in all_types
            srand(0)
            for i = 1:repeats
                r = rand(R)
                @test inv(r) == r'
                @test inv(r) == r.'
                @test inv(r)*r ≈ I
                @test r*inv(r) ≈ I
            end
        end
    end

    #########################################################################
    # Rotate some stuff
    #########################################################################

    # a random rotation of a random point
    @testset "Rotate Points" begin
        repeats = 100
        @testset "$(R)" for R in all_types
            srand(0)
            for i = 1:repeats
                r = rand(R)
                m = SMatrix(r)
                v = randn(SVector{3})

                @test r*v ≈ m*v
            end

            # Test Base.Vector also
            r = rand(R)
            m = SMatrix(r)
            v = randn(3)

            @test r*v ≈ m*v
        end
    end

    # compose two random rotations
    @testset "Compose rotations" begin
        repeats = 100
        @testset "$(R1) * $(R2)" for R1 in all_types, R2 in all_types
            srand(0)
            for i = 1:repeats
                r1 = rand(R1)
                m1 = SMatrix(r1)

                r2 = rand(R2)
                m2 = SMatrix(r2)

                @test r1*r2 ≈ m1*m2
            end
        end
    end

    #########################################################################
    # Test conversions between rotation types
    #########################################################################
    @testset "Convert rotations" begin
        repeats = 100
        @testset "convert $(R1) -> $(R2)" for R1 in all_types, R2 in rot_types
            srand(0)
            for i = 1:repeats
                r1 = rand(R1)
                m1 = SMatrix(r1)

                r2 = R2(r1)

                @test r2 ≈ m1
            end
        end
    end

    #########################################################################
    # Check angle and axis and inv work as expected
    #########################################################################

    @testset "Testing angle / axis extraction" begin
        repeats = 100
        @testset "$(R)" for R in rot_types
            srand(0)
            for i = 1:repeats
                r1 = rand(AngleAxis)

                angle = rotation_angle(r1)
                axis = rotation_axis(r1)

                r2 = R(r1)

                @test rotation_angle(r2) ≈ angle
                @test rotation_axis(r2) ≈ axis
            end
        end

        # TODO RotX, RotXY?
    end

    #########################################################################
    # Check roll, pitch, yaw constructors
    #########################################################################

    @testset "Testing roll / pitch / yaw constructors" begin
        repeats = 100
        s = 1e-4
        @testset "$(R)" for R in taitbyran_types
            srand(0)
            for i = 1:repeats
                roll = s*(rand()-0.5)
                pitch = s*(rand()-0.5)
                yaw = s*(rand()-0.5)

                # This tests whether the rotations are the same to first order
                # in roll, pitch and yaw. Second-order terms are small enough
                # to pass the isapprox() test, but not first-order terms.
                r1 = RotXYZ(roll, pitch, yaw)
                r2 = R(roll=roll, pitch=pitch, yaw=yaw)

                @test r1 ≈ r2
            end
        end
    end
end


text:
"""
    PersistenceInterval

Type for representing persistence intervals. It behaves exactly like a `Tuple{Float64,
Float64}`, but can have meta data attached to it. The metadata is accessible with
`getproperty` or the dot syntax.

# Example

```jldoctest
julia> interval = PersistenceInterval(1, Inf; meta1=:a, meta2=:b)
[1.0, ∞) with:
 meta1: Symbol
 meta2: Symbol

julia> birth(interval), death(interval), persistence(interval)
(1.0, Inf, Inf)

julia> isfinite(interval)
false

julia> propertynames(interval)
(:birth, :death, :meta1, :meta2)

julia> interval.meta1
:a
```
"""
struct PersistenceInterval
    birth::Float64
    death::Float64
    meta::NamedTuple
end
function PersistenceInterval(birth, death; kwargs...)
    meta = (; kwargs...)
    return PersistenceInterval(Float64(birth), Float64(death), meta)
end
function PersistenceInterval(t::Tuple{<:Any,<:Any}; kwargs...)
    meta = (; kwargs...)
    return PersistenceInterval(Float64(t[1]), Float64(t[2]), meta)
end
function PersistenceInterval(int::PersistenceInterval; kwargs...)
    meta = (; kwargs...)
    return PersistenceInterval(Float64(int[1]), Float64(int[2]), meta)
end

"""
    birth(interval)

Get the birth time of `interval`.
"""
birth(int::PersistenceInterval) = getfield(int, 1)
"""
    death(interval)

Get the death time of `interval`.
"""
death(int::PersistenceInterval) = getfield(int, 2)
"""
    persistence(interval)

Get the persistence of `interval`, which is equal to `death - birth`.
"""
persistence(int::PersistenceInterval) = death(int) - birth(int)

"""
    midlife(interval)

Get the midlife of the `interval`, which is equal to `(birth + death) / 2`.
"""
midlife(int::PersistenceInterval) = (birth(int) + death(int)) / 2

Base.isfinite(int::PersistenceInterval) = isfinite(death(int))

###
### Iteration
###
function Base.iterate(int::PersistenceInterval, i=1)
    if i == 1
        return birth(int), i + 1
    elseif i == 2
        return death(int), i + 1
    else
        return nothing
    end
end

Base.length(::PersistenceInterval) = 2
Base.IteratorSize(::Type{<:PersistenceInterval}) = Base.HasLength()
Base.IteratorEltype(::Type{<:PersistenceInterval}) = Base.HasEltype()
Base.eltype(::Type{<:PersistenceInterval}) = Float64

function Base.getindex(int::PersistenceInterval, i)
    if i == 1
        return birth(int)
    elseif i == 2
        return death(int)
    else
        throw(BoundsError(int, i))
    end
end

Base.firstindex(int::PersistenceInterval) = 1
Base.lastindex(int::PersistenceInterval) = 2

###
### Equality and ordering
###
function Base.:(==)(int1::PersistenceInterval, int2::PersistenceInterval)
    return birth(int1) == birth(int2) && death(int1) == death(int2)
end
Base.:(==)(int::PersistenceInterval, (b, d)::Tuple) = birth(int) == b && death(int) == d
Base.:(==)((b, d)::Tuple, int::PersistenceInterval) = birth(int) == b && death(int) == d

function Base.isless(int1::PersistenceInterval, int2::PersistenceInterval)
    return (birth(int1), death(int1)) < (birth(int2), death(int2))
end

###
### Printing
###
function Base.show(io::IO, int::PersistenceInterval)
    b = round(birth(int); sigdigits=3)
    d = isfinite(death(int)) ? round(death(int); sigdigits=3) : "∞"
    return print(io, "[$b, $d)")
end

function Base.show(io::IO, ::MIME"text/plain", int::PersistenceInterval)
    b = round(birth(int); sigdigits=3)
    d = isfinite(death(int)) ? round(death(int); sigdigits=3) : "∞"
    print(io, "[$b, $d)")
    if !isempty(int.meta)
        print(io, " with:")
        for (k, v) in zip(keys(int.meta), int.meta)
            print(io, "\n ", k, ": ", summary(v))
        end
    end
end

###
### Metadata
###
function Base.getproperty(int::PersistenceInterval, key::Symbol)
    if hasfield(typeof(int), key)
        return getfield(int, key)
    elseif haskey(int.meta, key)
        return int.meta[key]
    else
        error("interval $int has no $key")
    end
end
function Base.propertynames(int::PersistenceInterval, private::Bool=false)
    if private
        return tuple(:birth, :death, propertynames(int.meta)..., :meta)
    else
        return (:birth, :death, propertynames(int.meta)...)
    end
end

"""
    representative(interval::PersistenceInterval)

Get the representative (co)cycle attached to `interval`, if it has one.
"""
representative(int::PersistenceInterval) = int.representative

"""
    birth_simplex(interval::PersistenceInterval)

Get the critical birth simplex of `interval`, if it has one.
"""
birth_simplex(int::PersistenceInterval) = int.birth_simplex

"""
    death_simplex(interval::PersistenceInterval)

Get the critical death simplex of `interval`, if it has one.

!!! note
    An infinite interval's death simplex is `nothing`.
"""
death_simplex(int::PersistenceInterval) = int.death_simplex


text:
"""
	bin_by_date!(fp::FluPop; start=:auto, last=:auto, binwidth=Day(121), binspacing=Day(121))

Bin `fp` by dates, from `start` to `last`.
"""
function bin_by_date!(fp::FluPop{S};
	start=:auto,
	stop=:auto,
	binwidth=Day(121),
	binspacing=Day(121)) where S

	fp.datebin = Dict{Tuple{Date,Date}, Array{S,1}}()
	# Start and end dates
	if start == :auto
		startdate = findmin([x[:date] for x in values(fp.strains)])[1]
	else
		startdate = Date(start)
	end
	if stop == :auto
		stopdate = findmax([x[:date] for x in values(fp.strains)])[1]
	else
		stopdate = Date(stop)
	end
	#
	now = startdate + binwidth
	while now < stopdate # If dlat - dstart != 0 mod[binwidth], the last bin will be smaller than others and is not considered here
		fp.datebin[now-binwidth, now] = Array{S,1}(undef, 0)
		now += binspacing
	end

	binlist = sort(collect(keys(fp.datebin)))
	for s in values(fp.strains)
		date = s[:date]
		if date < binlist[end][2] && date > binlist[1][1]
			db = _find_datebin(date, binlist)
			push!(fp.datebin[db], s)
		end
	end

	nothing
end

"""
Assume sorted datebins
"""
function _find_datebin(date::Date, datebins::Array{Tuple{Date,Date}})
	N = length(datebins) # Max
	i0 = 1 # Min
	i = div(N,2)
	found = false
	out = datebins[i]
	c = 0
	while !found
		c += 1
		if date >= datebins[i][1]
			if date < datebins[i][2]
				found = true
				out = datebins[i]
			else
				i0 = i
				i = i0 + div(N-i0,2) + 1
			end
		else
			N = i
			i = i0 + div(N-i0,2)
		end
	end
	return out
end

"""
	filter_by_region!(fp::FluPop, r::Array{<:AbstractString,1})
	filter_by_region!(fp::FluPop, r::AbstractString)
"""
function filter_by_region!(fp::FluPop, r::Array{<:AbstractString,1})
	for S in values(fp.datebin)
		idx = findall(s-> !in(s[:region], r), S)
		deleteat!(S, idx)
	end
	for s in values(fp.strains)
		if !in(s[:region], r)
			# println(s[:region])
			# println(r[1])
			# @show s[:region] == r[1]
			delete!(fp.strains, s[:strain])
		end
	end
end
filter_by_region!(fp::FluPop, r::AbstractString) = filter_by_region!(fp, [r])
"""
	filter_by_region(fp::FluPop, r)
"""
filter_by_region(fp::FluPop, r) = begin out = deepcopy(fp); filter_by_region!(out, r); return out end

"""
	filter_by_country!(fp::FluPop, r::Array{<:AbstractString,1})
	filter_by_country(fp::FluPop, r::AbstractString)
"""
function filter_by_country!(fp::FluPop, r::Array{<:AbstractString,1})
	for S in values(fp.datebin)
		idx = findall(s-> !in(s[:country], r), S)
		for i in idx
			delete!(fp.strains, S[i][:strain])
		end
		deleteat!(S, idx)
	end
end
filter_by_country(fp::FluPop, r::AbstractString) = filter_by_country!(fp, [r])
"""
	filter_by_country(fp::FluPop, r)
"""
filter_by_country(fp::FluPop, r) = begin out = deepcopy(fp); filter_by_country!(out, r); return out end

"""
	find_strains(S::Array{<:AbstractStrain}, i::Int64, val; seqtype = :aa)

Find strains in `S` that carry value `val` at position `i`.
"""
function find_strains(S::Array{<:AbstractStrain}, i::Int64, val; seqtype = :aa)
	return S[findall(s->s.seq[i] == val, S)]
end

"""
	datebin_to_date(d)
"""
function datebin_to_date(d)
	if length(d) != 2 || d[2] < d[1]
		@error "Invalid date bin"
	end
	return d[1] + div(d[2] - d[1], 2)
end
"""
	date_to_datebin(date, binwidth)
"""
function date_to_datebin(date, binwidth)
	return (date-binwidth, date+binwidth)
end
"""
	find_datebin(date::Date, datebins)
	find_datebin(date::Date, fp::FluPop)
"""
function find_datebin(date::Date, datebins)
	found = false
	out = first(datebins)
	for db in datebins
		if db[1] <= date && db[2] > date
			out = db
			found = true
			break
		end
	end
	if !found
		@warn "date $date was not found"
	end
	return out
end
find_datebin(date::Date, fp::FluPop) = find_datebin(date, keys(fp.datebin))

function get_regions(S::Array{<:Strain,1})
	out = Dict{String,Int64}()
	for s in S
		if !in(s[:region], ["?",'?'])
			out[s[:region]] = get(out, s[:region], 0) + 1
		end
	end
	return out
end

function numdate(d::Date)
	return year(d) + (month(d)-1) /12 + day(d) / 365
end

# Filters for IO
"""
"""
is_flu_outlier(x::AbstractStrain, flulineage) = in(x[:strain], outliers[flulineage])
function is_flu_outlier(flulineage)
	return x->is_flu_outlier(x,flulineage)
end

"""
	has_mutation(x::AbstractStrain, mutation)
	has_mutation(mutation)

Mutation `m` is expected to be of the form `(m[1], m[2])` (position, symbol).
"""
has_mutation(x::AbstractStrain, mutation) = (x.seq[mutation[1]] == mutation[2])
function has_mutation(mutation)
	return x->has_mutation(x, mutation)
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
using Trilinos
using Test
using MPI
import Trilinos.Testing.arrayview_sum

a = collect(1.0:10.0)
s = sum(a)
s_half = sum(a[1:5])

@test arrayview_sum(Teuchos.ArrayView(a)) == s
@test arrayview_sum(a) == s
@test arrayview_sum(a[1:5]) == s_half
@test arrayview_sum(@view a[1:5]) == s_half
@test arrayview_sum(Teuchos.ArrayView(a,5)) == s_half

function time_av(a)
  cumsum = 0.0
  for i in 1:10000
    cumsum += arrayview_sum(Teuchos.ArrayView(a,5))
  end
  return cumsum
end

@time time_av(a)
@time time_av(a)
@time time_av(a)