text:
#' Read MODFLOW .hds File
#'
#' This function reads in a binary hds file and creates a data frame
#' composed of the following columns:
#' \describe{
#' \item{TRANS}{Transport Time Step}
#' \item{STP}{Flow Time Step}
#' \item{PERLEN}{Stress Period Length}
#' \item{TIME}{Elapsed Time}
#' \item{LAY}{Model Layer}
#' \item{ROW}{Model Row}
#' \item{COL}{Model Column}
#' \item{GWE}{Groundwater Elevation}
#' }
#' @param rootname This is the rootname of the flow model (i.e., "F95")
#' @param NLAY This is the number of layers assigned to the model and needs to be
#'             assigned
#' @param NPER This is the number of stress periods that are printed to the .hds file. This number
#' can be obtained from the .dis file
#' @export
#' @examples
#' readhds("F95", NLAY = 8, NTTS = 16)

readhds <- function(rootname = NA, NLAY, NPER){
    if(is.na(rootname)){
            rootname <- MFtools::getroot()
    }
    # NLAY IS THE NUMBER OF LAYERS IN THE MODEL
    # NPER IS THE NUMBER OF TRANSPORT TIME STEPS IN THE MODEL
    hdsname <- paste(rootname, ".hds", sep = "")
    to.read <- file(hdsname, "rb")
    TRANS <- c()
    STP <- c()
    PERLEN <- c()
    TIME <- c()
    TEXT <- c()
    LAY <- c()
    GWE <- c()
    dat <- c()
    readblock <- function(){
          TRANS <- readBin(to.read, integer(), n = 1)
          STP <- readBin(to.read, integer(), n = 1)
          PERLEN <- readBin(to.read, double(), size = 4, n = 1)
          TIME <- readBin(to.read, double(), size = 4, n = 1)
          TEXT <- readChar(to.read, 16)
          NC <- readBin(to.read, integer(), n = 1)
          NR <- readBin(to.read, integer(), n = 1)
          LAY <- readBin(to.read, integer(), n = 1)
          GWE = readBin(to.read, double(), size = 4, n = NR * NC, endian = "little")
    out <- list(TRANS, STP, PERLEN, TIME, TEXT, NC, NR, LAY, GWE)
    return(out)
    }
    for(Q in 1:NPER){
        for(K in 1:NLAY){
            dat[[length(dat) + 1]] <- readblock()
        }
    }
    close(to.read)
    TRANS  <- sapply(dat, "[[", 1)
    STP    <- sapply(dat, "[[", 2)
    PERLEN <- sapply(dat, "[[", 3)
    TIME   <- sapply(dat, "[[", 4)
    LAY    <- sapply(dat, "[[", 8)
    NC     <- dat[[1]][6] %>% as.integer
    NR     <- dat[[1]][7] %>% as.integer
    GWE    <- sapply(dat, "[[", 9)
    HDS    <- tibble::data_frame(
                      TRANS = rep(TRANS, each = (NC * NR)) %>% as.integer(),
                      STP = rep(STP, each = (NC * NR)) %>% as.integer(),
                      PERLEN = rep(PERLEN, each = (NC * NR)) %>% as.numeric(),
                      TIME = rep(TIME, each = (NC * NR)) %>% as.numeric(),
                      LAY = rep(LAY, each = (NC * NR)) %>% as.integer(),
                      ROW = rep(rep(rep(1:NR, each = NC), NLAY), NPER) %>% as.integer(),
                      COL = rep(rep(rep(seq(1, NC, 1), NR), NLAY), NPER) %>% as.integer(),
                      GWE = GWE %>% as.numeric()
           )
    rm(TRANS)
    rm(STP)
    rm(PERLEN)
    rm(TIME)
    rm(TEXT)
    rm(NC)
    rm(NR)
    rm(LAY)
    rm(GWE)
    gc()
    return(HDS)
    }

max_stars_repo_name: dpphat/MFtools

max_issues_repo_name: dpphat/MFtools

max_forks_repo_name: dpphat/MFtools


text:
#' GMO - General Multinomial Overdisperion Model
#'
#' Creates an General Multinomial Overdispersion Model for the given data matrix.
#' The data passed in should be in the form of multinomial clusters.
#'
#' The `corstr` option specifies the intra-cluster correlation structure:
#'    - `independence` is used when there is no intra-cluster correlation (i.e.: rho = 0) (default)
#'    - `exchangebale` is used when there is intra-cluster correlation
#'
#' The `link` option specifies the format of the output:
#'    - `glogit` option gives the estimates and variance matrix in terms of beta (real-valued) (default)
#'    - `identity` option gives the estimates and variance matrix in terms of pi's in [0, 1]
#'
#' @param dat matrix of multinomial clusters
#' @param corstr one of: independence, exchangeable
#' @param link one of: glogit, identity
#'
#' @examples
#' dat <- matrix(c(6, 4, 4, 3, 1, 2, 2, 2, 1, 3, 15, 6, 5, 8, 4, 2, 2, 1, 13, 10, 8, 1, 3, 2, 10, 22, 11, 31, 10, 14), nrow=10, ncol=3, byrow=TRUE)
#' gmo(dat, link='identity')
#'
#' @return a GMO object
#' @export
gmo <- function (dat, corstr = "independence", link = "glogit") {
  if (!is.matrix(dat))
    stop("dat must be a matrix")
  if (!(corstr %in% c("independence", "exchangeable")))
    stop("corstr must be one of: independence, exchangeable")
  if (!(link %in% c("glogit", "identity")))
    stop("link must be one of: glogit, identity")

  RHO <- (corstr == "exchangeable")
  ndim <- dim(dat)[[2]]
  num.clus <- dim(dat)[[1]]
  out <- list(ndim = ndim, num.clus = num.clus)

  estimates <- iterative_function(dat = dat, RHO = RHO)
  out$num.iter <- estimates$num.iter
  out$error <- estimates$error ## TEMP: FOR TESTING
  if (RHO) {
    out$rho <- estimates$rho
  }

  if (link == "glogit") {
    # return betas as estimates
    out$betas <- estimates$betas
    out$var.mat <- get_var_glogit(estimates$betas, estimates$rho, dat)

    class(out) <- "gmo.glogit"
  } else {
    # return p.vec as estimates
    p.vec <- exp(estimates$betas) / (1 + sum(exp(estimates$betas)))

    out$p.vec <- c(p.vec, 1 - sum(p.vec))
    out$var.mat <- get_var_ident(p.vec, estimates$rho, dat)

    class(out) <- "gmo.ident"
  }
  out
}

#' Print GMO glogit
#'
#' Prints GMO glogit object.
#'
#' @param x gom glogit object
#' @param ... additional params to \code{print(x, ...)}
#'
#' @return None
#' @export
print.gmo.glogit <- function(x, ...) {
  cat('Number of Iterations: ', x$num.iter, '\n')

  cat('\nEstimates: \n')
  names(x$betas) <- c(sprintf("Beta%d", seq(1:length(x$betas))))
  print(round(x$betas, 6), ...)

  if (exists("rho", where = x)) {
    names(x$rho) <- "Rho"
    print(round(x$rho, 6), ...)
  }

  cat('\nVar-Covar Matrix: \n')
  print(round(x$var.mat, 6), ...)
}

#' Print GMO identity
#'
#' Prints GMO identity object.
#'
#' @param x gom identity object
#' @param ... additional params to \code{print(x, ...)}
#'
#' @return None
#' @export
print.gmo.ident <- function(x, ...) {
  cat('Number of Iterations: ', x$num.iter, '\n')

  cat('\nEstimates: \n')
  names(x$p.vec) <- c(sprintf("P%d", seq(1:length(x$p.vec))))
  print(round(x$p.vec, 6), ...)

  if (exists("rho", where = x)) {
    names(x$rho) <- "Rho"
    print(round(x$rho, 6), ...)
  }

  cat('\nVar-Covar Matrix: \n')
  print(round(x$var.mat, 6), ...)
}

ext: r


text:
args<-commandArgs(T)
library("stringr")

mydata=read.table(args[1],header=F,sep="\t")
stat=read.table(args[2],header=F,sep="\t")

base = as.numeric(args[3])
h = as.numeric(args[4])
p = c(0,0,0,0)
j <- 1;
cur <- 0

rn <- nrow(stat)
names_len <- max(nchar(as.character(stat[,2])))
bar <- str_c(rep("-",as.integer(names_len / 2) * 1.3 ),collapse='')

colnames(mydata)<-c("x","y")
coln<-c(NA,"5'UTR length (log2)")
pdf(paste(args[1],".pdf",sep=""),width=12,height=10)
par(mar=c(20,8,2,2))
boxplot(y~x,data=mydata,ylab=coln[2],xlab=coln[1],xlim=c(1,rn),names=NA)
#text(0.5,base,labels=bar,srt=60,xpd=TRUE,adj=1,cex=1.5)

for(i in 1:rn){
	if(cur < stat[i,3]){
		text(i-0.5,base,labels=bar,srt=60,xpd=TRUE,adj=1,cex=1.5)
		p[j] <- i - 0.5
		j <- j + 1
	}
	text(i,base,labels=paste(stat[i,2]," (",stat[i,3],")"),srt=60,xpd=TRUE,adj=1)
	cur <- stat[i,3]
}
text(rn+0.5,base,labels=bar,srt=60,xpd=TRUE,adj=1,cex=1.5)
p[j] <- rn + 0.5

text(p[1]+(p[2]-p[1])/2-2,h,labels="biological_process",xpd=TRUE,adj=1)
text(p[2]+(p[3]-p[2])/2-2,h,labels="cellular_component",xpd=TRUE,adj=1)
text(p[3]+(p[4]-p[3])/2-2,h,labels="molecular_function",xpd=TRUE,adj=1)

dev.off()


text:
source("team_helper_functions.r")

data <- read.csv("output/hclust.csv", na.strings=c(""," ","NA"), header=TRUE) # Loading data
table(data$cluster, data$Is_Resigning)
table(data$cluster, data$JobLevel)
table(data$cluster, data$Gender)


text:
# set up locations and reused functions

#output.path
outfolder <- "out_dada"

# function to normalize subtables and re-merge
norm_esv <- function(st) {
    pt <- prop.table(st,margin=1)
    pt[is.nan(pt)] <- 0
    pt <- pt[!is.na(rownames(pt)),]
    pt
}

# function to add duplicates back in
# finds zero distances from distance matrix df.dist, adds de-duped rows from df.xy back into df.xy
addDups <- function(df.dist, df.xy) {
    #stackoverflow solution to insert row: https://stackoverflow.com/questions/11561856/add-new-row-to-dataframe-at-specific-row-index-not-appended/11562428#11562428
    insertRow2 <- function(existingDF, newrow, r) {
      existingDF <- rbind(existingDF,newrow)
      existingDF <- existingDF[order(c(1:(nrow(existingDF)-1),r-0.5)),]
      row.names(existingDF) <- 1:nrow(existingDF)
      return(existingDF)
    }

    mat <- as.matrix(df.dist)
    zeros <- which(mat==0,arr.ind=T)
    zeros <- zeros[zeros[,1]!=zeros[,2],]
    zg <- graph.data.frame(zeros,directed=FALSE)
    zc <- data.frame(components(zg)$membership)
    zc[,2] <- as.numeric(rownames(zc))
    zc <- zc[order(zc[,1],zc[,2]),]
    rownames(zc) <- NULL
    colnames(zc) <- c("group","row")

    dups <- NULL
    for(k in unique(zc$group)) {
        duped <- sort(zc$row[which(zc$group==k)])
        kept <- duped[1]
        duped <- duped[-1]
        dups <- rbind(dups, cbind(rep(kept,length(duped)),duped) )
    }
    dups <- as.matrix(dups[order(dups[,2]),])

    out.xy <- df.xy
    for(j in 1:nrow(dups)) {
        out.xy <- insertRow2(out.xy, out.xy[dups[j,1],], unname(dups[j,2]))
    }

    out.xy
}


text:
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#
#    This file is part of ICTP RegCM.
#
#    ICTP RegCM is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    ICTP RegCM is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY# without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with ICTP RegCM.  If not, see <http://www.gnu.org/licenses/>.
#
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

library("ncdf")
library("PCICt")
library("fields")

surface = open.ncdf("/scratch/ggiulian/test/output/verysmall_SRF.1999050100.nc")
lonmat  = get.var.ncdf(nc=surface,varid="xlon")
latmat  = get.var.ncdf(nc=surface,varid="xlat")

timearr   = get.var.ncdf(nc=surface,varid="time")
timeunits = att.get.ncdf(nc=surface,varid="time",attname="units")
calendar  = att.get.ncdf(nc=surface,varid="time",attname="calendar")

tpr       = get.var.ncdf(nc=surface,varid="pr")
tprunits  = att.get.ncdf(nc=surface,varid="pr",attname="units")

seconds_per_hour <- 3600.0
origin <- as.PCICt(substring(timeunits$value,13,31),cal=calendar$value)
times <- origin + timearr*seconds_per_hour

inds = (1:dim(timearr))
targettime = as.PCICt("1999-05-06 12:00:00",cal=calendar$value)
tind = inds[times == targettime]

image.plot(lonmat, latmat,tpr[,,tind])

close.ncdf(surface)

max_stars_repo_name: taobrienlbl/RegCM

max_issues_repo_name: taobrienlbl/RegCM

max_forks_repo_name: taobrienlbl/RegCM


text:
# Lesson 4: Filtering your data
# Objectives:
  #1.) Remove chr prefixes from chromosome name
  #2.) Order the chromosomes correctly
  #3.) Pick a subset of the types and rename them.
library(ggplot2)

filename <- "lesson_3/Encode_HMM_data.txt"

my_data <- read.csv(filename, sep="\t", header=FALSE)

head(my_data)

names(my_data)[1:4] <- c("chrom", "start", "stop", "type")

