text:
#' Write json pretty
#'
#' Re-write existing json file with pretty formatting
#'
#' @param path path to json file on disk
#' @keywords internal
#' @examples
#' \dontrun{
#' pretty_json("tests/testthat/api.coinbase.com/v2/accounts.json")
#' }
pretty_json <- function(path) {
  json <- jsonlite::read_json(path)
  jsonlite::write_json(json, path, pretty = TRUE, auto_unbox = TRUE)
}

#' Capture requests to tests directory
#'
#' Helper to build up mock fixtures for testing
#' @param expr R expresssion to execute
#' @examples
#' \dontrun{
#' capture_requests2(get_accounts())
#' }
capture_request2 <- function(expr) {
  httptest::capture_requests(expr, path = "tests/testthat")
}

max_stars_repo_name: zamorarr/gdaxr

max_issues_repo_name: zamorarr/rcoinbase

max_forks_repo_name: zamorarr/rcoinbase


text:
source("mtcars__dane.r")
source("machine_dane.r")
source("imports-85_dane.r")

#zobacz jak w SVMRegresja.r

ext: r

max_stars_repo_name: PLPOLAND/Projekt-R

max_issues_repo_name: PLPOLAND/Projekt-R

max_forks_repo_name: PLPOLAND/Projekt-R


text:
#!/usr/bin/env Rscript

# library(Biobase)
suppressPackageStartupMessages(library(pcaMethods, quietly = TRUE))
suppressPackageStartupMessages(library(colorspace, quietly = TRUE))
suppressPackageStartupMessages(library(RColorBrewer, quietly = TRUE))
suppressPackageStartupMessages(library(sva, quietly = TRUE))
suppressPackageStartupMessages(library(MASS, quietly = TRUE))
suppressPackageStartupMessages(library(NMF, quietly = TRUE))
suppressPackageStartupMessages(library(dendextend, quietly = TRUE))
suppressPackageStartupMessages(library(amap, quietly = TRUE))
suppressPackageStartupMessages(library(pvclust, quietly = TRUE))
suppressPackageStartupMessages(library(Rtsne, quietly = TRUE))

debug_mode = ifelse(length(commandArgs(trailingOnly = TRUE)) == 1, "debug", "batch")

log_prefix = "transcriptome_curation.r:"
cat(log_prefix, "mode =", debug_mode, "\n")
if (debug_mode == "debug") {
    infile = '/Users/s229181/MSN/cstmm/cross_species_tmm_normalized_counts/Helianthus_annuus_cstmm_counts.tsv'
    eff_file = '/Users/s229181/MSN/eff_length/Helianthus_annuus_eff_length.tsv'
   # infile = "/Users/kf/Dropbox (Personal)/collaborators/Ken Naito/20210509_Vigna/gfe_data/cross_species_tmm_normalized_counts/Vigna_angularis_cstmm_counts.tsv"
    #eff_file = "/Users/kf/Dropbox (Personal)/collaborators/Ken Naito/20210509_Vigna/gfe_data/merge/Vigna_angularis_eff_length.tsv"
    dist_method = "pearson"
    mapping_rate_cutoff = 0
    min_dif = 0
    plot_intermediate = 0
    selected_curate_group = c("root", "flower", "leaf")
    dir_count = "counts/"
    dir_eff_length = "eff_length/"
    mode = "msm"
    transform_method = "fpkm"
    # tmm norm debug
    correlation_threshold = 0.4
    tmm_norm = "yes"
    dir_work = '/Users/s229181/MSN/'
    srafile = '/Users/s229181/MSN/metadata/metadata_may2021_updated.tsv'
   # dir_work = "/Users/kf/Dropbox (Personal)/collaborators/Ken Naito/20210509_Vigna/gfe_data"
   # srafile = "/Users/kf/Dropbox (Personal)/collaborators/Ken Naito/20210509_Vigna/gfe_data/metadata/metadata/metadata_manual.tsv"

    dist_method = "pearson"
    min_dif = 0
    plot_intermediate = 1
    selected_curate_groups = c("root", "flower", "leaf")
    stop_after_tmm = TRUE

    # selected_curate_groups = strsplit('root|flower|leaf', '\\|')[[1]]
} else if (debug_mode == "batch") {
    args = commandArgs(trailingOnly = TRUE)
    print(args)
    infile = args[1]
    srafile = args[2]
    dir_work = args[3]
    eff_file = args[4]
    dist_method = args[5]
    mapping_rate_cutoff = as.numeric(args[6])
    min_dif = as.numeric(args[7])
    plot_intermediate = as.integer(args[8])
    selected_curate_groups = strsplit(args[9], "\\|")[[1]]
    transform_method = args[10]
    one_outlier_per_iteration = as.integer(args[11])
    correlation_threshold = as.numeric(args[12])

}
if (!endsWith(dir_work, "/")) {
    dir_work = paste0(dir_work, "/")
}

# set directory
cat(log_prefix, "dir_work =", dir_work, "\n")
dir.create(dir_work, showWarnings = FALSE)

dir_curate = file.path(dir_work,'curate')
if (!file.exists(dir_curate)) {
  dir.create(dir_curate)
}
setwd(dir_curate)

# create output directories
dir_pdf = file.path(dir_curate,'plots')
if (!file.exists(dir_pdf)) {
  dir.create(dir_pdf)
}
dir_rdata = file.path(dir_curate,'rdata')
if (!file.exists(dir_rdata)) {
  dir.create(dir_rdata)
}
dir_tsv = file.path(dir_curate,'tables')
if (!file.exists(dir_tsv)) {
  dir.create(dir_tsv)
}

tc_sra_intersect = function(tc, sra) {
    sra_run = sra[['run']]
    tc = tc[, colnames(tc) %in% sra_run]
    sra = sra[sra[['run']] %in% colnames(tc), ]
    return(list(tc = tc, sra = sra))
}

remove_nonexpressed_gene = function(tc) {
    gene_sum = apply(tc, 1, sum)
    tc_ex = tc[gene_sum > 0, ]
    tc_ne = tc[gene_sum == 0, ]
    return(list(tc_ex = tc_ex, tc_ne = tc_ne))
}

add_color_to_sra = function(sra, selected_curate_groups) {
    sra = sra[, (!names(sra) %in% c("bp_color", "sp_color", "curate_group_color"))]
    bioproject = as.character(sra[['bioproject']])
    bioproject_u = sort(unique(bioproject))
    scientific_name = as.character(sra[['scientific_name']])
    scientific_name_u = sort(unique(scientific_name))
    curate_group = as.character(sra[['curate_group']])
    curate_group_u = sort(unique(curate_group))
    if (length(selected_curate_groups) <= 8) {
        curate_group_color = brewer.pal(length(curate_group_u), "Dark2")
        bp_color = rainbow_hcl(length(bioproject_u), c = 50)
        sp_color = rainbow_hcl(length(scientific_name_u), c = 100)
    } else if (length(selected_curate_groups) <= 12) {
        curate_group_color = brewer.pal(length(unique(curate_group)), "Paired")
        bp_color = rainbow_hcl(length(bioproject_u), c = 50)
        sp_color = rainbow_hcl(length(scientific_name_u), c = 100)
    } else {
        curate_group_color = rainbow_hcl(length(selected_curate_groups), c = 100)
        bp_color = rainbow_hcl(length(bioproject_u), c = 50)
        sp_color = rainbow_hcl(length(scientific_name_u), c = 150)
    }
    df_curate_group = data.frame(curate_group=curate_group_u, curate_group_color=curate_group_color[1:length(curate_group_u)], stringsAsFactors = FALSE)
    df_bp = data.frame(bioproject=bioproject_u, bp_color=bp_color[1:length(bioproject_u)], stringsAsFactors = FALSE)
    df_sp = data.frame(scientific_name=scientific_name_u, sp_color=sp_color[1:length(scientific_name_u)], stringsAsFactors = FALSE)
    sra = merge(sra, df_bp, sort = FALSE, all.y = FALSE)
    sra = merge(sra, df_sp, sort = FALSE, all.y = FALSE)
    sra = merge(sra, df_curate_group, sort = FALSE, all.y = FALSE)
    return(sra)
}

sort_tc_and_sra = function(tc, sra, sort_columns = c("curate_group", "scientific_name", "bioproject")) {
    for (column in rev(sort_columns)) {
        sra = sra[order(sra[column]), ]
    }
    sra_intersection = sra[(sra[['run']] %in% colnames(tc)),'run']
    tc = tc[, sra_intersection]
    return(list(tc = tc, sra = sra))
}

sort_averaged_tc = function(tc) {
    split_colnames = strsplit(colnames(tc), "_")
    genus_names = c()
    specific_names = c()
    curate_group_names = c()
    for (i in 1:length(split_colnames)) {
        genus_names = c(genus_names, split_colnames[[i]][1])
        specific_names = c(specific_names, split_colnames[[i]][2])
        curate_group_names = c(curate_group_names, split_colnames[[i]][3])
    }
    colname_order = order(curate_group_names, genus_names, specific_names)
    tc = tc[, colname_order]
    return(tc)
}

cleanY = function(y, mod, svs) {
    X = cbind(mod, svs)
    Hat = solve(t(X) %*% X) %*% t(X)
    beta = (Hat %*% t(y))
    P = ncol(mod)
    return(y - t(as.matrix(X[, -c(1:P)]) %*% beta[-c(1:P), ]))
}

curate_group_mean = function(tc, sra, selected_curate_groups = NA, balance.bp = FALSE) {
    # if the data are SVA-corrected, balance.bp would not be necessary because project-specfic effects
    # were removed in SVA already.
    if (all(is.na(selected_curate_groups))) {
        sp_curate_groups = unique(sra[['curate_group']])
    } else {
        sp_curate_groups = selected_curate_groups[selected_curate_groups %in% unique(sra[['curate_group']])]
    }
    tc_ave = data.frame(matrix(rep(NA, length(sp_curate_groups) * nrow(tc)), nrow = nrow(tc)))
    colnames(tc_ave) = sp_curate_groups
    rownames(tc_ave) = rownames(tc)
    for (curate_group in sp_curate_groups) {
        run_curate_group = sra[(sra[['curate_group']] == curate_group),'run']
        run_curate_group = run_curate_group[run_curate_group %in% colnames(tc)]
        if (length(run_curate_group) == 1) {
            exp_curate_group = tc[, run_curate_group]
        } else {
            if (balance.bp) {
                is_run = (sra[['run']] %in% colnames(tc))
                is_curate_group = (sra[['curate_group']] == curate_group)
                is_no_exclusion = (sra[['exclusion']] == "no")
                bps = unique(sra[is_run & is_curate_group & is_no_exclusion, "bioproject"])
                df_tmp = data.frame(matrix(rep(NA, nrow(tc) * length(bps)), nrow = nrow(tc), ncol = length(bps)))
                colnames(df_tmp) = bps
                for (bp in bps) {
                    is_bp = (sra[['bioproject']] == bp)
                    is_curate_group = (sra[['curate_group']] == curate_group)
                    is_no_exclusion = (sra[['exclusion']] == "no")
                    sra_ids = sra[is_bp & is_curate_group & is_no_exclusion, "run"]
                    tc_bp = tc[, sra_ids]
                    if (class(tc_bp) == "numeric") {
                        df_tmp[bp] = tc_bp
                    } else {
                        df_tmp[bp] = rowMeans(tc_bp)
                    }
                }
                exp_curate_group = rowMeans(df_tmp)
            } else {
                exp_curate_group = rowMeans(tc[, run_curate_group])
            }
        }
        tc_ave[, curate_group] = exp_curate_group
    }
    return(tc_ave)
}

curate_group2tau = function(tc_curate_group, rich.annotation = TRUE, unlog = FALSE) {
    if (rich.annotation) {
        cols = c("tau", "highest", "order")
    } else {
        cols = c("tau")
    }
    df_tau = data.frame(matrix(rep(NA, length(cols) * nrow(tc_curate_group)), nrow = nrow(tc_curate_group)))
    colnames(df_tau) = cols
    rownames(df_tau) = rownames(tc_curate_group)
    if (unlog) {
        tc_curate_group = exp(tc_curate_group) - 1
        tc_curate_group[tc_curate_group < 0] = 0
    }
    xmax = apply(tc_curate_group, 1, max)
    df_tau[,'tau'] = apply((1 - (tc_curate_group/xmax))/(ncol(tc_curate_group) - 1), 1, sum)
    if (rich.annotation) {
        tc_curate_group[is.na(tc_curate_group)] = 0
        for (i in 1:nrow(tc_curate_group)) {
            is_nonzero = tc_curate_group[i, ] > 0
            if (sum(is_nonzero) > 0) {
                exp_order = order(tc_curate_group[i, is_nonzero], decreasing = TRUE)
                curate_group_ordered = colnames(tc_curate_group)[is_nonzero][exp_order]
                df_tau[i, "highest"] = curate_group_ordered[1]
                df_tau[i, "order"] = paste(curate_group_ordered, collapse = "|")
            }
        }
    }
    return(df_tau)
}

check_mapping_rate = function(tc, sra, mapping_rate_cutoff) {
    cat(paste0('Mapping rate cutoff: ', mapping_rate_cutoff*100, '%\n'))
    is_mapping_good = (sra[['mapping_rate']] >= mapping_rate_cutoff*100)
    is_mapping_good[is.na(is_mapping_good)] = TRUE
    if (any(!is_mapping_good)) {
        cat("Removed due to low mapping rate:\n")
        df_tmp = sra[!is_mapping_good,]
        for (i in rownames(df_tmp)) {
            sra_id = df_tmp[i,'run']
            mapping_rate = df_tmp[i,'mapping_rate']
            cat(paste0(sra_id, ': mapping rate = ', mapping_rate, '%\n'))
        }
        tc = tc[, colnames(tc) %in% sra[is_mapping_good, "run"]]
    } else {
        cat("No entry removed due to low mapping rate.\n")
    }
    sra[!is_mapping_good, "exclusion"] = "low_mapping_rate"
    return(list(tc = tc, sra = sra))
}

