text:
myTestRule {
#Input parameters are:
#  Session variable
#    Session variable without the $ sign
#    all      - output all of the defined variables
#  Output mode flag:
#    server   - log the output to the server log
#    client   - send the output to the client in rError
#    all      - send the output to both client and server
#Output from running the example is:
#  Variables are written to the log file
#Output in irods/server/log/rodsLog.2011.6.1 log file is:
#  msiGetSessionVarValue: userNameClient=rods
   msiGetSessionVarValue(*A,"server");
   writeLine("stdout","Variables are written to the log file");
}
INPUT *A="userNameClient"
OUTPUT ruleExecOut

ext: r

max_stars_repo_name: PlantandFoodResearch/irods

max_issues_repo_name: PlantandFoodResearch/irods

max_forks_repo_name: PlantandFoodResearch/irods


text:
### Add to history() and avoid repeatedly appending suffix.
addhistory <- function(read)
{
  tmp <- as.character(read)
  suffix <- paste0(" # ", getval(prompt))
  if (!grepl(x = tmp, pattern = paste0(suffix, "$"), perl = TRUE))
    tmp <- paste0(tmp, suffix)

  if (isWindows() && isRStudio())
  {
    ### This adds history to native R but is insufficient for Rtudio IDE.
    # .rs.registerReplaceHook("timestamp", "utils", function(original, ...)
    # {
    #   invisible(.Call("rs_timestamp", "bb"))
    # })
    # timestamp()
  }
  else
    utils::timestamp(stamp = tmp, prefix = "", suffix = "", quiet = TRUE)

  invisible()
}

max_stars_repo_name: cran/remoter

max_issues_repo_name: cran/remoter

max_forks_repo_name: cran/remoter


text:
## -------->>  [[file:../nstandr.src.org::*inset_formals][inset_formals:1]]
inset_formals <- function(formals_args, new_args) {
    checkmate::assert_subset(names(new_args)
                           , choices = names(formals_args)
                           , empty.ok = FALSE)
    for (i in 1:length(new_args)) {
        formals_args[names(new_args)[i]] <- new_args[i]
    }
    return(formals_args)
}

make_alias <- function(fun, ...) {
    .fun <- fun
    ## remove 1st and 2nd elemnts ('call' name and 'fun' arg)
    args <- as.list(sys.call())[-(1:2)]
    checkmate::assert_subset(names(args)
                           , choices = names(formals(.fun))
                           , empty.ok = FALSE)
    formals(.fun)[names(args)] <- args
    return(.fun)
}
## --------<<  inset_formals:1 ends here

max_stars_repo_name: stasvlasov/nstandr

max_issues_repo_name: stasvlasov/nstandr

max_forks_repo_name: stasvlasov/nstandr


text:
rm(list = ls())

gene_posSite<- function(input_file, output_file){
  bed_file = read.table(input_file, sep = '\t')
  posSite = bed_file[,-c(3,5)]
  posSite[,2] = posSite[,2] + 10
  colnames(posSite) = c('chr', 'coordinate', 'strand', 'IPD_ratio', 'methy_level', 'low_95', 'up_95')
  write.csv(file = output_file, posSite, row.names = F, quote = F)
}

# for D. melanogaster
gene_posSite(input_file = '825_m6A_25x_novar_dev03.bed', output_file = 'm6A_dm_pos.csv')

# for A. thaliana
gene_posSite(input_file = 'm6A_TAIR10.bed', output_file = 'm6A_ara_pos.csv')

# for A. thaliana
gene_posSite(input_file = 'm6A_ecoli.bed', output_file = 'm6A_ecoli_pos.csv')


text:
# Assignment: ASSIGNMENT 7
# Name: Angel, Daniel
# Date: 2021-08-01

## Set the working directory to the root of your DSC 520 directory
setwd("C:/Users/Danny/Documents")

## Load the `data/r4ds/heights.csv` to
heights_df <- read.csv("data/r4ds/heights.csv")

# Fit a linear model
earn_lm <-  lm(earn ~ ed + race + height + age + sex, data = heights_df)

# View the summary of your model
summary(earn_lm)

predicted_df <- data.frame(
  earn = predict(earn_lm, heights_df),
  ed=heights_df$ed, race=heights_df$race, height=heights_df$height,
  age=heights_df$age, sex=heights_df$sex
)

## Compute deviation (i.e. residuals)
mean_earn <- mean(heights_df$earn)
mean_earn
## Corrected Sum of Squares Total
sst <- sum((mean_earn - heights_df$earn)^2)
sst
## Corrected Sum of Squares for Model
ssm <- sum((mean_earn - predicted_df$earn)^2)
ssm
## Residuals
residuals <- heights_df$earn - predicted_df$earn
## Sum of Squares for Error
sse <- sum(residuals^2)
sse
## R Squared
r_squared <- ssm / sst
r_squared
## Number of observations
n <- nrow(heights_df)
n
## Number of regression parameters
p <- 8
p
## Corrected Degrees of Freedom for Model
dfm <- (p - 1)
dfm
## Degrees of Freedom for Error
dfe <- (n - p)
dfe
## Corrected Degrees of Freedom Total:   DFT = n - 1
dft <- (n - 1)
dft
## Mean of Squares for Model:   MSM = SSM / DFM
msm <- ssm / dfm
msm
## Mean of Squares for Error:   MSE = SSE / DFE
mse <- sse / dfe
mse
## Mean of Squares Total:   MST = SST / DFT
mst <- sst / dft
mst
## F Statistic
f_score <- msm / mse
f_score
## Adjusted R Squared R2 = 1 - (1 - R2)(n - 1) / (n - p)
adjusted_r_squared <- 1 - ((1-r_squared)*(n-1)/(n - p))
adjusted_r_squared


text:
## Example

## install.packages(c("moments"))
setwd("~/GitHub/SpatialStress")
source("R/functions.r")

# creating an initial landscape object

initial <- init_landscape(states = c("+","0","-"), cover = c(0.4,0.1,0.5))

plot(initial)
plot(initial, cols = c("darkgreen", "grey70", "white"))

summary(initial)

# change size of lattice

initial <- init_landscape(states = c("+","0","-"), cover = c(0.4,0.1,0.5), width = 25, height = 25)

plot(initial, cols = c("darkgreen", "grey70", "white"))

summary(initial)

# run simulation

parms_grazing <- list(
  del = 0.9,
  b = 0.3,
  c_ = 0.2,
  m0 = 0.05,
  g = 0.1,
  r = 0.01,
  f = 0.9,
  d = 0.2,
  protect = 0.9
)

grazingrun <- ca(initial, parms_grazing, model = grazing)

plot(grazingrun)

## specify what to plot
plot(grazingrun, plotstates = c(TRUE, TRUE, TRUE), cols = c("darkgreen", "grey50", "grey80"), lwd = 2)

## get plots of the snapshots
par(mfrow = c(2,5), mar = c(0,0,0,0)+0.4)
for(i in 1:9) {
  plot(grazingrun$timeseries[[i]])
}

dev.off()
# analyse output

## a simple summary
summary(grazingrun)

## a function that returns  indicators / early warning signs
indicators(grazingrun)

# select model

parms_mussel  <- list(
  r = 0.4, # recolonisation of empty sites dependent on local density
  d = 0.9, # probability of disturbance of occupied sites if at least one disturbed site
  delta = 0.01 # intrinsic disturbance rate
)

musselrun <- ca(initial, parms_mussel, model = musselbed)

## specify what to plot
plot(musselrun, plotstates = c(TRUE, FALSE, TRUE), cols = c("black", "grey50", "grey90"), lwd = 2)

## get plots of the snapshots
par(mfrow = c(2,5), mar = c(0,0,0,0)+0.4)
for(i in 1:9) {
  plot(musselrun$timeseries[[i]])
  mtext(musselrun$timeseries[[i]], outer = TRUE)
}
dev.off()

# simulate gradients

parms_mussel_gradient  <- list(
  r = 0.4, # recolonisation of empty sites dependent on local density
  d = seq(0,1,0.1), # probability of disturbance of occupied sites if at least one disturbed site
  delta = 0.01 # intrinsic disturbance rate
)

# provides parallel backend
library(foreach)
library(doSNOW)

workerlist <- c(rep("localhost", times = 7))

cl <- makeSOCKcluster(workerlist, outfile='out_messages.txt')

registerDoSNOW(cl)

musselgradient <- ca_array(initial, parms_mussel_gradient, model = musselbed)

stopCluster(cl)

musselgradient

plot(`+`~ d, data = musselgradient )

max_stars_repo_name: skefi/SpatialStress

max_issues_repo_name: skefi/SpatialStress

max_forks_repo_name: skefi/SpatialStress