ggplot(my_data, aes(x=chrom, fill=type)) + geom_bar()

dim(my_data)

summary(my_data$chrom)
summary(my_data$type)

# levels = possibilities/categories
# remove "chr" prefix
# gsub() allows you to replace a certain substring with another substring for any given levels.
# factor() allows you to create a data object that has those new categories
my_data$chrom <- factor(gsub("chr", "", my_data$chrom))
summary(factor(gsub("chr", "", my_data$chrom)))
summary(my_data$chrom)
ggplot(my_data, aes(x=chrom, fill=type)) + geom_bar()

# seq() works like range but is inclusive on both ends
seq(1,22)

c(seq(1,22), "X", "Y")

my_data$chrom <- factor(my_data$chrom, levels=c(seq(1,22), "X", "Y"))
factor(my_data$chrom)
summary(my_data$chrom)

# Run the plot to make sure that the prefixes have been renamed and reordered
ggplot(my_data, aes(x=chrom, fill=type)) + geom_bar()

# Filter the data to just only get plots of a few, interesting types
# The %in% returns a bool value for each entry in the data frame

my_data <- my_data[my_data$type %in% c("1_Active_Promoter", "4_Strong_Enhancer", "8_Insulator"), ]
summary(my_data$type)

# Rename the types
library(plyr) # We need the revalue() function from plyr

my_data$type <- revalue(my_data$type, c("1_Active_Promoter"="Promoter",
                                        "4_Strong_Enhancer"="Enhancer",
                                        "8_Insulator"="Insulator"))
# check the updated plot
png("lesson_4/updatedplot.png")
ggplot(my_data, aes(x=chrom, fill=type)) + geom_bar()
dev.off()


text:
require(httr)

data = '{"title":"china1"}'

res <- httr::POST(url = 'http://example.com/post', body = data)

max_stars_repo_name: verhovsky/curlconverter

max_issues_repo_name: verhovsky/curlconverter

max_forks_repo_name: verhovsky/curlconverter


text:
## **Importando Bibliotecas**