check_within_curate_group_correlation = function(tc, sra, dist_method, min_dif, selected_curate_groups, one_out_per_iter = TRUE, correlation_threshold) {
    out = tc_sra_intersect(tc, sra)
    tc = out[["tc"]]
    sra2 = out[["sra"]]

    sra2[,'num_other_run_same_bp_curate_group'] = 0
    selected_curate_groups = selected_curate_groups[selected_curate_groups %in% unique(sra2[['curate_group']])]
    num_curate_group = length(selected_curate_groups)
    exclude_runs = c()
    for (sra_run in colnames(tc)) {
        is_sra = (sra2[['run']] == sra_run)
        my_curate_group = sra2[is_sra, "curate_group"]
        my_bioproject = sra2[is_sra, "bioproject"]
        is_not_my_bp = (sra2[['bioproject']] != my_bioproject)
        is_my_curate_group = (sra2[['curate_group']] == my_curate_group)
        run_other_bp = sra2[(is_not_my_bp | !is_my_curate_group), "run"]
        run_other_bp = run_other_bp[run_other_bp %in% colnames(tc)]
        tc_other_bp = tc[, run_other_bp]
        num_other_run_same_bp_curate_group = length(unique(sra2[(is_not_my_bp & is_my_curate_group), "bioproject"]))
        sra2[is_sra, "num_other_run_same_bp_curate_group"] = num_other_run_same_bp_curate_group
        sra2_other_bp <- sra2[sra2[['run']] %in% run_other_bp,]

        # If one curate_group is completely sourced from the same bioproject, we can't remove the whole bioproject for tc_ave_other_bp

        num_other_bp_same_curate_group = sum(sra2_other_bp[['curate_group']] == my_curate_group, na.rm=TRUE)
        if (num_other_bp_same_curate_group == 0) {
            tc_ave_other_bp = curate_group_mean(tc, sra2, selected_curate_groups)
        } else {
            tc_ave_other_bp = curate_group_mean(tc_other_bp, sra2, selected_curate_groups)
        }
        tc_ave = curate_group_mean(tc, sra2, selected_curate_groups)
        coef = c()
        coef_other_bp = c()
        for (curate_group in selected_curate_groups) {
            tmp_coef = cor(tc[, sra_run], tc_ave[, curate_group], method = dist_method)

            if ((num_other_bp_same_curate_group == 0) & (curate_group == my_curate_group)) {
              tmp_coef_other_bp = cor(tc[, sra_run], tc_ave[, curate_group], method = dist_method)
            } else {
              tmp_coef_other_bp = cor(tc[, sra_run], tc_ave_other_bp[, curate_group], method = dist_method)
            }

            if (curate_group == my_curate_group) {
                tmp_coef = tmp_coef - min_dif
                tmp_coef_other_bp = tmp_coef_other_bp - min_dif
            }
            coef = c(coef, tmp_coef)
            coef_other_bp = c(coef_other_bp, tmp_coef_other_bp)
        }
        names(coef) = selected_curate_groups
        names(coef_other_bp) = selected_curate_groups
        if (max(coef) != coef[my_curate_group]) {
            cat('Registered as a candidate for exclusion. Better correlation to other categories:', sra_run, '\n')
            exclude_runs = c(exclude_runs, sra_run)
        }
        if (coef_other_bp[my_curate_group] < correlation_threshold) {
            cat('Registered as a candidate for exclusion. Low within-category correlation:', sra_run, '\n')
            exclude_runs = c(exclude_runs, sra_run)
        }
    }
    if (length(exclude_runs)) {
      if(one_out_per_iter == TRUE){
          cat("Excluding only one outlier per bioproject or same curate_group. \n")
          exclude_run_bps_and_curate_group = sra2[(sra2$run %in% exclude_runs), c("bioproject", "run", "curate_group")]
          first_bp_hit = exclude_run_bps_and_curate_group[match(unique(exclude_run_bps_and_curate_group$bioproject), exclude_run_bps_and_curate_group$bioproject),]
          first_same_curate_group_hit = exclude_run_bps_and_curate_group[match(unique(exclude_run_bps_and_curate_group$curate_group), exclude_run_bps_and_curate_group$curate_group),]
          # if a first_same_curate_group_hit is part of the same bioproject as the other removal candidates, ommit the same curate_group candidates
          if(any(first_same_curate_group_hit$bioproject %in% first_bp_hit$bioproject)){
              exclude_runs_tmp = c(first_bp_hit$run,first_same_curate_group_hit[!first_same_curate_group_hit$bioproject %in% first_bp_hit$bioproject]$run)
          }else{
              exclude_runs_tmp = c(first_bp_hit$run,first_same_curate_group_hit$run)
          }
          exclude_runs = unique(exclude_runs_tmp)
          # TODO This boolean vector should be all TRUE by definition. ???: exclude_run_bps_and_curate_group$run %in% exclude_runs
          exclude_bps = exclude_run_bps_and_curate_group[exclude_run_bps_and_curate_group$run %in% exclude_runs, "bioproject"]
      } else {
        exclude_run_bps = sra2[(sra2[['run']] %in% exclude_runs), c("bioproject", "run", "num_other_run_same_bp_curate_group")]
        exclude_bp_counts = data.frame(table(exclude_run_bps[['bioproject']]))
        exclude_run_bps = merge(exclude_run_bps, exclude_bp_counts, by.x = "bioproject", by.y = "Var1")
        exclude_run_bps = exclude_run_bps[order(exclude_run_bps[['num_other_run_same_bp_curate_group']], exclude_run_bps[['Freq']]),]
        rownames(exclude_run_bps) = 1:nrow(exclude_run_bps)
        min_other_run_same_bp_curate_group = exclude_run_bps[1, "num_other_run_same_bp_curate_group"]
        semimin_bp_count = exclude_run_bps[1, "Freq"]
        cat("minimum number of other BioProjects within curate_group:", min_other_run_same_bp_curate_group, "\n")
        cat("semi-minimum count of exclusion-candidate BioProjects:", semimin_bp_count, "\n")
        conditions = (exclude_run_bps[['Freq']] == semimin_bp_count)
        conditions = conditions & (exclude_run_bps[['num_other_run_same_bp_curate_group']] == min_other_run_same_bp_curate_group)
        exclude_bps = unique(exclude_run_bps[conditions, "bioproject"])
        exclude_runs = exclude_run_bps[(exclude_run_bps[['bioproject']] %in% exclude_bps), "run"]
      }

    }
    if (length(exclude_runs)) {
        cat("Partially removed BioProjects due to low within-curate_group correlation:\n")
        print(exclude_bps)
        cat("Removed Runs due to low within-curate_group correlation:\n")
        print(exclude_runs)
    }
    tc = tc[, !colnames(tc) %in% exclude_runs]
    sra[(sra[['run']] %in% exclude_runs), "exclusion"] = "low_within_curate_group_correlation"
    return(list(tc = tc, sra = sra))
}

sva_subtraction = function(tc, sra) {
    out = tc_sra_intersect(tc, sra)
    tc = out[["tc"]]
    sra = out[["sra"]]
    out = sort_tc_and_sra(tc, sra)
    tc = out[["tc"]]
    sra = out[["sra"]]
    out = remove_nonexpressed_gene(tc)
    tc = out[["tc_ex"]]
    tc_ne = out[["tc_ne"]]
    mod = model.matrix(~curate_group, data = sra)
    mod0 = model.matrix(~1, data = sra)
    set.seed(1)
    sva1 = try(sva(dat = as.matrix(tc), mod = mod, mod0 = mod0, B = 10))
    if (class(sva1) != "try-error") {
        cat("SVA correction was correctly performed.\n")
        tc = cleanY(y = tc, mod = mod, svs = sva1[['sv']])
        tc = rbind(tc, tc_ne)
    } else {
        cat("SVA correction failed. Returning the original transcriptome data.\n")
        tc = rbind(tc, tc_ne)
    }
    return(list(tc = tc, sva = sva1))
}

map_color = function(redundant_variables, c) {
    uniq_var = unique(redundant_variables)
    uniq_col = rainbow_hcl(length(uniq_var), c = c)
    df_unique = data.frame(var = uniq_var, col = uniq_col, stringsAsFactors = FALSE)
    df_redundant = data.frame(var = redundant_variables, order = seq(1, length(redundant_variables)),
                              stringsAsFactors = FALSE)
    df_redundant = merge(df_redundant, df_unique, by = "var", all.x = TRUE, stringsAsFactors = FALSE)
    df_redundant = df_redundant[order(df_redundant[['order']]), ]
    return(df_redundant[['col']])
}

draw_heatmap = function(sra, tc_dist_matrix, legend = TRUE, fontsize = 7) {
    bp_fac = factor(sub(";.*", "", sra[, c("bioproject")]))
    curate_group_fac = factor(sra[, c("curate_group")])
    ann_label = data.frame(bioproject = bp_fac, curate_group = curate_group_fac)
    bp_col_uniq = unique(sra[order(sra[['bioproject']]), 'bp_color'])
    curate_group_col_uniq = unique(sra[order(sra[['curate_group']]), 'curate_group_color'])
    ann_color = list(bioproject = bp_col_uniq, curate_group = curate_group_col_uniq)
    breaks = c(0, seq(0.3, 1, 0.01))
    aheatmap(tc_dist_matrix, color = "-RdYlBu2:71", Rowv = NA, Colv = NA, revC = TRUE, legend = TRUE,
             breaks = breaks, annCol = ann_label, annRow = ann_label, annColors = ann_color, annLegend = legend,
             fontsize = fontsize)
}

color_children2parent = function(node) {
    if (length(node) == 2) {
        child1_color = attributes(node[[1]])[['edgePar']][["col"]]
        child2_color = attributes(node[[2]])[['edgePar']][["col"]]
        if ((!is.null(child1_color)) & (!is.null(child2_color))) {
            if (child1_color == child2_color) {
                attributes(node)[['edgePar']][["col"]] = child1_color
            }
        }
    }
    return(node)
}

draw_dendrogram = function(sra, tc_dist_dist, fontsize = 7) {
    dend <- as.dendrogram(hclust(tc_dist_dist))
    dend_colors = sra[order.dendrogram(dend),'curate_group_color']
    labels_colors(dend) <- dend_colors
    dend_labels <- sra[order.dendrogram(dend), 'run']
    dend <- color_branches(dend, labels = dend_labels, col = dend_colors)
    dend <- set(dend, "branches_lwd", 1)
    for (i in 1:ncol(tc)) {
        dend = dendrapply(dend, color_children2parent)
    }
    cex.xlab = min(fontsize, max(0.2, 0.5/log10(nrow(sra))))
    par(cex = cex.xlab)
    plot(dend, las = 1, axes = FALSE)
    par(cex = 1)
    axis(side = 2, line = 0, las = 1)
    mtext("Distance", side = 2, line = 8.5, outer = FALSE)
    n = nrow(sra)
    symbols(
      1:n,
      rep(0, n),
      circles = rep(1, n),
      add = TRUE,
      inches = 0.02,
      xpd = TRUE,
      lwd = 1,
      bg = sra[order.dendrogram(dend), 'curate_group_color'],
      fg = sra[order.dendrogram(dend), 'bp_color']
    )
}

draw_dendrogram_ggplot = function(sra, tc_dist_dist, fontsize = 7) {

  cg_col = unique(sra[,c('curate_group', 'curate_group_color')])
  bp_col = unique(sra[,c('bioproject', 'bp_color')])
  colnames(cg_col) = c('Group', 'Color')
  colnames(bp_col) = c('Group', 'Color')
  sra_colors = rbind(cg_col,bp_col)

  group_colors <-  data.table(Group=sra_colors$Group, Color=sra_colors$Color, key="Group")
  group_colors <- transpose(group_colors, make.names = "Group")

  hc       <- hclust(tc_dist_dist)           # heirarchal clustering
  dendr    <- dendro_data(hc, type="rectangle") # convert for ggplot
  clust.df <- data.frame(label=sra$run, curate_group=factor(sra[sra$run %in% dendr$labels$label,'curate_group']), bioproject = factor(sra[sra$run %in% dendr$labels$label,'bioproject']))
  dendr[["labels"]] <- merge(dendr[["labels"]],clust.df, by="label")
  ggplot() +
    geom_segment(data=dendr$segments, aes(x=x, y=y, xend=xend, yend=yend), size = .8, show.legend = FALSE) +
    geom_segment(data=merge(dendr$segments[dendr$segments$yend==0,], dendr$labels[,c('label', 'curate_group', 'x')], by = 'x'), aes(x=x, y=y, xend=xend, yend=yend, color = curate_group), size = .8, show.legend = FALSE) +
    geom_text(data=dendr$labels, aes(x, y - .008, label=label, hjust=0, angle = 270, color = curate_group ), size=3, show.legend = FALSE) +
    scale_y_continuous(expand = c(.2, .1)) +
    geom_point(data = dendr$labels, aes(x,y, color = bioproject), size = 3, show.legend = FALSE) +
    geom_point(data = dendr$labels, aes(x,y, color = curate_group), size = 2, show.legend = FALSE) +
    theme(axis.line.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.text.x=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.ticks.y=element_blank(),
          axis.text.y=element_blank(),
          panel.background=element_rect(fill="white"),
          panel.grid=element_blank()
    ) + scale_color_manual(values=group_colors)
}

draw_dendrogram_pvclust = function(sra, tc, nboot, pvclust_file, fontsize = 7) {
    dist_fun = function(x) {
        Dist(t(x), method = "pearson")
    }
    sp = sub(" ", "_", sra[['scientific_name']][1])
    if (file.exists(pvclust_file)) {
        if (file.info(pvclust_file)[['size']]) {
            print("pvclust file found.")
            load(pvclust_file)
        }
    } else {
        print("no pvclust file found. Start bootstrapping.")
        result = pvclust(tc, method.dist = dist_fun, method.hclust = "average", nboot = nboot, parallel = FALSE)  # UPGMA
        save(result, file = pvclust_file)
    }
    dend = as.dendrogram(result)
    dend_colors = sra[order.dendrogram(dend), 'curate_group_color']
    labels_colors(dend) = dend_colors
    dend_labels = sra[order.dendrogram(dend), 'run']
    dend = color_branches(dend, labels = dend_labels, col = dend_colors)
    dend = set(dend, "branches_lwd", 2)
    for (i in 1:ncol(tc)) {
        dend = dendrapply(dend, color_children2parent)
    }
    cex.xlab = min(0.2 + 1/log10(fontsize))
    par(cex = cex.xlab)
    plot(dend, las = 1, ylab = "Distance", cex.axis = 1/cex.xlab, cex.lab = 1/cex.xlab)
    par(cex = 1)
    n = nrow(sra)
    symbols(
      1:n,
      rep(0, n),
      circles = rep(1, n),
      add = TRUE,
      inches = 0.04,
      xpd = TRUE,
      lwd = 2,
      bg = sra[order.dendrogram(dend), 'curate_group_color'],
      fg = sra[order.dendrogram(dend), 'bp_color']
    )
    text(result, print.num = FALSE, cex = 1, col.pv = "black")
}