text:
; ******** Source: suite-a.asm
     1                          !cpu 65816
     2
     3                          ; 2017-12-13 J.E. Klasek j+816 AT klasek DOT at
     4
     5
     6
     7                          videoram = $0400
     8                          colorram = $d800
     9
    10                          ;-------------------------------------------------------------------------------
    11                          	*=$07ff
    12  07ff 0108               	!word $0801
    13  0801 0b08               	!word bend
    14  0803 0a00               	!word 10
    15  0805 9e                 	!byte $9e
    16  0806 3230363100         	!text "2061", 0
    17  080b 0000               bend:       !word 0
    18                          ;-------------------------------------------------------------------------------
    19
    20  080d 78                 	sei
    21  080e a917               	lda #$17
    22  0810 8d18d0             	sta $d018
    23  0813 a935               	lda #$35
    24  0815 8501               	sta $01
    25  0817 a97f               	lda #$7f
    26  0819 8d0ddc             	sta $dc0d
    27  081c 8d0ddd             	sta $dd0d
    28  081f ad0ddc             	lda $dc0d
    29  0822 ad0ddd             	lda $dd0d
    30  0825 a200               	ldx #0
    31                          -
    32  0827 a920               	lda #$20
    33  0829 9d0004             	sta videoram,x
    34  082c 9d0005             	sta videoram+$0100,x
    35  082f 9d0006             	sta videoram+$0200,x
    36  0832 9d0007             	sta videoram+$0300,x
    37  0835 a901               	lda #1
    38  0837 9d00d8             	sta colorram,x
    39  083a 9d00d9             	sta colorram+$0100,x
    40  083d 9d00da             	sta colorram+$0200,x
    41  0840 9d00db             	sta colorram+$0300,x
    42  0843 e8                 	inx
    43  0844 d0e1               	bne -
    44
    45  0846 4c0010             	jmp start
    46                          theend:
    47  0849 e230               	sep #$30                        ; 8-bit for X/Y and A/M
    48                          	!as
    49                          	!rs
    50  084b af100204           	lda $040210
    51  084f c9ff               	cmp #$ff
    52  0851 d00d               	bne error
    53
    54  0853 a905                      	lda #5
    55  0855 8d20d0             	sta $d020
    56  0858 a200                      	ldx #0 ; success
    57  085a 8effd7             	stx $d7ff
    58  085d 4c5d08             	jmp *
    59                          error
    60  0860 8d0004             	sta $0400
    61  0863 af110204           	lda $040211			; failure map (which test failed)
    62  0867 8d0104             	sta $0401
    63  086a a90a               	lda #10
    64  086c 8d20d0             	sta $d020
    65  086f a2ff               	ldx #$ff ; failure
    66  0871 8effd7             	stx $d7ff
    67  0874 4c7408             	jmp *
    68
    69                          ;-------------------------------------------------------------------------------
    70
    71                          	* = $1000
    72                          start:
    73                          ; EXPECTED FINAL RESULTS: $0210 = FF
    74                          ; (any other number will be the
    75                          ;  test that failed)
    76
    77                          ; initialize:
    78  1000 a900               	lda #$00
    79  1002 8f100204           	sta $040210
    80  1006 8f110204           	sta $040211
    81
    82
    83                          test00:
    84
    85                          ; setup cpu
    86  100a 18                 	clc
    87  100b fb                 	xce                             ; native mode
    88  100c c230               	rep #$30                        ; 16-bit for X/Y and A/M
    89                          	!al
    90                          	!rl
    91
    92                          ; setup registers
    93  100e a90404             	lda #$0404                      ; Data Bank Register register
    94  1011 48                 	pha                             ; akku in 16 bit
    95  1012 ab                 	plb                             ; pull DBR twice
    96  1013 ab                 	plb
    97  1014 a08888             	ldy #$8888                      ; change marker
    98  1017 bb                 	tyx
    99  1018 98                 	tya
   100
   101                          ; setup memory
   102  1019 a95555             	lda #$5555                      ; wrap marker
   103  101c 8f878804           	sta $048887                     ; into bank 4, for LDX/LDY
   104  1020 a97777             	lda #$7777                      ; no-wrap marker
   105  1023 8f878805           	sta $058887			; into bank 5, for LDX/LDY
   106
   107                          ;---------------------------------------------------------------------
   108
   109  1027 9c0000             	stz $0000			; init wrap marker
   110  102a a97777             	lda #$7777                      ; no-wrap marker
   111  102d 8f000005           	sta $050000                     ; to start of bank 5
   112
   113  1031 8cffff             	sty $ffff                       ; high byte of Y is where?
   114  1034 ad0000             	lda $0000
   115  1037 d011               	bne +
   116  1039 adffff             	lda $ffff			; fetch, does not wrap
   117  103c c98888             	cmp #$8888
   118  103f d009               	bne +
   119  1041 af000005           	lda $050000
   120  1045 c98877             	cmp #$7788			; write to bank 5
   121  1048 f004               	beq ++
   122  104a ee1002             +	inc $0210			; fail counter
   123  104d 18                 	clc
   124                          ++
   125  104e 2e1102             	rol $0211			; update failure map
   126                          ;---------------------------------------------------------------------
   127
   128  1051 9c0000             	stz $0000			; init wrap marker
   129  1054 a97777             	lda #$7777                      ; no-wrap marker
   130  1057 8f000005           	sta $050000                     ; to start of bank 5
   131
   132  105b bb                 	tyx                             ; change marker
   133  105c 8effff             	stx $ffff                       ; high byte of Y is where?
   134  105f ad0000             	lda $0000
   135  1062 d011               	bne +
   136  1064 adffff             	lda $ffff			; fetch, does not wrap
   137  1067 c98888             	cmp #$8888
   138  106a d009               	bne +
   139  106c af000005           	lda $050000
   140  1070 c98877             	cmp #$7788			; write to bank 5
   141  1073 f004               	beq ++
   142  1075 ee1002             +	inc $0210			; fail counter
   143  1078 18                 	clc
   144                          ++
   145  1079 2e1102             	rol $0211			; update failure map
   146                          ;---------------------------------------------------------------------
   147
   148  107c bcffff             	ldy $ffff,x     ; Y=5555  Y=7777  value for Y comes from which bank?
   149  107f c07777             	cpy #$7777
   150  1082 f004               	beq +
   151  1084 ee1002             	inc $0210			; fail counter
   152  1087 18                 	clc
   153                          +
   154  1088 2e1102             	rol $0211			; update failure map
   155                          ;---------------------------------------------------------------------
   156
   157  108b 9b                 	txy				; reinitialize y
   158  108c beffff             	ldx $ffff,y     ; X=5555  X=7777  value for X comes from which bank?
   159  108f e07777             	cpx #$7777
   160  1092 f004               	beq +
   161  1094 ee1002             	inc $0210			; fail counter
   162  1097 18                 	clc
   163                          +
   164  1098 2e1102             	rol $0211			; update failure map
   165                          ;---------------------------------------------------------------------
   166
   167  109b 9c0000             	stz $0000			; init wrap marker
   168  109e a97777             	lda #$7777                      ; no-wrap marker
   169  10a1 8f000005           	sta $050000                     ; to start of bank 5
   170
   171  10a5 a98877             	lda #$7788
   172  10a8 ee0000             	inc $0000			; $0000 = 1
   173  10ab 1cffff             	trb $ffff                       ; 88 77 & ^(88 77) -> 00 00
   174  10ae ad0000             	lda $0000
   175  10b1 c90100             	cmp #$0001			; $0000 not reset by trb (does not wrap)
   176  10b4 d009               	bne +
   177  10b6 af000005           	lda $050000
   178  10ba c90077             	cmp #$7700			; $050001 reset by trb
   179  10bd f004               	beq ++
   180  10bf ee1002             +	inc $0210			; fail counter
   181  10c2 18                 	clc
   182                          ++
   183  10c3 2e1102             	rol $0211			; update failure map
   184                          ;---------------------------------------------------------------------
   185
   186  10c6 a98877             	lda #$7788
   187  10c9 8f000005           	sta $050000			; 00 88 | 88 77 -> 88 ff
   188  10cd 0cffff             	tsb $ffff                       ; set bits (which are already cleared)
   189  10d0 ad0000             	lda $0000
   190  10d3 c90100             	cmp #$0001			; $0000 not set by tsb (does not wrap!)
   191  10d6 d009               	bne +
   192  10d8 af000005           	lda $050000
   193  10dc c9ff77             	cmp #$77ff			; $050001 all bits set by tsb
   194  10df f004               	beq ++
   195  10e1 ee1002             +	inc $0210			; fail counter
   196  10e4 18                 	clc
   197                          ++
   198  10e5 2e1102             	rol $0211			; update failure map
   199                          ;---------------------------------------------------------------------
   200
   201
   202                          test00pass:
   203  10e8 ad1002             	lda $0210
   204  10eb 49003f             	eor #%0011111100000000		; invert failure map
   205  10ee 8d1002             	sta $0210
   206  10f1 d003               	bne +
   207  10f3 ce1002             	dec $0210			; 0 -> FF
   208                          +
   209
   210  10f6 a90000             	lda #$0000
   211  10f9 48                 	pha
   212  10fa ab                 	plb
   213  10fb ab                 	plb				; program bank = 0
   214  10fc 38                 	sec
   215  10fd fb                 	xce				; emulation mode
   216  10fe e230               	sep #$30			; a/m, x/y 8 bit
   217
   218  1100 4c4908                     jmp theend


text:
	maximize_wrench <- function(H_matrix_input, bounds_tuple_of_numeric, direction) {
		positive_fx_direction_constraint <- a_matrix_lhs_direction(H_matrix_input, direction = direction, bounds_tuple_of_numeric)
	    indices_for_muscles <- 1:ncol(H_matrix_input)
	    fmax_info <- lpsolve_force_in_dir("max", positive_fx_direction_constraint, indices_for_muscles, nrow(H_matrix_input))
	}

max_stars_repo_name: bc/stfeasibility

max_issues_repo_name: bc/stfeasibility

max_forks_repo_name: bc/stfeasibility


text:
#' Convert spatial data files to GeoJSON from various formats.
#'
#' You can use a web interface called Ogre, or do conversions locally using the
#' rgdal package.
#'
#' @export
#' @param input The file being uploaded, path to the file on your machine.
#' @param output Destination for output geojson file. Defaults to current
#' working directory, and gives a random alphanumeric file name
#' @param encoding (character) The encoding passed to
#' \code{\link[rgdal]{readOGR}}.  Default: CP1250
#' @param verbose (logical) Printing of \code{\link[rgdal]{readOGR}} progress.
#' Default: \code{FALSE}
#' @template read
#' @examples \dontrun{
#' file <- system.file("examples", "norway_maple.kml", package = "geojsonio")
#'
#' # KML type file - using the web method
#' file_to_geojson(input=file, method='web', output='kml_web')
#' ## read into memory
#' file_to_geojson(input=file, method='web', output = ":memory:")
#' file_to_geojson(input=file, method='local', output = ":memory:")
#'
#' # KML type file - using the local method
#' file_to_geojson(input=file, method='local', output='kml_local')
#'
#' # Shp type file - using the web method - input is a zipped shp bundle
#' file <- system.file("examples", "bison.zip", package = "geojsonio")
#' file_to_geojson(file, method='web', output='shp_web')
#'
#' # Shp type file - using the local method - input is the actual .shp file
#' file <- system.file("examples", "bison.zip", package = "geojsonio")
#' dir <- tempdir()
#' unzip(file, exdir = dir)
#' list.files(dir)
#' shpfile <- file.path(dir, "bison-Bison_bison-20130704-120856.shp")
#' file_to_geojson(shpfile, method='local', output='shp_local')
#'
#' # Neighborhoods in the US
#' ## beware, this is a long running example
#' # url <- 'http://www.nws.noaa.gov/geodata/catalog/national/data/ci08au12.zip'
#' # out <- file_to_geojson(input=url, method='web', output='cities')
#'
#' # geojson with .json extension
#' ## this doesn't work anymore, hmmm
#' # x <- gsub("\n", "", paste0('https://gist.githubusercontent.com/hunterowens/
#' # 25ea24e198c80c9fbcc7/raw/7fd3efda9009f902b5a991a506cea52db19ba143/
#' # wards2014.json', collapse = ""))
#' # res <- file_to_geojson(x)
#' # jsonlite::fromJSON(res)
#' # res <- file_to_geojson(x, method = "local")
#' # jsonlite::fromJSON(res)
#' }

file_to_geojson <- function(input, method = "web", output = ".", parse = FALSE,
                            encoding = "CP1250", verbose = FALSE, ...) {

  method <- match.arg(method, choices = c("web", "local"))
  if (!inherits(parse, "logical")) stop("parse must be logical", call. = FALSE)

  input <- handle_remote(input)
  mem <- ifelse(output == ":memory:", TRUE, FALSE)

  if (method == "web") {
    url <- "http://ogre.adc4gis.com/convert"
    tt <- httr::POST(url, body = list(upload = httr::upload_file(input)))
    if (tt$status_code > 201) {
      res <- tryCatch(
        httr::content(tt, as = "text", encoding = "UTF-8"),
        error = function(e) e
      )
      if (inherits(res, "error")) httr::stop_for_status(tt)
      res2 <- tryCatch(
        jsonlite::fromJSON(res),
        error = function(e) e
      )
      if (inherits(res2, "error")) httr::stop_for_status(tt)
      stop(paste0(res2[[1]], collapse = "\n"), call. = FALSE)
    }
    out <- httr::content(tt, as = "text", encoding = "UTF-8")
    if (mem) {
      jsonlite::fromJSON(out, parse)
    } else {
      if (output == ".") {
        temp <- tempfile()
        on.exit(unlink(temp))
        output <- basename(temp)
      }
      fileConn <- file(paste0(output, ".geojson"))
      writeLines(out, fileConn)
      close(fileConn)
      file_at(output)
      file_ret(path.expand(output))
    }
  } else {
    fileext <- ftype(input)
    if (output == ":memory:") {
      temp <- tempfile()
      output <- temp
      on.exit(unlink(temp))
      on.exit(unlink(output), add = TRUE)
    } else if (output == ".") {
      temp <- tempfile()
      output <- basename(temp)
      on.exit(unlink(temp))
      on.exit(unlink(output), add = TRUE)
    }

    output <- path.expand(output)
    if (fileext == "kml") {
      my_layer <- rgdal::ogrListLayers(input)
      x <- rgdal::readOGR(input, layer = my_layer[1],
                          drop_unsupported_fields = TRUE,
                          verbose = FALSE, stringsAsFactors = FALSE,
                          encoding = encoding, ...)
      write_ogr2(x, output)
      if (mem) {
        from_json(output, parse)
      } else {
        file_at(output)
        file_ret(output)
      }
    } else if (fileext == "shp") {
      x <- rgdal::readOGR(input, rgdal::ogrListLayers(input),
                          verbose = FALSE, stringsAsFactors = FALSE,
                          encoding = encoding, ...)
      write_ogr2(x, output)
      if (mem) {
        from_json(output, parse)
      } else {
        file_at(output)
        file_ret(output)
      }
    } else if (fileext %in% c("geojson", "json")) {
      unlink(paste0(output, ".geojson"))
      x <- rgdal::readOGR(input, rgdal::ogrListLayers(input),
                          verbose = FALSE, stringsAsFactors = FALSE,
                          encoding = encoding, ...)

      write_ogr2(x, output)
      if (mem) {
        from_json(output, parse)
      } else {
        file_at(output)
        file_ret(output)
      }
    } else {
      stop("only .shp, .kml, .topojson, and url's are supported",
           call. = FALSE)
    }
  }
}

write_ogr2 <- function(x, y) {
  unlink(paste0(y, ".geojson"))
  rgdal::writeOGR(x, paste0(y, ".geojson"), basename(y), driver = "GeoJSON",
                  check_exists = FALSE)
}

file_at <- function(x) {
  message(paste0("Success! File is at ", x, ".geojson"))
}

file_ret <- function(x) {
  invisible(paste0(x, ".geojson"))
}

from_json <- function(x, parse) {
  on.exit(unlink(paste0(x, ".geojson")))
  jsonlite::fromJSON(paste0(x, ".geojson"), parse)
}

ftype <- function(z) {
  if (is.url(z)) {
    z <- httr::parse_url(z)$path
  }
  fileext <- strsplit(z, "\\.")[[1]]
  fileext[length(fileext)]
}

# If given a url for a zip file, download it give back a path to the temporary file
handle_remote <- function(x){
  if (!is.url(x)) {
    return(x)
  } else {
    tfile <- tempfile(fileext = paste0(".", ftype(x)))
    res <- httr::GET(x, httr::write_disk(tfile))
    httr::stop_for_status(res)
    res$request$output$path
  }
}

max_stars_repo_name: pvictor/geojsonio

max_issues_repo_name: pvictor/geojsonio

max_forks_repo_name: pvictor/geojsonio


text:
#' @export
plotArrows = function (Y,X,add=F,colfun=rainbow,startend=T,...) {

  #//trace arrows around x,y points from start to end
  xl = c(min(X),max(X))
  yl = c(min(Y),max(Y))
  if(Y[length(Y)]==0) {Y=Y[-length(Y)]}
  if(length(Y) !=length(X)) {X=X[-length(X)]}
  n = length(Y)
  cl = colfun(n)
  if(add==F) plot(1,1,ylim=yl,xlim=xl,type='n')

  for(i in 1:n) {
    arrows(x0=X[i],x1=X[i+1],y0=Y[i],y1=Y[i+1],angle=45,col=cl[i],...)
  }
  points(X,Y,pch=16,cex=1,col='black')
  if(startend){
    points(X[1],Y[1],pch=17,cex=2,col='black')
    points(X[n],Y[n],pch=17,cex=2,col='black')
  }
}


