text:
# Looker API 3.0 Reference
#
# ### Authorization  The Looker API uses Looker **API3** credentials for authorization and access control. Looker admins can create API3 credentials on Looker's **Admin/Users** page. Pass API3 credentials to the **/login** endpoint to obtain a temporary access_token. Include that access_token in the Authorization header of Looker API requests. For details, see [Looker API Authorization](https://looker.com/docs/r/api/authorization)  ### Client SDKs  The Looker API is a RESTful system that should be usable by any programming language capable of making HTTPS requests. Client SDKs for a variety of programming languages can be generated from the Looker API's Swagger JSON metadata to streamline use of the Looker API in your applications. A client SDK for Ruby is available as an example. For more information, see [Looker API Client SDKs](https://looker.com/docs/r/api/client_sdks)  ### Try It Out!  The 'api-docs' page served by the Looker instance includes 'Try It Out!' buttons for each API method. After logging in with API3 credentials, you can use the \"Try It Out!\" buttons to call the API directly from the documentation page to interactively explore API features and responses.  ### Versioning  Future releases of Looker will expand this API release-by-release to securely expose more and more of the core power of Looker to API client applications. API endpoints marked as \"beta\" may receive breaking changes without warning. Stable (non-beta) API endpoints should not receive breaking changes in future releases. For more information, see [Looker API Versioning](https://looker.com/docs/r/api/versioning)
#
# OpenAPI spec version: 3.0.0
#
# Generated by: https://github.com/swagger-api/swagger-codegen.git

