text:
calcDeciles <- function(probabilities, labels) {
  probs.sorted <- probabilities[order(probabilities)]
  labels.sorted <- labels[order(probabilities)]

  deciles <- quantile(probs.sorted, probs=seq(0,1,.1))

  dec1.idx <- which(probs.sorted < deciles[2])
  dec2.idx <- which((probs.sorted < deciles[3]) & (probs.sorted >= deciles[2]))
  dec3.idx <- which((probs.sorted < deciles[4]) & (probs.sorted >= deciles[3]))
  dec4.idx <- which((probs.sorted < deciles[5]) & (probs.sorted >= deciles[4]))
  dec5.idx <- which((probs.sorted < deciles[6]) & (probs.sorted >= deciles[5]))
  dec6.idx <- which((probs.sorted < deciles[7]) & (probs.sorted >= deciles[6]))
  dec7.idx <- which((probs.sorted < deciles[8]) & (probs.sorted >= deciles[7]))
  dec8.idx <- which((probs.sorted < deciles[9]) & (probs.sorted >= deciles[8]))
  dec9.idx <- which((probs.sorted < deciles[10]) & (probs.sorted >= deciles[9]))
  dec10.idx <- which(probs.sorted >= deciles[10])

  #response
  dec1.resp <- probs.sorted[dec1.idx]
  dec2.resp <- probs.sorted[dec2.idx]
  dec3.resp <- probs.sorted[dec3.idx]
  dec4.resp <- probs.sorted[dec4.idx]
  dec5.resp <- probs.sorted[dec5.idx]
  dec6.resp <- probs.sorted[dec6.idx]
  dec7.resp <- probs.sorted[dec7.idx]
  dec8.resp <- probs.sorted[dec8.idx]
  dec9.resp <- probs.sorted[dec9.idx]
  dec10.resp <- probs.sorted[dec10.idx]

  dec1.labs <- labels.sorted[dec1.idx]
  dec2.labs <- labels.sorted[dec2.idx]
  dec3.labs <- labels.sorted[dec3.idx]
  dec4.labs <- labels.sorted[dec4.idx]
  dec5.labs <- labels.sorted[dec5.idx]
  dec6.labs <- labels.sorted[dec6.idx]
  dec7.labs <- labels.sorted[dec7.idx]
  dec8.labs <- labels.sorted[dec8.idx]
  dec9.labs <- labels.sorted[dec9.idx]
  dec10.labs <- labels.sorted[dec10.idx]

  predicted <- vector()
  observed <- vector()

  predicted <- 100* c(mean(dec1.resp),mean(dec2.resp),mean(dec3.resp),mean(dec4.resp),mean(dec5.resp),
                 mean(dec6.resp),mean(dec7.resp),mean(dec8.resp),mean(dec9.resp),mean(dec10.resp))

  dec1.rates <- 100* ((sum(dec1.labs))/(length(dec1.labs)))
  dec2.rates <- 100* ((sum(dec2.labs))/(length(dec2.labs)))
  dec3.rates <- 100* ((sum(dec3.labs))/(length(dec3.labs)))
  dec4.rates <- 100* ((sum(dec4.labs))/(length(dec4.labs)))
  dec5.rates <- 100* ((sum(dec5.labs))/(length(dec5.labs)))
  dec6.rates <- 100* ((sum(dec6.labs))/(length(dec6.labs)))
  dec7.rates <- 100* ((sum(dec7.labs))/(length(dec7.labs)))
  dec8.rates <- 100* ((sum(dec8.labs))/(length(dec8.labs)))
  dec9.rates <- 100* ((sum(dec9.labs))/(length(dec9.labs)))
  dec10.rates <- 100* ((sum(dec10.labs))/(length(dec10.labs)))

  observed <- c(dec1.rates,dec2.rates,dec3.rates,dec4.rates,dec5.rates,dec6.rates,dec7.rates,dec8.rates,dec9.rates,dec10.rates)

  calibration <- data.frame(predicted, observed)
  colnames(calibration) <- c('predicted rates', 'observed rates')
  rownames(calibration) <- c('decile 1', 'decile 2', 'decile 3', 'decile 4', 'decile 5',
                             'decile 6', 'decile 7', 'decile 8', 'decile 9', 'decile 10')
  list.res <- list()
  list.res[[1]] <- calibration
  list.res[[2]] <- probs.sorted
  list.res[[3]] <- labels.sorted
  return(list.res)
}

sortDeciles <- function(probabilities, labels, deciles) {
  probs.sorted <- probabilities[order(probabilities)]
  labels.sorted <- labels[order(probabilities)]

  #deciles <- quantile(probs.sorted, probs=seq(0,1,.1))

  dec1.idx <- which(probs.sorted < deciles[2])
  dec2.idx <- which((probs.sorted < deciles[3]) & (probs.sorted >= deciles[2]))
  dec3.idx <- which((probs.sorted < deciles[4]) & (probs.sorted >= deciles[3]))
  dec4.idx <- which((probs.sorted < deciles[5]) & (probs.sorted >= deciles[4]))
  dec5.idx <- which((probs.sorted < deciles[6]) & (probs.sorted >= deciles[5]))
  dec6.idx <- which((probs.sorted < deciles[7]) & (probs.sorted >= deciles[6]))
  dec7.idx <- which((probs.sorted < deciles[8]) & (probs.sorted >= deciles[7]))
  dec8.idx <- which((probs.sorted < deciles[9]) & (probs.sorted >= deciles[8]))
  dec9.idx <- which((probs.sorted < deciles[10]) & (probs.sorted >= deciles[9]))
  dec10.idx <- which(probs.sorted >= deciles[10])

  #response
  dec1.resp <- probs.sorted[dec1.idx]
  dec2.resp <- probs.sorted[dec2.idx]
  dec3.resp <- probs.sorted[dec3.idx]
  dec4.resp <- probs.sorted[dec4.idx]
  dec5.resp <- probs.sorted[dec5.idx]
  dec6.resp <- probs.sorted[dec6.idx]
  dec7.resp <- probs.sorted[dec7.idx]
  dec8.resp <- probs.sorted[dec8.idx]
  dec9.resp <- probs.sorted[dec9.idx]
  dec10.resp <- probs.sorted[dec10.idx]

  dec1.labs <- labels.sorted[dec1.idx]
  dec2.labs <- labels.sorted[dec2.idx]
  dec3.labs <- labels.sorted[dec3.idx]
  dec4.labs <- labels.sorted[dec4.idx]
  dec5.labs <- labels.sorted[dec5.idx]
  dec6.labs <- labels.sorted[dec6.idx]
  dec7.labs <- labels.sorted[dec7.idx]
  dec8.labs <- labels.sorted[dec8.idx]
  dec9.labs <- labels.sorted[dec9.idx]
  dec10.labs <- labels.sorted[dec10.idx]

  predicted <- vector()
  observed <- vector()
  std <- vector()
  lower <- vector()
  upper <- vector()

  predicted <- 100* c(mean(dec1.resp),mean(dec2.resp),mean(dec3.resp),mean(dec4.resp),mean(dec5.resp),
                      mean(dec6.resp),mean(dec7.resp),mean(dec8.resp),mean(dec9.resp),mean(dec10.resp))
  std <- 100* c(sd(dec1.resp),sd(dec2.resp),sd(dec3.resp),sd(dec4.resp),sd(dec5.resp),
                      sd(dec6.resp),sd(dec7.resp),sd(dec8.resp),sd(dec9.resp),sd(dec10.resp))
  lower <- 100* c(t.test(dec1.resp)$conf.int[1],t.test(dec2.resp)$conf.int[1],t.test(dec3.resp)$conf.int[1],t.test(dec4.resp)$conf.int[1],t.test(dec5.resp)$conf.int[1],
                t.test(dec6.resp)$conf.int[1],t.test(dec7.resp)$conf.int[1],t.test(dec8.resp)$conf.int[1],t.test(dec9.resp)$conf.int[1],t.test(dec10.resp)$conf.int[1])
  upper <- 100* c(t.test(dec1.resp)$conf.int[2],t.test(dec2.resp)$conf.int[2],t.test(dec3.resp)$conf.int[2],t.test(dec4.resp)$conf.int[2],t.test(dec5.resp)$conf.int[2],
                t.test(dec6.resp)$conf.int[2],t.test(dec7.resp)$conf.int[2],t.test(dec8.resp)$conf.int[2],t.test(dec9.resp)$conf.int[2],t.test(dec10.resp)$conf.int[2])
  dec1.rates <- 100* ((sum(dec1.labs))/(length(dec1.labs)))
  dec2.rates <- 100* ((sum(dec2.labs))/(length(dec2.labs)))
  dec3.rates <- 100* ((sum(dec3.labs))/(length(dec3.labs)))
  dec4.rates <- 100* ((sum(dec4.labs))/(length(dec4.labs)))
  dec5.rates <- 100* ((sum(dec5.labs))/(length(dec5.labs)))
  dec6.rates <- 100* ((sum(dec6.labs))/(length(dec6.labs)))
  dec7.rates <- 100* ((sum(dec7.labs))/(length(dec7.labs)))
  dec8.rates <- 100* ((sum(dec8.labs))/(length(dec8.labs)))
  dec9.rates <- 100* ((sum(dec9.labs))/(length(dec9.labs)))
  dec10.rates <- 100* ((sum(dec10.labs))/(length(dec10.labs)))

  observed <- c(dec1.rates,dec2.rates,dec3.rates,dec4.rates,dec5.rates,dec6.rates,dec7.rates,dec8.rates,dec9.rates,dec10.rates)

  calibration <- data.frame(predicted, std, lower, upper, observed)
  colnames(calibration) <- c('predicted rates', 'standard deviation', 'lower', 'upper', 'observed rates')
  rownames(calibration) <- c('decile 1', 'decile 2', 'decile 3', 'decile 4', 'decile 5',
                             'decile 6', 'decile 7', 'decile 8', 'decile 9', 'decile 10')
  list.res <- list()
  list.res[[1]] <- calibration
  list.res[[2]] <- probs.sorted
  list.res[[3]] <- labels.sorted
  return(list.res)
}

calcQuartiles <- function(probabilities, labels) {
  probs.sorted <- probabilities[order(probabilities)]
  labels.sorted <- labels[order(probabilities)]

  deciles <- quantile(probs.sorted, probs=seq(0,1,.25))

  dec1.idx <- which(probs.sorted < deciles[2])
  dec2.idx <- which((probs.sorted < deciles[3]) & (probs.sorted >= deciles[2]))
  dec3.idx <- which((probs.sorted < deciles[4]) & (probs.sorted >= deciles[3]))
  dec4.idx <- which(probs.sorted >= deciles[4])

  #response
  dec1.resp <- probs.sorted[dec1.idx]
  dec2.resp <- probs.sorted[dec2.idx]
  dec3.resp <- probs.sorted[dec3.idx]
  dec4.resp <- probs.sorted[dec4.idx]

  dec1.labs <- labels.sorted[dec1.idx]
  dec2.labs <- labels.sorted[dec2.idx]
  dec3.labs <- labels.sorted[dec3.idx]
  dec4.labs <- labels.sorted[dec4.idx]

  predicted <- vector()
  observed <- vector()

  predicted <- 100* c(mean(dec1.resp),mean(dec2.resp),mean(dec3.resp),mean(dec4.resp))

  dec1.rates <- 100* ((sum(dec1.labs))/(length(dec1.labs)))
  dec2.rates <- 100* ((sum(dec2.labs))/(length(dec2.labs)))
  dec3.rates <- 100* ((sum(dec3.labs))/(length(dec3.labs)))
  dec4.rates <- 100* ((sum(dec4.labs))/(length(dec4.labs)))

  observed <- c(dec1.rates,dec2.rates,dec3.rates,dec4.rates)

  calibration <- data.frame(predicted, observed)
  colnames(calibration) <- c('predicted rates', 'observed rates')
  rownames(calibration) <- c('decile 1', 'decile 2', 'decile 3', 'decile 4')
  list.res <- list()
  list.res[[1]] <- calibration
  list.res[[2]] <- probs.sorted
  list.res[[3]] <- labels.sorted
  return(list.res)
}

sortQuartiles <- function(probabilities, labels, deciles) {
  probs.sorted <- probabilities[order(probabilities)]
  labels.sorted <- labels[order(probabilities)]

  #deciles <- quantile(probs.sorted, probs=seq(0,1,.25))

  dec1.idx <- which(probs.sorted < deciles[2])
  dec2.idx <- which((probs.sorted < deciles[3]) & (probs.sorted >= deciles[2]))
  dec3.idx <- which((probs.sorted < deciles[4]) & (probs.sorted >= deciles[3]))
  dec4.idx <- which(probs.sorted >= deciles[4])

  #response
  dec1.resp <- probs.sorted[dec1.idx]
  dec2.resp <- probs.sorted[dec2.idx]
  dec3.resp <- probs.sorted[dec3.idx]
  dec4.resp <- probs.sorted[dec4.idx]

  dec1.labs <- labels.sorted[dec1.idx]
  dec2.labs <- labels.sorted[dec2.idx]
  dec3.labs <- labels.sorted[dec3.idx]
  dec4.labs <- labels.sorted[dec4.idx]

  predicted <- vector()
  observed <- vector()

  predicted <- 100* c(mean(dec1.resp),mean(dec2.resp),mean(dec3.resp),mean(dec4.resp))

  dec1.rates <- 100* ((sum(dec1.labs))/(length(dec1.labs)))
  dec2.rates <- 100* ((sum(dec2.labs))/(length(dec2.labs)))
  dec3.rates <- 100* ((sum(dec3.labs))/(length(dec3.labs)))
  dec4.rates <- 100* ((sum(dec4.labs))/(length(dec4.labs)))

  observed <- c(dec1.rates,dec2.rates,dec3.rates,dec4.rates)

  calibration <- data.frame(predicted, observed)
  colnames(calibration) <- c('predicted rates', 'observed rates')
  rownames(calibration) <- c('decile 1', 'decile 2', 'decile 3', 'decile 4')
  list.res <- list()
  list.res[[1]] <- calibration
  list.res[[2]] <- probs.sorted
  list.res[[3]] <- labels.sorted
  return(list.res)
}