text:
ver_crassy     = Sys.getenv("CRASSY_VER")
ver_cqlconnect = Sys.getenv("CQLCONNECT_VER")
ver_subprocess = Sys.getenv("SUBPROCESS_VER")

install_github = function(pkg, ver, repo, sep) {
    cat("Installing", pkg, ver, "\n")
    remotes::install_github(file.path(repo, ver, fsep = sep))
}

install_github("crassy", ver_crassy, "AkhilNairAmey/crassy", "@v")
install_github("CQLConnect", ver_cqlconnect, "AkhilNairAmey/CQLConnect", "@v")
install_github("subprocess", ver_subprocess, "lbartnik/subprocess", "@")

max_stars_repo_name: CipherAPJNP/cipher-base

max_issues_repo_name: CipherAPJNP/cipher-base

max_forks_repo_name: CipherAPJNP/cipher-base


text:
context("readFile")

test_that("file separator", {
  expect_equal(readFile("../../tests/shinytest/test-data/passes_everything.csv")$fileSeparator, ",")
})

max_stars_repo_name: cjrace/data-screening-shinyapp

max_issues_repo_name: cjrace/data-screening-shinyapp

max_forks_repo_name: cjrace/data-screening-shinyapp


text:
source("_common.r")
x = readLines("../runs/summit/dgemm/openblas.txt")
df = process(x, gpu=FALSE)

plotter(df, color=nthreads) +
  labs(color="Number of Threads") +
  ggtitle("Square Matrix Product", subtitle="OpenBLAS from R")

ggsave(last_plot(), file="openblas.pdf")

max_stars_repo_name: wrathematics/matprodbench

max_issues_repo_name: wrathematics/matprodbench

max_forks_repo_name: wrathematics/matprodbench


text:
<!DOCTYPE html>

  <meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">

<title>Make the UpSet plot — UpSet • ComplexHeatmap</title>

<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.3.7/css/bootstrap.min.css" integrity="sha256-916EbMg70RQy9LHiGkXzG8hSg9EdNy97GazNG/aiY1w=" crossorigin="anonymous" />

<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.7.1/css/all.min.css" integrity="sha256-nAmazAk6vS34Xqo0BSrTb+abbtFlgsFK7NKSi6o7Y78=" crossorigin="anonymous" />
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.7.1/css/v4-shims.min.css" integrity="sha256-6qHlizsOWFskGlwVOKuns+D1nB6ssZrHQrNj1wGplHc=" crossorigin="anonymous" />

<link href="../pkgdown.css" rel="stylesheet">

<meta property="og:title" content="Make the UpSet plot — UpSet" />

<meta property="og:description" content="Make the UpSet plot" />
<meta name="twitter:card" content="summary" />

      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        Toggle navigation

      </button>

        ComplexHeatmap
        2.1.0

  Reference

    Articles

      UNKNOWN TITLE

      UNKNOWN TITLE

Make the UpSet plot

UpSet.rd

Make the UpSet plot

UpSet(m,
    comb_col = "black",
    pt_size = unit(3, "mm"), lwd = 2,
    bg_col = "#F0F0F0", bg_pt_col = "#CCCCCC",
    set_order = order(set_size(m), decreasing = TRUE),
    comb_order = if(attr(m, "set_on_rows")) {
    order.comb_mat(m[set_order, ], decreasing = TRUE)
    } else {
    order.comb_mat(m[, set_order], decreasing = TRUE)
    },
    top_annotation = upset_top_annotation(m),
    right_annotation = upset_right_annotation(m),
    row_names_side = "left",
    ...)

Arguments

    <colgroup><col class="name" /><col class="desc" /></colgroup>

      m

A combination matrix returned by make_comb_mat. The matrix can be transposed to switch the position of sets and combination sets.

      comb_col

The color for the dots representing combination sets.

      pt_size

The point size for the dots representing combination sets.

      lwd

The line width for the combination sets.

      bg_col

Color for the background rectangles.

      bg_pt_col

Color for the dots representing the set is not selected.

      set_order

The order of sets.

      comb_order

The order of combination sets.

      top_annotation

A HeatmapAnnotation object on top of the combination matrix.

      right_annotation

A HeatmapAnnotation object on the right of the combination matrix.

      row_names_side

The side of row names.

      ...

Other arguments passed to Heatmap.

Details

By default, the sets are on rows and combination sets are on columns. The positions of the
two types of sets can be switched by transposing the matrix.

When sets are on rows, the default top annotation is the barplot showing the size of each
combination sets and the default right annotation is the barplot showing the size of the sets.
The annotations are simply constructed by HeatmapAnnotation and anno_barplot with some
parameters pre-set. Users can check the source code of upset_top_annotation and
upset_right_annotation to find out how the annotations are defined.

To change or to add annotations, users just need to define a new HeatmapAnnotation object.
E.g. if we want to change the side of the axis and name on top annotation:

    Upset(..., top_annotation =
        HeatmapAnnotation(
           "Intersection size" = anno_barplot(
               comb_size(m),
               border = FALSE,
               gp = gpar(fill = "black"),
               height = unit(2, "cm"),
               axis_param = list(side = "right")
           ),
           annotation_name_side = "right",
           annotation_name_rot = 0)
    )

To add more annotations on top, users just add it in HeatmapAnnotation:

    Upset(..., top_annotation =
        HeatmapAnnotation(
           "Intersection size" = anno_barplot(
               comb_size(m),
               border = FALSE,
               gp = gpar(fill = "black"),
               height = unit(2, "cm"),
               axis_param = list(side = "right")
           ),
           "anno1" = anno_points(...),
           "anno2" = some_vector,
           annotation_name_side = "right",
           annotation_name_rot = 0)
    )

And so is for the right annotations.

UpSet returns a Heatmap-class object, which means, you can add it with other heatmaps and annotations
by + or %v%.

Examples

set.seed(123)
lt = list(a = sample(letters, 10),
          b = sample(letters, 15),
          c = sample(letters, 20))
m = make_comb_mat(lt)
UpSet(m)

UpSet(t(m))

m = make_comb_mat(lt, mode = "union")
UpSet(m)

UpSet(m, comb_col = c(rep(2, 3), rep(3, 3), 1))

# compare two UpSet plots
set.seed(123)
lt1 = list(a = sample(letters, 10),
          b = sample(letters, 15),
          c = sample(letters, 20))
m1 = make_comb_mat(lt1)
set.seed(456)
lt2 = list(a = sample(letters, 10),
          b = sample(letters, 15),
          c = sample(letters, 20))
m2 = make_comb_mat(lt2)

max1 = max(c(set_size(m1), set_size(m2)))
max2 = max(c(comb_size(m1), comb_size(m2)))

UpSet(m1, top_annotation = upset_top_annotation(m1, ylim = c(0, max2)),
    right_annotation = upset_right_annotation(m1, ylim = c(0, max1)),
    column_title = "UpSet1") +
UpSet(m2, top_annotation = upset_top_annotation(m2, ylim = c(0, max2)),
    right_annotation = upset_right_annotation(m2, ylim = c(0, max1)),
    column_title = "UpSet2")

Contents

Arguments

Details

Examples

Developed by Zuguang Gu.

Site built with pkgdown 1.4.1.

ext: rd


text:
# Script to simulate control interventions in fitted dengue transmission model

library(rstan)
library(doParallel)
library(plyr)
library(magrittr)
library(lubridate)

# Loading and processing mosquito and case data

tseries <- read.csv("Data/Vitoria.data.csv")

data <- ddply(tseries, .(tot.week), summarise,
              obs = sum(Cases, na.rm = T),
              q = sum(Mosquitoes, na.rm = T),
              tau = sum(Trap, na.rm = T),
              year = unique(Year),
              week = unique(Week)
              )

ends <- daply(data, .(year), function(df) max(df$tot.week))

# Population size
pop <- 327801

# Wrangling weather data and computing EIP
weather <- read.csv("Data/Vitoria.weather.csv") %>%
  mutate(date = ymd(BRST), year = year(date), week = week(date))

weather <- subset(weather, date < date[1] + weeks(243))
weather$tot.week <- rep(c(1:243), each = 7)

covars <- ddply(weather, .(tot.week), summarise,
                temp = mean(Mean.TemperatureC, na.rm = T))

rov <- 7 * exp(0.21 * covars$temp - 7.9)

# Setting imports (if any)
firsts <- daply(data, .(year), function(df) min(df$tot.week))[2:5]
import <- rep(0, 243)
# import[firsts - 8] <- 10 / pop

model <- stan_model(file = "Code/gammaeip.stan", verbose = T, save_dso = TRUE, auto_write = TRUE)

#===============================================================================
# Loading mcmc samples

chain <- readRDS("Results/chain.rds")
samples <- rstan::extract(chain, permute = T)[1:18]

fitcases <- rstan::extract(chain, "state", permute = T)[[1]] %>%
  extract(, ,  11) %>%
  apply(1, diff)

nmcmc <- length(samples[[1]])

extract_sample <- function(x, k){
  if(length(dim(x)) == 1) return(x[k])
  else return(x[k, ])
}

#===============================================================================
# Adult control simulations

# Loops over each HMC iteration and simulates a 5% decresase in total mosquito
# population size in each of the first 210 weeks of the time series.
# Summarizes the difference in number of cases observed in each year between the
# control and uncontrolled time series.

# Setting up parallel
cl <- makeCluster(8, type = "SOCK")
registerDoParallel(cl)

# Parallel for-loop over mcmc iterations
reduction <- foreach(k = 1:nmcmc, .combine = "rbind", .packages = c("rstan", "magrittr", "lubridate")) %:%
  foreach(j = 1:191, .combine = "rbind") %dopar% {

    # Extract parameter values in kth iteration
    inits <- list(lapply(samples, extract_sample, k))

    # Set-up control structure
    control <- matrix(1.0, nrow = 243, ncol = 3)
    control[j, 3] <- 0.95
    control[j + 1, 2] <- 1.05

    # Stan data object
    dat.stan <- list(T = 243,
                     steps = 7,
                     y = data$obs,
                     q = data$q,
                     tau = data$tau,
                     rov = rov,
                     control = control,
                     import = import,
                     pop = pop)

    # Simulate trajectory with control
    sim <- sampling(model,
                    data = dat.stan,
                    init = inits,
                    iter = 1,
                    chains = 1,
                    warmup = 0,
                    algorithm = "Fixed_param")

    # Extract values of the state variables at the beginning of each year
    cases <- rstan::extract(sim, "state", permute = T)[[1]][1, , 11] %>% diff()

    # Assemble results
    data.frame(rep = k,
               control = j,
               dvmu = inits[[1]]$dv0,
               delta = inits[[1]]$delta_c,
               phi = inits[[1]]$phi_y,
               week = j:243,
               relweek = 0:(243 - j),
               year = year(ymd("2008-01-01") + weeks(j:243)),
               cases = cases[j:243],
               cases0 = fitcases[j:243, k])
}

stopCluster(cl)

saveRDS(reduction, "Results/control.rds")

# Effect of control in the number of cases in the following year
moving <- ddply(control, .(rep, control), summarise,
                tcases = sum(cases[relweek < 53]),
                tcases0 = sum(cases0[relweek < 53]),
                ratio = tcases / tcases0,
                diff = tcases0 - tcases,
                dvmu = mean(dvmu),
                delta = mean(delta),
                phi = mean(phi))

saveRDS(moving, "Results/moving_control.rds")

max_stars_repo_name: clint-leach/mosquito-recon

max_issues_repo_name: clint-leach/mosquito-recon

max_forks_repo_name: clint-leach/mosquito-recon


