text:
#' Legend guide
#'
#' Legend type guide shows key (i.e., geoms) mapped onto values.
#' Legend guides for various scales are integrated if possible.
#'
#' Guides can be specified in each `scale_*` or in [guides()].
#' `guide = "legend"` in `scale_*` is syntactic sugar for
#' `guide = guide_legend()` (e.g. `scale_color_manual(guide = "legend")`).
#' As for how to specify the guide for each scale in more detail,
#' see [guides()].
#'
#' @param title A character string or expression indicating a title of guide.
#'   If `NULL`, the title is not shown. By default
#'   ([waiver()]), the name of the scale object or the name
#'   specified in [labs()] is used for the title.
#' @param title.position A character string indicating the position of a
#'   title. One of "top" (default for a vertical guide), "bottom", "left"
#'  (default for a horizontal guide), or "right."
#' @param title.theme A theme object for rendering the title text. Usually the
#'   object of [element_text()] is expected. By default, the theme is
#'   specified by `legend.title` in [theme()] or theme.
#' @param title.hjust A number specifying horizontal justification of the
#'   title text.
#' @param title.vjust A number specifying vertical justification of the title
#'   text.
#' @param label logical. If `TRUE` then the labels are drawn. If
#'   `FALSE` then the labels are invisible.
#' @param label.position A character string indicating the position of a
#'   label. One of "top", "bottom" (default for horizontal guide), "left", or
#'   "right" (default for vertical guide).
#' @param label.theme A theme object for rendering the label text. Usually the
#'   object of [element_text()] is expected. By default, the theme is
#'   specified by `legend.text` in [theme()].
#' @param label.hjust A numeric specifying horizontal justification of the
#'   label text.
#' @param label.vjust A numeric specifying vertical justification of the label
#'   text.
#' @param keywidth A numeric or a [grid::unit()] object specifying
#'   the width of the legend key. Default value is `legend.key.width` or
#'   `legend.key.size` in [theme()].
#' @param keyheight A numeric or a [grid::unit()] object specifying
#'   the height of the legend key. Default value is `legend.key.height` or
#'   `legend.key.size` in [theme()].
#' @param direction  A character string indicating the direction of the guide.
#'   One of "horizontal" or "vertical."
#' @param default.unit A character string indicating [grid::unit()]
#'   for `keywidth` and `keyheight`.
#' @param override.aes A list specifying aesthetic parameters of legend key.
#'   See details and examples.
#' @param nrow The desired number of rows of legends.
#' @param ncol The desired number of column of legends.
#' @param byrow logical. If `FALSE` (the default) the legend-matrix is
#'   filled by columns, otherwise the legend-matrix is filled by rows.
#' @param reverse logical. If `TRUE` the order of legends is reversed.
#' @param order positive integer less than 99 that specifies the order of
#'   this guide among multiple guides. This controls the order in which
#'   multiple guides are displayed, not the contents of the guide itself.
#'   If 0 (default), the order is determined by a secret algorithm.
#' @param ... ignored.
#' @export
#' @family guides
#' @examples
#' \donttest{
#' df <- expand.grid(X1 = 1:10, X2 = 1:10)
#' df$value <- df$X1 * df$X2
#'
#' p1 <- ggplot(df, aes(X1, X2)) + geom_tile(aes(fill = value))
#' p2 <- p1 + geom_point(aes(size = value))
#'
#' # Basic form
#' p1 + scale_fill_continuous(guide = guide_legend())
#'
#' # Control styles
#'
#' # title position
#' p1 + guides(fill = guide_legend(title = "LEFT", title.position = "left"))
#'
#' # title text styles via element_text
#' p1 + guides(fill =
#'   guide_legend(
#'     title.theme = element_text(
#'       size = 15,
#'       face = "italic",
#'       colour = "red",
#'       angle = 0
#'     )
#'   )
#' )
#'
#' # label position
#' p1 + guides(fill = guide_legend(label.position = "left", label.hjust = 1))
#'
#' # label styles
#' p1 +
#'   scale_fill_continuous(
#'     breaks = c(5, 10, 15),
#'     labels = paste("long", c(5, 10, 15)),
#'     guide = guide_legend(
#'       direction = "horizontal",
#'       title.position = "top",
#'       label.position = "bottom",
#'       label.hjust = 0.5,
#'       label.vjust = 1,
#'       label.theme = element_text(angle = 90)
#'     )
#'   )
#'
#' # Set aesthetic of legend key
#' # very low alpha value make it difficult to see legend key
#' p3 <- ggplot(mtcars, aes(vs, am, colour = factor(cyl))) +
#'   geom_jitter(alpha = 1/5, width = 0.01, height = 0.01)
#' p3
#' # override.aes overwrites the alpha
#' p3 + guides(colour = guide_legend(override.aes = list(alpha = 1)))
#'
#' # multiple row/col legends
#' df <- data.frame(x = 1:20, y = 1:20, color = letters[1:20])
#' p <- ggplot(df, aes(x, y)) +
#'   geom_point(aes(colour = color))
#' p + guides(col = guide_legend(nrow = 8))
#' p + guides(col = guide_legend(ncol = 8))
#' p + guides(col = guide_legend(nrow = 8, byrow = TRUE))
#'
#' # reversed order legend
#' p + guides(col = guide_legend(reverse = TRUE))
#' }
guide_legend <- function(# title
                         title = waiver(),
                         title.position = NULL,
                         title.theme = NULL,
                         title.hjust = NULL,
                         title.vjust = NULL,

                         # label
                         label = TRUE,
                         label.position = NULL,
                         label.theme = NULL,
                         label.hjust = NULL,
                         label.vjust = NULL,

                         # key
                         keywidth = NULL,
                         keyheight = NULL,

                         # general
                         direction = NULL,
                         default.unit = "line",
                         override.aes = list(),
                         nrow = NULL,
                         ncol = NULL,
                         byrow = FALSE,
                         reverse = FALSE,
                         order = 0,
                         ...) {

  if (!is.null(keywidth) && !is.unit(keywidth)) {
    keywidth <- unit(keywidth, default.unit)
  }
  if (!is.null(keyheight) && !is.unit(keyheight)) {
    keyheight <- unit(keyheight, default.unit)
  }

  structure(
    list(
      # title
      title = title,
      title.position = title.position,
      title.theme = title.theme,
      title.hjust = title.hjust,
      title.vjust = title.vjust,

      # label
      label = label,
      label.position = label.position,
      label.theme = label.theme,
      label.hjust = label.hjust,
      label.vjust = label.vjust,

      # size of key
      keywidth = keywidth,
      keyheight = keyheight,

      # general
      direction = direction,
      override.aes = rename_aes(override.aes),
      nrow = nrow,
      ncol = ncol,
      byrow = byrow,
      reverse = reverse,
      order = order,

      # parameter
      available_aes = c("any"),
      ...,
      name = "legend"
    ),
    class = c("guide", "legend")
  )
}

#' @export
guide_train.legend <- function(guide, scale, aesthetic = NULL) {
  breaks <- scale$get_breaks()
  if (length(breaks) == 0 || all(is.na(breaks))) {
    return()
  }

  # in the key data frame, use either the aesthetic provided as
  # argument to this function or, as a fall back, the first in the vector
  # of possible aesthetics handled by the scale
  aes_column_name <- aesthetic %||% scale$aesthetics[1]
  key <- new_data_frame(setNames(list(scale$map(breaks)), aes_column_name))
  key$.label <- scale$get_labels(breaks)

  # Drop out-of-range values for continuous scale
  # (should use scale$oob?)
  if (!scale$is_discrete()) {
    limits <- scale$get_limits()
    noob <- !is.na(breaks) & limits[1] <= breaks & breaks <= limits[2]
    key <- key[noob, , drop = FALSE]
  }

  if (guide$reverse) key <- key[nrow(key):1, ]

  guide$key <- key
  guide$hash <- with(
    guide,
    digest::digest(list(title, key$.label, direction, name))
  )
  guide
}

#' @export
guide_merge.legend <- function(guide, new_guide) {
  new_guide$key$.label <- NULL
  guide$key <- cbind(guide$key, new_guide$key)
  guide$override.aes <- c(guide$override.aes, new_guide$override.aes)
  if (any(duplicated(names(guide$override.aes)))) {
    warn("Duplicated override.aes is ignored.")
  }
  guide$override.aes <- guide$override.aes[!duplicated(names(guide$override.aes))]
  guide
}

#' @export
guide_geom.legend <- function(guide, layers, default_mapping) {
  # arrange common data for vertical and horizontal guide
  guide$geoms <- lapply(layers, function(layer) {
    matched <- matched_aes(layer, guide, default_mapping)

    # check if this layer should be included
    include <- include_layer_in_guide(layer, matched)

    if (!include) {
      return(NULL)
    }

    if (length(matched) > 0) {
      # Filter out set aesthetics that can't be applied to the legend
      n <- vapply(layer$aes_params, length, integer(1))
      params <- layer$aes_params[n == 1]

      aesthetics <- layer$mapping
      modifiers <- aesthetics[is_scaled_aes(aesthetics) | is_staged_aes(aesthetics)]

      data <- tryCatch(
        layer$geom$use_defaults(guide$key[matched], params, modifiers),
        error = function(...) {
          warn("Failed to apply `after_scale()` modifications to legend")
          layer$geom$use_defaults(guide$key[matched], params, list())
        }
      )
    } else {
      data <- layer$geom$use_defaults(NULL, layer$aes_params)[rep(1, nrow(guide$key)), ]
    }

    # override.aes in guide_legend manually changes the geom
    data <- modify_list(data, guide$override.aes)

    list(
      draw_key = layer$geom$draw_key,
      data = data,
      params = c(layer$geom_params, layer$stat_params)
    )
  })

  # remove null geom
  guide$geoms <- compact(guide$geoms)

  # Finally, remove this guide if no layer is drawn
  if (length(guide$geoms) == 0) guide <- NULL
  guide
}