quartPoints <- function(list.probs, num.folds) {
  point.1 <- vector()
  point.2 <- vector()
  point.3 <- vector()
  point.4 <- vector()
  point.5 <- vector()

  for(f in 1:num.folds) {
    deciles <- quantile(list.probs[[f]], probs=seq(0,1,.25))
    point.1[f] <- deciles[1]
    point.2[f] <- deciles[2]
    point.3[f] <- deciles[3]
    point.4[f] <- deciles[4]
    point.5[f] <- deciles[5]
  }
  list.points <- list()
  list.points[[1]] <- point.1
  list.points[[2]] <- point.2
  list.points[[3]] <- point.3
  list.points[[4]] <- point.4
  list.points[[5]] <- point.5

  return(list.points)
}

decPoints <- function(list.probs, num.folds) {
  point.1 <- vector()
  point.2 <- vector()
  point.3 <- vector()
  point.4 <- vector()
  point.5 <- vector()
  point.6 <- vector()
  point.7 <- vector()
  point.8 <- vector()
  point.9 <- vector()
  point.10 <- vector()
  point.11 <- vector()

  for(f in 1:num.folds) {
    deciles <- quantile(list.probs[[f]], probs=seq(0,1,.1))
    point.1[f] <- deciles[1]
    point.2[f] <- deciles[2]
    point.3[f] <- deciles[3]
    point.4[f] <- deciles[4]
    point.5[f] <- deciles[5]
    point.6[f] <- deciles[6]
    point.7[f] <- deciles[7]
    point.8[f] <- deciles[8]
    point.9[f] <- deciles[9]
    point.10[f] <- deciles[10]
    point.11[f] <- deciles[11]
  }
  list.points <- list()
  list.points[[1]] <- point.1
  list.points[[2]] <- point.2
  list.points[[3]] <- point.3
  list.points[[4]] <- point.4
  list.points[[5]] <- point.5
  list.points[[6]] <- point.6
  list.points[[7]] <- point.7
  list.points[[8]] <- point.8
  list.points[[9]] <- point.9
  list.points[[10]] <- point.10
  list.points[[11]] <- point.11

  return(list.points)
}


text:
library(ggplot2)
library(dplyr)
library(hrbrthemes)
numModels <- c(12, 25, 50, 75, 100, 150, 200, 250, 300, 400)
rfModels <- data.frame(n = numeric(), time.t = numeric())
xgbModels <- data.frame(n = numeric(), time.t = numeric())
# read results
for(i in 1:length(numModels)){
  rfModels[i,] <- readRDS(file = paste0("./test/t", numModels[i], ".rds"))
  xgbModels[i,] <- readRDS(file = paste0("./test/xgboost", numModels[i], ".rds"))
}
# convert to minutes
rfModels[c(8:10),2] <- rfModels[c(8:10),2]*60
xgbModels[c(1:2),2] <- xgbModels[c(1:2),2]/60
# plot rf time
p1 <- ggplot(aes(x = n, y = time.t), data = rfModels)+
        geom_line(color = "purple4", size = 1) +
        geom_point(shape=20, color="grey69", fill="#69b3a2", size=4) +
        xlab("Number of models") +
        ylab("Time (in minutes)") +
        ggtitle("Random Forest Classifier - Computation Time.")
p1 + theme_bw()
# save plot
ggsave(filename = "./test/rfPlot.png", width = 15, height = 15, units = "cm")
# plot xgboost time
p2 <- ggplot(aes(x = n, y = time.t), data = xgbModels) +
        geom_line(color = "purple4", size = 1) +
        geom_point(shape=20, color="grey69", fill="#69b3a2", size=4) +
        xlab("Number of models") +
        ylab("Time (in minutes)") +
        ggtitle("xGBoost Classifier - Computation Time.")
p2 + theme_bw()
# save plot
ggsave(filename = "./test/xgbPlot.png", width = 15, height = 15, units = "cm")

ext: r

max_stars_repo_name: BiodataAnalysisGroup/metaDIEA

max_issues_repo_name: BiodataAnalysisGroup/metaDIEA

max_forks_repo_name: BiodataAnalysisGroup/metaDIEA


text:
##############################################################################
# graph_and_data_setup_for_DCI.r

# Edited by: G Oldford
    # Last modified: August, 2020
    # Inputs:
    #  adj_matrix, passability, lengths
    # Output:
    #  "summary_table_natural.csv" with passability values of each pathway
    #   if there were NO artificial barriers (the natural passability of the
    #   riverscape)
    # "summary_table_all.csv" with passability values
    #   of the pathways in the given riverscape (artificial and natural barriers
    #   included).  This way we can have an idea of how much the artificial barriers
    #   are really affecting the DCI of the riverscape.
    #
    # Output example:
    #  "start","end","path2","barrier_id","pathway_pass","start_section_length","finish_section_length"
    #  "55152_s","55152_s","55152_s","NA",1,120.1,120.1
    #  "55152_s","62689_s","55152_s,sink,62689_s","55152,62689",0.25,120.1,452.2
    #
    # Notes:
    #  this function calls sum_fx for natural only and natural / anthro barriers
    #  other than graph.fx
    #
    # to-do:

    #
    # Notes from previous coding work:
    #dci.fxs.r assigns adj.matrix and passability before it runs the graph.and.data.setup.for.DCI.r
    #output: get the summary tables of all possible pathways with their passability values,
    # start and end sections, barrier id, length of start and end segments
    #calls on "sum.fx.r" function

graph_and_data_setup_for_DCI<-function(adj_matrix,
                                       passability,
                                       lengths){
    # checks if any barriers are flagged as natural - GO, 2020
    NB<-sum(passability$nat_barrier==TRUE)>0

    #NB = natural barriers
    if (NB==TRUE){
        # Old note:
        # if there are natural barriers in the system, we want to know what the
        # overall DCInp and DCInd (n= natural) is for the system so we can compare
        # it to the DCIap and DCIad (a= all (anthropogenic + natural barriers))
        natural_passability<-passability
        # Old note:
        # only want to take into consideration passability values for natural barriers.
        # so we want to change the passability values for anthropogenic barriers to 1
        natural_passability$Pass[natural_passability$nat_barrier==FALSE]<-1

        sum_table_n<-sum_fx(adj_matrix=adj_matrix,
                            passability=natural_passability,
                            lengths=lengths)
        write.table(sum_table_n,
                    "summary_table_natural.csv",
                    row.names=F,
                    sep=",")
    }
    # Old note:
    #if you put the else statement here, then it will only do the below commands
    # only if the above statement is false.  if you take the else statement out,
    # it will do the above commands only if it's true, but it will ALWAYS do the
    # commands below.

    sum_table<-sum_fx(adj_matrix=adj_matrix,
                      passability=passability,
                      lengths=lengths)

    write.table(sum_table,
                "summary_table_all.csv",
                row.names=F,
                sep=",")

    return(NB)
}


text:
#' @export
boot.strata <- function (data, nless = 0, nresamp = 1, method = c("RESCALE",   "BWR", "NAIVE")) {
    if (!inherits(data, "strata"))
        stop("Not a legitimate strata object, need to run Stratify() first")
    call <- match.call(expand = FALSE)
    method <- match.arg(method)
    res1 <- data
    res3 <- summary(res1)
    if (nresamp == 0)
        stop("No resampling done.")
    out <- matrix(0, nrow = nresamp + 1, ncol = 3, dimnames = list(c("Actual",
        1:nresamp), c("Mean", "Variance",'gini')))
    out[1, ] <- c(res3$yst, (res3$se.yst)^2,res3$gini)
    if (method == "RESCALE") {
        res1$nh <- res1$nh - nless
        for (i in 1:nresamp) {
            yhib <- sapply(res1$yhi, rescale.boot, nless)
            out[i + 1, ] <- c(sum(res1$Wh * as.vector(sapply(yhib,
                mean)), na.rm = TRUE), (sum((((res1$Nh * (res1$Nh -
                res1$nh))/sum(res1$Nh)^2) * (as.vector(sapply(yhib,
                var))))/res1$nh, na.rm = TRUE)))
        }

    }
    if (method == "BWR") {
        fh <- res1$nh/res1$Nh
        kh <- (res1$nh - 1)/(1 - fh)
        ph <- ((1/kh) - (1/ceiling(kh)))/((1/floor(kh)) - (1/ceiling(kh)))
        for (i in 1:nresamp) {
            yhib <- bwr.boot(res1$yhi, kh, ph, sample, replace = TRUE,
                simplify = FALSE)
            yhib[res1$nh == 1] <- res1$yhi[res1$nh == 1]
            nhws = sapply(yhib,FUN=function(x) sum(x>0))

            out[i + 1, ] <- c(sum(res1$Wh * as.vector(sapply(yhib,
                mean)), na.rm = TRUE), (sum((((res1$Nh * (res1$Nh -
                as.vector(sapply(yhib, length))))/sum(res1$Nh)^2) *
                (as.vector(sapply(yhib, var))))/as.vector(sapply(yhib,
                length)), na.rm = TRUE)),gini(x=as.vector(sapply(yhib,
                mean)),y=res1$Nh))
        }
    }
    if (method == "NAIVE") {
        for (i in 1:nresamp) {
            yhib <- sapply(res1$yhi, sample, replace = TRUE,
                simplify = FALSE)
            out[i + 1, ] <- c(sum(res1$Wh * as.vector(sapply(yhib,
                mean)), na.rm = TRUE), (sum((((res1$Nh * (res1$Nh -
                as.vector(sapply(yhib, length))))/sum(res1$Nh)^2) *
                (as.vector(sapply(yhib, var))))/as.vector(sapply(yhib,
                length)), na.rm = TRUE)))
        }
    }
    res <- list(orig.mean = out[1, 1], orig.var = out[1, 2],
        boot.means = out[c(2:(nresamp + 1)), 1], boot.vars = out[c(2:(nresamp +
            1)), 2],gini = out[c(2:(nresamp + 1)), 3], accel = accel.str(res1), call = call, method = method)
    oldClass(res) <- "boot"
    return(res)
}


text:
REBOL []

name: 'Event
source: %event/mod-event.c
includes: [%prep/extensions/event]

depends: compose [
    %event/t-event.c

    (switch system-config/os-base [
        'Windows [
            [%event/event-windows.c]
        ]
    ] else [
        [%event/event-posix.c]
    ])
]

libraries: try switch system-config/os-base [
    'Windows [
        ;
        ; Needed for SetTimer(), GetMessage(), etc.
        ;
        [%user32]
    ]
]

max_stars_repo_name: hostilefork/rebol

max_issues_repo_name: hostilefork/rebol

max_forks_repo_name: hostilefork/rebol