text:
probitFit = function(y, X,  prior_var_vec, iters = 100){
	#' MCMC Samples for Probit Regression
	#'
	#' \code{probitFit} generates MCMCM samples of Probit regression, following
	#' the data-augmentation Gibbs sampler of Albert and Chib (1993).
	#' Rquires truncnorm package to sample from truncated normal.
	#' Assumes prior centered at 0, and diagonal prior covariance structure.
	#' Assumes response y is coded as 0-1.
	#'
	#' @param y <nX1> observed data vector
	#' @param X <nXp> design matrix (takes objects of class matrix or sparce Matrix)
	#' @param prior_var_vec <px1> prior variance vector
	#' @param iters Integer; number of desired output iterations

	#Number of coefficients
	n = dim(X)[1]
	p = dim(X)[2]
	B_inv = diag(1/prior_var_vec, ncol = p) #inverse prior variance mat

	#Truncation limits, 1 row per observation, colums are lower and upper limits.
	trunc_limits = matrix(0, nrow = n, ncol = 2)
	trunc_limits[y==0 ,1] = -Inf
	trunc_limits[y==1, 2] = Inf

	#Precompute Cholesky decomposition (to efficiently sample from multivariate normal)
	XtX = Matrix::crossprod(X)
	Xty = Matrix::crossprod(X,y)
	A = XtX + B_inv
	R = chol(A) #Cholesky decomposition: A = t(R) %*% R

	#Start MCMC algorithm at Least-Squares estimate (guarenteed to exist w/ nonzero prior variance)
	beta_ls = linCholSolver(R, Xty)

	#Generate the normal rancom variates outside of loop
	rnorm_mat = matrix( rnorm(p * iters, mean = 0, sd = 1),
				ncol = p)

	#Helper function to sample Z|... from truncated normal
	sampleZ = function(beta_vec){
		#Samples z given beta from truncated normal distribtion
		truncnorm::rtruncnorm(n, trunc_limits[ ,1], trunc_limits[ ,2],
							mean =  as.numeric(X %*% beta_vec))
	}

	#MCMC storage, each row is an iteration
	beta_mat = matrix(0, nrow = iters, ncol = p )
	beta_mat[1, ] = as.numeric(beta_ls)
	z_latent = sampleZ(beta_mat[1, ])
	for (i in 2:iters){
		#Gibbs sample beta given z
		Xtz = Matrix::crossprod(X, z_latent)
		beta_z = linCholSolver(R,  Xtz)
		beta_mat[i, ] =  as.numeric(beta_z +
			backsolve(R, rnorm_mat[i, ]))

		#Gibbs sample z
		z_latent = sampleZ(beta_mat[i, ])
	}
	return(beta_mat)
}

max_stars_repo_name: dhelkey/dgrank

max_issues_repo_name: dhelkey/dgrank

max_forks_repo_name: dhelkey/dgrank


text:
if (!exists("key")) message("You need to define a key variable containing the contents of your google maps API key")

geocode_url <- function(address) {
  geo <- "http://maps.google.com/maps/geo?"

  params <- c(
    key = key,
    q = address,
    output = "csv"
  )

  p <- paste(names(params), "=", laply(params,URLencode), sep="", collapse ="&")
  paste(geo, p, sep="")
}

geocode <- function(addresses) {
  paths <- laply(addresses, geocode_url)

  read <- function(path) {
    suppressWarnings(read.csv(url(path), header = FALSE))
  }
  df <- ldply(paths, read, .progress = "text")
  names(df) <- c("status", "accuracy", "lat", "long")
  df
}

addy <- function(row) {
  paste(row$street, ", ", row$city, " ", row$zip, sep="")
}

addresses <- read.csv("addresses.csv")

while(TRUE) {
  todo <- which(is.na(addresses$lat))[1:100]
  todo <- todo[!is.na(todo)]
  stopifnot(length(todo) > 0)
  loc <- geocode(addy(addresses[todo, ]))

  addresses[todo, names(loc)] <- loc
  write.table(addresses, "addresses.csv", sep=",", row=F)
}

addresses[!is.na(addresses$accuracy) & addresses$accuracy == 0, c("lat", "long")] <- NA

addresses[!is.na(addresses$lat) & addresses$lat < 30, c("lat", "long")] <- NA
qplot(long, lat, data=addresses)

max_stars_repo_name: hadley/sfhousing

max_issues_repo_name: hadley/sfhousing

max_forks_repo_name: hadley/sfhousing


text:
#READ.ME
#this script should be run AFTER the preselect script
#it samples the preselect data to create a test set that can be used for debugging
#user can set the percentage of the preselected script that they want

#THIS SCRIPT MUST BE SAVED IN THE CDMINSTANCES/PRESELECTION FOLDER

# default 10% of preselect_person_id will be sampled and used to generate a test set

rm(list=ls())
if(!require(rstudioapi)){install.packages("rstudioapi")}
if(!require(data.table)){install.packages("data.table")}
library(rstudioapi)
library(data.table)

#should be CDMInstances/preselect
projectFolder<-dirname(rstudioapi::getSourceEditorContext()$path)
setwd(projectFolder)

#set up test_sample folder
dir.create(paste0(projectFolder,"_sample_test"))

sample_folder<-(paste0(projectFolder,"_sample_test"))

my_perc<-0.1

tables_preselect<-list()
tables_preselect$EVENTS<-list.files(paste0(projectFolder,"/"), pattern="^EVENTS")
tables_preselect$MEDICAL_OBSERVATIONS<-list.files(paste0(projectFolder,"/"), pattern="^MEDICAL_OBSERVATIONS")
tables_preselect$SURVEY_OBSERVATIONS<-list.files(paste0(projectFolder,"/"), pattern="^SURVEY_OBSERVATIONS")
tables_preselect$MEDICINES<-list.files(paste0(projectFolder,"/"), pattern="^MEDICINES")
tables_preselect$VACCINES<-list.files(paste0(projectFolder,"/"), pattern="^VACCINES")
tables_preselect$SURVEY_ID<-list.files(paste0(projectFolder,"/"), pattern="^SURVEY_ID")
tables_preselect$EUROCAT<-list.files(paste0(projectFolder,"/"), pattern="^EUROCAT")
tables_preselect$PERSONS<-list.files(paste0(projectFolder,"/"), pattern="^PERSONS")

#get persons table
PERSONS<- fread(paste0(projectFolder,"/", tables_preselect$PERSONS[1]))
person_id<-PERSONS$person_id
num_sample<-round((length(person_id)*my_perc),0)
sample_id<-sample(person_id, num_sample)

tables_preselect<-unlist(tables_preselect)

for(i in 1:length(tables_preselect)){
  tablename<-(tables_preselect[i])
  mytable<-fread(paste0(projectFolder,"/",tablename))
  sample_table<-mytable[mytable$person_id%in%sample_id,]
  fwrite(sample_table, paste0(sample_folder, "/", tablename), row.names = F)
}


text:
.sp
.ce 1
\f(CHC\fR
.sp

max_stars_repo_name: gspu/Coherent

max_issues_repo_name: gspu/Coherent

max_forks_repo_name: gspu/Coherent


text:
#' Calculate the longitudinal dispersivity value based on the approach developed by Xu and Eckstein (1995)
#'
#' This function takes plume length (PLUME_L) and length units (UNITS) as arguments to calculate
#' longitudinal dispersivity (α_L). Length units can be meters or feet. The function will return the
#' longitudinal dispersivity in the same length units as plume length.
#'
#' disp(PLUME_L, UNITS)
#' @param PLUME_L is the plume length
#' @param UNITS is the length unit for plume length and longitudinal dispersivity
#' @export
#' @examples
#' disp(PLUME_L = 880, UNITS = "feet")
#' [1] 23.18827

disp <- function(PLUME_L = 1, UNITS = "feet"){
    CONV    <- ifelse(toupper(UNITS) %in% c("FEET", "FT"), 1 * 12 * 2.54 / 100,
               ifelse(toupper(UNITS) %in% c("METER", "METERS", "M"), 1, NA))
    ALPHA   <- 0.83 * log10(PLUME_L * CONV)^2.414 / CONV
    return(ALPHA)
}

max_stars_repo_name: dpphat/MFtools

max_issues_repo_name: dpphat/MFtools

max_forks_repo_name: dpphat/MFtools


text:
library(lmtest)
library(corrplot)

ejercicio1 <- read.csv2('../data/ejercicio1.csv', header = TRUE ,sep = ',')
ejercicio2 <- read.csv('../data/Advertising.csv', header = TRUE)

#Ejercicio 1

M <- cor(ejercicio1)
corrplot(M,method = 'number')

png(filename='../Images/Dispersion1.png')
pairs(ejercicio1)
dev.off()

e1 <- lm(Ventas~Gastos, data=ejercicio1)
summary(e1)

e2 <- lm(Ventas~Gastos+Mes, data=ejercicio1)
summary(e2)

mean(e1$residuals)
sum(e1$residuals)

mean(e2$residuals)
sum(e2$residuals)

png(filename='../Images/qqplot1.png')
qqnorm(e1$residuals)
qqline(e1$residuals)
dev.off()

png(filename='../Images/qqplot2.png')
qqnorm(e2$residuals)
qqline(e2$residuals)
dev.off()

png(filename='../Images/hist1.png')
hist(e1$residuals)
dev.off()

png(filename='../Images/hist2.png')
hist(e2$residuals)
dev.off()

dwtest(e1)
dwtest(e2)

png(filename='../Images/standard1.png')
plot(ejercicio1$Ventas, rstandard(e1), ylab = "Residuos estandarizados", xlab = "Ventas")
abline(h =0, lty = 2)
dev.off()

png(filename='../Images/standard2.png')
plot(ejercicio1$Ventas, rstandard(e2), ylab = "Residuos estandarizados", xlab = "Ventas")
abline(h =0, lty = 2)
dev.off()

#Ejercicio 2

M2 <- cor(ejercicio2)
corrplot(M2,method = 'number')

png(filename='../Images/Dispersion2.png')
pairs(ejercicio2)
dev.off()

e3 <- lm(sales~TV+radio+newspaper, data=ejercicio2)
summary(e3)

mean(e3$residuals)
sum(e3$residuals)

dwtest(e3)

png(filename='../Images/qqplot3.png')
qqnorm(e3$residuals)
qqline(e3$residuals)
dev.off()

png(filename='../Images/hist3.png')
hist(e3$residuals)
dev.off()

png(filename='../Images/standard3.png')
plot(ejercicio2$sales, rstandard(e3), ylab = "Residuos estandarizados", xlab = "sales")
abline(h =0, lty = 2)
dev.off()

e4 <- lm(sales~TV+radio, data=ejercicio2)
summary(e4)

mean(e4$residuals)
sum(e4$residuals)

dwtest(e4)

png(filename='../Images/qqplot4.png')
qqnorm(e4$residuals)
qqline(e4$residuals)
dev.off()

png(filename='../Images/hist4.png')
hist(e4$residuals)
dev.off()

png(filename='../Images/standard4.png')
plot(ejercicio2$sales, rstandard(e4), ylab = "Residuos estandarizados", xlab = "sales")
abline(h =0, lty = 2)
dev.off()

e5 <- lm(sales~TV, data=ejercicio2)
summary(e5)

mean(e5$residuals)
sum(e5$residuals)

dwtest(e5)

png(filename='../Images/qqplot5.png')
qqnorm(e5$residuals)
qqline(e5$residuals)
dev.off()

png(filename='../Images/hist5.png')
hist(e5$residuals)
dev.off()

png(filename='../Images/standard5.png')
plot(ejercicio2$sales, rstandard(e5), ylab = "Residuos estandarizados", xlab = "sales")
abline(h =0, lty = 2)
dev.off()

max_stars_repo_path: Second Assigment/Code/code.r

max_issues_repo_path: Second Assigment/Code/code.r

max_forks_repo_path: Second Assigment/Code/code.r


text:
s3 = import('s3')
test = s3$test

max_stars_repo_name: pooranis/rmodules

max_issues_repo_name: pooranis/rmodules

max_forks_repo_name: pooranis/rmodules


text:
# Prepare data for display on data table
# data table must display age adjusted cancer incidence rates of most recent 5 year average
# René Dario Herrera
# University of Arizona Cancer Center
# renedherrera at email dot arizona dot edu
# May 2021

# set up ----
# load packages ----
library(here)
library(tidyverse)
library(leaflet)
library(tigris)
library(sf)
library(ggplot2)

# load county shapefiles for AZ ----
county_spatial <- counties(
  state = "04"
)

# inspect
glimpse(county_spatial)
class(county_spatial)
st_crs(county_spatial)

# read data ----
incidence_table <- read_rds("communication/shiny_apps/dashboard_incidence_thematic_map/data/incidence_table.rds")

# inspect
glimpse(incidence_table)

# filter incidence table to show only counties
incidence_table <- incidence_table %>%
  filter(!(area == "US" | area == "Catchment" | area == "Arizona"),
         site == "All Cancer Sites Combined",
         sex == "Male and Female",
         race == "All Races")

# merge
incidence_spatial <- geo_join(
  spatial_data = county_spatial,
  data_frame = incidence_table,
  by_sp = "NAME",
  by_df = "area"
)

#inspect
class(incidence_spatial)
glimpse(incidence_spatial)
st_crs(incidence_spatial)

# transform CRS for use with Leaflet
incidence_spatial <- st_transform(incidence_spatial, '+proj=longlat +datum=WGS84')

# save to disk for use in shiny app
write_rds(incidence_spatial, "communication/shiny_apps/dashboard_incidence_thematic_map/data/incidence_spatial.rds")

# need to know bins to set color palette for leaflet choropleth
ggplot(data = incidence_table) +
  geom_histogram(mapping = aes(x = age_adjusted_rate), bins = 5)

# bins from histogram plot
bins <- c(0, 250, 300, 350, 400, 450)

# set palette
pal <- colorBin("viridis",
                domain = incidence_spatial$age_adjusted_rate,
                bins = bins)

# set labels
labels <- sprintf(
  "%s
%g per 100,000",
  incidence_spatial$NAME, incidence_spatial$age_adjusted_rate
) %>% lapply(htmltools::HTML)

# generate leaflet visualization
leaflet(data = incidence_spatial) %>%
  addProviderTiles("Stamen.TonerLite") %>%
  addPolygons(fillColor = ~pal(age_adjusted_rate),
              color = "white",
              fillOpacity = .7,
              weight = 1,
              label = labels,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "16px",
                direction = "auto"
              )) %>%
  addLegend(pal = pal, values = ~age_adjusted_rate, opacity = 0.7, title = NULL,
            position = "bottomright")