draw_pca = function(sra, tc_dist_matrix, fontsize = 7) {
    set.seed(1)
    pca = prcomp(tc_dist_matrix)
    xlabel = paste0("PC 1 (", round(summary(pca)[['importance']][2, 1] * 100, digits = 1), "%)")
    ylabel = paste0("PC 2 (", round(summary(pca)[['importance']][2, 2] * 100, digits = 1), "%)")
    plot(
      pca[['rotation']][, 1],
      pca[['rotation']][, 2],
      pch = 21,
      cex = 2,
      lwd = 1,
      bg = sra[['curate_group_color']],
      col = sra[['bp_color']],
      xlab = xlabel,
      ylab = ylabel,
      las = 1
    )
    # plot(pca$x[,1], pca$x[,2], pch=21, cex=2, lwd=2, bg=sra$curate_group_color, col=sra$bp_color, main=title,
    # xlab=xlabel, ylab=ylabel, las=1)
}

draw_mds = function(sra, tc_dist_dist, fontsize = 7) {
    set.seed(1)
    try_out = tryCatch({
        isoMDS(tc_dist_dist, k = 2, maxit = 100)
    }, error = function(a) {
        return("MDS failed.")
    })
    if (mode(try_out) == "character") {
        cat("MDS failed.\n")
        plot(c(0, 1), c(0, 1), ann = F, bty = "n", type = "n", xaxt = "n", yaxt = "n")
    } else {
        mds <- try_out
        plot(
          mds[['points']][, 1],
          mds[['points']][, 2],
          pch = 21,
          cex = 2,
          lwd = 1,
          bg = sra[['curate_group_color']],
          col = sra[['bp_color']],
          xlab = "MDS dimension 1",
          ylab = "MDS dimension 2",
          las = 1
        )
    }
}

draw_tsne = function(sra, tc, fontsize = 7) {
    perplexity = min(30, floor(nrow(sra)/4))
    set.seed(1)
    out_tsne = Rtsne(as.matrix(t(tc)), theta = 0, check_duplicates = FALSE, verbose = FALSE, perplexity = perplexity,
                     dims = 2)
    try_out = tryCatch({
        plot(
          out_tsne[['Y']][, 1],
          out_tsne[['Y']][, 2],
          pch = 21,
          cex = 2,
          lwd = 1,
          bg = sra[['curate_group_color']],
          col = sra[['bp_color']],
          xlab = "t-SNE dimension 1",
          ylab = "t-SNE dimension 2",
          las = 1
        )
    }, error = function(a) {
        return("t-SNE plot failed.")
    })
    if (mode(try_out) == "character") {
        cat("t-SNE failed.\n")
        plot(c(0, 1), c(0, 1), ann = F, bty = "n", type = "n", xaxt = "n", yaxt = "n")
    }
}

draw_sva_summary = function(sva_out, tc, sra, fontsize) {
    if ((is.null(sva_out)) | (class(sva_out) == "try-error")) {
        plot(c(0, 1), c(0, 1), ann = F, bty = "n", type = "n", xaxt = "n", yaxt = "n")
        df = NA
    } else {
        out = tc_sra_intersect(tc, sra)
        tc = out[["tc"]]
        sra = out[["sra"]]
        out = sort_tc_and_sra(tc, sra)
        tc = out[["tc"]]
        sra = out[["sra"]]
        cols = c("curate_group","bioproject","lib_selection","instrument","mapping_rate")
        label_cols = c("organ","BioProject","library selection","instrument","mapping rate")

        num_sv = sva_out[['n.sv']]
        df = data.frame(matrix(NA, num_sv, length(cols)))
        colnames(df) = cols
        rownames(df) = paste0("SV", 1:nrow(df))
        for (i in 1:length(cols)) {
            for (j in 1:num_sv) {
                if (length(unique(sra[, cols[i]])) == 1) {
                    df[j, i] = NA
                } else {
                    lm_summary = summary(lm(sva_out[['sv']][, j] ~ sra[, cols[i]]))
                    df[j, i] = lm_summary[['adj.r.squared']]
                }
            }
        }
        colnames(df) = label_cols
        breaks = seq(0, 1, 0.02)
        colors = colorRampPalette(c("blue", "yellow", "red"))(length(breaks))
        df2 = t(df)
        df2[df2 < 0] = 0
        aheatmap(df2, color = colors, Rowv = NA, Colv = NA, revC = TRUE, breaks = breaks, fontsize = fontsize)
    }
    return(df)
}

draw_boxplot = function(sra, tc_dist_matrix, fontsize = 7) {
    is_same_bp = outer(sra[['bioproject']], sra[['bioproject']], function(x, y) {
        x == y
    })
    is_same_curate_group = outer(sra[['curate_group']], sra[['curate_group']], function(x, y) {
        x == y
    })
    plot(c(0.5, 4.5), c(0, 1), type = "n", xlab = "", ylab = "Pearson's correlation\ncoefficient", las = 1,
         xaxt = "n")
    boxplot(tc_dist_matrix[(!is_same_bp) & (!is_same_curate_group)], at = 1, add = TRUE, col = "gray", yaxt = "n")
    boxplot(tc_dist_matrix[(is_same_bp) & (!is_same_curate_group)], at = 2, add = TRUE, col = "gray", yaxt = "n")
    boxplot(tc_dist_matrix[(!is_same_bp) & (is_same_curate_group)], at = 3, add = TRUE, col = "gray", yaxt = "n")
    boxplot(tc_dist_matrix[(is_same_bp) & (is_same_curate_group)], at = 4, add = TRUE, col = "gray", yaxt = "n")
    labels = c("bw\nbw", "bw\nwi", "wi\nbw", "wi\nwi")
    axis(side = 1, at = c(1, 2, 3, 4), labels = labels, padj = 0.5)
    axis(side = 1, at = 0.35, labels = "Organ\nBioProject", padj = 0.5, hadj = 1, tick = FALSE)

}

draw_tau_histogram = function(tc, sra, selected_curate_groups, fontsize = 7) {
    df_tau = curate_group2tau(curate_group_mean(tc, sra, selected_curate_groups), rich.annotation = FALSE, unlog = TRUE)
    hist_out = hist(df_tau[['tau']], breaks = seq(0, 1, 0.05), las = 1, xlab = "Tau (expression specificity)",
                    ylab = "Gene count", main = "", col = "gray")
    num_noexp = sum(is.na(df_tau[['tau']]))
    num_all = nrow(df_tau)
    # num_exp = nrow(df_tau) - num_noexp text_noexp = paste('Expressed genes:', num_exp,
    # '\nNon-expressed genes:', num_noexp)
    text_noexp = paste0("Excluded due to\nno expression:\n", num_noexp, "/", num_all, " genes")
    text(0, max(hist_out[['counts']]) * 0.85, text_noexp, pos = 4)
}

draw_exp_level_histogram = function(tc, sra, selected_curate_groups, fontsize = 7, transform_method) {
    tc_curate_group = curate_group_mean(tc, sra, selected_curate_groups)
    xmax = apply(tc_curate_group, 1, max)
    xmax[xmax < 0] = 0
    xmax[xmax > 15] = 15
    breaks = seq(0, 15, 1)
    hist_out = hist(xmax, breaks = breaks, las = 1, xlab = paste0("Max expression (log ",transform_method,"+1)"), ylab = "Gene count",
                    main = "", col = "gray")
}

draw_legend = function(sra, new = TRUE, pos = "center", fontsize = 7, nlabel.in.col) {
    if (new) {
        plot.new()
    }
    curate_group_unique = unique(sra[['curate_group']])
    bp_unique = unique(sub(";.*", "", sra[['bioproject']]))
    curate_group_color_unique = unique(sra[['curate_group_color']])
    bp_color_unique = unique(sra[['bp_color']])
    ncol = ceiling((length(curate_group_unique) + length(bp_unique) + 2)/nlabel.in.col)
    legend_text = c("Organ", as.character(curate_group_unique), "", "BioProject", as.character(bp_unique))
    legend_color = c(rgb(1, 1, 1, 0), rep(rgb(1, 1, 1, 0), length(curate_group_color_unique)), rgb(1, 1, 1,
                                                                                0), rgb(1, 1, 1, 0), bp_color_unique)
    legend_bg = c(rgb(1, 1, 1, 0), curate_group_color_unique, rgb(1, 1, 1, 0), rgb(1, 1, 1, 0), rep(rgb(1,
                                                                                1, 1, 0), length(bp_color_unique)))
    legend_font = c(2, rep(1, length(curate_group_color_unique)), 1, 2, rep(1, length(bp_color_unique)))
    legend(pos, legend = legend_text, pch = 21, lwd = 1, lty = 0, col = legend_color, pt.bg = legend_bg,
           text.font = legend_font, ncol = ncol, bty = "n")
}

save_plot = function(tc, sra, sva_out, dist_method, file, selected_curate_groups, fontsize = 7, transform_method) {
    out = tc_sra_intersect(tc, sra)
    tc = out[["tc"]]
    sra = out[["sra"]]
    sra = add_color_to_sra(sra, selected_curate_groups)
    out = sort_tc_and_sra(tc, sra)
    tc = out[["tc"]]
    sra = out[["sra"]]
    pdf(paste0(dir_pdf,'/',file, ".pdf"), height = 8, width = 7.2, fonts = "Helvetica", pointsize = fontsize)
    layout_matrix = matrix(c(2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1,
                             2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1,
                             1, 1, 1, 1, 1, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 3, 3,
                             3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 7, 7, 7, 8, 8, 8, 9, 9, 9,
                             9, 9, 9, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10,
                             10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10), 14, 12,
                           byrow = TRUE)
    layout(layout_matrix)
    tc_dist_matrix = cor(tc, method = dist_method)
    tc_dist_matrix[is.na(tc_dist_matrix)] = 0
    tc_dist_dist = Dist(t(tc), method = dist_method) + 1e-09
    tc_dist_dist[is.na(tc_dist_dist)] = 1
    par(mar = c(6, 6, 1, 0))
    draw_dendrogram(sra, tc_dist_dist, fontsize)
    par(mar = c(0, 0, 0, 0))
    draw_heatmap(sra, tc_dist_matrix, legend = FALSE)
    # draw_dendrogram(sra, tc, nboot=1000, cex.xlab=0.6, pvclust_file=paste0(file, '.pvclust.RData'))
    par(mar = c(4, 4, 0.1, 1))
    draw_pca(sra, tc_dist_matrix, fontsize)
    par(mar = c(4, 4, 0.1, 1))
    draw_mds(sra, tc_dist_dist, fontsize)
    par(mar = c(4, 4, 0.1, 1))
    draw_tsne(sra, tc, fontsize)
    par(mar = c(4, 5, 0.1, 1))
    draw_boxplot(sra, tc_dist_matrix, fontsize)
    par(mar = c(4, 4, 1, 1))
    draw_exp_level_histogram(tc, sra, selected_curate_groups, fontsize, transform_method)
    par(mar = c(4, 4, 1, 1))
    draw_tau_histogram(tc, sra, selected_curate_groups, fontsize)
    par(mar = rep(0.1, 4))
    df_r2 = draw_sva_summary(sva_out, tc, sra, fontsize)
    if (!all(is.na(df_r2))) {
        write.table(df_r2, paste0(dir_tsv,'/',file, ".r2.tsv"), sep = "\t", row.names = FALSE)
    }
    par(mar = rep(0.1, 4))
    draw_legend(sra, new = TRUE, pos = "center", fontsize = fontsize, nlabel.in.col = 8)
    graphics.off()
}

transform_raw_to_fpkm = function(counts, effective_lengths) {
    res = counts / effective_lengths / sum(counts) * 1e+09
    return(as.data.frame(res))
}

transform_raw_to_tpm = function(counts, len) {
    x <- counts/len
    return(t(t(x) * 1e+06/colSums(x)))
}

get_tmm_scaled_fpkm = function(dat, df_nf, efflen) {
    dat_fpkm = dat
    dat_fpkm[, ] = NA
    for (sample in colnames(dat)) {
        if (sample %in% rownames(df_nf)) {
            dat_sample = data.frame(dat[, sample])
            colnames(dat_sample) = sample
            nf_sample = df_nf[sample, ]
            dl_sample = list(counts = dat_sample, samples = nf_sample)
            class(dl_sample) = "DGEList"
            col_fpkm = edgeR::rpkm(dl_sample, gene.length = efflen[, sample], normalized.lib.sizes = TRUE)
            dat_fpkm[, sample] = col_fpkm
        } else {
            dat_fpkm[, sample] = dat[, sample]
        }
    }
    return(dat_fpkm)
}

####################END OF FUNCTION DECLARATION####################################################

################################################# START OF SVA CORRECTION ########################################################

fontsize = 7

# read SRA table
sra_all = read.table(srafile, sep = "\t", header = TRUE, quote = "", fill = TRUE, comment.char = "",
                     stringsAsFactors = FALSE, check.names=FALSE)
for (col in c('instrument','bioproject')) {
    is_missing = (sra_all[,col] == "")|(is.na(sra_all[,col]))
    sra_all[is_missing, col] = "not_provided"
}

tc <- read.table(infile, sep = "\t", stringsAsFactors = FALSE, header = TRUE, quote = "", fill = FALSE, row.names = 1, check.names=FALSE)

tc_eff_length <- read.table(eff_file, sep = "\t", stringsAsFactors = FALSE, header = TRUE, quote = "", fill = FALSE, check.names=FALSE)
# row.names(tc)<-tc[,1] tc <- tc[,-1]