text:
    # Title: R Margins, single figure example
    # Date: 2008-02-04
    # Project: R graphics

    # Cheerfully borrowed from::
	# http://cran.r-project.org/doc/manuals/R-intro.html#Graphics

    #rm(list = ls())     # clear objects
    #graphics.off()      # close graphics windows

    # Generate data
    x = 0:10;
    y = 0:10;

    # Plot the data

    # - Specify the layout parameters before any plotting
    #   If you don't specify them before any plotting, the
    #   results will be inconsistent and misaligned.
    #
    # - oma stands for 'Outer Margin Area', or the total margin space that is outside
    #   of the standard plotting region (see graph)
    #
    # - The vector is ordered, the first value corresponding to the bottom. The entire
    #   array is c(bottom, left, top, right)
    #
    # - All of the alternatives are:
    #   - oma: Specify width of margins in number of lines
    #   - omi: Specify width of margins in inches
    #   - omd: Specify width of margins in 'device coordinates'
    #       - Device coordinates place (0,0) in the upper left and (1,1) in the
    #         lower right corner

    par(oma=c(3,3,3,3))  # all sides have 3 lines of space
    #par(omi=c(1,1,1,1)) # alternative, uncomment this and comment the previous line to try

    # - The mar command represents the figure margins. The vector is in the same ordering of
    #   the oma commands.
    #
    # - The default size is c(5,4,4,2) + 0.1, (equivalent to c(5.1,4.1,4.1,2.1)).
    #
    # - The axes tick marks will go in the first line of the left and bottom with the axis
    #   label going in the second line.
    #
    # - The title will fit in the third line on the top of the graph.
    #
    # - All of the alternatives are:
    #   - mar: Specify the margins of the figure in number of lines
    #   - mai: Specify the margins of the figure in number of inches

    par(mar=c(5,4,4,2) + 0.1)
    #par(mai=c(2,1.5,1.5,.5)) # alternative, uncomment this and comment the previous line

    # Plot
    #plot(x, y, type="n", xlab="X", ylab="Y")    # type="n" hides the points
	plot(rnorm(n=20),ylim=c(-2.0,2.0),xlim=c(1,20))  # simple plot so show how lines are oriented in plot area

    # Place text in the plot and color everything plot-related red
    text(5,5, "Plot", col="red", cex=2)
    text(5,4, "text(5,5, \"Plot\", col=\"red\", cex=2)", col="red", cex=1)
    box("plot", col="red")

    # Place text in the margins and label the margins, all in green
    mtext("Margins", side=3, line=2, cex=1.5, col="darkgreen")
    mtext("par(mar=c(5,4,4,2) + 0.1)", side=3, line=1, cex=1, col="darkgreen")

	mtext("Line 0", side=4, line=0, adj=1.0, cex=1, col="darkgreen")
	mtext("Line 1", side=4, line=1, adj=1.0, cex=1, col="darkgreen")

    mtext("Line 0", side=3, line=0, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 1", side=3, line=1, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 2", side=3, line=2, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 3", side=3, line=3, adj=1.0, cex=1, col="darkgreen")

    mtext("Line 0", side=2, line=0, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 1", side=2, line=1, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 2", side=2, line=2, adj=1.0, cex=1, col="darkgreen")
    mtext("Line 3", side=2, line=3, adj=1.0, cex=1, col="darkgreen")

	mtext("Line 0", side=1, line=0, adj=1.0, cex=1, col="darkgreen")
	mtext("Line 1", side=1, line=1, adj=1.0, cex=1, col="darkgreen")
	mtext("Line 2", side=1, line=2, adj=1.0, cex=1, col="darkgreen")
	mtext("Line 3", side=1, line=3, adj=1.0, cex=1, col="darkgreen")
	mtext("Line 4", side=1, line=4, adj=1.0, cex=1, col="darkgreen")

    box("figure", col="darkgreen")

    # Label the outer margin area and color it blue
    # Note the 'outer=TRUE' command moves us from the figure margins to the outer
    # margins.
    mtext("Outer Margin Area", side=3, line=1, cex=1.5, col="blue", outer=TRUE)  # default for series is side=1, not 3
    mtext("par(oma=c(3,3,3,3))", side=3, line=2, cex=1, col="blue", outer=TRUE)
    mtext("Line 0", side=3, line=0, adj=0.0, cex=1, col="blue", outer=TRUE)
    mtext("Line 1", side=3, line=1, adj=0.0, cex=1, col="blue", outer=TRUE)
    mtext("Line 2", side=3, line=2, adj=0.0, cex=1, col="blue", outer=TRUE)
    box("outer", col="blue")

max_stars_repo_name: f-tonini/Intro-Spatial-R

max_issues_repo_name: f-tonini/Intro-Spatial-R

max_forks_repo_name: f-tonini/Intro-Spatial-R


text:
library(ggplot2)

plot_df <- read.csv("C:/Users/LinaL/OneDrive/Desktop/Tutorial essays, prac/Computing sessions/question-3/data/average-rainfall.csv")

g <- ggplot(
  data = plot_df,
  mapping = aes(
    x = as.integer(month_num),
    y = average_rainfall,
    colour = location)) +
  geom_point() +
  geom_line() +
  labs(
    x = NULL,
    y = "Average rainfall (mm)",
    colour = "City",
    title = "Rainfall",
    subtitle = "Monthly average (1855--2015)") +
  scale_x_continuous(
    breaks = plot_df$month_num[seq(2,12,2)],
    labels = plot_df$month[seq(2,12,2)]) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = -45))

ggsave(filename = "C:/Users/LinaL/OneDrive/Desktop/Tutorial essays, prac/Computing sessions/question-3/out/result.png",
       plot = g,
       height = 10.5, width = 14.8,
       units = "cm")

sink(file = "C:/Users/LinaL/OneDrive/Desktop/Tutorial essays, prac/Computing sessions/question-3/out/package-versions-make-plot.txt")
sessionInfo()
sink()


text:
args <- commandArgs(TRUE)
help(args[1])

max_stars_repo_name: ktisha/TheRPlugin

max_issues_repo_name: ktisha/TheRPlugin

max_forks_repo_name: ktisha/TheRPlugin


text:
# Example : 1   Chapter : 8.1     Page No: 413
# Find movements , tensions, elongations of spring
K<-matrix(c(2,-1,0,-1,2,-1,0,-1,2),ncol=3)
K1<-solve(K)
f<-c(1,1,1)#since all mi=m
u<-K1%*%f
print("Movements are given by mg/c * u and u  is")
print(u)
A<-matrix(c(1,-1,0,0,0,1,-1,0,0,0,1,-1),ncol=3)
e<-A%*%u
print("Elongations are given by mg/c * e and e is")
print(e)
print("Tensions are given by mg * e and e is")
print(e)


text:
# #############################################################################
# SPI.r
# Author: Boris Orlowsky (ETH, Switzerland)
# Contributor: Martin Evaldsson (SMHI, Sweden)
# #############################################################################
# Description
#    For each month, the precipitation over the preceding TIMESCALE months, x,
#    is summed. Then a two-parameter Gamma distribution of cumulative
#    probability, Gamma_\alpha,\beta , is fitted to the strictly positive
#    TIMESCALE month sums, such that the probability of a non-zero
#    precipitation sum being below a certain value x corresponds to
#    Gamma_\alpha,\beta (x). We estimate shape parameter \alpha and scale
#    parameter \beta with a maximum likelihood approach. If the estimation
#    does not converge, \alpha and \beta are approximated using empirical
#    relations (Bordi et al., 2001). Accounting for TIMESCALE month periods
#    of no precipitation, occurring at a frequency q, the total cumulative
#    probability distribution of a precipitation sum below x, H(x), becomes
#    H(x) = q + (1 - q)*Gamma_\alpha,\beta (x). In the last step, a
#    precipitation sum x is assigned to its corresponding SPI value by
#    computing the quantile q_N(0,1) of the standard normal distribution at
#    probability H(x). The SPI of a precipitation sum x, thus, corresponds
#    to the quantile of the standard normal distribution which is assigned
#    by preserving the probability of the original precipitation sum, H(x).
#
# Required diag_script_info attributes (diagnostics specific)
#
# Optional diag_script_info attributes (diagnostic specific)
#
# Required variable_info attributes (variable specific)
#
# Optional variable_info attributes (variable specific)
#
# Caveats
#    No output of processed files to log-file yet!
#
# Modification history
#    20151113-A_laue_ax: added header
#
# ############################################################################

source('diag_scripts/aux/SPI/SPI_auxiliary.r')
source('interface_data/r.interface')
source(diag_script_cfg)
source('diag_scripts/lib/R/info_output.r')

## Do not print warnings
options(warn=-1)

var0 <- variables[1]
field_type0 <- field_types[1]

info_output(paste0("<<<<<<<< Entering ", diag_script), verbosity, 4)
info_output("+++++++++++++++++++++++++++++++++++++++++++++++++", verbosity, 1)
info_output(paste0("plot - ", diag_script, " (var: ", variables[1], ")"), verbosity, 1)
info_output("+++++++++++++++++++++++++++++++++++++++++++++++++", verbosity, 1)

library(tools)
diag_base = file_path_sans_ext(diag_script)

## Create working dirs if they do not exist
dir.create(plot_dir, showWarnings = FALSE)
dir.create(file.path(plot_dir, diag_base), showWarnings = FALSE)
dir.create(climo_dir, showWarnings = FALSE)

##
## Run it all
##
for (model_idx in c(1:length(models_name))) {
    fullpath_filename <- interface_get_fullpath(var0, field_type0, model_idx)
    filename <- interface_get_infile(var0, field_type0, model_idx)
    spi.info <- generate.spi(fullpath_filename=fullpath_filename,
                             filename=filename,
                             start_year=models_start_year[model_idx],
                             end_year=models_end_year[model_idx],
                             timescale=timescale,
                             begin.ref.year=begin.ref.year,
                             end.ref.year=end.ref.year,
                             model_idx)

    evaluate.spi(spi.info)

    ##
    ## A simple map of the trend
    ##
    requested_data_type <- "_trends-n-pvals_wrt_"
    in.filename <- processed_spi_info(timescale, requested_data_type,
                                      begin.ref.year, end.ref.year,
                                      filename)

    suppressPackageStartupMessages(require("ncdf4"))
    eval.nc <- nc_open(file.path(climo_dir, in.filename))
    lon <- ncvar_get(eval.nc, "Lon")
    lat <- ncvar_get(eval.nc, "Lat")
    trends <- ncvar_get(eval.nc, "Decadal_Trend")
    p.values <- ncvar_get(eval.nc, "p_Value")

    ## Check if lats are ascending
    if (lat[1] > lat[length(lat)]) {
        lat <- rev(lat)
        trends <- trends[, dim(trends)[2]:1, ]
        p.values <- p.values[, dim(p.values)[2]:1, ]
    }

    ## Same range for all seasons
    max <- max(abs(trends), na.rm = TRUE)

    ## Choose the world
    {
        if (min(lon) < -175)
            my.world <- "world"
        else my.world <- "world2"
    }

    suppressPackageStartupMessages(require("fields"))
    suppressPackageStartupMessages(require("maps"))

    ## Output figure filename
    n_diag_script <- nchar(diag_script)
    diag_script_base <- substr(diag_script, 0, n_diag_script - 2)
    description <- "_figure_"
    aux_info <- processed_spi_info(timescale, description,
                                   begin.ref.year, end.ref.year,
                                   "")
    figure_filename <- interface_get_figure_filename(diag_script_base,
                                                     var0,
                                                     field_type0,
                                                     aux_info,
                                                     model_idx)
    figure_filename <- paste(figure_filename, output_file_type, sep=".")
    figure_filename <- file.path(plot_dir, diag_base, figure_filename)

    ## Chose output format for figure
    if (tolower(output_file_type) == "png") {
        png(filename = figure_filename,
            width <- png_width,
            height <- png_height,
            units <- png_units,
            pointsize <- png_pointsize,
            bg <- png_bg)
    } else if (tolower(output_file_type) == "pdf") {
        pdf(file <- figure_filename)
    } else if (tolower(output_file_type) == "eps") {
        setEPS()
        postscript(figure_filename)
    }
    par(mfrow = c(2, 3))
    for (i in 1:5) {
        ## Plot trends
        image(lon, lat, trends[, , i], main = paste(seasons[i],
              "SPI trends [1/decade] (", models_name[model_idx], ")"),
               xlab = "lon", ylab = "lat",
               zlim = c(-spi_colorbar_max, spi_colorbar_max),
               col = my.colors(11))

        ## Add significance level as contour
        contour(lon, lat, p.values[, , i], levels = 0.1, add = TRUE)

        ## Brave new world
        map(my.world, add = TRUE, int = FALSE)
    }

    ## Add colorbar
    plot.new()
    image.plot(legend.only = TRUE, zlim = c(-spi_colorbar_max, spi_colorbar_max),
               col = my.colors(11),
        smallplot <- c(0.1, 0.2, 0.3, 0.7))

    invisible(dev.off())
}
info_output(paste0(">>>>>>>> Leaving ", diag_script), verbosity, 4)


text:
# Environment for linear regression

loadLibs = function () {
  print("Loading libs...")
  library(dplyr)
  library(ggplot2)
  library(QuantPsyc)
  library(car)
  print("Libs loaded!")
}


text:
context("read_vcf")

v <- system.file("extdata", "test.vcf.gz", package = "valr")
x <- read_vcf(v)

test_that("colnames are set in vcf df", {
  expect_true("chrom" %in% colnames(x))
  expect_true("start" %in% colnames(x))
  expect_true("end" %in% colnames(x))
})

test_that("chrom names are set correctly", {
  expect_true(all(stringr::str_detect(x$chrom, "^chr")))
})

max_stars_repo_name: jimhester/valr

max_issues_repo_name: jimhester/valr

max_forks_repo_name: jimhester/valr


text:
\name{generate_BB}
\alias{generate_BB}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
  Simulating networks from the Bianconi-Barabasi model}
\description{
  This function generates networks from the Bianconi-\enc{Barabási}{Barabasi} model. It is a \sQuote{preferential attachment with fitness} model. In this model, the preferential attachment function is linear, i.e. \eqn{A_k = k}, and node fitnesses are sampled from some probability distribution.
}
\usage{
generate_BB(N              = 1000   ,
            num_seed       = 2      ,
            multiple_node  = 1      ,
            m              = 1      ,
            mode_f         = "gamma",
            s              = 10     )
}

%- maybe also 'usage' for other objects documented here.
\arguments{
  The parameters can be divided into two groups.

The first group specifies basic properties of the network:
    \item{N}{
      Integer. Total number of nodes in the network (including the nodes in the seed graph). Default value is \code{1000}.
    }
  \item{num_seed}{
    Integer. The number of nodes of the seed graph (the initial state of the network). The seed graph is a cycle. Default value is \code{2}.
  }
  \item{multiple_node}{
    Positive integer. The number of new nodes at each time-step. Default value is \code{1}.
  }
\item{m}{
  Positive integer. The number of edges of each new node. Default value is \code{1}.
}

The final group of parameters specifies the distribution from which node fitnesses are generated:
    \item{mode_f}{
      String. Possible values:\code{"gamma"}, \code{"log_normal"} or \code{"power_law"}. This parameter indicates the true distribution for node fitness. \code{"gamma"} = gamma distribution, \code{"log_normal"} = log-normal distribution. \code{"power_law"} = power-law (pareto) distribution. Default value is "gamma".
    }
\item{s}{
Non-negative numeric. The inverse variance parameter. The mean of the distribution is kept at \eqn{1} and the variance is \eqn{1/s} (since node fitnesses are only meaningful up to scale). This is achieved by setting shape and rate parameters of the Gamma distribution to \eqn{s}; setting mean and standard deviation in log-scale of the log-normal distribution to \eqn{-1/2*log (1/s + 1)} and \eqn{(log (1/s + 1))^{0.5}}; and setting shape and scale parameters of the pareto distribution to \eqn{(s+1)^{0.5} + 1} and \eqn{(s+1)^{0.5}/((s+1)^{0.5} + 1)}. If \code{s} is \code{0}, all node fitnesses \eqn{\eta} are fixed at \code{1} (i.e., \enc{Barabási}{Barabasi}-Albert model). The default value is \code{10}.
}
}

