text:
# Routine to download automatically SMN and gauge data from the DWD using their R library
# Daniel Wolfensberger, LTE - EPFL / MeteoSwiss, 2020
# The routine is meant to be called in command line
# Rscript retrieve_dwh_data.r t0 t1 threshold stations variables output_folder missing_value overwrite
# t0 : start time in YYYYMMDDHHMM format
# t1 : end time in YYYYMMDDHHMM format
# threshold : minimum value of hourly precipitation for the entire hour to be included in the database (i.e. all 6 10min timesteps)
# stations : list of station abbrevations, e.g. "PAY, GVE, OTL"
# variables : list of variables to retrieve, using the DWH names, for example "tre200s0,prestas0,ure200s0,rre150z0,dkl010z0,fkl010z0"
# output_folder : directory where to store the csv files with the retrieved data
# overwrite : whether or not to overwrite already existing data in the output_folder

# IMPORTANT: this function is called by the main routine in database.py
# so you should never have to call it manually

.libPaths("/store/msclim/share/CATs/cats/lib-R3.5.2/")
library('mchdwh')
.libPaths("/store/msrad/utils/anaconda3-wolfensb/envs/radardb/r_libraries/")
library('R.utils')
library('plyr')
library('lubridate')

# This function adds new data to a file, with/without overwriting already present data
append_to_file <- function(fname, df, overwrite = FALSE){
  if(!file.exists(fname)){
    df_merged = df
  }
  else{
    df_old = data.table::fread(fname, sep = ',')
    if(!overwrite){
      not_in_old = !df$TIMESTAMP %in% df_old$TIMESTAMP
      df_merged = bind_rows(df_old, df[not_in_old,])
    }
    else{
      not_in_new = !df_old$TIMESTAMP %in% df$TIMESTAMP
      df_merged = data.table::rbindlist(list(df_old[not_in_new,],df), fill = T)
    }
    df_merged = df_merged[order(df_merged$TIMESTAMP),]
  }
  return(as.data.frame(df_merged))
}
#
# Parsing command line arguments
initial.options <- commandArgs(trailingOnly = FALSE)
file.arg.name <- "--file="
script.name <- sub(file.arg.name, "", initial.options[grep(file.arg.name, initial.options)])
script.basename <- paste(dirname(script.name),'/',sep = '')
print(initial.options)
args <-commandArgs(TRUE)

t0_str <- args[1]
tend_str <- args[2]
threshold <- args[3]
stations <- args[4]
variables <- args[5]
output_folder <- args[6]
missing_value <- as.numeric(args[7])
overwrite <- as.numeric(args[8])
variables = gsub("[()]","",variables)
variables = strsplit(variables,',')[[1]]
stations = gsub("[()]","",stations)
stations = strsplit(stations,',')[[1]]

#For testing

# t0_str <- '201808101410'
# tend_str <- '201808251000'
# threshold <- 0.1
# stations = c('MRP')
# variables = c('rre150z0')
# output_folder = '/scratch/wolfensb/dbase_tests/gauge/'
# script.basename = '/users/wolfensb/RADAR_GAUGE_DB/'
# missing_value <- -9999
# overwrite <- TRUE

# Get directory of data_stations.csv, i.e rainforest/common/constants
station_info_path = file.path(dirname(script.basename), 'common', 'data', 'data_stations.csv')
station_info = data.table::fread(station_info_path, sep = ';')
station_info <- as.data.frame(station_info)

for(i in 1:length(stations)){
  print(paste('Retrieving station ',stations[i]))
  isvalid = TRUE
  tryCatch({
  if((station_info[["Type"]][ station_info[["Abbrev"]] == stations[i] ]) != 'SwissMetNet'){
    variables_var <- c('rre150z0')
  }
  else{
    variables_var <- variables
  }
  }, error=function(e){print(e); isvalid = FALSE})

  if(!isvalid){
    next
  }
  tryCatch({
  data <- dwhget_surface(nat_abbr = stations[i], param_short = variables_var, date = c(t0_str,tend_str))
  # Reorganize data into multidimensional array
  data <- reshape_data(x=data, reshape="stp")

  # Fill up missing columns
  for(j in 1:length(variables)){
    if(!(variables[j] %in% colnames(data))){
      data[variables[j]] = missing_value
    }
  }

  # Time of beginning of measurement
  tstamps <- as.POSIXct(strptime(data$datetime,'%Y%m%d%H%M',tz='UTC')) - 60 * 5
  stahours_all <- paste(data[["nat_abbr"]],format(ceiling_date(tstamps,'hour'), format="%Y%m%d%H"))
  data[["stahours"]] = stahours_all

  # Get all wet hours
  rain_h <- aggregate(data[["rre150z0"]], by=list(data[["stahours"]]), FUN=sum, na.rm=TRUE)
  hours_wet <- rain_h[["x"]] >= threshold
  stahours_wet <- rain_h[["Group.1"]][hours_wet]

  zhours_wet <- stahours_all %in% stahours_wet
  zdata_wet <- data[zhours_wet,c('datetime',variables)]

  name_out <- paste(output_folder, '/', stations[i], '.csv.gz', sep = '')
  # date format must be changed to linux timestamp
  zdata_wet$datetime = as.numeric(as.POSIXct(strptime(zdata_wet$datetime,'%Y%m%d%H%M',tz='UTC')))
  colnames(zdata_wet)[1] = 'timestamp'

  # Finally since the R format of hdf5 is not compatible with dask and pandas
  # I choose to stay with csv
  zdata_wet[is.na(zdata_wet)] = missing_value

  # Add station info
  zdata_wet <- cbind(station = stations[i], zdata_wet)

  # Reorder if needed to always get same col order
  right_order = c('station','timestamp', variables)
  zdata_wet = zdata_wet[,right_order]
  colnames(zdata_wet) <- toupper(colnames(zdata_wet))

  zdata_wet = append_to_file(name_out, zdata_wet, overwrite)
  write.table(zdata_wet, gzfile(name_out), sep=',',col.names = T, row.names = F)

  }, error=function(e){print(e);print(paste('Data not available for ',stations[i],sep = ''))})
  }

max_stars_repo_name: gugerlir/rainforest

max_issues_repo_name: gugerlir/rainforest

max_forks_repo_name: gugerlir/rainforest


text:
#' cDriver.
#'
#' @name cDriver
#' @docType package
#' @description
#' \code{cDriver} R package for finding candidate driver genes in cancers.
#' Useful functions:
#' \itemize{
#'      \item \code{CCF}
#'      \item \code{bcgr.lawrence}
#'      \item \code{bcgr}
#'      \item \code{bcgr.combine}
#'      \item \code{bayes.risk}
#'      \item \code{bayes.driver}
#'      \item \code{combine.ranking}
#'      \item \code{cut.off.risk}
#'      \item \code{plotSamplesMut}
#'      \item \code{plotMutChange}
#'      \item \code{boxplotCCF.mutations}
#'      \item \code{boxplotCCF.patients}
#'      \item \code{plotStaircase}
#'      \item \code{boxplotCCF.mutations} }
#' Useful data:
#' \itemize{
#'      \item \code{lawrence.df}
#'      \item \code{length.genes}
#'      \item \code{driver.genes.concensus}
#'      \item \code{sample.genes.mutect}}
NULL

ext: r

max_stars_repo_name: hanasusak/cDriver

max_issues_repo_name: hanasusak/cDriver

max_forks_repo_name: hanasusak/cDriver


text:
library(glue)
library(here)
library(purrr)

# Week vector
weeks <- c(1:13)

# Create directories
purrr::map(weeks, ~dir.create(here("lecture_notes", glue("week{.}"))))


text:
# wahpenayo at gmail dot com
# 2017-12-10
#-----------------------------------------------------------------
if (file.exists('e:/porta/projects/taigabench')) {
  setwd('e:/porta/projects/taigabench')
} else {
  setwd('c:/porta/projects/taigabench')
}
source('src/scripts/r/functions.r')
readr.show_progress <- FALSE
#-----------------------------------------------------------------
dataset <- 'ontime'
problem <- 'qcost'
response <- 'arrdelay'
dataf <- ontime.data
dtest <- dataf(test.file(dataset=dataset))
suffixes <- c('32768','131072','524288','2097152','8388608','33554432')
#-----------------------------------------------------------------
bench(
  dataset=dataset,problem=problem,dataf=dataf,dtest=dtest,
  response=response,
  suffixes=suffixes,
  trainf=qcost.randomForestSRC,
  prefix='randomForestSRC')

bench(
  dataset=dataset,problem=problem,dataf=dataf,dtest=dtest,
  response=response,
  suffixes=suffixes,
  trainf=qcost.quantregForest,
  prefix='quantregForest')
#-----------------------------------------------------------------

max_stars_repo_name: wahpenayo/taigabench

max_issues_repo_name: wahpenayo/taigabench

max_forks_repo_name: wahpenayo/taigabench


text:
# Numeric vector: 1 up to 10
my_vector <- 1:10

# Numeric matrix: 1 up to 9
my_matrix <- matrix(1:9, ncol = 3)

# Factor of sizes
my_factor <- factor(c("M","S","L","L","M"), ordered = TRUE, levels = c("S","M","L"))

# List containing vector, matrix and factor
my_list <- list(my_vector, my_matrix, my_factor)

# Construct my_super_list with the four data structures above
my_super_list <- list(my_vector,my_matrix,my_factor,my_list)

# Display structure of my_super_list
str(my_super_list)

# Construct my_list with these different elements
my_list <- list(vec=my_vector, mat= my_matrix, fac=my_factor)

# Print my_list to the console
print(my_list)

max_stars_repo_name: applecool/DataScience

max_issues_repo_name: applecool/DataScience

max_forks_repo_name: applecool/DataScience