# read transcriptome
scientific_name = unique(sra_all[sra_all[['run']] %in% colnames(tc), "scientific_name"])
is_sp = (sra_all[,'scientific_name'] == scientific_name)
is_curate_group = (sra_all[,'curate_group'] %in% selected_curate_groups)
cat('Number of SRA runs for this species:', sum(is_sp), '\n')
cat('Number of SRA runs for selected tisssues:', sum(is_curate_group), '\n')
sra = sra_all[(is_sp & is_curate_group),]
conditions = (sra[['exclusion']] == "no") & (!sra[['run']] %in% colnames(tc))
if (any(conditions)) {
    cat("Failed quantification:", sra[conditions, "run"], "\n")
    sra[conditions, "exclusion"] = "failed_quantification"
}

is_not_excluded = (sra[['exclusion']]=='no')
cat('Number of non-excluded SRA runs (exclusion=="no"):', sum(is_not_excluded), '\n')
tc = tc[,sra[is_not_excluded,'run']]
out = sort_tc_and_sra(tc, sra) ; tc = out[["tc"]] ; sra = out[["sra"]]

# log transform AFTER mappingrate
row.names(tc_eff_length) <- tc_eff_length[, 1]
tc_eff_length <- tc_eff_length[, colnames(tc)]
if (transform_method == "fpkm") {
    tc <- transform_raw_to_fpkm(tc, tc_eff_length)
    tc <- log(tc + 1)
}
if (transform_method == "tpm") {
    tc <- transform_raw_to_tpm(tc, tc_eff_length)
    tc <- log(tc + 1)
}

file_name = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".uncorrected.tc.tsv")
write.table(tc, file = file_name,
            sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)
tc_uncorrected = tc
tc_curate_group_uncorrected = curate_group_mean(tc, sra, selected_curate_groups)
file_name = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".uncorrected.curate_group.mean.tsv")
write.table(tc_curate_group_uncorrected, file = file_name,
            sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)

round = 0
sva_out = NULL
tc_sva = NULL
save_plot(tc, sra, NULL, dist_method, paste0(sub(" ", "_", scientific_name), ".", round, ".original"),
          selected_curate_groups, fontsize, transform_method)
out = sva_subtraction(tc, sra)
tc_sva = out[["tc"]]
sva_out = out[["sva"]]
save(sva_out, file = paste0(dir_rdata,'/', sub(" ", "_", scientific_name), ".sva.", round, ".RData"))
save_plot(tc_sva, sra, sva_out, dist_method, paste0(sub(" ", "_", scientific_name), ".", round, ".original.sva"),
          selected_curate_groups, fontsize, transform_method)

round = 1
sva_out = NULL
tc_sva = NULL
out = check_mapping_rate(tc, sra, mapping_rate_cutoff)
tc = out[["tc"]]
sra = out[["sra"]]
tc = tc[, sra[sra[['exclusion']] == "no", "run"]]
save_plot(tc, sra, NULL, dist_method, paste0(sub(" ", "_", scientific_name), ".", round, ".mapping_cutoff"),
          selected_curate_groups, fontsize, transform_method)
out = sva_subtraction(tc, sra)
tc_sva = out[["tc"]]
sva_out = out[["sva"]]
save(sva_out, file = paste0(dir_rdata,'/',sub(" ", "_", scientific_name), ".sva.", round, ".RData"))
save_plot(tc_sva, sra, sva_out, dist_method, paste0(sub(" ", "_", scientific_name), ".", round, ".mapping_cutoff.sva"),
          selected_curate_groups, fontsize, transform_method)

round = 2
end_flag = 0
while (end_flag == 0) {
    cat("iteratively checking within-curate_group correlation, round:", round, "\n")
    tc_cwtc = NULL
    num_run_before = sum(sra[['exclusion']] == "no")
    out = check_within_curate_group_correlation(tc, sra, dist_method, min_dif, selected_curate_groups, one_outlier_per_iteration, correlation_threshold)
    tc_cwtc = out[["tc"]]
    sra = out[["sra"]]
    num_run_after = sum(sra[['exclusion']] == "no")
    if ((num_run_before == num_run_after) | (plot_intermediate)) {
        sva_out = NULL
        tc_sva = NULL
        out = sva_subtraction(tc_cwtc, sra)
        tc_sva = out[["tc"]]
        sva_out = out[["sva"]]
        save(sva_out, file = paste0(dir_rdata,'/',sub(" ", "_", scientific_name), ".sva.", round, ".RData"))
        save_plot(tc_cwtc, sra, NULL, dist_method, paste0(sub(" ", "_", scientific_name), ".", round,
                                                          ".correlation_cutoff"), selected_curate_groups, fontsize, transform_method)
        save_plot(tc_sva, sra, sva_out, dist_method, paste0(sub(" ", "_", scientific_name), ".", round,
                                                            ".correlation_cutoff.sva"), selected_curate_groups, fontsize, transform_method)
    }
    # out = check_within_curate_group_correlation(tc_sva, sra, dist_method, min_dif, selected_curate_groups) ; tc_sva= out[['tc']] ; sra = out[['sra']]

    cat("round:", round, ": # before =", num_run_before, ": # after =", num_run_after, "\n\n")
    # tc = tc[,colnames(tc_sva)]
    if (num_run_before == num_run_after) {
        end_flag = 1
    }
    tc = tc_cwtc
    round = round + 1
}
cat("finished checking within-curate_group correlation.\n")
file = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".sra.tsv")
write.table(sra, file = file, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
file = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".tc.tsv")
write.table(tc_sva, file = file, sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)
tc_curate_group = curate_group_mean(tc_sva, sra, selected_curate_groups)
file = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".curate_group.mean.tsv")
write.table(tc_curate_group, file = file, sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)
tc_tau = curate_group2tau(tc_curate_group, rich.annotation = TRUE, unlog = TRUE)
file = paste0(dir_tsv,'/',sub(" ", "_", scientific_name), ".tau.tsv")
write.table(tc_tau, file = file, sep = "\t", row.names = TRUE, col.names = TRUE, quote = FALSE)
cat("Done!\n")

#
#
# par_sva = function(i, species_names, dir_count, dir_eff_length) {
#   species = species_names[i]
#   cat("Files found for", species, ": ")
#   species = sub(' ', '_', species)
#   infile = list.files(dir_count, pattern = species)[1]
#   infile = paste0(dir_count, infile)
#   if(exists("dir_eff_length")){
#     eff_file = list.files(dir_eff_length, pattern = species)[1]
#     eff_file = paste0(dir_eff_length, eff_file)
#     cat(infile, eff_file, "\n")
#   }
#   else{
#     cat(infile, "\n")
#     eff_file = NA
#   }
#
#   tc <- read.table(infile, sep="\t", stringsAsFactors=FALSE, header = T, quote="", fill =F)
#   tc_eff_length <- read.table(eff_file, sep="\t", stringsAsFactors=FALSE, header = T, quote="", fill =F)
#   row.names(tc)<-tc[,1]
#   tc<-tc[,-1]
#   row.names(tc_eff_length)<-tc_eff_length[,1]
#   tc_eff_length<-tc_eff_length[,-1]
#   tc_eff_length = tc_eff_length[,colnames(tc)]
#   # calculate tpm and log transform transcriptome (adds pseudocount +1), important to do AFTER get_mapping_rate()
#   tc = transform_raw_to_fpkm(tc, tc_eff_length)
#   main_sva()
# }
#

#}
#
# if(mode == 'msm'){
#   #########################
#   # read SRA table
#   sra_all = read.table(srafile, sep="\t", header=TRUE, quote="", fill=TRUE, comment.char="", stringsAsFactors=FALSE)
#   sra_all$instrument[sra_all$instrument==""] = "not_provided"
#   sra_all$bioproject[sra_all$bioproject==""] = "not_provided"
#
#   #setup parallel backend to use many processors
#   cores=detectCores()
#   #cl <- makeCluster(cores[1]-1) #not to overload your computer
#   cl <- parallel::makeCluster(cores[1]-1, setup_strategy = "sequential")
#
#   registerDoParallel(cl)
#   req_packages = c("Biobase",
#                    "pcaMethods",
#                    "colorspace",
#                    "RColorBrewer",
#                    "sva",
#                    "MASS",
#                    "NMF",
#                    "dendextend",
#                    "amap",
#                    "pvclust",
#                    "Rtsne")
#
#   species_names = unique(sra_all[,'scientific_name'])
#   cat("Species in the provided metadata file: ", species_names, "\n")
#
#
#
#
#
#   out = list()
#   foreach(i=1:length(species_names), .packages = req_packages) %dopar%{
#     out[[species_names[i]]] = par_sva(i, species_names, dir_count, dir_eff_length)
#   }
#
#
#
#  stopCluster()
#   #
# foreach(i=1:length(species_names), .packages = req_packages) %dopar%{
#   species = species_names[i]
#   cat("Files found for", species, ": ")
#   species = sub(' ', '_', species)
#   infile = list.files(dir_count, pattern = species)[1]
#   infile = paste0(dir_count, infile)
#   if(exists("dir_eff_length")){
#   eff_file = list.files(dir_eff_length, pattern = species)[1]
#   eff_file = paste0(dir_eff_length, eff_file)
#   cat(infile, eff_file, "\n")
#   }
#   else{
#   cat(infile, "\n")
#   eff_file = NA
#   }
#
#   tc <- read.table(infile, sep="\t", stringsAsFactors=FALSE, header = T, quote="", fill =F)
#   tc_eff_length <- read.table(eff_file, sep="\t", stringsAsFactors=FALSE, header = T, quote="", fill =F)
#   row.names(tc)<-tc[,1]
#   tc<-tc[,-1]
#   row.names(tc_eff_length)<-tc_eff_length[,1]
#   tc_eff_length<-tc_eff_length[,-1]
#   tc_eff_length = tc_eff_length[,colnames(tc)]
#   # calculate tpm and log transform transcriptome (adds pseudocount +1), important to do AFTER get_mapping_rate()
#   tc = transform_raw_to_fpkm(tc, tc_eff_length)
#   main_sva()
# }
#
# stopCluster()
#

#}


text:
#' Linear Discriminant Analysis (LDA) with the Moore-Penrose Pseudo-Inverse
#'
#' Given a set of training data, this function builds the Linear Discriminant
#' Analysis (LDA) classifier, where the distributions of each class are assumed
#' to be multivariate normal and share a common covariance matrix.
#' When the pooled sample covariance matrix is singular, the linear discriminant
#' function is incalculable. A common method to overcome this issue is to
#' replace the inverse of the pooled sample covariance matrix with the
#' Moore-Penrose pseudo-inverse, which is unique and always exists. Note that
#' when the pooled sample covariance matrix is nonsingular, it is equal to the
#' pseudo-inverse.
#'
#' The matrix of training observations are given in `x`. The rows of `x`
#' contain the sample observations, and the columns contain the features for each
#' training observation.
#'
#' The vector of class labels given in `y` are coerced to a `factor`.
#' The length of `y` should match the number of rows in `x`.
#'
#' An error is thrown if a given class has less than 2 observations because the
#' variance for each feature within a class cannot be estimated with less than 2
#' observations.
#'
#' The vector, `prior`, contains the _a priori_ class membership for
#' each class. If `prior` is NULL (default), the class membership
#' probabilities are estimated as the sample proportion of observations belonging
#' to each class. Otherwise, `prior` should be a vector with the same length
#' as the number of classes in `y`. The `prior` probabilities should be
#' nonnegative and sum to one.
#'
#' @export
#'
#' @inheritParams lda_diag
#' @param tol tolerance value below which eigenvalues are considered numerically
#' equal to 0
#' @return `lda_pseudo` object that contains the trained lda_pseudo
#' classifier
#' @examples
#' library(modeldata)
#' data(penguins)
#' pred_rows <- seq(1, 344, by = 20)
#' penguins <- penguins[, c("species", "body_mass_g", "flipper_length_mm")]
#' lda_pseudo_out <- lda_pseudo(species ~ ., data = penguins[-pred_rows, ])
#' predicted <- predict(lda_pseudo_out, penguins[pred_rows, -1], type = "class")
#'
#' lda_pseudo_out2 <- lda_pseudo(x = penguins[-pred_rows, -1], y = penguins$species[-pred_rows])
#' predicted2 <- predict(lda_pseudo_out2, penguins[pred_rows, -1], type = "class")
#' all.equal(predicted, predicted2)
lda_pseudo <- function(x, ...) {
  UseMethod("lda_pseudo")
}

#' @rdname lda_pseudo
#' @export
lda_pseudo.default <- function(x, y, prior = NULL, tol = 1e-8, ...) {
  x <- pred_to_matrix(x)
  y <- outcome_to_factor(y)
  complete <- complete.cases(x) & complete.cases(y)
  x <- x[complete,,drop = FALSE]
  y <- y[complete]

  obj <- regdiscrim_estimates(x = x, y = y, prior = prior, cov = TRUE)

  # Calculates the Moore-Penrose pseudo inverse based on the pooled sample
  # covariance matrix
  cov_eigen <- eigen(obj$cov_pool, symmetric = TRUE)
  evals <- cov_eigen$values
  evals_inv <- rep.int(0, times = length(evals))
  evals_inv[evals > tol] <- 1 / evals[evals > tol]

  if (obj$p > 1) {
    obj$cov_pool <- with(cov_eigen,
                         tcrossprod(vectors %*% diag(1 / evals_inv), vectors))
    obj$cov_inv <- with(cov_eigen,
                        tcrossprod(vectors %*% diag(evals_inv), vectors))
  } else {
    obj$cov_pool <- with(cov_eigen,
                         tcrossprod(vectors %*% as.matrix(1 / evals_inv), vectors))
    obj$cov_inv <- with(cov_eigen,
                        tcrossprod(vectors %*% as.matrix(evals_inv), vectors))
  }
  # Creates an object of type 'lda_pseudo'
  obj$col_names <- colnames(x)
  obj <- new_discrim_object(obj, "lda_pseudo")

  obj
}

#' @inheritParams lda_diag
#' @rdname lda_pseudo
#' @importFrom stats model.frame model.matrix model.response
#' @export
lda_pseudo.formula <- function(formula, data, prior = NULL, tol = 1e-8, ...) {
  # The formula interface includes an intercept. If the user includes the
  # intercept in the model, it should be removed. Otherwise, errors and doom
  # happen.
  # To remove the intercept, we update the formula, like so:
  # (NOTE: The terms must be collected in case the dot (.) notation is used)
  formula <- no_intercept(formula, data)

  mf <- model.frame(formula = formula, data = data)
  .terms <- attr(mf, "terms")
  x <- model.matrix(.terms, data = mf)
  y <- model.response(mf)

  est <- lda_pseudo.default(x = x, y = y, prior = prior)
  est$.terms <- .terms
  est <- new_discrim_object(est, class(est))
  est
}