\value{
  The output is a \code{PAFit_net} object, which is a List contains the following four fields:
    \item{graph}{a three-column matrix, where each row contains information of one edge, in the form of \code{(from_id, to_id, time_stamp)}. \code{from_id} is the id of the source, \code{to_id} is the id of the destination.}
    \item{type}{a string indicates whether the network is \code{"directed"} or \code{"undirected"}.}
    \item{PA}{a numeric vector contains the true PA function.}
   \item{fitness}{fitness values of nodes in the network. The name of each value is the ID of the node.}
}
\author{
  Thong Pham \email{thongphamthe@gmail.com}
}
\references{
  1. Bianconni, G. & \enc{Barabási}{Barabasi}, A. (2001). Competition and multiscaling in evolving networks. Europhys. Lett., 54, 436 (\doi{10.1209/epl/i2001-00260-6}).
}
\seealso{
  For subsequent estimation procedures, see \code{\link{get_statistics}}.

  For other functions to generate networks, see \code{\link{generate_net}}, \code{\link{generate_BA}}, \code{\link{generate_ER}} and \code{\link{generate_fit_only}}. }

\examples{
  library("PAFit")
  # generate a network from the BB model with alpha = 1, N = 100, m = 1
  # The inverse variance of the Gamma distribution of node fitnesses is s = 10
  net <- generate_BB(N = 100,m = 1,mode = 1, s = 10)
  str(net)
  plot(net)
}
\concept{Bianconi-Barabasi model}

ext: rd

max_stars_repo_name: thongphamthe/PAFit

max_issues_repo_name: thongphamthe/PAFit

max_forks_repo_name: thongphamthe/PAFit


text:
| pc = 0xc000 | a = 0xff | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc002 | a = 0xff | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x101] = 0xff |
| pc = 0xc005 | a = 0xff | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x1] = 0xff |
| pc = 0xc007 | a = 0x44 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc009 | a = 0x44 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x102] = 0x44 |
| pc = 0xc00c | a = 0x44 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x2] = 0x44 |
| pc = 0xc00e | a = 0x7c | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc010 | a = 0x7c | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x103] = 0x7c |
| pc = 0xc013 | a = 0x3 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc015 | a = 0x3 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x6] = 0x3 |
| pc = 0xc017 | a = 0x1 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc019 | a = 0x1 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x7] = 0x1 |
| pc = 0xc01b | a = 0x4 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc01d | a = 0x4 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x5] = 0x4 |
| pc = 0xc01f | a = 0x6 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc021 | a = 0x6 | x = 0x0 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 | MEM[0x4] = 0x6 |
| pc = 0xc023 | a = 0x6 | x = 0x1 | y = 0x0 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc025 | a = 0x6 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc027 | a = 0x57 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc029 | a = 0x56 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 00100000 |
| pc = 0xc02b | a = 0xa9 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x1] = 0xff |
| pc = 0xc02d | a = 0xed | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x2] = 0x44 |
| pc = 0xc02f | a = 0x12 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x101] = 0xff |
| pc = 0xc032 | a = 0x56 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x102] = 0x44 |
| pc = 0xc035 | a = 0x2a | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x103] = 0x7c |
| pc = 0xc038 | a = 0x56 | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x103] = 0x7c |
| pc = 0xc03a | a = 0x2a | x = 0x1 | y = 0x2 | sp = 0xff | p[NV-BDIZC] = 10100000 | MEM[0x103] = 0x7c |


text:
#' Translate an expression to sql.
#'
#' @section Base translation:
#' The base translator, `base_sql`,
#' provides custom mappings for `!` (to NOT), `&&` and `&` to
#' `AND`, `||` and `|` to `OR`, `^` to `POWER`,
#' \code{\%>\%} to \code{\%}, `ceiling` to `CEIL`, `mean` to
#' `AVG`, `var` to `VARIANCE`, `tolower` to `LOWER`,
#' `toupper` to `UPPER` and `nchar` to `LENGTH`.
#'
#' `c()` and `:` keep their usual R behaviour so you can easily create
#' vectors that are passed to sql.
#'
#' All other functions will be preserved as is. R's infix functions
#' (e.g. \code{\%like\%}) will be converted to their SQL equivalents
#' (e.g. `LIKE`). You can use this to access SQL string concatenation:
#' `||` is mapped to `OR`, but \code{\%||\%} is mapped to `||`.
#' To suppress this behaviour, and force errors immediately when dplyr doesn't
#' know how to translate a function it encounters, using set the
#' `dplyr.strict_sql` option to `TRUE`.
#'
#' You can also use [sql()] to insert a raw sql string.
#'
#' @section SQLite translation:
#' The SQLite variant currently only adds one additional function: a mapping
#' from `sd()` to the SQL aggregation function `STDEV`.
#'
#' @param ...,dots Expressions to translate. `translate_sql()`
#'   automatically quotes them for you.  `translate_sql_()` expects
#'   a list of already quoted objects.
#' @param con An optional database connection to control the details of
#'   the translation. The default, `NULL`, generates ANSI SQL.
#' @param vars Deprecated. Now call [partial_eval()] directly.
#' @param vars_group,vars_order,vars_frame Parameters used in the `OVER`
#'   expression of windowed functions.
#' @param window Use `FALSE` to suppress generation of the `OVER`
#'   statement used for window functions. This is necessary when generating
#'   SQL for a grouped summary.
#' @param context Use to carry information for special translation cases. For example, MS SQL needs a different conversion for is.na() in WHERE vs. SELECT clauses.  Expects a list.
#' @export
#' @examples
#' # Regular maths is translated in a very straightforward way
#' translate_sql(x + 1)
#' translate_sql(sin(x) + tan(y))
#'
#' # Note that all variable names are escaped
#' translate_sql(like == "x")
#' # In ANSI SQL: "" quotes variable _names_, '' quotes strings
#'
#' # Logical operators are converted to their sql equivalents
#' translate_sql(x < 5 & !(y >= 5))
#' # xor() doesn't have a direct SQL equivalent
#' translate_sql(xor(x, y))
#'
#' # If is translated into case when
#' translate_sql(if (x > 5) "big" else "small")
#'
#' # Infix functions are passed onto SQL with % removed
#' translate_sql(first %like% "Had%")
#' translate_sql(first %is% NULL)
#' translate_sql(first %in% c("John", "Roger", "Robert"))
#'
#' # And be careful if you really want integers
#' translate_sql(x == 1)
#' translate_sql(x == 1L)
#'
#' # If you have an already quoted object, use translate_sql_:
#' x <- quote(y + 1 / sin(t))
#' translate_sql_(list(x))
#'
#' # Windowed translation --------------------------------------------
#' # Known window functions automatically get OVER()
#' translate_sql(mpg > mean(mpg))
#'
#' # Suppress this with window = FALSE
#' translate_sql(mpg > mean(mpg), window = FALSE)
#'
#' # vars_group controls partition:
#' translate_sql(mpg > mean(mpg), vars_group = "cyl")
#'
#' # and vars_order controls ordering for those functions that need it
#' translate_sql(cumsum(mpg))
#' translate_sql(cumsum(mpg), vars_order = "mpg")
translate_sql <- function(...,
                          con = NULL,
                          vars = character(),
                          vars_group = NULL,
                          vars_order = NULL,
                          vars_frame = NULL,
                          window = TRUE) {
  if (!missing(vars)) {
    abort("`vars` is deprecated. Please use partial_eval() directly.")
  }

  translate_sql_(
    quos(...),
    con = con,
    vars_group = vars_group,
    vars_order = vars_order,
    vars_frame = vars_frame,
    window = window
  )
}

#' @export
#' @rdname translate_sql
translate_sql_ <- function(dots,
                           con = NULL,
                           vars_group = NULL,
                           vars_order = NULL,
                           vars_frame = NULL,
                           window = TRUE,
                           context = list()) {

  if (length(dots) == 0) {
    return(sql())
  }

  stopifnot(is.list(dots))

  if (!any(have_name(dots))) {
    names(dots) <- NULL
  }

  old_con <- set_current_con(con)
  on.exit(set_current_con(old_con), add = TRUE)

  if (length(context) > 0) {
    old_context <- set_current_context(context)
    on.exit(set_current_context(context), add = TRUE)
  }

  if (window) {
    old_group <- set_win_current_group(vars_group)
    on.exit(set_win_current_group(old_group), add = TRUE)

    old_order <- set_win_current_order(vars_order)
    on.exit(set_win_current_order(old_order), add = TRUE)

    old_frame <- set_win_current_frame(vars_frame)
    on.exit(set_win_current_frame(old_frame), add = TRUE)
  }

  variant <- sql_translate_env(con)
  pieces <- map(dots, function(x) {
    if (is_atomic(get_expr(x))) {
      escape(get_expr(x), con = con)
    } else {
      overscope <- sql_overscope(x, variant, con, window = window)
      on.exit(overscope_clean(overscope))
      escape(overscope_eval_next(overscope, x))
    }
  })

  sql(unlist(pieces))
}

sql_overscope <- function(expr, variant, con, window = FALSE,
                          strict = getOption("dplyr.strict_sql", FALSE)) {
  stopifnot(is.sql_variant(variant))

  # Default for unknown functions
  if (!strict) {
    unknown <- setdiff(all_calls(expr), names(variant))
    top_env <- ceply(unknown, default_op, parent = empty_env())
  } else {
    top_env <- child_env(NULL)
  }

  # Known R -> SQL functions
  special_calls <- copy_env(variant$scalar, parent = top_env)
  if (!window) {
    special_calls2 <- copy_env(variant$aggregate, parent = special_calls)
  } else {
    special_calls2 <- copy_env(variant$window, parent = special_calls)
  }

  # Existing symbols in expression
  names <- all_names(expr)
  name_env <- ceply(
    names, function(x) escape(ident(x), con = con),
    parent = special_calls2
  )

  # Known sql expressions
  symbol_env <- env_clone(base_symbols, parent = name_env)

  new_overscope(symbol_env, top_env)
}

is_infix_base <- function(x) {
  x %in% c("::", "$", "@", "^", "*", "/", "+", "-", ">", ">=", "<", "<=",
    "==", "!=", "!", "&", "&&", "|", "||", "~", "<-", "<<-")
}
is_infix_user <- function(x) {
  grepl("^%.*%$", x)
}

default_op <- function(x) {
  assert_that(is_string(x))

  if (is_infix_base(x)) {
    sql_infix(x)
  } else if (is_infix_user(x)) {
    x <- substr(x, 2, nchar(x) - 1)
    sql_infix(x)
  } else {
    sql_prefix(x)
  }
}

all_calls <- function(x) {
  if (!is.call(x)) return(NULL)

  fname <- as.character(x[[1]])
  unique(c(fname, unlist(lapply(x[-1], all_calls), use.names = FALSE)))
}

all_names <- function(x) {
  if (is.name(x)) return(as.character(x))
  if (!is.call(x)) return(NULL)

  unique(unlist(lapply(x[-1], all_names), use.names = FALSE))
}

# character vector -> environment
ceply <- function(x, f, ..., parent = parent.frame()) {
  if (length(x) == 0) return(new.env(parent = parent))
  l <- lapply(x, f, ...)
  names(l) <- x
  list2env(l, parent = parent)
}

max_stars_repo_name: tdsmith/dbplyr

max_issues_repo_name: tdsmith/dbplyr

max_forks_repo_name: tdsmith/dbplyr


text:
# Character input
 calories <- function(food) switch(food, apple=47, pizza=1500, stop("food not known"))
 calories("apple")  # 47
 calories("banana") # throws an error
# Numeric input
 alphabet <- function(number) switch(number, "a", "ab", "abc")
 alphabet(0) # null response
 alphabet(1) # "a"
 alphabet(2) # "ab"
 alphabet(3) # "abc"
 alphabet(4) # null response
# Note that no 'otherwise' option is allowed when the input is numeric.

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
r=359.96
https://sandbox.dams.library.ucdavis.edu/fcrepo/rest/collection/sherry-lehmann/catalogs/d7r592/media/images/d7r592-012/svc:tesseract/full/full/359.96/default.jpg Accept:application/hocr+xml

max_stars_repo_name: ucd-library/wine-price-extraction

max_issues_repo_name: ucd-library/wine-price-extraction

max_forks_repo_name: ucd-library/wine-price-extraction


text:
#' Likelihood ratio test for model fits
#'
#' @inheritParams  fit
#' @param  reduced  character vector of parameters used in the reduced model
#' @return  a data.frame with terms, their fitted values, and LRT result
lrt = function(x, y, reduced, weights=rep(1,nrow(x)), offset=0, control=list(),
                lower=-Inf, upper=Inf, verbose=FALSE) {
    full = fit(x, y, weights, offset, control, lower, upper, verbose=verbose)
    red = fit(x[,reduced, drop=FALSE], y, weights, offset, control,
              lower, upper, verbose=verbose)

    df = length(y) - ncol(x)
    stat = 2 * (red$value - full$value)
    pval = stats::pchisq(stat, df=df, lower.tail=FALSE)

    # only report the term difference between full and reduced model
    data.frame(statistic=stat, p.value=pval)
}

max_stars_repo_name: mschubert/boundie

max_issues_repo_name: mschubert/boundie