text:
#!/usr/bin/env r

usePG <- TRUE

if (usePG) {
    cat("Using Pg\n")
    tempdb <- "pgdatetime"
    system(paste("createdb", tempdb))   # create a temp database

    suppressMessages(library(RPostgreSQL))
    drv <- dbDriver("PostgreSQL")
    con <- dbConnect(drv, dbname=tempdb)

} else {
    cat("Using SQLite\n")
    tempdb <- "/tmp/tempdb.sqlite"          # assumed to not exist
    suppressMessages(library(RSQLite))
    drv <- dbDriver("SQLite")
    if (file.exists(tempdb))
        unlink(tempdb)
    con <- dbConnect(drv, tempdb)
}

sql <- "create table foo (i integer, r real, t text)"
res <- dbSendQuery(con, sql)
cat("Created table\n")

i <- as.integer(11)
r <- as.numeric(22.22)
txt <- as.character("blim blom")

sql <- paste("insert into foo ",
             "values (",
             i, ",",
             r, ", '",
             txt, "') ", sep="")
res <- dbSendQuery(con, sql)
cat("Wrote values\n")

df <- dbReadTable(con, "foo")
cat("Read values\n")

## now test the types of the colums we got
stopifnot( class(df[,1]) == "integer" )
stopifnot( class(df[,2]) == "numeric" )
stopifnot( class(df[,3]) == "character" )
cat("GOOD -- all types are as expected\n")

## and test the values
stopifnot( identical( df[1,1], i))
stopifnot( identical( df[1,2], r))
stopifnot( identical( df[1,3], txt))
cat("GOOD -- all values are as expected\n")

if (usePG) {
    dbDisconnect(con)
    dbUnloadDriver(drv)
    system(paste("dropdb", tempdb))
} else {
    dbDisconnect(con)
    dbUnloadDriver(drv)
    unlink(tempdb)
}

cat("DONE\n")

max_stars_repo_name: geaviation/shiny

max_issues_repo_name: geaviation/shiny

max_forks_repo_name: geaviation/shiny


text:
halfhelmert_contrasts = function(...){
	contr.helmert(...)*.5
}


text:
#Assign spatial products to Ramirez synthesis data.
#Climate from worldclim2. N-dep from NADP.
#These are based on functions written by Colin.
rm(list = ls())
source('paths.r')
source('paths_fall2019.r')
source('NEFI_functions/extract_ndep.r')
source('NEFI_functions/worldclim2_grab.r')
source('NEFI_functions/arid_extract.r')
source('NEFI_functions/extract_pH.r')
source('NEFI_functions/extract_npp.r')
source('NEFI_functions/extract_ndep_global.r')
source('NEFI_functions/extract_EM.r')

# set output path
output.path <- ramirez_clean_map.path

#load 16S meta data file with lat/long
d <- readRDS(ramirez_raw_mapping_and_abundance.path)

# removing taxonomic data
#metadata.cols <- colnames(d)[!grepl("^[phylum|class|order|family|genus]", colnames(d))]
metadata.cols <- colnames(d)[!grepl("^[dpcofgs]__[bacteria]", colnames(d))]
metadata.cols <- metadata.cols[!grepl("^[dpcofgs]__[vu]", metadata.cols)]
#d <- d[,..metadata.cols]
d <- d[,metadata.cols]

#rownames(d) <- d$sampleID
# d$latitude <- as.numeric(d$latitude.x)
# d$longitude <- as.numeric(d$longitude.x)
#d$altitude.elevation..meter. <- as.numeric(d$altitude.elevation..meter.)
# remove studies without lat/lon (would be cool to get this netherlands one back later)
d <- d[!is.na(d$latitude),]
#extract worldclim.

# on pecan:
# clim <- worldclim2_grab(d$latitude,d$longitude,d$altitude.elevation..meter.)
# clim$aridity <- arid_extract(d$latitude,d$longitude)
# on SCC:
clim <- worldclim2_grab(d$latitude,d$longitude,
                        #d$altitude.elevation..meter.,
                        worldclim2_folder = "/projectnb/talbot-lab-data/spatial_raster_data/WorldClim2/")
clim$aridity <- arid_extract(d$latitude,d$longitude, folder = "/projectnb/talbot-lab-data/spatial_raster_data/Global_Aridity/aridity/")

#extract N deposition. Half of these are NA because our Ndep product only covers the United States.
#ignore the warnings.
#ndep <- extract_ndep(d$longitude,d$latitude, folder = "/projectnb/talbot-lab-data/spatial_raster_data/CASTNET_Ndep/")

#extract pH
#ph_estim <- extract_pH(d$latitude,d$longitude, folder = "/projectnb/talbot-lab-data/spatial_raster_data/SoilGrids_uncertainty/")

#NPP <- extract_npp(d$latitude, d$longitude)
NPP <- extract_npp(d$latitude, d$longitude, path="/projectnb/talbot-lab-data/NEFI_data/covariate_data/NPP/w001001.adf")

# from Ackerman et al.: https://conservancy.umn.edu/handle/11299/197613
#ndep.glob <- extract.ndep.global(d$latitude, d$longitude)

# extract % basal area of ecto trees (using the not-masked raster)
relEM <- extract_EM(d$latitude, d$longitude, path = "/projectnb/talbot-lab-data/NEFI_data/covariate_data/EM_AM_maps/em_1degree.nc")

#put together all spatial products
spatial <- cbind(clim, NPP, relEM)
#spatial <- cbind(ndep, ph_estim,clim, NPP, ndep.glob, relEM)

#remove any columns from previous spatial extractions, update with new extraction.
d <- cbind(d,spatial)

# add NPP for a couple studies. these are the values if the longitudes was .1 different
#extract_npp(latitude = 40.77552, longitude = -74.01,path="/projectnb/talbot-lab-data/NEFI_data/covariate_data/NPP/w001001.adf")
#extract_npp(latitude = 50.0, longitude = -5,path="/projectnb/talbot-lab-data/NEFI_data/covariate_data/NPP/w001001.adf")
d[round(d$latitude, 1)==40.8 & round(d$longitude, 1)==-74.0|round(d$longitude, 1)==-73.9,]$NPP <- .736
d[round(d$latitude, 1)==50.0 & round(d$longitude, 1)==-5.2,]$NPP <- 0.974

#save
saveRDS(d,output.path)


text:
library(RJSONIO)
langUrl <- "http://rosettacode.org/mw/api.php?action=query&format=json&cmtitle=Category:Solutions_by_Programming_Language&list=categorymembers&cmlimit=500"

languages <- fromJSON(langUrl)$query$categorymembers
languages <- sapply(languages, function(x) sub("Category:", "", x$title))

# fails if there are more than 500 users per language
user <- function (lang) {
  userBaseUrl <- "http://rosettacode.org/mw/api.php?action=query&format=json&list=categorymembers&cmlimit=500&cmtitle=Category:"
  userUrl <- paste(userBaseUrl, URLencode(paste(lang, " User", sep="")),sep="")
  length(fromJSON(userUrl)$query$categorymembers)
}

users <- sapply(languages, user)
head(sort(users, decreasing=TRUE),15)

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
library(XML)
library(RCurl)
library(httr)

rm(list=ls(all.names=TRUE))

fpUrl <- 'http://m.click108.com.tw/astro/index.php?astroNum='

zodiac <- c('Aries','Taurus','Gemini','Cancer','Leo','Virgo','Libra','Scorpio','Sagittarius'
            ,'Capricorn','Aquarius','Pisces')
factor <- c('Love','Finance','Career')

data <- matrix(0,3,12)
rownames(data) <- factor
colnames(data) <- zodiac

xmlPathLove <- '//div[@id="astroDailyScore_love"]/@style'
xmlPathCareer <- '//div[@id="astroDailyScore_career"]/@style'
xmlPathMoney <- '//div[@id="astroDailyScore_money"]/@style'

for( i in 1:12){
  url <- paste0(fpUrl,i)

  tmp <- getURL(url,encoding='big5')

  xmldoc <- htmlParse(tmp)

  love <- xpathSApply(xmldoc,xmlPathLove)
  money <- xpathSApply(xmldoc,xmlPathMoney)
  career <-xpathSApply(xmldoc,xmlPathCareer)

  score <- c( as.numeric(substr(love,84,84)),
              as.numeric(substr(money,85,85)),
              as.numeric(substr(career,84,84))
              )
  data[,i] <- score
}

date <- Sys.Date()

filename <- paste0('data\\httpm.click108.com\\httpm.click108.com',date,'.csv')

write.csv(data,filename)

max_stars_repo_name: ivorydaae/TuesClass

max_issues_repo_name: ivorydaae/TuesClass

max_forks_repo_name: ivorydaae/TuesClass


text:
#BoxPlot (or WhiskerPlot)
# ЧАСТЬ 1: делаем data.frame
	# шаг-1. вчитываем таблицу. делаем из нее датафрейм.
MDepths <- read.csv("Depths.csv", header=TRUE, sep = ",")
	# шаг-2. чистим датафрейм от NA значений
MDepths <- na.omit(MDepths)
row.has.na <- apply(MDF, 1, function(x){any(is.na(x))}) # проверяем, удалил ли все NA
sum(row.has.na) # суммируем все NA, должно получиться: [1] 0
head(MDepths) # смотрим очищенный датафрейм. теперь с ним работаем.

# ЧАСТЬ 2: делаем диаграмму распределения данных "whisker boxplot"
	# шаг-3. создаем лист названий профилей	(что будет отображаться по оси X)
profile_names <- paste(c("profile"), seq(1:25), sep="")
	# шаг-4. строим боксплот по указанным аргументам (здесь: 25 профилей)
p<- boxplot(MDepths$profile1, MDepths$profile2, MDepths$profile3, MDepths$profile4, MDepths$profile5, MDepths$profile6,MDepths$profile7, MDepths$profile8, MDepths$profile9, MDepths$profile10, MDepths$profile11, MDepths$profile12, MDepths$profile13, MDepths$profile14, MDepths$profile15, MDepths$profile16, MDepths$profile17, MDepths$profile18, MDepths$profile19, MDepths$profile20, MDepths$profile21, MDepths$profile22, MDepths$profile23, MDepths$profile24, MDepths$profile25,
	main = "Mariana Trench Depths Boxplot",
	outline = TRUE,
	outlier.color = "seagreen", outlier.shape = 8, outlier.size = 2,
#	xlab = "Profiles",
	ylab = "Depths",
	las = 2,
	col = viridis(25, alpha=.2),
	names = profile_names)
p

max_stars_repo_name: paulinelemenkova/R-Boxplot-Whiskerplot-

max_issues_repo_name: paulinelemenkova/R-Boxplot-Whiskerplot-

max_forks_repo_name: paulinelemenkova/R-Boxplot-Whiskerplot-


text:
type 'MErr' {
	array Strs {
    	pstring;
    };
};

/*
resource 'MErr' (16091, "File Type Wrong Length") {
	{
	"The file type is an incorrect length.  I"
	"t must be exactly 4 characters long.",
	"This is a second string"
	}
};
*/

max_stars_repo_name: mickaelcala/TNTbasic

max_issues_repo_name: mctully/tntbasic

max_forks_repo_name: mctully/tntbasic


text:
#  James Rekow

rSubgraph = function(g, v){

  #  ARGS: g - input graph
  #        v - number of vertices from g to include in the random subgraph
  #
  #  RETRURNS: rSub - a connected subgraph of g with v vertices and all of the edges between them

  i = 0
  numVertices = vcount(g)
  rSub = induced.subgraph(g, sample(numVertices, v))

  while(!is.connected(rSub)){

    rSub = induced.subgraph(g, sample(numVertices, v))
    i = i + 1

  } #  end while
  print(i)
  return(rSub)

} #  end rSubgraph function


text:
get_sequence <- function (object, observed = FALSE, length_out = 30) {
  assert_that(is.flag(observed))
  assert_that(is.count(length_out))
  assert_that(length_out > 0)

  UseMethod("get_sequence", object)
}

get_sequence.continuous <- function(object, observed = FALSE, length_out = 30) {

  if(length_out == 1)
    return (get_mean(object, observed = observed))

  if (observed) {
    if(length_out >= length(object$obs))
      return (object$obs)
    x <- seq(from = 1, to = length(object$obs), length.out = length_out)
    x <- as.integer(round(x))
    return (object$obs[x])
  }
  x <- seq(from = get_min(object), to = get_max(object), length.out = length_out)
  x <- unique(x)
  return (x)
}

get_sequence.categorical <- function (object, observed = FALSE, length_out = 30) {

  if(length_out == 1)
    return (get_min(object, observed = observed))

  if(observed) {
    if (length_out >= length(object$obs))
      return (object$obs)
    return (object$obs[1:length_out])
  }
  if (length_out >= length(object$levels))
    return (object$levels)
  return (object$levels[1:length_out])
}