library(readr)
library(tidyverse)
library(rpart)
library(fastDummies)
library(caret)
library(rattle)
```
## **Importando dados**

dados <- read_csv("acessov2.csv")
dados

## **Tratando os dados**

table(dados$busca)

### Metodo de separaçao em Dummies, para poder trabalhar com a coluna busca.

dados <- dummy_cols(dados, remove_selected_columns = T)
dados

set.seed(17)
idx <- createDataPartition(dados$comprou, p = 0.80, list = F)
dados.train <- dados[idx,]
dados.train <- dados.train[c(1,2,4,5,6,3)]
dados.test <- dados[-idx,]
dados.test.r <- dados.test[3]
dados.test <- dados.test[-3]

### **Checando proporçoes**

# Proporçao dos dados inicias
prop.table(table(dados$comprou))
# Proporçao do treino
prop.table(table(dados.train$comprou))
# Proporçao do teste
prop.table(table(dados.test.r$comprou))

## **Modelo Decision Tree Rpart**
### Funçao para categorizar em 0 e 1.

# 0.375 0.611111111111111 0.704697986577181 0.795275590551181 1
bins <- function(df)map_dbl(df, function (x){
  if (x <= 0.704697986577181){
    return(0)
  }
  else{
    return(1)
  }})

modelo1 <- rpart(comprou~., dados.train, method = "anova")
previsao_rpart <- predict(modelo1, dados.test)
dados.test.r["previsao_rpart"] <- bins(previsao_rpart)

## **Grafico do modelo rpart**

fancyRpartPlot(modelo1, main = "Arvore de Decisoes", sub = "")

## **Matriz de confusao**

confusionMatrix(table(dados.test.r$comprou,
                      dados.test.r$previsao_rpart),
                positive = "1")

max_stars_repo_path: rpart e carat.r

max_stars_repo_name: Oaccioly/R-Machine-Learning

max_issues_repo_path: rpart e carat.r

max_issues_repo_name: Oaccioly/R-Machine-Learning

max_forks_repo_path: rpart e carat.r

max_forks_repo_name: Oaccioly/R-Machine-Learning


text:
#' as.banded
#'
#' Convert a regular R matrix into a banded matrix.
#'
#' @details
#' The values \code{kl} and \code{ku} describe the bandwidth of the band
#' matrix.  If they are different, they will produce a general banded
#' matrix.  If they are both 0, they will produce a diagonal matris.
#'
#' @param x
#' A matrix containing numeric, integer, or logical data.
#' @param kl,ku
#' The size of the lower and upper bands.  At least one value \code{kl} or
#' \code{ku} must be supplied; however, if one is missing, the other will
#' be assumed to be the same value. Each supplied value should be a
#' non-negative integer. See the details section for more information
#' about these values.
#'
#' @return
#' A banded matrix.
#'
#' @examples
#' x <- matrix(1:25, nrow=5)
#'
#' # A diagonal matrix
#' as.banded(x, 0)
#'
#' # A general banded matrix
#' as.banded(x, 2, 1)
#'
#' @references \url{http://www.netlib.org/lapack/lug/node124.html}
#' @seealso \code{\link{kdim}} and \code{\link{classes}}
#' @keywords Casters
#' @export
as.banded <- function(x, kl, ku)
{
  check.is.matrix(x)
  if (!is.numeric(x) && !is.logical(x))
    stop("argument 'x' must contain numeric, integer, or logical data")

  if (missing(kl) && missing(ku))
    stop("must supply at least one of 'kl' or 'ku'")
  else if (missing(kl))
    kl <- ku
  else if (missing(ku))
    ku <- kl

  dim <- as.integer(dim(x))

  if (ku > dim[1L] || kl > dim[1L])
    stop("'ku' and 'kl' must each be smaller than nrow(x)")
  if (max(ku, kl) > min(dim))
    warning("banded storage schemes are only practical if kl, ku << min(dim(x))")

  ku <- as.integer(ku)
  kl <- as.integer(kl)

  out <- .Call(R_tobanded, x, kl, ku)

  if (kl == 0L && ku == 0L)
    ret <- diagmat(Data=out, dim=dim)
  else
    ret <- genbandmat(Data=out, dim=dim, kl=kl, ku=ku)

  ret
}

max_stars_repo_name: wrathematics/band

max_issues_repo_name: wrathematics/band

max_forks_repo_name: wrathematics/band


text:
Using the dplyr package, use the 6 different operations to analyze/transform the data -
GroupBy, Summarize, Mutate, Filter, Select, and Arrange – Remember this isn’t just
modifying data, you are learning about your data also – so play around and start to understand your dataset in more detail

data_house <- read_xlsx("week-7-housing.xlsx")

#GroupBy:-

> data_house %>%
+ group_by(square_feet_total_living) %>%
+ summarise(n = n())
# A tibble: 654 x 2
   square_feet_total_living     n
                      <dbl> <int>
 1                      240     1
 2                      310     1
 3                      340     1
 4                      410     1
 5                      430     1
 6                      480     2
 7                      530     1
 8                      540     1
 9                      550     2
10                      570     1
# ... with 644 more rows
>

#Summarize:-

summarise(data_house, square_feet_total_livings = mean(data_house$square_feet_total_living,na.rm = TRUE ))
# A tibble: 1 x 1
  square_feet_total_livings
                      <dbl>
1                     2540.summarise(data_house, square_feet_total_livings = mean(data_house$square_feet_total_living,na.rm = TRUE ))
# A tibble: 1 x 1
  square_feet_total_livings
                      <dbl>
1                     2540.

#Mutate:-

data_house1 <- mutate(data_house, years_old = as.numeric (format(as.Date(data_house$`Sale Date`), "%Y")) - year_built)

#Filter:-

filter(data_house1 ,years_old > 10)

#Select:-

select(data_house1,`Sale Date`, years_old)

#Arrange:-

> arrange(data_house1,desc(years_old), `Sale Date`)

Using the purrr package – perform 2 functions on your dataset.  You could use zip_n, keep, discard, compact, etc.

keep:-

data_house1$years_old %>% keep(~ .x  > 10)

discard :-

data_house1$years_old %>% discard(~ .x  > 10)

compact :-

data_house1$years_old %>% compact(~ .x)

Use the cbind and rbind function on your dataset

#cbind:-

> cbind_hdata <- cbind(data_house,data_house1)
> colnames(cbind_hdata)
 [1] "Sale Date"                "Sale Price"
 [3] "sale_reason"              "sale_instrument"
 [5] "sale_warning"             "sitetype"
 [7] "addr_full"                "zip5"
 [9] "ctyname"                  "postalctyn"
[11] "lon"                      "lat"
[13] "building_grade"           "square_feet_total_living"
[15] "bedrooms"                 "bath_full_count"
[17] "bath_half_count"          "bath_3qtr_count"
[19] "year_built"               "year_renovated"
[21] "current_zoning"           "sq_ft_lot"
[23] "prop_type"                "present_use"
[25] "Sale Date"                "Sale Price"
[27] "sale_reason"              "sale_instrument"
[29] "sale_warning"             "sitetype"
[31] "addr_full"                "zip5"
[33] "ctyname"                  "postalctyn"
[35] "lon"                      "lat"
[37] "building_grade"           "square_feet_total_living"
[39] "bedrooms"                 "bath_full_count"
[41] "bath_half_count"          "bath_3qtr_count"
[43] "year_built"               "year_renovated"
[45] "current_zoning"           "sq_ft_lot"
[47] "prop_type"                "present_use"
[49] "years_old"
>
rbind:-

rbind(new_bindhouse,data_house1)

Split a string, then concatenate the results back together

Split String:-

df=data.frame(do.call("rbind", strsplit(as.character(rs), " ", fixed = TRUE)))

> df
    X1 X2    X3 X4     X5      X6
1 This is First  R String Example

concatenate (Paste):-

totstr <- paste(df$X1,df$X2, df$X3, df$X4, df$X5, df$X6, sep=" ")

> totstr
[1] "This is First R String Example"


text:
a = 10
class(a)
a = as.integer(a)
class(a)
b = "red"
b = 'red'
b
d = TRUE
class(d)
d = T
class(d)
print(a)
print(d)
d = 0
class(d)
d = as.logical(d)
d
a = c(1,2,3,4,5)
class(a)
b = c('vijay', 'kumar')
class(b)
print(a)
print(b)
a
b
c
d
a[0]
print(a[0])
a[1]
a[1:3]
a[3:5]
a[-3]
a[6]
length(a)
a[2,4]
a[-1,-3,5]
a[-1,-3,-5]
a[2::4]
a[c(2,4)]
a[2] = 10
a
a[6] = 6
a
a = 12
a
a = c(1,2,3,4,5,6)
a
a = c(a,13,12,15)
a
a = c(a[1:2],11,12,a[3:6])
a
[print(a[1])]
print(a[1])
print(a[1:2])
prin(a[1:2],[-7,-8])
print(a[1:2],[-7:-8])
print(a[-7:-8])
print(c(a[1:2],a[8:9]))
a
a[9] = 15
a
sort(a)
a
d = sort(a)
d
a = sort(a,decreasing = TRUE)
a
a = 1:100
a
a = seq(1,100,by = 2)
a
a = seq(1,100, by = 3)
a
a = seq(2,100, by = 2)
a
a = c(1,3, seq(5:11, by = 2))
a = c(1,3, seq(5,11, by = 2))
a
a = rep(0,3)
a
a = c()
a
a = c(a,1)
a
a = c(1)
a
a = c(a,1)
a = c(2)
a = c(a,1)
a
a = c(a, 2:3)
a
a = c("red")
a = c("red", "white", "yellow")
a
class(a)
a[1]
a[4] = c("violet")
a
a = 1:10
b = 11:20
c = a + b
a* b
b = 1:5
a * b
a[2] * 4
a = 10
b = 2
a / b
b = 3
a / b
b = 2
a % b
a %% b
a // b
a %/% b
c = 3
a %/% b
a %/% c
a / c
a = red
a = "red"
b = "white"
c(a+b)
paste(a,b)
paste(c,d)
a = factor(c("M", "F", "M", "M", "M", "F"))
a
class(a)
b = factor(c("vijay", "vijay"))
b
a
table(a)
table(b)
a = c("E", "W", "N", "S", "E", "W")
a
class(a)
a = as.factor(a)
a
table(a)
class(a)
savehistory("D:/Study/3rd Sem/IDS/Programs/r-objects.r")
a = matrix(1:10, nrow = 5)
a
a = matrix(1:10, ncol = 5)
a
a = matrix(0,3,3)
a
a[1,1]
a
a = matrix(1:10, 5)
a
a[1, ]
a[3, -1]
a[3, 1]
a[3, -2]
a
a*2
a * b
a*b
b = c(1,2,NA,5)
b
class(b)
a = 1:5
a
mean(a)
median(a)
mode(a)
b = NA.omit(b)
b = na.omit(b)
b
mean(b)
a * b
sqrt(16)
sqrt(16)
sqrt(729)
b = 1947.0
b = as.character(b)
class(b)
d = c(1:6)
class(d)
e = c(1, 'a', 2, 'b')
class(e)
f = c(a:z)
f = c(1:26)
f = as.character(f)
class(f)
length(f)
f[1] = "a"
f
names = c('Ayush', 'Vijay', 'Sudhanva', 'Karthik', 'Shreedhar')
names
class(names)
length(names)
names[1:2]
names[2:3]
sort(names)
sort(names, decreasing = TRUE)
rep(1,11)
a = sort(rep('a',2))
a
a = sort(rep('a',2), (1,1))
a = sort(rep('a',2)(1,1))
a = sort(rep('a',2,1:8))
a = c(1,2,NA,4)
a
a = na.omit(a)
a
rep('a',2,seq(1:5),seq(7,11))
a = c(rep("a", 2),1:5,seq(7,11,by = 2))
a
a = Letters[1:26]
a
a
z
z = LETTERS[1:26]
LETTERS[1:26]
a = sort(rep('a',2)(1,1))
LETTERS[1:26]
a = 1:10000
a
b = sample(a,100)
b
y = sample(y,100)
y = 1:10000
y = sample(y,100)
y
b = sample(a,100, replace = TRUE)
b
c = sample(a,7)
c
mean(b)
mean(c)
summary(b)
summary(c)
savehistory("D:/Study/3rd Sem/IDS/Programs/r-objects.r")

max_stars_repo_name: nsudhanva/mca-code

max_issues_repo_name: nsudhanva/mca-code

max_forks_repo_name: nsudhanva/mca-code


text:
#' Metabolite set enrichment analysis (MSEA) (using a hypergeometric test) using pathway knowledge curated by Metabolon
#'
#' A function that returns the pathway enrichment score for all perturbed metabolites in a patient's full metabolomic profile.
#' @param met.profile - A character vector of a patient's metabolomic profile, including KEGG IDs
#' and the associated z-score or p-value describing the level of the metabolite compared to controls.
#' @param threshold - A cutoff to select metabolites with a zscore > threshold or < -1*threshold.
#' @param type - Either "p-value" or "z-score".
#' @param gene.profile - Default set to NULL, meaning the default enrichment analysis only considers metabolites.
#'                       However, if you have gene data, too, set this parameter to a character vector of the gene
#'                      names with found variants in the patient's record. Gene IDs must be converted to Entrez Identifiers.
#' @export stat.getORA_Metabolon
#' @examples
#' pathway.data = stat.getORA_Metabolon(met.profile, threhold=3, "z-score", NULL)
stat.getORA_Metabolon = function(met.profile, threshold=3, type="zscore", gene.profile=NULL) {
  met.profile = met.profile[which(!(is.na(met.profile)))]
  if (type=="zscore") {
    perturbed.mets = met.profile[which(abs(met.profile) > threshold)]
  } else {
    perturbed.mets = met.profile[which(met.profile < 0.05)]
  }
  nms.perturbed.mets = unique(as.character(unlist(sapply(names(perturbed.mets), function(i) strsplit(i, split=";")))))

  # The size of the population of total possible metabolites to draw from
  population = names(met.profile)
  print(sprintf("Total number of metabolites in profile. = %d.", length(population)))
  paths.hsa = list.dirs(path=system.file("extdata", package="CTDext"), full.names = FALSE)
  paths.hsa = paths.hsa[-which(paths.hsa %in% c("", "RData", "allPathways", "MSEA_Datasets", "Pathway_GMTs"))]
  row = 1
  pathway.data = data.frame(Pathway=character(), FDR=numeric(), Pvalue=numeric(), Hits=integer(), Size=integer(), stringsAsFactors = FALSE)
  for (pathway in 1:length(paths.hsa)) {
    load(system.file(sprintf("extdata/RData/%s.RData", paths.hsa[pathway]), package="CTDext"))

    pathway.compounds = tolower(V(ig)$label[which(V(ig)$shape=="circle")])
    pathCompIDs = unique(tolower(pathway.compounds[which(pathway.compounds %in% population)]))
    print(sprintf("Total number of metabolites in profile also in current Metabolon pathway. = %d.", length(pathCompIDs)))

    # q (sample successes), m (population successes), n (population failures), k (sample size)
    sampleSuccesses = length(which(nms.perturbed.mets %in% pathCompIDs))
    populationSuccesses = length(intersect(pathCompIDs, population))
    N = length(population)
    populationFailures=N-populationSuccesses
    numDraws=length(perturbed.mets)
    if (populationSuccesses>0) {
      pathway.data[row, "Pathway"] = paths.hsa[pathway]
      pathway.data[row, "Pvalue"] = phyper(q=sampleSuccesses-1, m=populationSuccesses, n=populationFailures, k=numDraws, lower.tail=FALSE)
      pathway.data[row, "Hits"] = sampleSuccesses
      pathway.data[row, "Size"] = populationSuccesses
      row = row + 1
    }
  }
  pathway.data[,"FDR"] = p.adjust(pathway.data[,"Pvalue"], method="fdr")
  # Sort by FDR, then Pvalue
  pathway.data = pathway.data[order(pathway.data[,"FDR"], pathway.data[,"Pvalue"]),]

  # Then, only return rows that have Pvalue < 0.25
  pathway.data = pathway.data[which(pathway.data[,"Pvalue"]<0.25),]

  return(pathway.data)
}

max_stars_repo_name: BRL-BCM/CTDext

max_issues_repo_name: BRL-BCM/CTDext

max_forks_repo_name: BRL-BCM/CTDext


text:
n <- 68
Fmax <- function(y, theta, log = FALSE){
  ans <- n*(log(y)-log(theta))
  if(y > theta) ans <- 0
  if(!log) ans <- exp(ans)
  return(ans)
}
poder <- function(theta){
  return(Fmax(2.9, theta) + 1 - Fmax(4, theta))
}
poder <- Vectorize(poder)

curve(poder, 1, 5, xlab = expression(theta),
      ylab = expression(pi(theta*"|"*delta)), lwd = 2)
curve(poder, 2.9, 4, lwd = 6, col = 2, add = TRUE)
abline(v = c(3, 4), lwd = 2, lty = 2)

curve(poder, 3, 4, xlab = expression(theta),
 ylab = expression(pi(theta*"|"*delta)), lwd = 2)


text:
#' Column Operations
#'
#' Column operations (currently sums/means) for shaq objects.
#'
#' @section Communication:
#' The operation consists of a local column sum operation, followed by an
#' \code{allreduce()} call, quadratic on the number of columns.
#'
#' @param x
#' A shaq.
#' @param na.rm
#' Should \code{NA}'s be removed?
#' @param dims
#' Ignored.
#'
#' @return
#' A regular vector.
#'
#' @examples
#' \dontrun{
#' library(kazaam)
#' x = ranshaq(runif, 10, 3)
#' cs = colSums(x)
#' comm.print(cs)
#'
#' finalize()
#' }
#'
#' @name col_ops
#' @rdname col_ops
NULL

colsum.shaq = function(x, na.rm = FALSE, dims = 1L)
{
  check.is.flag(na.rm)

  cs = colSums(DATA(x), na.rm=na.rm)
  MPI_Allreduce(cs)
}

colmean.shaq = function(x, na.rm = FALSE, dims = 1L)
{
  check.is.flag(na.rm)

  cs = colSums(DATA(x), na.rm=na.rm)
  MPI_Allreduce(cs) / nrow(x)
}

#' @rdname col_ops
#' @export
setMethod("colSums", signature(x="shaq"), colsum.shaq)

#' @rdname col_ops
#' @export
setMethod("colMeans", signature(x="shaq"), colmean.shaq)

max_stars_repo_name: snoweye/kazaam

max_issues_repo_name: snoweye/kazaam

max_forks_repo_name: snoweye/kazaam


text:
ExtractMostRelevantCoefficient <- function(singleModel) {
  # This function extracts for each country the most relevant coefficient in
  # terms of absolute value, and save it into a data table.
  #
  # ARGS
  # - singleModel [model]  the model from which one wants to extract the
  #                        coefficients
  #
  # RETURN
  # - dt  [data table]  a data table with the name of the coefficient, its value,
  #                     and the associated country

  dt <- data.table()

  # For each country
  for(country in 1:length(singleModel)) {
    precisModel <- precis(singleModel[[country]])

    # Select intervals not containing 0
    valueWithZero <- precisModel[[3]][-1] <= 0 & precisModel[[4]][-1] >= 0
    selectValue <- precisModel[[1]][-1][valueWithZero == FALSE]

    if (length(selectValue) > 0) {
      # Find the maximum among the coefficients (no intercept)
      maxCoef <- max(abs(selectValue))

      # Find the row of the associated coefficent
      mainCoeff <- precisModel[abs(precisModel$mean) == maxCoef, ]

      tmp <- data.table(name = row.names(mainCoeff),
                        mainCoef = mainCoeff[[1]],
                        country = names(singleModel)[[country]])
    } else {
      tmp <- data.table(name = "NA",
                        mainCoef = NA,
                        country = names(singleModel)[[country]])
    }

    dt <- rbind(dt, tmp)
  }

  return(dt)
}

max_stars_repo_name: formozov/Global-Preferences-Survey

max_issues_repo_name: formozov/Global-Preferences-Survey

max_forks_repo_name: formozov/Global-Preferences-Survey


text:
Begin fired first
Begin fired second
Begin fired third
Begin fired fourth
Begin fired fifth
End1 fired after exit
End2 fired after exit
End3 fired after exit
End4 fired after exit

max_stars_repo_name: tjfontaine/dtrace-utils

max_issues_repo_name: tjfontaine/dtrace-utils

max_forks_repo_name: tjfontaine/dtrace-utils


text:
#
#           script meteodusol_EM50.r
#           version 20220323
#           Philippe Choler
#

# Ce script gÃ©nÃ¨re deux fichiers csv contenant (1) les mÃ©tadonnÃ©es et (2) les donnÃ©es de teneur en eau du sol pour leur bancarisation dans OSUG-DC dans le cadre du projet climaplant.

# Les mesures sont rÃ©alisÃ©es avec les enresgistreurs EM50

# PrÃ©-requis
# 1. un fichier contenant les attributs des sites (mysite_EM50.csv)
#  SourceID = nom du site attribuÃ© par le fournisseur de donnÃ©es
#  X_WGS84  = longitude en degrÃ©s dÃ©cimaux
#  Y_WGS84  = latitude en degrÃ©s dÃ©cimaux
#  loggerID   = numÃ©ro de sÃ©rie de l'enregistreur

# 2. un rÃ©pertoire (DIR.CSV) contenant les fichiers export des enregistreurs TOMST-TMS
# le nom des fichiers est du type "loggerID_xxxxx.csv", le loggerID ne doit pas contenir le caractÃ¨re "_"
# colonne 1 : numÃ©ro de ligne
# colonne 2 : DATA HEURE GMT+00:00" (par exemple 2017-05-18 11:00:00)
# colonne 4 Ã  6 : tempÃ©rature en Â°C. (par exemple 10,2 et non pas 10.2)
# colonne 7 : valeurs brutes (raw values) de teneur en eau
# ...
# le sÃ©parateur de colonnes est ;

# IMPORTANT : ParamÃ¨tres Ã  ajuster par chaque utilisateur

# You need first to set working directory (WD) to source file location
setwd("WD")

# a. DIR.CSV le chemin du rÃ©pertoire contenant les fichiers csv des enregistreurs Hobo
DIR.CSV    <- "./CSV_EM50/"
# b. mysite = le chemin d'accÃ¨s au fichier mysite
mysite     <- "./mysite_EM50.csv"

# c. dataset  = le nom du lot de donnÃ©es
dataset    <- "pne"

# d. DIR.OUT = le chemin du rÃ©pertoire contenant les fichiers de sortie
DIR.EXPORT  <- "./EXPORT_EM50/"

# e. initialisation du lot pour incrÃ©mentation
INIT <- 1

# vÃ©rification des fichiers d'entrÃ©e
FILES.CSV  <- list.files(DIR.CSV,pattern=".csv",full=T)
FILES.CSVs <- list.files(DIR.CSV,pattern=".csv",full=F)
NAMES      <- unlist(lapply(strsplit(FILES.CSVs,"_"),function(x) x[1]))

(SITE           <- read.table(mysite,sep=";",stringsAsFactors = F,header=T))
colnames(SITE) <- c("SourceID","X_WGS84","Y_WGS84") # renomme les colonnes

# prepare data and metadata
md         <- NULL
data       <- NULL
DOY        <- gsub("-","",Sys.Date())   # date

set.seed(2067)

# import files and produce data and metadata
for (i in 1:length(FILES.CSV)){
  print(i)
  TS        <- read.csv(FILES.CSV[i],skip=2,sep=",",header=T)
  TIME      <- as.POSIXct(TS[,1],format = '%d/%m/%Y %H:%M:%S',tz="UTC")
  DAY       <- as.Date(TIME)
  SourceID  <- NAMES[i]
  META      <- SITE[match(SourceID,SITE$SourceID),]
  NL        <- nrow(TS)
  ID        <- round(runif(1,10^9,10^10))

  # metadata
  TMPmd <- data.frame(
    INCREMENT     = paste0(dataset,INIT+i-1),
    ID            = ID,
    X             = round(META$"X_WGS84",5),
    Y             = round(META$"Y_WGS84",5),
    sensor_type   = "DECAGON",
    date_start    = DAY[1],
    date_end      = DAY[NL],
    subdataset    = dataset,
    sourceID      = META$SourceID
  )
  md <- rbind(md,TMPmd)

  # data
  for (j in 2:ncol(TS)) TS[,j] <- as.numeric(gsub(",",".",TS[,j]))
  TMPdata <- data.frame(
    INCREMENT     = paste0(dataset,INIT+i-1),
    ID            = ID,
    date_time     = format(TIME,"%Y-%m-%d %H:%M"),
    eau           = round(TS[,-1],2)
  )

  nam  <- expand.grid(c("VWC","TEMP"),1:4)
  colnames(TMPdata)[-(1:3)] <- paste0(nam[,1],nam[,2])
  data  <- rbind(data,TMPdata)

} # end of CSV files

# sauvegarde des fichiers
write.csv(md,paste0(DIR.EXPORT,DOY,"-",dataset,"-metadata-eau-soil.csv"),row.names=F,quote=F)
write.csv(data,paste0(DIR.EXPORT,DOY,"-",dataset,"-data-eau-soil.csv"),row.names=F,quote=F)

max_stars_repo_name: pcholer/meteodusol

max_issues_repo_name: pcholer/meteodusol

max_forks_repo_name: pcholer/meteodusol


text:
library(data.table)
library(dplyr)

indir = "~/git/R0t"

#
# Contact matrices from Zhang et al. Science https://science.sciencemag.org/content/368/6498/1481

pre_W = matrix(nrow = 14, ncol = 14, byrow = T,
               c(0.11, 0.35, 0.42, 0.55, 0.88, 0.35, 0.58, 0.12, 0.00, 0.00, 0.00, 0.00, 0.11, 0.33,
                 1.69, 2.31, 2.27, 5.84, 0.61, 0.58, 0.62, 0.55, 0.21, 0.23, 0.29, 0.17, 0.16, 0.18,
                 0.06, 3.76, 4.08, 8.39, 0.65, 0.75, 0.43, 0.31, 0.47, 0.60, 0.20, 0.31, 0.03, 0.16,
                 0.00, 1.35, 1.62, 4.41, 0.56, 0.43, 0.40, 0.24, 0.61, 0.47, 0.24, 0.28, 0.23, 0.14,
                 0.31, 0.38, 0.16, 0.82, 2.66, 1.65, 1.67, 1.37, 0.81, 0.98, 0.61, 0.67, 0.18, 0.41,
                 1.41, 0.70, 0.37, 0.99, 3.08, 1.67, 1.66, 1.43, 0.75, 0.97, 0.43, 0.48, 0.20, 0.45,
                 1.33, 1.29, 0.35, 1.11, 2.65, 1.82, 1.60, 1.26, 0.81, 1.15, 0.76, 0.62, 0.29, 0.48,
                 0.84, 0.99, 0.53, 1.36, 3.10, 1.92, 1.83, 1.70, 1.33, 1.45, 0.83, 0.83, 0.50, 0.86,
                 0.84, 0.60, 0.33, 0.83, 2.37, 1.20, 1.21, 1.23, 0.78, 1.30, 0.76, 0.58, 0.76, 1.60,
                 0.43, 0.88, 0.42, 1.39, 1.41, 0.76, 0.94, 0.91, 1.03, 1.18, 0.97, 0.79, 0.49, 0.88,
                 0.37, 0.58, 0.22, 0.75, 1.26, 0.94, 0.80, 0.60, 0.92, 1.20, 0.89, 0.66, 1.00, 2.03,
                 0.55, 0.69, 0.24, 1.07, 1.92, 0.98, 1.02, 1.34, 1.16, 0.62, 0.95, 0.84, 1.64, 4.10,
                 0.40, 0.24, 0.13, 0.39, 1.33, 0.73, 0.78, 0.61, 1.01, 1.76, 0.89, 0.76, 1.18, 2.76,
                 1.06, 0.68, 0.31, 0.86, 2.06, 1.50, 1.31, 0.69, 1.11, 1.17, 0.74, 1.22, 0.98, 1.91))

post_W = matrix(nrow = 14, ncol = 14, byrow = T,
                c(0.00, 0.11, 0.00, 0.00, 0.00, 0.00, 1.11, 0.44, 0.23, 0.00, 0.00, 0.00, 0.21, 0.11,
                  0.00, 0.00, 0.00, 0.00, 0.06, 0.05, 0.45, 0.87, 0.20, 0.10, 0.07, 0.09, 0.22, 0.09,
                  0.00, 0.00, 0.00, 0.04, 0.00, 0.00, 0.03, 0.78, 0.59, 0.32, 0.16, 0.03, 0.00, 0.07,
                  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.43, 0.47, 0.72, 0.15, 0.00, 0.05, 0.20,
                  0.04, 0.02, 0.01, 0.05, 0.13, 0.05, 0.06, 0.02, 0.19, 0.43, 0.42, 0.21, 0.11, 0.10,
                  0.24, 0.16, 0.03, 0.04, 0.10, 0.31, 0.25, 0.10, 0.10, 0.04, 0.36, 0.21, 0.10, 0.03,
                  0.29, 0.38, 0.16, 0.01, 0.00, 0.07, 0.39, 0.24, 0.10, 0.02, 0.06, 0.11, 0.28, 0.19,
                  0.14, 0.30, 0.29, 0.06, 0.00, 0.02, 0.22, 0.31, 0.28, 0.04, 0.03, 0.04, 0.17, 0.19,
                  0.05, 0.07, 0.32, 0.27, 0.08, 0.01, 0.02, 0.11, 0.40, 0.27, 0.07, 0.02, 0.05, 0.30,
                  0.00, 0.08, 0.03, 0.24, 0.27, 0.16, 0.02, 0.08, 0.17, 0.37, 0.28, 0.09, 0.02, 0.41,
                  0.04, 0.01, 0.03, 0.03, 0.24, 0.23, 0.09, 0.05, 0.03, 0.13, 0.32, 0.25, 0.05, 0.21,
                  0.26, 0.00, 0.05, 0.00, 0.11, 0.20, 0.16, 0.05, 0.20, 0.00, 0.15, 0.31, 0.10, 0.36,
                  0.06, 0.10, 0.01, 0.01, 0.05, 0.01, 0.08, 0.08, 0.02, 0.06, 0.08, 0.14, 0.53, 0.23,
                  0.00, 0.06, 0.06, 0.00, 0.00, 0.00, 0.12, 0.12, 0.06, 0.06, 0.06, 0.06, 0.23, 0.47))

pre_S = matrix(nrow = 14, ncol = 14, byrow = T,
               c(3.17, 0.24,  0.06, 0.03,  0.07, 0.82, 1.21, 0.68, 0.17, 0.11, 0.26, 0.53, 0.63, 0.62,
                 0.29, 13.03, 1.76, 0.04,  0.13, 0.36, 0.93, 1.41, 0.26, 0.16, 0.04, 0.30, 0.57, 0.86,
                 0.00, 1.14,  17.5, 2.51,  0.25, 0.29, 0.47, 1.04, 0.92, 0.49, 0.23, 0.13, 0.16, 0.69,
                 0.00, 0.05,  0.97, 17.64, 0.85, 1.03, 1.64, 2.57, 1.65, 1.59, 0.94, 0.71, 0.28, 0.92,
                 0.10, 0.00,  0.08, 0.07,  2.15, 2.67, 2.44, 3.04, 2.95, 3.17, 1.72, 1.27, 0.61, 0.82,
                 0.16, 0.07,  0.07, 0.21,  1.14, 2.76, 2.56, 3.29, 2.05, 2.41, 1.86, 1.74, 0.87, 1.27,
                 0.42, 0.41,  0.06, 0.32,  1.44, 2.65, 4.11, 4.53, 2.81, 3.38, 1.98, 1.86, 1.05, 1.17,
                 0.11, 0.29,  0.33, 0.18,  0.76, 1.57, 2.30, 3.07, 2.17, 2.32, 1.31, 1.33, 0.69, 0.63,
                 0.04, 0.08,  0.37, 0.45,  0.76, 1.55, 2.20, 3.23, 2.21, 2.82, 1.34, 0.90, 0.37, 1.51,
                 0.04, 0.09,  0.25, 0.27,  1.42, 1.95, 2.64, 2.81, 2.78, 2.91, 2.03, 1.18, 0.58, 1.53,
                 0.09, 0.06,  0.12, 0.24,  1.04, 1.77, 1.48, 2.70, 2.15, 2.57, 2.10, 1.68, 1.48, 2.56,
                 0.10, 0.09,  0.04, 0.07,  0.61, 1.24, 1.66, 1.96, 1.32, 1.75, 1.46, 2.37, 1.08, 2.79,
                 0.05, 0.03,  0.10, 0.05,  0.26, 0.72, 0.96, 1.10, 0.82, 0.93, 1.03, 1.57, 1.78, 2.46,
                 0.03, 0.12,  0.06, 0.06,  0.13, 0.27, 0.25, 0.56, 0.69, 0.74, 0.68, 0.88, 1.48, 5.56))

post_S = matrix(nrow = 14, ncol = 14, byrow = T,
                c(0.00, 0.00, 0.00, 0.00, 0.00, 0.29, 0.55, 1.01, 0.14, 0.00, 0.14, 0.00, 0.00, 0.00,
                  0.04, 0.00, 0.00, 0.00, 0.00, 0.26, 0.51, 0.75, 0.11, 0.08, 0.00, 0.00, 0.00, 0.00,
                  0.00, 0.00, 0.00, 0.00, 0.14, 0.05, 0.38, 0.47, 0.34, 0.54, 0.09, 0.00, 0.00, 0.09,
                  0.12, 0.06, 0.19, 0.18, 0.25, 0.43, 0.11, 0.20, 0.72, 0.52, 0.43, 0.23, 0.09, 0.46,
                  0.02, 0.02, 0.06, 0.04, 0.08, 0.09, 0.05, 0.09, 0.25, 0.32, 0.42, 0.26, 0.13, 0.18,
                  0.05, 0.06, 0.00, 0.03, 0.09, 0.22, 0.11, 0.08, 0.08, 0.07, 0.35, 0.27, 0.04, 0.06,
                  0.18, 0.20, 0.05, 0.01, 0.08, 0.12, 0.57, 0.28, 0.16, 0.07, 0.10, 0.17, 0.25, 0.19,
                  0.15, 0.26, 0.19, 0.08, 0.05, 0.10, 0.15, 0.43, 0.35, 0.11, 0.04, 0.03, 0.31, 0.28,
                  0.16, 0.23, 0.31, 0.24, 0.11, 0.08, 0.15, 0.23, 0.54, 0.34, 0.14, 0.09, 0.21, 0.66,
                  0.00, 0.00, 0.28, 0.29, 0.23, 0.12, 0.11, 0.17, 0.25, 0.39, 0.33, 0.14, 0.00, 0.10,
                  0.02, 0.03, 0.07, 0.05, 0.24, 0.36, 0.07, 0.06, 0.00, 0.10, 0.43, 0.19, 0.22, 0.12,
                  0.03, 0.04, 0.00, 0.00, 0.16, 0.08, 0.22, 0.00, 0.23, 0.08, 0.16, 0.26, 0.35, 0.15,
                  0.02, 0.06, 0.02, 0.03, 0.06, 0.08, 0.21, 0.12, 0.13, 0.00, 0.02, 0.10, 0.40, 0.38,
                  0.00, 0.04, 0.01, 0.05, 0.00, 0.05, 0.10, 0.09, 0.06, 0.00, 0.00, 0.05, 0.28, 0.58))

# Define age categories

age_category = c("0-4", "5-9", "10-14", "15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60-64", "65+")

# Save

# Shanghai
tmp = data.table(age_index = rep(age_category, length(age_category)),
                 age_contact = rep(age_category, each = length(age_category) ),
                 count_pre = as.vector(pre_S),
                 count_post = as.vector(post_S) )
tmp[, city := "Shanghai"]

# Wuhan
tmp1 = data.table(age_index = rep(age_category, length(age_category)),
                  age_contact = rep(age_category, each = length(age_category) ),
                  count_pre = as.vector(pre_W),
                  count_post = as.vector(post_W) )
tmp1[, city := "Wuhan"]
tmp = rbind(tmp, tmp1)

tmp[, age_index := factor(age_index, levels = age_category)]
tmp[, age_contact := factor(age_contact, levels = age_category)]

saveRDS(tmp, file = file.path(indir, "usa", "data", "meta", "estimate_contact_intensities_outbreak_Zhang.rds"))

#
# Contact matrices from Jarvis et al. BMC https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-020-01597-8

library('socialmixr')

pre_UK = contact_matrix(polymod, countries = "United Kingdom", age.limits = c(0, 5, 18, 30, 40, 50, 60, 70))$matrix

post_UK = matrix(ncol = 8, nrow = 8, byrow = T,
                 c(0.525697998, 0.402583748, 0.368722035, 0.538440709, 0.194780472, 0.173757609, 0.098934288,	0.047254293,
                   0.14245221,	 2.430861989, 0.42609208	, 0.514546595, 0.48340504	, 0.188256272, 0.12084115	, 0.111607242,
                   0.137760561, 0.449900339, 1.194867894, 0.483464468, 0.480124784,	0.346225317, 0.184960293,	0.121814092,
                   0.209770231, 0.566523235, 0.504132528, 0.757773889, 0.52903225	, 0.346035983, 0.239739093,	0.08818435,
                   0.077825866, 0.545854342, 0.51346029	, 0.542568686, 0.795064526, 0.362278814, 0.233184095,	0.237601009,
                   0.079827148, 0.244423681, 0.425735426, 0.408058164, 0.416553864, 0.509931872, 0.264286413,	0.204649283,
                   0.058580903, 0.202213955, 0.293131116, 0.364369583, 0.345564886, 0.340625601, 0.321611885,	0.271976601,
                   0.023550101, 0.157192176, 0.162488722, 0.11280726 , 0.296361125, 0.22200101,	 0.228914771,	0.369024457))

# Define age categories

age_category = c("0-4", "5-17", "18-29", "30-39", "40-49", "50-59", "60-69", "70+")

tmp = data.table(age_index = rep(age_category, length(age_category)),
                 age_contact = rep(age_category, each = length(age_category) ),
                 count_pre = as.vector(pre_UK),
                 count_post = as.vector(post_UK) )

tmp[, age_index := factor(age_index, levels = age_category)]
tmp[, age_contact := factor(age_contact, levels = age_category)]

saveRDS(tmp, file = file.path(indir, "usa", "data", "meta", "estimate_contact_intensities_outbreak_UK.rds"))


text:
#' Search Bills
#'
#' Use this request to search the title and full text of legislation by keyword
#' to get the 20 most recent bills. Searches cover House and Senate bills from
#' the 113th Congress through the current Congress (116th). If multiple words
#' are given (e.g. query=health care) the search is treated as multiple keywords
#' using the OR operator. Quoting the words (e.g. query="health care") makes
#' it a phrase search. Search results can be sorted by date (the default) or by
#' relevance, and in ascending or descending order.
#'
#' @param query search text
#' @param sort sort by lastest major action date or relevance score
#' @param dir sort in descending or ascending order
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bills_search("health care")
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#search-bills}
bills_search <- function(query, sort = c("date", "_score"),
                         dir = c("desc", "asc"), page = 1L) {
  sort <- match.arg(sort)
  dir <- match.arg(dir)
  path <- "bills/search.json"
  params <- list(
    query = query,
    sort = sort,
    dir = dir,
    offset = offset_from_page(page)
  )

  # get results
  results <- congress_api(path, params, class = "bills_search")
  results
}

#' Get Recent Bills
#'
#' Use this request type to get summaries of the 20 most recent bills by type.
#' For the current Congress, “recent bills” can be one of four types
#' (see the descriptions below). For previous Congresses, “recent bills” means
#' the last 20 bills of that Congress. In the responses, an \code{active} value of true
#' means that the bill has seen action beyond introduction and committee
#' referral. Requests include a \code{chamber} value; to get recent bills from the
#' House and Senate, use both as the value. You can paginate through bills using
#' the \code{offset} querystring parameter that accepts multiples of 20. Bills before
#' the 113th Congress (prior to 2013) have fewer attribute values than those
#' from the 113th Congress onward, because the more recent bill data comes from
#' the bulk data provided by the Government Publishing Office. Details for the
#' older bills came from scraping Thomas.gov, the former congressional site of
#' the Library of Congress.
#'
#' @param chamber house or senate
#' @param congress 105-116
#' @param type introduced, updated, active, passed, enacted, or vetoed
#' @param page page of results
#'
#' @examples
#' \dontrun{
#' bills_recent("house", 116, "enacted")
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#search-bills}
bills_recent <- function(chamber = c("house", "senate"),
                         congress = 116,
                         type = c("introduced", "updated", "active", "passed", "enacted", "vetoed"),
                         page = 1L) {

  chamber <- match.arg(chamber)
  type <- match.arg(type)

  path <- sprintf("%s/%s/bills/%s.json", congress, chamber, type)
  params <- list(
    offset = offset_from_page(page)
  )

  # get results
  results <- congress_api(path, params, class = "bills_recent")
  results
}

#' Recent Bills by a Specific Subject
#'
#' Use this request type to get the 20 most recently updated bills for a specific
#' legislative subject. Results can include more than one Congress.
#'
#' @param subject A slug version of a legislative subject
#' @param page page of results
#' @keywords internal
#'
#' @examples
#' \dontrun{
#' bills_recent_subject("meat")
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-recent-bills-by-a-specific-subject}
bills_recent_subject <- function(subject, page = 1L) {
  path <- sprintf("bills/subjects/%s.json", subject)
  #params <- list(offset = offset_from_page(page))

  # get results
  results <- congress_api(path, class = "bills_recent_subject", extract = FALSE)
  results
}

#' Search a Specific Bill Subject
#'
#' Use this request type to search for bill subjects that contain a specified
#' term.
#'
#' @param query a word or phrase to search
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bills_search_subjects("climate")
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-a-specific-bill-subject}
bills_search_subjects <- function(query, page = 1L) {
  path <- "bills/subjects/search.json"
  params <- list(
    query = query,
    offset = offset_from_page(page)
  )

  # get results
  results <- congress_api(path, params, class = "bills_search_subjects")
  results
}

#' Get Upcoming Bills
#'
#' Use this request type to get details on bills that may be considered by the
#' House or Senate in the near future, based on scheduled published or announced
#' by congressional leadership. The bills and their potential consideration are
#' taken from the \href{http://docs.house.gov/floor/}{House Majority Leader} and
#' floor updates from
#' \href{http://www.republicans.senate.gov/public/index.cfm/floor-updates}{Senate Republicans}.
#' The responses include a \code{legislative_day} attribute which is the earliest the
#' bills could be considered, and a \code{range} attribute that indicates whether the
#' bill information comes from a weekly schedule or a daily one. Combine the two
#' for the best sense of when a bill might come up for consideration. For Senate
#' bills, the response includes a \code{context} attribute reproducing the sentence that
#' includes mention of the bill. These responses omit bills that have not yet been
#' assigned a bill number or introduced, and additional bills may be considered
#' at any time.
#'
#' @param chamber house or senate
#' @export
#'
#' @examples
#' \dontrun{
#' bills_upcoming("house")
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-upcoming-bills}
bills_upcoming <- function(chamber = c("house", "senate")) {
  chamber <- match.arg(chamber)
  path <- sprintf("bills/upcoming/%s.json", chamber)

  # get results
  results <- congress_api(path, class = "bills_upcoming")
  results
}

#' Get a Specific Bill
#'
#' Use this request type to get details about a particular bill, including
#' actions taken and votes. The attributes \code{house_passage_vote} and
#' \code{senate_passage_vote} are populated (with the date of passage) only upon
#' successful passage of the bill. Bills before the 113th Congress (prior to 2013)
#' have fewer attribute values than those from the 113th Congress onward,
#' because the more recent bill data comes from the bulk data provided by the
#' Government Publishing Office. Details for the older bills came from scraping
#' Thomas.gov, the former congressional site of the Library of Congress.
#'
#' @param bill_id a bill slug, for example hr4881 - these can be found in the recent bill response.
#' @param congress 105-116
#' @export
#'
#' @examples
#' \dontrun{
#' bill("hr502", 116)
#' bill("s1436", 116)
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-a-specific-bill}
bill <- function(bill_id, congress = 116L) {
  path <- sprintf("%s/bills/%s.json", congress, bill_id)

  # get results
  results <- congress_api(path, class = "bill")
  results
}

#' Get Amendments for a Specific Bill
#'
#' Use this request type to get Library of Congress-assigned subjects about a
#' particular bill. This request returns the 20 most recent results and supports
#' paginated requests.
#'
#' @param bill_id a bill slug, for example hr4881 - these can be found in the recent bill response.
#' @param congress 105-116
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bill_amendments("hr1628", 115)
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-amendments-for-a-specific-bill}
bill_amendments <- function(bill_id, congress = 116L, page = 1L) {
  path <- sprintf("%s/bills/%s/amendments.json", congress, bill_id)
  params <- list(offset = offset_from_page(page))

  # get results
  results <- congress_api(path, class = "bill_amendments")
  results
}

#' Get Subjects for a Specific Bill
#'
#' Use this request type to get Library of Congress-assigned subjects about a
#' particular bill. This request returns the 20 most recent results and supports
#' paginated requests.
#'
#' @param bill_id a bill slug, for example hr4881 - these can be found in the recent bill response.
#' @param congress 105-116
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bill_subjects("hr2810", 115)
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-subjects-for-a-specific-bill}
bill_subjects <- function(bill_id, congress = 116L, page = 1L) {
  path <- sprintf("%s/bills/%s/subjects.json", congress, bill_id)
  params <- list(offset = offset_from_page(page))

  # get results
  results <- congress_api(path, class = "bill_subjects")
  results
}

#' Get Related Bills for a Specific Bill
#'
#' Use this request type to get Library of Congress-identified related bills for
#' a particular bill. This request returns the 20 most recent results and
#' supports paginated requests.
#'
#' @param bill_id a bill slug, for example hr4881 - these can be found in the recent bill response.
#' @param congress 105-116
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bill_related("hr3219", 115)
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-related-bills-for-a-specific-bill}
bill_related <- function(bill_id, congress = 116L, page = 1L) {
  path <- sprintf("%s/bills/%s/related.json", congress, bill_id)
  params <- list(offset = offset_from_page(page))

  # get results
  results <- congress_api(path, class = "bill_related")
  results
}

#' Get Cosponsors for a Specific Bill
#'
#' Get information about the cosponsors of a particular bill
#'
#' @param bill_id a bill slug, for example hr4881 - these can be found in the recent bill response.
#' @param congress 105-116
#' @param page page of results
#' @export
#'
#' @examples
#' \dontrun{
#' bill_cosponsors("hr4249", 114)
#' }
#'
#' @references \url{https://projects.propublica.org/api-docs/congress-api/bills/#get-cosponsors-for-a-specific-bill}
bill_cosponsors <- function(bill_id, congress = 116L, page = 1L) {
  path <- sprintf("%s/bills/%s/cosponsors.json", congress, bill_id)
  params <- list(offset = offset_from_page(page))

  # get results
  results <- congress_api(path, class = "bill_cosponsors")
  results
}

max_stars_repo_name: zamorarr/ppcongress

max_issues_repo_name: zamorarr/ppcongress

max_forks_repo_name: zamorarr/ppcongress


text:
if (length(commandArgs(TRUE)) == 0) {
  system("Rscript barplot.r -h", intern = F)
  q("no")
}

# load packages that are provided in the conda env
options(show.error.messages = F,
       error = function() {
           cat(geterrmessage(), file = stderr())
                           q("no", 1, F)
       }
)
loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
warnings()
library(optparse)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)

option_list <- list(
  make_option(
    c("-i", "--input"),
    default = NA,
    type = "character",
    help = "Input file that contains count data (no header)"
  ),
  make_option(
    c("-o", "--barplot"),
    default = NA,
    type = "character",
    help = "PDF output file"
  )
)

opt <- parse_args(OptionParser(option_list = option_list),
                 args = commandArgs(trailingOnly = TRUE))

annotations <- read.delim(opt$input, header = F)
colnames(annotations) <- c("sample", "class", "percent_of_reads", "total")
annotations$percent <- round(annotations$percent_of_reads / annotations$total * 100,
                            digits = 2)
## ggplot2 plotting

# Define the number of colors you want
mycolors <- c("#e6194b", "#3cb44b", "#ffe119", "#4363d8", "#f58231",
                              "#911eb4", "#46f0f0", "#f032e6", "#bcf60c",
                              "#008080", "#e6beff", "#9a6324", "#fffac8", "#800000",
                              "#aaffc3", "#808000", "#ffd8b1", "#000075", "#808080")
ggtitle("Class proportions")
ggplot(annotations, aes(x = total / 2, y = percent_of_reads, fill = class, width = total)) +
       geom_bar(position = "fill", stat = "identity") +
       facet_wrap(~sample, ncol = 3) +
       geom_label_repel(aes(label = percent), position = position_fill(vjust = 0.5),
                        size = 2, show.legend = F) +
       coord_polar(theta = "y") +
       labs(x = "Class fractions (%)") +
       scale_fill_manual(values = mycolors) +
       theme(axis.text = element_blank(),
             axis.ticks = element_blank(),
             panel.grid  = element_blank(),
             axis.title.y = element_blank(),
             legend.position = "bottom") +
       geom_text(aes(x = total / 2, y = .5, label = paste(round(total / 1000000, digits = 3),
                     "M"),
                 vjust = 4, hjust = -1), size = 2)
ggsave(file = opt$barplot, device = "pdf")

max_stars_repo_name: ARTbio/tools-artbio

max_issues_repo_name: ARTbio/tools-artbio

max_forks_repo_name: ARTbio/tools-artbio


text:
#
# trajectory with selected 1000 cells
#
x=VariableFeatures(d)
y=d$RNA@data[x,]
e=apply(y,2,mean)

j=match( data.table(id=colnames(y), e=e, c=d$orig.ident)[order(-e,c),head(.SD,100),by=c]$id, colnames(y) )
w=Embeddings(d,"harmony")[j,]
c=d$orig.ident[j]

my_color <- brewer.pal(8,"Set2")
names(my_color) <- unique(names(table(c)))

        ggplot()+ geom_point(aes(x = x, y = y,col=color)) +
        scale_color_manual(labels=names(my_color),values=as.character(my_color))

x = data.frame(x=w[,1],y=w[,2],color=my_color[c])
plot(x$x,x$y, col = my_color[c], pch=16, asp = 1)
legend("topright",legend = names(my_color), fill = my_color)
lines(sds, lwd=2, type = 'lineages', col = c("black"))

plot(reducedDim(sds), col = my_color[c], pch = 16, cex = 0.5)
lines(sds, lwd = 2, type = 'lineages', col = 'black')

library(slingshot)
sds <- slingshot(w, clusterLabels = c, start.clus = "Lee-1", stretch = 0)

library(ggbeeswarm)
data.frame(x=slingPseudotime(sds)[,1],y=factor(c),c=my_color[c]) %>%
        filter(!is.na(x)) %>%
        ggplot(aes(x=x,y=y,colour=c)) +
         scale_color_manual(labels=names(my_color),values=as.character(my_color)) +
         geom_quasirandom(groupOnX = FALSE) + theme_classic()

max_stars_repo_name: hmgene/mudcookies

max_issues_repo_name: hmgene/mudcookies

max_forks_repo_name: hmgene/mudcookies


text:
KmeansGap=function(dat,nozeros=FALSE,multiD=0,numnulls=100,verbose=1,mink=1,maxk=NULL,nulltype='shuffle',nstartingpoints=100,YanYe=0,plot=0,line90=1){
	set.seed(0)
	if(class(dat)!='data.frame') stop('Valid input to KmeansGap must be a data frame')
	if(!multiD & !'trait'%in%names(dat)) stop('Input data frame for the 1d case must have a "trait" column')
	if(!'N'%in%names(dat)) stop('Input data frame must have a "N" column')

	## prepare data for analysis
	dat=dat[apply(dat,1,function(v) all(!is.na(v))),]
	if(nozeros) dat=dat[dat$N>0,]
	dat$N=dat$N/min(1,min(dat$N[dat$N>0]))
	if(!multiD) dat=dat[order(dat$trait),]

	traitdata=if(multiD) subset(dat,select=-N) else data.frame(trait=dat$trait)

	if(is.null(maxk)) maxk=round(min(20,sum(dat$N>0)/2))
	kvec=mink:maxk

	library(plyr)
	rep.row=function(r,a) colwise(function(x) rep(x,a))(r)
	kms=function(mat,n,k){
		if(ncol(mat)>1) com=rep.row(r=mat,a=n) else com=rep(mat$trait,n)
		if(k==1){
			z=kmeans(com,centers=1,iter.max=100)
		}else{
			nullcenters=apply(unique(mat),2,function(v) sample(seq(min(v),max(v),l=k)))
			z=try(kmeans(com,centers=nullcenters,iter.max=100),silent=TRUE)
			if(!is.null(attr(z,'condition'))) z=kmeans(com,centers=k,nstart=nstartingpoints,iter.max=100)
		}
		return(z)
	}

	Dispersion=function(dtf,weight){
		if(weight==0) return(dtf$tot.withinss)							## Unweighted sum of dist^2
		if(weight==1) return(with(dtf,sum(withinss/(2*size))))			## Tibshirani et al 2001
		if(weight==2) return(with(dtf,sum(withinss/(2*size*(size-1)))))	## Correction from Yan & Ye 2007
	}

	Wk=sapply(kvec,function(k){ 	## Wk = within-cluster sum of squares (ie squared distances); D in Tibshirani et al 2001
		mod=kms(mat=traitdata,n=dat$N,k)
		Dispersion(dtf=mod,weight=YanYe)
	})

	nullWk=sapply(kvec,function(k){
		if(verbose) cat(".",if(k%%10==0 | k==maxk) paste(k,"\n"))
		if(nulltype=='shuffle'){
			ndat=sapply(seq(numnulls),function(run) sample(dat$N))
			return(apply(ndat,2,function(N){
				mod=kms(mat=traitdata,n=N,k)
				Dispersion(dtf=mod,weight=YanYe)
			}))
		}else if(nulltype=='spline'){
			ndat=rnorm(length(dat$N),mean=dat$N,sd=sd(dat$N)); ndat[ndat<0]=0
			return(apply(ndat,2,function(N){
				mod=kms(mat=traitdata,n=N,k)
				Dispersion(dtf=mod,weight=YanYe)
			}))
		}else{
			return(apply(ndat,2,function(N){
				tmp=apply(as.matrix(traitdata),2,function(v) runif(length(v),min=min(v),max=max(v)))
				mod=kms(mat=as.data.frame(tmp),n=N,k)
				Dispersion(dtf=mod,weight=YanYe)
			}))
		}
	})

	logWk=log(Wk)										## vector length(kvec)
	lognullWk=log(nullWk)								## matrix numnulls x length(kvec)
	ElogWk=apply(lognullWk,2,mean,na.rm=TRUE)			## vector length(kvec)
	Gapk=ElogWk-logWk									## vector length(kvec)
	nullGapk=t(ElogWk-t(lognullWk))						## matrix numnulls x length(kvec)

	maxgap=max(Gapk,na.rm=TRUE)							## scalar
	maxnullgap=apply(nullGapk,1,max,na.rm=TRUE)			## vector numnulls
	kmaxnullgap=apply(nullGapk,1,which.max)				## vector numnulls
	kmaxnullquant=merge(data.frame(kmax=kvec),ddply(data.frame(kmax=kmaxnullgap,max=maxnullgap),.(kmax),function(v) quantile(v$max,.95)),by='kmax',all.x=TRUE)
	kmaxnullquant90=merge(data.frame(kmax=kvec),ddply(data.frame(kmax=kmaxnullgap,max=maxnullgap),.(kmax),function(v) quantile(v$max,.9)),by='kmax',all.x=TRUE)
	khat=kvec[which.max(Gapk)]							## scalar
	z.score=(maxgap-mean(maxnullgap))/sd(maxnullgap)	## scalar
	p.value=sum(maxnullgap>maxgap)/length(maxnullgap)	## scalar

	mngap=apply(nullGapk,2,mean,na.rm=TRUE); 			## mean  of the null gaps --- vector length(kvec)
	sdgap=apply(nullGapk,2,sd,na.rm=TRUE); 				## stdev of the null gaps --- vector length(kvec)
	nullquant=apply(nullGapk,2,function(vec) as.numeric(quantile(vec,.95,na.rm=TRUE)))	## 95% quantile of the null gaps --- vector length(kvec)
	nullquant90=apply(nullGapk,2,function(vec) as.numeric(quantile(vec,.9,na.rm=TRUE)))	## 90% quantile of the null gaps --- vector length(kvec)
	maxnullquant=quantile(maxnullgap,.95)				## 95% quantile of the max null gaps --- scalar
	maxnullquant90=quantile(maxnullgap,.9)				## 90% quantile of the max null gaps --- scalar

	if(plot){
		plot(0,t='n',xlim=c(min(kvec),max(kvec)),ylim=c(0,max(Gapk,maxnullquant)),las=1,xlab='No. clusters',ylab='Gap')
		polygon(x=c(kvec,rev(kvec)),y=c(mngap,rev(nullquant)),col='grey90',border=NA)
		polygon(x=c(kvec,rev(kvec)),y=c(mngap,rev(nullquant90)),col='grey80',border=NA)
		lines(kvec,rep(maxnullquant,length(kvec)),col=2)
		if(line90) lines(kvec,rep(maxnullquant90,length(kvec)),col=2,lty=2)
		points(kvec,Gapk,t='o',pch=20,lwd=2)
		polygon(x=c(kvec,rev(kvec)),y=c(mngap,rep(-maxnullquant,length(kvec))),col='white',border=NA)
		box()
	}

	if(YanYe==0) disp='D'; if(YanYe==1) disp='D/2n'; if(YanYe==2) disp='D/(2n(n-1))'

	return(list(
		data=data.frame(
			k=kvec,
			gap=Gapk,
			Egap=mngap,
			sdgap=sdgap,
			nullquant=nullquant,
			nullquant90=nullquant90,
			kmaxnullquant=kmaxnullquant[,2],
			kmaxnullquant90=kmaxnullquant90[,2],
			logWk=logWk,
			ElogWk=ElogWk
		),
		khat=khat,
		maxgap=maxgap,
		maxnullquant=maxnullquant,
		maxnullquant90=maxnullquant90,
		mink=mink,
		maxk=maxk,
		z.score=z.score,
		p.value=p.value,
		dispersion=disp
	))
}


text:
####################################################################################################################################
#                                                        Crime data                                                                #
####################################################################################################################################

library(data.table)
library(sf)
library(RCurl)
library(RPostgreSQL)

source("~/Dropbox/git/postgis/postgis_functions.r")

####### Attention! In the end I didn't use crime rate per catchment but per cluster, so the following was not used #################

##### 1. Download data #############################################################################################################
data_url <- "https://data.police.uk/data/archive/"
# down_dir <- "~/Dropbox/liverpool/retail_typology/data/crime"
# dir.create(down_dir)
year <- 2016
for (m in 1:12){
  m <- ifelse(m < 10, paste0("0", m), as.character(m))
  full_data_url <- paste0(data_url, year, "-", m, ".zip")
  temp <- tempfile()
  download.file(full_data_url, temp)
  csv_names <- unzip(temp, list = T)
  unzip(temp, junkpaths = T, exdir = "~/temp", overwrite = T)
  unlink(temp)
}

csvs <- list.files("~/temp")
con <- dbConnect(PostgreSQL(), dbname = "crime", user = "user", password = "pwd")
dbExecute(con, "CREATE TABLE crime(crime_id TEXT,	month TEXT,	reported_by TEXT,	falls_within TEXT, longitude NUMERIC,
                latitude NUMERIC,	location TEXT, lsoa_code TEXT,	lsoa_name TEXT,	crime_type TEXT, last_outcome_category TEXT, context TEXT)")
for (csv in csvs){
  system(paste0("psql -U postgres -c \"copy crime FROM '~/temp/", csv, "' CSV HEADER\" crime"))
}

dbExecute(con, "ALTER TABLE crime ADD COLUMN geog geography(POINT,4326)")
dbExecute(con, "UPDATE crime SET geog = ST_GeogFromText('SRID=4326;POINT(' || longitude || ' ' || latitude || ')')")
dbExecute(con, "CREATE INDEX ON crime USING gist(geog)")
dbExecute(con, "CREATE INDEX ON crime(crime_type)")
# import_or_append(con = con, working_dirs = "~/Dropbox/liverpool/GIS_Data/LADs_2011", shp_names = "Local_UnitaryAuthority.shp", srid = 27700, table_name = "lads")
import_or_append(con = con, working_dirs = "/home/michalis/Dropbox/liverpool/GIS_Data/OAs_2011", shp_names = "oas11.shp", srid = 27700, table_name = "oas")
dbExecute(con, "create index on oas(code)")

dbExecute(con, "CREATE TABLE oas_84 (gid serial PRIMARY KEY, oa11 TEXT, geog geography(MultiPolygon, 4326))")
dbExecute(con, "INSERT INTO oas_84(oa11, geog) SELECT code, ST_Transform(geom, 4326) FROM oas")
dbExecute(con, "CREATE INDEX ON oas_84(code)")
dbExecute(con, "CREATE INDEX ON oas_84 USING gist(geog)")

import_or_append(con = con, working_dirs = "/home/michalis/Dropbox/liverpool/boundaries/new_boundaries", shp_names = "all_boundaries_fin.shp", srid = 27700, table_name = "boundaries")
dbExecute(con, "CREATE INDEX ON boundaries(id)")

dbExecute(con, "CREATE TABLE boundaries84_buf100(id TEXT, geog geography(Polygon, 4326))")
dbExecute(con, "INSERT INTO boundaries84_buf100(id, geog) SELECT id, ST_Transform(ST_Buffer(geom, 100), 4326) FROM boundaries")
dbExecute(con, "CREATE INDEX ON boundaries84_buf100(id)")
dbExecute(con, "CREATE INDEX ON boundaries84_buf100 USING gist(geog)")

dbExecute(con, paste("CREATE TABLE count_crimes AS",
                     "SELECT t1.id, ST_Area(t1.geog), count(t2.*) FROM boundaries84_buf100 t1 LEFT JOIN crime t2 ON ST_Covers(t1.geog, t2.geog)",
                     "WHERE t2.crime_type NOT IN ('Burglary', 'Other theft') GROUP BY t1.id, t1.geog"))

# dbExecute(con, paste("CREATE TABLE count_crimes AS",
#                      "SELECT t1.id, count(t2.*) FROM boundaries84 t1 LEFT JOIN crime t2 ON ST_DWithin(t1.geog, t2.geog, 100, false)",
#                      "WHERE t2.crime_type NOT IN ('Burglary', 'Other theft') GROUP BY id"))
count_crimes <- dbGetQuery(con, "select * from count_crimes")
names(count_crimes)[2:3] <- c("buffer_retail_area", "count_crime")
count_crimes$buffer_retail_area <- count_crimes$buffer_retail_area / 10^4
count_crimes$crime_density <- count_crimes$count_crime / count_crimes$buffer_retail_area
write.csv(count_crimes, "~/Dropbox/liverpool/retail_typology/data/crime/count_crime.csv")

#### 2. Crime/pop per cluster id #########################################################################################################
crimes <- fread("~/Dropbox/liverpool/retail_typology/data/crime/count_crime.csv")[,2:3]
all_probs <- fread("~/Dropbox/liverpool/retail_typology/data/catchments_oa_lookup.csv")

# sanity check
sum(unique(all_probs$oa11) %in% crimes$code) # 172636
sum(!unique(all_probs$oa11) %in% crimes$code) # 46710

table(substr(unique(all_probs[oa11 %in% crimes$code, ]$oa11), 1, 1))
#      E      S      W
# 163489    283   8864
table(substr(unique(all_probs[!oa11 %in% crimes$code, ]$oa11), 1, 1))
#    E     S     W
# 4890 41676   144

probs_crime <- crimes[all_probs, on = c(code = "oa11")]
probs_crime[is.na(count), count := 0]

crime_rate <- catchments_crime[, list(sum(count) / sum(pop), sum(count)), keyby = "cluster_id"]
setnames(crime_rate, c("V1", "V2"), c("crime_per_pop", "total_crime"))

fwrite(crime_rate, "~/Dropbox/liverpool/retail_typology/analysis/indicators/crime_rate.csv")


text:
 ZPPCON Example Program Results

 Estimate of condition number =  1.51E+02


text:
## -------->>  [[file:../nstandr.src.org::*Package documentation][Package documentation:1]]
#' @details
#' nstandr package standardizes organizational names
#'     mainly using procedures described in Thoma et al. (2010) and
#'     Magerman, Looy, Bart, & Song (2006) but not only.  This is work
#'     in progress. Please, file an issues or suggestion if you have
#'     any.  The main function is [standardize()].
#' @keywords internal
"_PACKAGE"
## --------<<  Package documentation:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_toupper][standardize_toupper:1]]
##' @eval attr(standardize_toupper, "description")
##'
##' @param x data
##'
##' @inheritDotParams standardize_options
##'
##' @return updated data (as data.table)
##' @export
standardize_toupper <- function(x, ...) {
    get_target(x) |>
        toupper() |>
        inset_target(x)
}

attr(standardize_toupper, "description") <-
"Uppercases vector of interest in the object (table)"
## --------<<  standardize_toupper:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_remove_brackets][standardize_remove_brackets:1]]
##' @eval attr(standardize_remove_brackets, "@title")
##' @param x object (table)
##' @inheritDotParams standardize_options
##' @return updated object
##'
##' @export
standardize_remove_brackets  <- function(x, ...) {
    get_target(x) |>
        stringi::stri_replace_all_regex("<[^<>]*>|\\([^()]*\\)|\\{[^{}]*\\}|\\[[^\\[\\]]*\\]", "") |>
        inset_target(x)
}

attr(standardize_remove_brackets, "@title") <- "Removes brackets and content in brackets"
## --------<<  standardize_remove_brackets:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_remove_quotes][standardize_remove_quotes:1]]
##' Removes double quotes
##'
##' @param x an object
##' @inheritDotParams standardize_options
##' @return updated object
##' @export
standardize_remove_quotes <- function(x, ...) {
        get_target(x) |>
          stringi::stri_replace_all_regex("\"", "") |>
          inset_target(x)
}
## --------<<  standardize_remove_quotes:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_dehtmlize][standardize_dehtmlize:1]]
#' Converts HTML characters to UTF-8
#'
#' The method is about 1/3 faster than htmlParse but it is still quite slow
#' @param x object (table)
#' @param as_single_string If set then collapse characters in the main column of the `x` (i.e., `x.col`) as to a single string. It will increase performance (at least for relatively short tables). Default is FALSE
#' @param as_single_string_sep delimiter for collapsed strings to uncollapse it later. Default is "#_|".
#' @param use_read_xml If set the it will parse XML. Default is FALSE which means it parses HTML
#' @inheritDotParams standardize_options
#' @return updated object
#' @references http://stackoverflow.com/questions/5060076
#'
#' @export
standardize_dehtmlize <- function(x
                              , as_single_string = FALSE
                              , as_single_string_sep = "#_|"
                              , use_read_xml = FALSE
                              , ...) {
    x_vector <- get_target(x)
    if(as_single_string) {
        x_vector <- paste0(x_vector, collapse = as_single_string_sep)
        x_vector <- paste0("<x>", x_vector, "</x>")
        x_vector <-
            (if(use_read_xml) {
                 xml2::read_xml(x_vector)
             } else {
                 xml2::read_html(x_vector)
             }) |> xml2::xml_text()
        strsplit(x_vector, as_single_string_sep, fixed = TRUE)[[1]]
    } else {
        sapply(x_vector, \(str) {
            str <- paste0("<x>", str, "</x>")
            (if(use_read_xml) {
                 xml2::read_xml(str)
             } else {
                 xml2::read_html(str)
             }) |> xml2::xml_text()
        }, USE.NAMES = FALSE)
    } |> inset_target(x)
}
## --------<<  standardize_dehtmlize:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_detect_enc][standardize_detect_enc:1]]
#' Detects string encoding
#' @param x object
#' @param output_codes_col_name Same as in [detect_patterns()]
#' @param return_only_codes Same as in [detect_patterns()]
#' @param ...
#' @inheritDotParams standardize_options
#' @return updated object
#'
#' @export
standardize_detect_enc <- function(x
                               , output_codes_col_name = "{col_name_}encoding"
                               , return_only_codes = FALSE
                               , ...) {
    available_enc_list <- iconvlist()
    x_vector <- get_target(x) |>
        stringi::stri_enc_detect() |>
        lapply(function(enc) {
            enc <- enc[["Encoding"]]
            first_ok_enc <- which(enc %in% available_enc_list)[1]
            if(length(first_ok_enc) == 0) ""
            else enc[[first_ok_enc]]
        }) |> unlist()
    if(return_only_codes) {
        x_vector
    } else {
        inset_target(x_vector
                   , x
                   , output_placement = "omit"
                   , output_copy_col_name = output_codes_col_name
                   , append_output_copy = TRUE)
    }
}
## --------<<  standardize_detect_enc:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_toascii][standardize_toascii:1]]
#' Translates non-ascii symbols to its ascii equivalent
#'
#' @param x String to translate
#' @param detect_encoding Detect encoding of individual elements (slower). Allows to work with mixed encodings.
#' @inheritDotParams standardize_options
#'
#' @export
standardize_toascii <- function(x
                            , detect_encoding = FALSE
                            , ...) {
  str <- get_target(x)
  utf <- nstandr_patterns_ascii$utf |> paste(collapse = "")
  ascii <- nstandr_patterns_ascii$ascii |> paste(collapse = "")
  (if(detect_encoding) {
       mapply(
           \(name, enc) chartr(utf, ascii, iconv(name, from = enc, to = "UTF-8", sub = ""))
         , str
         , standardize_detect_enc(str, return_only_codes = TRUE)
         , SIMPLIFY = FALSE, USE.NAMES = FALSE) |>
           unlist() |>
           iconv(to = "ASCII", sub = "")
   } else {
       chartr(utf, ascii, enc2utf8(str)) |>
           iconv(to = "ASCII", sub = "")
   }) |> inset_target(x)
}
## --------<<  standardize_toascii:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize_x_split][standardize_x_split:1]]
##' Splits the object (table) in chunks by rows
##'
##' Convenient to apply some function to the table in chunks, e.g., if you want to add display of progress.
##'
##' @param x object or table
##' @param by number of rows to split by
##' @param len length of the table (nrow). If it is NULL then use x_length(x)
##'
##' @return List of (sub)tables
standardize_x_split <- function(x, by, len = NULL) {
    if(is.null(len)) len <- x_length(x)
    split(x, rep(seq(1, len %/% by +1)
               , each = by
               , length.out = len))
}
## --------<<  standardize_x_split:1 ends here

## -------->>  [[file:../nstandr.src.org::*standardize][standardize:1]]
#' Standardizes organizational names. Takes either vector or column in the table.
#'
#' @param x object (table)
#' @param procedures Named list of procedures (closures) to apply to x. If we need to pass arguments to some of the procedures it can be done by specifying sub-list where the first element is procedure and the rest its arguments. Names of the list elements are used for progress messages. Procedures can also be passed as data.frame in which case it will be converted to list of procedures with `standardize_make_procedures_list` (see its help for the correct format of data.frame with procedures). Default is `nstandr_default_procedures_table`
#' @param show_progress Whether to report progress percentage. Default is TRUE
#' @param nrows_min_to_show_progress The minimum number of rows the x should have for automatic progress estimation. If x has less rows no progress will be shown. Default is 10^5
#' @param progress_step_nrows If set it will divide the x into chunk of this amount of rows. Default is NULL.
#' @param progress_step_in_percent Number of percents that represent one step in progress. Value should be between 0.1 and 50. Default is 1 which means it will try to chunk the x into 100 pieces.
#' @param progress_message_use_names Should we use names from `procedures` list to report progress. Default is TRUE.
#' @param quite Suppress all messages. Default is FALSE.
#' @param save_intermediate_x_to_var For debuging of standartization procedures. Saves intermediate results to this variable. If procedures finish without errors then the variable will be removed.
#' @inheritDotParams standardize_options
#'
#' @aliases make_standard_names make_std_names
#'
#' @return standardized data
#'
#' @export
standardize <- function(x
                    , procedures = nstandr_default_procedures_table
                    , show_progress = TRUE
                    , nrows_min_to_show_progress = 10^5
                    , progress_step_nrows = NULL
                    , progress_step_in_percent = 1
                    , progress_message_use_names = TRUE
                    , quite = FALSE
                    , save_intermediate_x_to_var = NULL
                    , ...) {
    checkmate::assert_string(save_intermediate_x_to_var, null.ok = TRUE)
    checkmate::assert_flag(show_progress)
    checkmate::assert_flag(quite)
    checkmate::assert_flag(progress_message_use_names)
    checkmate::assert_multi_class(procedures, classes = c("list", "data.frame"))
    if(is.data.frame(procedures)) {
        procedures <- standardize_make_procedures_list(procedures)
    }
    ## make format of the massages for procedures
    message_delimiter <- paste(c("\n", rep("-", 65), "\n"), collapse = "")
    message_init <- paste0("\nApplying standardization procedures:", message_delimiter)
    message_done  <- "\b\b\b\bDONE"
    progress_format <- "\b\b\b\b%3.0f%%"
    message_format <- "* %-60.60s...."
    message_fin <- paste0(message_delimiter, "Standardization is done!\n")
    ## ensure that x is either vector or data.table
    if(missing(x)) return(NULL)
    checkmate::assert_multi_class(x
                                , classes = c("character"
                                            , "numeric"
                                            , "integer"
                                            , "logical"
                                            , "data.frame"
                                            , "data.table"))
    ## check progress_step_in_percent
    checkmate::assert_number(progress_step_in_percent, lower = 0.1, upper = 50)
    checkmate::assert_number(progress_step_nrows
                           , lower = x_length(x)/1000
                           , upper = x_length(x)/2
                           , null.ok = TRUE)
    x_len <- x_length(x)
    ## Set progress_step_nrows
    progress_step_nrows <-
        if (show_progress && !quite) {
            if(x_len < nrows_min_to_show_progress) {
                NULL
            } else if(!is.null(progress_step_nrows)) {
                progress_step_nrows
            } else {
                round(x_len / (100 / progress_step_in_percent))
            }
        }
    ## Apply Procedures
    if(!quite) message(message_init)
    for(p in 1:length(procedures)) {
        ## get procedure function
        procedure_fun <- procedures[[p]][[1]]
        ## get procedure arguments (remove show_progressarg if it is there)
        procedure_args <- procedures[[p]][
            -c(1, which(names(procedures[[p]]) == "show_progress"))
        ]
        ## get procedure names
        procedure_name <-
            if(standardize_is_data_empty(names(procedures)[p]) | !progress_message_use_names) {
                procedure_fun
            } else {
                names(procedures)[p]
            }
        ## Anounce Procedure Name
        if(!quite) packageStartupMessage(
                       sprintf(message_format, procedure_name)
                     , appendLF = FALSE)
        ## Check if we need report progress:
        ## progress is set & progress = FALSE is absent in the arguments
        if(!is.null(progress_step_nrows) &
           !isFALSE(procedures[[p]]["show_progress"][[TRUE]])) {
            ## check if we need to split..
            if(!isTRUE(class(x) == "list")) {
                x <- standardize_x_split(x, progress_step_nrows)
            }
            ## set progress counter
            i <- 0; env <- environment()
            ## Apply procedure to list!
            x <- lapply(x, \(x_by) {
                if(!is.null(save_intermediate_x_to_var))
                    assign(save_intermediate_x_to_var, x, pos = 1)
                ## apply procedure fun with args
                x_by <- do.call(procedure_fun
                              , c(list(x_by), procedure_args))
                ## Increment progress counter
                assign("i", i + 100 * progress_step_nrows / x_len, envir = env)
                ## Anounce progress
                packageStartupMessage(
                    sprintf(progress_format, i)
                  , appendLF = FALSE)
                return(x_by)
            })
        } else {
            ## check if we need to rbindlist..
            if(isTRUE(class(x) == "list")) {
                if(is.atomic(x[[1]])) {
                    x <- unlist(x, use.names = FALSE)
                } else {
                    x <- rbindlist(x)
                }
            }
            if(!is.null(save_intermediate_x_to_var))
                assign(save_intermediate_x_to_var, x, pos = 1)
            ## Apply procedure fun with args!
            x <- do.call(procedure_fun, c(list(x), procedure_args))
        }
        ## Anounce DONE
        if(!quite) packageStartupMessage(message_done)
    }
    if(!quite) message(message_fin)
    ## return x
    if(isTRUE(class(x) == "list")) {
        if(is.atomic(x[[1]])) {
            x <- unlist(x, use.names = FALSE)
        } else {
            x <- rbindlist(x)
        }
    }
    ## remove intermediate saves if procedures finished without error
    remove
    if(!is.null(save_intermediate_x_to_var))
        rm(list = save_intermediate_x_to_var, pos = 1)
    return(x)
}

## Alias

#' @rdname standardize
#' @export
make_std_names <- standardize

#' @rdname standardize
#' @export
make_standard_names <- standardize
## --------<<  standardize:1 ends here

max_stars_repo_name: stasvlasov/nstandr

max_issues_repo_name: stasvlasov/nstandr

max_forks_repo_name: stasvlasov/nstandr


text:
load_install_packages <- function(list.of.packages) {
    new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
    if(length(new.packages)) install.packages(new.packages,repos="http://cran.cnr.Berkeley.edu/")
    for (i in list.of.packages) {
        print(sprintf("loading %s",i))
        library(i,character.only=T)
    }
}

package_lst <- c(
    "coin",
    "missForest",
    "knockoff",
    "glmnet",
    "ranger",
    "Xmisc",
    "doMC",
    "CoxBoost",
    "randomForestSRC",
    "party",
    "softImpute",
    "ranger",
    'IRkernel'
    )

load_install_packages(package_lst)

max_stars_repo_name: loramf/mlforhealthlabpub

max_issues_repo_name: loramf/mlforhealthlabpub

max_forks_repo_name: loramf/mlforhealthlabpub


text:
#source: tls.s
#as: -a64
#ld: -shared
#readelf: -WSsrl
#target: powerpc64*-*-*

There are [0-9]+ section headers, starting at offset 0x[0-9a-f]+:

Section Headers:
 +\[Nr\] Name +Type +Address +Off +Size +ES Flg Lk Inf Al
 +\[[ 0-9]+\] +NULL +0+ 0+ 0+ 0+ +0 +0 +0
 +\[[ 0-9]+\] \.hash .*
 +\[[ 0-9]+\] \.dynsym .*
 +\[[ 0-9]+\] \.dynstr .*
 +\[[ 0-9]+\] \.rela\.dyn .*
 +\[[ 0-9]+\] \.rela\.plt .*
 +\[[ 0-9]+\] \.text .*
 +\[[ 0-9]+\] \.tdata +PROGBITS .* 0+38 0+ WAT +0 +0 +8
 +\[[ 0-9]+\] \.tbss +NOBITS .* 0+38 0+ WAT +0 +0 +8
 +\[[ 0-9]+\] \.dynamic .*
 +\[[ 0-9]+\] \.opd .*
 +\[[ 0-9]+\] \.got .*
 +\[[ 0-9]+\] \.plt .*
 +\[[ 0-9]+\] \.symtab .*
 +\[[ 0-9]+\] \.strtab .*
 +\[[ 0-9]+\] \.shstrtab .*
#...

Elf file type is DYN \(Shared object file\)
Entry point 0x[0-9a-f]+
There are [0-9]+ program headers, starting at offset [0-9]+

Program Headers:
 +Type +Offset +VirtAddr +PhysAddr +FileSiz +MemSiz +Flg Align
 +LOAD .* R E 0x10000
 +LOAD .* RW +0x10000
 +DYNAMIC .* RW +0x8
 +TLS .* 0x0+38 0x0+70 R +0x8

 Section to Segment mapping:
 +Segment Sections\.\.\.
 +0+ +\.hash \.dynsym \.dynstr \.rela\.dyn \.rela\.plt \.text *
 +01 +\.tdata \.dynamic \.opd \.got \.plt *
 +02 +\.dynamic *
 +03 +\.tdata \.tbss *

Relocation section '\.rela\.dyn' at offset .* contains 20 entries:
 +Offset +Info +Type +Symbol's Value +Symbol's Name \+ Addend
[0-9a-f ]+R_PPC64_RELATIVE +51c
[0-9a-f ]+R_PPC64_RELATIVE +18800
[0-9a-f ]+R_PPC64_TPREL16 +0+60 le0 \+ 0
[0-9a-f ]+R_PPC64_TPREL16_HA +0+68 le1 \+ 0
[0-9a-f ]+R_PPC64_TPREL16_LO +0+68 le1 \+ 0
[0-9a-f ]+R_PPC64_TPREL16_DS +28
[0-9a-f ]+R_PPC64_TPREL16_HA +30
[0-9a-f ]+R_PPC64_TPREL16_LO +30
[0-9a-f ]+R_PPC64_DTPMOD64 +0
[0-9a-f ]+R_PPC64_DTPREL64 +0
[0-9a-f ]+R_PPC64_DTPREL64 +18
[0-9a-f ]+R_PPC64_DTPMOD64 +0
[0-9a-f ]+R_PPC64_DTPMOD64 +0+ gd \+ 0
[0-9a-f ]+R_PPC64_DTPREL64 +0+ gd \+ 0
[0-9a-f ]+R_PPC64_DTPMOD64 +0+40 ld0 \+ 0
[0-9a-f ]+R_PPC64_DTPMOD64 +0+ ld \+ 0
[0-9a-f ]+R_PPC64_DTPREL64 +0+50 ld2 \+ 0
[0-9a-f ]+R_PPC64_DTPMOD64 +0+38 gd0 \+ 0
[0-9a-f ]+R_PPC64_DTPREL64 +0+38 gd0 \+ 0
[0-9a-f ]+R_PPC64_TPREL64 +0+58 ie0 \+ 0

Relocation section '\.rela\.plt' at offset .* contains 1 entry:
 +Offset +Info +Type +Symbol's Value +Symbol's Name \+ Addend
[0-9a-f ]+R_PPC64_JMP_SLOT +0+ __tls_get_addr \+ 0

Symbol table '\.dynsym' contains [0-9]+ entries:
 +Num: +Value +Size +Type +Bind +Vis +Ndx +Name
.* NOTYPE +LOCAL +DEFAULT +UND *
.* SECTION +LOCAL +DEFAULT +6.*
.* SECTION +LOCAL +DEFAULT +10.*
.* TLS +GLOBAL +DEFAULT +UND gd
.* TLS +GLOBAL +DEFAULT +8 le0
.* NOTYPE +GLOBAL +DEFAULT +UND __tls_get_addr
.* TLS +GLOBAL +DEFAULT +8 ld0
.* TLS +GLOBAL +DEFAULT +8 le1
.* TLS +GLOBAL +DEFAULT +UND ld
.* FUNC +GLOBAL +DEFAULT +10 _start
.* TLS +GLOBAL +DEFAULT +8 ld2
.* TLS +GLOBAL +DEFAULT +8 ld1
.* TLS +GLOBAL +DEFAULT +8 gd0
.* TLS +GLOBAL +DEFAULT +8 ie0

Symbol table '\.symtab' contains [0-9]+ entries:
 +Num: +Value +Size +Type +Bind +Vis +Ndx +Name
.* NOTYPE +LOCAL +DEFAULT +UND *
.* FILE +LOCAL +DEFAULT +ABS .*
.* NOTYPE +LOCAL +DEFAULT +ABS TLSMARK
.* TLS +LOCAL +DEFAULT +7 gd4
.* TLS +LOCAL +DEFAULT +7 ld4
.* TLS +LOCAL +DEFAULT +7 ld5
.* TLS +LOCAL +DEFAULT +7 ld6
.* TLS +LOCAL +DEFAULT +7 ie4
.* TLS +LOCAL +DEFAULT +7 le4
.* TLS +LOCAL +DEFAULT +7 le5
.* FILE +LOCAL +DEFAULT +ABS
.* OBJECT +LOCAL +DEFAULT +9 _DYNAMIC
.* NOTYPE +LOCAL +DEFAULT +6 __glink_PLTresolve
.* NOTYPE +LOCAL +DEFAULT +6 .*\.plt_call\.__tls_get_addr
.* TLS +GLOBAL +DEFAULT +UND gd
.* TLS +GLOBAL +DEFAULT +8 le0
.* NOTYPE +GLOBAL +DEFAULT +UND __tls_get_addr
.* TLS +GLOBAL +DEFAULT +8 ld0
.* TLS +GLOBAL +DEFAULT +8 le1
.* TLS +GLOBAL +DEFAULT +UND ld
.* FUNC +GLOBAL +DEFAULT +10 _start
.* TLS +GLOBAL +DEFAULT +8 ld2
.* TLS +GLOBAL +DEFAULT +8 ld1
.* TLS +GLOBAL +DEFAULT +8 gd0
.* TLS +GLOBAL +DEFAULT +8 ie0


text:
###############################################################################
#
# plothascountsvalues.r - plot inertion and unique values bar graphs
#
# File:    plothascountsvalues.r
# Author:  Alex Stivala
# Created: January 2013
#
#
#
# Plot hashcount (insertion count) and unique value inserted (entries in
# hash table at end) frmo instrumentatino output summarized by
# summarize_instrument.sh
#
# Usage:
#       R --vanilla --slave -f plothascountsvalues.r
#
#  Reads the rtab files teslam2070.instrumentrtab,  etc.
#  Output is  PostScript file knapsack_hashandvalucounts.eps
#
# Requires the gplots package from CRAN to draw error bars.
#
# $Id: plothascountsvalues.r 2526 2009-06-15 06:23:41Z astivala $
#
###############################################################################

library(gplots)

#
# globals
#

colorvec=c('deepskyblue4','brown','red','turquoise','blue','purple','green','cyan','gray20','magenta','darkolivegreen2','midnightblue','magenta3','darkseagreen','violetred3','darkslategray3')
ltyvec=c(1,2,4,5,6,1,2,1,5,6,1,2,4,5,6,1,2)
pchvec=c(20,21,22,23,24,25)
# these must line up so the GPU output in gpu_fileprefixvec[i] corresponds

maxthreads <- 8 # number of CPU threads to use

filenamevec <- c('edwardgpunode.instrument.rtab',
                 'edwardgpunode_randomstart.instrument.rtab',
                 'teslam2070_norandomstart.instrument.rtab',
                 'teslam2070.instrument.rtab')
namevec=c('CPU', 'CPU random starts', 'GPU', 'GPU random starts')

count_types=c('total computations (h/h1)','subproblems computed (v/v1)')
basefilename <- filenamevec[1] # contains threads==1 as baseline

plotcountsandvaluesincrease <- function()
{
  epsfile <- 'knapsack_hashandvaluecounts.eps'

  # EPS suitable for inserting into LaTeX
  postscript(epsfile, onefile=FALSE,paper="special",horizontal=FALSE,
             width = 9, height = 6)

  basertab <- read.table(basefilename, header=T)
  basevaluecount <- mean(subset(basertab, threads==1)$hn)
  basehashcount <- mean(subset(basertab, threads==1)$hc)

  means <- c()
  stdev <- c()
  ciw <- c()
  for (j in 1:length(filenamevec)) {
    rtab <- read.table(filenamevec[j], header=TRUE)
    if (!is.na(pmatch('GPU', namevec[j])))
      isgpu <- T
    else
      isgpu <- F
    if (isgpu)  {
      means <- c(means,  mean(subset(rtab, is.na(threads))$hc/basehashcount) )
      means <- c(means,  mean(subset(rtab, is.na(threads))$hn/basevaluecount) )
    }
    else {
      means <- c(means,  mean(subset(rtab, threads == maxthreads)$hc/basehashcount) )
      means <- c(means,  mean(subset(rtab, threads == maxthreads)$hn/basevaluecount) )
    }
    this_hashstdev <-  sqrt(var(subset(rtab, is.na(threads))$hc/basehashcount))
    this_valuestdev <-  sqrt(var(subset(rtab, is.na(threads))$hn/basevaluecount))
    stdev <- c(stdev , this_hashstdev, this_valuestdev)
    n <- length(subset(rtab, is.na(threads))$hn)
    ciw <- c(ciw,  qt(0.975, n) * this_hashstdev / sqrt(n))
    ciw <- c(ciw,  qt(0.975, n) * this_valuestdev / sqrt(n))
   }
   mean_increase_matrix <- matrix(means, nrow = length(count_types), byrow=F)
   ci_matrix <- matrix(ciw,  nrow = length(count_types)   , byrow=F)

   ylim <- c(0,110)
   mp <- barplot2(mean_increase_matrix, #names.arg = namevec,
            beside = TRUE,
            ylab="computation increase ratio",  ylim=ylim,
            plot.ci = TRUE,
            ci.u = mean_increase_matrix + ci_matrix,
            ci.l = mean_increase_matrix - ci_matrix,
            col=colorvec[1:2],
            names.arg=namevec
            )
   text(mp, mean_increase_matrix+05, format(means, digits=2), xpd=TRUE)
   legend('topright', fill=colorvec, legend=count_types,bty='n')
   box()
   dev.off()
}

#
# main
#

plotcountsandvaluesincrease()


text:
#' Basic Data Import for Water Flow Data
#'
#' Imports data from user-supplied data file. Specifically used to import water flow data for use in the EGRET package.
#' For EGRET usage, the first column is expected to be dates. If the data is daily data, then next column is
#' expected to be the measured values. If the data is sampled data, the next column is remark codes, and the third
#' column is values.
#'
#' @param filePath character specifying the path to the file
#' @param fileName character name of file to open
#' @param hasHeader logical true if the first row of data is the column headers
#' @param separator character character that separates data cells
#' @keywords data import file
#' @return retval dataframe with dateTime, value, and code columns
#' @export
#' @examples
#' filePath <- system.file("extdata", package="EGRET")
#' fileName <- 'ChoptankRiverFlow.txt'
#' ChopData <- readDataFromFile(filePath,fileName, separator="\t")
readDataFromFile <- function (filePath,fileName,hasHeader=TRUE,separator=","){
  totalPath <- file.path(filePath,fileName)
  if(file.exists(totalPath)){
    retval <- read.delim(
      totalPath,
      header = hasHeader,
      sep=separator,
      colClasses=c('character'),
      fill = TRUE,
      comment.char="#")

    if(ncol(retval) == 2){
      numCol <- 2
    } else {
      numCol <- seq(from = 3,to = ncol(retval), by = 2)
    }

    if(dateFormatCheck(retval[,1])){
      retval[,1] <- as.Date(retval[,1])
    } else {
      retval[,1] <- as.Date(retval[,1],format="%m/%d/%Y")
    }

    retval[,numCol] <- sapply(retval[,numCol], as.numeric)

    return (retval)
  } else {

    stop("File does not exist")
  }
}

max_stars_repo_name: ldecicco-USGS/EGRET

max_issues_repo_name: ldecicco-USGS/EGRET

max_forks_repo_name: ldecicco-USGS/EGRET


text:
\name{circos.genomicLabels}
\alias{circos.genomicLabels}
\title{
Add labels to specified genomic regions
}
\description{
Add labels to specified genomic regions
}
\usage{
circos.genomicLabels(bed, labels = NULL, labels.column = NULL,
    facing = "clockwise", niceFacing = TRUE,
    col = par("col"), cex = 0.8, font = par("font"), padding = 0.4,
    connection_height = convert_height(5, "mm"),
    line_col = par("col"), line_lwd = par("lwd"), line_lty = par("lty"),
    labels_height = min(c(convert_height(1.5, "cm"),
    max(strwidth(labels, cex = cex, font = font)))),
    side = c("inside", "outside"), track.margin = circos.par("track.margin"))
}
\arguments{

  \item{bed}{a data frame in bed format
}
  \item{labels}{a vector of labels corresponding to rows in \code{bed}
}
  \item{labels.column}{if the label column is already in \code{bed}, the index for this column in \code{bed}
}
  \item{facing}{facing of the labels. The value can only be 'clockwise' or 'reverse.clockwise'.
}
  \item{niceFacing}{whether automatically adjust the facing of the labels.
}
  \item{col}{color for the labels
}
  \item{cex}{size of the labels
}
  \item{font}{font of the labels
}
  \item{padding}{padding of the labels, the value is the ratio to the height of the label
}
  \item{connection_height}{height of the connection track
}
  \item{line_col}{color for the connection lines
}
  \item{line_lwd}{line width for the connection lines
}
  \item{line_lty}{line type for the connectioin lines
}
  \item{labels_height}{height of the labels track
}
  \item{side}{side of the labels track, is it in the inside of the track where the regions are marked?
}
  \item{track.margin}{bottom and top margins
}

}
\details{
The function adds labels for the specified regions. The positions of labels are arranged
so that they are not overlapping to each other.
}
\author{
Zuguang Gu <z.gu@dkfz.de>
}
\examples{
\dontrun{
circos.initializeWithIdeogram(plotType = c("labels", "axis"))
bed = generateRandomBed(nr = 100, fun = function(k) sample(letters, k, replace = TRUE))
bed[1, 4] = "aaaaaaaa"
circos.genomicLabels(bed, labels.column = 4, side = "inside",
    col = as.numeric(factor(bed[[1]])))
circos.genomicLabels(bed, labels.column = 4, side = "outside",
    line_col = as.numeric(factor(bed[[1]])))
}
}

ext: rd

max_stars_repo_name: calpan/circlize

max_issues_repo_name: Nexller/circlize

max_forks_repo_name: Nexller/circlize