max_forks_repo_name: mschubert/boundie


text:
context("DESCRIPTION checks")

test_that("Parse DESCRIPTION file", {
  pkg <- pkg_desc("testNamespace")

  expect_identical(package_version("0.1"), pkg$get_version())
  expect_identical("testNamespace", pkg_name("testNamespace"))
})


text:
library(benchmarkMetrics)
source("../gitProjectExtras/gitBasedProjects/R/sourceAllLibs.r")
sourceAllLibs('../gitProjectExtras/gitBasedProjects/R/')

## uncomment on windows
#source('../gitProjectExtras/package_git2r.r')
#config(repository(), user.name="Douglas Kelley", user.email="douglas.i.kelley@gmail.com")

library(raster)
library(ncdf4)
library(rasterExtras)
library(rasterPlot)
library(plotrix)
library(mapdata)
library(mapplots)
library(ellipse)
library(vegan)
library(RcppEigen)
library(parallel)
library(snow)
library(reldist)
library(MASS)
data(worldHiresMapEnv)

sourceAllLibs('../rasterextrafuns/rasterPlotFunctions/R/')
sourceAllLibs('../rasterextrafuns/rasterExtras/R/')

setupProjectStructure(dirn = c("outputs", "data", "temp", "figs"))
ens_dir = paste(outputs_dir, '/ensembles_noTempC_negDist/', sep = '')
makeDir(ens_dir)

pr_datasets  = c('MSWEP', 'CRU', 'GPCC', 'CMORPH')
drought_vars = c('MADD', 'MADM', 'MConc', 'MDDM')

pr_datasets  = c('MSWEP')
drought_vars = c('MADD')

ensemble_no = round(seq(1,240000, length.out = 101) )
ensemble_no = round(seq(1,1000, length.out = 6) )

grab_cache_default = TRUE

sourceAllLibs('libs/')
sourceAllLibs('libs/LimTREE_r/')
sourceAllLibs('libs/plotting/')

extent = c(-180, 180, -30, 30)

fire_datasets = c("GFED_four_s")
pr_datasets  = c('MSWEP')#, 'CRU', 'GPCC', 'CMORPH')
drought_vars = c('MADD')#, 'MADM', 'MConc', 'MDDM')

paramFile = 'outputs/params-fullnewModel'

conFire_params = read.csv('../LimFIRE/data/params.csv')

plot_title = c('MAP', 'MAT', 'Stress', 'Exclusion')

Jules_fire_off_LU_off_fname = '../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-S2/'
Jules_fire_off_LU_on_fname = '../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-mort0/'
Jules_fire_on_LU_off_fname = '../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3/'
Jules_fire_on_LU_on_fname = '../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-2/'

Jules_fire_on_LU_on_fnames = c(Jules_fire_on_LU_on_fname,
			       "../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-mortHlf/",
			       "../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-mort10th/",
			       "../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-mort100th/",
			       "../fireMIPbenchmarking/data/ModelOutputs/JULES-INFERNO-SF3-fracHalf/")


text:
LC_Linear_Analysis_nlme_lmList <- function(Input_FileName) {

# R Functions to process dataset from Load Cell Array
# It calcuates slope and intercept.
#
# Author : Kuo-Hsine Chang, Ph.D. Dept. of Plant Agricuture, University of Guelph
# Contact: changks.github.io
# Version: 04162018

  # load library
    library(zoo)
    library(data.table)
    library(tidyverse)
    library(lubridate)
    library(stats)

  # read the csv file
    dt  <- fread(Input_FileName, header = T, sep = ',');

    df  <- data.frame(dt)

  # reshape data to get Time repeated for each variable
    df2 <- df %>%
    gather(var, value, -Time) %>%
         mutate(Time=ymd_hms(Time),                ## convert to date-time variable
         date=date(Time),                          ## date info only
         timeval=Time-floor_date(Time, "day"),      ## time since beginning of day
         datevar=interaction(date, var))            ## date/var combo

   # Linear Regression Analysis

   # na.omit: observations are removed if they contain any missing values
   # na.exclude: functions will pad residuals and predictions to the correct length by inserting NAs for omitted cases.
   # na.pass: keep all data, including NAs
   # na.fail: returns the object only if it contains no missing values

   # LM <- nlme::lmList(value ~ timeval | datevar, df2, na.action=na.pass)
    LM <- nlme::lmList(value ~ timeval | datevar, df2, na.action=na.exclude)

   #  LM <- nlme::lmList(value ~ timeval | datevar, df2, na.action = na.omit)

   # export to a csv file
     result <- data.frame(Subject = rownames(coef(LM)), coef(LM), check.names = FALSE)
     write.table(result, paste0(tools::file_path_sans_ext(Input_FileName), "_nlme_lmList.csv"), row.names = FALSE, col.names = TRUE, sep = ",")
}


text:
# Example : 1   Chapter : 7.3     Page No: 400
# Diagonaliazation of matrix
A<-matrix(c(0.5,-0.5,-0.5,0.5),ncol=2)
ev<-eigen(A)$values
D<-matrix(c(ev[1],0,0,ev[2]),ncol=2)
print("The diagonialized matrix")
print(D)


text:
library(devtools)
source_url("https://raw.githubusercontent.com/fdrtd/fdrtd/main/fdrtd/clients/r/fdrtd.r")

api <- Api("http://127.0.0.1:55502")
lst <- api$report()
print(lst)


text:
#' Coinbase API
#'
#' Make request to Coinbase API.
#'
#' @param method appropriate web method. These map to httr::GET, httr::POST, etc...
#' @param ... extra arguments passed to httr::GET, httr::POST, etc...
#' @keywords internal
#' @examples
#' \dontrun{
#' coinbase_get("prices/BTC-USD/spot", query = list(date = "2017-01-01"))
#' }
coinbase_api <- function(endpoint, method = c("get", "post", "put", "delete"), ...) {
  stopifnot(length(endpoint) == 1L)

  baseurl <- "https://api.coinbase.com"
  path = paste("/v2", endpoint, sep = "/")
  url <- httr::modify_url(baseurl, path = path)

  # get api key
  api_key <- Sys.getenv("COINBASE_KEY")

  # extract and validate arguments
  method <- match.arg(method)
  extra <- list(...)
  request_body <- get0("body", list2env(extra), inherits = FALSE, ifnotfound = "")

  # generate signature and timestamp
  timestamp <- as.character(as.integer(Sys.time()))
  signature <- access_signature(timestamp, path, request_body, method = method)

  # build headers
  # Please supply API version (YYYY-MM-DD) as CB-VERSION header
  headers <- httr::add_headers(
    "CB-ACCESS-TIMESTAMP" = timestamp,
    "CB-ACCESS-KEY" = api_key,
    "CB-ACCESS-SIGN" = signature,
    "CB-VERSION" = "2017-11-28",
    "Content-Type" = "application/json",
    "Accept" = "application/json",
    "User-Agent" = "github.com/zamorarr/rcoinbase")

  # get response
  method_fun <- switch(
    method,
    "get" = httr::GET,
    "post" = httr::POST,
    "put" = httr::PUT,
    "delete" = httr::DELETE
  )
  response <- method_fun(url, headers, ...)

  # check response type
  if (httr::http_type(response) != "application/json") {
    stop("API did not return json", call. = FALSE)
  }

  # parse content
  content <- httr::content(response, "text")
  json <- jsonlite::fromJSON(content, simplifyVector = FALSE)

  # check errors
  if (httr::http_error(response)) {
    stop(
      sprintf(
        "Coinbase API request failed [%s]\n%s\n",
        httr::status_code(response),
        json$errors[[1]]$message
        ),
      call. = FALSE
    )
  }

  # return object
  structure(
    list(
      content = json,
      path = path,
      response = response
    ),
    class = "coinbase_api"
    )
}

#' @rdname coinbase_api
#' @export
coinbase_get <- function(endpoint, ...) coinbase_api(endpoint, "get", ...)

#' @rdname coinbase_api
#' @export
coinbase_post <- function(endpoint, ...) coinbase_api(endpoint, "post", ..., encode = "json")

#' @rdname coinbase_api
#' @export
coinbase_put <- function(endpoint, ...) coinbase_api(endpoint, "put", ..., encode = "json")

#' @rdname coinbase_api
#' @export
coinbase_delete <- function(endpoint, ...) coinbase_api(endpoint, "delete", ..., encode = "json")

#' @export
print.coinbase_api <- function(x, ...) {
  cat("<Coinbase ", x$path, ">\n", sep = "")
  utils::str(x$content, 1)
  invisible(x)
}

#' Generate an api access signature
#'
#' The CB-ACCESS-SIGN header is generated by creating a sha256 HMAC using the
#' secret key on the prehash string timestamp + method + requestPath + body
#' (where + represents string concatenation). The timestamp value is the same as
#'  the CB-ACCESS-TIMESTAMP header.
#'
#' @param timestamp The CB-ACCESS-TIMESTAMP header MUST be number of seconds since Unix Epoch.
#' @param request_path path requesting
#' @param request_body The body is the request body string or omitted if there is no request body (typically for GET requests).
#' @param method The method should be UPPER CASE.
#'
#' @keywords internal
access_signature <- function(timestamp, request_path, request_body = "",
                             method = c("get", "post", "put", "delete") ) {
  method <- toupper(match.arg(method))

  if (any(nchar(request_body) > 0)) {
    body <- paste(jsonlite::toJSON(request_body, auto_unbox = TRUE), collapse = "\n")
  } else {
    body <- ""
  }

  msg <- paste0(timestamp, method, request_path, body)

  api_secret <- Sys.getenv("COINBASE_SECRET")
  digest::hmac(api_secret, msg, "sha256", raw = FALSE)
}

max_stars_repo_name: zamorarr/rcoinbase

max_issues_repo_name: zamorarr/rcoinbase

max_forks_repo_name: zamorarr/rcoinbase


text:
\name{tinamous}
\alias{tinamous}
\docType{data}
\title{Example phylota object (tinamous)}
\description{
  Final phylota object produced by running pipeline for the tinamous, Tinamiformes (8802).
}
\usage{
   data(tinamous)
}
\format{
  \code{tinamous} is generated from read_phylota().
}
\examples{
data(tinamous) # load object
str(tinamous)
}
\keyword{datasets}

max_stars_repo_name: hettling/phylotaR

max_issues_repo_name: haithamsghaier/phylotaR

max_forks_repo_name: haithamsghaier/phylotaR


text:
#' Consolidate a set of similarly-cohorted data sets
#'
#' @param data_list the list of similarly-cohorted data sets.
#' @param on the variable to consolidate on.
#' @importFrom purrr reduce
#' @export
consolidate <- function(data_list, on) {
  reduce(data_list, full_join, by = on)
}

max_stars_repo_name: kaneplusplus/forceps

max_issues_repo_name: kaneplusplus/foreceps

max_forks_repo_name: kaneplusplus/foreceps


text:
io = import('../io')

#' Creates a table of different identifiers and caches it
#'
#' @param dset  Ensembl data set, e.g. '{hsapiens,mmusculus}_gene_ensembl'
#' @param force  Re-generate table if it already exists
#' @return       A data.frame with gene and transcript-level information
ortho_table = function(dset="hsapiens_gene_ensembl", to="mmusculus", force=FALSE) {
    fname = sprintf("orthologue_table-%s-%s.RData", dset, to)
    cache = file.path(module_file(), "cache", fname)
    if (file.exists(cache) && !force)
        return(io$load(cache))

    mart = biomaRt::useMart(biomart="ensembl", dataset=dset)
    ids = c('external_gene_name', 'ensembl_gene_id',
            paste0(to, '_homolog_ensembl_gene'),
            paste0(to, '_homolog_associated_gene_name'))
    mapping = biomaRt::getBM(attributes=ids, mart=mart)
    for (col in colnames(mapping)) {
        is_empty = nchar(as.character(mapping[,col])) == 0
        mapping[[col]][is_empty] = NA
    }

    dir.create(dirname(cache), showWarnings=FALSE)
    save(mapping, file=cache)
    mapping
}

max_stars_repo_name: mschubert/ebits

max_issues_repo_name: mschubert/ebits

max_forks_repo_name: mschubert/ebits


text:
isaterm <- function()
{
  gui <- .Platform$GUI

  if (!isatty(stdin()))
    return(FALSE)
  # ban emacs: here and everywhere else in life
  else if (Sys.getenv("EMACS") == "t" || identical(getOption("STERM"), "iESS"))
    return(FALSE)
  else if (gui == "RTerm" || gui == "X11")
    return(TRUE)
  else if (gui == "unknown" && .Platform$OS.type == "unix" && Sys.getenv("RSTUDIO") != 1 && Sys.getenv("R_GUI_APP_VERSION") == "")
    return(TRUE) # I think?
  else
    return(FALSE)
}

hastcltk <- function()
{
  test <- tryCatch(requireNamespace("tcltk", quietly=TRUE), warning=identity)
  if (!is.logical(test))
    test <- FALSE

  test
}

max_stars_repo_name: cran/getPass

max_issues_repo_name: cran/getPass

max_forks_repo_name: cran/getPass


text:
r=0.23
https://sandbox.dams.library.ucdavis.edu/fcrepo/rest/collection/sherry-lehmann/catalogs/d77c7q/media/images/d77c7q-047/svc:tesseract/full/full/0.23/default.jpg Accept:application/hocr+xml


text:
#Analyses for Tutorial 04: PCMs

library(phytools)

#A quick simulation
mytree<- pbtree(n=50, scale=1) #one way to simulate a tree
  #note: in R, one has functions to simulate BD trees, random splits trees, etc. using different functions