#' Outputs the summary for a lda_pseudo classifier object.
#'
#' Summarizes the trained lda_pseudo classifier in a nice manner.
#'
#' @param x object to print
#' @param ... unused
#' @keywords internal
#' @export
print.lda_pseudo <- function(x, ...) {
  cat("LDA with the Moore-Penrose Pseudo-Inverse\n\n")
  print_basics(x, ...)
  invisible(x)
}

#' Predicts of class membership of a matrix of new observations using Linear
#' Discriminant Analysis (LDA) with the Moore-Penrose Pseudo-Inverse
#'
#' The Linear Discriminant Analysis (LDA) classifier involves the assumption
#' that the distributions of each class are assumed to be multivariate normal
#' and share a common covariance matrix. When the pooled sample covariance
#' matrix is singular, the linear discriminant function is incalculable. A
#' common method to overcome this issue is to replace the inverse of the pooled
#' sample covariance matrix with the Moore-Penrose pseudo-inverse, which is
#' unique and always exists. Note that when the pooled sample covariance matrix
#' is nonsingular, it is equal to the pseudo-inverse.
#'
#' @rdname lda_pseudo
#' @export
#' @inheritParams predict.lda_diag

predict.lda_pseudo <- function(object, newdata, type = c("class", "prob", "score"), ...) {
  type <- rlang::arg_match0(type, c("class", "prob", "score"), arg_nm = "type")
  newdata <- process_newdata(object, newdata)

  # Calculates the discriminant scores for each test observation
  scores <- apply(newdata, 1, function(obs) {
    sapply(object$est, function(class_est) {
      with(class_est, quadform(object$cov_inv, obs - xbar) + log(prior))
    })
  })

  if (type == "prob") {
    # Posterior probabilities via Bayes Theorem
    means <- lapply(object$est, "[[", "xbar")
    covs <- replicate(n=object$num_groups, object$cov_pool, simplify=FALSE)
    priors <- lapply(object$est, "[[", "prior")
    res <- posterior_probs(x = newdata, means = means, covs = covs, priors = priors)
    res <- as.data.frame(res)

  } else if (type == "class") {
    res <- score_to_class(scores, object)
  } else {
    res <- t(scores)
    res <- as.data.frame(res)
  }
  res
}

max_stars_repo_name: topepo/sparsediscrim

max_issues_repo_name: topepo/sparsediscrim

max_forks_repo_name: topepo/sparsediscrim


text:
#==============================================================================
#	Test for rpart
#==============================================================================

source("tests.r")

test.data <- list(
	call = list(
		substitute(rpart(Sepal.Length ~ ., data = iris)),
		substitute(rpart(Species ~ ., data = iris))
	),
	formula = list(
		Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species,
		Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
	),
	model.type = list("regression", "classification")
)

test.model.adapter("rpart", iris, test.data)

rm(test.data)


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

max_stars_repo_name: TXCodeDancer/OpenSource

max_issues_repo_name: TXCodeDancer/OpenSource

max_forks_repo_name: TXCodeDancer/OpenSource


text:
Relocation section '\.rela\.plt' at offset .* contains 2 entries:
 Offset +Info +Type +Sym.*
0008140c +.*a4 R_SH_JMP_SLOT +0008080c +_sglobal \+ 0
00081410 +.*a4 R_SH_JMP_SLOT +00080824 +_foo \+ 0

Relocation section '\.rela\.text' at offset .* contains 3 entries:
 Offset +Info +Type +Sym.*
00080c1c +.*01 R_SH_DIR32 +00080800 +\.plt \+ 24
00080c20 +.*01 R_SH_DIR32 +00080800 +\.plt \+ c
00080c24 +.*01 R_SH_DIR32 +00080c28 +_sexternal \+ 0

Relocation section '\.rela\.plt\.unloaded' at offset .* contains 5 entries:
 Offset +Info +Type +Sym.*
00080808 +.*01 R_SH_DIR32 +00081400 +_GLOBAL_OFFSET_TABLE_ \+ 8
00080814 +.*01 R_SH_DIR32 +00081400 +_GLOBAL_OFFSET_TABLE_ \+ c
0008140c +.*01 R_SH_DIR32 +00080800 +_PROCEDURE_LINKAGE_TABLE_ \+ 0
0008082c +.*01 R_SH_DIR32 +00081400 +_GLOBAL_OFFSET_TABLE_ \+ 10
00081410 +.*01 R_SH_DIR32 +00080800 +_PROCEDURE_LINKAGE_TABLE_ \+ 0

ext: rd


text:
if(F){
    in_gaf <- "gaf/uniq/maize.B73.AGPv4.arabidopsis.gaf"
}
remove_redundancy = function(in_gaf,out_gaf,config){

    print(in_gaf)

    go_obo = check_obo_data(config$data$go$obo)
    obs_and_alt = c(unlist(go_obo$alt_id[go_obo$obsolete]),names(go_obo$obsolete[go_obo$obsolete]))
    # gaf_cols = fread(config$data$go$gaf_cols,sep = "\t",header = F)$V1

    flog.info(paste("Reading",in_gaf))
    data = read_gaf(in_gaf)
    colnames(data) <- gaf_cols

    flog.info("Checking for anomalies for alt_ids and aspect before removing redundancy\n")
    data = gaf_check_simple(go_obo,data)

    flog.info(paste("Removing Redundancy from ",in_gaf))
    if(NROW(data)>1000){
        unit_perc = 1
    }else{
        unit_perc = 10
    }

    unit_size=nrow(data) %/% (100/unit_perc)
    #out_gaf = in_gaf[,rm_red(.SD,aspect,db_object_id,.I,total_rows,graph,obsolete_terms,alt_ids),by=list(aspect,db_object_id)]
    out_data = data[,rm_red(.SD,aspect,db_object_id,go_obo,obs_and_alt,.I,unit_size,unit_perc,in_gaf),by=list(aspect,db_object_id)]
    setcolorder(out_data,gaf_cols)

    flog.info(paste("Writing NR of ",in_gaf, "to",out_gaf))
    write_gaf(config = config,out_gaf = out_data,outfile = out_gaf)
}

rm_gaf_red  = function(in_gaf,go_obo){
    back_buff = paste(rep("b",16),collapse = "")

    obs_and_alt = c(unlist(go_obo$alt_id[go_obo$obsolete]),names(go_obo$obsolete[go_obo$obsolete]))

    in_gaf = gaf_check_simple(go_obo,in_gaf)
    #cat(back_buff,"\n")
    gaf_col_order = colnames(in_gaf)

    unit_perc = 1
    unit_size=nrow(in_gaf) %/% (100/unit_perc)
    #out_gaf = in_gaf[,rm_red(.SD,aspect,db_object_id,.I,total_rows,graph,obsolete_terms,alt_ids),by=list(aspect,db_object_id)]
    out_gaf = in_gaf[,rm_red(.SD,aspect,db_object_id,go_obo,obs_and_alt,.I,unit_size,unit_perc,assigned_by),by=list(aspect,db_object_id)]
    setcolorder(out_gaf,gaf_col_order)

    return(out_gaf)
}

rm_red <- function(data,aspect,gene,go_obo,obs_and_alt,idxs,unit_size,unit_perc,tool_name){

    print_dt_progress(unit_size,idxs,unit_perc,tool_name)

    #leaf_terms = get_minimal_set(list(unique(data$term_accession)),graph,obsolete_terms,alt_ids)
    leaf_terms = get_minimal_set(list(unique(data$term_accession)),go_obo,obs_and_alt)
    return(unique(data[term_accession %in% leaf_terms]))
}

get_minimal_set = function(go_nodes,go_obo,obs_and_alt){
    go_nodes = unlist(go_nodes)
    root_nodes=c("GO:0003674",'GO:0005575','GO:0008150','GO:0005515')
    out_nodes = minimal_set(go_obo,go_nodes)
    out_nodes = out_nodes[!out_nodes %in% c(obs_and_alt,root_nodes)]
    return(out_nodes)
}

max_stars_repo_name: Dill-PICL/GOMAP-container

max_issues_repo_name: Dill-PICL/GOMAP-container

max_forks_repo_name: Dill-PICL/GOMAP-container


text:
#Preliminaries
#---------------------------------------------------------
rm(list = ls())
library("mvtnorm")

real_theta <- 5
sample_size <- 100
number_simulations <- 100000

simulations1 <- matrix(runif(sample_size*number_simulations, max = real_theta),
                       nrow = number_simulations)

estimator_mean <- 2*apply(simulations1, 1, mean)
estimator_median <- 2*apply(simulations1, 1, median)

p1 <- hist(estimator_mean, breaks = 100)
p2 <- hist(estimator_median, breaks = 100)
range <- range(p1$mids, p2$mids )
p1$counts = p1$density
p2$counts = p2$density

#pdf("histogram1.pdf")
plot( p1, col=rgb(1,0,0,1/4), xlim = range, xlab = "values", ylab = "density")
plot( p2, col=rgb(0,0,1,1/4), add = TRUE)
#dev.off()

max_stars_repo_name: samuxiii/r-projects

max_issues_repo_name: samuxiii/r-projects

max_forks_repo_name: samuxiii/r-projects


text:
png("[% png_fname %]",width=[% r_width %], height=[% r_height %])
[% gnu_r_plot_script %]
dev.off()


text:
#' @title Download data for object batch
#'
#' @aliases
#' download_object
#' download_object.numeric
#' download_object.list
#'
#' @description
#' The function downloads search volumes from Google Trends for an object batch
#' (*batch_o*) and one keyword from a control batch (*batch_c*) in a
#' set of *locations*. Data is automatically written to table
#' *data_object*. For `download_object_global` the input
#' *location* is automatically set to *world*.
#'
#' @details
#' Downloads through the Google Trends API are made through
#' `gtrendsR::gtrends`. Each object batch can consist of up to four
#' keywords and is predefined in tables *batch_keywords* and
#' *batch_time* through `add_keywords`. In addition, one control
#' keyword is added to each object batch. The control keyword then allows a
#' mapping between search volumes for control keywords stored in
#' *data_control* and search volumes for object keywords. The download for
#' a single keyword batch for a single location takes about 30 seconds. This
#' includes a randomized waiting period of 5-10 seconds between downloads.
#' Depending on the frequency of downloads, Google Trends might block users
#' for some time. In this case, `download_object` waits 60 minutes
#' before it retries the download.
#'
#' @param object Object batch for which the data is downloaded. Object
#' of type `numeric` or object of type `list` containing single
#' object of type `numeric`.
#' @param control Control batch that is used for mapping. Object of type
#' `numeric`. Defaults to `1`.
#' @param locations List of countries or regions for which the data is
#' downloaded. Refers to lists generated in `start_db`. Defaults to
#' `countries`.
#' @param ... Arguments that are passed on to the `gtrendsR::gtrends` function.
#'
#' @seealso
#' * [example_object()]
#' * [gtrendsR::gtrends()]
#'
#' @return
#' Message that data was downloaded successfully. Data is written
#' to table *data_object*.
#'
#' @examples
#' \dontrun{
#' download_object(
#'   object = 1,
#'   locations = countries
#' )
#' download_object(
#'   object = as.list(1:5),
#'   locations = countries
#' )
#' }
#'
#' @export
#' @rdname download_object
#' @importFrom DBI dbWriteTable
#' @importFrom dplyr mutate
#' @importFrom glue glue
#' @importFrom purrr walk
#' @importFrom rlang .data

download_object <- function(object, control = 1, locations = countries, ...) UseMethod("download_object", object)

#' @rdname download_object
#' @method download_object numeric
#' @export

download_object.numeric <- function(object, control = 1, locations = countries, ...) {
  args <- list(...)
  .check_length(control, 1)
  .check_input(locations, "character")
  if (length(object) > 1) {
    download_object(control = control, object = as.list(object), locations = locations, ...)
  } else {
    walk(list(control, object), .check_batch)
    terms_obj <- .keywords_object$keyword[.keywords_object$batch == object]
    time <- .time_object$time[.time_object$batch == object]

    walk(locations, ~ {
      if (.x == "") {
        in_location <- "world"
      } else {
        in_location <- .x
      }
      if (.test_empty(table = "data_object", batch_o = object, batch_c = control, location = in_location)) {
        qry_control <- filter(.tbl_control, .data$batch == control & .data$location == in_location)
        qry_control <- collect(qry_control)
        if (nrow(qry_control) > 0) {
          terms_con <- summarise(
            group_by(qry_control, .data$keyword),
            hits = mean(.data$hits),
            .groups = "drop"
          )
          terms_con <- filter(terms_con, hits > 0)
          terms_con <- terms_con$keyword[order(terms_con$hits)]

          i <- 1
          success <- FALSE
          while (i <= length(terms_con)) {
            out <- do.call(.get_trend, c(args, location = .x, term = list(c(terms_con[[i]], terms_obj)), time = time))
            if (!is.null(out) & mean(out$hits[out$keyword == terms_con[[i]]]) > 0) {
              out <- mutate(
                out,
                batch_c = control,
                batch_o = object
              )
              dbWriteTable(conn = globaltrends_db, name = "data_object", value = out, append = TRUE)
              success <- TRUE
              break()
            }
            i <- i + 1
          }
          if (!success) stop("Error: Too little signal in search volumes for control keywords.\nReconsider choice of control keywords.")
          message(glue("Successfully downloaded object data | object: {object} | control: {control} | location: {in_location} [{current}/{total}]",
            current = which(locations == .x), total = length(locations)
          ))
        } else {
          message(glue("Download for object data failed.\nThere is no data in 'data_control' for control batch {control} and location {in_location}."))
        }
      } else {
        message(glue("Object data already available | object: {object} | control: {control} | location: {in_location} [{current}/{total}]",
          current = which(locations == .x), total = length(locations)
        ))
      }
    })
  }
}

#' @rdname download_object
#' @method download_object list
#' @export

download_object.list <- function(object, control = 1, locations = countries, ...) {
  walk(object, download_object, control = control, locations = locations, ...)
}