#' @title Project operations
#' @description looker.Project
#'
#' @field path Stores url path of the request.
#' @field apiClient Handles the client-server communication.
#' @field userAgent Set the user agent of the request.
#'
#' @importFrom R6 R6Class
#'
#' @section Methods:
#' \describe{
#'
#' all_git_branches Get All Git Branchs
#'
#'
#' all_git_connection_tests Get All Git Connection Tests
#'
#'
#' all_project_files Get All Project Files
#'
#'
#' all_projects Get All Projects
#'
#'
#' create_git_deploy_key Create Deploy Key
#'
#'
#' create_project Create Project
#'
#'
#' git_deploy_key Git Deploy Key
#'
#'
#' project Get Project
#'
#'
#' project_file Get Project File
#'
#'
#' project_validation_results Cached Project Validation Results
#'
#'
#' project_workspace Get Project Workspace
#'
#'
#' reset_project_to_production Reset To Production
#'
#'
#' reset_project_to_remote Reset To Remote
#'
#'
#' run_git_connection_test Run Git Connection Test
#'
#'
#' update_project Update Project
#'
#'
#' validate_project Validate Project
#'
#' }
#'
#' @export
ProjectApi <- R6::R6Class(
  'ProjectApi',
  public = list(
    userAgent = "Swagger-Codegen/1.0.0/r",
    apiClient = NULL,
    initialize = function(apiClient){
      if (!missing(apiClient)) {
        self$apiClient <- apiClient
      }
      else {
        self$apiClient <- ApiClient$new()
      }
    },
    all_git_branches = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/git_branches"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- GitBranch$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    all_git_connection_tests = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/git_connection_tests"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- GitConnectionTest$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    all_project_files = function(project_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}/files"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- ProjectFile$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    all_projects = function(fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects"
      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Project$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    create_git_deploy_key = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/git/deploy_key"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "POST",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Character$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    create_project = function(body, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`body`)) {
        body <- `body`$toJSONString()
      } else {
        body <- NULL
      }

      urlPath <- "/projects"
      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "POST",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Project$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    git_deploy_key = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/git/deploy_key"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Character$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    project = function(project_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Project$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    project_file = function(project_id, file_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`file_id`)) {
        queryParams['file_id'] <- file_id
      }

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}/files/file"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- ProjectFile$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    project_validation_results = function(project_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}/validate"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- ProjectValidationCache$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    project_workspace = function(project_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}/current_workspace"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- ProjectWorkspace$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    reset_project_to_production = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/reset_to_production"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "POST",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Character$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    reset_project_to_remote = function(project_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/reset_to_remote"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "POST",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Character$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    run_git_connection_test = function(project_id, test_id, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      urlPath <- "/projects/{project_id}/git_connection_tests/{test_id}"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      if (!missing(`test_id`)) {
        urlPath <- gsub(paste0("\\{", "test_id", "\\}"), `test_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "GET",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- GitConnectionTestResult$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    update_project = function(project_id, body, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      if (!missing(`body`)) {
        body <- `body`$toJSONString()
      } else {
        body <- NULL
      }

      urlPath <- "/projects/{project_id}"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "PATCH",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- Project$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    },
    validate_project = function(project_id, fields, ...){
      args <- list(...)
      queryParams <- list()
      headerParams <- character()

      if (!missing(`fields`)) {
        queryParams['fields'] <- fields
      }

      urlPath <- "/projects/{project_id}/validate"
      if (!missing(`project_id`)) {
        urlPath <- gsub(paste0("\\{", "project_id", "\\}"), `project_id`, urlPath)
      }

      resp <- self$apiClient$callApi(url = paste0(self$apiClient$basePath, urlPath),
                                 method = "POST",
                                 queryParams = queryParams,
                                 headerParams = headerParams,
                                 body = body,
                                 ...)

      if (httr::status_code(resp) >= 200 && httr::status_code(resp) <= 299) {
        returnObject <- ProjectValidation$new()
        result <- returnObject$fromJSON(httr::content(resp, "text", encoding = "UTF-8"))
        Response$new(returnObject, resp)
      } else if (httr::status_code(resp) >= 400 && httr::status_code(resp) <= 499) {
        Response$new("API client error", resp)
      } else if (httr::status_code(resp) >= 500 && httr::status_code(resp) <= 599) {
        Response$new("API server error", resp)
      }

    }
  )
)

max_stars_repo_name: adamlouis/lookr

max_issues_repo_name: adamlouis/lookr

max_forks_repo_name: adamlouis/lookr


text: step(aov(log(Days+2.5) ~ .^4, data=Quine))

ext: r

max_stars_repo_name: DeagleGross/Rplugin

max_issues_repo_name: DeagleGross/Rplugin

max_forks_repo_name: DeagleGross/Rplugin


text:
#' @export
set_plot_options <- function(...) {
    .Deprecated('options', msg = 'use the `repr.plot.*` options from the repr package instead')
    opts <- list(...)
    names(opts) <- paste0('repr.plot.', names(opts))
    do.call(options, opts)
}

#' @export
get_plot_options <- function() {
    .Deprecated('getOption', msg = 'use the `repr.plot.*` options from the repr package instead')
    all.opts <- options()
    plot.opt.idx <- grep('^repr.plot', names(all.opts))
    plot.opts <- all.opts[plot.opt.idx]
    names(plot.opts) <- gsub('^repr\\.plot\\.', '', names(plot.opts))
    return(as.list(plot.opts))
}

opt.defaults <- list(
    jupyter.plot_mimetypes = c(
        'text/plain',
        'image/png',
        'image/svg+xml'))

.onLoad <- function(libname = NULL, pkgname = NULL) {
    for (option in names(opt.defaults)) {
        if (is.null(getOption(option))) {
            do.call(options, opt.defaults[option])  # single []: preserve name
        }
    }
}

max_stars_repo_name: ChinaQuants/IRkernel

max_issues_repo_name: ChinaQuants/IRkernel

max_forks_repo_name: ChinaQuants/IRkernel


text:
library(affy)

basedir<- "/Users/Tammy/rawdata-master/"
file.path(basedir, "celfiles/sampleinfo.txt")

sample<- read.delim(file.path(basedir, "celfiles/sampleinfo.txt"), header=TRUE)

sample[sample$filenames=="1521a99hpp_av06.CEL.gz", "X36311_at"]

basedir<- "/Users/Tammy/rawdata-master/celfiles"
fns<- list.celfiles(basedir)
fns %in% sample[,1]
ab<- ReadAffy(filenames= file.path(basedir, sample[,1]), phenoData=sample)

table(probeNames(ab)=="36311_at")
## there are 16 probes for the same feature 36311_at

## gene expression for all the probes
pm(ab)[1:10,1:6]

## note it is different from exprs(ab) which summarise the gene expression at the feature level

which(probeNames(ab)=="36085_at")

## extract all the probes for this feature 36085_at
par(mfrow=c(1,1))
mat<- pm(ab)[which(probeNames(ab)=="36085_at"),]
matplot(mat, type="l")
matplot(log2(mat[,"1532a99hpp_av04.CEL.gz"]), type="l")
matplot(log2(mat[,"1532b99hpp_av04.CEL.gz"]), type="l")
matplot(mat[,"1532a99hpp_av04.CEL.gz"], mat[,"1532b99hpp_av04.CEL.gz"], type="l")

eset<- rma(ab)
g<-  factor(pData(eset)[,2])

library(genefilter)
rowttests(eset, g)["36085_at",]

### spike-ins
sig<- colnames(pData(ab))[-1]

##get rid of the leading X for each feature
sig

rowttests(eset,g)[sig,]
library(limma)
datadir = "/Users/Tammy/rawdata-master/"
basedir = file.path(datadir, "/agilent")
setwd(basedir)
targets = readTargets("TargetBeta7.txt")
RG = read.maimages(targets$FileName, source="genepix")
MA<- MA.RG(RG, bc.method="none")

###################### NGS data
biocLite("pasillaBamSubset")
library(pasillaBamSubset)
library(Rsamtools)

# the path to the bam file
filename <- untreated1_chr4()

# add parathensis will evaluate the expression
(bf <- BamFile(filename))

seqinfo(bf)
(sl <- seqlengths(bf))

# summary of bam file
quickBamFlagSummary(bf)

?scanBamParam

## genomic ranges for chr4
(gr <- GRanges("chr4",IRanges(1, sl["chr4"])))

countBam(bf, param=ScanBamParam(which = gr))

?ScanBamParam
ScanBamParam(which = gr)

reads<- scanBam(BamFile(filename, yieldSize=5))
?scanBam

class(reads)
names(reads[[1]])

reads[[1]]$pos
reads[[1]]$rname
reads[[1]]$strand
reads[[1]]$seq

library(GenomicAlignments)
(ga <- readGAlignments(bf))

length(ga)
granges(ga[1])

gr <- GRanges("chr4", IRanges(700000, 800000))
(fo <- findOverlaps(ga, gr)) # which reads over this range

countOverlaps(gr, ga) # count overlaps of range with the reads

table(ga %over% gr) # logical vector of read overlaps with the range

## If we had run countOverlaps(ga, gr) it would return an integer vector
## with the number of overlaps for each read with the range in gr.
countOverlaps(ga, gr)

gr <- GRanges("chr4", IRanges(440000, 470000))
countOverlaps(gr, ga)

reads <- scanBam(bf, param=ScanBamParam(what=c("pos","strand","seq"), which=gr))
## access sequences of the reads by
reads$`chr4:440000-470000`$seq

## GC content of the sequences
apply(letterFrequency(x=reads$`chr4:440000-470000`$seq, letters=c("C","G"), as.prob=T), 1, sum)

mean (apply(letterFrequency(x=reads$`chr4:440000-470000`$seq, letters=c("C","G"), as.prob=T), 1, sum))

biocLite("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
g<- genes(TxDb.Dmelanogaster.UCSC.dm3.ensGene)

## Let's focus on the genes at a particular location:

g2<- g[g %over% GRanges("chr4",IRanges(200000, 300000))]

g2

g2['FBgn0039890']
ranges(g2['FBgn0039890'])

## how many reads mapped on to this gene?
countOverlaps(g2['FBgn0039890'], ga, ignore.strand=FALSE)

##### count reads from a bam file

library(pasillaBamSubset)
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
library(GenomicAlignments)

## genes on chr4
g<- genes(txdb)
g<- g[seqnames(g) == "chr4"]

txdb<- TxDb.Dmelanogaster.UCSC.dm3.ensGene
grl<- exonsBy(txdb, by="gene")

#GRangesList for chr4
grl<- grl[names(g)]

# the 10th gene
grl[10]

# access the GRanges
grl[[10]]

#first exon
grl[[10]][1]

# Test for the same names:
all.equal(names(g), names(grl))

bf<- BamFile(untreated1_chr4())
so1<- summarizeOverlaps(features=grl,
                        reads= bf,
                        ignore.strand=TRUE)
head(assay(so1))
colSums(assay(so1))
rowData(so1)
colData(so1)
colData(so1) <- c("one", "two")
metadata(rowData(so1)) # rowData depricated
metadata(rowRanges(so1))

so2<- summarizeOverlaps(features=g,
                        reads= bf,
                        ignore.strand=TRUE)

head(assay(so2))

plot(log2(assay(so2)[,1]), log2(assay(so1)[,1]))
abline(a=0, b=1)

## some reads mapped to exon-exon junctions, we see more reads fro g compared with grl
mean( assay(so1)[assay(so2)[,1] >0,] / assay(so2)[assay(so2)[,1] >0,] )

head(assay(so1))
## FPM (fragments per million counts for each gene (exons))
head( assay(so1)[,1]/colSums(assay(so1)) * 10^6 )

### FPKM

## reduce can directly work on GRangelist,
## try reduce on Grange object as well.

## look at the first gene
reduce(grl[[1]])
width(reduce(grl[[1]]))
width(grl[[1]])

## compare with this
width(reduce(grl))

## sum of exons for each gene
sapply(width(reduce(grl)), sum)
exons_length_per_gene<-sapply(width(reduce(grl)), sum)

##
summary(sapply(width(reduce(grl)), sum))

### calculate RPKM
FPKM<- (assay(so1)[,1]/colSums(assay(so1)) * 10^6)/exons_length_per_gene * 1000
FPKM["FBgn0002521"]

#### visualize
# coverage of the GenomicAlignment object ga
fl1<- untreated1_chr4()
fl2<- untreated3_chr4()

x<- BamFile(fl1)
y<- BamFile(fl2)

xcov<- coverage(x)
xcov
z<- GRanges("chr4", IRanges(456500,466000))
# subset the coverage using z
xcov[z]  # it is a Rlelist
## or
xcov$chr4[ranges(z)]  # it is integer-Rle

xnum<- as.numeric(unlist(xcov[z]))
plot(xnum)

ycov<- coverage(y)
ynum<- as.numeric(ycov$chr4[ranges(z)])
plot(xnum, type="l", col="blue", lwd=2)
lines(ynum, col="red", lwd=2)

## zoom in
plot(xnum, type="l", col="blue", lwd=2, xlim=c(6200,6600))
lines(ynum, col="red", lwd=2)

### check gene lgs using biomart
library(biomaRt)
listMarts()
m<- useMart("ensembl", dataset="dmelanogaster_gene_ensembl")
lf<- listFilters(m)
lf[grep("name", lf$description, ignore.case=TRUE),]

map<- getBM(mart=m,
            attributes=c("ensembl_gene_id", "flybasename_gene"),
            filters="flybasename_gene",
            values="lgs")

map

grl<- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, by="gene")
gene<- grl[[map$ensembl_gene_id[1]]]
gene

### visualize by Gviz
library(Gviz)
gtrack<- GenomeAxisTrack()
atrack<- AnnotationTrack(gene, name="Gene Model")
plotTracks(list(gtrack, atrack))

xgr<- as(xcov, "GRanges")
ygr<- as(ycov, "GRanges")
dtrack1<- DataTrack(xgr[xgr %over% z], name="sample 1")
dtrack2<- DataTrack(ygr[ygr %over% z], name="sample 2")

plotTracks(list(gtrack, atrack, dtrack1, dtrack2))
plotTracks(list(gtrack, atrack, dtrack1, dtrack2), type="polygon")

## use ggbio  built on ggplot2
library(ggbio)
autoplot(gene)
autoplot(fl1, which=z)
autoplot(fl2, which=z)

max_stars_repo_name: crazyhottommy/scripts-general-use

max_issues_repo_name: crazyhottommy/scripts-general-use

max_forks_repo_name: crazyhottommy/scripts-general-use


text:
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

#' Print MVP Banner
#'
#' Build date: Aug 30, 2017
#' Last update: Dec 12, 2018
#'
#' @author Lilin Yin, Haohao Zhang, and Xiaolei Liu
#'
#' @param width the width of the message
#' @param verbose whether to print detail.
#'
#' @return version number.
#'
#' @examples
#' MVP.Version()
MVP.Version <- function(width=60, verbose=TRUE) {
    welcome <- "Welcome to MVP"
    title   <- "A Memory-efficient, Visualization-enhanced, and Parallel-accelerated Tool For GWAS"
    authors <- c("Designed and Maintained by Lilin Yin, Haohao Zhang, and Xiaolei Liu",
                 "Contributors: Zhenshuang Tang, Jingya Xu, Dong Yin, Zhiwu Zhang, Xiaohui Yuan, Mengjin Zhu, Shuhong Zhao, Xinyun Li")
    contact <- "Contact: xiaoleiliu@mail.hzau.edu.cn"
    logo_s  <- c(" __  __  __   __  ___",
                 "|  \\/  | \\ \\ / / | _ \\",
                 "| |\\/| |  \\ V /  |  _/",
                 "|_|  |_|   \\_/   |_|")

    version <- print_info(welcome = welcome, title = title, logo = logo_s, authors = authors, contact = contact, linechar = '=', width = width, verbose = verbose)
    return(invisible(version))
}

#' Print progress bar
#'
#' @param i the current loop number
#' @param n max loop number
#' @param type type1 for "for" function
#' @param symbol the symbol for the rate of progress
#' @param tmp.file the opened file of "fifo" function
#' @param symbol.head the head for the bar
#' @param symbol.tail the tail for the bar
#' @param fixed.points whether use the setted points which will be printed
#' @param points the setted points which will be printed
#' @param symbol.len the total length of progress bar
#'
#' @keywords internal
print_bar <- function(i,
                    n,
                    type = c("type1", "type3"),
                    symbol = "-",
                    tmp.file = NULL,
                    symbol.head = ">>>",
                    symbol.tail = ">" ,
                    fixed.points = TRUE,
                    points = seq(0, 100, 1),
                    symbol.len = 48,
                    verbose = TRUE
) {
    switch(
        match.arg(type),
        "type1"={
            if(fixed.points){
                point.index <- points
                point.index <- point.index[point.index > floor(100*(i-1)/n)]
                if(floor(100*i/n) %in% point.index){
                    if(floor(100*i/n) != max(point.index)){
                        print.len <- floor(symbol.len*i/n)
                        logging.log(
                            paste("\r",
                                  paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                                  paste(rep(" ", symbol.len-print.len), collapse=""),
                                  sprintf("%.2f%%", 100*i/n)
                                  , sep=""),
                            verbose = verbose
                        )
                    }else{
                        print.len <- floor(symbol.len*i/n)
                        logging.log(
                            paste("\r",
                                  paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                                  sprintf("%.2f%%", 100*i/n), "\n"
                                  , sep=""),
                            verbose = verbose
                        )
                    }
                }
            }else{
                if(i < n){
                    print.len <- floor(symbol.len*i/n)
                    logging.log(
                        paste("\r",
                              paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                              paste(rep(" ", symbol.len-print.len), collapse=""),
                              sprintf("%.2f%%", 100*i/n)
                              , sep=""),
                        verbose = verbose
                    )
                }else{
                    print.len <- floor(symbol.len*i/n)
                    logging.log(
                        paste("\r",
                              paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                              sprintf("%.2f%%", 100*i/n), "\n"
                              , sep=""),
                        verbose = verbose
                    )
                }
            }
        },
        # "type2"={
        #     if(inherits(parallel:::mcfork(), "masterProcess")) {
        #         progress <- 0.0
        #         while(progress < n && !isIncomplete(tmp.file)){
        #             msg <- readBin(tmp.file, "double")
        #             progress <- progress + as.numeric(msg)
        #             print.len <- round(symbol.len * progress / n)
        #             if(fixed.points){
        #                 if(progress %in% round(points * n / 100)){
        #                     logging.log(paste("\r",
        #                               paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
        #                               paste(rep(" ", symbol.len-print.len), collapse=""),
        #                               sprintf("%.2f%%", progress * 100 / n), sep=""))
        #                 }
        #             }else{
        #                 logging.log(paste("\r",
        #                           paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
        #                           paste(rep(" ", symbol.len-print.len), collapse=""),
        #                           sprintf("%.2f%%", progress * 100 / n), sep=""))
        #             }
        #         }
        #         parallel:::mcexit()
        #     }
        # },
        "type3"={
            progress <- readBin(tmp.file, "double") + 1
            writeBin(progress, tmp.file)
            print.len <- round(symbol.len * progress / n)
            if(fixed.points){
                if(progress %in% round(points * n / 100)){
                    logging.log(
                        paste("\r",
                              paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                              paste(rep(" ", symbol.len-print.len), collapse=""),
                              sprintf("%.2f%%", progress * 100 / n)
                              , sep=""),
                        verbose = verbose
                    )
                }
            }else{
                logging.log(
                    paste("\r",
                          paste(c(symbol.head, rep("-", print.len), symbol.tail), collapse=""),
                          paste(rep(" ", symbol.len-print.len), collapse=""),
                          sprintf("%.2f%%", progress * 100 / n)
                          , sep=""),
                    verbose = verbose
                )
            }
        }
    )
}

print_accomplished <- function(width = 60, verbose = TRUE) {
    logging.log(make_line("MVP ACCOMPLISHED", width = width, linechar = '='), "\n", verbose = verbose)
}

#' Print R Package information, include title, short_title, logo, version, authors, contact
#'
#' Build date: Oct 22, 2018
#' Last update: Oct 22, 2018
#'
#' @keywords internal
#' @author Haohao Zhang
#'
#' @param welcome welcome text, for example: "Welcom to <Packagename>"
#' @param title long text to introduct package
#' @param short_title short label, top-left of logo
#' @param logo logo
#' @param version short label, bottom-right of logo
#' @param authors
#' @param contact email or website
#' @param line 1, 2, or char
#' @param width banner width
#'
#' @examples
#' welcome <- "Welcome to MVP"
#' title   <- "A Memory-efficient, Visualization-enhanced, and Parallel-accelerated Tool For GWAS"
#' authors <- "Authors: Lilin Yin, Haohao Zhang, and Xiaolei Liu"
#' contact <- "Contact: xiaoleiliu@mail.hzau.edu.cn"
#' logo_s  <- c(" __  __  __   __  ___",
#'              "|  \/  | \ \ / / | _ \",
#'              "| |\/| |  \ V /  |  _/",
#'              "|_|  |_|   \_/   |_|")
#' print_info(welcome = welcome, title = title, logo = logo_s, authors = authors, contact = contact, linechar = '=', width = width)
print_info <- function(welcome=NULL, title=NULL, short_title=NULL, logo=NULL, version=NULL, authors=NULL, contact=NULL, linechar = '=', width=NULL, verbose=TRUE) {
    msg <- c()
    # width
    if (is.null(width)) { width <- getOption('width') }
    # version
    if (is.null(version)) {
        if (getPackageName() == ".GlobalEnv") {
            version <- "devel"
        } else {
            version <- as.character(packageVersion(getPackageName()))
        }
    }
    # welcome
    if (is.null(welcome)) {
        if (getPackageName() == ".GlobalEnv") {
            welcome <- ""
        } else {
            welcome <- paste0("Welcome to ", getPackageName())
        }
    }
    msg <- c(msg, make_line(welcome, linechar = linechar, width = width))
    # title
    if (!is.null(title)) {
        msg <- c(msg, rule_wrap(string = title, width = width, align = "center"))
    }

    # align logo
    logo_width <- max(sapply(logo, nchar))
    for (i in 1:length(logo)) {
        l <- paste0(logo[i], paste(rep(" ", logo_width - nchar(logo[i])), collapse = ""))
        l <- make_line(l, width)
        msg <- c(msg, l)
    }

    # paste short_title label to logo top-left
    if (!is.null(short_title)) {
        i <- length(msg) - length(logo) + 1
        msg[i] <- paste_label(msg[i], paste0(short_title), side = "left")
    }

    # paste version label to logo bottom-right
    msg[length(msg)] <- paste_label(msg[length(msg)], paste0("Version: ", version), side = "right")

    # authors
    if (!is.null(authors)) {
        msg <- c(msg, rule_wrap(string = authors, align = "left", linechar = " ", width = width))
    }
    # contact
    if (!is.null(contact)) {
        msg <- c(msg, rule_wrap(string = contact, align = "left", linechar = " ", width = width))
    }
    # bottom line
    msg <- c(msg, paste0(rep(linechar, width), collapse = ''))

    logging.log(msg, sep = "\n", verbose = verbose)

    return(version)
}

#' make line
#'
#' Build date: Dec 12, 2018
#' Last update: Dec 12, 2018
#'
#' @keywords internal
#' @author Haohao Zhang
make_line <- function(string, width, linechar = " ", align = "center", margin = 1) {
    string <- paste0(paste0(rep(" ", margin), collapse = ""),
                     string,
                     paste0(rep(" ", margin), collapse = ""))

    if (align == "center") {
        if (width > nchar(string)) {
            left_width <- (width - nchar(string)) %/% 2
            right_width <- width - nchar(string) - left_width
            string <-
                paste0(paste0(rep(linechar, left_width), collapse = ""),
                       string,
                       paste0(rep(linechar, right_width), collapse = ""))
        }
    } else if (align == "left") {
        if (width > nchar(string)) {
            string <-
                paste0(linechar,
                       string,
                       paste0(rep(linechar, width - nchar(string) - 1), collapse = ""))
        }
    }
    return(string)
}

#' wrap text to multiple line, align left, right or center.
#'
#' Build date: Oct 22, 2018
#' Last update: Dec 12, 2018
#' by using base::strwarp.
#'
#' @keywords internal
#' @author Haohao Zhang
rule_wrap <- function(string, width, align = "center", linechar = " ") {
    # define
    msg <- c()
    lines <- strwrap(string, width = width - 4)

    # wrap
    for (i in 1:length(lines)) {
        l <- make_line(lines[i], width = width, linechar = linechar, align = align)
        msg <- c(msg, l)
    }
    return(msg)
}

#' Paste label to a line
#'
#' Build date: Oct 22, 2018
#' Last update: Oct 22, 2018
#'
#' @param line long text
#' @param label short label
#' @param side "right" or "left"
#' @param margin default 2
#'
#' @keywords internal
#' @author Haohao Zhang
paste_label <- function(line, label, side = "right", margin = 2) {
    if (side == "right") {
        end   <- nchar(line) - margin
        start <- end - (nchar(label) - 1)
    } else {
        start <- 1 + margin
        end   <- start + (nchar(label) - 1)
    }
    substr(line, start, end) <- label
    return(line)
}

#' format time
#'
#' @param x seconds
#'
#' @keywords internal
#' @author Lilin Yin
format_time <- function(x) {
    h <- x %/% 3600
    m <- (x %% 3600) %/% 60
    s <- ((x %% 3600) %% 60)
    index <- which(c(h, m, s) != 0)
    num <- c(h, m, s)[index]
    num <- round(num, 0)
    char <- c("h", "m", "s")[index]
    return(paste0(num, char, collapse = ""))
}

load_if_installed <- function(package) {
    if (!identical(system.file(package = package), "")) {
        do.call('library', list(package))
        return(TRUE)
    } else {
        return(FALSE)
    }
}

mkl_env <- function(exprs, threads = 1) {
    if (load_if_installed("RevoUtilsMath")) {
        math.cores <- eval(parse(text = "getMKLthreads()"))
        eval(parse(text = "setMKLthreads(threads)"))
    }
    result <- exprs
    if (load_if_installed("RevoUtilsMath")) {
        eval(parse(text = "setMKLthreads(math.cores)"))
    }
    return(result)
}

remove_bigmatrix <- function(x, desc_suffix=".geno.desc", bin_suffix=".geno.bin") {
    name <- basename(x)
    path <- dirname(x)

    descfile <- paste0(x, desc_suffix)
    binfile  <- paste0(x, bin_suffix)

    remove_var <- function(binfile, envir) {
        for (v in ls(envir = envir)) {
            if (class(get(v, envir = envir)) == "big.matrix") {
                desc <- describe(get(v, envir = envir))@description
                if (desc$filename == binfile) {
                    rm(list = v, envir = envir)
                    gc()
                }
            }
        }
    }

    # remove_var(binfile, globalenv())
    remove_var(binfile, as.environment(-1L))

    if (file.exists(descfile)) {
        file.remove(descfile)
    }
    if (file.exists(binfile)) {
        file.remove(binfile)
    }
}

max_stars_repo_name: xiaolei-lab/rMVP

max_issues_repo_name: hxxonly/rMVP

max_forks_repo_name: hxxonly/rMVP


text:
#  James Rekow

disconnectedSubgraphs = function(g, v, numSubgraphs, numSample = 100, maxIters = 100){

  #  ARGS: g - input graph
  #        v - number of vertices in the output graphs
  #
  #  RETURNS: vList - a list of numSample disconnected subgraphs of g with v vertices, or NULL
  #                   if no such subgraphs exist. That is, each element of v is a vector of
  #                   vertices corresponding to a subgraph whose adjacency matrix is the
  #                   v x v zero matrix

  #  create adjacency matrix for connectivity graph
  adjMat = as.matrix(get.adjacency(g, type = "both"))

  #  get number of vertices in g
  numVertices = vcount(g)

  pairMat = combn(numVertices, 2)

  produceSubgraphs = function(x = NULL){

    #  create a list of all pairs of vertices
    vList = lapply(as.list(1:choose(numVertices, 2)), function(i) pairMat[ , i])

    isDisconnected = function(vVec){

      #  ARGS: vVec - vector of vertices contained in subgraph, where it is known that vVec without
      #               it's last element forms a completely disconnected subgraph
      #
      #  RETURNS: TRUE/FALSE depending on whether or not the vertices in vVec form a completely
      #           disconnected subgraph

      #  number of vertices in vVec
      numV = length(vVec)

      #  create vector of vertices which are already known to form a completely disconnected subgraph
      disconnectedV = vVec[-numV]

      #  check whether there are any edges connecting the newest vertex to the previous disconnected
      #  vertices
      disconnected = !any(adjMat[ , vVec[numV]][disconnectedV] == 1)

      return(disconnected)

    } #  end isDisconnected function

    #  remove all vertex vectors which correspond to subgraphs containing at least one edge from vList
    vList = vList[unlist(lapply(vList, isDisconnected))]

    createNewVList = function(vList, numVertices){

      #  ARGS: vList - input vList
      #
      #  RETURNS: newVList - a list of all of the elements in vList with each element of 1:M not in
      #           a given vList appended to vList separately

      expandVVec = function(vVec, numVertices){

        #  create vector of potential vertices
        baseV = 1:numVertices

        #  select elements from baseV which are not already present in vVec
        #newV = baseV[-which(baseV %in% vVec)] - deprecated
        newV = baseV[-vVec]

        #  create a list of new vertex vectors, where each element is a vector that is vVec with one
        #  element of newV appended to it, and each element in newV corresponds to exactly one vector
        #  in this newly created list
        vAppendedList = lapply(as.list(newV), function(x) append(vVec, x))

        return(vAppendedList)

      } #  end expandVVec function

      #  create new list of vertex vectors
      newVList = unlist(lapply(vList, function(vVec) expandVVec(vVec = vVec,
                                                      numVertices = numVertices)), recursive = FALSE)

      return(newVList)

    } #  end createNewVList function

    #  return NULL if every pair of vertices has an edge between them (i.e. g is completely
    #  connected)
    if(length(vList) == 0){
      return(NULL)
    } #  end if

    #  track the number of vertices in each vector of vertices in vList
    numCurrentV = 2

    #  grow the size of the vertex vectors in vList iteratively
    while(numCurrentV < v){

      #  if there are no vectors in vList, return NULL
      if(length(vList) == 0){
        return(NULL)
      } #  end if

      #  append new values to each vector in vList
      vList = createNewVList(vList = vList, numVertices = numVertices)

      #  remove duplicate vectors (vectors corresponding to the same set of vertices)
      #  unnecessary and very slow
      #vList = uniqueSet(vList = vList, select = "ran")

      #  remove vectors from vList that correspond to subgraphs with at least one edge
      vList = vList[unlist(lapply(vList, isDisconnected))]

      #  if a value for numSample is specified, randomly select that many vectors from vList at each step
      #  this
      if(!is.null(numSample)){
        vList = sample(vList, numSample, replace = TRUE)
      } #  end if

      if(length(vList) == 0){
        vList = NULL
        break
      } #  end if

      numCurrentV = numCurrentV + 1

    } #  end while loop

    return(vList)

  } #  end produceSubgraphs function

  #  list to store vertex vectors
  vList = list()

  #  count the number of iterations so the function can terminate if it iterates too many times
  iter = 0

  #  produce the vertex vectors for vList
  while(length(vList) < numSubgraphs && iter <= maxIters){

    newSubs = produceSubgraphs()

    vList = append(vList, newSubs)

    iter = iter + 1

  } #  end while loop

  #  randomly select the desired number of vertex vectors from the created vList (which may
  #  be too long before this step)
  vList = sample(vList, numSubgraphs)

  return(vList)

} #  end disconnectedSubgraphs function


text:
##
##  Branch cut algorithm.
## This solves a variant of the assignment problem to minimize
## the total length of all branch cuts.
##
## note: parameter pen is a penalty added to distances from charges to edges.
##		 I'm not sure this is really useful. Seems to work fine with the default
##		 of 0.
##

#' Linear program based phase unwrapping
#'
#' `brcutpuw` implements a branch cut algorithm for phase unwrapping
#'
#' @details This package implements two distinct algorithms for two dimensional
#'   phase unwrapping that can be set up and solved as general linear programs.
#'   The linear programs are solved using the LP solver `CBC` from [COIN-OR](https://github.com/coin-or/Cbc)
#'   using the interface provided by the package [rcbc](https://dirkschumacher.github.io/rcbc/index.html).
#'
#' @param phase matrix with phase map to be unwrapped (phases are in radians)
#' @param pen penalty for making a branch cut from a residue to an edge (`brcutpuw` only)
#' @param wts matrix of weights for cost function with same dimension as `phase` (`netflowpuw` only)
#' @param details Return some details of the solution?
#' @param trace Send some info from the LP solver to the console?
#'
#' @return The unwrapped wavefront in units of fringes.
#'   If `details` is `TRUE` additional details of the solution
#'   are returned in a named list with the first member
#'   `puw` containing the unwrapped wavefront.
#'
#' @seealso There is a function with the same name `brcutpuw` in package [zernike].
#'
#' @describeIn brcutpuw Branch cut algorithm for phase unwrapping
#'
#' @section Note:
#'   According to the documentation for `rcbc` different levels of detail
#'   from the LP solver can be printed with trace levels up to 15,
#'   however the same output seems to be returned for all values.
#'   Setting `trace=0` will produce silent output, which may not
#'   be advisable since these can take some time to run.
#'
#'   The value of trace doesn't matter when running in a Windows "GUI"
#'   console window or Rstudio because the CBC output isn't passed. If
#'   you want to see some output run in a terminal window instead.
#'
#' @examples
#'   data("phasemaps", package="lppuw")
#'   mtext(zernike::rmap(phi, plot=TRUE))
#'   wf.bc <- brcutpuw(phi)
#'   wf.nf <- netflowpuw(phi, mod)
#'   X11()
#'   zernike::plot.pupil(wf.nf, col=zernike::rygcb(400))
#'   cat("Summary of the difference between the two unwrapped wavefronts:\n")
#'   zernike::summary.pupil(wf.nf - wf.bc)
#'
brcutpuw <- function(phase, pen=0, details=FALSE, trace=1) {

  ## distance between points specified by their x,y coordinates.

  dmat <- function(p1, p2) {
    fn <- function(x,y) abs(y-x)
    xleg <- outer(p1[,1], p2[,1], FUN=fn)
    yleg <- outer(p1[,2], p2[,2], FUN=fn)
    xleg+yleg
  }

  ## Make a branch cut

  brcut <- function(startp, endp, mask) {
    xl <- endp[1]-startp[1]
    yl <- endp[2]-startp[2]
    if (abs(yl) <= abs(xl)) {
      slope <- yl/xl
      ix <- startp[1]:endp[1]
      iy <- startp[2] + round(slope*(ix-startp[1]))
    } else {
      slope <- xl/yl
      iy <- startp[2]:endp[2]
      ix <- startp[1] + round(slope*(iy-startp[2]))
    }
    mask[cbind(ix,iy)] <- NA
    mask
  }

  res <- zernike::rmap(phase)
  ## no residues, so a single call to idiffpuw is all we need
  if (sum(abs(res),na.rm=TRUE) == 0) {
    return(zernike::idiffpuw(phase,ucall=TRUE))
  }

  nr <- nrow(phase)
  nc <- ncol(phase)
  dx <- zernike::wrap(zernike::.fdiff(phase))
  dy <- zernike::wrap(t(zernike::.fdiff(t(phase))))

  cutlen <- NULL

  mask <- matrix(1,nr,nc)
  mask[is.na(phase)] <- NA

  ## positive charges

  chp <- which(res==1, arr.ind=TRUE)
  ncp <- nrow(chp)
  ## negative charges

  chm <- which(res== -1, arr.ind=TRUE)
  ncm <- nrow(chm)

  ## get the boundary points - ones one pixel outside the area of defined phase.
  ## step in

  ptsb <- which(is.na(phase) & (!is.na(rbind(phase[-1,], rep(NA, nc)))), arr.ind=TRUE)
  ptsb <- rbind(ptsb, which(is.na(phase) & (!is.na(cbind(phase[,-1], rep(NA, nr)))), arr.ind=TRUE))

  ## step out

  ptsc <- which((!is.na(phase)) & is.na(rbind(phase[-1,], rep(NA, nc))), arr.ind=TRUE)
  ptsc[,1] <- ptsc[,1]+1
  ptsb <- rbind(ptsb, ptsc)
  ptsc <- which((!is.na(phase)) & is.na(cbind(phase[,-1], rep(NA, nr))), arr.ind=TRUE)
  ptsc[,2] <- ptsc[,2]+1
  ptsb <- rbind(ptsb, ptsc)

  if ((ncp>0) & (ncm>0)) {

    ## distances between residues

    dpm <- dmat(chp, chm)

    ## distance from each positive charge to nearest edge

    dce <- dmat(chp, ptsb)
    ipb <- apply(dce, 1, which.min)
    dpe <- apply(dce, 1, min)+pen
    pe <- matrix(ptsb[ipb,], ncp, 2)

    ## distance from each negative charge to nearest edge

    dce <- dmat(chm, ptsb)
    ipb <- apply(dce, 1, which.min)
    dme <- apply(dce, 1, min)+pen
    me <- matrix(ptsb[ipb,], ncm, 2)

    ## set up the LP. This is almost the assignment problem.

    obj <- c(as.vector(dpm), dpe, dme)
    nrow <- ncp + ncm
    ncol <- ncp*ncm + ncp + ncm
    row_lb <- rep(1, nrow)
    i_ind <- c(rep(1:ncp, each=ncm+1), ncp+rep(1:ncm, each=ncp+1))
    j_mp <- t(matrix(1:ncp, ncp, ncm)) + matrix(ncp*(0:(ncm-1)), ncm, ncp)
    j_p <- rbind(j_mp, (ncm*ncp)+(1:ncp))
    j_p <- as.vector(j_p)
    j_m <- rbind(t(j_mp), (ncm*ncp + ncp)+(1:ncm))
    j_m <- as.vector(j_m)
    j_ind <- c(j_p, j_m)
    if (packageVersion("Matrix") >= "1.3.0") {
      mat <- Matrix::sparseMatrix(i=i_ind, j=j_ind, dims=c(nrow, ncol), repr="T")
    } else {
      mat <- Matrix::sparseMatrix(i=i_ind, j=j_ind, dims=c(nrow, ncol), giveCsparse=FALSE)
    }
    lpout <- rcbc::cbc_solve(obj=obj, mat=mat,
                            row_lb = row_lb,
                            row_ub = row_lb,
                            col_lb = rep(0, ncol),
                            col_ub = rep(1, ncol),
                            cbc_args = list(log=trace, verbose=trace)
                          )
    if (rcbc::solution_status(lpout) != "optimal") {
      warning("non-convergence detected")
      return(lpout)
    }

    lpsol <- lpout$column_solution
    cutlen <- lpout$objective_value
    dpcuts <- matrix(round(lpsol[1:(ncp*ncm)]), ncp, ncm)
    pecuts <- round(lpsol[(ncp*ncm+1):(ncp*ncm+ncp)])
    mecuts <- round(lpsol[(ncp*ncm+ncp+1):(ncp*ncm+ncp+ncm)])
    dpcuts <- which(dpcuts==1, arr.ind=TRUE)
    i <- 1
    while (i <= nrow(dpcuts)) {
      startp <- chp[dpcuts[i,1],]
      endp <- chm[dpcuts[i,2],]
      mask <- brcut(startp, endp, mask)
      i <- i+1
    }
    pecuts <- which(pecuts==1)
    i <- 1
    while (i <= length(pecuts)) {
      startp <- chp[pecuts[i],]
      endp <- pe[pecuts[i],]
      mask <- brcut(startp, endp, mask)
      i <- i+1
    }
    mecuts <- which(mecuts==1)
    i <- 1
    while (i <= length(mecuts)) {
      startp <- chm[mecuts[i],]
      endp <- me[mecuts[i],]
      mask <- brcut(startp, endp, mask)
      i <- i+1
    }
  } else if (ncp > 0) {
      ## if there are only plus or minus charges connect to nearest edge (should be rare case, but...)

      ## distance from each positive charge to nearest edge

      dce <- dmat(chp, ptsb)
      ipb <- apply(dce, 1, which.min)
      dpe <- apply(dce, 1, min)+pen
      pe <- matrix(ptsb[ipb,],ncp,2)
      for (i in 1:ncp) {
        startp <- chp[i,]
        endp <- pe[i,]
        mask <- brcut(startp, endp, mask)
      }
  } else if (ncm > 0) {

    ## distance from each negative charge to nearest edge

      dce <- dmat(chm, ptsb)
      ipb <- apply(dce, 1, which.min)
      dme <- apply(dce, 1, min)+pen
      me <- matrix(ptsb[ipb,],ncm,2)
      for (i in 1:ncm) {
        startp <- chm[i,]
        endp <- me[i,]
        mask <- brcut(startp, endp, mask)
      }
  }

  ## call idiffpuw to unwrap all unmasked pixels

  uwum <- zernike::idiffpuw(phase, mask, ucall=FALSE, dx=dx, dy=dy)
  puw <- uwum$puw
  uw <- uwum$uw
  uw[is.na(phase)] <- NA

  ## fill in whatever was left wrapped
  ## This inverts the logic of fill algorithm in idiffpuw
  ## todo initially contains a list of pixels that haven't been unwrapped
  ## that should be. We look for neighbors that have been unwrapped
  ## and unwrap from them, removing the newly unwrapped pixels from
  ## the todo list.
  ## This should work because branch cuts are just a pixel wide and should
  ## adjoin unmasked regions. It's possible that cuts could completely
  ## enclose an area, in which case it will still be filled
  ## with most likely erroneous values as we fill through the cuts.

  todo <- which(!uw)
  kn <- c(-1, 1, -nr, nr)
  while (length(todo) > 0) {
    for (i in 1:4) {
      b <- todo + kn[i]
      subs <- which(uw[b])
      puw[todo[subs]] <- puw[b[subs]] + switch(i, dx[b[subs]], -dx[todo[subs]],
                                                dy[b[subs]], -dy[todo[subs]])
      uw[todo[subs]] <- TRUE
      todo <- todo[-subs]
      if (length(todo)==0) break
    }
  }
  if (details) {
    bcuts <- matrix(NA, nr, nc)
    bcuts[is.na(mask) & (!is.na(phase))] <- 1
    list(puw=puw/(2*pi), cutlen=cutlen, bcuts=bcuts,
         obj=obj, mat=mat, lpout=lpout)
  } else {
    puw/(2*pi)
  }
}

max_stars_repo_name: mlpeck/lppuw

max_issues_repo_name: mlpeck/lppuw

max_forks_repo_name: mlpeck/lppuw


text:
for (x in 1:10) {
    print(x)
}

max_stars_repo_name: teddythinh/R-language-Practice

max_issues_repo_name: teddythinh/R-language-Practice

max_forks_repo_name: teddythinh/R-language-Practice


text:
y <- rnorm(n=36,mean=250,sd=2)
x <- sapply(y, function(y) rnorm(1,1.1*y,1))
write.table(round(x, digits=0), row.names = FALSE)


text:
# Import and format growth parameters
# then define growth function

# growth parameters
Gpars=list(intcpt=rep(NA,n_spp),intcpt.yr=matrix(0,Nyrs,n_spp), intcpt.gr=matrix(0,6,n_spp),
  slope=rep(NA,n_spp),slope.yr=matrix(0,Nyrs,n_spp),
  nb=matrix(0,n_spp,n_spp),alpha=matrix(NA,n_spp,n_spp),
  sigma2.a=rep(NA,n_spp),sigma2.b=rep(NA,n_spp))
for(i in 1:n_spp){
  infile <- paste("../extdata/Growth_params_",spp_list[i],".csv",sep="")
  Gdata=read.csv(infile)
  Gpars$intcpt[i]=Gdata$Intercept[1]
  tmp=which(names(Gdata)=="Group")
  if(length(tmp)>0) Gpars$intcpt.gr[,i]=Gdata$Group[!is.na(Gdata$Group)]
  Gpars$intcpt.yr[,i]=Gdata$Intercept.yr
  Gpars$slope[i]=Gdata$logarea.t0[1]
  # random effects on slope
  tmp=which(names(Gdata)=="logarea.t0.yr")
  if(length(tmp)>0) Gpars$slope.yr[,i]=Gdata[,tmp]
  # get competition coefficients
  tmp=paste("crowd",1:length(spp_list),sep="")
  tmp=which(is.element(names(Gdata),tmp))
  if(length(tmp)>0) Gpars$nb[i,]=as.numeric(Gdata[1,tmp])
  Gpars$alpha[i,]=Gdata$alpha[1:length(spp_list)]
  Gpars$sigma2.a[i]=Gdata$sigma.a[1]
  Gpars$sigma2.b[i]=Gdata$sigma.b[1]
} # next i
rm(Gdata)

# growth function
G=function(v,u,W,Gpars,doYear,doSpp){
  mu=Gpars$intcpt[doSpp]+Gpars$intcpt.yr[doYear,doSpp]+(Gpars$slope[doSpp]+Gpars$slope.yr[doYear,doSpp])*u+
    W%*%(Gpars$nb[doSpp,])
  sigma2=Gpars$sigma2.a[doSpp]*exp(Gpars$sigma2.b[doSpp]*mu)
  out=dnorm(v,mu,sqrt(sigma2))
  out
}

max_stars_repo_name: atredennick/IPMdoit

max_issues_repo_name: atredennick/IPMdoit

max_forks_repo_name: atredennick/IPMdoit


text:
#  James Rekow

marginalDissDensity = function(M = 40, numSubgroups = 2, N = 20, iStrength = 1, univ = 1, sigmaMax = 0.1,
                                    thresholdMult = 10 ^ (-1), maxSteps = 10 ^ 4, tStep = 10 ^ (-2),
                                    intTime = 1, interSmplMult = 0.01, lambda = 0, returnParams = FALSE,
                                    conGraph = NULL, overRange = 0.05, highOver = 0.9){

  #  ARGS: numSubgroups - number of subgroups in the cohort
  #        rollRange - the maximum difference in overlap between two DOC points for which the two points
  #                    are still considered close enough that they are included in each others local
  #                    dissimilarity profiles (e.g. in computing rolling mean, etc.)
  #
  #  RETURNS:

  library(ggplot2)
  library(grid)
  library(gridExtra)
  source("subgroupsAbdListCreator.r")
  source("DOCProcedure.r")

  #  default indices of the 2 subgroups created from subgroupsAbdListCreator
  subgroup1Ix = 1:20
  subgroup2Ix = 21:40

  #  create abdList with desired number of subgroups
  abdList = subgroupsAbdListCreator(M = M, numSubgroups = numSubgroups, N = N, iStrength = iStrength,
                                    univ = univ, sigmaMax = sigmaMax, thresholdMult = thresholdMult,
                                    maxSteps = maxSteps, tStep = tStep, intTime = intTime,
                                    interSmplMult = interSmplMult, lambda = lambda,
                                    returnParams = returnParams, conGraph = conGraph)

  #  compute DOC
  doc = DOCProcedure(abdList)

  #  make DOC a data frame (by default a doc is a named list)
  doc = data.frame(over = doc$x, diss = doc$y)

  #  create matrix whose ith column is the pair of indices corresponding to the ith DOC point in the doc
  #  data frame
  ixMat = combn(1:M, 2)

  #  check whether either of the indices corresponding to a given DOC point come from a sample in
  #  subgroup 1(2)
  ixInSubgroup1 = apply(ixMat, 2, function(x) {x[1] %in% subgroup1Ix} || {x[2]} %in% subgroup1Ix)
  ixInSubgroup2 = apply(ixMat, 2, function(x) {x[1] %in% subgroup2Ix} || {x[2]} %in% subgroup2Ix)

  #  check which DOC points come from samples in the same or different subgroups
  diffSubgroup = ixInSubgroup1 & ixInSubgroup2
  sameSubgroup = !diffSubgroup

  #  add column to data frame
  doc$sameSubgroup = sameSubgroup

  #  create ggplot coloring DOC points based on classification as same or different subgroups (of
  #  corresponding samples)
  p1 = ggplot(doc, aes(x = over, y = diss, col = sameSubgroup)) + geom_point()

  ggplotSubgroupDissDens = function(over){

    #  ARGS: over - numeric overlap value
    #
    #  RETURNS: densPlot - ggplot object plotting the dissimilarity density of DOC points with overlap
    #                      similar to over, grouped by sameSubgroup

    #  identify indices of DOC points with similar overlap as input overlap value
    similarOver = abs(doc$over - over) < overRange

    #  create DOC data frame with only the points that have similar overlap to the input over value
    similarOverDOC = subset(doc, similarOver == TRUE)

    #  create ggplot object plotting the two marginal dissimilarity densities
    densPlot = ggplot(similarOverDOC, aes(diss, col = sameSubgroup)) + geom_density() + xlab("") + ylab("")

    return(densPlot)

  } #  end ggplotSubgroupDissDens

  #  compute mean overlap
  meanOver = mean(doc$over)

  #  create ggplot objects for the marginal dissimilarity distributions
  meanOverDissDensPlot = ggplotSubgroupDissDens(over = meanOver)
  highOverDissDensPlot = ggplotSubgroupDissDens(over = highOver)

  titleText = "Marginal Dissimilarity Densities"
  p2 = grid.arrange(arrangeGrob(meanOverDissDensPlot, top = "Mean Overlap"),
                    arrangeGrob(highOverDissDensPlot, top = paste("High Overlap: ", highOver)),
                    ncol = 2, left = "Density", bottom = "Dissimilarity", top = titleText)

  #  plot first ggplot object in new graphics device
  dev.new()
  plot(p1)

  #  plot second ggplot object in new graphics device
  dev.new()
  plot(p2)

  return(doc)

} #  end marginalDissDensity function


text:
#' rel_acclimatized
#'
#' Calculate treshshold of risk realated to wbgt index for acclimatized people
#'
#' @param met numeric Metabolic rate in Watt following level of ISO FDIS 7243.#
#' @return REL
#'
#' @author Istituto per la Bioeconomia CNR Firenze Italy  Alfonso Crisci \email{alfonso.crisci@@ibe.cnr.it}
#' @references  Ergonomics of the thermal environment – Assessment of heat stress using the WBGT (wet bulb globe temperature)   index, ISO FDIS 7243 (2016)
#' @keywords  REL
#' @export
#'
#'
#'
#'

rel_acclimatized=function(met)  {
  ct$assign("met", as.array(met))
  ct$eval("var res=[]; for(var i=0, len=met.length; i < len; i++){ res[i]=rel_acclimatized(met[i])};")
  res=ct$get("res")
  return(res)
}

max_stars_repo_name: alfcrisci/rBiometeo

max_issues_repo_name: alfcrisci/rBiometeo

max_forks_repo_name: alfcrisci/rBiometeo


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

max_stars_repo_name: gspu/Coherent

max_issues_repo_name: gspu/Coherent

max_forks_repo_name: gspu/Coherent


text:
# A gibbs sampler for the DP-GMM model presented in
# https://www.seas.harvard.edu/courses/cs281/papers/rasmussen-1999a.pdf
# author: Alberto Lumbreras
#
# There is a matlab implementation:
# https://github.com/mim/igmm/blob/master/igmm_mv.m

library(coda)
library(mixtools)
source('conditionals.r')

gibbs <- function(A, z_init, iters=100, plots=TRUE, plots.freq=10){
  # Arguments:
  #   A: matrix with user features (one column per user)
  #   z_init: initial assignments (one per user)
  ######################################################
  par(mfrow=c(1,1))

  D <- dim(A)[1]
  N <- dim(A)[2]

  # traces
  traces.z <- matrix(NA, iters, N)
  traces <- matrix(NA, iters, 6)
  colnames(traces) <- c("n_clusters", "alpha", "mean_mu_a0", "det(R_a0)", "det(W_a0)", "beta_a0")

  # data mean and covariance
  Sigma_a <- cov(t(A))
  Lambda_a <- solve(Sigma_a)
  mu_a <- rowMeans(A)

  # init state
  alpha = 2
  z <- z_init
  mu_a0 <- mu_a
  R_a0 <- solve(Sigma_a)
  W_a0 <- Sigma_a
  beta_a0 <- dim(A)[1]
  S_ar <- array(NA, dim=c(D, D, 100))
  mu_ar <- matrix(NA, D, 100)
  for (k in 1:length(unique(z))){
    mu_ar[,k] <- rowMeans(A[,z==k])
  }

  for(i in 1:iters){
    cat('\n', i)

    # Active components
    K <- length(unique(z))

    # Sample component parameters
    for (k in 1:length(unique(z))){
      mask <- (z==k)
      S_ar[,,k] <- sample_S_ar(A[,mask, drop=FALSE], mu_ar[,k, drop=FALSE], beta_a0, W_a0)
      mu_ar[,k] <- sample_mu_ar(A[,mask, drop=FALSE], S_ar[,,k], mu_a0, R_a0)
    }

    # Sample hyperpriors
    mu_a0 <- sample_mu_a0(Lambda_a, mu_a, mu_ar[,1:K, drop=FALSE], R_a0)
    R_a0 <- sample_R_a0(Sigma_a, mu_ar[,1:K, drop=FALSE], mu_a0)
    W_a0 <- sample_W_a0(Lambda_a, S_ar[,,1:K, drop=FALSE], beta_a0)
    beta_a0 <- sample_beta_a0(S_ar[,,1:K, drop=FALSE], W_a0, init=beta_a0)

    # Sample assignments
    for (n in 1:N){
      res <- sample_z(n, A, alpha, z, mu_ar, S_ar, mu_a0, R_a0,  beta_a0, W_a0)
      z <- res$z
      mu_ar <- res$mu_ar
      S_ar <- res$S_ar
    }
    cat('\n', table(z))

    # Sample concentration parameter
    K <- length(unique(z))
    alpha <- sample_alpha(K = K, U = N)

    traces[i,] <- c(length(unique(z)), alpha, mean(mu_a0), det(R_a0), det(W_a0), beta_a0)
    traces.z[i,] <- z

    if(plots && (i%%plots.freq==0)){
      plot(t(A[c(1,2),]), col=palette()[z+1])
      title(paste0("Sample #", i))
      for(k in 1:K){
        ellipse(mu=mu_a0[c(1,2)], sigma=W_a0[c(1,2),c(1,2)], alpha = .25, lwd=5, npoints = 250, col="black")
        ellipse(mu=mu_ar[c(1,2),k], sigma=solve(S_ar[c(1,2),c(1,2),k]),  alpha = .25, lwd=3, npoints = 250, col=palette()[k+1])
      }
    }

  }
  res <- list(traces=traces,
              traces.z=traces.z)

  return(res)
}

max_stars_repo_name: alumbreras/DP-GMM

max_issues_repo_name: alumbreras/DP-GMM

max_forks_repo_name: alumbreras/DP-GMM