plot(mytree)
X<-fastBM(tree=mytree) #simulates a continuous trait on phylogeny under Brownian motion
Y<-fastBM(tree=mytree)
cor(X,Y) #Not zero

#Now, what happens if we do this many times?
nsim=1000
trees<-pbtree(n=50, scale=1, nsim = nsim) #100 PB trees
X<-lapply(1:nsim, function(j) fastBM(tree=trees[[j]]))
Y<-lapply(1:nsim, function(j) fastBM(tree=trees[[j]]))
cor.raw<-unlist(lapply(1:nsim, function(j) cor(X[[j]],Y[[j]])))
hist(cor.raw, xlim=c(-1,1))

#Now condition on phylogeny:
X.pic<-lapply(1:nsim, function(j) pic(x = X[[j]],phy = trees[[j]]))
Y.pic<-lapply(1:nsim, function(j) pic(x = Y[[j]],phy = trees[[j]]))
cor.pic<-unlist(lapply(1:nsim, function(j) cor(X.pic[[j]],Y.pic[[j]])))
hist(cor.pic, xlim=c(-1,1)) #mean near zero, but LOOK at variation!

library(ape)
library(geiger)

#Here is a large time-dated molecular phylogeny (a chronogram):
ManderTree<-read.tree("TutorialData/Mander.tre",tree.names=T)
plot(ManderTree)

#Now read in data for species in the genus Hydromantes
Hyd.dat<-read.csv("TutorialData/HydromantesMns.csv", header=TRUE, row.names = 1)
Hyd.dat  #Notice we read in the first row as row.names. These MUST match the names in the phylogeny

#Now, prune the tree to match the data and vice-versa:
Hydro.new<-treedata(phy=ManderTree,data = Hyd.dat, warnings=FALSE)
plot(Hydro.new$phy)  #We have matched data to the tree!

#### 1: Phylogenetically Independent Contrasts
#Read in data, phylogeny, and match/prune them to one another
ManderTree<-read.tree("TutorialData/Mander.tre",tree.names=T)
   plot(ManderTree)
Mander.dat<-read.csv("TutorialData/PlethodonMns.csv", header=TRUE, row.names = 1)
Gps<-Mander.dat[,12]; names(Gps)<-row.names(Mander.dat)
Mander.dat<-Mander.dat[,-12]  #must remove prior to 'treedata', which does not like factors

Pleth.new<-treedata(phy=ManderTree,data = Mander.dat, warnings=FALSE)
   plot(Pleth.new$phy)
SVL.new<-Pleth.new$data[,1] #grabs appropriate column from data matrix after treedata
HL.new<-Pleth.new$data[,3]
  match(names(Gps),row.names(Pleth.new$data))  #double check that the factor matches

#Phylogenetically-naive analyses
plot(SVL.new,HL.new)
cor.test(SVL.new,HL.new)  #correlation of two traits: pretty high
anova(lm(HL.new~SVL.new))

#Generate PICs and test while conditioning on phylogeny
X.pic<-pic(x = SVL.new, phy = Pleth.new$phy)
Y.pic<-pic(x = HL.new, phy = Pleth.new$phy)
cor.test(X.pic,Y.pic)
plot(X.pic,Y.pic)
anova(lm(Y.pic~X.pic+0))   #regression through the origin (b/c order of taxa for contrast is arbitrary [can 'spin' on node])
summary(lm(Y.pic~X.pic+0))   #coefficients of the model

#### 2: Phylogenetic Generalized Least Squares
#To perform PGLS in R, we must first estimate the phylogenetic covariance matrix V:
V<-corBrownian(phy=Pleth.new$phy)

#Now we run the analysis:
library(nlme)
bm.gls<-gls(HL.new~SVL.new,correlation = V, data=data.frame(SVL.new, HL.new))
summary(bm.gls)
anova(bm.gls)

#### 3: Phylogenetic Transform
library(geomorph)
gdf<-geomorph.data.frame(SVL=SVL.new,HL=HL.new,Gps=Gps, tree=Pleth.new$phy)
res.PhyT<-procD.pgls(HL~SVL, data = gdf, phy = Pleth.new$phy, print.progress = FALSE)
summary(res.PhyT)    #Same F-ratio as before (note: p-value found using permutation in this case)
res.PhyT$pgls.coefficients  #Same parameter estimates as before

#### Phylogenetic ANOVA
#Phylogenetic anova
aov.phylo(SVL.new~Gps, phy = Pleth.new$phy)
procD.pgls(SVL~Gps, data = gdf, phy = Pleth.new$phy,print.progress = FALSE)
anova(gls(SVL.new~Gps,correlation = corBrownian(phy=Pleth.new$phy), data=data.frame(SVL.new, Gps)))  #identical to Phylo.transform

####Phylogenetic Signal
#Phylogenetic Signal
phylosig(tree=Pleth.new$phy, x=SVL.new, method="K", test=T, nsim=1000)  #phytools
  res<-physignal(A=SVL.new,phy = Pleth.new$phy, print.progress = FALSE)   #geomorph
summary(res)  #the same, but the latter method (from geomorph) can accomodate multivariate data
plot(res)

####Additional Simulation Approaches

#Simulate correlated data
library(MASS)
R<-matrix(0.7, nrow=2,ncol=2); diag(R)<-1
dat.sim<-mvrnorm(n=10000,Sigma = R,mu = c(0,0))
plot(dat.sim, asp=1)

#Simulate BM correlated data on phylogeny
mytree<- pbtree(n=100, scale=1) #one way to simulate a tree
R<-matrix(0.7, nrow=2,ncol=2); diag(R)<-1
dat.BM<-sim.char(phy = mytree,par = R,nsim = 1)[,,1]
plot(dat.BM, asp=1)

##Practical Questions

#To apply what you have learned, please do the following:

#1: Simulate a phylogeny (or read one into R)
#2: Simultate two continuous traits with some known correlation between them
#3: Estimate the evolutionary association of the traits, conditioned on the phylogeny, using one of the methods you have learned

#Upload your code to the assignment center on Canvas


text:
# Decision Lens API
#
# No description provided (generated by Swagger Codegen https://github.com/swagger-api/swagger-codegen)
#
# OpenAPI spec version: 1.0
#
# Generated by: https://github.com/swagger-api/swagger-codegen.git

#' ValueCount Class
#'
#' @field value
#' @field count
#'
#' @importFrom R6 R6Class
#' @importFrom jsonlite fromJSON toJSON
#' @export
ValueCount <- R6::R6Class(
  'ValueCount',
  public = list(
    `value` = NULL,
    `count` = NULL,
    initialize = function(`value`, `count`){
      if (!missing(`value`)) {
        stopifnot(is.character(`value`), length(`value`) == 1)
        self$`value` <- `value`
      }
      if (!missing(`count`)) {
        stopifnot(is.numeric(`count`), length(`count`) == 1)
        self$`count` <- `count`
      }
    },
    toJSON = function() {
      ValueCountObject <- list()
      if (!is.null(self$`value`)) {
        ValueCountObject[['value']] <- self$`value`
      }
      if (!is.null(self$`count`)) {
        ValueCountObject[['count']] <- self$`count`
      }

      ValueCountObject
    },
    fromJSON = function(ValueCountJson) {
      ValueCountObject <- dlensFromJSON(ValueCountJson)
      if (!is.null(ValueCountObject$`value`)) {
        self$`value` <- ValueCountObject$`value`
      }
      if (!is.null(ValueCountObject$`count`)) {
        self$`count` <- ValueCountObject$`count`
      }
    },
    toJSONString = function() {
       sprintf(
        '{
           "value": %s,
           "count": %d
        }',
        self$`value`,
        self$`count`
      )
    },
    fromJSONString = function(ValueCountJson) {
      ValueCountObject <- dlensFromJSON(ValueCountJson)
      self$`value` <- ValueCountObject$`value`
      self$`count` <- ValueCountObject$`count`
    }
  )
)

max_stars_repo_name: dlens/dlxapi

max_issues_repo_name: dlens/dlxapi

max_forks_repo_name: dlens/dlxapi


text:
A -> D (3)
C -> D (1)
D -> E (2)
E -> B (1)
E -> G (1)
D -> F (6)

max_stars_repo_name: TXCodeDancer/OpenSource

max_issues_repo_name: TXCodeDancer/OpenSource

max_forks_repo_name: TXCodeDancer/OpenSource


text:
library(covid19AgeModel)

args = list()
args$indir = "~/git/R0t/covid19AgeModel/inst"

args$file_age_ifr <- file.path(args$indir,'data','ifr-by-age-prior_Verity_raw_200429.csv')

ifr.by.age <- as_tibble(read.csv(args$file_age_ifr))
stopifnot( c('age','ifr_mean')%in%colnames(ifr.by.age) )
ifr.by.age <- ifr.by.age %>%
  select(age, ifr_mean, ifr_cl, ifr_cu)

#
# Fit smoth spline to find ifr by 1 year age band
library(splines)
x = seq(mean(c(0,9)),mean(c(80,89)),10)
y.mean = unique(ifr.by.age$ifr_mean)
y.cl = unique(ifr.by.age$ifr_cl)
y.cu = unique(ifr.by.age$ifr_cu)
fit.smooth.mean = smooth.spline(x = x, y = y.mean, cv = TRUE)
fit.smooth.cl = smooth.spline(x = x, y = y.cl, cv = TRUE)
fit.smooth.cu = smooth.spline(x = x, y = y.cu, cv = TRUE)
ifr.by.age$ifr_mean = predict(fit.smooth.mean, x = 0:100)$y
ifr.by.age$ifr_cl = predict(fit.smooth.cl, x = 0:100)$y
ifr.by.age$ifr_cu = predict(fit.smooth.cu, x = 0:100)$y + 2*predict(fit.smooth.mean, x = 0:100)$y

#
# Aggregate by 5 year age band
GFNAME_us_population = file.path(args$indir,"data","us_population_withnyc.rds")
ifr.by.age = as.data.table(read_ifr_data_by_age(path_to_file_pop_us = GFNAME_us_population, ifr.by.age = ifr.by.age, max_age = 95))

#
# Find hyperparameters log normal distribution
library(lognorm)
ifr.by.age[, ln_mu := getParmsLognormForLowerAndUpper(lower = ifr_cl, upper = ifr_cu, sigmaFac=2)[1], by = c("age")]
ifr.by.age[, ln_sd := getParmsLognormForLowerAndUpper(lower = ifr_cl, upper = ifr_cu, sigmaFac=2)[2], by = c("age")]

#
# Find hyperparameters beta distribution
library("prevalence")
ifr.by.age = as.data.table(ifr.by.age)
ifr.by.age[, ifr_beta_alpha := betaExpert(ifr_mean, ifr_cl, ifr_cu, p = 0.95, method = "mean")$alpha, by = c("age")]
ifr.by.age[, ifr_beta_beta := betaExpert(ifr_mean, ifr_cl, ifr_cu, p = 0.95, method = "mean")$beta, by = c("age")]

#
# clean
setnames(ifr.by.age, "age", "age_cat")
ifr.by.age$age_cat_label = c(paste0(seq(0, 80, 5), "-", seq(4,84,5)), "85+")

write.csv(ifr.by.age, file =file.path(args$indir,'data','ifr-by-age-prior_Verity_200624.csv'), row.names = F)


text:
#' Extract the information from the simulation data frame to analyse the individual conditional expectation
#'
#' @param allsim dataset with all simulations values
#' @param dataset dataset with all variables
#' @param dr a vector with dose response values
#' @param squem squeme of the values of the prediction values
#' @param remove_extrem boolean parameter to remove the extrem values
#' @return a data frame with interactions
#' @examples
#' data(expose_data)
#' data(simu)
#' data(gen)
#' delta=c(1,0)
#' seku <- seq(0,1,0.05)
#' Exposures<- c('Var1','Var2','Var3','Var4','Var5')
#' summary_table_lines <- gen[[2]]
#' ice_res <- ice(allsim = simu[[1]], dataset = expose_data, dr = seku,
#' squem = summary_table_lines, remove_extrem = FALSE)
#' @export

ice <- function(allsim, dataset, dr = seq(0, 1, 0.1), squem, remove_extrem = FALSE) {
  dataset <- data.frame(dataset)
    ext <- squem[grep("DR_", squem$Group), ]

    hm_pol <- unique(ext$Case)

    n <- nrow(dataset)

    ext2 <- ext[ext$Case == hm_pol[1], ]

    df <- data.frame(allsim)

    ini <- ext2$From[1]
    fin <- ext2$To[nrow(ext2)]

    allpcb <- data.frame(df[ini:fin, 1])

    allpcb$id <- rep(1:nrow(dataset), length(dr))

    mm <- NULL
    len2 <- length(dr) - 1
    for (g in 0:len2) {
        mm <- c(mm, rep(g, n))
    }

    allpcb$X <- mm
    names(allpcb)[1] <- paste0(hm_pol[1], "_pred")

    if (length(hm_pol) > 1) {
        hm_pol2 <- hm_pol[-1]
        for (h in hm_pol2) {
            # extracting again for each pollutant
            ext2 <- ext[ext$Case == h, ]

            df <- data.frame(allsim)

            ini <- ext2$From[1]
            fin <- ext2$To[nrow(ext2)]

            allpcb2 <- data.frame(df[ini:fin, 1])

            allpcb2$id <- rep(1:nrow(dataset), length(dr))

            allpcb2$X <- mm
            names(allpcb2)[1] <- paste0(h, "_pred")

            allpcb <- merge(allpcb, allpcb2, by = c("id", "X"))

        }
    }
    names(dataset)[1] <- "id"
    ice <- merge(dataset, allpcb, by = "id")

    seku_breaks <- as.character(sort(as.numeric(as.character(unique(ice$X)))))

    ice$X <- factor(ice$X, levels = seku_breaks)

    if (remove_extrem) {
        sekut <- dr[c(1, length(dr))]
        ice <- ice[!ice$X %in% sekut, ]
    }

    return(ice)
}