get_sequence.vinteger <- function(object, observed = FALSE, length_out = 30) {
  if (observed) {
    if(length_out >= length(object$obs))
      return (object$obs)
    x <- seq(from = 1, to = length(object$obs), length.out = length_out)
    x <- as.integer(round(x))
    return (object$obs[x])
  }
  x <- seq(from = get_min(object), to = get_max(object), length.out = length_out)
  x <- unique(as.integer(round(x)))
  return (x)
}

max_stars_repo_name: poissonconsulting/jaggernaut

max_issues_repo_name: poissonconsulting/jaggernaut

max_forks_repo_name: poissonconsulting/jaggernaut


text:
set.seed(7)
#TAMANHO
size <- 800

#INDICES
indices <- 1:size

#CATEGORIAS
genero <- c("Masculino",
            "Feminino")

formacao <- c("Médio",
              "Superior incompleto",
              "Superior completo",
              "Pós-graducação",
              "Mestrado",
              "Doutorado")

religiao <- c("Nd",
              "Cristianismo",
              "Budismo",
              "Hinduísmo",
              "Islamismo")

data <- data.frame(indice = indices,
                   genero = sample(genero, size, TRUE),
                   formacao = sample(formacao, size, TRUE),
                   religiao = sample(religiao, size, TRUE),
                   peso = round(rnorm(size, 80, 10), digits = 1),
                   altura = round(rnorm(size, 1.75, c(0.05)), digits = 2),
                   salario = round(rnorm(size, 2000, 500), digits = 2),
                   filho = floor(runif(size, 0, 5)),
                   automovel =  floor(runif(size, 0, 3)))

#Salvando dataframe
arquivo <- "r_file" # Define nome do arquivo
caminho <- paste(getwd(), arquivo, sep="/") # Define diretorio
write.table(data, file = caminho, sep = ";") # Salva dataframe

#GERANDO INFORMAÇÕES

#Histogramas
hist(data[,5], main = "Histograma de peso")
hist(data[,6], main = "Histograma de altura")

#Gráfico de peso x altura
plot(x = data[,5], y = data[,6])

#Gráfico de salário x filhos
plot(x = data[,7], y = data[,8])

# Split por categoria
data_genero <- split(data[,5:9], data$genero)
data_formacao <- split(data[,5:9], data$formacao)
data_religiao <- split(data[,5:9], data$religiao)

# Médias por genero
medias_por_genero <- data.frame(
  media_geral = sapply(data[, 5:9], mean),
  media_feminino = sapply(data_genero$Feminino[,1:5], mean),
  media_masculino = sapply(data_genero$Masculino[,1:5], mean))

#Visualização
View(medias_por_genero)

#Soma por religião
filhos_e_automoveis_por_religiao <- data.frame(
    soma_total = sapply(data[,8:9], sum),
    soma_budismo = sapply(data_religiao$Budismo[,4:5], sum),
    soma_cristianismo = sapply(data_religiao$Cristianismo[,4:5], sum),
    soma_nd = sapply(data_religiao$Nd[,4:5], sum),
    soma_hinduismo = sapply(data_religiao$Hinduísmo[,4:5], sum),
    soma_islamismo = sapply(data_religiao$Islamismo[,4:5], sum))

View(filhos_e_automoveis_por_religiao)

# Quantidade por categoria
categorias <- apply(data[2:4], 2, table)
categorias$genero
categorias$religiao
categorias$formacao

max_stars_repo_name: abrantesasf/ApoemaTraining

max_issues_repo_name: abrantesasf/ApoemaTraining

max_forks_repo_name: abrantesasf/ApoemaTraining


text:
StanDDM_NCEN_St <- function(simulation, modal, data, control, seed, warmup, num_cores, iter){

    model_name <- paste('StanDDM_NCEN_St', modal, sep="") #modal assigned in main
    dir.create(paste(model_name, '/', sep=''), showWarnings = FALSE)
    cat('\n///////////////Computing model: ', model_name, '///////////////\n')

    if(simulation){
        cat('\nCreating random parameters and simulating Data...\n')
        fakeParams <- makeFakeParams(nsub = 9, include = c('st'))    #'include' is model dependent!
        fakedat <- simulDat(fakeParams)
        fakedat <- fakedat[[1]]
        data <- fake_data_processing(fakedat$data)
    }

    model_definition <-
    "
    data {
      int<lower=1> N;      // Number of subjects
      int<lower=0> Nu_max; // Max (across subjects) number of upper boundary responses
      int<lower=0> Nl_max; // Max (across subjects) number of lower boundary responses
      int<lower=0> Nu[N];  // Number of upper boundary responses for each subj
      int<lower=0> Nl[N];  // Number of lower boundary responses for each subj
      real RTu[N, Nu_max];  // upper boundary response times
      real RTl[N, Nl_max];  // lower boundary response times
      real minRT[N];       // minimum RT for each subject of the observed data
      real RTbound;        // lower bound or RT across all subjects (e.g., 0.1 second)
    }

    parameters {
    //Group level parameters in vecotized form, one for each parameter.
      vector[4] mu_mu; //group level mean means
      vector<lower=0>[4] mu_sigma; //group level mean variances

      vector[4] sigma_mu; //group level variance means
      vector<lower=0>[4] sigma_sigma; //group level variance variances

      // Subject-level raw parameters (for Matt trick)
      vector[N] alpha_mu_pr;
      vector[N] beta_mu_pr;
      vector[N] delta_mu_pr;

      vector[N] tau_mu_pr;
      vector[N] tau_sigma_pr;

      // Trial-level raw parameters (for Matt trick)
      vector[(sum(Nu)+sum(Nl))] tau_pr;
    }

    transformed parameters {
      // Transform  raw parameters (primes)

      // Subject-level declarations: MEANS
      vector<lower=0>[N]         alpha_mu; // boundary separation mean subject level
      vector<lower=0, upper=1>[N] beta_mu;  // initial bias mean subject level
      vector<lower=0>[N]         delta_mu; // drift rate mean subject level
      vector<lower=RTbound, upper=max(minRT)>[N] tau_mu; // nondecision time mean subject level

      // Subject-level declarations: VARIANCES
      vector<lower=0>[N]         tau_sigma; // drift rate variance subject level

      // Trial-level declarations
      vector[(sum(Nu)+sum(Nl))] tau;

    //-------------------------------------

      //Subject-level MEANS transformations
      for (i in 1:N) {
        beta_mu[i] = Phi_approx(mu_mu[2] + mu_sigma[2] * beta_mu_pr[i]);
        tau_mu[i]  = Phi_approx(mu_mu[4] + mu_sigma[4] * tau_mu_pr[i]) * (minRT[N]-RTbound) + RTbound;
      }

      alpha_mu = exp(mu_mu[1] + mu_sigma[1] * alpha_mu_pr); //reparametrization as in  Gelman manual second ed pg 313 and Kruschkes manual pg. 281
      delta_mu = exp(mu_mu[3] + mu_sigma[3] * delta_mu_pr);

      //Subject-level VARIANCES transformations
      tau_sigma = exp(sigma_mu[4] + sigma_sigma[4] * tau_sigma_pr);//already vectorized for subject level in declaration

    //Trial-level transformations: MUST BE SLICED AS IN MODEL BLOCK

        for (i in 1:N) {//begin subject loop

            //loop-scoped definitions for index calcs
            int NuL; int NuU; int NlL; int NlU;

            NuL = sum(Nu[1:i]) -Nu[i] +1;
            NuU = sum(Nu[1:i]);
            NlL = sum(Nl[1:i]) -Nl[i] +1 + sum(Nu);
            NlU = sum(Nl[1:i]) + sum(Nu);

    		tau[NuL:NuU] = Phi_approx(tau_mu[i] + tau_sigma[i] * tau_pr[NuL:NuU])* (minRT[i]-RTbound) + RTbound;
    		tau[NlL:NlU] = Phi_approx(tau_mu[i] + tau_sigma[i] * tau_pr[NlL:NlU])* (minRT[i]-RTbound) + RTbound;

        }//end subject loop

    }

    model {
        // Group level parameters (all vectorized!)
        mu_mu  ~ normal(0, 1); //Group mean mean
        mu_sigma ~ cauchy(0, 5); //Group mean variance

        sigma_mu ~ cauchy(0, 5); //Group variance mean
        sigma_sigma ~ cauchy(0, 5); //Group variance variance

        // Subject-level sampling statements: MEANS
        alpha_mu_pr ~ normal(0, 1);
        beta_mu_pr ~ normal(0, 1);
        delta_mu_pr ~ normal(0, 1);
        tau_mu_pr ~ normal(0, 1);

        // Subject-level sampling statements: VARIANCES
        tau_sigma_pr ~ cauchy(0, 5);

        // Trial-level sampling statements: PRIMES
        //alpha_pr ~ normal(0, 1);
        //beta_pr  ~ normal(0, 1);
        //delta_pr ~ normal(0, 1);
        tau_pr   ~ normal(0, 1);

        // Begin subject loop
        for (i in 1:N) {

            int NuL; int NuU; int NlL; int NlU;

            NuL = sum(Nu[1:i]) -Nu[i] +1;
            NuU = sum(Nu[1:i]);

            RTu[i, 1:Nu[i]] ~ wiener(alpha_mu[i], tau[NuL:NuU], beta_mu[i], delta_mu[i]);

            NlL = sum(Nl[1:i]) -Nl[i] +1 + sum(Nu);
            NlU = sum(Nl[1:i]) + sum(Nu);

            RTl[i, 1:Nl[i]] ~ wiener(alpha_mu[i], tau[NlL:NlU], 1-beta_mu[i], -delta_mu[i]);

        } // end of subject loop
    }// end model block

    generated quantities {

        vector[(sum(Nu)+sum(Nl))] log_lik;

        for (i in 1:N) {

            int NuL; int NuU; int NlL; int NlU;

            NuL = sum(Nu[1:i]) -Nu[i] +1;
            NuU = sum(Nu[1:i]);

            for (j in NuL:NuU){
                log_lik[j] = wiener_lpdf(RTu[i, 1:Nu[i]] | alpha_mu[i], tau[j], beta_mu[i], delta_mu[i]);
            }

            NlL = sum(Nl[1:i]) -Nl[i] +1 + sum(Nu);
            NlU = sum(Nl[1:i]) + sum(Nu);

            for (j in NlL:NlU){
                log_lik[j] = wiener_lpdf(RTl[i, 1:Nl[i]] | alpha_mu[i], tau[j], 1-beta_mu[i], -delta_mu[i]);
            }
        }
    }
    "

    inits_fixed <- c(0.5, 0.5, 0.5, 0.15)
    init <- lapply(1:num_cores, function(i) {
        list(
            mu_mu     = c(log(inits_fixed[1]), qnorm(inits_fixed[2]), log(inits_fixed[3]), qnorm(inits_fixed[4])),
            mu_sigma    = c(1.0, 1.0, 1.0, 1.0),
            sigma_mu     = c((inits_fixed[1]), qnorm(inits_fixed[2]), log(inits_fixed[3]), qnorm(inits_fixed[4])),
            sigma_sigma    = c(0.1, 0.1, 0.1, 0.1),
            alpha_mu_pr = rep(log(inits_fixed[1]), data$forstan$N),
            beta_mu_pr  = rep(qnorm(inits_fixed[2]), data$forstan$N),
            tau_mu_pr   = rep(qnorm(inits_fixed[4]), data$forstan$N),
            tau_sigma_pr   = rep(0.1, data$forstan$N),
            delta_mu_pr   = rep(qnorm(inits_fixed[3]), data$forstan$N),
            tau_pr = rep(log(inits_fixed[4]), sum(data$forstan$Nl)+sum(data$forstan$Nu))
        )
    })

    pars <- c("tau", 'tau_sigma', "alpha_mu", "beta_mu", "delta_mu", "tau_mu", 'log_lik')

    cat('\nCompiling and fitting the model...\n')
    #-----------------------------------------------------------------
        fit <- stan(model_code = model_definition,
                    seed = seed,
                    data = data$forstan,
                    pars = pars,
                    warmup = warmup,
                    cores = num_cores,
                    control = control,
                    init = init,
                    iter = iter,
                    chains = num_cores)

    cat('\nSaving Stanfit object...\n')
    #-----------------------------------------------------------------
        saveit(fit = fit,
               string = model_name,
               file=paste(model_name, '/', model_name, '.RData', sep = ''))

    cat('\nExtracting and saving parameters...\n')
    #-----------------------------------------------------------------
        extracted_params <- parameter_extraction(
            fit,
            numsub=data$forstan$N,
            names=data$forsim$names)

        names(extracted_params) <- c('st','a', 'z', 'v', 't') #dependent on model
        write.csv(extracted_params, file = paste(model_name, '/',model_name,".csv", sep = ''), row.names=FALSE)

        if(simulation){ #needs model_name, else that would be in main
            write.csv(fakeParams, file = paste(model_name, '/',model_name,"_template_params.csv", sep = ''), row.names=FALSE)

        }
        #TODO: name parsing?

    cat(paste('\nSimulating', data$forstan$N, 'subjects with fitted parameters:\n'))
    #-----------------------------------------------------------------
        sims <- simulDat(extracted_params)
        sims <- sims[[1]]

    # cat('\nSimulated incorrect responses:', mean(mean(cor$value==0)),'%',
    #     '\nSimulated correct responses:', mean(mean(cor$value==1)), '%', '\n')

    cat(paste('\nSaving Plots...\n'))
    #-----------------------------------------------------------------
        models_plots(experim_dat = data$forsim$rawdat,
                     simul_dat = sims,
                     model_name = model_name)

        fit_quality(experim_dat = data$forsim$rawdat,
                 simul_dat = sims,
                 model_name = model_name)

    if(simulation){
        param_recovery(fit, data$forstan$N, fakeParams, model_name = model_name)
    }

    cat(paste('\nComputing/plotting model diagnostics:\n'))
    #-----------------------------------------------------------------
        model_diagnostics(fit)

    cat(paste('\nFor Reproducibility:\n'))
    #-----------------------------------------------------------------
        print(devtools::session_info())
        cat(paste('Seed:', seed, sep = ' '))

    cat('\n///////////////', model_name, 'has been computed.///////////////\n')
}# end of function