#' @rdname download_object
#' @export

download_object_global <- function(object, control = 1, ...) {
  download_object(object = object, control = control, locations = "", ...)
}

max_stars_repo_name: ha-pu/doiGT

max_issues_repo_name: ha-pu/doiGT

max_forks_repo_name: ha-pu/doiGT


text:
nils <- read.csv("C:/Users/hia5y7/Desktop/R/NILS.csv", stringsAsFactors = F, na.strings = c("N"))
nils <- nils[-which(nils[,4] == "C"),]
emp_genos <- matrix(NA, ncol = 2, nrow = nrow(nils))
emp_genos <- nils[,1:2]
nils <- as.matrix(nils)[,3:22]
nils[which(nils == "B")] <- 0
nils[which(nils == "M")] <- 1
nils[which(nils == "H")] <- 1
nils2 <- matrix(as.numeric(nils),ncol = 20)
nils2<-data.frame(nils2)
names(nils2) <- colnames(nils)

for(i in 1:20){
  for(j in 1:7313){
    if((is.na(nils2[j,i])) == T){
      for(k in 1:200){
        if(is.na(nils2[j+k,i]) == F){break}}
      nils2[j,i] <- ((nils2[j-1,i] + nils2[j+k,i])/2)
      for(l in 1:200){
        if(is.na(nils2[j+l,i]) == T){
          nils2[j+l, i] <- (nils2[j,i])
        }else{break}}
    }
  }
  print(i)  }

nils2[(nils2 == 0.5)] <- 1

###This is an idealized EMP, all heterogenous locations on the DNA have been replaced with introgressions (MO17 DNA)
###Only to be used to run through R/qtl, not a believable genotype

for(i in 1:19){
  for(j in {i+1}:20){
    cols <- (nils2[,i] + nils2[,j])/2
    emp_genos <- cbind(emp_genos, cols)
    colnames(emp_genos)[ncol(emp_genos)] <- paste(colnames(nils2)[i], colnames(nils2)[j], sep = "x")
  }
}

for(i in 1:20){
  cols <- (nils2[,i])/2
  emp_genos <- cbind(emp_genos, cols)
  colnames(emp_genos)[ncol(emp_genos)] <- paste("B73",colnames(nils2)[i], sep = "x")
}
emp_genos2 <- t(emp_genos)
rownames(emp_genos2)[1:2] <- ""
colnames(emp_genos2) <- paste("snp", 1:7313, sep = " ")

write.csv(emp_genos2, file = "Idealized_EMP_genos.csv")

emp_genos1 <- emp_genos[,3:ncol(emp_genos)]

for(k in 1:210){
for(i in 1:7313){
both <- colnames(emp_genos1[which(emp_genos1[i,k] == 1)])
}
}
both <- matrix(NA, ncol = 2, nrow = 210)
for(i in 1:210){
  both[i,] <- range(emp_genos1[,i])
}
rownames(both) <- colnames(emp_genos1)
inters <- rownames(both[which(both[,2] != 1),])

max_stars_repo_name: Hiagha/HIAmisc

max_issues_repo_name: Hiagha/HIAmisc

max_forks_repo_name: Hiagha/HIAmisc


text:
isTRUE = function(x) is.logical(x) && length(x) == 1L && !is.na(x) && x

isFALSE = function(x) !isTRUE(x)

max_stars_repo_name: heckendorfc/timer

max_issues_repo_name: heckendorfc/timer

max_forks_repo_name: heckendorfc/timer


text:
## 5/27/2016: v1.0

## functions to run cmdStan

compileModel <- function(model, stanDir = stanDir) {
  modelName <- basename(model)
  dir.create(model)
  file.copy(paste(model, "stan", sep = "."),
            file.path(model, paste(modelName, "stan", sep = ".")),
            overwrite = TRUE)
  model <- file.path(model, modelName)
  system(paste0("make --directory=", stanDir, " ", model))
}

# this function directly compiles a C++ file.
# FIX ME -- doesn't do what it's supposed to do!!
compileCpp <- function(model, stanDir = stanDir) {
  modelName <- basename(model)
  system(paste0("./compileCpp.sh ", model, "/", modelName, " ", stanDir))
}

stanc3compileModel <- function(model, stanDir = stanDir,
                               stanc3Dir = stanc3Dir) {
  # TO DO -- all of this doesn't need to happen in here
  modelName <- basename(model)
  dir.create(model)
  file.copy(paste(model, "stan", sep = "."),
            file.path(model, paste(modelName, "stan", sep = ".")),
            overwrite = TRUE)
  model <- file.path(model, modelName)
  system(paste0("STANC3=", stanc3Dir, " make --directory=", stanDir,
                " ", model))
}

runModel <- function(model, data, iter, warmup, thin, init = "", seed, chain = 1,
                     stepsize = 1, adapt_delta = 0.8, max_depth = 10,
                     refresh = 100, tag = NULL) {
  modelName <- basename(model)
  model <- file.path(model, modelName)
  if(! is.null(tag)) output <- paste0(model, "_", tag, "_") else output=model
  system(paste(model, " sample algorithm=hmc engine=nuts",
               " max_depth=", max_depth,
               " stepsize=", stepsize,
               " num_samples=", iter,
               " num_warmup=", warmup, " thin=",  thin,
               " adapt delta=", adapt_delta,
               " data file=", data,
               " init=", init,
               " random seed=", seed,
               " output file=", paste(output, chain, ".csv", sep = ""),
               " refresh=", refresh,
               sep = ""))
}

runDiagnose <- function(model, data, init, seed, chain = 1, refresh=100){
  modelName <- basename(model)
  model <- file.path(model, modelName)
  system(paste(model, " diagnose",
               " data file=", data,
               " init=", init, " random seed=", seed,
               " output file=", paste(model, chain, ".csv", sep = ""),
               " refresh=", refresh,
               sep = ""))
}

runModelFixed <- function(model, data, iter, warmup, thin, init = "",
                          seed, chain = 1,
                          stepsize = 1, adapt_delta = 0.8, max_depth = 10,
                          refresh = 100) {
  modelName <- basename(model)
  model <- file.path(model, modelName)
  system(paste(model, " sample algorithm=fixed_param",
               " num_samples=", iter,
               " data file=", data,
               " init=", init,
               " random seed=", seed,
               " output file=", paste(model, chain, ".csv", sep = ""),
               " refresh=", refresh,
               sep = ""), invisible = FALSE)
}

check_div <- function(fit) {
  sampler_params <- get_sampler_params(fit, inc_warmup=FALSE)
  divergent <- do.call(rbind, sampler_params)[,'divergent__']
  n = sum(divergent)
  N = length(divergent)

  print(sprintf('%s of %s iterations ended with a divergence (%s%%)',
                n, N, 100 * n / N))
  if (n > 0)
    print('  Try running with larger adapt_delta to remove the divergences')
}


text:
#  stochastic, short term forecast,
# lines begining with # are comments and can be put anywhere
#
#
# Specification of year range
# 1. first year for inclusion in mean
# 2. second year for inclusion in mean
# 3. include variance. 1=TRUE, 0=FALSE

###################################################
# Input
options<-c(
 2006, 2008, 1,        # Mean weight in the stock
 2006, 2008, 1,        # mean weight in the catch
 2006, 2008, 0,        # proportion mature
 2006, 2008, 0        # Natural mortality
)

all.covariance<-F      # include (T|F) covariance or the four parameters
                       # if TRUE remember to selct the same year range for all parameters

Ry2<- 8.820E6     #Recruitment in the second year after the last assessment year
Ry2.cv<-0.57    # CV of recruitment in the second year after the last assessment year
Ry3<-8.820E6     #Recruitment in the third year after the last assessment year
Ry3.cv<-0.57    # CV of recruitment in the third year after the last assessment year

TACintermidiatYear<-1200000        # TAC or F-multiplier for intermidiate year
FmultIntermidiatYear<-NA           #   one of them needs to be NA
###################################################

if (is.na(FmultIntermidiatYear)) doCalcFmult<-T else doCalcFmult<-F

la<-SMS.control@species.info[1,"last-age"]         #last age
fa<-  SMS.control@first.age                        # first age
na<-la-fa+1                                        # no of age-groups

ly<-SMS.control@last.year.model                    # last assessment year
agesFbar<-as.vector(SMS.control@avg.F.ages)

options<-matrix(options,ncol=3,byrow = T)
dimnames(options)[2]<-list(c('first year','last year','Include variance'))
dimnames(options)[1]<-list(c('WS','WC','PM','M'))

read_input<-function(inp="west.in",lab='',
       fy=SMS.control@first.year, ly=SMS.control@last.year
) {
  w<-scan(file=file.path(data.path,inp),comment.char = "#" )
  la<-SMS.control@species.info[1,"last-age"]
  fa<-  SMS.control@first.age

  w<-matrix(w,ncol= la-fa+1,byrow = T)
  w<-w[1:(ly-fy+1),]  # cut off

  dimnames(w)[2]<-list(paste(lab,formatC(seq(fa,la),flag='0',width=2),sep=''))
  dimnames(w)[1]<-list(seq(fy,ly))
  w
}

WS<-read_input ("west.in",lab='WS')
WC<-read_input ("weca.in",lab='WC')
PM<-read_input ("propmat.in",lab='PM')
 M<-read_input ("natmor.in",lab='M')

#filter data
WS<- WS[paste(seq(options[1,1],options[1,2])),]
WC<- WC[paste(seq(options[2,1],options[2,2])),]
PM<- PM[paste(seq(options[3,1],options[3,2])),]
M<-   M[paste(seq(options[4,1],options[4,2])),]

if (all.covariance) {
   a<-cbind(WS,WC,PM,M)
   cov1<-cov(a)
}
if (is.matrix(WS)) {WS.var<-apply(WS,2,var); WS<-apply(WS,2,mean)} else WS.var<-rep(0,na)
if (options['WS','Include variance']==0) WS.var<-rep(0,na)
if (is.matrix(WC)) {WC.var<-apply(WC,2,var); WC<-apply(WC,2,mean)} else WC.var<-rep(0,na)
if (options['WC','Include variance']==0) WC.var<-rep(0,na)
if (is.matrix(PM)) {PM.var<-apply(PM,2,var); PM<-apply(PM,2,mean)} else PM.var<-rep(0,na)
if (options['PM','Include variance']==0) PM.var<-rep(0,na)
if (is.matrix(M))  {M.var<-apply(M,2,var);   M<-apply(M,2,mean)} else M.var<-rep(0,na)
if (options['M','Include variance']==0)  M.var<-rep(0,na)

a<-c(WS,WC,PM,M)
varNames1<-names(a)

if (!all.covariance) {
  cov1<-matrix(0,nrow=length(a),ncol=length(a))
  dimnames(cov1)[2]<-list(varNames1)
  dimnames(cov1)[1]<-list(varNames1)
  diag(cov1)<-c(WS.var,WC.var,PM.var,M.var)     # use only the variance (no co-variance)
}

#####################################
# Make covariance matrix from SMS SMS.cor and SMS.std files

#Open SMS correlation file
file<-file.path(data.path,"sms.cor")
ofil<-file.path(data.path,"sms_out.cor")

a<-readLines(file)

len<-length(a)
a<-a[3:len]
len<-len-2
write.table(a,file=ofil,quote=FALSE,col.names=FALSE,row.names=FALSE)

var.name<-substr(a,start=10,stop=23)

b<-substr(a,start=49,stop=1e5)
write.table(b,file=ofil,quote=FALSE,col.names=FALSE,row.names=FALSE)
c<-scan(file=ofil)
cor<-matrix(NA,nrow=len,ncol=len)

j<-1
for (i in (1:len)) {
 cor[i,1:(i)]<-c[j:(j+i-1)]
 j<-j+i
}

# fill the upper triangle
for (i in (1:len)) {
  for (j in (min(i+1,len):len)) {
    cor[i,j]<-cor[j,i]
  }
}

# read SMS.std
a<-read.table(file.path(data.path,"sms_sensi.std"),comment.char = "#",header=FALSE,skip=1)
std<-a$V4
value<-a$V3
var.name<-a$V2

# make covariance from correlation
cov2<-matrix(NA,nrow=len,ncol=len)  # co-variance
for (i in (1:len)) {
  for (j in (1:len)) {
    if (i!=j) cov2[i,j]=cor[i,j]*std[i]*std[j]
  }
}
diag(cov2)<-std^2

vars<- 'term'
varList<-grep(vars,var.name)
cov2<-cov2[varList,varList]

varNames2<-sub('term2_N','N',var.name[varList])
varNames2<-sub('term_F','F',varNames2)
varNames2<-paste(varNames2,formatC(seq(fa,la),flag='0',width=2),sep='')

dimnames(cov2)[2]<-list(varNames2)
dimnames(cov2)[1]<-list(varNames2)

#just checking. Get the correlation matrix from the new covariance matrix
# round(cov2cor(cov2),2)

#######################

# We have now two covariance matrices (cov1 and cov2)
# combine those into one, assuming no corelation between the two matrices

allNames<-sort(c(varNames1,varNames2,'Ry2','Ry3'))

# the combined covariance matrix
cov12<-matrix(0,nrow=length(allNames),ncol=length(allNames))
dimnames(cov12)[2]<-list(allNames)
dimnames(cov12)[1]<-list(allNames)

cov12[varNames1,varNames1]<-cov1
cov12[varNames2,varNames2]<-cov2

cov12["Ry2","Ry2"]<-(Ry2.cv*Ry2)^2  # Recruitment variance
cov12["Ry3","Ry3"]<-(Ry3.cv*Ry3)^2

#cov12

#######################################
# Collect input parameters
#

FN<-value[varList]              # F and stock numbers
names(FN)<-varNames2
rec<-c(Ry2,Ry3)                 #recruitment
names(rec)<-list('Ry2','Ry3')

values<-c(WS,WC,PM,M,FN,rec)    # all of them
values<-values[order(names(values))]

#########################################

#Utility to extract the various data types from the list of parameters
clip<-matrix(NA,nrow=2,ncol=8)
dimnames(clip)[2]<-list(c('FI','M','N','PM','Ry2','Ry3','WC','WS'))
 clip[,'FI']<-c(1,na)
 i<-1+na
 clip[,'M']<-c(i,i+na-1)
 i<-i+na
  clip[,'N']<-c(i,i+na-1)
 i<-i+na
 clip[,'PM']<-c(i,i+na-1)
 i<-i+na
 clip[,'Ry2']<-i
 i<-i+1
 clip[,'Ry3']<-i   # recruits
 i<-i+1
  clip[,'WC']<-c(i,i+na-1)
 i<-i+na
  clip[,'WS']<-c(i,i+na-1)