text:
fortyFivePlot <- function(input, output,values){
    prod = seq(0,input$ymax,length.out = 100)
    to_plot = data.frame(revenu = prod)
    fig = plot_ly(to_plot, x = ~revenu)

    fig = fig %>% add_trace(y = prod, type = "scatter", mode = "lines", name = "$$Y^d=y$$", color = I("blue"),
      hovertemplate = "y=yd=%{y:.0f}<extra></extra>")
    fig = fig %>% add_trace(y = (input$alpha + input$iy) * prod + input$id*input$D/values$p + input$ir*values$eq$r + input$cpi + input$bari + input$g, type = "scatter", mode = "lines", name = "$$Y_1^d$$", color = I("red"),
      hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
    if(values$eq$y > 0){
      fig = fig %>% add_segments(0,values$eq$y,values$eq$y,values$eq$y, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
        hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
      fig = fig %>% add_segments(values$eq$y,0,values$eq$y,values$eq$y, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
        hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
    }
    # fig = fig %>% add_trace(x = values$POTY, y = values$POTY, mode = "markers", marker = list(color = "red"), showlegend = F)

    if(values$shock & !is.na(input$new_value_IS)){
      fig = fig %>% add_trace(y = (values$shocked_params_IS["alpha"] + values$shocked_params_IS["iy"]) * prod + values$shocked_params_IS["id"]*input$D/values$shocked_params_LM["p"] + values$shocked_params_IS["ir"]*values$shocked_params_IS["r"] + values$shocked_params_IS["cpi"] + values$shocked_params_IS["bari"] + values$shocked_params_IS["g"], type = "scatter", mode = "lines", name = "$$Y_2^d$$", color = "rgb(0,200,20)",
        hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))

      if(values$new_eq$y > 0){
        fig = fig %>% add_segments(0,values$new_eq$y,values$new_eq$y,values$new_eq$y, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend=F,
          hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
        fig = fig %>% add_segments(values$new_eq$y,0,values$new_eq$y,values$new_eq$y, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend=F,
          hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))

        if(input$show_path){
          eq_tmp = values$eq$y
          for(i in 1:10){
            new_eq_tmp = unname((values$shocked_params_IS["alpha"]+values$shocked_params_IS["iy"])*eq_tmp + values$shocked_params_IS["id"]*input$D/values$shocked_params_LM["p"] + values$shocked_params_IS["ir"]*values$shocked_params_IS["r"]+values$shocked_params_IS["cpi"]+values$shocked_params_IS["bari"]+values$shocked_params_IS["g"])
            fig = fig %>% add_segments(eq_tmp,eq_tmp,eq_tmp,new_eq_tmp, line = list(color = 'rgb(150,150,150', width = 1, dash = "dash"), showlegend=F,
              hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
            fig = fig %>% add_segments(eq_tmp,new_eq_tmp,new_eq_tmp,new_eq_tmp, line = list(color = 'rgb(150,150,150', width = 1, dash = "dash"), showlegend=F,
              hovertemplate = paste("y=%{x:.0f}","
yd=%{y:.0f}","<extra></extra>"))
            eq_tmp = new_eq_tmp
          }
        }

        new_eq = list(
          x = values$new_eq$y,
          y = values$new_eq$y,
          text = paste0("y*=",round(values$new_eq$y))
        )
        fig = fig %>% layout(annotations = new_eq)
      }
    } else if(values$eq$y > 0){
      eq = list(
        x = values$eq$y,
        y = values$eq$y,
        text = paste0("y*=",round(values$eq$y))
      )
      fig = fig %>% layout(annotations = eq)
    }

    f <- list(
      family = "Courier New, monospace",
      size = 18,
      color = "#7f7f7f"
    )
    revenu <- list(
      title = "Revenu, y",
      titlefont = f
    )
    demand <- list(
      title = "Demande, Yd",
      titlefont = f
    )
    fig = fig %>% layout(xaxis = revenu, yaxis = demand)
}

ISPlot <- function(input,output,values){
  prod = seq(0,input$ymax,length.out = 100)
  to_plot = data.frame(revenu = prod)

  fig = plot_ly(to_plot, x = ~revenu)

  fig = fig %>% add_trace(y = ((1-input$alpha-input$iy)*prod-input$cpi-input$id*input$D/values$p -input$bari-input$g)/input$ir, type = "scatter", mode = "lines", name = "$$IS_1$$", color = I("red"),
                            hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
  if(values$eq$y > 0){
    fig = fig %>% add_segments(0,values$eq$r,values$eq$y,values$eq$r, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                 hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
    fig = fig %>% add_segments(values$eq$y,values$eq$r,values$eq$y,0, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                 hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
  }
  if(values$shock & !is.na(input$new_value_IS)){
    fig = fig %>% add_trace(y = ((1-values$shocked_params_IS["alpha"]-values$shocked_params_IS["iy"])*prod-values$shocked_params_IS["cpi"]-values$shocked_params_IS["bari"]-values$shocked_params_IS["g"])/values$shocked_params_IS["ir"], type = "scatter", mode = "lines", name = "$$IS_2$$", color = "rgb(0,150,20)",
                              hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
    if(values$new_eq$y > 0){
      fig = fig %>% add_segments(0,values$shocked_params_IS["r"],values$new_eq$y,values$shocked_params_IS["r"], line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                   hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
      fig = fig %>% add_segments(values$new_eq$y,values$shocked_params_IS["r"],values$new_eq$y,0, line = list(color = 'rgb(150, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                   hovertemplate = paste("y=%{x:.0f}","
r=%{y:.2f}","<extra></extra>"))
      is_eq = list(
        x = values$new_eq$y,
        y = values$shocked_params_IS["r"],
        text = paste0("y*=",round(values$new_eq$y))
      )
      fig = fig %>% layout(annotations = is_eq)
    }
  } else if(values$eq$y > 0){
    is_eq = list(
      x = values$eq$y,
      y = values$eq$r,
      text = paste0("y*=",round(values$eq$y))
    )
    fig = fig %>% layout(annotations = is_eq)
  }

  f <- list(
    family = "Courier New, monospace",
    size = 18,
    color = "#7f7f7f"
  )
  revenu <- list(
    title = "Revenu, y",
    titlefont = f
  )
  interest <- list(
    title = "Taux d'intérêt, r",
    titlefont = f,
    range = c(0,10)
  )
  fig = fig %>% layout(xaxis = revenu, yaxis = interest)

}

IeqSPlot <- function(input,output,values){
    prod = seq(0,input$ymax,length.out = 100)
    to_plot = data.frame(revenu = prod)

    invEpEq = (1-input$alpha)*values$eq$y-input$cpi
    fig = plot_ly(to_plot, x = ~revenu)
    fig = fig %>% add_trace(y = (1-input$alpha)*prod-input$cpi, type = "scatter", mode = "lines", name = "$$S_1(y)$$", line = list(color = "#00b70c"),
                              hovertemplate = paste("y=%{x:.0f}","
S=%{y:.0f}","<extra></extra>"))
    fig = fig %>% add_trace(y = input$iy*prod+input$id*input$D/values$p+input$ir*values$eq$r+input$bari+input$g, type = "scatter", mode = "lines", name = "$$I_1(y,r_0)$$", line = list(color = "#ab3bf2"),
                              hovertemplate = paste("y=%{x:.0f}","
I=%{y:.0f}","<extra></extra>"))
    if(invEpEq > 0){
      fig = fig %>% add_segments(0,invEpEq,values$eq$y,invEpEq, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                   hovertemplate = paste("y=%{x:.0f}","
I=S=%{y:.0f}","<extra></extra>"))
      fig = fig %>% add_segments(values$eq$y,invEpEq,values$eq$y,0, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                   hovertemplate = paste("y=%{x:.0f}","
I=S=%{y:.0f}","<extra></extra>"))
    }

    if(values$shock & !is.na(input$new_value_IS)){
      new_invEpEq = (1-values$shocked_params_IS["alpha"])*values$new_eq$y-values$shocked_params_IS["cpi"]
      if(input$shocked_var_IS == "cpi" | input$shocked_var_IS == "alpha"){
        fig = fig %>% add_trace(y = (1-values$shocked_params_IS["alpha"])*prod-values$shocked_params_IS["cpi"], type = "scatter", mode = "lines", name = "$$S_2(y)$$", line = list(color = "#006806"),
                                  hovertemplate = paste("y=%{x:.0f}","
S=%{y:.0f}","<extra></extra>"))
      } else{
        fig = fig %>% add_trace(y = values$shocked_params_IS["iy"]*prod+values$shocked_params_IS["ir"]*values$shocked_params_IS["r"]+values$shocked_params_IS["bari"]+values$shocked_params_IS["g"], type = "scatter", mode = "lines", name = "$$I_2(y,r_0)$$", line = list(color = "#581e7f"),
                                  hovertemplate = paste("y=%{x:.0f}","
I=%{y:.0f}","<extra></extra>"))
      }
      if(new_invEpEq > 0){
        fig = fig %>% add_segments(0,new_invEpEq,values$new_eq$y,new_invEpEq, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                     hovertemplate = paste("y=%{x:.0f}","
I=S=%{y:.0f}","<extra></extra>"))
        fig = fig %>% add_segments(values$new_eq$y,new_invEpEq,values$new_eq$y,0, line = list(color = 'rgb(200, 0, 0)', width = 1, dash = 'dash'), showlegend = FALSE,
                                     hovertemplate = paste("y=%{x:.0f}","
I=S=%{y:.0f}","<extra></extra>"))

        is_eq = list(
          x = values$new_eq$y,
          y = new_invEpEq,
          text = paste0("y*=",round(values$new_eq$y), ", I=S=",round(new_invEpEq))
        )
        fig = fig %>% layout(annotations = is_eq)
      }
    } else{
      invEp_eq = list(
        x = values$eq$y,
        y = invEpEq,
        text = paste0("y*=",round(values$eq$y),", I=S=",round(invEpEq))
      )
      fig = fig %>% layout(annotations = invEp_eq)
    }

    f <- list(
      family = "Courier New, monospace",
      size = 18,
      color = "#7f7f7f"
    )
    revenu <- list(
      title = "Revenu, y",
      titlefont = f
    )
    invEp <- list(
      title = "Investissement, Épargne, I,S",
      titlefont = f
    )
    fig = fig %>% layout(xaxis = revenu, yaxis = invEp)
}


text:
library(scales)
library(zoo)
library(dplyr)

# lighthouse, lodestar, nimbus, prysm, teku data
lh = read.csv("../dat/210128164829-lighthouse.csv", header = TRUE)
ld = read.csv("../dat/210128164829-lodestar.csv", header = TRUE)
nb = read.csv("../dat/210128164829-nimbus.csv", header = TRUE)
pr = read.csv("../dat/210128164829-prysm.csv", header = TRUE)
tk = read.csv("../dat/210128164829-teku.csv", header = TRUE)

# mutate outbound traffic to be negative
lh = lh %>% mutate(nout = out * -1)
ld = ld %>% mutate(nout = out * -1)
nb = nb %>% mutate(nout = out * -1)
pr = pr %>% mutate(nout = out * -1)
tk = tk %>% mutate(nout = out * -1)

# sum up inb and out traffic
lh = lh %>% mutate(trff = out + inc)
ld = ld %>% mutate(trff = out + inc)
nb = nb %>% mutate(trff = out + inc)
pr = pr %>% mutate(trff = out + inc)
tk = tk %>% mutate(trff = out + inc)

# orange, black, blue, purple, green gradients
col_lh = c("#7f2704","#a63603","#d94801","#f16913","#fd8d3c","#fdae6b","#fdd0a2","#fee6ce","#fff5eb")
col_ld = c("#000000","#252525","#525252","#737373","#969696","#bdbdbd","#d9d9d9","#f0f0f0","#ffffff")
col_nb = c("#08306b","#08519c","#2171b5","#4292c6","#6baed6","#9ecae1","#c6dbef","#deebf7","#f7fbff")
col_pr = c("#3f007d","#54278f","#6a51a3","#807dba","#9e9ac8","#bcbddc","#dadaeb","#efedf5","#fcfbfd")
col_tk = c("#00441b","#006d2c","#238b45","#41ab5d","#74c476","#a1d99b","#c7e9c0","#e5f5e0","#f7fcf5")

# time, slot, slots per second, disk usage, memory usage, peer count, traffic out/in, cpu usage
pairs(lh[,c(2, 3, 12, 5, 6, 7, 8, 9, 10)], main="Lighthouse", pch=".", col=alpha(col_lh[4], 0.3))
pairs(ld[,c(2, 3, 12, 5, 6, 7, 8, 9, 10)], main="Lodestar", pch=".", col=alpha(col_ld[4], 0.3))
pairs(nb[,c(2, 3, 12, 5, 6, 7, 8, 9, 10)], main="Nimbus", pch=".", col=alpha(col_nb[4], 0.3))
pairs(pr[,c(2, 3, 12, 5, 6, 7, 8, 9, 10)], main="Prysm", pch=".", col=alpha(col_pr[4], 0.3))
pairs(tk[,c(2, 3, 12, 5, 6, 7, 8, 9, 10)], main="Teku", pch=".", col=alpha(col_tk[4], 0.3))

# slot over time
sync <- plot(lh$time, lh$slot, main="Synchronization Progress over Time", xlab="Time [s]", ylab="Slot [1]", pch=".", xlim=c(0, 30027+5000), ylim=c(0, 424334+10000), col=alpha(col_lh[5], 0.3))
sync <- lines(pr$time, pr$slot, type="p", pch=".", col=alpha(col_pr[5], 0.3))
sync <- lines(tk$time, tk$slot, type="p", pch=".", col=alpha(col_tk[5], 0.3))
sync <- lines(nb$time, nb$slot, type="p", pch=".", col=alpha(col_nb[5], 0.3))
sync <- lines(ld$time, ld$slot, type="p", pch=".", col=alpha(col_ld[5], 0.3))
sync <- lines(lh$time, rollmean(lh$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
sync <- lines(pr$time, rollmean(pr$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
sync <- lines(tk$time, rollmean(tk$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
sync <- lines(nb$time, rollmean(nb$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
sync <- lines(ld$time, rollmean(ld$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
sync <- text(9029-1000, 419493+10000, "2h 30min", cex=0.8, pos=4, col=col_lh[2])
sync <- text(14759-1000, 419970+10000, "4h 5min", cex=0.8, pos=4, col=col_pr[2])
sync <- text(7739-0, 0+10000, "Prysm fetching genesis state (#8209)", cex=0.8, pos=4, col=col_pr[2])
sync <- text(18059-1500, 420243+10000, "5h 59sec", cex=0.8, pos=4, col=col_tk[2])
sync <- text(19213+0, 420341+10000, "5h 20min", cex=0.8, pos=4, col=col_nb[2])
sync <- text(30027-1000, 421243+10000, "8h 20min", cex=0.8, pos=4, col=col_ld[2])
sync <- text(24243+1000, 328097+0, "Lodestar JS-heap out-of-memory (#2005)", cex=0.8, pos=4, col=col_ld[2])
sync <- legend("bottomright", pch=19,
               legend = c("Lighthouse", "Prysm", "Teku", "Nimbus", "Lodestar"),
               col=c(col=col_lh[4],col=col_pr[4],col=col_tk[4],col=col_nb[4],col=col_ld[4]))

# slot over vtime (adjusted)
vsync <- plot(lh$vtime, lh$slot, main="Synchronization Progress with Adjusted Start-Time", xlab="Time [s]", ylab="Slot [1]", pch=".", xlim=c(0, 30027+5000), ylim=c(0, 424334+10000), col=alpha(col_lh[5], 0.3))
vsync <- lines(pr$vtime, pr$slot, type="p", pch=".", col=alpha(col_pr[5], 0.3))
vsync <- lines(tk$vtime, tk$slot, type="p", pch=".", col=alpha(col_tk[5], 0.3))
vsync <- lines(nb$vtime, nb$slot, type="p", pch=".", col=alpha(col_nb[5], 0.3))
vsync <- lines(ld$vtime, ld$slot, type="p", pch=".", col=alpha(col_ld[5], 0.3))
vsync <- lines(lh$vtime, rollmean(lh$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
vsync <- lines(pr$vtime, rollmean(pr$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
vsync <- lines(tk$vtime, rollmean(tk$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
vsync <- lines(nb$vtime, rollmean(nb$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
vsync <- lines(ld$vtime, rollmean(ld$slot, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
vsync <- text(14759-7739-2000, 419970+10000, "1h 57min", cex=0.8, pos=4, col=col_pr[2])
vsync <- text(0+1000, 0+10000, "Start time at first slot synced (#8209)", cex=0.8, pos=4, col=col_pr[2])
vsync <- text(9029-0-1000, 419493+10000, "2h 30min", cex=0.8, pos=4, col=col_lh[2])
vsync <- text(18059-38-1500, 420243+10000, "5h 21sec", cex=0.8, pos=4, col=col_tk[2])
vsync <- text(19213-5+0, 420341+10000, "5h 20min", cex=0.8, pos=4, col=col_nb[2])
vsync <- text(30027-6-1000, 421243+10000, "8h 20min", cex=0.8, pos=4, col=col_ld[2])
vsync <- text(24243+1000, 328097+0, "Lodestar JS-heap out-of-memory (#2005)", cex=0.8, pos=4, col=col_ld[2])
vsync <- legend("bottomright", pch=19,
               legend = c("Prysm", "Lighthouse", "Teku", "Nimbus", "Lodestar"),
               col=c(col=col_pr[4],col=col_lh[4],col=col_tk[4],col=col_nb[4],col=col_ld[4]))

# sps60 over time
sps <- plot(lh$time, lh$sps60, main="Synchronization Speed over Time", xlab="Time [s]", ylab="Slots per Second [1/s]", log="y", xlim=c(0, 30027+5000), ylim=c(1, 150), pch=".", col=alpha(col_lh[5], 0.3))
sps <- lines(pr$time, pr$sps60, type="p", pch=".", col=alpha(col_pr[5], 0.3))
sps <- lines(tk$time, tk$sps60, type="p", pch=".", col=alpha(col_tk[5], 0.3))
sps <- lines(nb$time, nb$sps60, type="p", pch=".", col=alpha(col_nb[5], 0.3))
sps <- lines(ld$time, ld$sps60, type="p", pch=".", col=alpha(col_ld[5], 0.3))
sps <- lines(lh$time[0:8954], rollmean(lh$sps60[0:8954], 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
sps <- lines(pr$time[6902:13809], rollmean(pr$sps60[6902:13809], 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
sps <- lines(tk$time[38:17955], rollmean(tk$sps60[38:17955], 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
sps <- lines(nb$time[4:19113], rollmean(nb$sps60[4:19113], 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
sps <- lines(ld$time[6:22501], rollmean(ld$sps60[6:22501], 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.2))
sps <- lines(ld$time[22614:28279], rollmean(ld$sps60[22614:28279], 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.2))
# > 419970/(14759-7739)
# [1] 59.82478632478632478632
sps <- text(14759-1500, 59.825, "59.825/s", cex=0.8, pos=4, col=col_pr[2])
# > 419493/(9029-0)
# [1] 46.46062686897773839849
sps <- text(9029-1500, 46.461, "46.461/s", cex=0.8, pos=4, col=col_lh[2])
# > 420243/(18059-38)
# [1] 23.31962710171466622274
sps <- text(18059-1500, 23.320, "23.320/s", cex=0.8, pos=4, col=col_tk[2])
# > 420341/(19213-5)
# [1] 21.88364223240316534777
sps <- text(19213-0, 21.884, "21.884/s", cex=0.8, pos=4, col=col_nb[2])
# > 357502/(24145-6)
# [1] 14.81014126517254235884
sps <- text(24145-1500, 14.810, "14.810/s", cex=0.8, pos=4, col=col_ld[2])
# > (421243-328097)/(30027-24243)
# [1] 16.10408022130013831259
sps <- text(30027-0, 16.104, "16.104/s (after restart)", cex=0.8, pos=4, col=col_ld[2])
sps <- legend("bottomright", pch=19,
                legend = c("Prysm", "Lighthouse", "Teku", "Nimbus", "Lodestar"),
                col=c(col=col_pr[4],col=col_lh[4],col=col_tk[4],col=col_nb[4],col=col_ld[4]))

# sps60 over slot
spsl <- plot(lh$slot, lh$sps60, main="Synchronization Speed over Sync Progress", xlab="Slot [1]", ylab="Slots per Second [1/s]", log="y", xlim=c(0, 424334+10000), ylim=c(1, 150), pch=".", col=alpha(col_lh[5], 0.3))
spsl <- lines(pr$slot, pr$sps60, type="p", pch=".", col=alpha(col_pr[5], 0.3))
spsl <- lines(tk$slot, tk$sps60, type="p", pch=".", col=alpha(col_tk[5], 0.3))
spsl <- lines(nb$slot, nb$sps60, type="p", pch=".", col=alpha(col_nb[5], 0.3))
spsl <- lines(ld$slot, ld$sps60, type="p", pch=".", col=alpha(col_ld[5], 0.3))
spsl <- lines(lh$slot[0:8954], rollmean(lh$sps60[0:8954], 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
spsl <- lines(pr$slot[6902:13809], rollmean(pr$sps60[6902:13809], 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
spsl <- lines(tk$slot[38:17955], rollmean(tk$sps60[38:17955], 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
spsl <- lines(nb$slot[4:19113], rollmean(nb$sps60[4:19113], 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
spsl <- lines(ld$slot[6:22501], rollmean(ld$sps60[6:22501], 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.1))
spsl <- lines(ld$slot[22614:28279], rollmean(ld$sps60[22614:28279], 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.1))
# > 419970/(14759-7739)
# [1] 59.82478632478632478632
spsl <- text(424334, 59.825, "59.825/s", cex=0.8, pos=4, col=col_pr[2])
# > 419493/(9029-0)
# [1] 46.46062686897773839849
spsl <- text(424334, 46.461, "46.461/s", cex=0.8, pos=4, col=col_lh[2])
# > 420243/(18059-38)
# [1] 23.31962710171466622274
spsl <- text(424334, 24, "23.320/s", cex=0.8, pos=4, col=col_tk[2])
# > 420341/(19213-5)
# [1] 21.88364223240316534777
spsl <- text(424334, 21, "21.884/s", cex=0.8, pos=4, col=col_nb[2])
# > 357502/(24145-6)
# [1] 14.81014126517254235884
spsl <- text(328097-10000, 14.810, "14.810/s", cex=0.8, pos=4, col=col_ld[2])
# > (421243-328097)/(30027-24243)
# [1] 16.10408022130013831259
spsl <- text(424334, 16.104, "16.104/s", cex=0.8, pos=4, col=col_ld[2])
spsl <- legend("bottomright", pch=19,
              legend = c("Prysm", "Lighthouse", "Teku", "Nimbus", "Lodestar"),
              col=c(col=col_pr[4],col=col_lh[4],col=col_tk[4],col=col_nb[4],col=col_ld[4]))

# db over time
dbt <- plot(tk$time, tk$db, main="Disk Usage over Time", xlab="Time [s]", ylab="Database Size [B]", pch=".", log="y", xlim=c(0, 30027+5000), ylim=c(50000000, 7000000000), col=alpha(col_tk[5], 0.3))
dbt <- lines(pr$time, pr$db, type="p", pch=".", col=alpha(col_pr[5], 0.3))
dbt <- lines(lh$time, lh$db, type="p", pch=".", col=alpha(col_lh[5], 0.3))
dbt <- lines(nb$time, nb$db, type="p", pch=".", col=alpha(col_nb[5], 0.3))
dbt <- lines(ld$time, ld$db, type="p", pch=".", col=alpha(col_ld[5], 0.3))
dbt <- lines(lh$time, rollmean(lh$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
dbt <- lines(pr$time, rollmean(pr$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
dbt <- lines(tk$time, rollmean(tk$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
dbt <- lines(nb$time, rollmean(nb$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
dbt <- lines(ld$time, rollmean(ld$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
dbt <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
dbt <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
dbt <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
dbt <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
dbt <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
# > max(tk$db[0:30027])
# [1] 5001367200
# > max(tk$db[0:30027])/1024/1024/1024
# [1] 4.657886
dbt <- text(15962, 5001367200+200000000, "4.658 GiB", cex=0.8, pos=4, col=col_tk[2])
# > max(nb$db[0:30027])
# [1] 4822584504
# > max(nb$db[0:30027])/1024/1024/1024
# [1] 4.491382
dbt <- text(30098, 4822584504+300000000, "4.491 GiB", cex=0.8, pos=4, col=col_nb[2])
# > max(pr$db[0:30027])
# [1] 3393183744
# > max(pr$db[0:30027])/1024/1024/1024
# [1] 3.160149
dbt <- text(29494, 3393183744+200000000, "3.160 GiB", cex=0.8, pos=4, col=col_pr[2])
# > max(lh$db[0:30027])
# [1] 3194929693
# > max(lh$db[0:30027])/1024/1024/1024
# [1] 2.97551
dbt <- text(8936, 3194929693-50000000, "2.976 GiB", cex=0.8, pos=4, col=col_lh[2])
# > max(ld$db[0:30027])
# [1] 1599823563
# > max(ld$db[0:30027])/1024/1024/1024
# [1] 1.489952
dbt <- text(30996, 1599823563+50000000, "1.490 GiB", cex=0.8, pos=4, col=col_ld[2])
dbt <- legend("bottomright", pch=19,
             legend = c("Teku", "Nimbus", "Prysm", "Lighthouse", "Lodestar"),
             col=c(col=col_tk[4],col=col_nb[4],col=col_pr[4],col=col_lh[4],col=col_ld[4]))

# db over slot
dbs <- plot(tk$slot, tk$db, main="Disk Usage over Sync Progress", xlab="Slot [1]", ylab="Database Size [B]", pch=".", log="y", xlim=c(0, 424334+20000), ylim=c(50000000, 7000000000), col=alpha(col_tk[5], 0.3))
dbs <- lines(pr$slot, pr$db, type="p", pch=".", col=alpha(col_pr[5], 0.3))
dbs <- lines(lh$slot, lh$db, type="p", pch=".", col=alpha(col_lh[5], 0.3))
dbs <- lines(nb$slot, nb$db, type="p", pch=".", col=alpha(col_nb[5], 0.3))
dbs <- lines(ld$slot, ld$db, type="p", pch=".", col=alpha(col_ld[5], 0.3))
dbs <- lines(lh$slot, rollmean(lh$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
dbs <- lines(pr$slot, rollmean(pr$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
dbs <- lines(tk$slot, rollmean(tk$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
dbs <- lines(nb$slot, rollmean(nb$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
dbs <- lines(ld$slot, rollmean(ld$db, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.1))
# > max(tk$db[0:30027])
# [1] 5001367200
# > max(tk$db[0:30027])/1024/1024/1024
# [1] 4.657886
dbs <- text(424334, 5001367200+300000000, "4.658 GiB", cex=0.8, pos=4, col=col_tk[2])
# > max(nb$db[0:30027])
# [1] 4822584504
# > max(nb$db[0:30027])/1024/1024/1024
# [1] 4.491382
dbs <- text(424334, 4822584504-100000000, "4.491 GiB", cex=0.8, pos=4, col=col_nb[2])
# > max(pr$db[0:30027])
# [1] 3393183744
# > max(pr$db[0:30027])/1024/1024/1024
# [1] 3.160149
dbs <- text(424334, 3393183744+100000000, "3.160 GiB", cex=0.8, pos=4, col=col_pr[2])
# > max(lh$db[0:30027])
# [1] 3194929693
# > max(lh$db[0:30027])/1024/1024/1024
# [1] 2.97551
dbs <- text(424334, 3194929693-50000000, "2.976 GiB", cex=0.8, pos=4, col=col_lh[2])
# > max(ld$db[0:30027])
# [1] 1599823563
# > max(ld$db[0:30027])/1024/1024/1024
# [1] 1.489952
dbs <- text(424334, 1599823563+0, "1.490 GiB", cex=0.8, pos=4, col=col_ld[2])
dbs <- legend("bottomright", pch=19,
             legend = c("Teku", "Nimbus", "Prysm", "Lighthouse", "Lodestar"),
             col=c(col=col_tk[4],col=col_nb[4],col=col_pr[4],col=col_lh[4],col=col_ld[4]))

# mem over time
memt <- plot(tk$time, tk$mem, main="Memory Usage over Time", xlab="Time [s]", ylab="Resident Memory [B]", pch=".", log="y", xlim=c(0, 30027+5000), ylim=c(100000000, 13000000000), col=alpha(col_tk[5], 0.3))
memt <- lines(pr$time, pr$mem, type="p", pch=".", col=alpha(col_pr[5], 0.3))
memt <- lines(lh$time, lh$mem, type="p", pch=".", col=alpha(col_lh[5], 0.3))
memt <- lines(nb$time, nb$mem, type="p", pch=".", col=alpha(col_nb[5], 0.3))
memt <- lines(ld$time, ld$mem, type="p", pch=".", col=alpha(col_ld[5], 0.3))
memt <- lines(lh$time, rollmean(lh$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
memt <- lines(pr$time, rollmean(pr$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
memt <- lines(tk$time, rollmean(tk$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
memt <- lines(nb$time, rollmean(nb$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
memt <- lines(ld$time, rollmean(ld$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
memt <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
memt <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
memt <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
memt <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
memt <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
# > max(tk$mem[0:30027])
# [1] 12298780672
# > max(tk$mem[0:30027])/1024/1024/1024
# [1] 11.45413
memt <- text(30027, 12298780672+500000000, "11.454 GiB", cex=0.8, pos=4, col=col_tk[2])
# > max(pr$mem[0:30027])
# [1] 6228918272
# > max(pr$mem[0:30027])/1024/1024/1024
# [1] 5.801132
memt <- text(30027, 6228918272+300000000, "5.801 GiB", cex=0.8, pos=4, col=col_pr[2])
# > max(ld$mem[0:22501])
# [1] 3126464512
# > max(ld$mem[0:22501])/1024/1024/1024
# [1] 2.911747
memt <- text(24145, 3126464512+300000000, "2.912 GiB", cex=0.8, pos=4, col=col_ld[2])
# > max(ld$mem[22614:28279])
# [1] 2311073792
# > max(ld$mem[22614:28279])/1024/1024/1024
# [1] 2.152355
memt <- text(30027, 2311073792+0, "2.152 GiB", cex=0.8, pos=4, col=col_ld[2])
# > max(lh$mem[0:30027])
# [1] 2900750336
# > max(lh$mem[0:30027])/1024/1024/1024
# [1] 2.701534
memt <- text(30027, 2900750336+200000000, "2.702 GiB", cex=0.8, pos=4, col=col_lh[2])
# > max(nb$mem[0:30027])
# [1] 1092087808
# > max(nb$mem[0:30027])/1024/1024/1024
# [1] 1.017086
memt <- text(30027, 1092087808+100000000, "1.017 GiB", cex=0.8, pos=4, col=col_nb[2])
memt <- legend("bottomright", pch=19,
             legend = c("Teku", "Prysm", "Lodestar", "Lighthouse", "Nimbus"),
             col=c(col=col_tk[4],col=col_pr[4],col=col_ld[4],col=col_lh[4],col=col_nb[4]))

# mem over slot
mems <- plot(tk$slot, tk$mem, main="Memory Usage over Sync Progress", xlab="Slot [1]", ylab="Resident Memory [B]", pch=".", log="y", xlim=c(0, 424334+20000), ylim=c(100000000, 13000000000), col=alpha(col_tk[5], 0.3))
mems <- lines(pr$slot, pr$mem, type="p", pch=".", col=alpha(col_pr[5], 0.3))
mems <- lines(lh$slot, lh$mem, type="p", pch=".", col=alpha(col_lh[5], 0.3))
mems <- lines(nb$slot, nb$mem, type="p", pch=".", col=alpha(col_nb[5], 0.3))
mems <- lines(ld$slot, ld$mem, type="p", pch=".", col=alpha(col_ld[5], 0.3))
mems <- lines(lh$slot, rollmean(lh$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
mems <- lines(pr$slot, rollmean(pr$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
mems <- lines(tk$slot, rollmean(tk$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
mems <- lines(nb$slot, rollmean(nb$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
mems <- lines(ld$slot, rollmean(ld$mem, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.1))
# > max(tk$mem[0:30027])
# [1] 12298780672
# > max(tk$mem[0:30027])/1024/1024/1024
# [1] 11.45413
mems <- text(424334, 12298780672, "11.454 GiB", cex=0.8, pos=4, col=col_tk[2])
# > max(pr$mem[0:30027])
# [1] 6228918272
# > max(pr$mem[0:30027])/1024/1024/1024
# [1] 5.801132
mems <- text(424334, 6228918272, "5.801 GiB", cex=0.8, pos=4, col=col_pr[2])
# > max(ld$mem[0:22501])
# [1] 3126464512
# > max(ld$mem[0:22501])/1024/1024/1024
# [1] 2.911747
mems <- text(357502, 3126464512, "2.912 GiB", cex=0.8, pos=4, col=col_ld[2])
# > max(ld$mem[22614:28279])
# [1] 2311073792
# > max(ld$mem[22614:28279])/1024/1024/1024
# [1] 2.152355
mems <- text(424334, 2311073792, "2.152 GiB", cex=0.8, pos=4, col=col_ld[2])
# > max(lh$mem[0:30027])
# [1] 2900750336
# > max(lh$mem[0:30027])/1024/1024/1024
# [1] 2.701534
mems <- text(424334, 2900750336, "2.702 GiB", cex=0.8, pos=4, col=col_lh[2])
# > max(nb$mem[0:30027])
# [1] 1092087808
# > max(nb$mem[0:30027])/1024/1024/1024
# [1] 1.017086
mems <- text(424334, 1092087808, "1.017 GiB", cex=0.8, pos=4, col=col_nb[2])
mems <- legend("bottomright", pch=19,
             legend = c("Teku", "Prysm", "Lodestar", "Lighthouse", "Nimbus"),
             col=c(col=col_tk[4],col=col_pr[4],col=col_ld[4],col=col_lh[4],col=col_nb[4]))

# pc over time
pert <- plot(nb$time, nb$pc, main="P2P Connections over Time", xlab="Time [s]", ylab="Peers [1]", pch=".", ylim=c(-20,170), xlim=c(0, 30027+5000), col=alpha(col_nb[5], 0.3))
pert <- lines(pr$time, pr$pc, type="p", pch=".", col=alpha(col_pr[5], 0.3))
pert <- lines(lh$time, lh$pc, type="p", pch=".", col=alpha(col_lh[5], 0.3))
pert <- lines(tk$time, tk$pc, type="p", pch=".", col=alpha(col_tk[5], 0.3))
pert <- lines(ld$time, ld$pc, type="p", pch=".", col=alpha(col_ld[5], 0.3))
pert <- lines(lh$time, rollmean(lh$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
pert <- lines(pr$time, rollmean(pr$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
pert <- lines(tk$time, rollmean(tk$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
pert <- lines(nb$time, rollmean(nb$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
pert <- lines(ld$time, rollmean(ld$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
pert <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
pert <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
pert <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
pert <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
pert <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
pert <- legend("bottomright", pch=19,
               legend = c("Nimbus", "Teku", "Lighthouse", "Prysm", "Lodestar"),
               col=c(col=col_nb[4],col=col_tk[4],col=col_lh[4],col=col_pr[4],col=col_ld[4]))

# pc over slot
pers <- plot(nb$slot, nb$pc, main="P2P Connections over Sync Progress", xlab="Slot [1]", ylab="Peers [1]", pch=".", ylim=c(-20,170), xlim=c(0, 424334+20000), col=alpha(col_nb[5], 0.3))
pers <- lines(pr$slot, pr$pc, type="p", pch=".", col=alpha(col_pr[5], 0.3))
pers <- lines(lh$slot, lh$pc, type="p", pch=".", col=alpha(col_lh[5], 0.3))
pers <- lines(tk$slot, tk$pc, type="p", pch=".", col=alpha(col_tk[5], 0.3))
pers <- lines(ld$slot, ld$pc, type="p", pch=".", col=alpha(col_ld[5], 0.3))
pers <- lines(lh$slot, rollmean(lh$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
pers <- lines(pr$slot, rollmean(pr$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
pers <- lines(tk$slot, rollmean(tk$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
pers <- lines(nb$slot, rollmean(nb$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
pers <- lines(ld$slot, rollmean(ld$pc, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.1))
pers <- legend("bottomright", pch=19,
               legend = c("Nimbus", "Teku", "Lighthouse", "Prysm", "Lodestar"),
               col=c(col=col_nb[4],col=col_tk[4],col=col_lh[4],col=col_pr[4],col=col_ld[4]))

# out and inbount traffic
trff <- plot(nb$time, nb$inc, main="Traffic over Time (In/Out)", xlab="Time [s]", ylab="Traffic [B/s]", pch=".", ylim=c(-80000, 80000), xlim=c(0, 30027+5000), col=alpha(col_nb[5], 0.1))
trff <- lines(pr$time, pr$inc, type="p", pch=".", col=alpha(col_pr[5], 0.1))
trff <- lines(lh$time, lh$inc, type="p", pch=".", col=alpha(col_lh[5], 0.1))
trff <- lines(tk$time, tk$inc, type="p", pch=".", col=alpha(col_tk[5], 0.1))
trff <- lines(ld$time, ld$inc, type="p", pch=".", col=alpha(col_ld[5], 0.1))
trff <- lines(lh$time, rollmean(lh$inc, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
trff <- lines(pr$time, rollmean(pr$inc, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
trff <- lines(tk$time, rollmean(tk$inc, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
trff <- lines(nb$time, rollmean(nb$inc, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
trff <- lines(ld$time, rollmean(ld$inc, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
trff <- lines(pr$time, pr$nout, type="p", pch=".", col=alpha(col_pr[5], 0.1))
trff <- lines(pr$time, pr$nout, type="p", pch=".", col=alpha(col_pr[5], 0.1))
trff <- lines(lh$time, lh$nout, type="p", pch=".", col=alpha(col_lh[5], 0.1))
trff <- lines(tk$time, tk$nout, type="p", pch=".", col=alpha(col_tk[5], 0.1))
trff <- lines(ld$time, ld$nout, type="p", pch=".", col=alpha(col_ld[5], 0.1))
trff <- lines(lh$time, rollmean(lh$nout, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
trff <- lines(pr$time, rollmean(pr$nout, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
trff <- lines(tk$time, rollmean(tk$nout, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
trff <- lines(nb$time, rollmean(nb$nout, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
trff <- lines(ld$time, rollmean(ld$nout, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
trff <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
trff <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
trff <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
trff <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
trff <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
trff <- text(33000, +70000, "+ in", cex=0.8, pos=4, col=col_ld[1])
trff <- text(33000, -70000, "- out", cex=0.8, pos=4, col=col_ld[1])
trff <- legend("bottomleft", pch=19,
               legend = c("Nimbus", "Teku", "Prysm", "Lodestar", "Lighthouse"),
               col=c(col=col_nb[4],col=col_tk[4],col=col_pr[4],col=col_ld[4],col=col_lh[4]))

# out and inbount traffic
trfs <- plot(nb$time, nb$trff, main="Traffic over Time (Total)", xlab="Time [s]", ylab="Traffic [B/s]", pch=".", log="y", xlim=c(0, 30027+5000), col=alpha(col_nb[5], 0.1))
trfs <- lines(pr$time, pr$trff, type="p", pch=".", col=alpha(col_pr[5], 0.1))
trfs <- lines(lh$time, lh$trff, type="p", pch=".", col=alpha(col_lh[5], 0.1))
trfs <- lines(tk$time, tk$trff, type="p", pch=".", col=alpha(col_tk[5], 0.1))
trfs <- lines(ld$time, ld$trff, type="p", pch=".", col=alpha(col_ld[5], 0.1))
trfs <- lines(lh$time, rollmean(lh$trff, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
trfs <- lines(pr$time, rollmean(pr$trff, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
trfs <- lines(tk$time, rollmean(tk$trff, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
trfs <- lines(nb$time, rollmean(nb$trff, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
trfs <- lines(ld$time, rollmean(ld$trff, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
trfs <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
trfs <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
trfs <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
trfs <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
trfs <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
trfs <- legend("bottomright", pch=19,
               legend = c("Nimbus", "Teku", "Prysm", "Lodestar", "Lighthouse"),
               col=c(col=col_nb[4],col=col_tk[4],col=col_pr[4],col=col_ld[4],col=col_lh[4]))

# cpu over time
cput <- plot(lh$time, lh$cpu, main="CPU Usage over Time", xlab="Time [s]", ylab="User [%]", pch=".", ylim=c(0, 100), xlim=c(0, 30027+5000), col=alpha(col_lh[5], 0.3))
cput <- lines(pr$time, pr$cpu, type="p", pch=".", col=alpha(col_pr[5], 0.3))
cput <- lines(nb$time, nb$cpu, type="p", pch=".", col=alpha(col_nb[5], 0.3))
cput <- lines(tk$time, tk$cpu, type="p", pch=".", col=alpha(col_tk[5], 0.3))
cput <- lines(ld$time, ld$cpu, type="p", pch=".", col=alpha(col_ld[5], 0.3))
cput <- lines(lh$time, rollmean(lh$cpu, 600, fill=list(NA, NULL, NA)), col=alpha(col_lh[3], 0.5))
cput <- lines(pr$time, rollmean(pr$cpu, 600, fill=list(NA, NULL, NA)), col=alpha(col_pr[3], 0.5))
cput <- lines(tk$time, rollmean(tk$cpu, 600, fill=list(NA, NULL, NA)), col=alpha(col_tk[3], 0.5))
cput <- lines(nb$time, rollmean(nb$cpu, 600, fill=list(NA, NULL, NA)), col=alpha(col_nb[3], 0.5))
cput <- lines(ld$time, rollmean(ld$cpu, 600, fill=list(NA, NULL, NA)), col=alpha(col_ld[3], 0.5))
cput <- text(-1000, 0.0/8.0*100.0, "0 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 1.0/8.0*100.0, "1 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 2.0/8.0*100.0, "2 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 3.0/8.0*100.0, "3 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 4.0/8.0*100.0, "4 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 5.0/8.0*100.0, "5 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 6.0/8.0*100.0, "6 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 7.0/8.0*100.0, "7 -", cex=0.8, pos=4, col=col_ld[1])
cput <- text(-1000, 8.0/8.0*100.0, "8 - cores", cex=0.8, pos=4, col=col_ld[1])
cput <- abline(v = 9029, col=col_lh[3], lwd=1, lty=3)
cput <- abline(v = 14759, col=col_pr[3], lwd=1, lty=3)
cput <- abline(v = 18059, col=col_tk[3], lwd=1, lty=3)
cput <- abline(v = 19213, col=col_nb[3], lwd=1, lty=3)
cput <- abline(v = 30027, col=col_ld[3], lwd=1, lty=3)
cput <- legend("topright", pch=19,
               legend = c("Lighthouse", "Prysm", "Teku", "Lodestar", "Nimbus"),
               col=c(col=col_lh[4],col=col_pr[4],col=col_tk[4],col=col_ld[4],col=col_nb[4]))


text:
#!/usr/bin/env Rscript

# a script to process the output from hulk smash
# And vizualize them as a heatmap and a ordination plot

#functions
# check.packages function: install and load multiple R packages.
# (downloaded from: https://gist.github.com/smithdanielle/9913897)
# Check to see if packages are installed. Install them if they are not, then load them into the R session.
# NOTE CRAN MIRROR SELECT IS norway MIRROR !!!
check.packages <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE, repos="https://cran.uib.no/")
  sapply(pkg, require, character.only = TRUE)
}

#Libraries to load and install if needed.
packages<-c("heatmap3")
check.packages(packages)

# loading and processing the data
all_samples <- read.csv("all_samples.Weighted_Jaccard.hulk-matrix.csv", head=TRUE)
colnames(all_samples) <- sub(".clean.fq.gz.json","",colnames(all_samples))
rownames(all_samples) <- colnames(all_samples)

samples <- ncol(all_samples)

#transforming the data to a matrix
all_samples_matrix <- as.matrix(all_samples)

#creating a heatmap using the Jaccard distances, and saving it.
colors <- colorRampPalette(c("steelblue", "white","orange2"))(100)

pdf(file= "Heatmap.all_samples.WeightedJaccard.pdf",
    width = 12, height = 9, bg="white")

heatmap3(all_samples_matrix, scale="none", col=colors,cexRow=(0.8), cexCol=0.8, margins=c(10,6))
dev.off()


text:
# create a dialog box to select a file

newProgressBar = function(text = 'Progress bar')
{
  # adapted from http://metabolomics-forum.com/viewtopic.php?f=9&t=411

  pwindow <<- gwindow(container = TRUE, width = 400, height = 50, title = sprintf("%s (0%% done)", text))
  bp <<- gtkProgressBarNew()
  bp$'width-request' = 400
  bp$'height-request' = 50
  add(pwindow, bp)
}

updateProgressBar = function(p, text = 'Progress bar')
{
  bp$setFraction(p)
  svalue(pwindow) = sprintf("%s (%.1f%% done)", text, p*100)
}

fileChoose <- function(action="print", text = "Select a file...", type="open", ...)
{
  gfile(text=text, type=type, ..., action = action, # create dialog box
    handler =function(h,...) { # with handler for when action performed
      if(h$action != "") # if the action has not been set to blank
      {
        print(h$action) # print the action name to R
        setwd(dirname(h$file)) # set directory of file to new working directory
        do.call(h$action, list(h$file)) # perform the action on the file (usually print)
      }
    })
}

# File save dialog
fileSaveChoose <- function(action="print", text = "Select a file...", type="save", ...) {
  gfile(text=text, type=type, ..., action = action, # create dialog box
    handler =function(h,...) { # with handler for when action performed
      if(h$action != "")
      {
        print(h$action) # print the action name to R
        setwd(dirname(h$file)) # set directory of file to new working directory
        do.call(h$action, list(h$file)) # perform the action on the file (usually print)
      }
    })
}

## Simple confirmation dialog
confirmDialog <- function(message, handler=NULL) {
   window <- gwindow("Confirm")
   group <- ggroup(container = window)
   gimage("info", dirname="stock", size="dialog", container=group)

   ## A group for the message and buttons
   inner.group <- ggroup(horizontal=FALSE, container = group)
   glabel(message, container=inner.group, expand=TRUE)

   ## A group to organize the buttons
   button.group <- ggroup(container = inner.group)
   ## Push buttons to right
   addSpring(button.group)
   gbutton("OK", handler = function(h,...) dispose(window), container=button.group)
   return()
}

# routine to load a user-selected raster file,
# outputs a list of the necessary outputs
# the.file = name of loaded file
# the.rast = data in raster form
# the.data = data in matrix form
raster.load.button.handler = function()
{
  the.file=fileChoose("print") # select file using dialog
  if (!is.na(the.file)) # if a file has been selected
  {
    ext = extension(the.file) # work out its extension
    if (ext == '.csv')  # if it is a csv file, use read.csv (use try to catch errors)
      DATA = try(read.csv(the.file, header = FALSE), silent = TRUE)
    else   # otherwise, use raster (should work on all standard raster formats)
      DATA = try(raster(the.file), silent = TRUE)
    try.error = inherits(DATA, 'try-error')
    if (try.error) # if an error has been caught, complain in a dialog box (TODO: need further action?)
      confirmDialog("File not compatible! Please select file in either .csv or raster format")
    else # otherwise, convert data for output
    {
      if (ext == '.csv')
        the.rast = raster(as.matrix(DATA)) # make a raster out of the csv file
      else
        # read raster into memory (note: this is necessary as DATA itself
        # doesn't actually contain the raster data - only a reference
        # to the filename of the loaded raster)
        the.rast = readAll(DATA)
      the.data = raster::as.matrix(the.rast) # TODO: default version didn't work when SPADE made package for some reason...
    }
    # return list of required outputs
    list(the.file = the.file, the.rast = the.rast,  the.data = the.data)
  }
  else # if no valid file has been selected, return NULL
    NULL
}

# set global (i.e. terrain-based) dispersal and fencing information to defaults
# and plot
default.rasters = function(h,...)
{
  svalue(input.pad) <<- 1 # choose the Rasters tab
  # set EL.rast (raster information) to 1 everywhere
  # - using the first species' carrying capacity as a template
  EL.all.rast <<- K.rast[[1]]*0 + 1
  EL.all <<- raster::as.matrix(EL.all.rast) # set matrix version
  svalue(EL.all.label)='Default (1 for all cells)' # show default on label
  visible(EL.all.plot.raster)=TRUE  # select the plot device
  svalue(raster.pad)=1
  plotdata(EL.all) # plot data

  # set fenced.rast (raster information) to 1 everywhere
  # - using the first species' carrying capacity as a template
  fenced.rast <<- K.rast[[1]]*0
  fenced <<- raster::as.matrix(fenced.rast) # set matrix version
  svalue(fenced.label)='Default (0 for all cells)' # show default on label
  visible(fenced.plot.raster)=TRUE # select the plot device
  svalue(raster.pad)=2
  plotdata(fenced)  # plot data
}

# routine to update the strategy.params list based on the current values of the management text boxes
# (the spin buttons update the variables automatically)
update.strategy.sliders = function()
{
  if (!is.null(curr.strategy)) # only if has already been defined
  {
    if (curr.strategy <= N.strategies) # and is an existent option (added because of remove strategy)
    {
      # stored param values into strategy.params from numbers currently on gtext objects (using previous droplist place)
      strategy.params$init.cull.abs[curr.strategy] <<- as.double(svalue(init.cull.abs.slider))
      strategy.params$maint.cull.abs[curr.strategy] <<- as.double(svalue(maint.cull.abs.slider))
      strategy.params$target.density[curr.strategy] <<- as.double(svalue(target.density.slider))
      strategy.params$cb.a[curr.strategy] <<- as.double(svalue(cb.a.slider))
      strategy.params$cb.b[curr.strategy] <<- as.double(svalue(cb.b.slider))
      strategy.params$cb.c[curr.strategy] <<- as.double(svalue(cb.c.slider))
  }}
}

# function to generate Management window
create.man.window = function()
{
  strategies.group = ggroup(container=man.group) # group for Strategies selector
  a.label = glabel('Strategies:', container = strategies.group) # label
  addSpring(strategies.group)# make sure droplist is on right
  strategies.droplist <<- gdroplist(items = 'None', container = strategies.group, handler = sdh) # droplist to select current strategy

  sar.group = ggroup(container = man.group, horizontal = TRUE) # strategies add/remove buttons group

  # Add strategy button with handler for when button is pressed
  add.strategy = gbutton("Add strategy", container = sar.group,
    handler = add.strategy.handler) # handler in handlers.r

  # Remove strategy button with handler for when button is pressed
  rem.strategy = gbutton("Remove strategy", container = sar.group,
    handler = rem.strategy.handler) # handler in handlers.r

  # define group to contain all the parameters for each strategy
  strategy.params.group <<- ggroup(container = man.group, horizontal = FALSE)
  # set visibility depending on whether there are any strategies defined
  visible(strategy.params.group) = (N.strategies > 0)

  culling.group=ggroup(container=strategy.params.group) # create space for label and droplist
  a.label=glabel("Species to manage:",container=culling.group) # label
  addSpring(culling.group) # make sure droplist is aligned to the right
  culling.droplist <<- gdroplist(items = species, container=culling.group, selected = 1, # droplist for selecting species to manage with handler
    handler=function(h,...) # handles when option is changed
    {
      if (length(svalue(h$obj, index=TRUE))>0) # if it has a value
        strategy.params$cull.species[curr.strategy] <<- svalue(h$obj, index=TRUE) # set relevant cull.species option to current selected species
    }
  )

  seasons.group = ggroup(container = strategy.params.group) # create space for label and checkbox(es)
  a.label = glabel('Seasons to cull:', container = seasons.group) # label
  seasons.boxes <<- vector('list', seasons) # create list of checkbox objects (one checkboxbox for each season)
  for (i in 1:seasons) # for each season, make a checkbox
  {
    # create checkbox with handler for when checkbox is checked/unchecked
    seasons.boxes[[i]] <<- gcheckbox(i, checked = strategy.params$cull.seasons[curr.strategy, i], container = seasons.group,
      handler = function(h,...) strategy.params$cull.seasons[curr.strategy, id(h$obj)] <<- svalue(h$obj)) # handler sets relevant cull.seasons option depending on checked status
    id(seasons.boxes[[i]]) <<- i # sets the id tag for the checkbox so that above handler works correctly
  }

  area.target.group=ggroup(container=strategy.params.group) # create space for label and spinbutton
  area.target.label=glabel("Target proportional population in region (used for optimisation)",container=area.target.group) # label
  addSpring(area.target.group) # make sure spinbutton is aligned to the right
  # slider with handler to set value in strategy.params to spinbutton value when changed
  area.target.slider <<- gspinbutton(from=0.01,to=1,by=.01,value=strategy.params$area.target[curr.strategy],container=area.target.group,
    handler=function(h,...) strategy.params$area.target[curr.strategy] <<-svalue(area.target.slider))

  # Create a radio group to select the method of culling
  culling.frame = gframe("Removal method:", container = strategy.params.group, horizontal = FALSE) # label/frame
  culling.method <<- gradio(c("Proportional", "Absolute", "Capped proportional"), container = culling.frame,
    handler = function(h,...){ # handler for when method is changed
      a = svalue(h$obj, index = TRUE) # extract newly selected method
      strategy.params$culling.choice[curr.strategy] <<- a # add information to strategy.params
      visible(init.cull.group) = (a == 1 | a == 3)    # visible if "Proportional" is selected
      visible(maint.cull.group) = (a == 1 | a == 3)
      visible(init.cull.abs.group) = (a == 2 | a == 3) # visible if "Absolute" is selected
      visible(maint.cull.abs.group) = (a == 2 | a == 3)
  })

  init.cull.group<<-ggroup(container=culling.frame) # create space for label and spinbutton
  init.cull.label=glabel("Removal rate in first year",container=init.cull.group) # label
  addSpring(init.cull.group) # make sure spinbutton is aligned to the right
  # slider with handler to set value in strategy.params to spinbutton value when changed
  init.cull.slider <<- gspinbutton(from=0,to=0.99,by=.01,value=strategy.params$init.cull[curr.strategy],container=init.cull.group,
    handler=function(h,...) strategy.params$init.cull[curr.strategy] <<-svalue(init.cull.slider))

  maint.cull.group<<-ggroup(container=culling.frame) # create space for label and spinbutton
  maint.cull.label=glabel("Removal rate after first year",container=maint.cull.group) # label
  addSpring(maint.cull.group) # make sure spinbutton is aligned to the right
  # slider with handler to set value in strategy.params to spinbutton value when changed
  maint.cull.slider <<- gspinbutton(from=0,to=0.99,by=.01,value=strategy.params$maint.cull[curr.strategy],container=maint.cull.group,
    handler=function(h,...) strategy.params$maint.cull[curr.strategy] <<-svalue(maint.cull.slider))

  init.cull.abs.group<<-ggroup(container=culling.frame) # create space for label and edit box
  init.cull.abs.label=glabel("Total removal in first year",container=init.cull.abs.group) # label
  addSpring(init.cull.abs.group) # make sure edit box is aligned to the right
  init.cull.abs.slider<<-gedit(text=strategy.params$init.cull.abs[curr.strategy],container=init.cull.abs.group) # edit box
  visible(init.cull.abs.group) = FALSE # set to originally be invisible (as Proportional mode is selected by default)

  maint.cull.abs.group<<-ggroup(container=culling.frame) # create space for label and edit box
  maint.cull.abs.label=glabel("Total removal after first year",container=maint.cull.abs.group)# label
  addSpring(maint.cull.abs.group) # make sure edit box is aligned to the right
  maint.cull.abs.slider<<-gedit(text=strategy.params$maint.cull.abs[curr.strategy],container=maint.cull.abs.group) # edit box
  visible(maint.cull.abs.group) = FALSE # set to originally be invisible

  # create group for spatially implicit definitions
  # so that we can set its visibility depending on whether we're using rasters or not
  # (they start off as visible)
  spatially.implicit.group1 <<- ggroup(container = culling.frame, horizontal = FALSE)

  target.density.group=ggroup(container=spatially.implicit.group1) # create space for label and edit box
  target.density.label=glabel("Species density (per km^2) at which removal stops in a cell",container=target.density.group) # label
  addSpring(target.density.group) # make sure edit box is aligned to the right
  target.density.slider<<-gedit(text=strategy.params$target.density[curr.strategy],container=target.density.group) # edit box

  # create groups for spatially explicit definitions
  # so that we can set their visibility depending on whether we're using rasters or not
  # (they start off as invisible = value of use.cost.rasters)

  # first group for buttons (stacked horizontally to save space)
  cost.raster.group1 <<- ggroup(container=culling.frame, horizontal = TRUE)

  # buttons to load raster from file (spatially explicit)
  TD.button <- gbutton("Load Target Density Raster", container=cost.raster.group1,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata)) # if data loaded
      {
        # store data in strategy.params
        strategy.params$TD.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$TD[[curr.strategy]] <<- allthedata$the.data
        visible(TD.plot.raster)=TRUE # select plot window
        svalue(strats.raster.pad)=3
        plotdata(strategy.params$TD[[curr.strategy]]) # plot raster
      }
    })

  # Create a radio group to select the cost per animal method
  cb.frame = gframe("Cost per animal method:", container = strategy.params.group, horizontal = FALSE) # label/frame
  cb.method <<- gradio(c("Power", "Exponential"), container = cb.frame,
    handler = function(h,...){ # handler for when method is changed
      a = svalue(h$obj, index = TRUE) # extract newly selected method
      strategy.params$cb.choice[curr.strategy] <<- a # add information to strategy.params
      if (a == 1)
      {
        svalue(cb.a.label) = "Intercept"
        svalue(cb.b.label) = "Slope"
        svalue(cb.c.label) = "Hourly cost"
        svalue(CB.A.button) = "Load Intercept Raster"
        svalue(CB.B.button) = "Load Slope Raster"
        svalue(CB.C.button) = "Load Hourly Cost Raster"
      }
      else
      {
        svalue(cb.a.label) = "Intercept"
        svalue(cb.b.label) = "Coefficient"
        svalue(cb.c.label) = "Exponential rate"
        svalue(CB.A.button) = "Load Intercept Raster"
        svalue(CB.B.button) = "Load Coefficient Raster"
        svalue(CB.C.button) = "Load Exponential Rate Raster"
      }
  })

  # create group for spatially implicit definitions
  # so that we can set its visibility depending on whether we're using rasters or not
  # (they start off as visible)
  spatially.implicit.group2 <<- ggroup(container = cb.frame, horizontal = FALSE)

  cb.a.group=ggroup(container=spatially.implicit.group2) # create space for label and edit box
  cb.a.label<<-glabel("Intercept",container=cb.a.group) # label
  addSpring(cb.a.group) # make sure edit box is aligned to the right
  cb.a.slider<<-gedit(text=strategy.params$cb.a[curr.strategy],container=cb.a.group) # edit box

  cb.b.group=ggroup(container=spatially.implicit.group2) # create space for label and edit box
  cb.b.label<<-glabel("Slope",container=cb.b.group) # label
  addSpring(cb.b.group) # make sure edit box is aligned to the right
  cb.b.slider<<-gedit(text=strategy.params$cb.b[curr.strategy],container=cb.b.group) # edit box

  cb.c.group=ggroup(container=spatially.implicit.group2) # create space for label and edit box
  cb.c.label<<-glabel("Hourly cost",container=cb.c.group) # label
  addSpring(cb.c.group) # make sure edit box is aligned to the right
  cb.c.slider<<-gedit(text=strategy.params$cb.c[curr.strategy],container=cb.c.group)  # edit box

  # buttons to load rasters from file

  # second group for buttons (stacked horizontally to save space)

  cost.raster.group2 <<- ggroup(container=cb.frame, horizontal = TRUE)

  CB.A.button <<- gbutton("Load Cost Intercept Raster", container=cost.raster.group2,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata)) # if data loaded
      {
        # store data in strategy.params
        strategy.params$CB.B.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$CB.B[[curr.strategy]] <<- allthedata$the.data
        visible(CB.B.plot.raster)=TRUE  # select plot window
        svalue(strats.raster.pad)=4
        plotdata(strategy.params$CB.B[[curr.strategy]]) # plot raster
      }
    })

  CB.B.button <<- gbutton("Load Cost Slope Raster", container=cost.raster.group2,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata)) # if data loaded
      {
        # store data in strategy.params
        strategy.params$CB.A.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$CB.A[[curr.strategy]] <<- allthedata$the.data
        visible(CB.A.plot.raster)=TRUE # select plot window
        svalue(strats.raster.pad)=5
        plotdata(strategy.params$CB.A[[curr.strategy]]) # plot raster
      }
    })

  CB.C.button <<- gbutton("Load Hourly Cost Raster", container=cost.raster.group2,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata)) # if data loaded
      {
        # store data in strategy.params
        strategy.params$CB.C.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$CB.C[[curr.strategy]] <<- allthedata$the.data
        visible(CB.C.plot.raster)=TRUE # select plot window
        svalue(strats.raster.pad)=6
        plotdata(strategy.params$CB.C[[curr.strategy]]) # plot raster
      }
    })

  cost.raster.label.group = ggroup(container = strategy.params.group, horizontal = TRUE) # create space for label and checkbox
  a.label = glabel('Spatially explicit cost and target density:', container = cost.raster.label.group)
  # set a variable which defines whether we're using rasters for management or not (initialise to false)
  use.cost.rasters <<- FALSE
  # create a checkbox that defines whether the strategy uses rasters or spatially homogeneous values
  cost.raster.check <<- gcheckbox('', checked = use.cost.rasters, container = cost.raster.label.group,
    handler = cost.raster.check.handler) # handler in handlers.r

  # set visibility of raster group (initially false)
  visible(cost.raster.group1) = use.cost.rasters
  visible(cost.raster.group2) = use.cost.rasters

  # buttons to load Priority and Culling Mask rasters from file
  # (always visible as opposed to above four buttons)

  # label to show filename of loaded raster
  # TODO: do we need this? if so, should we have it for the above four rasters too?
  PR.label <<- glabel("",container=strategy.params.group)
  PR.button <- gbutton("Load Priority Raster", container=strategy.params.group,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata))  # if data loaded
      {
        # store data in strategy.params
        strategy.params$PR.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$PR[[curr.strategy]] <<- allthedata$the.data
        svalue(PR.label)=allthedata$the.file # set label to filename
        visible(PR.plot.raster)=TRUE # select plot window
        svalue(strats.raster.pad)=1
        plotdata(strategy.params$PR[[curr.strategy]]) # plot raster
      }
    })

  C.label <<- glabel("",container=strategy.params.group) # label to show filename of loaded raster
  C.button <- gbutton("Load Culling Mask Raster", container=strategy.params.group,
    handler=function(h,...){ # handler for when button pressed
      allthedata = raster.load.button.handler() # function to load user-selected data
      if (!is.null(allthedata))  # if data loaded
      {
        # store data in strategy.params
        strategy.params$C.rast[[curr.strategy]] <<- allthedata$the.rast
        strategy.params$C.mask[[curr.strategy]] <<- allthedata$the.data
        svalue(C.label)=allthedata$the.file # set label to filename
        visible(C.plot.raster)=TRUE  # select plot window
        svalue(strats.raster.pad)=2
        plotdata(strategy.params$C.mask[[curr.strategy]]) # plot raster
      }
    })

  # create notebook to display raster plots
  strats.raster.pad <<- gnotebook(container=strategy.params.group,expand=TRUE)

  # TODO: any way to limit repetition in the raster plots below?
  # Identical process, different only in that performed on different variables

  # create a tab in the notebook for the Priority raster
  PR.plot.group <- ggroup(container=strats.raster.pad,label="Priority",horizontal=FALSE)
  PR.plot.raster <<- ggraphics(container=PR.plot.group) # create plot window for tab
  PR.plot.buttons <- ggroup(container = PR.plot.group, horizontal = TRUE) # create space for buttons
  PR.plot.edit <- gbutton("Edit",container=PR.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$PR[[curr.strategy]] <<- edit(strategy.params$PR[[curr.strategy]])
      visible(PR.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=1
      strategy.params$PR.rast[[curr.strategy]] <<- raster(strategy.params$PR[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$PR[[curr.strategy]]) # plot raster
    })
  PR.plot.save <- gbutton("Save",container=PR.plot.buttons, # Save button
    handler=function(h,...){ # handler for when button pressed
      filename=fileSaveChoose("print") # choose filename to save raster file to
      if(!is.na(filename)){ # if valid filename chosen
        strategy.params$PR.rast[[curr.strategy]]=raster(strategy.params$PR[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
        writeRaster(strategy.params$PR.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
      }
    })

  # create a tab in the notebook for the Culling Mask raster
  C.plot.group <- ggroup(container=strats.raster.pad,label="Culling Mask",horizontal=FALSE)
  C.plot.raster <<- ggraphics(container=C.plot.group) # create plot window for tab
  C.plot.buttons <- ggroup(container = C.plot.group, horizontal = TRUE) # create space for buttons
  C.plot.edit <- gbutton("Edit",container=C.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$C.mask[[curr.strategy]] <<- edit(strategy.params$C.mask[[curr.strategy]])
      visible(C.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=2
      strategy.params$C.rast[[curr.strategy]]<<- raster(strategy.params$C.mask[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$C.mask[[curr.strategy]]) # plot raster
    })
  C.plot.save <- gbutton("Save",container=C.plot.buttons, # Save button
    handler=function(h,...){ # handler for when button pressed
      filename=fileSaveChoose("print") # choose filename to save raster file to
      if(!is.na(filename)){ # if valid filename chosen
        strategy.params$C.rast[[curr.strategy]] <<- raster(strategy.params$C.mask[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
        writeRaster(strategy.params$C.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
      }
    })

  # create a tab in the notebook for the Culling Mask raster
  TD.plot.group <- ggroup(container=strats.raster.pad,label="Target Density",horizontal=FALSE)
  TD.plot.raster <<- ggraphics(container=TD.plot.group) # create plot window for tab
  TD.plot.buttons <- ggroup(container = TD.plot.group, horizontal = TRUE) # create space for buttons
  TD.plot.edit <- gbutton("Edit",container=TD.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$TD[[curr.strategy]] <<- edit(strategy.params$TD[[curr.strategy]])
      visible(TD.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=3
      strategy.params$TD.rast[[curr.strategy]] <<- raster(strategy.params$TD[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$TD[[curr.strategy]]) # plot raster
    })
  TD.plot.save <- gbutton("Save",container=TD.plot.buttons, # Save button
    handler=function(h,...){ # handler for when button pressed
      filename=fileSaveChoose("print") # choose filename to save raster file to
      if(!is.na(filename)){ # if valid filename chosen
        strategy.params$TD.rast[[curr.strategy]] <<- raster(strategy.params$TD[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
        writeRaster(strategy.params$TD.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
      }
    })

  # create a tab in the notebook for the Cost Intercept raster
  CB.B.plot.group <- ggroup(container=strats.raster.pad,label="Cost Intercept",horizontal=FALSE)
  CB.B.plot.raster <<- ggraphics(container=CB.B.plot.group) # create plot window for tab
  CB.B.plot.buttons <- ggroup(container = CB.B.plot.group, horizontal = TRUE) # create space for buttons
  CB.B.plot.edit <- gbutton("Edit",container=CB.B.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$CB.B[[curr.strategy]] <<- edit(strategy.params$CB.B[[curr.strategy]])
      visible(CB.B.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=4
      strategy.params$CB.B.rast[[curr.strategy]] <<- raster(strategy.params$CB.B[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$CB.B[[curr.strategy]]) # plot raster
    })
  CB.B.plot.save <- gbutton("Save",container=CB.B.plot.buttons, # Save button
    handler=function(h,...){ # handler for when button pressed
      filename=fileSaveChoose("print") # choose filename to save raster file to
      if(!is.na(filename)){ # if valid filename chosen
        strategy.params$CB.B.rast[[curr.strategy]] <<- raster(strategy.params$CB.B[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
        writeRaster(strategy.params$CB.B.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
      }
    })

  # create a tab in the notebook for the Cost Slope raster
  CB.A.plot.group <- ggroup(container=strats.raster.pad,label="Cost Slope",horizontal=FALSE)
  CB.A.plot.raster <<- ggraphics(container=CB.A.plot.group) # create plot window for tab
  CB.A.plot.buttons <- ggroup(container = CB.A.plot.group, horizontal = TRUE) # create space for buttons
  CB.A.plot.edit <- gbutton("Edit",container=CB.A.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$CB.A[[curr.strategy]] <<- edit(strategy.params$CB.A[[curr.strategy]])
      visible(CB.A.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=5
      strategy.params$CB.A.rast[[curr.strategy]]<<-raster(strategy.params$CB.A[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$CB.A[[curr.strategy]]) # plot raster
    })
  CB.A.plot.save <- gbutton("Save",container=CB.A.plot.buttons, # Save button
    handler=function(h,...){ # handler for when button pressed
      filename=fileSaveChoose("print") # choose filename to save raster file to
      if(!is.na(filename)){ # if valid filename chosen
        strategy.params$CB.A.rast[[curr.strategy]]<<-raster(strategy.params$CB.A[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
        writeRaster(strategy.params$CB.A.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
      }
    })

  # create a tab in the notebook for the Hourly Cost raster
  CB.C.plot.group <- ggroup(container=strats.raster.pad,label="Hourly Cost",horizontal=FALSE)
  CB.C.plot.raster <<- ggraphics(container=CB.C.plot.group) # create plot window for tab
  CB.C.plot.buttons <- ggroup(container = CB.C.plot.group, horizontal = TRUE) # create space for buttons
  CB.C.plot.edit <- gbutton("Edit",container=CB.C.plot.buttons, # Edit button
    handler=function(h,...){ # handler for when button pressed
      # open R's inbuilt editor, save results
      strategy.params$CB.C[[curr.strategy]] <<- edit(strategy.params$CB.C[[curr.strategy]])
      visible(CB.C.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad)=6
      strategy.params$CB.C.rast[[curr.strategy]]<<-raster(strategy.params$CB.C[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      plotdata(strategy.params$CB.C[[curr.strategy]]) # plot raster
    })
  CB.C.plot.save <- gbutton("Save",container=CB.C.plot.buttons, # Save button
  handler=function(h,...){ # handler for when button pressed
    filename=fileSaveChoose("print") # choose filename to save raster file to
    if(!is.na(filename)){ # if valid filename chosen
      strategy.params$CB.C.rast[[curr.strategy]]<<-raster(strategy.params$CB.C[[curr.strategy]]) # save raster version (TODO: use first species carrying capacity as template)
      writeRaster(strategy.params$CB.C.rast[[curr.strategy]],filename,format="GTiff") # write raster to file in GTiff format
    }
  })

  # if the current strategy is defined, plot Priority and Culling Mask rasters
  # (the other four are dealt with by the checkbox)
  if (!is.null(curr.strategy))
  {
    if (length(strategy.params$PR.rast[[curr.strategy]]) > 0) # if opening window with carrying capacity already present
    {
      DATA = raster::as.matrix(strategy.params$PR.rast[[curr.strategy]]) # extract data from strategy.params
      visible(PR.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad) = 1
      plotdata(DATA) # plot data
    }

    if (length(strategy.params$C.rast[[curr.strategy]]) > 0) # if opening window with carrying capacity already present
    {
      DATA = raster::as.matrix(strategy.params$C.rast[[curr.strategy]])  # extract data from strategy.params
      visible(C.plot.raster)=TRUE # select plot window
      svalue(strats.raster.pad) = 2
      plotdata(DATA) # plot data
    }
  }
}

max_stars_repo_name: nickbeeton/SPADE

max_issues_repo_name: nickbeeton/SPADE

max_forks_repo_name: nickbeeton/SPADE


text:
hoge <- read.csv("synthesized_info.csv")
dead <- subset(hoge, status == "dead")
arrive <- subset(hoge, status == "arrive")
par(mfrow=c(2,1))

plot(dead$lifetime, dead$average.bug, xlim = c(0,2500), ylim = c(0,10), xlab = "lifetime [day]" ,ylab = "# of average bug per month", main = "Removed test codes")
abline(v=median(dead$lifetime), lty=2)
plot(arrive$lifetime, arrive$average.bug, xlim = c(0,2500), ylim = c(0,10), xlab = "lifetime [day]" ,ylab = "# of average bug per month", main = "Arriving test codes")
abline(v=median(arrive$lifetime), lty=2)
dev.copy2eps(file="../plot_lifetime_bug.eps")

plot(dead$commit.interval, dead$average.bug, xlim = c(0,950), ylim = c(0,10),  xlab = "commit interval [day]" ,ylab = "# of average bug per month", main = "Removed test codes")
abline(v=median(dead$commit.interval), lty=2)
plot(arrive$commit.interval, arrive$average.bug, xlim = c(0,950), ylim = c(0,10),  xlab = "commit interval [day]" ,ylab = "# of average bug per month", main = "Arriving test codes")
abline(v=median(arrive$commit.interval), lty=2)
dev.copy2eps(file="../plot_interval_bug.eps")

plot(dead$co.evolution.rate, dead$average.bug, xlim = c(0,1), ylim = c(0,10), xlab = "co-evolution rate" ,ylab = "# of average bug per month", main = "Removed test codes")
abline(v=median(dead$co.evolution.rate), lty=2)
plot(arrive$co.evolution.rate, arrive$average.bug, xlim = c(0,1), ylim = c(0,10), xlab = "co-evolution rate" ,ylab = "# of average bug per month", main = "Arriving test codes")
abline(v=median(arrive$co.evolution.rate), lty=2)
dev.copy2eps(file="../plot_co_evo_bug.eps")

> median(dead$lifetime)
[1] 469
> median(arrive$lifetime)
[1] 1289.5
> median(dead$commit.interval)
[1] 28.5
> median(arrive$commit.interval)
[1] 33.75
> median(dead$co.evolution.rate)
[1] 0.1431154
> median(arrive$co.evolution.rate)
[1] 0.05882774


text:
Every step +1 Char

$_SESSION['captcha_length']

Login Fail Count

Wrong captcha

captcha Image Must be class for refresh all together.

each captcha input must be has separated parent

max_stars_repo_name: Milad-Abooali/safego

max_issues_repo_name: Milad-Abooali/safego

max_forks_repo_name: Milad-Abooali/safego


text:
function [lia,locb] = ismember(a,b,flag1,flag2)
%ISMEMBER True for set member.
%   LIA = ISMEMBER(A,B) for cell arrays A and B returns an array of the same
%   size as A containing true where the elements of A are in B and false
%   otherwise.
%
%   [LIA,LOCB] = ISMEMBER(A,B) also returns an array LOCB containing the
%   lowest absolute index in B for each element in A which is a member of
%   B and 0 if there is no such index.
%
%   Inputs A and B must be strings or cell arrays of strings.
%
%     The behavior of ISMEMBER has changed.  This includes:
%       - occurrence of indices in LOCB switched from highest to lowest
%     If this change in behavior has adversely affected your code, you may
%     preserve the previous behavior with:
%        [LIA,LOCB] = ISMEMBER(A,B,'legacy')
%
%   Example:
%
%      a = {'red','green'}
%      b = {'gray','blue','red','orange'}
%      [lia,locb] = ismember(a,b)
%
%      % returns  lia = [1 0]  and  locb = [3 0]
%
%   See also UNIQUE, UNION, INTERSECT, SETDIFF, SETXOR, SORT, SORTROWS.

%   Copyright 1984-2014 The MathWorks, Inc.

if nargin <= 2
    % The only A and B allowed are character arrays or cellstr.
    if ~((ischar(a) || iscellstr(a)) && (ischar(b) || iscellstr(b)))
        error(message('MATLAB:ISMEMBER:InputClass',class(a),class(b)));
    end
    % Scalar A: no sort needed
    if (ischar(a) && isrow(a)) || isscalar(a)
        match = strcmp(a,b);
        lia = any(match(:));
        if nargout > 1
            if lia
                locb = find(match,1,'first');
            else
                locb = 0;
            end
        end
    % Scalar B: no sort needed
    elseif (ischar(b) && isrow(b)) || isscalar(b)
        lia = strcmp(a,b);
        if nargout > 1
            locb = double(lia);
        end
    else
        % If A or B is char, convert it to a cellstr and remove trailing spaces
        if ischar(a)
            a = cellstr(a);
        end
        if ischar(b)
            b = cellstr(b);
        end
        % Duplicates within the sets are eliminated
        [uA,~,icA] = unique(a(:),'sorted');
        if nargout <= 1
            uB = unique(b(:),'sorted');
        else
            [uB,ib] = unique(b(:),'sorted');
        end
        % Sort the unique elements of A and B, duplicate entries are adjacent
        [sortuAuB,IndSortuAuB] = sort([uA;uB]);
        % d indicates the indices matching entries
        d = strcmp(sortuAuB(1:end-1),sortuAuB(2:end));
        indReps = IndSortuAuB(d);                  % Find locations of repeats
        if nargout <= 1
            lia = ismember(icA,indReps);           % Find repeats among original list
        else
            szuA = size(uA,1);
            d = find(d);
            [lia,locb] = ismember(icA,indReps);    % Find locb by using given indices
            newd = d(locb(lia));                   % NEWD is D for non-unique A
            where = ib(IndSortuAuB(newd+1)-szuA);
            locb(lia) = where;
        end
        lia = reshape(lia,size(a));
        if nargout > 1
            locb = reshape(locb,size(a));
        end
    end
elseif nargin == 3
    % ISMEMBER(a,b,'legacy'), ISMEMBER(a,b,'R2012a')
    if nargout < 2
        lia = cellismemberlegacy(a,b,flag1);
    else
        [lia,locb] = cellismemberlegacy(a,b,flag1);
    end
else
    % ISMEMBER(a,b,'rows','legacy'), ISMEMBER(a,b,'rows','R2012a')
    if nargout < 2
        lia = cellismemberlegacy(a,b,flag1,flag2);
    else
        [lia,locb] = cellismemberlegacy(a,b,flag1,flag2);
    end
end

max_stars_repo_name: Bondzio/lego

max_issues_repo_name: Bondzio/lego

max_forks_repo_name: Bondzio/lego


text:
#' No labels in ggplot2
#'
#' @description Object, that clears ggplot2 labels
#'
#'
#' @return ggplot2 object
  #' @export
#'
#' @examples
#'
#' library(ggplot2)
#'
#' plot <-
#'   ggplot(mtcars,aes(hp,drat))+
#'   geom_point()
#'
#' plot
#'
#' plot + no_labels
#'

no_labels <-
  ggplot2::labs(
    x = "",
    y = "",
    col = "",
    fill = "",
    size = "",
    alpha = ""
    )

max_stars_repo_name: larissabf/relper

max_issues_repo_name: larissabf/relper

max_forks_repo_name: larissabf/relper


text:
# SS(Total) = Syy
SStotal=6066.1667
rxy=0.606
SSregression = (rxy^2)*SStotal
print(SSregression)

SSresidual=SStotal-SSregression
print(SSresidual)


text:
#' Simple function to extract dimensions of EQ-5D profile
#'
#' Extracts a specific one-digit number from a five-digit number.
#' @param input EQ-5D profile, usually a five-digit number.
#' @param place Number from 1 to 5 indicating which number to extract.
#' @return Number if input is a five-digit number, else NA.
#' @examples
#' eq5d_3l(input = 12312, 3)
#' # [1] 3
#'
#' eq5d_3l(input = NA, 3)
#' # [1] NA
#'
#' eq5d_3l(input = 1231, 3)
#' # [1] NA
#'
#' eq5d_3l(input = c(12312, NA, 1231, 11121), 3)
#' # [1]  3 NA NA  1
#' @export
eq5d_3l <- function(input, place) {
  if (is.factor(input)) {
    input <- as.character(input)
  }
  ret <-
    ifelse(nchar(input) != 5 | is.na(input), NA, substr(input, place, place))
  ret <- as.numeric(ret)
  return(ret)
}

#' Calculate Swedish EQ-5D-index using all five values extracted from eq5d_3l()
#'
#' Calculates Swedish EQ-ED-index using a known linear model (Model 4).
#' @param Mobility First number extracted from eq5d_3l().
#' @param Selfcare Second number extracted from eq5d_3l().
#' @param Usual Third number extracted from eq5d_3l().
#' @param Pain Fourth number extracted from eq5d_3l().
#' @param Anxiety Fifth number extracted from eq5d_3l().
#' @return Number, EQ-5D-index.
#' @examples
#' # Returns intercept
#' eqindSwed(1,1,1,1,1)
#' # [1] 0.9694
#'
#' eqindSwed(NA,1,1,1,1)
#' # [1] NA
#'
#' eqindSwed(5,1,1,1,1)
#' # [1] NA
#' @export

eqindSwed <- function(Mobility, Selfcare, Usual, Pain, Anxiety) {
  Mobility[!Mobility %in% 1:3] <- NA
  Selfcare[!Selfcare %in% 1:3] <- NA
  Usual[!Usual %in% 1:3] <- NA
  Pain[!Pain %in% 1:3] <- NA
  Anxiety[!Anxiety %in% 1:3] <- NA

  mobility <- ifelse(Mobility == 1, 0, ifelse(Mobility == 2,-0.0666,-0.1247))
  selfcare <- ifelse(Selfcare == 1, 0, -0.0276)
  usual    <- ifelse(Usual    == 1, 0, ifelse(Usual    == 2,-0.1012,-0.1355))
  pain     <- ifelse(Pain     == 1, 0, ifelse(Pain     == 2,-0.0345,-0.0904))
  anxiety  <- ifelse(Anxiety  == 1, 0, ifelse(Anxiety  == 2,-0.0552,-0.2077))
  N3       <-
    ifelse(3 %in% c(Mobility, Selfcare, Usual, Pain, Anxiety), -0.0433, 0)

  0.9694 + mobility + selfcare + usual + pain + anxiety + N3
}

#' Calculate Brittish EQ-5D-index using all five values extracted from eq5d_3l()
#'
#' Calculates Brittish EQ-ED-index using a known linear model (Dolan).
#' @param Mobility First number extracted from eq5d_3l().
#' @param Selfcare Second number extracted from eq5d_3l().
#' @param Usual Third number extracted from eq5d_3l().
#' @param Pain Fourth number extracted from eq5d_3l().
#' @param Anxiety Fifth number extracted from eq5d_3l().
#' @return Number, EQ-5D-index.
#' @examples
#' # Returns intercept
#' eqindBrit(1,1,1,1,1)
#' eqindBrit(NA,1,1,1,1)
#' eqindBrit(5,1,1,1,1)
#' @export

eqindBrit <- function(Mobility, Selfcare, Usual, Pain, Anxiety) {
  Mobility[!Mobility %in% 1:3] <- NA
  Selfcare[!Selfcare %in% 1:3] <- NA
  Usual[!Usual %in% 1:3] <- NA
  Pain[!Pain %in% 1:3] <- NA
  Anxiety[!Anxiety %in% 1:3] <- NA

  mobility <- ifelse(Mobility == 1, 0, ifelse(Mobility == 2,-0.069,-0.314))
  selfcare <- ifelse(Selfcare == 1, 0, ifelse(Selfcare == 2,-0.104,-0.214))
  usual    <- ifelse(Usual == 1,    0, ifelse(Usual    == 2,-0.036,-0.094))
  pain     <- ifelse(Pain == 1,     0, ifelse(Pain     == 2,-0.123,-0.386))
  anxiety  <- ifelse(Anxiety == 1,  0, ifelse(Anxiety  == 2,-0.071,-0.236))
  N2       <-  ifelse((
      Mobility %in% c(2, 3) |
        Selfcare %in% c(2, 3) |
        Usual    %in% c(2, 3) |
        Pain     %in% c(2, 3) |
        Anxiety  %in% c(2, 3)
      ), -0.081, 0)

  N3 <- ifelse(3 %in% c(Mobility, Selfcare, Usual, Pain, Anxiety), -0.269, 0)

  1 + mobility + selfcare + usual + pain + anxiety + N2 + N3
}

max_stars_repo_name: gforge/shprplotfun

max_issues_repo_name: gforge/shprplotfun

max_forks_repo_name: gforge/shprplotfun


text:
GetGeneList_IntPath <- function(gene = NA, label = NA, direction = 1, steps = 20) {
  # This function gets certain downstream/upstream genes in given kegg pathways

  ## === input ===
  ### gene = NA: default as return all the genes
  ### label: list of pathway ids (some problem with the space!)
  ### direction: 1 for downstream; -1 for upstream
  ### step = 20: default as all the upstream/downstream genes; if specified, iterate given steps

  ## == output ===
  ### genelist: certain downstream/upstream genes in given kegg pathways

  load("./repertoire/IntPathRepertoire.RData")

  if (!is.na(label)) {
    filter = rep(F, dim(IntPathGenePairs)[1])
    for (i in 1:length(label)) {
      filter = filter | (as.character(IntPathGenePairs$label) == label[i])
    }
    filter = filter & !(as.character(IntPathGenePairs$reaction_type) == "GPrel") # here we don't consider binding as downstream
    trimmed.file = IntPathGenePairs[filter,]
  } else {
    trimmed.file = IntPathGenePairs
  }

  if (is.na(gene)) {
    result = unique(c(as.character(trimmed.file$node1), as.character(trimmed.file$node2)))
  } else {
    # generate matrix representing linkages
    node1.list = as.character(trimmed.file$node1)
    node2.list = as.character(trimmed.file$node2)
    all_nodes = unique(c(node1.list, node2.list))
    linkage.matrix = matrix(0, nrow = length(all_nodes), ncol = length(all_nodes), dimnames = list(all_nodes, all_nodes))
    for (j in 1:dim(trimmed.file)[1]) {
      linkage.matrix[node1.list[j], node2.list[j]] = 1
    }
    # get required genes
    if (direction == -1) {
      linkage.matrix = t(linkage.matrix)
    }
    result = c()
    current_number = 0
    for (cnt in 1:steps) {
      if (length(gene) > 1) {
        gene = colnames(linkage.matrix)[colSums(linkage.matrix[gene,]) > 0]
      } else {
        gene = colnames(linkage.matrix)[linkage.matrix[gene,] == 1]
      }
      result = c(result, gene)
      result = unique(result)
      tmp_cnt = length(result)
      if(current_number == tmp_cnt) {
        break
      } else {
        current_number = tmp_cnt
      }
    }
  }

  list(gene.list = result, steps =cnt)
}

max_stars_repo_name: ShaoGroup/PICheM

max_issues_repo_name: ShaoGroup/PICheM

max_forks_repo_name: ShaoGroup/PICheM


text:
#Usage:
#Rscript cell_classified_genes.r gene_array
library("monocle")
args <- commandArgs(trailingOnly=TRUE)
#print(args)
if(length(args)==0)
{
	print("Usage: Rscript cell_classified_genes.r gene_array")
}
HSMM <- readRDS("classified_cells.rds")
selected_genes <- row.names(subset(fData(HSMM),gene_short_name %in% args[1:length(args)]))
print(selected_genes)
cds_subset <- HSMM[selected_genes,]
print(plot_genes_jitter(cds_subset, grouping='CellType',color_by="CellType"))
ggsave("cell_classified_genes_plot.png")
ggsave("cell_classified_genes_plot.tiff", width = 7, height = 7,dpi=600)
dev.off()

max_stars_repo_name: nikhil/RAS

max_issues_repo_name: nikhil/RAS

max_forks_repo_name: nikhil/RAS


text:
###################################################################
#
# Mark Cembrowski, Janelia Farm, Jan 29 2014
#
# This function takes a string of gene symbols and returns all
# associated gene IDs.
#
# INPUT:
# syms: character.  one or more gene symbols to convert.
#
# OPTIONAL ARGUMENTS:
# keepMissed=F: logical.  if supplied with symbol(s) that do not
# match any gene_short_name, keep these as 'NA' in the output?
# quiet=F: logical. suppress printout of any missed conversions.
#
###################################################################

symToId <- function(syms,keepMissed=F,quiet=F){
	# Load up vectors corresponding to gene short names and
	# gene ids.
	numSamps <- length(unique(fpkmPool$sample_name))

	numInds <- nrow(fpkmPool)/numSamps
	allGeneIds <- as.vector( (fpkmPool$gene_id)[1:numInds] )
	allGsns <- as.vector( (fpkmPool$gene_short_name)[1:numInds] )

	# Identify indices associated with supplied short names, and
	# convert to gene IDs.
	theInds <- which(allGsns%in%syms)
	theIds <- allGeneIds[theInds]
	names(theIds) <- allGsns[theInds]

	theIds <- theIds[syms] # Order according to order of input symbols.

	# Screen for any missed symbols; announce.
	missed <- syms[!(syms%in%allGsns)]

	if(length(missed)>0.1){
		if(!quiet){
			print("Note: missed the following during a sym -> id conversion:")
			print(missed)
		}

		# Remove from list, if option selected.
		if(!keepMissed){
			theIds <- theIds[-which(!(syms%in%allGsns))]
		}
	}

	theOut <- list(theIds,missed)

	invisible(theOut)
}

max_stars_repo_name: cembrowskim/hippXValidate

max_issues_repo_name: cembrowskim/hippXValidate

max_forks_repo_name: cembrowskim/hippXValidate


text:
#' @export
#'
#' @title Compute average carcass persistence
#'
#' @description This computes average carcass persistence
#' during a set of searches, given a cummulative hazard distribution
#'
#' @param searches A vector containing julian dates of all
#' carcass searches. Only the number of searches and difference
#' between min and max matter.  Average search interval is
#' computed and average persistence probability is computed
#' between 0 and average interval.  Hence, placement of searches
#' within \code{min(searches)} and \code{max(searches)} does
#' not matter.
#'
#' @param dist Character string naming the survival distribution
#' to use. Valid values are "exponential" and "weibull".
#'
#' @param shape Shape parameter of the survival distribution.
#' If \code{dist} == "exponential", this parameter is ignored.
#'
#' @param scale Scale parameter of the survival distribution.
#' If \code{dist} == "exponential", this parameter is the mean,
#' which is 1/rate, where rate is the parameter of \code{pexp} in R.
#'
#' @return A single number.  The average probability of a carcass
#' surviving from 0 to the average search interval. i.e.,
#' average survival under the survival curve from 0 to
#' \code{diff(range(searches))/2}.
#'
#' @author Trent McDonald
#'
#' @seealso \code{\link{eoar}}
#'
#' @examples
#' avgPersist(c(0,15,30),"exponential",scale=8.6251)
#' avgPersist(c(0,15,30),"weibull",0.57792,7.1683)
#'
#' @export
#'
avgPersist <- function(searches, dist="exponential", shape, scale){

 Ir <- diff(range(searches))/2

 tmp <- seq(0,Ir,length=100)
 tmp2 <- tmp[-length(tmp)] + diff(tmp)[1]/2

 if( dist=="exponential"){
    ans <- mean(1-pexp(tmp2,1/scale))
 } else if( dist == "weibull"){
    ans <- mean(1-pweibull(tmp2,shape,scale))
 } else {
    stop("Unknown carcass persistence distribution.")
 }

 ans

}