text:
# loads the library urca, which contains the command ur.df()
# and applies Unit-Root Testing based on Dickey-Fuller test
# install.packages(c('fGarch', 'urca','tseries','MASS'))
library(urca)
library(tseries)
library(MASS)
library(fGarch)

# Importing data into R can be carried out in various ways. Below, the command read.table is used:
data_all = read.table('/home/socital1/Desktop/george/AUEB/TimeSeries/time_series_in_r/Data-Assignment-Feb-2018.txt')

xdata_cols <- c('V5','V6', 'V7','V8','V9','V10','V11','V12','V13','V14','V15','V16','V17','V18')

xdata <- data_all[c(xdata_cols)]

head(xdata)

y2 <- data_all$V2

j2 = ts(y2, frequency=12, start = c(1991, 1), end = c(2004, 12))
xfactors=ts(xdata, frequency=12,  start = c(1991, 1), end = c(2004, 12))

train_data = data_all[1: 168,]
test_data = data_all[169:180,]

attach(train_data)

fitnull <- lm(V2 ~ 1)

stepSR <- step(fitnull,
               scope=list(lower = ~ 1,
                          upper = ~ V5 + V6 + V7 + V8 + V9 + V10 + V11 + V12 + V13 + V14 + V15 + V16 + V17 + V18),
               direction="both",
               data=train_data)

stepSR$anova

summary(stepSR)

detach(train_data)

stepSRresidualsTS <- ts(stepSR$residuals,
                        frequency=12,
                        start = c(1991,1),
                        end = c(2004, 12))
stepSRresidualsTS

par(mfrow=c(3,2))

# set up the graphics
acf(ts(stepSRresidualsTS,freq=1), 48, main="ACF of residuals")
pacf(ts(stepSRresidualsTS,freq=1), 48, main="PACF of residuals")
acf(ts(stepSRresidualsTS^2,freq=1), 48, main="ACF of squared residuals")
pacf(ts(stepSRresidualsTS^2,freq=1), 48, main="PACF of squared residuals")
qqnorm(stepSRresidualsTS,main="Normal QQplot of residuals")
qqline(stepSRresidualsTS)
plot(j2, type="l", col='red', lwd=1, main="Time Series plot of Y2", ylab="Monthly returns")

ar2fit <- arima(stepSRresidualsTS, order=c(2,0,0))
ar2fit

par(mfrow=c(3,2))
ar2fit_residualsTS <- ts(ar2fit$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(ar2fit_residualsTS,freq=1), 48, main="ACF on residuals of AR(2)")
pacf(ts(ar2fit_residualsTS,freq=1), 48, main="PACF on residuals of AR(2)")
acf(ts(ar2fit_residualsTS^2,freq=1), 48, main="ACF on squared residuals of AR(2)")
pacf(ts(ar2fit_residualsTS^2,freq=1), 48, main="PACF on squared residuals of AR(2)")
qqnorm(ar2fit_residualsTS,main="Normal QQplot on residuals of AR(2)")
qqline(ar2fit_residualsTS)

Box.test(ar2fit_residualsTS, lag=12, type="Ljung")
Box.test(ar2fit_residualsTS^2, lag=12, type="Ljung")

ma2fit <- arima(stepSRresidualsTS, order=c(0,0,2))
ma2fit

par(mfrow=c(3,2))
ma2fit_residualsTS <- ts(ma2fit$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(ma2fit_residualsTS,freq=1), 48, main="ACF on residuals of MA(2)")
pacf(ts(ma2fit_residualsTS,freq=1), 48, main="PACF on residuals of MA(2)")
acf(ts(ma2fit_residualsTS^2,freq=1), 48, main="ACF on squared residuals of MA(2)")
pacf(ts(ma2fit_residualsTS^2,freq=1), 48, main="PACF on squared residuals of MA(2)")
qqnorm(ma2fit_residualsTS,main="Normal QQplot on residuals of MA(2)")
qqline(ma2fit_residualsTS)

Box.test(ma2fit_residualsTS, lag=12, type="Ljung")
Box.test(ma2fit_residualsTS^2, lag=12, type="Ljung")

arma22fit <- arima(stepSRresidualsTS, order=c(2, 0,2))
arma22fit

par(mfrow=c(3,2))
arma22fit_residualsTS <- ts(arma22fit$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(arma22fit_residualsTS,freq=1), 48, main="ACF on residuals of ARMA(2, 2)")
pacf(ts(arma22fit_residualsTS,freq=1), 48, main="PACF on residuals of ARMA(2, 2)")
acf(ts(arma22fit_residualsTS^2,freq=1), 48, main="ACF on squared residuals of ARMA(2, 2)")
pacf(ts(arma22fit_residualsTS^2,freq=1), 48, main="PACF on squared residuals of ARMA(2, 2)")
qqnorm(arma22fit_residualsTS, main="Normal QQplot on residuals of ARMA")
qqline(arma22fit_residualsTS)

Box.test(arma22fit_residualsTS, lag=12, type="Ljung")
Box.test(arma22fit_residualsTS^2, lag=12, type="Ljung")

# First, we will fit an ARCH(1) model:
m1arch <- garchFit(~garch(1,0), data=stepSR$residuals, trace=F)
# trace = F reduces the summary
summary(m1arch)

predict(m1arch, 6)

m1arch_student=garchFit(~garch(1,0),data=stepSR$residuals, cond.dist="std",trace=F)
summary(m1arch_student)

m2garch=garchFit(~garch(1,1), data=stepSR$residuals,trace=F)
summary(m2garch)

arma11_garch11fit <- garchFit(formula= ~arma(1,1) + garch(1,1),
                              data=stepSR$residuals,
                              trace=F)
summary(arma11_garch11fit)

forecast <- predict(arma11_garch11fit, 12)

# plot of forecasts with 1 standard error
par(mfrow=c(1,1))
UL <- forecast$meanForecast + forecast$meanError
LL <- forecast$meanForecast - forecast$meanError

predTS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))
UL_TS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))
LL_TS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))

min_x <- min(stepSRresidualsTS, LL)
max_x <- max(stepSRresidualsTS, UL)

ts.plot(stepSRresidualsTS, predTS)
lines(predTS, col="red", type="o")
lines(UL_TS, col="blue", lty="dashed")
lines(LL_TS, col="blue", lty="dashed")

library(hydroGOF)
MSE <- mse(forecast$meanForecast, forecast$meanError)
MSE

max_issues_repo_name: gperakis/Time-Series-Analysis-Using-R

max_forks_repo_name: gperakis/Time-Series-Analysis-Using-R


text:
fdr.plot <- function(tp, tn, data.set, main="", limit, legx, legy){
    plot(1, type="n", xlab="", ylab="", xlim=c(1, limit), ylim=c(0, 1))
    med.mat <- matrix(data=NA, nrow=limit, ncol=5)
    mtext(side=3,line=0.2,font=1, cex=1,main)
    mtext(side=2,line=2.5,font=1, cex=0.8,"Proportion")
    mtext(side=1,line=2.5,font=1, cex=0.8, "N samples")

   # this needs to be a for loop or nested lapply when the full dataset is collected
   for(i in 2:limit){
        TP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tp))) #/ length(tp)
        FP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tn))) #/ length(tn)

        FN.dist <- sapply(data.set[[i]], function(x) length( setdiff(tp,  (intersect(x,tp)) )) )

        TN.dist <- length(c(tp, tn)) - (TP.dist+FP.dist+FN.dist)

        tpr.dist <- (TP.dist) / (TP.dist + FN.dist)
        fpr.dist <- (FP.dist) / (FP.dist + TN.dist)
        fnr.dist <- (FN.dist) / (TP.dist + FN.dist)
        fdr.dist <- (FP.dist) / (TP.dist + FP.dist )

        acc.dist <- (TP.dist + TN.dist) / (TP.dist + FP.dist + TN.dist + FN.dist)

        segments(i,quantile(tpr.dist, probs=c(0.025, 0.975)[1], na.rm=T), i, quantile(tpr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0,0,1,0.3))
        segments(i,quantile(fpr.dist, probs=c(0.025, 0.975)[1], na.rm=T), i, quantile(fpr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(1,0,0,0.3))
        segments(i,quantile(fdr.dist, probs=c(0.025, 0.975)[1], na.rm=T), i, quantile(fdr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0,0,0,0.3))

        segments(i,quantile(acc.dist, probs=c(0.025, 0.975)[1], na.rm=T), i, quantile(acc.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0.8,0.4,0.1,0.3))

        med.mat[i,] <- c(median(tpr.dist), median(fnr.dist), median(fpr.dist),median(fdr.dist),median(acc.dist))
    }
    points(2:limit,med.mat[2:limit,1],  pch=21, type="b", cex=0.5, col="blue",bg="white") # tpr
    points(2:limit,med.mat[2:limit,3],  pch=21, type="b", cex=0.5, col="red",bg="white") # fpr
    points(2:limit,med.mat[2:limit,4],  pch=21, type="b", cex=0.5   , col="black",bg="white") # fdr
    points(2:limit,med.mat[2:limit,5],  pch=21, type="b", cex=0.5   , col="orange",bg="white") # fdr

    legend(legx, legy, legend=c("tpr","fpr", "fdr", "acc"), pch=21, col=c("blue","red", "black","orange"))
    abline(h=0,lty=2, col="grey")

}

analyze_data <- function(tp,mp,tn,data.set, main="", size, limit) {
    med.mat <- matrix(data=NA, nrow=limit, ncol=5)

    plot(1, type="n", xlab="", ylab="", xlim=c(1, limit), ylim=c(0, size), main=main)

    for(i in 2:limit){
        # this needs to be a for loop or nested lapply when the full dataset is collected
        TP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tp))) #/ length(tp)
        MP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,mp))) #/ length(mp)
        FP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tn))) #/ length(tn)

        FN.dist <- sapply(data.set[[i]], function(x) length( setdiff(c(tp,mp),  (intersect(x,c(tp,mp))) )) )

        TN.dist <- length(c(tp, mp, tn)) - (TP.dist+MP.dist+FP.dist+FN.dist)

        segments(i,quantile(TP.dist, probs=c(0.025, 0.975)[1]), i, quantile(TP.dist, probs=c(0.025, 0.975)[2]), lwd=2, col=rgb(0,0,0,0.3))
        segments(i,quantile(MP.dist, probs=c(0.025, 0.975)[1]), i, quantile(MP.dist, probs=c(0.025, 0.975)[2]), lwd=2, col=rgb(0,0,1,0.3))
        segments(i,quantile(FP.dist, probs=c(0.025, 0.975)[1]), i, quantile(FP.dist, probs=c(0.025, 0.975)[2]), lwd=2, col=rgb(1,0,0,0.3))
        segments(i,quantile(TN.dist, probs=c(0.025, 0.975)[1]), i, quantile(TN.dist, probs=c(0.025, 0.975)[2]), lwd=2, col=rgb(0,1,0,0.3))
        segments(i,quantile(FN.dist, probs=c(0.025, 0.975)[1]), i, quantile(FN.dist, probs=c(0.025, 0.975)[2]), lwd=2, col=rgb(1,0,1,0.3))

        med.mat[i,] <- c(median(TP.dist), median(MP.dist), median(FP.dist),median(TN.dist),median(FN.dist))

    }
    points(2:limit,med.mat[2:limit,1],  pch=21, type="b", col="black",bg="white")
    points(2:limit,med.mat[2:limit,2],  pch=21, type="b", col="blue",bg="white")
    points(2:limit,med.mat[2:limit,3],  pch=21, type="b", col="red",bg="white")
    points(2:limit,med.mat[2:limit,4],  pch=21, type="b", col="green",bg="white")
    points(2:limit,med.mat[2:limit,5],  pch=21, type="b", col="magenta",bg="white")

    legend(1,(0.95 * size), legend=c("TP", "MP","TN","FP","FN"), pch=21, col=c("black","blue","green","red","magenta"))
    abline(h=0,lty=2, col="grey")

}