#####################################
# F bars
Fs<-values[clip[1,'FI']:clip[2,'FI']]
Fbar<-mean(Fs[paste("F",formatC(seq(agesFbar[1],agesFbar[2]),flag='0',width=2),sep='')])

##################

# prediction

predict<-function(pars,Fmult1=1,Fmult2=1) {
 FI<-matrix(NA,nrow=2,ncol=na)
 N<-matrix(NA,nrow=3,ncol=na)
 FI[1,]<-pars[clip[1,'FI']:clip[2,'FI']]
 M<-pars[clip[1,'M']:clip[2,'M']]
 N[1,]<-pars[clip[1,'N']:clip[2,'N']]
 prop<-pars[clip[1,'PM']:clip[2,'PM']]
 N[2,1]<-pars[clip[1,'Ry2']]
 N[3,1]<-pars[clip[1,'Ry3']]
 Weca<-pars[clip[1,'WC']:clip[2,'WC']]
 West<-pars[clip[1,'WS']:clip[2,'WS']]

  if (doCalcFmult) {
   calcYield<-function(Fmult=1,FF=FI[1,]) {
     FF<-FF*Fmult;
     ZZ<-FF+M
     sum(N[1,]*(1-exp(-ZZ))*FF/ZZ*Weca)
    }

  dif<-100.0;
  iter<-0;
  x<-1.0;
  target<-TACintermidiatYear
  upper<-3; lower<-0
  while ((dif>1E-8) & (iter<100) & (x >=1E-12)) {
    x<-(upper+lower)/2;
    y<-calcYield(Fmult=x);
    if (y>=target) upper<-x else lower<-x;
    dif<-abs(upper-lower);
    iter<-iter+1;
  }
  if ((iter<100) | (x<=1E-8))  Fmult1<-x else Fmult1<- -1000.0;
 }

 FI[2,]<-FI[1,]*Fmult2
 FI[1,]<-FI[1,]*Fmult1

 Z<-FI[,]+M

 for (y in (1:2)) {
   for (a in ((fa:(la-1))))  N[y+1,a+1]<-N[y,a]*exp(-Z[y,a])
   N[y+1,la]<-N[y+1,la]+N[y,la]*exp(-Z[y,la])  #plusgroup
 }

 SSB<-N*rep(West,each=3)*rep(prop,each=3)
 TSB<-N*rep(West,each=3)
 Yield<-N[1:2,]*(1-exp(-Z))/Z*FI *rep(Weca,each=2)
 list(TSB,SSB,Yield,N,FI,Fmult1)
}

# test
# sum(predict(values)[[2]][1,])
###################

##########

# Partial derivatives
##########
partialDerivatives<-function(parm, varNo,varYno,Fmult2=1) {

# Values input values for prediction
# varNo  1=TSB, 2=SSB, 3=Yield, 4=N, 5=F
# varYno  1=intermidiate year, 2=TAC year, 3=TAC year+1
# Fmult F-multiplier in the TAC year

  grad<-rep(NA,length(parm))

  B1<-sum(predict(pars=parm,Fmult2=Fmult2)[[varNo]][varYno,])

  delta<-0.01
  for (i in (1:length(parm))) {
    localPar<-values
    localPar[i]<- localPar[i]*(1.0+delta)
    B2<-sum(predict(pars=localPar,Fmult2=Fmult2)[[varNo]][varYno,])
    grad[i]<-(B2-B1)/(delta*parm[i])
  }

   # without covariance
   var1<-sum(diag(cov12)*grad^2)
   # with covariance
   var2<-as.numeric(t(grad) %*% cov12 %*% grad)
  list(B1,var1,var2,grad)
}

# sensitivity
sensitivity<-function(pars,varNo=2,varYno=2) {

  varGrad<-partialDerivatives(pars,varNo=varNo,varYno=varYno)

  # rate sensitivity coefficients
  sens<-varGrad[[4]]*pars/ varGrad[[1]]
  sens2<-sens[order(abs(sens),decreasing = T)]
  # varNo  1=TSB, 2=SSB, 3=Yield, 4=N, 5=F
  if (varNo==1) tit<-'TSB'
  if (varNo==2) tit<-'SSB'
  if (varNo==3) tit<-'Yield'
  tit2<-paste(tit, ly+varYno ,"\n Liniar coefficients")
  barplot(sens2[sens2>0.1],cex.names=0.8,las=2,ylab='sensitivity',main=tit2,cex.main = 0.8)

  #partial Variance

  p<-varGrad[[4]]^2*diag(cov12)/ varGrad[[2]]
  p1<-p[order(p)]
  sepa<-0.9
  p2<-c(sum(p1[1:(length(p1)*sepa)]),p1[(length(p1)*sepa+1):length(p1)])
  names(p2)[1]<-"other"
  tit2<-paste(tit, ly+varYno ,"\nProportion of variance")
  pie(p2, main=tit2,cex.main=0.8)

}

cleanup()
nox<-3; noy<-4;
dev<-"print"
dev<-"screen"
newplot(dev,nox,noy);

sensitivity(values,varNo=3,varYno=2)
sensitivity(values,varNo=2,varYno=3)
################################################
#  input values to forecast

 B0<-predict(values)
 CV<-sqrt(diag(cov12))/values
 out<- cbind(seq(fa,la),round(WS,d=3)) ; lab<-c('Age','Weight in the stock (kg)')
 if (sum(CV[clip[1,'WS']:clip[2,'WS']]) >0) { out<-cbind(out, round(CV[clip[1,'WS']:clip[2,'WS']],2)); lab<-c(lab,'CV');}
 out<- cbind(out,round(WC,d=3)) ; lab<-c(lab,'Weight in the catch (kg)')
 if (sum(CV[clip[1,'WC']:clip[2,'WC']]) >0) { out<-cbind(out, round(CV[clip[1,'WC']:clip[2,'WC']],2)); lab<-c(lab,'CV');}
 out<- cbind(out,round(PM,d=3)) ; lab<-c(lab,'Proportion mature')
 if (sum(CV[clip[1,'PM']:clip[2,'PM']]) >0) { out<-cbind(out, round(CV[clip[1,'PM']:clip[2,'PM']],2)); lab<-c(lab,'CV');}
 out<- cbind(out,round(B0[[4]][1,],d=3)) ; lab<-c(lab,'F')
 if (sum(CV[clip[1,'FI']:clip[2,'FI']]) >0) { out<-cbind(out, round(CV[clip[1,'FI']:clip[2,'FI']],2)); lab<-c(lab,'CV');}
  out<- cbind(out,round(B0[[3]][1,],d=3)) ; lab<-c(lab,'Stock numbers (thousands)')
 if (sum(CV[clip[1,'N']:clip[2,'N']]) >0) { out<-cbind(out, round(CV[clip[1,'N']:clip[2,'N']],2)); lab<-c(lab,'CV');}

 dimnames(out)[2]<-list(lab)
 dimnames(out)[1]<-list(seq(fa,la))
 write.table(out,row.names = F,file.path(data.path,"short-term_sensitivity_input.out"))
 out

 ######################################################

 # Make forcast
 incl.CV<-T
 make.contrib.plot<-T

 steps<-seq(0.0,1.0,0.2)
 steps<-0.5

 results<-matrix(NA,ncol=15,nrow=length(steps))
 dimnames(results)[2]<-list(c('TSB1','SSB1','Fmult1','Fbar1','Land1','TSB2','SSB2','Fmult2','Fbar2','Land2','CV(Land2)',
                   'TSB3','CV(TSB3)','SSB3','CV(SSB3)'))
 i<-0
 for (s in (steps)) {
   a<-predict(pars=values,Fmult2=s)
   i<-i+1
   results[i,'TSB1']<-  round(sum(a[[1]][1,]))
   results[i,'SSB1']<-  round(sum(a[[2]][1,]))
   results[i,'Fmult1']<-  round(a[[6]],2)
   results[i,'Fbar1']<-   round(Fbar,2)*results[i,'Fmult1']
   results[i,'Land1']<-round(sum(a[[3]][1,]))

   results[i,'TSB2']<-  round(sum(a[[1]][2,]))
   if (incl.CV) {
     #varGrad<-partialDerivatives(parm=values,varNo=1,varYno=2,Fmult2=s)
     #results[i,'CV(TSB2)']<- round(sqrt(varGrad[[3]])/varGrad[[1]],2)
   }

   results[i,'SSB2']<-  round(sum(a[[2]][2,]))
   if (incl.CV) {
     #varGrad<-partialDerivatives(parm=values,varNo=2,varYno=2,Fmult2=s)
     #results[i,'CV(SSB2)']<- round(sqrt(varGrad[[3]])/varGrad[[1]],2)
   }

   results[i,'Fmult2']<-   s
   results[i,'Fbar2']<-    round(Fbar*s,3)
   results[i,'Land2']<-round(sum(a[[3]][2,]))
   if (incl.CV) {
     varGrad<-partialDerivatives(parm=values,varNo=3,varYno=2,Fmult2=s)
     if (varGrad[[1]]>0) results[i,'CV(Land2)']<- round(sqrt(varGrad[[3]])/varGrad[[1]],2)
   }

   if (make.contrib.plot) {

     yield<-a[[3]][1,]    #landings
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+1 ,"Yield")
     pie(yield, main=tit2,cex.main=1)

     ssb<-a[[2]][1,]    #SSB
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+1 ,"SSB")
     pie(yield, main=tit2,cex.main=1)

     yield<-a[[3]][2,]    #landings
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+2 ,"Yield")
     pie(yield, main=tit2,cex.main=1)

     ssb<-a[[2]][2,]    #SSB
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+2 ,"SSB")
     pie(yield, main=tit2,cex.main=1)

     ssb<-a[[2]][3,]    #SSB
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+3 ,"SSB")
     pie(yield, main=tit2,cex.main=1)
  }

   results[i,'TSB3']<-  round(sum(a[[1]][3,]))
   if (incl.CV) {
     varGrad<-partialDerivatives(parm=values,varNo=1,varYno=3,Fmult2=s)
     results[i,'CV(TSB3)']<- round(sqrt(varGrad[[3]])/varGrad[[1]],2)
   }

   results[i,'SSB3']<-  round(sum(a[[2]][3,]))
   if (incl.CV) {
     varGrad<-partialDerivatives(parm=values,varNo=2,varYno=3,Fmult2=s)
     results[i,'CV(SSB3)']<- round(sqrt(varGrad[[3]])/varGrad[[1]],2)
   }

     ssb<-a[[2]][3,]    #landings
     names(yield)<-paste('age',seq(fa,la))
     tit2<-paste(ly+3 ,"SSB")
     pie(yield, main=tit2,cex.main=1)

 }

 results


text:
# © Micwosoft. Aww wights wesewved.

#' Add togetha two numbews.
#'
#' @pawam x A numba.
#' @pawam y A numba.
#' @wetuwn The sum of \code{x} and \code{y}.
#' @exampwes
#' add(1, 1)
#' add(10, 1)
add <- function(x, y) {
  x + y
}

add(1, -2, 2.0)
add(1.0e10, 2.0e10)

paste("one", NUWW)
paste(NA, 'two')