#' @export
guide_gengrob.legend <- function(guide, theme) {

  # default setting
  label.position <- guide$label.position %||% "right"
  if (!label.position %in% c("top", "bottom", "left", "right"))
    abort(glue("label position `{label.position}` is invalid"))

  nbreak <- nrow(guide$key)

  # obtain the theme for the legend title. We need this both for the title grob
  # and to obtain the title fontsize.
  title.theme <- guide$title.theme %||% calc_element("legend.title", theme)

  title.hjust <- guide$title.hjust %||% theme$legend.title.align %||% title.theme$hjust %||% 0
  title.vjust <- guide$title.vjust %||% title.theme$vjust %||% 0.5

  grob.title <- ggname("guide.title",
    element_grob(
      title.theme,
      label = guide$title,
      hjust = title.hjust,
      vjust = title.vjust,
      margin_x = TRUE,
      margin_y = TRUE
    )
  )

  title_width <- width_cm(grob.title)
  title_height <- height_cm(grob.title)
  title_fontsize <- title.theme$size %||% calc_element("legend.title", theme)$size %||%
    calc_element("text", theme)$size %||% 11

  # gap between keys etc
  # the default horizontal and vertical gap need to be the same to avoid strange
  # effects for certain guide layouts
  hgap <- width_cm(theme$legend.spacing.x  %||% (0.5 * unit(title_fontsize, "pt")))
  vgap <- height_cm(theme$legend.spacing.y %||% (0.5 * unit(title_fontsize, "pt")))

  # Labels

  # first get the label theme, we need it below even when there are no labels
  label.theme <- guide$label.theme %||% calc_element("legend.text", theme)

  if (!guide$label || is.null(guide$key$.label)) {
    grob.labels <- rep(list(zeroGrob()), nrow(guide$key))
  } else {
    # get the defaults for label justification. The defaults are complicated and depend
    # on the direction of the legend and on label placement
    just_defaults <- label_just_defaults.legend(guide$direction, label.position)
    # don't set expressions left-justified
    if (just_defaults$hjust == 0 && any(is.expression(guide$key$.label))) just_defaults$hjust <- 1

    # We break inheritance for hjust and vjust, because that's more intuitive here; it still allows manual
    # setting of hjust and vjust if desired. The alternative is to ignore hjust and vjust altogether, which
    # seems worse
    if (is.null(guide$label.theme$hjust) && is.null(theme$legend.text$hjust)) label.theme$hjust <- NULL
    if (is.null(guide$label.theme$vjust) && is.null(theme$legend.text$vjust)) label.theme$vjust <- NULL

    # label.theme in param of guide_legend() > theme$legend.text.align > default
    hjust <- guide$label.hjust %||% theme$legend.text.align %||% label.theme$hjust %||%
      just_defaults$hjust
    vjust <- guide$label.vjust %||% label.theme$vjust %||%
      just_defaults$vjust

    grob.labels <- lapply(guide$key$.label, function(label, ...) {
      g <- element_grob(
        element = label.theme,
        label = label,
        hjust = hjust,
        vjust = vjust,
        margin_x = TRUE,
        margin_y = TRUE
      )
      ggname("guide.label", g)
    })
  }

  label_widths <- width_cm(grob.labels)
  label_heights <- height_cm(grob.labels)

  # Keys
  key_width <- width_cm(
    guide$keywidth %||% theme$legend.key.width %||% theme$legend.key.size
  )
  key_height <- height_cm(
    guide$keyheight %||% theme$legend.key.height %||% theme$legend.key.size
  )

  key_size_mat <- do.call("cbind", lapply(guide$geoms, function(g) g$data$size / 10))
  if (nrow(key_size_mat) == 0 || ncol(key_size_mat) == 0) {
    key_size_mat <- matrix(0, ncol = 1, nrow = nbreak)
  }
  key_sizes <- apply(key_size_mat, 1, max)

  if (!is.null(guide$nrow) && !is.null(guide$ncol) &&
      guide$nrow * guide$ncol < nbreak) {
    abort("`nrow` * `ncol` needs to be larger than the number of breaks")
  }

  # If neither nrow/ncol specified, guess with "reasonable" values
  if (is.null(guide$nrow) && is.null(guide$ncol)) {
    if (guide$direction == "horizontal") {
      guide$nrow <- ceiling(nbreak / 5)
    } else {
      guide$ncol <- ceiling(nbreak / 20)
    }
  }
  legend.nrow <- guide$nrow %||% ceiling(nbreak / guide$ncol)
  legend.ncol <- guide$ncol %||% ceiling(nbreak / guide$nrow)

  key_sizes <- matrix(
    c(key_sizes, rep(0, legend.nrow * legend.ncol - nbreak)),
    legend.nrow,
    legend.ncol,
    byrow = guide$byrow
  )

  key_widths <- pmax(key_width, apply(key_sizes, 2, max))
  key_heights <- pmax(key_height, apply(key_sizes, 1, max))

  label_widths <- apply(
    matrix(
      c(label_widths, rep(0, legend.nrow * legend.ncol - nbreak)),
      legend.nrow,
      legend.ncol,
      byrow = guide$byrow
    ),
    2,
    max
  )
  label_heights <- apply(
    matrix(
      c(label_heights, rep(0, legend.nrow * legend.ncol - nbreak)),
      legend.nrow,
      legend.ncol,
      byrow = guide$byrow
    ),
    1,
    max
  )

  if (guide$byrow) {
    vps <- new_data_frame(list(
      R = ceiling(seq(nbreak) / legend.ncol),
      C = (seq(nbreak) - 1) %% legend.ncol + 1
    ))
  } else {
    vps <- mat_2_df(arrayInd(seq(nbreak), dim(key_sizes)), c("R", "C"))
  }

  # layout of key-label depends on the direction of the guide
  if (guide$byrow == TRUE) {
    switch(
      label.position,
      "top" = {
        kl_widths <- pmax(label_widths, key_widths)
        kl_heights <- utils::head(
          interleave(label_heights, vgap, key_heights, vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 4 - 1,
          key.col = C,
          label.row = R * 4 - 3,
          label.col = C
        )
      },
      "bottom" = {
        kl_widths <- pmax(label_widths, key_widths)
        kl_heights <- utils::head(
          interleave(key_heights, vgap, label_heights, vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 4 - 3,
          key.col = C,
          label.row = R * 4 - 1,
          label.col = C
        )
      },
      "left" = {
        kl_widths <- utils::head(
          interleave(label_widths, hgap, key_widths, hgap),
          -1
        )
        kl_heights <- utils::head(
          interleave(pmax(label_heights, key_heights), vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 2 - 1,
          key.col = C * 4 - 1,
          label.row = R * 2 - 1,
          label.col = C * 4 - 3
        )
      },
      "right" = {
        kl_widths <- utils::head(
          interleave(key_widths, hgap, label_widths, hgap),
          -1
        )
        kl_heights <- utils::head(
          interleave(pmax(label_heights, key_heights), vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 2 - 1,
          key.col = C * 4 - 3,
          label.row = R * 2 - 1,
          label.col = C * 4 - 1
        )
      })
  } else {
    switch(
      label.position,
      "top" = {
        kl_widths <- utils::head(
          interleave(pmax(label_widths, key_widths), hgap),
          -1
        )
        kl_heights <- utils::head(
          interleave(label_heights, vgap, key_heights, vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 4 - 1,
          key.col = C * 2 - 1,
          label.row = R * 4 - 3,
          label.col = C * 2 - 1
        )
      },
      "bottom" = {
        kl_widths <- utils::head(
          interleave(pmax(label_widths, key_widths), hgap),
          -1
        )
        kl_heights <- utils::head(
          interleave(key_heights, vgap, label_heights, vgap),
          -1
        )
        vps <- transform(
          vps,
          key.row = R * 4 - 3,
          key.col = C * 2 - 1,
          label.row = R * 4 - 1,
          label.col = C * 2 - 1
        )
      },
      "left" = {
        kl_widths <- utils::head(
          interleave(label_widths, hgap, key_widths, hgap),
          -1
        )
        kl_heights <- pmax(key_heights, label_heights)
        vps <- transform(
          vps,
          key.row = R,
          key.col = C * 4 - 1,
          label.row = R,
          label.col = C * 4 - 3
        )
      },
      "right" = {
        kl_widths <- utils::head(
          interleave(key_widths, hgap, label_widths, hgap),
          -1
        )
        kl_heights <- pmax(key_heights, label_heights)
        vps <- transform(
          vps,
          key.row = R,
          key.col = C * 4 - 3,
          label.row = R,
          label.col = C * 4 - 1
        )
      })
  }

  # layout the title over key-label
  switch(guide$title.position,
         "top" = {
           widths <- c(kl_widths, max(0, title_width - sum(kl_widths)))
           heights <- c(title_height, vgap, kl_heights)
           vps <- transform(
             vps,
             key.row = key.row + 2,
             key.col = key.col,
             label.row = label.row + 2,
             label.col = label.col
           )
           vps.title.row = 1; vps.title.col = 1:length(widths)
         },
         "bottom" = {
           widths <- c(kl_widths, max(0, title_width - sum(kl_widths)))
           heights <- c(kl_heights, vgap, title_height)
           vps <- transform(
             vps,
             key.row = key.row,
             key.col = key.col,
             label.row = label.row,
             label.col = label.col
           )
           vps.title.row = length(heights); vps.title.col = 1:length(widths)
         },
         "left" = {
           widths <- c(title_width, hgap, kl_widths)
           heights <- c(kl_heights, max(0, title_height - sum(kl_heights)))
           vps <- transform(
             vps,
             key.row = key.row,
             key.col = key.col + 2,
             label.row = label.row,
             label.col = label.col + 2
           )
           vps.title.row = 1:length(heights); vps.title.col = 1
         },
         "right" = {
           widths <- c(kl_widths, hgap, title_width)
           heights <- c(kl_heights, max(0, title_height - sum(kl_heights)))
           vps <- transform(
             vps,
             key.row = key.row,
             key.col = key.col,
             label.row = label.row,
             label.col = label.col
           )
           vps.title.row = 1:length(heights); vps.title.col = length(widths)
         })

  # grob for key
  key_size <- c(key_width, key_height) * 10

  draw_key <- function(i) {
    bg <- element_render(theme, "legend.key")
    keys <- lapply(guide$geoms, function(g) {
      g$draw_key(g$data[i, ], g$params, key_size)
    })
    c(list(bg), keys)
  }
  grob.keys <- unlist(lapply(seq_len(nbreak), draw_key), recursive = FALSE)

  # background
  grob.background <- element_render(theme, "legend.background")

  ngeom <- length(guide$geoms) + 1
  kcols <- rep(vps$key.col, each = ngeom)
  krows <- rep(vps$key.row, each = ngeom)

  # padding
  padding <- convertUnit(theme$legend.margin %||% margin(), "cm", valueOnly = TRUE)
  widths <- c(padding[4], widths, padding[2])
  heights <- c(padding[1], heights, padding[3])

  # Create the gtable for the legend
  gt <- gtable(widths = unit(widths, "cm"), heights = unit(heights, "cm"))
  gt <- gtable_add_grob(
    gt,
    grob.background,
    name = "background",
    clip = "off",
    t = 1,
    r = -1,
    b = -1,
    l = 1
  )
  gt <- gtable_add_grob(
    gt,
    justify_grobs(
      grob.title,
      hjust = title.hjust,
      vjust = title.vjust,
      int_angle = title.theme$angle,
      debug = title.theme$debug
    ),
    name = "title",
    clip = "off",
    t = 1 + min(vps.title.row),
    r = 1 + max(vps.title.col),
    b = 1 + max(vps.title.row),
    l = 1 + min(vps.title.col)
  )
  gt <- gtable_add_grob(
    gt,
    grob.keys,
    name = paste("key", krows, kcols, c("bg", seq(ngeom - 1)), sep = "-"),
    clip = "off",
    t = 1 + krows,
    r = 1 + kcols,
    b = 1 + krows,
    l = 1 + kcols
  )
  gt <- gtable_add_grob(
    gt,
    justify_grobs(
      grob.labels,
      hjust = hjust,
      vjust = vjust,
      int_angle = label.theme$angle,
      debug = label.theme$debug
    ),
    name = paste("label", vps$label.row, vps$label.col, sep = "-"),
    clip = "off",
    t = 1 + vps$label.row,
    r = 1 + vps$label.col,
    b = 1 + vps$label.row,
    l = 1 + vps$label.col
  )
  gt
}

#' Calculate the default hjust and vjust settings depending on legend
#' direction and position.
#'
#' @noRd
label_just_defaults.legend <- function(direction, position) {
  if (direction == "horizontal") {
    switch(
      position,
      "top" = list(hjust = 0.5, vjust = 0),
      "bottom" = list(hjust = 0.5, vjust = 1),
      "left" = list(hjust = 1, vjust = 0.5),
      list(hjust = 0, vjust = 0.5)
    )
  }
  else {
    switch(
      position,
      "top" = list(hjust = 0.5, vjust = 0),
      "bottom" = list(hjust = 0.5, vjust = 1),
      "left" = list(hjust = 1, vjust = 0.5),
      list(hjust = 0, vjust = 0.5)
    )

  }

}

utils::globalVariables(c("C", "R", "key.row", "key.col", "label.row", "label.col"))


text:
##' Provides access to arguments of nested functions. Sort of an alterative mechanism to passing `...` arguments but with more features.
##'
##' Provides access to higher level call's arguments (including `...` dots arguments) without explicitly passing it through calling stack and allows updating default values that are explicitly set throughout calling stack (i.e., lower calls take prevalence).
##'
##' @param function_or_arg_list The end function that meant to accept dots arguments (default arguments accessed with `formals(function_or_arg_list)`) or just explicit list of default dots arguments that will be searched up in calling stack and updated if set explicitly in higher calls. If set to NULL then use formals of the parent call (assessed with `sys.function(-1L)`).
##' @param select_args Which arguments to select from `function_or_arg_list`. Ether character or numeric vector.
##' @param search_calls_with_formals Formals (parameters, arguments) that should be present in each upper call to continue looking up the call stack for updates in dots arguments.
##' @param search_calls_of_env Environment/package name (character string) to which each function in upper calls to should belong to continue looking up the call stack for updates in dots arguments.
##' @param search_calls_regexp Regular expression that each function name in upper calls to should matched to continue looking up the call stack for updates in dots arguments.
##' @param search_depth Number of frames (aka environments) down in calling stack to look up arguments.
##' @param search_up_to_call The name of the call before which to continue looking up the call stack for updates in dots arguments.
##' @param skip_checks_for_parent_call Whether to skip checking `search_calls_with_formals` `search_calls_of_env` `search_calls_regexp`
##' @param eval_default_args Whether to evaluate default arguments. Default is do not evaluate (FALSE) assuming that all argument are simple values (i.e., evaluates to itself)
##' @param return_unlisted_if_single_arg Toggle wether unlist when returning a single argument. Default is TRUE
##' @examples
##' # Make get_dots available for following examples
##' get_dots <- nstandr:::get_dots
##' # Basic usage
##' util <- function(foo = 0, bar = 0) {
##'     # get dots and bind updated arguments into environment
##'     dots <- get_dots()
##'     for (v in names(dots)) assign(v, dots[[v]])
##'     # util just reports it arguments
##'     message("foo: ", foo, ", bar: ", bar)
##' }
##'
##' util()
##' #> foo: 0, bar: 0
##'
##' main <- function (...) {
##'     util()
##'     util(foo = 1)
##'     util(bar = 1)
##' }
##'
##' main(foo = 2, bar = 2)
##' #> foo: 2, bar: 2
##' #> foo: 1, bar: 2  # THIS WORKS NOW!
##' #> foo: 2, bar: 1  # THIS WORKS NOW!
##'
##' # Usage in nested calls
##' util <- function(foo = 0, bar = 0) {
##'     # get dots and bind updated arguments into environment
##'     dots <- get_dots(search_depth = 3L)
##'     for (v in names(dots)) assign(v, dots[[v]])
##'     # util just reports it arguments
##'     message("foo: ", foo, ", bar: ", bar)
##' }
##'
##' main <- function (...) {
##'     util()
##'     sub_main(foo = 1)
##' }
##'
##' sub_main <- function (...) {
##'     util()
##'     sub_sub_main(bar = 2)
##' }
##'
##' sub_sub_main <- function (...) {
##'     util()
##' }
##'
##' main()
##' #> foo: 0, bar: 0
##' #> foo: 1, bar: 0
##' #> foo: 0, bar: 2
##'
##' @return List of updated dots arguments
##'
##' @md
get_dots <- function(function_or_arg_list = NULL
                   , select_args = NULL
                   , search_calls_with_formals = "..."
                   , search_calls_of_env = NULL
                   , search_calls_regexp = NULL
                   , search_depth = 1L
                   , search_up_to_call = NULL
                   , skip_checks_for_parent_call = TRUE
                   , eval_default_args = FALSE
                   , return_unlisted_if_single_arg = TRUE) {

    ## check arguments with checkmate (optionally)
    if (requireNamespace("checkmate", quietly = TRUE)) {
        checkmate::assert(checkmate::check_function(function_or_arg_list, null.ok = TRUE)
                        , checkmate::check_list(function_or_arg_list))
        checkmate::assert_character(select_args, null.ok = TRUE)
        checkmate::assert_character(search_calls_with_formals, null.ok = TRUE)
        checkmate::assert_character(search_calls_of_env, null.ok = TRUE)
        checkmate::assert_character(search_calls_regexp, null.ok = TRUE)
        checkmate::assert_integer(search_depth)
        checkmate::assert_character(search_up_to_call, null.ok = TRUE)
        checkmate::assert_flag(eval_default_args)
        checkmate::assert_flag(return_unlisted_if_single_arg)
    }
    ## set default_args
    if (is.null(function_or_arg_list)) {
        function_or_arg_list <- sys.function(-1L)
    }
    if (is.null(function_or_arg_list)) {
        stop("get_dots -- supposed to be called inside function (nframe >= 1)")
    }
    if (is.function(function_or_arg_list)) {
        default_args <- formals(function_or_arg_list) |> as.list()
    } else if(is.list(function_or_arg_list)) {
        default_args <- function_or_arg_list
    }
    ## filter ... argument out (just in case)
    default_args <- default_args[!(default_args %in% "...")]
    if (length(select_args) > 0) {
        default_args <- default_args[select_args]
        if (length(default_args) == 0) stop("get_dots -- 'select_args' are not in 'formals(function_or_arg_list)'")
    }
    ## eval default args just in case
    if(eval_default_args &&
       length(default_args) != 0 &&
       is.function(function_or_arg_list)) {
        default_args <-
            lapply(default_args, eval, environment(function_or_arg_list))
    }
    ## collect explicit args in parents
    explicit_args <- list()
    sp <- sys.parent()
    for (fr in unique(rev(sys.parents()))) {
        ## stop searching frames stack deeper than search_depth
        if (fr < 1 || (sp - fr) > search_depth) break()
        ## check if we are searching only in 'friendly' functions:
        ## meaning that at least search_calls_with_formals should exist in calls
        parent_fun <- sys.function(fr)
        parent_default_args <- formals(parent_fun) |> as.list()
        if (!(skip_checks_for_parent_call && fr == sp) &&
            !all(search_calls_with_formals %in% names(parent_default_args))) next()
        ## check if call belongs to an env (package)
        if (!(skip_checks_for_parent_call && fr == sp) &&
            !is.null(search_calls_of_env) &&
            !(environmentName(environment(parent_fun)) %in% search_calls_of_env)) next()
        ## check if call matches regexp
        parent_call <- sys.call(fr) |> as.list()
        if (!(skip_checks_for_parent_call && fr == sp) &&
            !is.null(search_calls_regexp) &&
            !grepl(search_calls_regexp, as.character(parent_call[[1]]), perl = TRUE)) next()
        ## update defautls other parent defaults
        if(any(default_args_update <-
                   names(parent_default_args) %in% names(default_args))) {
            default_args_update <- parent_default_args[default_args_update]
            if(eval_default_args) {
                default_args_update <- lapply(default_args_update, eval, environment(parent_fun))
            }
            default_args <-
                c(default_args[!(names(default_args) %in% names(parent_default_args))]
                ,  default_args_update)
        }
        ## if explicit arg is in args list and not already added then add it
        parent_args <- parent_call[-1]
        if (length(parent_args) > 0) {
            args_to_add <-
                (names(parent_args) %in% names(default_args)) &
                !(names(parent_args) %in% names(explicit_args))
            if (any(args_to_add)) {
                ## evaluate args in the frame that was created before 'fr'
                ## the direct parent of 'fr' might be wrong environment to look for args
                args_to_add <-
                    lapply(parent_args[args_to_add]
                         , eval
                         , envir = sys.parents()[which(sys.parents() == fr)[1] - 1])
                explicit_args <- c(explicit_args, args_to_add)
            }
        }
        ## stop searching frames stack at search_up_to_call call
        if (parent_call[1] %in% search_up_to_call) break()
    }
    ## merge default and explicit args
    args_updated <- default_args
    if (length(explicit_args) != 0) {
        args_updated <-
            c(explicit_args
            , default_args[!(names(default_args) %in% names(explicit_args))])
    }
    ## just return argument if it is single argument
    if(return_unlisted_if_single_arg && length(args_updated) == 1) {
        args_updated <- args_updated[[1]]
    }
    return (args_updated)
}

ext: r

max_stars_repo_name: stasvlasov/nstandr

max_issues_repo_name: stasvlasov/nstandr

max_forks_repo_name: stasvlasov/nstandr


text:
  observer.get.counts.by.class = function( p, x ) {
       # sex codes
    male = 0
    female = 1
    sex.unknown = 2

    # maturity codes
    immature = 0
    mature = 1
    mat.unknown = 2

    #good.data = which(  x$cw>=95 & x$sex==male & ( x$shell>=2 | x$durometer >= 68 ) &  x$prodcd_id==0 & is.finite(x$mat) & is.finite(x$cw)  )
    good.data = which(  x$cw>=95 & x$sex==male & ( x$shell>=2 | x$durometer >= 68 ) &  x$prodcd_id==0 & is.finite(x$mat) & is.finite(x$cw)  )

    x = x[ good.data ,]
    x$dummy =1

    classnames = p$nodes
    classes = p$nodelabels

    nclasses = length(classes)
    x$class = NA
    for (c in 1:nclasses) x$class[ filter.class(x, classes[c]) ] = classnames[c]

    s = tapply2( x=x, indices=c("fishyr", "class"), var="dummy", newvars=c("count"), func=function(x){ sum(x) } )

    names(s) = c("yr", "vars", "n")
    s = factor2number(s, c("yr", "n" ) )
    return(s)
  }


text:
#' Symmetrises a scale
#'
#' also works in a faceted grid
#' from: http://stackoverflow.com/q/37635016
#'
#' @param p     A ggplot2 object
#' @param axis  Which axis to symmetrise, either "x" or "y"
symmetrise_scale = function(p, axis = "x"){
    gb = ggplot_build(p)
    type = switch(axis, "x" = "x.range", "y" = "y.range")
    lims = sapply(gb$layout$panel_ranges, "[[", type)
    fname = as.character(p$facet$params$facets)
    facets = gb$layout$panel_layout[ ,fname]
    lims2 = as.vector(t(tcrossprod(apply(abs(lims), 2, max), c(-1,1))))
    dummy = setNames(data.frame(rep(facets, each = 2), lims2), c(fname, axis))
    switch(axis,
           "x" = p + geom_blank(data=dummy, aes(x=x, y=Inf)),
           "y" = p + geom_blank(data=dummy, aes(x=Inf, y=y)))
}

max_stars_repo_name: mschubert/ebits

max_issues_repo_name: mschubert/ebits

max_forks_repo_name: mschubert/ebits


text:
source('topics.r')

library("argparser")

p <- arg_parser("Find the best fit of a topic model thing")

p <- add_argument(p, "--plot", help="create a plot", flag=TRUE)
p <- add_argument(p, "--method", help="choose method in use", default = "logll")
p <- add_argument(p, "--dtm", help="dtm used" )
p <- add_argument(p, "folder", help="folders to analyse")

args <- parse_args( p )

value_method <- check_fitness_ll; ## default to logll

if( args$method == 'perplexity' ) {

  library('topicmodels')

  ## there should be additional parameter for dtm
  load( '/Users/mnelimar/OneDrive - Aalto-yliopisto/projects/2016-topicmodelkritiikki/dtm_perdocument.rdata' )
  value_method <- function( x ) { perplexity( x, dtm3 ); } ## make similar to logkelihood

}

##for( path in commandArgs(trailingOnly=TRUE) ) {

  df = data.frame( k = integer(), value =integer() )

  path <- args$folder;

  for( f in list.files( path, pattern = 'topic*') ){
        print( f )
  	load( paste(path, f, sep = '') )
  	k <- model@k
    value <- value_method( model )
  	row = c(k, value)
  	df[ nrow(df)+1,] <- row
  }

  print("Examinging")
  print( path )

  print("Best fit k" )
  print( df$k[ which.max( df$value ) ] )

  print("Best fit value")
  print( df$value[ which.max( df$value ) ] )
  print("") ## Empty line

  if( args$plot ) {

    library('ggplot2')

    g <- ggplot( df , aes(k , value ) ) +
    geom_line() +
    xlab('k') + ylab('measure') +
    theme_minimal() +
    xlim(2,250) ## TODO: make args

    ggsave( file = 'plot.pdf' , g)

  }

##}

max_stars_repo_name: HIIT/hybra-core

max_issues_repo_name: HIIT/hybra-core

max_forks_repo_name: HIIT/hybra-core


text:
#' Linear simple pipeline for turning a matrix to a MC and plotting std figs
#'
#' The simplification here is that all objects (mat, gstat, marker gset, cgraph, coclust, mc, mc2d) will have the same ID and will be produced linerly without much questions asked. Also the parameters are simplified for typical use in small datasets.
#'
#'	@param mat_id the id to start from
#' @param T_vm threshold on normalized varmin, for selecting markers. Lower this is you don't get enough info from the selected markers. (0.2)
#' @param T_tot threshold on minimal total expression for markers (200)
#' @param T_top3 threshold on minimal number of genes in the top3 expressing cells for markers (3)
#' @param Knn the raw Knn parameter used for constructing a balanced Knn graph. This should be set to roughly the size of desired metacells (120)
#' @param n_resamp number of bootstraps (500). Decrease to save time.
#' @param T_weight number of co-clustering occurences in bootstrap to add an edge in the phase of generating metacells from the bootstrap co-clustering matrix.
#'
#' @export

mcell_pipe_mat2mc2d = function(mat_id,
	T_vm=0.2,
	T_tot=200,
	T_top3=3,
	Knn=120,
	n_resamp=500,
	T_weight=round(n_resamp*0.75/8))
{
	id = mat_id
	id_plot_marks = sprintf("%s.plot_marks", id)

	mcell_add_gene_stat(id, id)

	mcell_gset_filter_varmean(id, id, T_vm=T_vm, force_new=T)
	mcell_gset_filter_cov(id, id, T_tot=T_tot, T_top3=T_top3)

	mcell_add_cgraph_from_mat_bknn(id, id, id, K=Knn)

	mcell_coclust_from_graph_resamp(id, id,
	                                min_mc_size=round(Knn/6),
	                                p_resamp=0.75,
	                                n_resamp=n_resamp)

	mcell_mc_from_coclust_louv_sub(id, id, id,
			 max_clust_size=4000,
			 max_mc_size=100,
			 min_mc_size=round(Knn/6),
			 T_weight = T_weight)

	mcell_gset_from_mc_markers(gset_id=id_plot_marks, mc_id=id)

	if(!is.null(mark_fn)) {
	  marks_colors = read.table(mark_fn, sep="\t", h=T, stringsAsFactors=F)
	  mc_colorize(mc_id=id, marker_colors=marks_colors)
	}

	mcell_mc_plot_marks(mc_id=id,
						gset_id=id,
						mat_id=id)

	mcell_mc2d_force_knn(id, id, id)
	mcell_mc2d_plot(mc2d_id=id)
}

max_stars_repo_name: echomsky/metacell


text:
#!/usr/bin/env R

#!R
Args <- commandArgs()
NullModel=Args[6]
outputname=Args[7]
boxcolor=Args[8]
linecolor=Args[9]
outputname2=Args[10]
Window_Size=Args[11]

outputFileType=strsplit(outputname,'[.]')[[1]][length(strsplit(outputname,'[.]')[[1]])]
dataNull=read.table(file=NullModel,header=F)
dataNull=dataNull[order(dataNull[,1],decreasing=F),]
if(sum(dataNull[,2])>2^18-1){
    shrink=as.integer(sum(dataNull[,2])/2^18-1)+3
    dataNull[,2]=dataNull[,2]/shrink
    }
#MedianData=median(dataNull[,1])
#for(i in nrow(dataNull):1){
#    if(dataNull[i,1]>3*MedianData){
#        dataNull=dataNull[-i,]
#    }
#    }
StatType=strsplit(strsplit(NullModel,'/')[[1]][length(strsplit(NullModel,'/')[[1]])],'[.]')[[1]][1]

if (StatType=='RDNull'){
    data_RD=c()
    for(i in 1:nrow(dataNull)){    data_RD=c(data_RD,rep(dataNull[i,1],dataNull[i,2]))}
    RD_Mean=mean(data_RD)
    RD_Median=median(data_RD)
    RD_STD=sd(data_RD)
    total_length=length(data_RD)
    data_RD=data_RD[as.integer(total_length/40):as.integer(total_length-total_length/40)]
    p=RD_Mean/RD_STD^2
    r=RD_Mean*p/(1-p)
    xrange=range(min(data_RD),max(data_RD))
    yrange=range(0,max(dataNull[,2]/sum(dataNull[,2])))
    if(outputFileType=='pdf'){    pdf(outputname)}
    if(outputFileType=='jpg'){    jpeg(outputname)}
    if(outputFileType=='png'){    png(outputname)}
    par(mfrow=c(1,2))
    par(mar=rep(3,4))
    par(fig=c(0.1,0.6,0,1))
    cex_main=0.9
    break_num=as.integer(nrow(dataNull)/2)
    if (nrow(dataNull)>500){        break_num=as.integer(nrow(dataNull)/5)}
    #plot(xrange,yrange,xlab='Read Depth',ylab='Frequency',cex.main=0.8,main='Fitted by NB',type='n',xaxt='n',yaxt='n',frame.plot=F)
    data_RD3=data_RD[data_RD<3*RD_Mean]
    hist(data_RD3,xaxt='n',xlab='Read Depth',ylab='',freq=F,cex.main=cex_main,main='Fitted by NB',col=boxcolor,breaks=break_num)
    lines(seq(range(data_RD3)[1],range(data_RD3)[2],1),dnbinom(seq(range(data_RD3)[1],range(data_RD3)[2],1),r,p),type='l',lwd=2,col=linecolor)
    axis(1,c(xrange[1],xrange[1]+(xrange[2]-xrange[1])/4,xrange[1]+(xrange[2]-xrange[1])/2,xrange[1]+3*(xrange[2]-xrange[1])/4,xrange[2]),round(c(xrange[1],xrange[1]+(xrange[2]-xrange[1])/4,xrange[1]+(xrange[2]-xrange[1])/2,xrange[1]+3*(xrange[2]-xrange[1])/4,xrange[2])/as.double(Window_Size),digits=2))

    par(fig=c(0,0.15,0,1),new=T)
    par(mar=rep(1,4))
    plot(c(0,1),c(0,1),type='n',xlab='',ylab='',axes=FALSE)
    text(0.5,0.5,'Frequency',srt=90)

    par(fig=c(0.5,1,0,1),new=T)
    par(mar=rep(3,4))
    hist( data_RD3,xaxt='n',freq=F,col=boxcolor,breaks=break_num,xlab='Read Depth',ylab='',cex.main=0.8,main='Fitted by Normal')
    lines(seq(range(data_RD3)[1],range(data_RD3)[2],1),dnorm(seq(range(data_RD3)[1],range(data_RD3)[2],1),RD_Mean,RD_STD),type='l',lwd=2,col=linecolor)
    axis(1,c(xrange[1],xrange[1]+(xrange[2]-xrange[1])/4,xrange[1]+(xrange[2]-xrange[1])/2,xrange[1]+3*(xrange[2]-xrange[1])/4,xrange[2]),round(c(xrange[1],xrange[1]+(xrange[2]-xrange[1])/4,xrange[1]+(xrange[2]-xrange[1])/2,xrange[1]+3*(xrange[2]-xrange[1])/4,xrange[2])/as.double(Window_Size),digits=2))
    #axis(2,c(yrange[1],yrange[1]+(yrange[2]-yrange[1])/3,yrange[1]+2*(yrange[2]-yrange[1])/3,yrange[2]),round(c(yrange[1],yrange[1]+(yrange[2]-yrange[1])/3,yrange[1]+2*(yrange[2]-yrange[1])/3,yrange[2]),digits=2))
    dev.off()
    StatMatrix=data.frame(Mean=RD_Mean,Median=RD_Median,STD=RD_STD)
    write.table(file=outputname2,StatMatrix,quote=F,row.name=F)
    }

if (StatType=='ILNull'){
    library(mixtools)
    data_IL=c()
    for(i in 1:nrow(dataNull)){
    data_IL=c(data_IL,rep(dataNull[i,1],dataNull[i,2]))
    }
    total_length=length(data_IL)
    data_IL2=data_IL[as.integer(total_length/40): as.integer(total_length-total_length/40)]
    IL_Mean=mean(data_IL2)
    IL_Median=median(data_IL2)
    IL_STD=sd(data_IL2)
    mixmdl = normalmixEM(data_IL2)
    IL_Normal_Mean=IL_Mean
    IL_Normal_STD=IL_STD
    if(outputFileType=='pdf'){    pdf(outputname)}
    if(outputFileType=='jpg'){    jpeg(outputname)}
    if(outputFileType=='png'){    png(outputname)}
    xin=c(1:dataNull[nrow(dataNull),1])
    yin=mixmdl$lambda[1] * dnorm(xin,mean=mixmdl$mu[1],sd=mixmdl$sigma[1])+mixmdl$lambda[2] * dnorm(xin,mean=mixmdl$mu[2],sd=mixmdl$sigma[2])
    xrange=range(min(data_IL)*0.9,min(max(data_IL)*1.1,IL_Mean*3))
    yrange=range(0,max(yin)*1.2)
    par(mfrow=c(1,2))
    par(mar=rep(3,4))
    par(fig=c(0.1,0.6,0,1))
    cex_main=0.9
    break_num=as.integer(nrow(dataNull)/2)
    if (nrow(dataNull)>500){        break_num=as.integer(nrow(dataNull)/5)}
    plot(xrange,yrange,type='n',frame.plot=F,xlab='Insert Length',cex.main=cex_main,main='Fitted By Bimodal')
    hist(data_IL2,freq=F,col=boxcolor,breaks=break_num,add=T)
    lines(xin,yin,pch='l',col=linecolor)

    par(fig=c(0,0.15,0,1),new=T)
    par(mar=rep(1,4))
    plot(c(0,1),c(0,1),type='n',xlab='',ylab='',axes=FALSE)
    text(0.5,0.5,'Frequency',srt=90)

    par(fig=c(0.5,1,0,1),new=T)
    par(mar=rep(3,4))
    plot(xrange,yrange,type='n',frame.plot=F,xlab='Insert Length',cex.main=cex_main,main='Fitted By Normal')
    hist(data_IL2,freq=F,col=boxcolor,breaks=break_num,xlim=xrange,add=T)
    lines(seq(range(data_IL)[1],range(data_IL)[2],1),dnorm(seq(range(data_IL)[1],range(data_IL)[2],1),IL_Normal_Mean,IL_Normal_STD),type='l',lwd=2,col=linecolor)
    dev.off()

    StatMatrix1=data.frame(Mean=IL_Mean, Median=IL_Median, STD=IL_STD)
    StatMatrix2=data.frame(Bimodal1=mixmdl$lambda[1],Mean1=mixmdl$mu[1], STD1=mixmdl$sigma[1])
    StatMatrix3=data.frame(Bimodal2=mixmdl$lambda[2],Mean2=mixmdl$mu[2], STD2=mixmdl$sigma[2])
    StatMatrix4=data.frame(Normal=1, Mean=IL_Normal_Mean, STD=IL_Normal_STD)
    StatMatrix=rbind(colnames(StatMatrix1),c(StatMatrix1[1,]),
    colnames(StatMatrix2),c(StatMatrix2[1,]),
    colnames(StatMatrix3),c(StatMatrix3[1,]),
    colnames(StatMatrix4),c(StatMatrix4[1,]))
    write.table(file=outputname2,StatMatrix,quote=F,col.name=F,row.name=F)
    }

if (StatType=='TBNull'){
    library(mixtools)
    data_TB=c()
    for(i in 1:nrow(dataNull)){
    data_TB=c(data_TB,rep(dataNull[i,1],dataNull[i,2]))
    }
    total_length=length(data_TB)
    data_TB2=data_TB[as.integer(total_length/40): as.integer(total_length-total_length/40)]
    TB_Mean=mean(data_TB2)
    TB_Median=median(data_TB2)
    TB_STD=sd(data_TB2)
    mixmdl = normalmixEM(data_TB2)
    TB_Normal_Mean=TB_Mean
    TB_Normal_STD=TB_STD
    if(outputFileType=='pdf'){    pdf(outputname)}
    if(outputFileType=='jpg'){    jpeg(outputname)}
    if(outputFileType=='png'){    png(outputname)}
    xin=c(1:dataNull[nrow(dataNull),1])
    yin=mixmdl$lambda[1] * dnorm(xin,mean=mixmdl$mu[1],sd=mixmdl$sigma[1])+mixmdl$lambda[2] * dnorm(xin,mean=mixmdl$mu[2],sd=mixmdl$sigma[2])
    xrange=range(min(data_TB)*0.9,min(max(data_TB)*1.1,TB_Mean*3))
    yrange=range(0,max(yin)*1.2)
    par(mfrow=c(1,2))
    par(mar=rep(3,4))
    par(fig=c(0.1,0.6,0,1))
    cex_main=0.9
    break_num=as.integer(nrow(dataNull)/2)
    if (nrow(dataNull)>500){        break_num=as.integer(nrow(dataNull)/5)}
    plot(xrange,yrange,type='n',frame.plot=F,xlab='Insert Length',cex.main=cex_main,main='Fitted By Bimodal')
    hist(data_TB2,freq=F,col=boxcolor,breaks=break_num,add=T)
    lines(xin,yin,pch='l',col=linecolor)

    par(fig=c(0,0.15,0,1),new=T)
    par(mar=rep(1,4))
    plot(c(0,1),c(0,1),type='n',xlab='',ylab='',axes=FALSE)
    text(0.5,0.5,'Frequency',srt=90)

    par(fig=c(0.5,1,0,1),new=T)
    par(mar=rep(3,4))
    plot(xrange,yrange,type='n',frame.plot=F,xlab='Insert Length',cex.main=cex_main,main='Fitted By Normal')
    hist(data_TB2,freq=F,col=boxcolor,breaks=break_num,xlim=xrange,add=T)
    lines(seq(range(data_TB)[1],range(data_TB)[2],1),dnorm(seq(range(data_TB)[1],range(data_TB)[2],1),TB_Normal_Mean,TB_Normal_STD),type='l',lwd=2,col=linecolor)
    dev.off()

    StatMatrix1=data.frame(Mean=TB_Mean, Median=TB_Median, STD=TB_STD)
    StatMatrix2=data.frame(Bimodal1=mixmdl$lambda[1],Mean1=mixmdl$mu[1], STD1=mixmdl$sigma[1])
    StatMatrix3=data.frame(Bimodal2=mixmdl$lambda[2],Mean2=mixmdl$mu[2], STD2=mixmdl$sigma[2])
    StatMatrix4=data.frame(Normal=1, Mean=TB_Normal_Mean, STD=TB_Normal_STD)
    StatMatrix=rbind(colnames(StatMatrix1),c(StatMatrix1[1,]),
    colnames(StatMatrix2),c(StatMatrix2[1,]),
    colnames(StatMatrix3),c(StatMatrix3[1,]),
    colnames(StatMatrix4),c(StatMatrix4[1,]))
    write.table(file=outputname2,StatMatrix,quote=F,col.name=F,row.name=F)
    }

max_stars_repo_name: mills-lab/svelter

max_issues_repo_name: mills-lab/svelter

max_forks_repo_name: mills-lab/svelter


text:
.fp \np gr
.cs gr 20
.lg 0


text:
stdin <- file("stdin")
open(stdin)

invisible(readLines(stdin, n = 1, warn = FALSE)) # Consumes line to ignore `n`
x <- as.integer(strsplit(trimws(
  readLines(stdin, n = 1, warn = FALSE), which = "both"), " ")[[1]])
w <- as.integer(strsplit(trimws(
  readLines(stdin, n = 1, warn = FALSE), which = "both"), " ")[[1]])

close(stdin)

weighted_mean <- function(x, w) {
  sum(x * w) / sum(w)
}
cat(sprintf("%.1f", weighted_mean(x, w)), "\n")

max_stars_repo_name: vitorgt/problem-solving

max_issues_repo_name: vitorgt/problem-solving

max_forks_repo_name: vitorgt/problem-solving


text:
myRule {
       msiPidDeleteHandle(*handle);
       writeLine("stdout", "*handle");
}

INPUT *handle="21.T12995/78bea5ee-21a9-40e6-94a1-8374c8cbcf93"
OUTPUT ruleExecOut

max_stars_repo_name: sara-nl/irods-msi-persistent-id

max_issues_repo_name: sara-nl/irods-msi-persistent-id

max_forks_repo_name: sara-nl/irods-msi-persistent-id


text:
# RICHIESTE
prima_1 <- 9.999
prima_2 <- 13.328
seconda_1 <- 6.001
seconda_2 <- 8.999
terza <- 5.176

regions <- c(1, 6, 9, 10, 11.382, 13.329, 15.348, 18)
values <- c(0.004, 0.167, 0.168, 0.093, 0.161, 0.155, 0.109, 0.143)

# PROGRAMMA PRINCIPALE
#
#
regions <- c(c(-1), regions, c(30))
values <- c(c(0), values)
plot.new()
EPSILON <- 0.001
f_values <- c(0)

for(i in 2:length(values)) {
  f_values <- cbind(f_values,f_values[i-1] + values[i])
}

x <- seq(0, 20, by = EPSILON)

groups <- cut(x, regions)
fx <- f_values[groups]

par(new=TRUE)
plot(fx, x = x, type="l", xlab = "x", ylab = "F(x)", cex = EPSILON, col = "Blue4", bg = "transparent")

for(i in 2:length(regions)) {
  abline(v = regions[i], lty="dotted", col = "Gray")
  points(regions[i], f_values[i], col = "Blue4", bg = "Blue4", pch = 21, cex = 0.8)
  points(regions[i], f_values[i-1], col = "Blue4", pch = 21, cex = 0.8)
}

f <- function(x) {
  ifelse(x < regions[2], f_values[1],
    ifelse(x >= regions[2] && x < regions[3], f_values[2],
      ifelse(x >= regions[3] && x < regions[4], f_values[3],
        ifelse(x >= regions[4] && x < regions[5], f_values[4],
          ifelse(x >= regions[5] && x < regions[6], f_values[5],
            ifelse(x >= regions[6] && x < regions[7], f_values[6],
               ifelse(x >= regions[7] && x < regions[8], f_values[7],
                  ifelse(x >= regions[8] && x < regions[9], f_values[8],
                    ifelse(x >= regions[9], fvalues[9], NaN)))))))))
}

print(f(prima_2)-f(prima_1))
print(f(seconda_2))
print(1-f(terza))

max_stars_repo_path: code/Esercizio 24.r

max_stars_repo_name: mfranzil/PSUniTN

max_issues_repo_path: code/Esercizio 24.r

max_issues_repo_name: mfranzil/PSUniTN

max_forks_repo_path: code/Esercizio 24.r

max_forks_repo_name: mfranzil/PSUniTN


text:
library(gtools)
library(tibble)
library(dplyr)
library(stringr)
inp <- paste(readLines('~/Documents/CodeRepositories/AdventofCode/inputs2020/Advent2020_01.txt'), collapse = ",")
inp <- strsplit(inp, split = ",")
inp <- as.integer(inp[[1]])

combos <- as_tibble(gtools::combinations(n = 200, r = 2, v = inp)) %>%
  rowwise() %>%
  filter(sum(c(V1, V2)) == 2020) %>%
  mutate(product = V1 * V2)

print(stringr::str_interp("AoC 2020 Day 1, Part 1 answer is (${V1}, ${V2}) is ${product}", combos[1,]))

combos <- as_tibble(gtools::combinations(n = 200, r = 3, v = inp)) %>%
  rowwise() %>%
  filter(sum(c(V1, V2, V3)) == 2020) %>%
  mutate(product = V1 * V2 * V3)

print(stringr::str_interp("AoC 2020 Day 1, Part 2 answer is (${V1}, ${V2}, ${V3}) is ${product}", combos[1,]))


text:
str <- "Lorem ipsum dolor sit amet,
consectetur adipiscing elit,
sed do eiusmod tempor incididunt
ut labore et dolore magna aliqua."

print(str)

print(cat(str))

hi <- "Hello"
print(nchar(hi))

str1 <- "Hello World!"
print(grepl("H", str1))
print(grepl("Hello", str1))
print(grepl("X", str1))

str2 <- "Hello"
str3 <- "World"
print(paste(str2, str3))

max_stars_repo_name: teddythinh/R-language-Practice

max_issues_repo_name: teddythinh/R-language-Practice

max_forks_repo_name: teddythinh/R-language-Practice


text:
###############################
###  2010, David Ellinghaus ###
###############################

rm(list=ls())

file.het   = commandArgs()[4]
file.imiss = commandArgs()[5]
fraction_miss.thresh = as.numeric(commandArgs()[6])

# heterozygosity +- 5SD way from mean
sd.times = 5

imiss = read.table(file.imiss,header=T)
het   = read.table(file.het,header=T)

imiss.log10F_MISS = log10(imiss$F_MISS)
het = 1 - het$O.HOM./het$N.NM.

png(file=paste(file.het, ".logscale.2.png", sep=""), width=960, height=960)

dcols = densCols(imiss.log10F_MISS,het)
plot(imiss.log10F_MISS,het, xlim=c(-3,0), ylim=c(0,0.5), xlab="Individual missingness", ylab="Heterozygosity rate", col=dcols, pch=20, axes=F)
axis(1, at=c(-3,-2,log10(fraction_miss.thresh),-1,0), labels=c(0.001,0.01,fraction_miss.thresh,0.1,1))
axis(2, at=seq(0,0.5,0.05), tick=T)
abline(h=mean(het)-(sd.times*sd(het)), lty=2, col="red")
abline(h=mean(het)+(sd.times*sd(het)), lty=2, col="red")
abline(v=log10(fraction_miss.thresh), col="red", lty=2)

dev.off()

max_stars_repo_name: jkaessens/gwas-assoc

max_issues_repo_name: jkaessens/gwas-assoc

max_forks_repo_name: jkaessens/gwas-assoc


text:
setClass("FLIndex.SMS",
    contains="FLIndex",
    representation(
        range.SMS="vector"
    ) ,
     prototype=prototype(range.SMS=list(season=1, power.age=-1, q.age=0,
        var.age.group=as.vector(0,mode="list"),minCV=0.3))
)

FLIndex.SMS <- function(name=character(0), desc=character(0), distribution=character(0),
    type=character(0), startf=NA, endf=NA, plusgroup=NA, season=NA, power.age=NA,
    q.age=NA,var.age.group=NA,minCV=0.3, ...) {

    args <- list(...)
    if(length(args)==0)
        args <- list(index=FLQuant())

    dimnames <- dimnames(args[[names(lapply(args, is.FLQuant)==TRUE)[1]]])
    sdimnames <- dimnames
    sdimnames[1] <- "all"

    if(!is.FLQuant(args['index']))
        index <- FLQuant(dimnames=dimnames)

    dims <- dims(index)

    new <- new("FLIndex.SMS", name = name, desc = desc, distribution = distribution,
        type=type,
        index = index, index.var = FLQuant(dimnames=dimnames),
        index.q = FLQuant(dimnames=dimnames), sel.pattern = FLQuant(dimnames=dimnames),
        catch.n = FLQuant(dimnames=dimnames), catch.wt = FLQuant(dimnames=dimnames),
        effort = FLQuant(dimnames=sdimnames),
        range = unlist(list(min=dims$min, max=dims$max,
        plusgroup=NA, minyear=dims$minyear, maxyear=dims$maxyear, startf=startf, endf=endf)))
        range.SMS = unlist(list(season=season, power.age=power.age,
                 q.age=q.age,var.age.group=var.age.group, minCV=minCV))

    # Add extra arguments
    for(i in names(args)[names(args)!='iniFLQuant'])
        slot(new, i) <- args[[i]]

    return(new)
}

SMS2FLIndices<-function(control,path=NULL,fleet.inf="fleet_info.dat",fleet.index="fleet_catch.in",
                        fleet.name="fleet_names.in") {

    old.wd<-getwd()
    if (!is.null(path)) setwd(path)

    nsp<-slot(control,"no.species")
    nq<-slot(control,"last.season")

    #count number of other predators
    info<-slot(control,"species.info")[,"predator"]
    no.oth<-sum(info==2)
    nsp<-nsp-no.oth

    s<-readLines(fleet.name, n=1000)
    s<-gsub('_',' ',s)
    fl.names<-sub('[[:space:]]+$', '', s)

    info<-scan(fleet.inf,comment.char = "#",quiet=TRUE)
    minCV<-info[1]
    i<-2
    n.fleet<-as.vector(info[i:(i-1+nsp)])
    i<-i+nsp
    sum.fleet<-sum(n.fleet)
    fl.info<-matrix(info[i:(i-1+sum.fleet*10)],ncol=10,nrow=sum.fleet,byrow=TRUE)
    i<-i+sum.fleet*10

    sum.var.age<-sum(fl.info[,10])
    fl.var<-as.vector(info[i:(i-1+sum.var.age)])

    CE<-scan(fleet.index,comment.char = "#",quiet=TRUE)

    # creates empty FLIndices object
    FLIndices. <- FLIndices()

    i<-1
    v<-1
    sp.fl<-0
    for (sp in 1:nsp) {
     for (fl in 1:n.fleet[sp]) {
      sp.fl<-sp.fl+1
      fy<-fl.info[sp.fl,1]
      ly<-fl.info[sp.fl,2]
      alfa<-fl.info[sp.fl,3]
      beta<-fl.info[sp.fl,4]
      fa<-fl.info[sp.fl,5]
      la<-fl.info[sp.fl,6]
      la.q<-fl.info[sp.fl,7]
      la.p<-fl.info[sp.fl,8]
      seas<-fl.info[sp.fl,9]
      n.var<-fl.info[sp.fl,10]

      nyr<-ly-fy+1
      nages<-la-fa+1

       # template for input to quant

      dim<-c(nages,nyr,1,1,1,1)
      dim2<-c(1,nyr,1,1,1,1)
      dimnames<-list(age=fa:la,year=fy:ly,unit="all",season=seas,area="all",iter="none")
      dimnames2<-list(age="all",year=fy:ly,unit="all",season=seas,area="all",iter="none")

      tmp<-matrix(CE[i:(nyr*(nages+1)+i-1)],ncol=nages+1,nrow=nyr,byrow=TRUE)

      effort<-array(tmp[,1],dim=dim2,dimnames=dimnames2)
      catch<-matrix(tmp[,2:(nages+1)],ncol=nyr,nrow=nages,byrow=TRUE)
      #print(catch)
      catch<-array(catch,dim=dim,dimnames=dimnames)
      #print(catch)
      index<-catch/rep(effort,each=nages)

      indx<-FLIndex.SMS(index=as.FLQuant(index),effort=as.FLQuant(effort),catch.n=as.FLQuant(catch),
                 name = fl.names[sp.fl], desc = fl.names[sp.fl])

      indx@range<-unlist(list("min"=fa,"max"=la,"plusgroup"=NA,"minyear"=fy,"maxyear"=ly,"startf"=alfa,"endf"=beta))

      indx@range.SMS<-list("season"=seas, "power.age"=la.p,
                            "q.age"=la.q,"var.age.group"=as.vector(fl.var[v:(v-1+n.var)]),"minCV"=minCV)

      v<-v+n.var
      i<-i+nyr*(nages+1)
      FLIndices.[[sp.fl]]<-indx
    }
   }
   setwd(old.wd)
   FLIndices.
 }

FLIndices2SMS<-function(out.path=NULL,indices=NULL,control=NULL,fleet.inf="fleet_info.dat",
                        fleet.index="fleet_catch.in",fleet.name="fleet_names.in") {
    old.wd<-getwd()
    if (!is.null(out.path)) setwd(out.path)

    if (is.null(indices))
        stop("A 'FLIndices' must be given")

    if (!inherits(indices, "FLIndices"))
        stop("indices must be an 'FLIndices' object!")

    for (i in 1:length(indices)) {
     if (is.na(indices[[i]]@range["startf"]) || is.na(indices[[i]]@range["endf"]))
         stop(paste("Must supply startf & endf for range in FLIndex",i))

      if (!all(names(indices[[i]]@range) == c("min","max","plusgroup","minyear","maxyear","startf","endf")))
         stop("Range must have names 'min','max','plusgroup','minyear','maxyear','startf','endf'")
    }

    if (!inherits(control, "FLSMS.control"))
    stop("control must be an 'FLSMS.control' object!")

    if (!validObject(control)) stop("control is not valid!")

    nsp<-slot(control,"no.species")
    #count number of other predators
    info<-slot(control,"species.info")[,7]
    no.oth<-sum(info==2)
    nsp<-nsp-no.oth  # no of VPA species

    first.year<-slot(control,"first.year")
    last.year<-slot(control,"last.year.model")
    last.season<-slot(control,"last.season")
    n.season<-last.season

    no.indices<-rep(0,nsp)
    info<-matrix(0,ncol=10,nrow=length(indices))
    fl.name<-rep('',length(indices))
    v.age<-list()

    sp<-1; n<-1
    old.sp<-substr(indices[[1]]@desc,1,3)

    cat("# file fleet_catch.in\n",file=fleet.index)
    for (idc in indices) {
        fl.name[n]<-idc@name
        sp.name<-substr(idc@desc,1,3)
        if (nsp>1 & sp.name!=old.sp) {sp<-sp+1; old.sp<-sp.name}
        cat("# ",sp.name,",",fl.name[n],"\n",file=fleet.index,append=TRUE)
        no.indices[sp]=no.indices[sp]+1
        range<-idc@range
        info[n,1]<-range["minyear"]
        info[n,2]<-range["maxyear"]
        info[n,3]<-range["startf"]
        info[n,4]<-range["endf"]
        info[n,5]<-range["min"]
        info[n,6]<-range["max"]
        info[n,7]<-idc@range.SMS$q.age
        info[n,8]<-idc@range.SMS$power.age
        info[n,9]<-idc@range.SMS$season
        info[n,10]<-length(idc@range.SMS$var.age.group)
        v.age<-c(v.age,list(idc@range.SMS$var.age.group))
        minCV<-idc@range.SMS$minCV
        write.table(cbind(as.vector(idc@effort),
                    t(matrix(idc@catch.n,ncol=info[n,2]-info[n,1]+1,byrow=FALSE))),
                    file=fleet.index,row.names=FALSE,col.names=FALSE,quote=FALSE,append=TRUE)
        n<-n+1
    }
    cat("-999 #  Check value\n",file=fleet.index,append=TRUE)

    cat(paste("# file: fleet_info.dat\n",minCV," #min CV of CPUE observations\n"),file=fleet.inf)
    cat("# number of fleets by species\n",file=fleet.inf,append=TRUE)
    write(no.indices,file=fleet.inf,ncolumns=nsp,append=TRUE)
    cat("#############",
        "\n# 1-2, First year last year,",
        "\n# 3-4. Alpha and beta - the start and end of the fishing period for the fleet given as fractions of the season (or year if annual data are used),",
        "\n# 5-6   first and last age,",
        "\n# 7.   last age with age dependent catchability,",
        "\n# 8.   last age for stock size dependent catchability (power model), -1 indicated no ages uses power model,",
        "\n# 9.   season for survey,",
        "\n# 10.  number of variance groups for estimated cathability,",
        "\n# by species and fleet",
        "\n#############\n", file=fleet.inf,append=TRUE)
    i<-0
    for (s in (1:nsp)) {
      cat("# ",control@species.names[s+no.oth],"\n",file=fleet.inf,append=TRUE)
      for (id in (1:no.indices[s])) {
        i<-i+1
        cat("# ",fl.name[i],"\n",file=fleet.inf,append=TRUE)
        write(file=fleet.inf ,info[i,],ncolumns=10,append=TRUE)
      }
    }

#          write.table(file=fleet.inf ,info,row.names=FALSE,col.names=FALSE,quote=FALSE,append=TRUE)

    cat("# variance groups\n",file=fleet.inf,append=TRUE)
    #for (a in v.age) write(a,file=fleet.inf,append=TRUE)
    i<-0
    for (s in (1:nsp)) {
      cat("# ",control@species.names[s+no.oth],"\n",file=fleet.inf,append=TRUE)
      for (id in (1:no.indices[s])) {
        i<-i+1
        cat("# ",fl.name[i],"\n",file=fleet.inf,append=TRUE)
        write(file=fleet.inf ,v.age[[i]],append=TRUE)
      }
    }
    cat("-999 #  Check value\n",file=fleet.inf,append=TRUE)

    fl.name<-substr(paste(fl.name,"__________________________",sep=''),1,26)
    fl.name<-gsub(' ','_',fl.name)
    write(fl.name,file=fleet.name)

    setwd(old.wd)
}

   # test FLIndices2SMS(out.path=out.path,indices=SMS.indices,control=SMS.dat)


text:
error_foo <- function() {
  cat(geterrmessage(), file = stderr());
  q("no", 1, F)
}
options(show.error.messages = F, error = error_foo)

# we need that to not crash galaxy with an UTF8 error on German LC settings.
loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")

suppressPackageStartupMessages({
    library("pathview")
    library("optparse")
})

sessionInfo()

option_list <- list(
  make_option(c("--pathway_id"), type = "character", default = NULL, help = "Path to tabular file with gene data"),
  make_option(c("--pathway_id_fp"), type = "character", default = NULL, help = "Path to tabular file with pathway ids"),
  make_option(c("--pathway_id_header"), type = "logical", default = FALSE, help = "Header for tabular file with pathway ids"),
  make_option(c("--species"), type = "character", default = "hsa", help = "KEGG code, scientific name or the common name of the species"),
  make_option(c("--gene_data"), type = "character", default = NULL, help = "Path to tabular file with gene data"),
  make_option(c("--gd_header"), type = "logical", default = FALSE, help = "Header for the gene data file"),
  make_option(c("--gene_idtype"), type = "character", default = "entrez", help = "ID type used for the gene data"),
  make_option(c("--cpd_data"), type = "character", default = NULL, help = "Path to tabular file with compound data"),
  make_option(c("--cpd_header"), type = "logical", default = FALSE, help = "Header for the compound data file"),
  make_option(c("--cpd_idtype"), type = "character", default = "kegg", help = "ID type used for the compound data"),
  make_option(c("--multi_state"), type = "logical", default = TRUE, help = "Are the gene and compound data paired?"),
  make_option(c("--match_data"), type = "logical", default = TRUE, help = "Are the gene and compound data paired?"),
  make_option(c("--kegg_native"), type = "logical", default = TRUE, help = "Render pathway graph as native KEGG grap? Alternative is the Graphviz layout"),
  make_option(c("--same_layer"), type = "logical", default = TRUE, help = "Plot on same layer?"),
  make_option(c("--map_null"), type = "logical", default = TRUE, help = "Map the NULL gene or compound data to pathway?"),
  make_option(c("--split_group"), type = "logical", default = FALSE, help = "Split node groups into individual nodes?"),
  make_option(c("--expand_node"), type = "logical", default = FALSE, help = "Expand multiple-gene nodes into single-gene nodes?"),
  make_option(c("--sign_pos"), type = "character", default = "bottomright", help = "Position of pathview signature")
)

parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
args <- parse_args(parser)
print(args)

read_table <- function(fp, header, rownames = 1, colclasses = NA) {
  table <- read.table(fp, header = header, sep = "\t", row.names = rownames, colClasses = colclasses)
  # transform to vector if only one column
  if (dim(table)[2] == 1) {
    names <- rownames(table)
    table <- table[, 1]
    names(table) <- names
  }
  return(table)
}

get_table <- function(fp, header) {
  table <- NULL
  if (!is.null(fp)) {
    table <- read_table(fp, header, rownames = 1)
  }
  return(table)
}

# load gene_data file
gene_data <- get_table(args$gene_data, args$gd_header)

# load compound data file
cpd_data <- get_table(args$cpd_data, args$cpd_header)

run_pathview <- function(pathway_id) {
  pathview(
    pathway.id = pathway_id,
    gene.data = gene_data,
    gene.idtype = args$gene_idtype,
    cpd.data = cpd_data,
    cpd.idtype = args$cpd_idtype,
    species = args$species,
    multi.state = args$multi_state,
    match.data = args$match_data,
    kegg.native = args$kegg_native,
    same.layer = args$same_layer,
    split.group = args$split_group,
    expand.node = args$expand_node,
    sign.pos = args$sign_pos,
    map.null = args$map_null)
}

# get pathway ids
if (!is.null(args$pathway_id)) {
  run_pathview(args$pathway_id)
} else {
  pthws <- read_table(args$pathway_id_fp, args$pathway_id_header, rownames = NULL, colclasses = "character")
  for (p in pthws) {
    run_pathview(p)
  }
}

max_stars_repo_name: crashfrog/tools-iuc

max_issues_repo_name: Delphine-L/tools-iuc

max_forks_repo_name: Delphine-L/tools-iuc


text:
#' Load a matrix from a simple dense table with genes in rows and cells in columns. First column is the gene name
#'
#' @param mat_nm - the name of the new matrix in scdb
#' @param fn  path of the matrix to read (columns names should represent Cell.IDs, rownames should represent gene names, values are expected to represent number of UMIs per gene per Cell)
#' @param meta_fn meta data table to read. If this is null, all cells will be onsidred as part of a single Amp.Batch and Dataset.ID. The table must include at least the fields "Cell.ID", "Amp.Batch.ID", "Seq.Batch.ID", "Batch.Set.ID", where Cell.ID should match the names of cells in the data matrix.
#' @param force - if true, will import from 10x files even when the matrix is present in the DB
#'
#' @export

mcell_import_scmat_tsv = function(mat_nm, fn, genes_fn = NULL, meta_fn = NULL, dset_nm = NULL, force = FALSE, paralogs_policy = "sum")
{
	if(!scdb_is_valid()) {
		stop("MCERR - scdb is not initialized, cannot import")
	}
	if(!force & !is.null(scdb_mat(mat_nm))) {
		return(TRUE)
	}
	if(!file.exists(fn)) {
		stop("MCERR - umi matrix file ", fn, " does not exist")
	}
	umis = fread(fn)
	if(!is.null(genes_fn)) {
		all_ids = unlist(umis[,1])

		genes = as.data.frame(fread(genes_fn, header=F, stringsAsFactors = F))

		colnames(genes) = c('id', 'name')
		rownames(genes) = genes$id
		gname = genes[all_ids, 'name']
		umis = umis[!is.na(gname),]
		gname = gname[!is.na(gname)]

		unique_genes = names(which(table(gname) == 1))
		umis1 = umis[gname %in% unique_genes, -1]
		rownames(umis1) = gname[gname %in% unique_genes]
		if (paralogs_policy == 'remove') {
			umis = umis1
		}
		else if (paralogs_policy == 'sum') {
			non_unique_genes = setdiff(gname, unique_genes)
			umis2 = as.matrix(umis[gname %in% non_unique_genes, -1])

			message(sprintf("summing up total of %d paralog genes into %d unique genes", nrow(umis2), length(non_unique_genes)))

		   dup_genes = gname[gname %in% non_unique_genes]

			umis2s = apply(umis2, 2, function(x) {
					tapply(x, INDEX=dup_genes, FUN=sum) } )

			umis = rbind(umis1, umis2s)
			rownames(umis) = c(rownames(umis1), rownames(umis2s))
		}
	} else {
		rownames(umis) = umis[,1]
		umis = umis[,-1]
	}
	if(!is.null(meta_fn)) {
		md = fread(meta_fn, sep="\t")
		mandatory = c("Cell.ID", "Amp.Batch.ID", "Seq.Batch.ID", "Batch.Set.ID")
		miss_f = setdiff(mandatory, colnames(md))
		if(length(miss_f)>0) {
			stop("MC-ERR: missing fields in matrix metadata : ", miss_f)
		}
		rownames(md) = md$Cell.ID
		md = md %>% select(-Cell.ID)
	} else {
		if(is.null(dset_nm)) {
			stop("either provide a meta data file, or supply a dataset name (parameter dset_nm) and we'll fake it for you...")
		}
		dsets = data.frame(batch_set_id = dset_nm,
							 	 amp_batch_id = dset_nm,
							    seq_batch_id = dset_nm)
		md = matrix(c('import', dset_nm, dset_nm, dset_nm),
				nrow=ncol(umis), ncol=4, byrow=T,
				dimnames=list(colnames(umis),
							c('type', 'batch_set_id', 'amp_batch_id', 'seq_batch_id')))
		md = as.data.frame(md)
	}
	sparse_m = Matrix(as.matrix(umis))
	rownames(sparse_m) = rownames(umis)
	scmat = scm_new_matrix(sparse_m, md)
	scdb_add_mat(mat_nm, scmat)
	return(TRUE)
}


text:
#'@title calclcb
#'
#'@description Calculate longitudinal position of center of buoyancy
#'(\code{lcb}).
#'
#'@param lwl Waterline length (vector of numericals, m) (see \code{\link{calclwl}})
#'
#'@details
#'The longitudinal position of center of buoyancy describes center of buoyancy
#'forward of half of the waterline length (\code{lwl}) (decimal percentage of
#'waterline length).
#'
#'@return \code{lcb} (vector of numericals)
#'
#'@references
#'\href{http://hdl.handle.net/11250/2410741}{Rakke, S. G. 2016. "Ship Emissions
#'Calculation from AIS." NTNU.}
#'
#'@seealso
#'\code{\link{calclwl}}
#'
#'@examples
#'calclcb(c(218.75,400,262))
#'
#'@export

calclcb<- function(lwl){

  lcb<- -0.75*(lwl/2)/100

  return(lcb)
}


text:
test_that("multiplication works", {
  expect_equal(2 * 2, 4)
})

max_stars_repo_name: apsteinmetz/spotfuzz

max_issues_repo_name: kkkkkkkh/khlearn

max_forks_repo_name: kkkkkkkh/khlearn


text:
validate.catch <- TRUE
validate.fisheries <- TRUE
validate.mark.info <- TRUE

validate.stocks <- TRUE
validate.escapment.flags <- TRUE

max_stars_repo_name: methiess/PSC-FRAM-Admin

max_issues_repo_name: methiess/PSC-FRAM-Admin

max_forks_repo_name: methiess/PSC-FRAM-Admin


text:
rulemsiCollReplAsynch {
#Input parameters are:
#  Data Object path
#  Optional flags in form keyword=value
#    destRescName - the target resource for the replica
#    backupRescName - specifies use of the resource for the backup mode
#    rescName - the resource holding the source data
#    updateRepl= - specifies all replicas will be updated
#    replNum - specifies the replica number to use as the source
#    numThreads - specifies the number of threads to use for transmission
#    all - specifies to replicate to all resources in a resource group
#    irodsAdmin - enables administrator to replicate other users' files
#    verifyChksum - verify the transfer using checksums
#    rbudpTransfer - use Reliable Blast UDP for transport
#    rbudpSendRate - the transmission rate in kbits/sec, default is 600 kbits/sec
#    rbudpPackSize - the packet size in bytes, default is 8192
#Output parameter is:
#  Status
#Output from running the example is:
#  The file /tempZone/home/rods/sub1/foo3 is replicated onto resource testResc

	delay (*DelayInfo) {
		msiCollRepl(*SourceFile,"destRescName=*Resource",*Status);
		writeLine("stdout","The collection *SourceFile is replicated onto resource *Resource");
	}

  }
INPUT  *DelayInfo="<PLUSET>1m</PLUSET><EF>24h</EF>",*SourceFile="/tempZone/home/rods/sub1/foo3",*Resource="testResc"
OUTPUT ruleExecOut

max_stars_repo_name: korydraughn/jargon

max_issues_repo_name: korydraughn/jargon

max_forks_repo_name: korydraughn/jargon


text:
set.seed(22)

## Generate N random integers that sum to M
rand_vect <- function(N, M, sd = 1, pos.only = TRUE) {
  vec <- rnorm(N, M/N, sd)
  if (abs(sum(vec)) < 0.01) vec <- vec + 1
  vec <- round(vec / sum(vec) * M)
  deviation <- M - sum(vec)
  for (. in seq_len(abs(deviation))) {
    vec[i] <- vec[i <- sample(N, 1)] + sign(deviation)
  }
  if (pos.only) while (any(vec < 0)) {
    negs <- vec < 0
    pos  <- vec > 0
    vec[negs][i] <- vec[negs][i <- sample(sum(negs), 1)] + 1
    vec[pos][i]  <- vec[pos ][i <- sample(sum(pos ), 1)] - 1
  }
  return(vec)
}

## silumale missing indicators under MCAR
## missing completely at random; omega is the missing indicator: 1: missing; omega is a matrix of m x n; this function controls that their is not a whole row/col missing
f.omega.mcar = function(m, n, prob){
  r = round( min(m,n) / 2)
  ind = matrix(rbinom(r^2, size=1, prob),r,r)
  repcol = rand_vect(r,m)
  reprow = rand_vect(r,n)
  omega = ind[rep(1:r,repcol), rep(1:r,reprow)]
  return(omega)
}

#### mcar: silumale missing indicators
## missing completely at random
## omega is the missing indicator: 1: missing
## omega is a matrix of m x n
## when missing prob is hing, it's likely to have a whole row/col missing
f.omega.mcar = function(m, n, mprob){
  sim=TRUE
  while(sim){
    omega=matrix(rbinom(m*n,size=1,mprob),m,n)
    sim = sum(colSums(omega)==m) | sum(rowSums(omega)==n)
  }
  return(omega)
}

## silumale missing indicators under MAR
f.omega.mar = function(Y, prob=0.2, method){
  m_Y= dim(Y)[1]; n_Y = dim(Y)[2]
  ## missing at random
  if(method=="sv"){
    ## missing according to magnatitude of singular values
    ## more likely to miss when the biggest dj contributes more to a point?
    ## omega is the missing indicator: 1: missing
    ## omega is a matrix of m x n
    ## Ysvd = svd(Y)
    Ysvd = svd(Y)
    U = Ysvd$u
    D = Ysvd$d
    V = Ysvd$v

    notmaxD = which(D < max(D))
    Dnew = D
    Dnew[notmaxD] = 0
    Ynew = U %*% diag(Dnew) %*% t(V)

    beta = 10
    f.alpha = function(alpha, beta, Y){
      m_Y = dim(Y)[1]; n_Y = dim(Y)[2]
      a = sum(apply(Y, c(1,2), function(x){1/(1+exp(-alpha-beta*x))} ))
      return(a - prob * m_Y * n_Y)
    }

    alpha = uniroot(f.alpha, beta=beta, Y=Ynew, c(-100000000, 100000000))$root
    omega=apply(Ynew,c(1,2),function(x){ p=1/(1+exp(-alpha-beta*x)); rbinom(1, size=1, prob=p)})

  }else if(method=="corr"){
    ## each column of the data Y missing depends on the column that's most highly correlated with that column
    corY = cor(Y)
    diag(corY) = 0
    corYmax = apply(corY, 2, function(x){ x=abs(x);maxX = which(x == max(x));return(maxX)})
    omega = matrix(NA, m_Y,n_Y)
    for(j in 1:n_Y){
      beta = 10
      f.alpha = function(alpha, beta, Y){
        n_Y = dim(Y)[2]
        a = sum(sapply(Y[,corYmax[j]], function(x){1/(1+exp(-alpha-beta*x))} ))
        return(a - prob * m_Y)
      }
      alpha = uniroot(f.alpha,beta=beta, Y=Y, c(-100000000, 100000000))$root
      omega[,j]=sapply(Y[,corYmax[j]],function(x){ p=1/(1+exp(-alpha-beta*x)); rbinom(1, size=1, prob=p)})
    }
  }else{print("please specify method as 'sv'(singular value) or 'corr'")}
  return(omega)
}

args = commandArgs(TRUE)
m = as.numeric(args[1])
n = as.numeric(args[2])
K = as.numeric(args[3])
df = as.numeric(args[4])

print(m)
print(n)
print(K)
print(df)

omega = f.omega.mcar(m=m, n=n, mprob=0.15)
dir.create("./output/omega")
write.table(omega, file=paste("./output/omega/omega_mcar", "_df_", df, "_m_", m, "_n_", n, "_K_", K, ".csv", sep = ""), row.names = FALSE, sep=",", col.names = FALSE)

M = read.csv(paste("./output/simulated_data/M",  "_df_", df, "_m_", m, "_n_", n, "_K_", K, ".csv", sep = ""), header = FALSE)
omega = f.omega.mar(Y=M, prob=0.15, method="sv")
write.table(omega, file=paste("./output/omega/omega_sv", "_df_", df, "_m_", m, "_n_", n, "_K_", K, ".csv", sep = ""), row.names = FALSE, sep=",", col.names = FALSE)

omega = f.omega.mar(Y=M, prob=0.15, method="corr")
write.table(omega, file=paste("./output/omega/omega_corr",  "_df_", df, "_m_", m, "_n_", n, "_K_", K, ".csv", sep = ""), row.names = FALSE, sep=",", col.names = FALSE)

q()


text:
set.seed.if.possible <- function(object) {
	if (!is.null(object$seed)) {
		set.seed(object$seed)
	}
}

determine.positive.class <- function(object) {
	response <- object$adapter$y.vars[[1]]
	if (is.numeric(response)) {
		return("1")
	}
	if (is.logical(response)) {
		return("TRUE")
	}
	if (!is.null(object$positive.class)) {
		return(object$positive.class)
	}
	if (is.factor(response)) {
		return(levels(response)[1])
	}
	if (is.character(response)) {
		return(levels(as.factor(response))[1])
	}
	stop(
		"Response of the model should be integer, logical, factor or character."
	)
}


text:
 # var1s.r

 nn <- 100
 x <- rnorm(nn, sd =2); vx <- var(x)
 chi0 <- (nn-1)*vx / 2^2
 chi0
 p.val <- 1-pchisq(chi0, nn-1)
 p.val
 ci <- c( (nn-1)*vx /qchisq(0.975, nn-1), (nn-1)*vx /qchisq(0.025, nn-1) )
 ci

max_stars_repo_name: tolkien/misc

max_issues_repo_name: tolkien/misc

max_forks_repo_name: tolkien/misc


text:
# EWAS of left-handedness: basic model, NTR adults
# V.Odintsova, J.van Dongen

library(gee)

# Load DNA methylation data (DNAm)
# Load phenotype and covariate data (covariates)

# Variabels for basic model
# familynumber
# sex
# age (age at DNA bloodsampling)
# Neut_Perc (neutrophils percentage)
# Mono_Perc (monocytes percentage)
# Eos_Perc (eosinophils percentage)
# Array_rownum (array rownumber)
# Sample_Plate (sample plate)
# hand01 (left-handedness variable: right-handed = 0, left-handed = 1)

mydata <- data.frame(familynumber, sex,age,Mono_Perc,Eos_Perc,Neut_Perc,Array_rownum,Sample_Plate,hand01)
mydata <- mydata[order(mydata$familynumber),]

# GEE model for a continuous outcome (DNA methylation level) correcting for clustering of samples
#------------------------------------------------------------------------------------------------------------------------------------------------------
geemodel <- function (data)
{
gee(CpGi~hand01+sex+age+Mono_Perc+Eos_Perc+Neut_Perc+Array_rownum+Sample_Plate,data=data, id=familynumber, family=gaussian, corstr="exchangeable", maxiter=100, na.action=na.omit,silent=TRUE)
}
#------------------------------------------------------------------------------------------------------------------------------------------------------

Nestimates <-nrow(coeff)
Ncpgs      <-ncol(DNAm)
pval       <-matrix(NA,Ncpgs,Nestimates)
estimate   <-matrix(NA,Ncpgs,Nestimates)
RobustSE   <-matrix(NA,Ncpgs,Nestimates)
RobustZ    <-matrix(NA,Ncpgs,Nestimates)
N          <- matrix(NA,Ncpgs,1)
Ncpgs
Nestimates

for (i in 1:Ncpgs)
{
    mydata$CpGi <- DNAm[,i]
    N[i,] <- length(which(!is.na(mydata$CpGi)))
    coeff <- summary(geemodel(mydata))$coefficients
    for (j in 1: Nestimates)
{
    estimate[i,j] <- coeff[j,1]
    RobustSE[i,j] <- coeff[j,4]
    RobustZ[i,j]  <- coeff[j,5]
    pval[i,j]     <- 2*pnorm(-abs(coeff[j,5]))
}
}

colnames(estimate)   <- rownames(coeff)
colnames(RobustSE)   <- rownames(coeff)
colnames(RobustZ)    <- rownames(coeff)
colnames(pval)       <- rownames(coeff)
rownames(estimate)   <- colnames(DNAm)
rownames(RobustSE)   <- colnames(DNAm)
rownames(RobustZ)    <- colnames(DNAm)
rownames(pval)       <- colnames(DNAm)
rownames(N)          <- colnames(DNAm)

save(estimate,pval,RobustSE,RobustZ,N, file="EWASresults.RData")

max_stars_repo_path: ntr/ewas_basic model_NTRadults.r

max_stars_repo_name: MRCIEU/handedness-ewas

max_issues_repo_path: ntr/ewas_basic model_NTRadults.r

max_issues_repo_name: MRCIEU/handedness-ewas

max_forks_repo_path: ntr/ewas_basic model_NTRadults.r

max_forks_repo_name: MRCIEU/handedness-ewas


text:
wtd.means = by(firms, size.category,
  function(piece) (sum(piece$mktcap*piece$spread)/sum(piece$mktcap)))

max_stars_repo_name: DeagleGross/Rplugin

max_issues_repo_name: DeagleGross/Rplugin

max_forks_repo_name: DeagleGross/Rplugin


text:
REBOL [
    Title: "Pe-info"
    Date: 10-May-2013/10:14:56+2:00
    Name: none
    Version: none
    File: none
    Home: none
    Author: "A Rebol"
    Owner: none
    Rights: none
    Needs: none
    Tabs: none
    Usage: none
    Purpose: none
    Comment: none
    History: none
    Language: none
    Type: none
    Content: none
    Email: none
	require: [
    	rs-project 'stream-io
	]
]

PE: make stream-io [
	PE-port: none
	PE-header: none

	MSDOS-stub:
	file-header:	none
	section-header: none
	rsrc-table:     none
	rsrc-ptr:       none

	section-headers: copy []

	seekToBuffer: func[bytes][
		;print ["seekToBuffer bytes" bytes]
		setStreamBuffer copy/part PE-port bytes
		PE-port: skip PE-port bytes
		inBuffer
	]
	open: func["Reads PE file into buffer" PE-file [file! url!]][
		;AVImainHeader: AVIstreamHeader: none
		probe PE-file
		PE-port: system/words/open/read/binary/seek PE-file
		seekToBuffer size? PE-file
	]

	read-file-header: does [
		file-header: context [
			machine:          readUI16
			sections-nb:      readUI16
			timestamp:        readUI32		;-- file creation timestamp
			symbols-ptr:      readUI32		;-- OBJ only
			symbols-nb:       readUI32		;-- OBJ only
			opt-headers-size: readUI16		;-- zero for OBJ
			flags:            readUI16
		]
	]
	read-optional-header: does [
		optional-header: context[			;-- optional header for image only
			magic:          readUI16			;-- different for 32/64-bit
			major-link-version: readUI8
			minor-link-version: readUI8
			code-size:       readUI32
			initdata-size:       readUI32
			uninitdata-size:       readUI32
			entry-point-addr:       readUI32
			code-base:       readUI32
			data-base:       readUI32		;-- 32-bit only (remove field for 64-bit)
			image-base:       readUI32		;-- 8 bytes for 64-bit
			memory-align:       readUI32
			file-align:       readUI32
			major-OS-version: readUI16			;-- should be 4.0
			minor-OS-version: readUI16
			major-img-version: readUI16
			minor-img-version: readUI16
			major-sub-version: readUI16
			minor-sub-version: readUI16
			win32-ver-value:       readUI32		;-- reserved, must be zero
			image-size:       readUI32
			headers-size:       readUI32
			checksum:        readUI32	;-- for drivers and DLL only
			sub-system: readUI16
			dll-flags: readUI16		;-- DLL only
			stack-res-size:  readUI32		;-- 8 bytes for 64-bit
			stack-com-size:  readUI32		;-- 8 bytes for 64-bit
			heap-res-size:   readUI32		;-- 8 bytes for 64-bit
			heap-com-size:   readUI32		;-- 8 bytes for 64-bit
			loader-flags:    readUI32		;-- reserved, must be zero
			data-dir-nb:     readUI32
			;-- Data Directory
			export-addr:     readUI32
			export-size:     readUI32
			import-addr:     readUI32
			import-size:     readUI32
			rsrc-addr:       readUI32
			rsrc-size:       readUI32
			except-addr:     readUI32
			except-size:     readUI32
			cert-addr:       readUI32
			cert-size:       readUI32
			reloc-addr:      readUI32
			reloc-size:      readUI32
			debug-addr:      readUI32
			debug-size:      readUI32
			arch-addr:       readUI32		;-- reserved, must be zero
			arch-size:       readUI32		;-- reserved, must be zero
			gptr-addr:       readUI32
			gptr-size:       readUI32
			TLS-addr:        readUI32
			TLS-size:        readUI32
			config-addr:     readUI32
			config-size:     readUI32
			b-import-addr:   readUI32
			b-import-size:   readUI32
			IAT-addr:        readUI32
			IAT-size:        readUI32
			d-import-addr:   readUI32
			d-import-size:   readUI32
			CLR-addr:        readUI32
			CLR-size:        readUI32
			reserved:        readUI32	;-- reserved, must be zero
			reserved2:       readUI32		;-- reserved, must be zero
		]
	]
	read-section-header: does [
		section-header: context [
			name:             trim/with as-string readBytes 8 #"^@"		;-- placeholder for an 8 bytes string
			virtual-size:     readUI32
			virtual-address:  readUI32
			raw-data-size:    readUI32
			raw-data-ptr:     readUI32
			relocations-ptr:  readUI32
			line-num-ptr:     readUI32
			relocations-nb:   readUI16
			line-num-nb:      readUI16
			flags:            readBytes 4
		]
		repend section-headers [section-header/name section-header]
		section-header
	]
	read-rsrc-table: func[
		section-header
	   /local
		num-name-entries
		num-id-entries
		rsrc-table ptr
	][
		ptr: index? inBuffer
		rsrc-table: context [
			flags:            readBytes 4
			timestamp:        readUI32
			major-version:    readUI16
			minor-version:    readUI16
			name-entries: copy []
			id-entries:   copy []
		]
		num-name-entries: readUI16
		num-id-entries:   readUI16
		loop num-name-entries [
			append rsrc-table/name-entries read-rsrc-name-entry
		]
		loop num-id-entries [
			append rsrc-table/id-entries read-rsrc-id-entry
		]

		foreach dir rsrc-table/id-entries [
			either dir/subdirectory [
				inBuffer: at head inBuffer (rsrc-ptr + (dir/rva and 2147483647))
				dir/subdirectory: read-rsrc-table section-header
			][
				inBuffer: at head inBuffer (rsrc-ptr + dir/rva)
				dir/data-entry: read-rsrc-data-entry section-header/virtual-address
			]
		]
		inBuffer: at head inBuffer ptr
		rsrc-table
	]
	read-rsrc-name-entry: does [
		context [
			name-rva:         readUI32
			rva:              readUI32
		]
	]
	read-rsrc-id-entry: does [
		context [
			id:               readUI32
			rva:              readUI32
			subdirectory:    0 <> (rva and -2147483648) ;high bit = 0
			data-entry:       none
		]
	]
	read-rsrc-data-entry: func[virtual-addres /local result] [
		result: context [
			data-rva: readUI32 - virtual-addres
			size:     readUI32
			codepage: readUI32
			reserved: readUI32
			data:     none
		]
		inBuffer: at head inBuffer (rsrc-ptr + result/data-rva)
		result/data: readBytes result/size
		result
	]
	parse: func[
		"Parses PE file and prints info about it's content"
		PE-file [file! url!]
		/local sh
	][
		clear section-headers
		open PE-file
		either #{4D5A} = readBytes 2 [
			probe readBytes 58
			probe PESignatureOffset: readUI32
			skipBytes PESignatureOffset - 64
			if #{50450000} <> probe readBytes 4 [
				print "PE Signature not found!"
				close PE-port
				return false
			]
			probe read-file-header
			probe file-header/opt-headers-size
			probe read-optional-header
			loop file-header/sections-nb [
				probe read-section-header
			]
			if any [
				sh: select section-headers ".rsrc"
				sh: select section-headers "rsrc"
			]	[
				print "RSRC"
				inBuffer: at head inBuffer rsrc-ptr: sh/raw-data-ptr + 1
				probe length? inBuffer
				probe read-rsrc-table sh
				;probe inBuffer

				;probe xxx: inBuffer
			]
		][
			print "Not MSDOS file"
			close PE-port
			return false
		]
		close PE-port
	]
]
echo %/d/t1.txt PE/parse %/d/hx.exe      echo none
echo %/d/t2.txt PE/parse %/d/hello2s.EXE echo none

max_stars_repo_name: Oldes/rs

max_issues_repo_name: Oldes/rs

max_forks_repo_name: Oldes/rs