sens_spec <- function(tp,mp,tn,data.set, main="", limit){
    med.mat <- matrix(data=NA, nrow=limit, ncol=5)

    plot(1, type="n", xlab="", ylab="", xlim=c(1, limit), ylim=c(0, 1), main=main)

    for(i in 2:limit){
        # this needs to be a for loop or nested lapply when the full dataset is collected
        TP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tp))) #/ length(tp)
        MP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,mp))) #/ length(mp)
        FP.dist <- sapply(data.set[[i]], function(x) length(intersect(x,tn))) #/ length(tn)

        FN.dist <- sapply(data.set[[i]], function(x) length( setdiff(c(tp,mp),  (intersect(x,c(tp,mp))) )) )

        TN.dist <- length(c(tp, mp, tn)) - (TP.dist+MP.dist+FP.dist+FN.dist)

        tpr.dist <- (TP.dist) / (TP.dist + FN.dist)
        fpr.dist <- (FP.dist) / (FP.dist + TN.dist)
        fnr.dist <- (FN.dist) / (TP.dist + FN.dist)
        fdr.dist <- (FP.dist) / (TP.dist + FP.dist + 0.1)

        segments(i,quantile(tpr.dist, probs=c(0.025, 0.975)[1]), i, quantile(tpr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0,0,0,0.3))
        segments(i,quantile(fnr.dist, probs=c(0.025, 0.975)[1]), i, quantile(fnr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0,0,1,0.3))
        segments(i,quantile(fpr.dist, probs=c(0.025, 0.975)[1]), i, quantile(fpr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(1,0,0,0.3))
        segments(i,quantile(fdr.dist, probs=c(0.025, 0.975)[1]), i, quantile(fdr.dist, probs=c(0.025, 0.975)[2], na.rm=T), lwd=2, col=rgb(0,1,0,0.3))

        med.mat[i,] <- c(median(tpr.dist), median(fnr.dist), median(fpr.dist),median(fdr.dist),median(fdr.dist))
    }
    points(2:limit,med.mat[2:limit,1],  pch=21, type="b", col="black",bg="white")
    points(2:limit,med.mat[2:limit,2],  pch=21, type="b", col="blue",bg="white")
    points(2:limit,med.mat[2:limit,3],  pch=21, type="b", col="red",bg="white")
    points(2:limit,med.mat[2:limit,4],  pch=21, type="b", col="green",bg="white")

    legend(1,1, legend=c("tpr", "fnr","fpr", "fdr"), pch=21, col=c("black","blue", "red", "green"))
    abline(h=0,lty=2, col="grey")

}

max_stars_repo_name: ggloor/effect

max_issues_repo_name: ggloor/effect

max_forks_repo_name: ggloor/effect


text:
context("resolve mixture")

test_that("100 samples, 2 components, same as glm.nb", {
    f1 = runif(100)
    f2 = 1 - f1
    c1 = stats::rnbinom(n=100, size=10, mu=1000)
    c2 = stats::rnbinom(n=100, size=20, mu=10)
    both = round(f1 * c1 + f2 * c2)

    ctl = MASS::glm.nb(both ~ 0 + f1 + f2, link=identity)

    mf = model.frame(both ~ 0 + f1 + f2)
    terms = attr(mf, "terms")
    x = model.matrix(terms, mf)
    y = model.response(mf)
    res = fit(x, y, lower=c(1e-5,1e-5,-Inf))

    # allow 20/50% deviance
    expect_true(all(abs(res$par[1:2] - coef(ctl)) < c(200, 5)))
})

max_stars_repo_name: mschubert/boundie

max_issues_repo_name: mschubert/boundie

max_forks_repo_name: mschubert/boundie


text:
hello_world <- "Hello, world!"

if (hello_world == "Good morning") {
  print("morning")
} else if (grepl("Hello", hello_world, fixed = TRUE)) {
  print("noon")
} else {
  print("night")
}

time_slot <-
switch(hello_world,
  "Hello, world!" = "nonn",
  "Good morning, world!" = "morning",
  "night"
)

print(time_slot)

max_stars_repo_name: RyoWakabayashi/hello-worlds

max_issues_repo_name: RyoWakabayashi/hello-worlds

max_forks_repo_name: RyoWakabayashi/hello-worlds


text:
gro <- " !\"#$%&'()*+,-./0123456789:<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"

cat(sapply(strsplit(readLines(tail(commandArgs(), n=1)), ";"), function(s) {
  i <- 1
  paste(lapply(strsplit(s[2], "")[[1]], function(t) {
    g <- (nchar(gro) + regexpr(t, gro, fixed = TRUE)[1] - 1 - as.integer(substr(s[1], i, i))) %% nchar(gro) + 1
    i <<- i %% nchar(s[1]) + 1
    substr(gro, g, g)
  }), collapse="")
}), sep="\n")

max_stars_repo_name: IlkhamGaysin/ce-challenges

max_issues_repo_name: IlkhamGaysin/ce-challenges

max_forks_repo_name: IlkhamGaysin/ce-challenges


text:
lfm <- function(T0=.01, TE=.006, TP=.004, TA=.0002, AV=1, Fs=16000)
{
  id = 3
  Ur = 0.0
  Xe = 0.0
  itot = round(T0 * Fs)
  dug = rep(0,itot)

#  avs = 10.0^((AV - 49.5)/20.0)
  avs = AV

  T = 1.0/Fs;
  TE = round(TE * Fs) * T

#
# These adjustments to Tp were not part of the Lin (1990)
# Fortran code, but have been added to silently force Tp
# to have reasonable values.
#
  if (TP >= .99*TE)
    TP = .99*TE
  if (TP / TE < .5)
     TP = .5 * TE
########################

  tb = T0 - TE
  if (TA >= tb)
	  TA = 0.98 * tb
  if (TA == 0.0) {
	  id = 1
  }

  if (id == 3) {
  	Xe0 = 1.0 / TA
	  for (i in 0:199) {
	    xyz = exp(-tb * Xe0)
	    q = xyz - 1.0 + TA * Xe0;
	    dq = -tb * xyz + TA
	    xyz = q / dq
	    Xe = Xe0 - xyz
	    xyz = abs(-xyz / Xe0)
	    if (xyz <= 1.0e-4)
		    break
	    Xe0 = Xe
	  }

	  Xe = Xe0
	  ex = 1.0 / (Xe * TA)
	  Ur = ex * (exp(-Xe * tb) * (1.0 / Xe + tb) - 1.0 / Xe)
  }

  og = pi / TP

  if (TP / TE > 0.5) { # Lin had Te/Tp
    og2 = og * og
    ogte = og * TE
    sd1 = sin(ogte)
    cd1 = cos(ogte)
    eesd = abs(1.0 / sd1)
    ogcd = og * cd1

    tw = (TE - TP) + TA
    al = 200.0
    if (tw/TP >= 0.85)
	    al = 40.0

    for (i in 0:199) {
	    ed1 = exp(al * TE)
      e0 = eesd / ed1
      gg = al * al + og2
      ue = e0 / gg * (ed1 * (al * sd1 - ogcd) + og)
      ut = ue + Ur
      utprime = -TE*ue + (-2.0*al*ue + e0*(ed1 * (TE * (al*sd1 - ogcd) + sd1))) / gg
      xyz = ut / utprime
      aa = al - xyz
      if (abs(xyz/al) <= 1.0e-2)
        break
      al = aa
    }

    if (avs < 0) {
	    e0 = abs(avs) * e0
    } else {
	    e0 = avs / og / (1.0 + exp(al / og*pi)) * gg
    }
  } else {
	  if (avs <= 0.0) {
	    cat("AV must be > 0.0 when TP/TE == 0.5\n")
	    return(dug)
	  }
	  id = 1
	  al = 0.0
	  e0 = avs * og / 2.0
  }

  ite = round(TE * Fs)

  afilter = exp(al / Fs)
  xinput = e0 * afilter * sin(og / Fs)
  bfilter = 2.0 * afilter * cos(og / Fs)
  afilter = afilter * afilter
  y1filter = 0.0
  yout = xinput

  for (i in 1:ite) {
	  dug[i] = yout
	  y2filter = y1filter
	  y1filter = yout
	  yout = bfilter * y1filter - afilter * y2filter
  }

  if (id == 3) {
	  edconst = exp(-tb*Xe)
	  gain = -dug[ite] / (1.0 - edconst)
	  edconst = gain * edconst
	  afilter = exp(-Xe/Fs)
	  y1filter = -gain
	  for (i in (ite+1):itot) {
	    yout = afilter * y1filter
	    y1filter = yout
	    dug[i] = (yout + edconst)
	  }
  } else {
	  for (i in (ite+1):itot)
	    dug[i] = 0.0
  }

  return(dug)
}

max_stars_repo_name: NemoursResearch/FormantTracking

max_issues_repo_name: NemoursResearch/FormantTracking

max_forks_repo_name: NemoursResearch/FormantTracking


text:
#Bands

library(RCurl)
library(rjson)
library(plyr)

#set URL variables
url1 <- "http://www.metal-archives.com/browse/ajax-letter/l/"
url2 <- "/json/1?sEcho=2&iColumns=4&sColumns=&iDisplayStart="
url3 <- "&iDisplayLength=2000&mDataProp_0=0&mDataProp_1=1&mDataProp_2=2&mDataProp_3=3&iSortCol_0=0&sSortDir_0=asc&iSortingCols=1&bSortable_0=true&bSortable_1=true&bSortable_2=true&bSortable_3=false&_=1462052260436"

j<-1
i<-1
		#First pass - establish df
		x <- fromJSON(getURL(paste0(url1, i, url2, j, url3)))

		x <- lapply(x, function(aaData)
		{
		data.frame(matrix(unlist(aaData), ncol=4, byrow=TRUE))
		})

		x <- do.call(rbind, x)
		assign("bands", x)

#loop through all letters
for (i in letters)
	{
		j<-1
		x <- fromJSON(getURL(paste0(url1, i, url2, j, url3)))

		x <- lapply(x, function(aaData)
		{
		data.frame(matrix(unlist(aaData), ncol=4, byrow=TRUE))
		})

		x <- do.call(rbind, x)

		#get count of records for particular letter
		maxrecord <- x[1,1]

		j<-j+499
		maxrecord <- as.numeric(maxrecord)
		maxrecord<- round_any(maxrecord, 500, f=floor)
		bands <-rbind(bands, x)
		#loop over all records for particular letter
		while (j<=maxrecord)
			{
			x <- fromJSON(getURL(paste0(url1, i, url2, j, url3)))

			x <- lapply(x, function(aaData)
			{
			data.frame(matrix(unlist(aaData), ncol=4, byrow=TRUE))
			})

			x <- do.call(rbind, x)
			bands <- rbind(bands, x)
			j<-j+500
			}
	}

#whoops!  Forgot bands whose name begins with a number.  Could expand this to make it more resilient - currently not many bands starting with number

x <-fromJSON(getURL("http://www.metal-archives.com/browse/ajax-letter/l/NBR/json/1?sEcho=1&iColumns=4&sColumns=&iDisplayStart=0&iDisplayLength=500&mDataProp_0=0&mDataProp_1=1&mDataProp_2=2&mDataProp_3=3&iSortCol_0=0&sSortDir_0=asc&iSortingCols=1&bSortable_0=true&bSortable_1=true&bSortable_2=true&bSortable_3=false&_=1463181637758"))
x <- lapply(x, function)
		x <- lapply(x, function(aaData)
		{
		data.frame(matrix(unlist(aaData), ncol=4, byrow=TRUE))
		})

		x <- do.call(rbind, x)
bands <-rbind(bands, x)

max_stars_repo_path: R/metal - working script.r

max_issues_repo_path: R/metal - working script.r

max_forks_repo_path: R/metal - working script.r


text:
library(dplyr)

## requires targets df
region_table <- function(df) {
  by_target_region <- df %>%
    filter(!is.na(Region)) %>%
    group_by(Region) %>%
    summarise(total = n()) %>%
    arrange(-total)

  return(by_target_region)
}


text:
\name{draw_annotation_legend-HeatmapList-method}
\alias{draw_annotation_legend,HeatmapList-method}
\alias{draw_annotation_legend}
\title{
Draw legends for All Annotations
}
\description{
Draw legends for All Annotations
}
\usage{
\S4method{draw_annotation_legend}{HeatmapList}(object, legend_list = list(), ...)
}
\arguments{

  \item{object}{A \code{\link{HeatmapList-class}} object.}
  \item{legend_list}{A list of self-defined legends, should be wrapped into \code{\link[grid:grid.grob]{grob}} objects. It is normally constructed by \code{\link{Legend}}.}
  \item{...}{Other arguments.}

}
\details{
We call the "annotation legends" as the secondary legends.
For horizontal heamtap list, the legends are those from all top/bottom annotations, and for vertical heatmap list,
the legends are those from all left/right annotations.

A viewport is created which contains annotation legends.

This function is only for internal use.
}
\value{
This function returns no value.
}
\author{
Zuguang Gu <z.gu@dkfz.de>
}
\examples{
# There is no example
NULL

}

max_stars_repo_name: zhongmicai/complexHeatmap