max_stars_repo_name: itamuria/expose

max_issues_repo_name: itamuria/expose

max_forks_repo_name: itamuria/expose


text:
# Convert to coordinate system specified by EPSG code ----------------------------------------------------------
convert_coords = function(lat, long, from = CRS("+init=epsg:4326"), to) {
  data = data.frame(long = long, lat = lat)
  coordinates(data) = ~ long + lat
  proj4string(data) = from
  #
  xy = data.frame(sp::spTransform(data, to))
  colnames(xy) = c("x", "y")
  return(unlist(xy))
}
# Calculate route points as line ------------------------------------------
points2line_trajectory = function(p) {
  c = st_coordinates(p)
  i = seq(nrow(p) - 2)
  l = purrr::map(i, ~ sf::st_linestring(c[.x:(.x + 1),]))
  s = purrr::map_dbl(i, function(x) {
    geosphere::distHaversine(c[x,], c[(x + 1),]) /
      as.numeric(p$time[x + 1] - p$time[x])
  })
  lfc = sf::st_sfc(l)
  a = seq(length(lfc)) + 1 # sequence to subset
  p_data = cbind(sf::st_set_geometry(p[a,], NULL), s)
  sf::st_sf(p_data, geometry = lfc)
}

max_stars_repo_name: natesheehan/flashman-tracks-and-trails

max_issues_repo_name: natesheehan/flashman-tracks-and-trails

max_forks_repo_name: natesheehan/flashman-tracks-and-trails


text:
read.admb <-
function(ifile)
{
	ret=read.fit(ifile)

	fn=paste(ifile,'.rep', sep='')
	A=read.rep(fn)
	A$fit=ret

	pfn=paste(ifile,'.psv',sep='')
	if(file.exists(pfn))
		A$post.samp=read.psv(pfn)

	return(A)
}

read.fit <-
function(ifile)
{
	# __Example:
	#	file <-("~/admb/simple")
	#	A <- reptoRlist(file)
	#	Note there is no extension on the file name.

	## The following is a contribution from:
	## Anders Nielsen that reads the par & cor files.
	ret<-list()
	parfile<-as.numeric(scan(paste(ifile,'.par', sep=''),
	 what='', n=16, quiet=TRUE)[c(6,11,16)])
	ret$nopar<-as.integer(parfile[1])
	ret$nlogl<-parfile[2]
	ret$maxgrad<-parfile[3]
	file<-paste(ifile,'.cor', sep='')
	lin<-readLines(file)
	ret$npar<-length(lin)-2
	ret$logDetHess<-as.numeric(strsplit(lin[1], '=')[[1]][2])
	sublin<-lapply(strsplit(lin[1:ret$npar+2], ' '),function(x)x[x!=''])
	ret$names<-unlist(lapply(sublin,function(x)x[2]))
	ret$est<-as.numeric(unlist(lapply(sublin,function(x)x[3])))
	ret$std<-as.numeric(unlist(lapply(sublin,function(x)x[4])))
	ret$cor<-matrix(NA, ret$npar, ret$npar)
	corvec<-unlist(sapply(1:length(sublin), function(i)sublin[[i]][5:(4+i)]))
	ret$cor[upper.tri(ret$cor, diag=TRUE)]<-as.numeric(corvec)
	ret$cor[lower.tri(ret$cor)] <- t(ret$cor)[lower.tri(ret$cor)]
	ret$cov<-ret$cor*(ret$std%o%ret$std)
	return(ret)
}

read.rep <-
function(fn)
{
	# The following reads a report file
	# Then the 'A' object contains a list structure
	# with all the elemements in the report file.
	# In the REPORT_SECTION of the AMDB template use
	# the following format to output objects:
	#  	report<<"object \n"<<object<<endl;
	#
	# The part in quotations becomes the list name.
	# Created By Steven Martell
	options(warn=-1)  #Suppress the NA message in the coercion to double

	ifile=scan(fn,what="character",flush=TRUE,blank.lines.skip=FALSE,quiet=TRUE)
	idx=sapply(as.double(ifile),is.na)
	vnam=ifile[idx] #list names
	nv=length(vnam) #number of objects
	A=list()
	ir=0
	for(i in 1:nv)
	{
		ir=match(vnam[i],ifile)
		if(i!=nv) irr=match(vnam[i+1],ifile) else irr=length(ifile)+1 #next row
		dum=NA
		if(irr-ir==2) dum=as.double(scan(fn,skip=ir,nlines=1,quiet=TRUE,what=""))
		if(irr-ir>2) dum=as.matrix(read.table(fn,skip=ir,nrow=irr-ir-1,fill=TRUE))

		if(is.numeric(dum))#Logical test to ensure dealing with numbers
		{
			A[[vnam[i]]]=dum
		}
	}
	options(warn=0)

	return(A)
}

read.psv <-
function(fn, nsamples=10000)
{
	#This function reads the binary output from ADMB
	#-mcsave command line option.
	#fn = paste(ifile,'.psv',sep='')
	filen <- file(fn, "rb")
	nopar <- readBin(filen, what = integer(), n = 1)
	mcmc <- readBin(filen, what = numeric(), n = nopar * nsamples)
	mcmc <- matrix(mcmc, byrow = TRUE, ncol = nopar)
	close(filen)
	return(mcmc)
}


text:
# Watson Translate
#
# IBM SPSS Modeler Extension
#
# Install function for packages
packages <- function(x){
  x <- as.character(match.call()[[2]])
  if (!require(x,character.only=TRUE)){
    install.packages(pkgs=x,repos="http://cran.r-project.org")
    require(x,character.only=TRUE)
  }
}

packages(httr)
packages(XML)

# This function is used to generate automatically the dataModel for SPSS Modeler
getMetaData <- function (data) {
  if( is.null(dim(data)))
    stop("Invalid data received: not a data.frame")
  if (dim(data)[1]<=0) {
    print("Warning : modelerData has no line, all fieldStorage fields set to strings")
    getStorage <- function(x){return("string")}
  } else {
    getStorage <- function(x) {
      x <- unlist(x)
      res <- NULL
      #if x is a factor, typeof will return an integer so we treat the case on the side
      if(is.factor(x)) {
        res <- "string"
      } else {
        res <- switch(typeof(x),
                      integer="integer",
                      double = "real",
                      "string")
      }
      return (res)
    }
  }
  col = vector("list", dim(data)[2])
  for (i in 1:dim(data)[2]) {
    col[[i]] <- c(fieldName=names(data[i]),
                  fieldLabel="",
                  fieldStorage=getStorage(data[i]),
                  fieldMeasure="",
                  fieldFormat="",
                  fieldRole="")
  }
  mdm<-do.call(cbind,col)
  mdm<-data.frame(mdm)
  return(mdm)
}

version <- R.Version()$major #Get R version for Post call
print(version)

data <- data.frame() #Data frame to be populated with results
u <- "%%user%%"   #Username and Password provided in CDB
p <- "%%password%%"
source_language =%%source%%  #From CDB Combo Boxes
  target_language =%%target%%
  base <- "https://gateway.watsonplatform.net/language-translator/api/v2/translate?source="

if (source_language != target_language){
  for(i in 1:nrow(modelerData)) {
    text <-modelerData$%%text%%[i] #Load text
    if (!(is.element("en", c(source_language, target_language)))){
      language<-"Unable to convert non-english to non-english at this time" }
    else if(source_language != target_language) {
      full_url <- paste(base,source_language,"&target=", target_language, sep = '')
      pars <- list(txt=iconv(text,to="UTF-8"))
      if (version == 2){
        r<-POST(full_url,authenticate(u, p), body = pars)
      } else{
        r<-POST(full_url,authenticate(u, p),add_headers(Accept = "text/plain"), body = pars)
      }
      stop_for_status(r)
      language<-content(r)
    }
    data[i,1] <- language
  }

  modelerData <-cbind(modelerData,data)
  colnames(modelerData)[ncol(modelerData)] <- "Translated Text" #Name new column
  modelerDataModel <- getMetaData(modelerData)
} else {
  print("Choose different languages for source and target values!")
}

max_stars_repo_path: Source code/script.r

max_stars_repo_name: IBMPredictiveAnalytics/Watson-Translate

max_issues_repo_path: Source code/script.r

max_issues_repo_name: IBMPredictiveAnalytics/Watson-Translate

max_forks_repo_path: Source code/script.r

max_forks_repo_name: IBMPredictiveAnalytics/Watson-Translate


text:
#'vms 2 viirs
#'
#'@param vmstooviirshp is the output of the 'vms2viirs2shp' function.
#'@param viirsfile is a geospatial dataframe with viirs data obtained from https://eogdata.mines.edu/vbd/
#'@param vmsfile is a geospatial dataframe with VMS data over a certain period of time obtained via KKP (the Fisheries Ministrie of Indonesia)
#'@param bufferdistance is the maximum distance the VMS and VIIRS data are allowed to be apart
#'@return a geopsatial dataframe with all detected vessels through VIIRS satellite that could be connected to VMS vessel data.
#'@export
vms2viirsanalysis <- function(vmstoviirshp,viirsfile,vmsfile,bufferdistance) {
  dagen <- unique(vmstoviirshp$DAY)
  viirsfile$Date_Mscan <- as.POSIXct(viirsfile$Date_Mscan)
  viirsfile$day <- as.Date(viirsfile$Date_Mscan)
  vmsfile$PING_TIME <- as.POSIXct(vmsfile$PING_TIME)
  vmsfile$day <- as.Date(vmsfile$PING_TIME)
  outfile <- viirsfile[0,]
  outfile <- as.data.frame(outfile)
  for (i in dagen) {
    dayn <- i
    vmstoviirshperday <- subset(vmstoviirshp, DAY == dayn)
    viirsfileperday <- subset(viirsfile, day == dayn)
    vmsfileperday <- subset(vmsfile, day == dayn)
    transmitters <- unique(vmstoviirshperday$TRANSMITTE)
    for (i in transmitters) {
      transnr <- i
      vmstoviirshperdaytransmitte <- subset(vmstoviirshperday, TRANSMITTE == transnr)
      if (nrow(vmstoviirshperdaytransmitte)<1) {
        next
      }
      #bufferdistance <- 500 #knot2kph/2
      buvver <- function(x,y) {
        requireNamespace("sp")
        buv1 <- x[,-(2:8)]
        buv2 <- sp::spTransform(buv1, sp::CRS('+init=epsg:3857'))
        requireNamespace("rgeos")
        buv3 <- rgeos::gBuffer(buv2, byid=TRUE, widt=y)
        buv4 <- sp::spTransform(buv3, sp::CRS('+init=epsg:4326'))
        return(buv4)
      }
      vmstoviirsbuffer <- buvver(vmstoviirshperdaytransmitte,bufferdistance)
      overlayer <- function(y,z) {
        x <- sp::over(y,z,returnlist = TRUE)
        x[x == "list(id = numeric(0))"] <- "NA"
        x[x == "list(id = integer(0))"] <- "NA"
      }
      raster::crs(viirsfileperday) <- raster::crs(vmstoviirsbuffer)
      planka <- viirsfileperday[vmstoviirsbuffer, ]
      if (nrow(planka)<1) {
        next
      }
      planka$TRANSMITTE <- rep(transnr,length(planka))
      planka <- as.data.frame(planka)
      outfile <- rbind2(outfile,planka)
    }

  }
  haron <- dplyr::select(outfile, dplyr::contains("id_Key"), dplyr::contains("TRANSMITTE"))
  texigo <- sp::merge(viirsfile, haron, by = 'id_Key' , duplicateGeoms = TRUE)
  texigo$TRANSMITTE[is.na(texigo$TRANSMITTE)] <- "unidentified"
  return(texigo)
}

max_stars_repo_name: Cglobemaps/LLFI

max_issues_repo_name: Cglobemaps/LLFI

max_forks_repo_name: Cglobemaps/LLFI


text:
## netcontext.R - Using ggnetwork library to place gene models onto their network context
##
## by Artem Sokolov

library( optparse )

## Define accepted command-line arguments
opts <- list(
    make_option(
        c( "-n", "--network" ), type="character", default=NULL,
        help="A .sif file defining the interaction graph. Can be gzipped." ),
    make_option(
        c( "-w", "--weights" ), type="character", default=NULL,
        help="A tab-delimited two-column file that maps genes to their weights. Can be gzipped." ),
    make_option(
        c( "-s", "--seed" ), type="character", default=NULL,
        help="(Optional) A file that contains the seed genes, one per line. If not provided, the seed genes are taken to be the first column of the -w argument." ),
    make_option(
        c( "-o", "--output" ), type="character", default="out.pdf",
        help="Output filename. Default: out.pdf" ) )

## Parse the provided arguments and verify the presence of mandatory ones
op <- OptionParser( option_list = opts )
argv <- parse_args(op)
if( is.null( argv$network ) || is.null( argv$weights ) )
{
    print_help(op)
    stop( "Network and gene weights must be provided." )
}

## Load the remaining libraries
suppressMessages( library( dplyr ) )
suppressMessages( library( igraph ) )
suppressMessages( library( sna ) )
suppressMessages( library( network ) )
library( ggplot2 )
library( ggrepel )
library( ggnetwork )