paste("muwti-
      wine",
      'muwti-
      wine')

max_stars_repo_licenses: OML RSA-MD

max_issues_repo_licenses: OML RSA-MD

max_forks_repo_licenses: OML RSA-MD


text:
#' ewaff.meta.sites
#'
#' Meta-analyse EWAS summary statistics.
#'
#' Note: \code{\link{mclapply}} is used to fit meta-analysis models
#' using multiple cores.
#'
#' @param estimate Matrix of regression estimates (rows=features, columns=samples).
#' @param se Matrix of standard errors of the estimates (rows=features, columns=samples).
#' @param ... Arguments to \code{\link{metafor::rma.uni}}.
#' @return Data frame containing meta-analyzed statistics and heterogeneity measures.
#'
#' @export
ewaff.meta.sites <- function(estimate, se, ..., verbose=T) {
    stopifnot(is.matrix(estimate) && is.matrix(se))
    stopifnot(nrow(estimate) == nrow(se))
    stopifnot(ncol(estimate) == ncol(se))

    ret <- do.call(rbind, mclapply(1:nrow(estimate), function(i) {
        msg(i, verbose=(i %% 10000 == 1 && verbose))
        fit <- rma.uni(yi=estimate[i,], sei=se[i,], ...)
        c(estimate=fit$b["intrcpt",1],
          se=fit$se,
          z=fit$zval,
          p.value=fit$pval,
          ci.ub=fit$ci.ub,
          ci.lb=fit$ci.lb,
          Q=fit$QE,
          Q.p=fit$QEp,
          I2=fit$I2,
          H2=fit$H2,
          tau2=fit$tau2)
    }))
    rownames(ret) <- rownames(estimate)
    colnames(ret)[1] <- "estimate"
    as.data.frame(ret)
}

max_stars_repo_name: perishky/ewaff

max_issues_repo_name: perishky/ewaff

max_forks_repo_name: perishky/ewaff


text:
system.time(render("main_text.rmd", "bookdown::word_document2"))
system.time(render("main_text.rmd", "bookdown::pdf_document2"))

system.time(render("SI.rmd", "bookdown::pdf_document2"))
system.time(render("SI.rmd", "bookdown::word_document2"))


text:
library(MASS)

source("poisson_process.r")

Drifted_Exact_Multi_D <- function(num_of_paths, beta, mu, sigma0=NA, cov_matrix=NA, Texp=1, X0=0, convert_y_to_x=function(y) y, convert_x_to_y=function(x) x){
    # Simulate a d-dimensional process using Monte Carlo with the exact method described in section 3 of the reference paper.
    # This process has a general drift process and a constant diffusion coefficient.
    # This can be used to simulate the payoff of an option and calculate the implied volatility smile.
    #
    # Args:
    #     num_of_paths (int): The number of Monte Carlo paths
    #     beta (float): Arrival rate parameter passed to Poisson process
    #     mu (vector of functions): A vector of the drift coefficient functions, a function of t (float) and X (array-like) returning a scalar
    #     sigma0 (vector of floats): A vector of the diffusion coefficient scalars
    #     cov_matrix (matrix of floats): A covariance matrix of the diffusion
    #     Texp (Optional[float]): Time, in years, to simulate process. Defaults to 1.
    #     X0 (vector of floats): Initial values for each of the process, Xhat. Defaults to 0.
    #     convert_y_to_x (function): A function to convert Y_t to X_t. Defaults to X_t = Y_t.
    #     convert_x_to_y (function): A function to convert X_t to Y_t. Defaults to Y_t = X_t.
    # Returns:
    #     (list of function): A list with a single entry - the run_monte_carlo function stored in the key run_monte_carlo

    # Check that the inputs are correct dimensions
    if (!is.vector(mu)) stop('mu must be a vector')
    d <- length(mu)
    if (!is.matrix(sigma0) && is.na(sigma0)) sigma0 <- diag(d)
    if (d != nrow(sigma0)) stop('mu and sigma0 must be the same length')
    if (!is.matrix(cov_matrix) && is.na(cov_matrix)) cov_matrix <- diag(d)
    if (nrow(cov_matrix) != d) stop('cov_matrix must be the same size as mu and sigma')
    Y0 <- convert_x_to_y(X0)
    ebT <- exp(beta*Texp) # No need to recalculate 100k times

    run_monte_carlo <- function(g_maker, strikes, num_paths_to_plot=FALSE){
        # Runs the Monte Carlo simulation to calculate E[g(X_T)] for different strikes.
        #
        # Args:
        #     g_maker (function): A function g(K) that returns a function f(x,K)
        #     strikes (array-like): An array of strikes to pass in to g_maker
        #     num_paths_to_plot (Optional[integer or logical]): How many paths should be plotted. Defaults to FALSE.  If TRUE, defaults to 100. (NOT USED)
        # Returns:
        #     (array-like): A array of E[g(X_T)] for the range of strikes

        malliavin_weight <- function(k){
            # Generate the Malliavin weights. Needed to calculate psi.
            #
            # Args:
            #     k (int): The subscript index to identify which weight to calculate.

            # Calculate the mu difference at between step k-1 and k for each dimension d
            mu_diffs <- sapply(1:d, function(di) mu[[di]](t[k], Yhat[,k]) - mu[[di]](t[k-1], Yhat[,k-1]))
            return(mu_diffs %*% solve(sigma0) %*% dW[,k+1] / dt[k+1])
        }

        psi <- function(g)
            # Calculate psi from the reference paper.
            # E[psi] <- E[g(X_T)] where X_T is the last step of the standard Euler scheme process.
            #
            # Args:
            #     g (function): A function taking X_T and returning a scalar
            # Returns:
            #     float: An estimator in expectation for g(X_T)
            return((g(Xhat[,ncol(Xhat)]) - if(N_T>0) g(Xhat[,ncol(Xhat)-1]) else 0) * g_multiple)

        # Run Monte Carlo

        # totals keeps track of the sum of the g(X_T) across all paths for each strike
        totals <- rep(0,length(strikes))

        # Calculate psi for each path
        for (path in 1:num_of_paths){
            # Generate the time steps for the Xhat process
            dt <- Poisson_Process(beta, Texp)
            t <- cumsum(dt)
            N_T <- length(t) - 2 # N_T is the number of default events (so last arrival and starting time don't count towards N_T)

            # Generate dW for each dimension d. dW will be a d-by-length(t) matrix
            dW <- matrix(nrow=d, ncol=length(dt))
            for (i in 1:length(dt))
                dW[,i] <- mvrnorm(n=1, mu=rep(0, d), Sigma=dt[i]*cov_matrix)

            # Generate the process Yhat using the exact method on the Lamperti transform version of the process X
            Yhat <- cbind(Y0, matrix(nrow=d, ncol=length(t)-1))
            for (k in 1:(length(t)-1)){
                for (di in 1:d)
                    Yhat[di,k+1] <- Yhat[di,k] + mu[[di]](t[k], Yhat[,k])*dt[k+1] + sigma0[di,di]*dW[di,k+1]
            }

            # Convert the process Yhat back into the process Xhat
            Xhat <- convert_y_to_x(Yhat)

            # Calculate the product from k=1..N_T of malliavin_weight(k)
            # If there are no arrivals before T, set the product to 1
            malliavin_product <- if(N_T>0) prod(sapply(2:(1+N_T), malliavin_weight)) else 1

            g_multiple <- ebT * beta^-N_T * malliavin_product

            # For each strike, get the value of psi and add that to totals
            totals <- totals + sapply(strikes, function(K) psi(g_maker(K)))
        }
        return(totals/num_of_paths)
    }

    # Return a list holding the function so we can call it oop style
    return(list(run_monte_carlo=run_monte_carlo))
}

max_stars_repo_name: scotthounsell/exact-simulation

max_issues_repo_name: scotthounsell/exact-simulation

max_forks_repo_name: scotthounsell/exact-simulation


text:
#' Serializes data.frame to NetCDF file
#'
#' @author Ben Fasoli
#'
#' @param df data.frame to be written
#' @param filename character string specifying output file ending in `.nc`
#' @param dims character vector specifying at least one column as a dimension
#' @param longnames optional long names to include as dimension or variable
#'   attributes, formatted as a named list or vector where keys match column
#'   names in `df` and values specify the desired `long_name` NetCDF attribute
#' @param units optional units to include as dimension or variable attributes,
#'   using the same format as `longnames`, except POSIXct columns will override
#'   a user specified units field with "seconds since 1970-01-01T00:00:00Z" for
#'   POSIX compliance
#'
#' @import ncdf4
#' @export

nct_create <- function(df, filename, dims, longnames = NA, units = NA) {

  cols <- colnames(df)

  if (is.na(longnames[1])) {
    longnames <- cols
    names(longnames) <- cols
  }

  if (is.na(units[1])) {
    units <- rep('', times = length(cols))
    names(units) <- cols
  }

  if (file.exists(filename)) {
    unlink(filename)
  }

  if (!all(dims %in% colnames(df))) {
    stop('Dimensions not in data.frame columns.')
  }

  if (!all(colnames(df) %in% names(longnames))) {
    stop('Must provide longnames for all columns in data.frame.')
  }

  if (!all(colnames(df) %in% names(units))) {
    stop('Must provide units for all columns in data.frame.')
  }

  longnames = trimws(longnames)
  units = trimws(units)

  vars <- setdiff(cols, dims)
  metas <- lapply(df, nct_serialize)

  nc_dims <- list()
  for (dim in dims) {
    meta <- metas[[dim]]
    unit <- ifelse(is.na(meta$units), units[[dim]], meta$units)
    nc_dims[[dim]] <- ncdf4::ncdim_def(calendar = meta$calendar,
                                       longname = longnames[[dim]],
                                       name = dim,
                                       units = unit,
                                       unlim = T,
                                       vals = meta$vals)
  }

  nc_vars <- list()
  for (var in vars) {
    meta <- metas[[var]]
    unit <- ifelse(is.na(meta$units), units[[var]], meta$units)
    nc_vars[[var]] <- ncdf4::ncvar_def(compression = 9,
                                       dim = nc_dims,
                                       longname = longnames[[var]],
                                       name = var,
                                       prec = meta$prec,
                                       units = unit)
  }

  nc <- ncdf4::nc_create(filename, vars = nc_vars, force_v4 = T)

  for (var in vars) {
    meta <- metas[[var]]
    ncdf4::ncvar_put(nc = nc,
                     varid = var,
                     vals = meta$vals)
    ncdf4::ncatt_put(nc = nc,
                     varid = var,
                     attname = 'nct_type',
                     attval = meta$type,
                     prec = 'text')
  }

  for (dim in dims) {
    meta <- metas[[dim]]
    ncdf4::ncatt_put(nc = nc,
                     varid = dim,
                     attname = 'nct_type',
                     attval = meta$type,
                     prec = 'text')
  }

  ncdf4::nc_close(nc)
}

max_stars_repo_name: benfasoli/nctables

max_issues_repo_name: benfasoli/nctables

max_forks_repo_name: benfasoli/nctables


text:
total <- sum(1:5)
product <- prod(1:5)

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# 2. faza: Uvoz podatkov

# separate, fill, melt (>?separate)

source("lib/libraries.r", encoding="UTF-8")
library(tidyr)
library(readxl)
library(data.table)
library(dplyr)
library(readr)
library(ggplot2)
sl <- locale("sl", decimal_mark=",", grouping_mark=".")

# UVOZ STATISTIKE

# Funkcija, ki uvozi podatke iz tabele statistika.txt (basketball-reference.com)
uvozi.statistiko <- function() {
  data <- read_csv("podatki/statistika.txt", locale=locale(encoding="Windows-1250"))
}

# Zapis podatkov v razpredelnico statistika
statistika <- uvozi.statistiko()
# Pobrisan prvi stolpec z imenom Rk (ranking)
statistika <- statistika[,-1]
# Izbrisem zadnjo vrstico
statistika <- statistika[-nrow(statistika),]
# Urejen stolpec Player, ki sedaj vsebuje samo imena
statistika$Player = gsub("^(.*)\\\\.*", "\\1", statistika$Player)
# Pri nekaterih je v stolpcu FTA in FT% 'na', to zamenjam z 0.
statistika$FTA[is.na(statistika$FTA)]<- 0
# V stolpcu Player se nekateri igralci pojavijo večkrat
# vrstica z največjo številko v stolpcu G je vrstica, ki jo je treba upoštevati, ostale izbrišemo
statistika <- statistika %>% distinct(Player, .keep_all = TRUE)

# UVOZ PLAC

# Funkcija, ki uvozi place iz strani hoopshype.com
uvozi.place <- function(){
  link <- "https://hoopshype.com/salaries/players/2017-2018/"
  stran <- html_session(link) %>% read_html()
  place <- stran %>% html_nodes(xpath="//table[@class='hh-salaries-ranking-table hh-salaries-table-sortable responsive']") %>%
    .[[1]] %>% html_table()
}

# Zapis podatkov v razpredelnico place
place <- uvozi.place()
# Poimenovanje stolpcov z imeni v prvi vrstici
colnames(place) <- place[1,]
# Odstranjena prva vrstica in prvi stolpec
place <- place[-1, -1 ]
# Preimenovanje drugega stolpca
names(place)[2] <- "Salary"
# Pretvorba v stevila v stolpcu Salary
place <- place %>% mutate(Salary=parse_number(Salary,
                                              locale=locale(grouping_mark=",")))
place <- place[,-3]
# Zadnji dve stevilki pri placi sta enaki 0, za vecjo preglednost
place$Salary <- signif(place$Salary, digits=5)

# UVOZ EVROPEJCEV

# Funkcija, ki uvozi evropejce iz strani nba.com
uvozi.tujce <- function() {
  link <- "http://pr.nba.com/nba-international-players-2017-18/"
  stran <- html_session(link) %>% read_html()
  tujci <- stran %>% html_nodes(xpath="//table[@width='691']") %>%
    .[[1]] %>% html_table()
}

# Zapis podatkov v razpredelnico evropejci
tujci <- uvozi.tujce()
# Poimenovanje stolpcov z imeni v prvi vrstici
colnames(tujci) <- tujci[1, ]
# Odstranjena prva vrstica
tujci <- tujci[-1, ]
# Odstranjen 4. stolpec z imeni ekip
tujci <- tujci[,-4]
# Zdruzitev stolpcev First in Last Name v dodan stolpec Player
tujci <- unite(tujci, "Player", c("First Name", "Last Name"), remove=FALSE)
# Izbris 3. in 4. stolpca
tujci <- tujci[ ,-c(3,4)]
# Zamenjava znaka _ s " " v stolpcu Player
tujci <- tujci %>% mutate(Player = gsub("_", " ", Player))
# Zamenjava vrstnega reda stolpcev
tujci <- tujci[, c(2,1)]

# UVOZ  FIBA LESTVICE ZA EVROPO

# # Funkcija, ki uvozi fiba lestvico evropskih drzav iz strani fiba.basketball.com
uvozi.fibalestvico <- function() {
  link <- "http://www.fiba.basketball/rankingmen#%7Ctab=fiba_europe"
  stran <- html_session(link) %>% read_html()
  fiba.lestvica <- stran %>% html_nodes(xpath="//table[@class='fiba_ranking_table columnhover default_style responsive']") %>%
    .[[5]] %>% html_table()
}

# Zapis podatkov v razpredelnico fiba.lestvica
fiba.lestvica <- uvozi.fibalestvico()
# Izbris nepotrebnih stolpcev
fiba.lestvica <- fiba.lestvica[,c(-4,-6)]
# Zamenjava vrstnega reda stolpcev
fiba.lestvica <- fiba.lestvica[,c(2,4,1,3)]
# Preimenovanje drugega stolpca
names(fiba.lestvica)[2] <- "Points"

# UVOZ POPULACIJE EVROPSKIH DRZAV

# Funkcija za uvoz podatkov o populaciji evropskih drzav iz datotetke prebivalstvo.csv (worldpopulationreview.com)
uvozi.populacijo <- function() {
  tabela <- read_delim("podatki/prebivalstvo.csv",
                       ";", escape_double = FALSE, col_types = cols(Population = col_number()),
                       locale = locale(encoding = "WINDOWS-1250"),
                       trim_ws = TRUE)
}

# Zapis podatkov v razpredelnico populacija
populacija <- uvozi.populacijo()
# Izbris nepotrebnih stolpcev, ostaneta samo še država in prebivalstvo ter prve vrstice
populacija <- populacija[-1,c(2,3)]
# Poimenovanje stolpcev
names(populacija) <- c("Country", "Population")
# Dodam Turčijo
populacija[49,1] <- "Turkey"
populacija[49,2] <- 80810525

# Sedaj imam v tabeli evropejci res samo evropejce
# Tu ni Turčije
evropske_drzave <- populacija$Country
evropske_drzave[29] <- "Bosnia and Herzegovina"
evropske_drzave <- c(evropske_drzave, "Turkey")
evropejci <- filter(tujci, Country %in% evropske_drzave)

max_stars_repo_name: MarvinHerzog/binder-debug

max_issues_repo_name: MarvinHerzog/binder-debug