## "Projects" a set of selected genes (seed) onto a gene-gene interaction network:
##   w - a named vector of gene weights (can be differential expression scores, model weights, etc.)
##   g - an igraph object that encapsulates the interaction network
##   seed - a character vector specifying genes to use as a plot seed
## Returns a ggplot object
network.context <- function( w, g, seed )
{
    ## Generate a connected subgraph over the selected genes
    sp <- shortest_paths( g, seed, seed ) %>% unlist() %>% unique()
    gs <- induced_subgraph( g, sp )

    ## Set up the associated weight vector and populate it with known values
    ww <- setNames( rep(0, length(V(gs))), names(V(gs)) )
    j <- intersect( names(ww), names(w) )
    if( length(j) > 0 )
        ww[j] <- w[j]

    ## Associate the new weight vector with the graph
    V(gs)$Weight <- ww

    ## Distinguish seed and linker nodes
    V(gs)$Type <- c( "Linker", "Seed" )[ names(V(gs)) %in% seed + 1 ]

    ## Plot all the things
    ggplot( gs, aes(x = x, y = y, xend = xend, yend = yend ) ) +
        geom_edges( size=1, aes(linetype=Interaction), color="gray" ) +
        geom_nodes( aes( fill = Weight, size=Type ), color="black", shape=21 ) +
        geom_nodetext_repel( aes( label = vertex.names ), size=4, fontface="bold" ) +
        scale_fill_gradient2( low = "steelblue", high = "tomato", mid="white", midpoint = 0 ) +
        scale_size_manual( values=c( "Linker" = 5, "Seed" = 10 ) ) +
        theme_blank()
}

main <- function( argv )
{
    ## Load the provided network
    cat( "Loading the provided .sif\n" )
    X <- read.delim( argv$network, header=FALSE, as.is=TRUE )

    ## Generate a graph from the list of edges
    cat( "Generating a network\n" )
    g <- X %>% select( V1, V3 ) %>% as.matrix() %>% graph_from_edgelist( directed = FALSE )
    E(g)$Interaction <- X$V2

    ## Handle duplicate edges by defining a new "composite-edge" interaction class
    g <- simplify( g, edge.attr.comb=function(v) {ifelse( length(v) > 1, "composite-edge", v )} )

    ## Isolate the largest connected component of the graph
    if( is_connected( g ) == FALSE )
    {
        cat( "Reducing the graph to the largest connected component\n" )
        lg <- decompose( g )
        g <- lg[[ lapply( lg, vcount ) %>% which.max() ]]
    }

    ## Load the weights
    w <- read.delim( argv$weights, header=FALSE ) %>%
        tibble::column_to_rownames( "V1" ) %>% as.matrix() %>% drop()

    ## Load the seed
    if( is.null( argv$seed ) )
        s <- names(w)
    else
        s <- scan( argv$seed, what=character() )

    ## Reduce the seed set to the available vertex set
    s <- V(g) %>% names() %>% intersect( s )
    if( length(s) < 1 )
        stop( "None of the requested genes are present in the graph." )
    if( length(s) > 50 )
        stop( 'The tool is not meant to plot more than 50 nodes, as this creates \"hairball\" networks that are difficult to interpret.' )

    ## Generate and save the plot
    gg <- network.context( w, g, s )
    ggsave( argv$output, gg )
}

main( argv )

max_stars_repo_name: smkartik/msda

max_issues_repo_name: smkartik/msda

max_forks_repo_name: datarail/msda


text:
Data_Cars <- mtcars

print(dim(Data_Cars))

print(names(Data_Cars))

print(row.names(Data_Cars))

print(summary(Data_Cars))

print(which.max(Data_Cars$hp))
print(which.min(Data_Cars$hp))

print(quantile(Data_Cars$wt))

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:
 ### TEST
        cell.type <- "Monocytes"
        tr <-  "Ustekinumab"
        ts <- "NF-?B"
        stimulation <- "TNF?"

        block<-paste(cell.type, stimulation, ts, sep="__")
        print(block)
        ds.t<-ds[ds$Cell.type==cell.type & ds$Stimulation==stimulation & ds$Transcription.factor==ts,]

        ###TEST
        ds.t<-ds[ds$Cell.type==cell.type & ds$Stimulation==stimulation & ds$Transcription.factor==ts,]
        #t4<-lme(Rd.score ~  1 + Gender +  Ethnicity + Age + Treatment*tp, random = ~ 1|Patient.ID, data = ds.t, na.action="na.omit")
        t4<-lm(Rd.score ~  1 + Gender +  Ethnicity + Age + factor(tp)*Treatment, data = ds.t, na.action="na.omit")

        par(mfcol=c(1,2))
        plot(t4,1:2)
        Anova(t4)

        t4.means<-interactionMeans(t4)
        plot(t4.means)

        interactionMeans(t4, factors="Treatment")

        plot(t4.means, multiple=FALSE)
        testInteractions(t4, fixed="Treatment", across="tp")

        #t1=summary(lme(Rd.score ~  1 + Gender +  Ethnicity + Age + Treatment + tp + Treatment*tp, random = ~ 1|Patient.ID, data = ds.t, na.action="na.omit"))
        #t2=summary(lme(Rd.score ~  1 + Gender +  Ethnicity + Age + Treatment*tp, random = ~ 1|Patient.ID, data = ds.t, na.action="na.omit"))
        #t3=summary(lme(Rd.score ~  1 + Gender +  Ethnicity + Age + Treatment + tp + Treatment:tp, random = ~ 1|Patient.ID, data = ds.t, na.action="na.omit"))

        t4<-lm(relativePASI.log10 ~  1 + Gender +  Ethnicity + Age + Rd.score, data = ds, na.action="na.omit")

        par(mfcol=c(1,2))
        plot(t4,1:2)
        Anova(t4)

        tmp <- expand.grid(tp = factor(unique(ds.t$tp)), Treatment = unique(ds.t$Treatment))

X <- model.matrix(~ Treatment * factor(tp), data = tmp)
glht(mod, linfct = X)
Tukey <- contrMat(table(ds.t$tp), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(ds.t$Treatment)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(ds.t$Treatment)[2], rownames(K2), sep = ":")
K <- rbind(K1, K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))
summary(glht(mod, linfct = K %*% X))
glht(mod, linfct = mcp(x = "Tukey"))

#z3 <- lme(Rd.score ~ GenderId +  Ethnicity + Age + Block2*Treatment, random = ~ 1|Patient.ID, data = ds, na.action="na.omit")
#b=summary(z3)
#interaction.plot(ds$Treatment, ds$Block2, ds$Rd.score, ds$)
#b <- summary(z3)
#b$tTable
#print(contrast(z3, list(Treatment="Adalimumab", Block2="1_Monocytes_NF-?B_IL-17", GenderId ="1", Ethnicity="3", Age=23),       ##AGE, GENDER AND ETHNICITY ARE JUST DUMMY VALUS AND ARE NOT USED IN THE STATISTICS
#                   list(Treatment="Adalimumab",  Block2="1_Monocytes_NF-?B_TNF?", GenderId ="1", Ethnicity="3",Age=23), X=TRUE))
#
#ix1=ds$Treatment=="Adalimumab" & ds$Block2=="1_Monocytes_NF-?B_IL-17"
#ix2=ds$Treatment=="Adalimumab" & ds$Block2=="1_Monocytes_NF-?B_TNF?"
#boxplot(ds$Rd.score[ix1],ds$Rd.score[ix2], notch=TRUE)
#
#fig1 = xyplot(ds$Rd.score ~ ds$Treatment | ds$Block1, pch = 4)
#print(fig1)


text:
# +--------------------------------------------------------+
# | Test of secondary metabolites. Not an official FOCUS   |
# | dummy substance and metabolites. Called here GW-M      |
# |                                                        |
# | The results are compared with those initially obtained |
# | with macrounchained, in order to detect accidental     |
# | changes. No benchmark currently available to compare   |
# | with in the case of secondary metabolites              |
# +--------------------------------------------------------+

rm(list=ls(all=TRUE))

prj_name  <- "macrounchained"
pkg_name  <- "macrounchained"
test_name <- "test_GW_05_secondary-metabolite"

if( as.logical( as.integer( Sys.getenv( "macrounchained_developer", unset = "0" ) ) ) ){
    rPackagesDir <- Sys.getenv( "rPackagesDir", unset = NA_character_ )

    if( !is.na( rPackagesDir ) ){
         output_dir <- file.path( rPackagesDir, prj_name,
            pkg_name, "inst", "more_tests", "results" )
    }else{
        output_dir <- getwd()
    }
}else{
    output_dir <- getwd()
}

library( "macrounchained" )

param <- data.frame(
    "soil"              = "chateaudun",
    "crop"              = "cereals, winter",
    "id"                = c(     1L,       2L,       3L ),
    "name"              = c( "GW-M", "Met-M1", "Met-M2" ),
    "kfoc"              = c(     12,     114,        46 ),
    "nf"                = c(      1,     0.9,       0.9 ),
    "dt50"              = c(     23,      19,       136 ),
    "dt50_ref_temp"     = c(     20,      20,        20 ),
    "dt50_pf"           = c(      2,       2,         2 ),
    "exp_temp_resp"     = c( 0.0948,  0.0948,    0.0948 ),
    "exp_moist_resp"    = c(   0.49,    0.49,      0.49 ),
    "crop_upt_f"        = c(    0.5,     0.5,       0.5 ),
    "diff_coef"         = c(  5E-10,   5E-10,     5E-10 ),
    "parent_id"         = c(     NA,      1L,        2L ),
    "g_per_mol"         = c(    381,     183,       140 ),
    "transf_f"          = c(     NA,     0.2,       1.0 ), # mol met formed / mol parent degraded
    "g_as_per_ha"       = c(     6,        0,         0 ),
    "app_j_day"         = c(   309L,    309L,      309L ),
    stringsAsFactors    = FALSE
)

res <- macrounchainedFocusGW(
    s         = param,
    overwrite = TRUE,
    run       = TRUE )

#   Results obtained on 2019/01/23 with MACRO In FOCUS 5.5.4
#   macrounchained 0.9.0 (git revision: cbfe075),
#   rmacrolite 0.9.2 (git revision: b792510) and
#   macroutils2 2.2.1 (git revision: fafcc2b)
expected_results_s <- data.frame(
    "name"           = c( "GW-M", "Met-M1", "Met-M2" ),
    "ug_per_L_rnd"   = c(  0.140,  0.00269,  0.04150 ),
    "index_period1"  = c(     14,        7,        8 ),
    "index_period2"  = c(      7,       15,        2 ),
    stringsAsFactors        = FALSE )

expected_results_w <- data.frame(
    "name"                  = c(    "GW-M",  "Met-M1",  "Met-M2" ),
    "perc_period1_mm"       = c(  234.0759,  256.3556,  237.0148 ),
    "perc_period2_mm"       = c( 256.35547,  38.68237,  19.49554 ),
    stringsAsFactors        = FALSE )

modelVar  <- rmacrolite::rmacroliteGetModelVar()
log_file  <- file.path( output_dir,
    sprintf( "%s_results.txt", test_name ) )
path      <- rmacrolite::rmacroliteGetModelVar()[["path"]]
log_width <- rmacrolite::getRmlPar( "log_width" )
verbose   <- TRUE

file.copy(
    from = file.path( path, res[["extra_files"]]["log_file"] ),
    to   = log_file,
    overwrite = TRUE )

macrounchained:::.muc_logMessage( m = "Benchmark (results)",
    frame = "*", verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

obtained_results <- res[[ "analyse_summary_output" ]]

macrounchained:::.muc_logMessage( m = "Expected results (water):",
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

macrounchained:::.muc_print( x = expected_results_w,
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

macrounchained:::.muc_logMessage( m = "Expected results (solute):",
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

macrounchained:::.muc_print( x = expected_results_s,
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

macrounchained:::.muc_logMessage( m = "Obtained results (water and solute):",
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

macrounchained:::.muc_print( x = obtained_results,
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

obtained_results_s <- obtained_results[, colnames( obtained_results )
    %in% colnames( expected_results_s ) ]

obtained_results_w <- obtained_results[, colnames( obtained_results )
    %in% colnames( expected_results_w ) ]

diff_w_percent <- ( (expected_results_w[, c( "perc_period1_mm", "perc_period2_mm" ) ] -
                obtained_results[, c( "perc_period1_mm", "perc_period2_mm" ) ]) /
            expected_results_w[, c( "perc_period1_mm", "perc_period2_mm" ) ] ) * 100

macrounchained:::.muc_logMessage( m = "Differences in percolation: %s to %s %% of the expected percolation",
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE,
    values = list( round( min( diff_w_percent ), 2L ), round( max( diff_w_percent ), 2L ) ) )

test_s <- all.equal( obtained_results_s, expected_results_s,
    check.attributes = FALSE )

if( all( is.logical( test_s ) ) ){
    macrounchained:::.muc_logMessage( m = "Test conclusions (solute): ALL EQUAL",
        verbose = verbose, log_width = log_width,
        logfiles = log_file, append = TRUE )

}else{
    macrounchained:::.muc_logMessage( m = "Test conclusions (solute): DIFFERENCES DETECTED",
        verbose = verbose, log_width = log_width,
        logfiles = log_file, append = TRUE )
}

macrounchained:::.muc_logMessage( m = "END",
    verbose = verbose, log_width = log_width,
    logfiles = log_file, append = TRUE )

max_stars_repo_name: julienmoeys/macrounchained

max_issues_repo_name: julienmoeys/macrounchained

max_forks_repo_name: julienmoeys/macrounchained