text:
#########################################################
## script prepares data and exports files for analysis ##
#########################################################

# working dir
setwd("~/Dropbox/data/rollcall/ife_cg/ife-update/")

# some libraries
library(lubridate) # easier dates

# read data
d <- read.csv(file = "data/base_ife_eric_feb2021.csv", stringsAsFactors = FALSE)
str(d)

# fix logroll w vote at end not start
d$folio [d$folio==2299] <- 99999
d$folio [d$folio==2315] <- 2299
d$folio[d$folio==99999] <- 2315
# sort
d <- d[order(d$folio),]
# fix another logroll
d$dlogroll[d$folio==2623] <- 0
d$nconjunto[d$folio>=2503 & d$folio<=2514] <- 1

# drop recent consejeros without coded votes (feb 2021)
sel <- which(colnames(d) %in%  c("favela", "murayama", "ruiz", "sanmatin", "castillo", "galindo"))
d <- d[,-sel]

# turn NAs to vote=0
sel <- which(colnames(d) %in%  c("segob", "senpri", "senprd", "dippri", "dippan", "creel", "granados", "pinchetti", "pozas", "zertuche", "woldenberg", "barragan", "cantu", "cardenas", "lujambio", "merino", "molinar", "peschard", "zebadua", "luken", "rivera", "albo", "alcantar", "andrade", "glezluna", "latapi", "lopezflores", "morales", "sanchez", "ugalde", "banos", "nacif", "valdes", "elizondo", "figueroa", "guerrero", "cordova", "garcia", "marvan"))
tmp <- d[,sel]
tmp[is.na(tmp)] <- 0
d[,sel] <- tmp

# compute total ayes
d$tmp00 <- d$tmp01 <- d$tmp02 <- d$tmp03 <- d$tmp04 <- d$tmp05 <- d$tmp06 <- d$tmp07 <- d$tmp08 <- d$tmp09 <- d$tmp10 <- d$tmp11 <- d$tmp12 <- d$tmp13 <- d$tmp14 <- d$tmp15 <- d$tmp16 <- d$tmp17 <- d$tmp18 <- d$tmp19 <- d$tmp20 <- d$tmp21 <- d$tmp22 <- d$tmp23 <- d$tmp24 <- d$tmp25 <- d$tmp26 <- d$tmp27 <- d$tmp28 <- d$tmp29 <- d$tmp30 <- d$tmp31 <- d$tmp32 <- d$tmp33 <- d$tmp34 <- d$tmp35 <- d$tmp36 <- d$tmp37 <- d$tmp38 <- 0
#
d$tmp00[d$segob==1]     <- 1
d$tmp01[d$senpri==1]    <- 1
d$tmp02[d$senprd==1]    <- 1
d$tmp03[d$dippri==1]    <- 1
d$tmp04[d$dippan==1]    <- 1
d$tmp05[d$creel==1]     <- 1
d$tmp06[d$granados==1]  <- 1
d$tmp07[d$pinchetti==1] <- 1
d$tmp08[d$pozas==1]     <- 1
d$tmp09[d$zertuche==1]  <- 1
#
d$tmp10[d$woldenberg==1] <- 1
d$tmp11[d$barragan==1]   <- 1
d$tmp12[d$cantu==1]      <- 1
d$tmp13[d$cardenas==1]   <- 1
d$tmp14[d$lujambio==1]   <- 1
d$tmp15[d$merino==1]     <- 1
d$tmp16[d$molinar==1]    <- 1
d$tmp17[d$peschard==1]   <- 1
d$tmp18[d$zebadua==1]    <- 1
#
d$tmp19[d$luken==1]  <- 1
d$tmp20[d$rivera==1] <- 1
#
d$tmp21[d$albo==1]        <- 1
d$tmp22[d$alcantar==1]    <- 1
d$tmp23[d$andrade==1]     <- 1
d$tmp24[d$glezluna==1]    <- 1
d$tmp25[d$latapi==1]      <- 1
d$tmp26[d$lopezflores==1] <- 1
d$tmp27[d$morales==1]     <- 1
d$tmp28[d$sanchez==1]     <- 1
d$tmp29[d$ugalde==1]      <- 1
#
d$tmp30[d$banos==1]    <- 1
d$tmp31[d$nacif==1]    <- 1
d$tmp32[d$valdes==1]   <- 1
d$tmp33[d$elizondo==1] <- 1
d$tmp34[d$figueroa==1] <- 1
d$tmp35[d$guerrero==1] <- 1
#
d$tmp36[d$cordova==1] <- 1
d$tmp37[d$garcia==1]  <- 1
d$tmp38[d$marvan==1]  <- 1
#
d <- within(d, ayes <- tmp00 + tmp01 + tmp02 + tmp03 + tmp04 + tmp05 + tmp06 + tmp07 + tmp08 + tmp09 + tmp10 + tmp11 + tmp12 + tmp13 + tmp14 + tmp15 + tmp16 + tmp17 + tmp18 + tmp19 + tmp20 + tmp21 + tmp22 + tmp23 + tmp24 + tmp25 + tmp26 + tmp27 + tmp28 + tmp29 + tmp30 + tmp31 + tmp32 + tmp33 + tmp34 + tmp35 + tmp36 + tmp37 + tmp38)
#
# uncount logroll elements
tmp <- c(NA, d$nconjunto) # lag one obs (assumes numconj consecutive)
d$tmp <- tmp[-length(tmp)]
sel <- which(d$dlogroll==1 & d$nconjunto==d$tmp)
d$ayes[sel] <- NA
# and drop them (making uncounting redundant)
d <- d[-sel,]
# clean
sel <- grep("^tmp.*", colnames(d))
d <- d[,-sel] # clean

# drops 4 votes (with logrolls) where all 9 abstained
sel <- which(d$folio %in% c(5013, 5024, 5031, 5037, 2503))
d <- d[,-sel] # clean
rm(tmp,sel)

# total votes cast
d <- within(d, {
    segob1       <- as.numeric(      segob==1 | segob==2);
    senpri1      <- as.numeric(     senpri==1 | senpri==2);
    senprd1      <- as.numeric(     senprd==1 | senprd==2);
    dippri1      <- as.numeric(     dippri==1 | dippri==2);
    dippan1      <- as.numeric(     dippan==1 | dippan==2);
    creel1       <- as.numeric(      creel==1 | creel==2);
    granados1    <- as.numeric(   granados==1 | granados==2);
    pinchetti1   <- as.numeric(  pinchetti==1 | pinchetti==2);
    pozas1       <- as.numeric(      pozas==1 | pozas==2);
    woldenberg1  <- as.numeric( woldenberg==1 | woldenberg==2);
    zertuche1    <- as.numeric(   zertuche==1 | zertuche==2);
    barragan1    <- as.numeric(   barragan==1 | barragan==2);
    cardenas1    <- as.numeric(   cardenas==1 | cardenas==2);
    cantu1       <- as.numeric(      cantu==1 | cantu==2);
    lujambio1    <- as.numeric(   lujambio==1 | lujambio==2);
    merino1      <- as.numeric(     merino==1 | merino==2);
    molinar1     <- as.numeric(    molinar==1 | molinar==2);
    peschard1    <- as.numeric(   peschard==1 | peschard==2);
    zebadua1     <- as.numeric(    zebadua==1 | zebadua==2);
    luken1       <- as.numeric(      luken==1 | luken==2);
    rivera1      <- as.numeric(     rivera==1 | rivera==2);
    albo1        <- as.numeric(       albo==1 | albo==2);
    andrade1     <- as.numeric(    andrade==1 | andrade==2);
    alcantar1    <- as.numeric(   alcantar==1 | alcantar==2);
    glezluna1    <- as.numeric(   glezluna==1 | glezluna==2);
    latapi1      <- as.numeric(     latapi==1 | latapi==2);
    lopezflores1 <- as.numeric(lopezflores==1 | lopezflores==2);
    morales1     <- as.numeric(    morales==1 | morales==2);
    sanchez1     <- as.numeric(    sanchez==1 | sanchez==2);
    ugalde1      <- as.numeric(     ugalde==1 | ugalde==2);
    banos1       <- as.numeric(      banos==1 | banos==2);
    nacif1       <- as.numeric(      nacif==1 | nacif==2);
    valdes1      <- as.numeric(     valdes==1 | valdes==2);
    elizondo1    <- as.numeric(   elizondo==1 | elizondo==2);
    figueroa1    <- as.numeric(   figueroa==1 | figueroa==2);
    guerrero1    <- as.numeric(   guerrero==1 | guerrero==2);
    cordova1     <- as.numeric(    cordova==1 | cordova==2);
    garcia1      <- as.numeric(     garcia==1 | garcia==2);
    marvan1      <- as.numeric(     marvan==1 | marvan==2);
})
#
d$vtot <- NA
tmp <- d[d$term==1,]
tmp <- within(tmp, vtot <- segob1 + senpri1 + senprd1 + dippri1 + dippan1 + creel1 + granados1 + pinchetti1 + pozas1 + woldenberg1 + zertuche1)
d[d$term==1,] <- tmp
#
tmp <- d[d$term==2,]
tmp <- within(tmp, vtot <- woldenberg1 + peschard1 + lujambio1 + cardenas1 + cantu1 + merino1 + barragan1 + zebadua1 + molinar1)
d[d$term==2,] <- tmp
#
tmp <- d[d$term==3,]
tmp <- within(tmp, vtot <- woldenberg1 + peschard1 + lujambio1 + cardenas1 + cantu1 + merino1 + barragan1 + rivera1 + luken1)
d[d$term==3,] <- tmp
#
tmp <- d[d$term==4,]
tmp <- within(tmp, vtot <- ugalde1 + albo1 + andrade1 + alcantar1 + glezluna1 + latapi1 + lopezflores1 + morales1 + sanchez1)
d[d$term==4,] <- tmp
#
tmp <- d[d$term==5,]
tmp <- within(tmp, vtot <- albo1 + andrade1 + alcantar1 + glezluna1 + latapi1 + lopezflores1 + morales1 + sanchez1)
d[d$term==5,] <- tmp
#
tmp <- d[d$term==6,]
tmp <- within(tmp, vtot <- albo1 + andrade1 + alcantar1 + glezluna1 +         lopezflores1 +          sanchez1 + valdes1 + banos1 + nacif1)
d[d$term==6,] <- tmp
#
tmp <- d[d$term==7,]
tmp <- within(tmp, vtot <-       andrade1 + alcantar1 +                                         sanchez1 + valdes1 + banos1 + nacif1 + elizondo1 + figueroa1 + guerrero1)
d[d$term==7,] <- tmp
#
tmp <- d[d$term==8,]
tmp <- within(tmp, vtot <- valdes1 + banos1 + nacif1 + elizondo1 + figueroa1 + guerrero1)
d[d$term==8,] <- tmp
#
tmp <- d[d$term==9,]
tmp <- within(tmp, vtot <- valdes1 + banos1 + nacif1 + elizondo1 + figueroa1 + guerrero1 + marvan1 + cordova1 + garcia1)
d[d$term==9,] <- tmp
#
tmp <- d[d$term==10,]
tmp <- within(tmp, vtot <- valdes1 + banos1 + nacif1 + elizondo1 + figueroa1 + guerrero1 + marvan1 + cordova1)
d[d$term==10,] <- tmp
#
tmp <- d[d$term==11,]
tmp <- within(tmp, vtot <- banos1 + nacif1 + marvan1 + cordova1)
d[d$term==11,] <- tmp
#
d$absten <- NA
d$absten[d$term==1]           <- 11 - d$vtot[d$term==1]
d$absten[d$term>1 & d$term<5] <-  9 - d$vtot[d$term>1 & d$term<5]
d$absten[d$term==5]           <-  8 - d$vtot[d$term==5]
d$absten[d$term>5 & d$term<8] <-  9 - d$vtot[d$term>5 & d$term<8]
d$absten[d$term==8]           <-  6 - d$vtot[d$term==8]
d$absten[d$term==9]           <-  9 - d$vtot[d$term==9]
d$absten[d$term==10]          <-  8 - d$vtot[d$term==10]
d$absten[d$term==11]          <-  4 - d$vtot[d$term==11]
#
d <- within(d, nays <- vtot - ayes)
#
# verify
## table(d$ayes, useNA = "always")
## table(d$nays, useNA = "always")
## table(d$absten, useNA = "always")
#
# recompute result
d$tmp <- 1 - d$result # for comparison (was coded 0 pass 1 fail)
d <- within(d, result <- as.numeric(ayes - nays > 0))
# there were lots of inconsistent codings
## table(d$result, d$tmp, useNA = "always")
## sel <- which(d$result!=d$tmp)
## d[sel[2],]
# clean
sel <- grep(".*1$", colnames(d))
d <- d[,-sel]

# unanimous v contested votes
#
###############################################
## version up to US-Mex presentation in 2010 ##
###############################################
## replace unanime=1 if absten==0 & (ayes==9 | ayes==0) & term>1
## replace unanime=1 if absten==0 & (ayes==11 | ayes==0) & term==1
## replace unanime=1 if absten==1 & segob==. & (ayes==10 | ayes==0) & term==1
## replace unanime=1 if absten==9 & term>1
## replace unanime=1 if absten==0 & (ayes==8 | ayes==0) & term==5
## replace unanime=1 if absten==0 & (ayes==6 | ayes==0) & term==8
d$dunan.old <- 0
#with(d, table(ayes, nays, absten))
d$dunan.old[d$absten==0 & (d$ayes==0 | d$nays==0)]  <- 1
#
#############################################################
## versión que excluye ausentes de la cuenta de unanimidad ##
#############################################################
# falta para carpizo
d <- within(d, {
    barragangone <-    as.numeric(barragan==4    | barragan==5);
    cantugone <-       as.numeric(cantu==4       | cantu==5);
    cardenasgone <-    as.numeric(cardenas==4    | cardenas==5);
    molinargone <-     as.numeric(molinar==4     | molinar==5);
    lujambiogone <-    as.numeric(lujambio==4    | lujambio==5);
    merinogone <-      as.numeric(merino==4      | merino==5);
    peschardgone <-    as.numeric(peschard==4    | peschard==5);
    woldenberggone <-  as.numeric(woldenberg==4  | woldenberg==5);
    zebaduagone <-     as.numeric(zebadua==4     | zebadua==5);
    lukengone <-       as.numeric(luken==4       | luken==5);
    riveragone <-      as.numeric(rivera==4      | rivera==5);
    albogone <-        as.numeric(albo==4        | albo==5);
    alcantargone <-    as.numeric(alcantar==4    | alcantar==5);
    andradegone <-     as.numeric(andrade==4     | andrade==5);
    glezlunagone <-    as.numeric(glezluna==4    | glezluna==5);
    latapigone <-      as.numeric(latapi==4      | latapi==5);
    lopezfloresgone <- as.numeric(lopezflores==4 | lopezflores==5);
    moralesgone <-     as.numeric(morales==4     | morales==5);
    sanchezgone <-     as.numeric(sanchez==4     | sanchez==5);
    ugaldegone <-      as.numeric(ugalde==4      | ugalde==5);
    banosgone <-       as.numeric(banos==4       | banos==5);
    nacifgone <-       as.numeric(nacif==4       | nacif==5);
    valdesgone <-      as.numeric(valdes==4      | valdes==5);
    elizondogone <-    as.numeric(elizondo==4    | elizondo==5);
    figueroagone <-    as.numeric(figueroa==4    | figueroa==5);
    guerrerogone <-    as.numeric(guerrero==4    | guerrero==5);
    cordovagone <-     as.numeric(cordova==4     | cordova==5);
    garciagone <-      as.numeric(garcia==4      | garcia==5);
    marvangone <-      as.numeric(marvan==4      | marvan==5);
})
d <- within(d, noshow <- woldenberggone + molinargone + lujambiogone + peschardgone + merinogone + cardenasgone + barragangone + cantugone + zebaduagone + lukengone + riveragone + albogone + glezlunagone + sanchezgone + moralesgone + ugaldegone + latapigone + andradegone + lopezfloresgone + alcantargone + valdesgone + banosgone + nacifgone + elizondogone + figueroagone + guerrerogone + marvangone + cordovagone + garciagone)
# check
#table(d$noshow)
# clean
sel <- grep(".*gone$", colnames(d))
d <- d[,-sel]
d$tmp <- NULL
#
# noshows must be deducted from abstentions
d$absten <- d$absten - d$noshow
#
# "A vote qualifies as contested when, ignoring absences, at least one councilor voted contrary to the rest or abstained" ifedyn06
d$dunan <- 0
d$dunan[d$absten==0 & (d$ayes==0 | d$nays==0)]  <- 1
#
########################################################
## COULD GET OTHER UNANIMITY MEASURE FROM ife.do HERE ##
########################################################
#
# inspect
table(factor(d$dunan, labels = c("contested","not")), factor(d$dunan.old, labels = c("old-contested","not")), useNA = "always")
table(factor(d$dunan, labels = c("contested","not")), d$term, useNA = "ifany")

# sort agg vote count columns
sel <- which(colnames(d) %in% c("vtot","ayes","nays","absten","noshow"))
tmp <- d[,c("vtot","ayes","nays","absten","noshow")]
d[,sel] <- tmp
d <- d[,-sel]
d <- cbind(d, tmp)

# add date
d <- within(d, date <- ymd(yr*10000+mo*100+dy))
# add quarter/semester in calendar year (ojo: no longer our october-to-march version)
d <- within(d, {
    qtr <- year(date) + quarter(date)/10;
    sem <- year(date) + semester(date)/10;
})

#################
## export data ##
#################
#info.cols <- c("folio", "date", "qtr", "sem", "term")
info.cols <- c("folio", "date", "yr", "mo", "dy", "qtr", "sem", "term", "dunan")
#
sel.r <- which(d$term %in% c(2,3))
sel.c <- c("woldenberg", "barragan", "cantu", "cardenas", "lujambio", "merino", "molinar", "peschard", "zebadua", "rivera", "luken", info.cols)
tmp <- d[sel.r, sel.c]
write.csv(tmp, file = "data/v23.csv", row.names = FALSE)
#
## sel.r <- which(d$term %in% 4:8)
## c("ugalde", "albo", "andrade", "alcantar", "glezluna", "latapi", "lopezflores", "morales", "sanchez", "valdes", "banos", "nacif", "elizondo", "figueroa", "guerrero", info.cols)
## tmp <- d[sel.r, sel.c]
## write.csv(tmp, file = "data/v45678.csv", row.names = FALSE)
## #
## sel.r <- which(d$term %in% 4:9)
## c("ugalde", "albo", "andrade", "alcantar", "glezluna", "latapi", "lopezflores", "morales", "sanchez", "valdes", "banos", "nacif", "elizondo", "figueroa", "guerrero", "marvan", "cordova", "garcia", info.cols)
## tmp <- d[sel.r, sel.c]
## write.csv(tmp, file = "data/v456789.csv", row.names = FALSE)
#
sel.r <- which(d$term %in% 4:11)
sel.c <- c("ugalde", "albo", "andrade", "alcantar", "glezluna", "latapi", "lopezflores", "morales", "sanchez", "valdes", "banos", "nacif", "elizondo", "figueroa", "guerrero", "cordova", "garcia", "marvan", info.cols)
tmp <- d[sel.r, sel.c]
write.csv(tmp, file = "data/v456789ab.csv", row.names = FALSE)

# summarize contested votes
with(d[d$dunan==0,], plot(as.factor(year(date)+quarter(date)/10), main = "N monthly contested votes"))
with(d[d$dunan==0,], plot(as.factor(date), main = "N contested votes by session"))

cuts <- c(
    ymd("19961031"), # 1  to 2
    ymd("20001211"), # 2  to 3
    ymd("20031105"), # 3  to 4
    ymd("20071217"), # 4  to 5
    ymd("20080215"), # 5  to 6
    ymd("20080821"), # 6  to 7
    ymd("20101031"), # 7  to 8
    ymd("20111215"), # 8  to 9
    ymd("20130220"), # 9  to 10
    ymd("20131031")  # 10 to 11
)

max_stars_repo_name: grosasballina/ife-update

max_issues_repo_name: grosasballina/ife-update

max_forks_repo_name: grosasballina/ife-update


text:
#!/usr/bin/env Rscript
library(data.table)
library(tidyverse)
library(stringr)
library(stringi)
library(tidygraph)
library(zoo)
library(igraph)
library(parallel)

args <- commandArgs(TRUE)

pfam_cnames <- read_tsv("data/DBs/pfam_files/Pfam-31_names_mod_01122019.tsv", col_names = TRUE)
annot_kept  <- fread(paste0("gzip -dc ", args[1]), colClasses = c(V1 = "character", V8="character"),
                     stringsAsFactors = F, header = F) %>%
  dplyr::select(V1,V2,V3,V5,V6,V7,V8) %>%
  setNames(c("cl_name","repres","rep_annot","memb","memb_annot","clan","partial")) %>%
  as_tibble()

rename_pfam <- function(X){
  paste0(plyr::mapvalues(as.vector(str_split(X, "\\|", simplify = TRUE)), from = pfam_cnames$pfam, to = pfam_cnames$cname, warn_missing = FALSE),
         collapse = "|")
}
expand_repeats <- function(X, rep){
  lapply(X, function(X){
    if (is.element(X, rep$from)) {
      rep(X, rep %>% filter(from == X) %>% .$n + 1)
    }else{
      X
    }
  }) %>% unlist %>% paste(collapse = "|")
}

encode_repeats <- function(X, rep){
  lapply(X, function(X){
    if (is.element(X, rep$from)) {
      rep %>% filter(from == X) %>% .$n + 1
    }else{
      1
    }
  }) %>% unlist %>% paste(collapse = "|")
}

get_cluster_da_red <- function(X) {
  df_multi_process <- NULL
  df_mono_process <- NULL
  # Conditions
  # 1. First we check if we have complete ORFs
  # 2. We always will try to use the complete data, if not we will use the partial

  #X <- cls_mu[["1281746"]]
  df <- X %>% filter(!is.na(memb_annot))
  memb_annot <- df$memb_annot
  df <- as.data.table(df)
  df <- df[, memb_annot_reduced := rename_pfam(memb_annot),by=seq_len(nrow(df))] %>% as.tibble() %>%
    mutate(number_multi=stri_count_fixed(memb_annot_reduced,"|"))
  df %>% group_by(memb_annot_reduced,partial) %>% count() %>% View()
  if(max(df$number_multi)>1){
  df_rep <- df %>% filter(number_multi==max(number_multi)) %>% select(memb_annot_reduced) %>% distinct() %>%
    mutate(numb_rep=sapply(strsplit(memb_annot_reduced,"\\|"),function(x){length(unique(x))}))
    if(max(df$number_multi)>2 && df_rep$numb_rep == 1){
      df_rep <- df_rep %>% mutate(type_rep=sapply(strsplit(memb_annot_reduced,"\\|"), function(x){unique(x)}))
      df <- df %>% mutate(memb_annot_reduced=ifelse(grepl(df_rep$type_rep,memb_annot_reduced),df_rep$memb_annot_reduced,memb_annot_reduced)) %>%
        mutate(partial=ifelse(grepl(df_rep$type_rep,memb_annot_reduced),"00",partial))
    }
  }
  df_multi <- df %>% filter(number_multi > 0)

  df_multi_complete <- df_multi %>% filter(partial == "00")
  df_multi_partial <- df_multi %>% filter(partial != "00")

  df_multi_complete_n <- nrow(df_multi_complete)
  df_multi_partial_n <- nrow(df_multi_partial)

  df_mono <- df %>% mutate(number_multi = stri_count_fixed(memb_annot_reduced, "|")) %>% filter(number_multi == 0)
  df_mono_complete <- df_mono %>% filter(partial == "00")
  df_mono_partial <- df_mono %>% filter(partial != "00")

  df_mono_complete_n <- nrow(df_mono_complete)
  df_mono_partial_n <- nrow(df_mono_partial)

  # can we use the complete data:
  if (df_mono_complete_n < 1) {
    df_mono_process <- df_mono_partial %>% select(memb_annot_reduced) %>%
      group_by(memb_annot_reduced) %>% count() %>% rename(archit = memb_annot_reduced) %>% ungroup()
  }else{
    df_mono_process <- df_mono_complete %>% select(memb_annot_reduced) %>%
      group_by(memb_annot_reduced) %>% count() %>% rename(archit = memb_annot_reduced) %>% ungroup()
  }

# can we use the complete data:
  if (df_multi_complete_n < 1) {
    if (df_multi_partial_n > 0){
      df_multi_process <- df_multi_partial
    }
  }else{
    df_multi_process <- df_multi_complete
  }

  if (!is.null(df_multi_process)){
    c <- map(df_multi_process %>% .$memb_annot_reduced %>% strsplit("\\|") %>% unique(),
             function(x) {
               g <- rollapply(data = x, 2, by=1, c) %>%
                 graph_from_data_frame(directed = TRUE)

             })

    c <- graph.union(c) %>% as_tbl_graph() %>% mutate(degree = centrality_degree(loops = TRUE))

    archits <- do.call("rbind",lapply(
      df_multi_process %>% .$memb_annot_reduced %>% strsplit("\\|"),
      function(Y) {
        Y_g <- rollapply(data = Y, 2, by=1, c)  %>% graph_from_data_frame(directed = TRUE)
        g1 <- graph.intersection(c, Y_g) %>% as_tbl_graph() %>% activate(nodes) %>% filter(!(node_is_isolated()))
        #induced_subgraph(c,vids = Y,)#c %>% filter(name %in% Y)
        if (is_connected(g1)){
          # Identify loops
          E(g1)$is_loop <- which_loop(g1)
          # Identify repeats
          repeats <- g1 %>%
            igraph::as_data_frame() %>%
            as_tibble() %>%
            filter(is_loop == TRUE) %>%
            group_by(from) %>%
            count()
          arch <- topo_sort(g1) %>% as_ids()
          exp <- expand_repeats(arch, repeats)
          enc <- encode_repeats(arch, repeats)
          da <- data_frame(archit = paste0(topo_sort(g1) %>% as_ids(),collapse = "|"),
          exp_rep = exp,
          enc_rep = enc)
        }
      })) %>%
      as_tibble() %>% group_by(archit) %>% add_count(sort = TRUE) %>% distinct() %>%
      rename(n_in_path = n) %>% ungroup() %>%
      dplyr::top_n(n = 1, wt = n_in_path) %>% select(-n_in_path)
  }
  #path seen more time in the big cluster (c) add colum with number of reduced ex: AAA|P.... 2|1 etc...original..original one
  if(archits$exp_rep != ""){
  archits <- archits %>%
    mutate(class = "multi") %>%
    bind_rows(df_mono_process %>% select(-n) %>% mutate(class = "mono")) %>%
    #group_by(archit) %>%
    #summarise(n = sum(n)) %>%
    bind_cols(map_df(.$exp_rep, function(X) {
      d <- df %>% filter(memb_annot_reduced == X)
      tibble(complete = d %>% filter(partial == "00") %>% nrow,
             partial = d %>% filter(partial != "00") %>% nrow)
    }
    )) %>%
    arrange(desc(complete), desc(partial))
  }else{
    archits <- archits %>%
      mutate(class = "multi") %>%
      bind_rows(df_mono_process %>% select(-n) %>% mutate(class = "mono")) %>%
      #group_by(archit) %>%
      #summarise(n = sum(n)) %>%
      bind_cols(map_df(.$archit, function(X) {
        d <- df %>% filter(memb_annot_reduced == X)
        tibble(complete = d %>% filter(partial == "00") %>% nrow,
               partial = d %>% filter(partial != "00") %>% nrow)
      }
      )) %>%
      arrange(desc(complete), desc(partial))
  }
  list(g = c, domains = c %>% as_tibble() %>% .$name, n_comps = components(c)$no, archit = archits %>% ungroup(),
       mono_complete = df_mono_complete_n, mono_partial = df_mono_partial_n,
       multi_complete = df_multi_complete_n, multi_partial = df_multi_partial_n)
}

get_cluster_da_or <- function(X) {
  df_multi_process <- NULL
  df_mono_process <- NULL
  # Conditions
  # 1. First we check if we have complete ORFs
  # 2. We always will try to use the complete data, if not we will use the partial

  # X <- "10000404"
  df <- X %>% filter(!is.na(memb_annot))
  memb_annot <- df$memb_annot
  df <- df %>% rowwise() %>% mutate(memb_annot_reduced = rename_pfam(memb_annot)) %>% ungroup() %>%
    mutate(number_multi=stri_count_fixed(memb_annot,"|"))
  df_rep <- df %>% filter(number_multi==max(number_multi)) %>% select(memb_annot) %>% distinct() %>%
    mutate(numb_rep=sapply(strsplit(memb_annot,"\\|"),function(x){length(unique(x))}),type_rep=sapply(strsplit(memb_annot,"\\|"), unique))
  if(max(df$number_multi)>2 && df_rep$numb_rep == 1){
    df <- df %>% mutate(memb_annot=ifelse(grepl(df_rep$type_rep,memb_annot),df_rep$memb_annot,memb_annot))
   }
  df_multi <- df %>% mutate(number_multi = stri_count_fixed(memb_annot, "|")) %>% filter(number_multi > 0)
  df_multi_complete <- df_multi %>% filter(partial == "00")
  df_multi_partial <- df_multi %>% filter(partial != "00")

  df_multi_complete_n <- nrow(df_multi_complete)
  df_multi_partial_n <- nrow(df_multi_partial)

  df_mono <- df %>% mutate(number_multi = stri_count_fixed(memb_annot, "|")) %>% filter(number_multi == 0)
  df_mono_complete <- df_mono %>% filter(partial == "00")
  df_mono_partial <- df_mono %>% filter(partial != "00")

  df_mono_complete_n <- nrow(df_mono_complete)
  df_mono_partial_n <- nrow(df_mono_partial)

  # can we use the complete data:
  if (df_mono_complete_n < 1) {
    df_mono_process <- df_mono_partial %>% select(memb_annot) %>%
      group_by(memb_annot) %>% count() %>% rename(archit = memb_annot) %>% ungroup()
  }else{
    df_mono_process <- df_mono_complete %>% select(memb_annot) %>%
      group_by(memb_annot) %>% count() %>% rename(archit = memb_annot) %>% ungroup()
  }

  # can we use the complete data:
  if (df_multi_complete_n < 1) {
    if (df_multi_partial_n > 0){
      df_multi_process <- df_multi_partial
    }
  }else{
    df_multi_process <- df_multi_complete
  }

  if (!is.null(df_multi_process)){
    c <- map(df_multi_process %>% .$memb_annot %>% strsplit("\\|") %>% unique(),
             function(x) {
               g <- rollapply(data = x, 2, by=1, c) %>%
                 graph_from_data_frame(directed = TRUE)

             })

    c <- graph.union(c) %>% as_tbl_graph() %>% mutate(degree = centrality_degree(loops = TRUE))

    archits_or <- do.call("rbind",lapply(
      df_multi_process %>% .$memb_annot %>% strsplit("\\|"),
      function(Y) {
        g1 <- c %>% filter(name %in% Y)
        if (is_connected(g1)){
          arch <- paste0(topo_sort(g1) %>% as_ids(), collapse = "|")
          if(!grepl("\\|", arch)){
            dg <- g1 %>% activate(nodes) %>% as.tibble() %>% filter(name %in% arch) %>% .$degree
            if(dg>1){
              paste(rep(arch, dg+1), collapse = "|")
            }else{
              arch
            }
          }else{
            arch
          }
        }
      })) %>%
      as_tibble() %>% group_by(V1) %>% count(sort = TRUE) %>%
      rename(archit = V1, n_in_path = n) %>%
      top_n(n = 1, wt = n_in_path) %>% select(-n_in_path)
  }
  # archits <- df_multi_process %>% as_tibble() %>% group_by(cl_name,memb_annot) %>%
  #   count(sort = TRUE) %>% rename(archit=memb_annot, n_in_path = n) %>% ungroup() %>%
  #   dplyr::top_n(1, wt = n_in_path) %>% select(-n_in_path)
  archits <- archits_or %>%
    mutate(class = "multi") %>%
    bind_rows(df_mono_process %>% select(-n) %>% mutate(class = "mono")) %>%
    #group_by(archit) %>%
    #summarise(n = sum(n)) %>%
    bind_cols(map_df(.$archit, function(X) {
      d <- df %>% filter(memb_annot == X)
      tibble(complete = d %>% filter(partial == "00") %>% nrow,
             partial = d %>% filter(partial != "00") %>% nrow)
    }
    )) %>%
    arrange(desc(complete), desc(partial)) #%>% dplyr::slice(1)

  list(g = c, domains = c %>% as_tibble() %>% .$name, n_comps = components(c)$no, archit = archits %>% ungroup(),
        mono_complete = df_mono_complete_n, mono_partial = df_mono_partial_n,
        multi_complete = df_multi_complete_n, multi_partial = df_multi_partial_n)
}

get_cluster_da_mono <- function(X) {
  df_multi_process <- NULL
  df_mono_process <- NULL
  # Conditions
  # 1. First we check if we have complete ORFs
  # 2. We always will try to use the complete data, if not we will use the partial

  # X <- "10000404"
  df <- X %>% filter(!is.na(memb_annot))
  memb_annot <- df$memb_annot
  df <- df %>% rowwise() %>% mutate(memb_annot_reduced = rename_pfam(memb_annot)) %>% ungroup()
  if(length(unique(df$memb_annot))==1){
    homog = TRUE
  } else {
    homog = FALSE
  }
  df_mono <- df
  df_mono_complete <- df_mono %>% filter(partial == "00")
  df_mono_partial <- df_mono %>% filter(partial != "00")

  df_mono_complete_n <- nrow(df_mono_complete)
  df_mono_partial_n <- nrow(df_mono_partial)

  # can we use the complete data:
  if (df_mono_complete_n < 1) {
    df_mono_process <- df_mono_partial %>% select(memb_annot) %>%
      group_by(memb_annot) %>% count() %>% rename(archit = memb_annot) %>% ungroup()
  }else{
    df_mono_process <- df_mono_complete %>% select(memb_annot) %>%
      group_by(memb_annot) %>% count() %>% rename(archit = memb_annot) %>% ungroup()
  }

  archits <- df_mono_process %>% select(-n) %>% mutate(class = "mono") %>%
    #group_by(archit) %>%
    #summarise(n = sum(n)) %>%
    bind_cols(map_df(.$archit, function(X) {
      d <- df %>% filter(memb_annot == X)
      tibble(complete = d %>% filter(partial == "00") %>% nrow,
             partial = d %>% filter(partial != "00") %>% nrow)
    }
    )) %>% mutate(homog=homog) %>%
    arrange(desc(complete), desc(partial))

  list(archit = archits %>% ungroup(),
       mono_complete = df_mono_complete_n, mono_partial = df_mono_partial_n)
}

# we get the ones that are multidomain
h <- annot_kept %>% filter(grepl("\\|", memb_annot)) %>% mutate(number_multi = stri_count_fixed(memb_annot, "|"))
multi <- h %>% select(cl_name,number_multi) %>% distinct() %>% left_join(annot_kept)
# we create a list out of the different groups
cls_mu <- multi %>% split(.$cl_name)

## Reduced Pfam names
da_results_red <- mclapply(cls_mu_left, get_cluster_da_red, mc.cores = 20)
names(da_results_red) <- names(da_results_red)
# we select the first archtecture entry, it has been sorted the number of complete and partial oRFs
da_best_red <- map_df(da_results_red, function(X) X$archit %>% dplyr::slice(1), .id = "cl_name")

## Original Pfam names
da_results_or <- mclapply(cls_mu, get_cluster_da_or, mc.cores = 18)
names(da_results_or) <- names(da_results_or)
# we select the first archtecture entry, it has been sorted the number of complete and partial oRFs
da_best_or <- map_df(da_results_or, function(X) X$archit %>% dplyr::slice(1), .id = "cl_name")
da_best <- da_best_red1 %>% select(cl_name,archit) %>% left_join(da_best_or %>% rename(original=archit))

# we get the ones that are monodomain
multi_cl_name <- multi %>% select(cl_name) %>% distinct
m <- annot_kept %>% filter(!cl_name %in% multi_cl_name$cl_name)
# we create a list out of the different groups
cls_mo <- m %>% split(.$cl_name)
da_results_m1 <- mclapply(cls_mo[1:200000], get_cluster_da_mono, mc.cores = 18)
names(da_results_m1) <- names(da_results_m1)
da_best_m1 <- map_df(da_results_m1, function(X) X$archit %>% slice(1), .id = "cl_name")

da_results_m2 <- mclapply(cls_mo[200001:400000], get_cluster_da_mono, mc.cores = 18)
names(da_results_m2) <- names(da_results_m2)
da_best_m2 <- map_df(da_results_m2, function(X) X$archit %>% slice(1), .id = "cl_name")

da_results_m3 <- mclapply(cls_mo[400001:600000], get_cluster_da_mono, mc.cores = 18)
names(da_results_m3) <- names(da_results_m3)
da_best_m3 <- map_df(da_results_m3, function(X) X$archit %>% slice(1), .id = "cl_name")

da_results_m4 <- mclapply(cls_mo[600001:800000], get_cluster_da_mono, mc.cores = 18)
names(da_results_m4) <- names(da_results_m4)
da_best_m4 <- map_df(da_results_m4, function(X) X$archit %>% slice(1), .id = "cl_name")

da_results_m5 <- mclapply(cls_mo[800001:855345], get_cluster_da_mono, mc.cores = 18)
names(da_results_m5) <- names(da_results_m5)
da_best_m5 <- map_df(da_results_m5, function(X) X$archit %>% slice(1), .id = "cl_name")

# we select the first archtecture entry, it has been sorted the number of complete and partial oRFs
da_best_m_or <- rbind(da_best_m1,da_best_m2, da_best_m3, da_best_m4, da_best_m5)
da_best_m_or <- da_best_m_or[-6]
da_best_m <- da_best_m_red %>% select(cl_name,archit) %>% left_join(da_best_m_or %>% rename(original=archit))
#da_best_m <- da_best_m %>% filter(!cl_name %in% da_orig$cl_name)

da_best_all <- rbind(da_best_m,da_best)

fix_or <- function(X){
  df <- X %>% filter(!is.na(memb_annot))
  memb_annot <- df$memb_annot
  arch <- da_best_all %>% filter(cl_name %in% unique(df$cl_name)) %>% .$exp_rep
  n_multi <- stri_count_fixed(arch,"|")
  df <- df %>% rowwise() %>% mutate(memb_annot_reduced = rename_pfam(memb_annot)) %>% ungroup() %>%
    mutate(number_multi=stri_count_fixed(memb_annot,"|"))
  df %>% group_by(memb_annot) %>%
    filter(grepl(arch,memb_annot) & number_multi==n_multi) %>%
    mutate(complete=length(which(partial == "00")),partial1=length(which(partial != "00"))) %>%
    select(cl_name,memb_annot,complete,partial1) %>% rename(original=memb_annot,partial=partial1) %>% distinct() %>% ungroup() %>%
    arrange(desc(complete),desc(partial)) %>% dplyr::slice(1) %>% mutate(archit=arch)
}

da_res_or_mo <- mclapply(cls_mo, fix_or, mc.cores = 16)
da_res_or_mo <- plyr::ldply(da_res_or_mo,as.data.frame)
da_res_or_mo <- da_res_or_mo[-1]

da_res_or_mu <- mclapply(cls_mu_left, fix_or, mc.cores = 16)
da_res_or_mu <- plyr::ldply(da_res_or_mu,as.data.frame)
da_res_or_mu <- da_res_or_mu[-1]

da_res_mono <- da_best_all %>% filter(cl_name %in% names(cls_mo)) %>%
  left_join(da_res_or_mo %>% select(cl_name,archit,original), by="cl_name")
da_res_multi <- da_best_all %>% filter(cl_name %in% names(cls_mu)) %>%
  left_join(da_res_or_mu %>% select(cl_name,archit,original), by="cl_name")
da_best_all <- rbind(da_best_mono,da_best_multi)

#Define DUFs and PFs clusters (A cluster is consider annotated to DUFs if 100% of its annotated members are DUFs)
#Combine results and define DUFs and PFs clusters:
# A cluster is consider annotated to DUFs if the original consensus DA is 100% DUF
da_best_all <- rbind(da_best_multi,da_best_mono %>% mutate(exp_rep=NA,enc_rep=NA)) %>%
  separate_rows(original,sep="\\|") %>% group_by(cl_name) %>%
  mutate(type=ifelse(all(!grepl('^DUF',original)),"PF","DUF")) %>%
  group_by(cl_name,type) %>% mutate(original=paste(original,collapse = "|")) %>% distinct()
write.table(da_best_all, paste(args[2],"cluster_pfam_domain_architect.tsv",sep="/"), col.names = T, row.names = F, quote = F, sep = "\t")

#write.table(da_best_all %>% filter(type=="DUF"), paste(args[2],"kept_DUF.tsv",sep="/"), col.names = T, row.names = F, quote = F, sep = "\t")
#write.table(da_best_all %>% filter(type=="PF"), paste(args[2],"kept_PF.tsv",sep="/"), col.names = T, row.names = F, quote = F, sep = "\t")

ext: r


text:
title = "Fac (time)"

masterdata <- read.table("master-times-Fac", header=T)
recdata  <- read.table("recursive-times-Fac", header=T)
cpsdata <- read.table("cps-times-Fac", header=T)
origdata <- read.table("orig-times-Fac", header=T)

all <-  c(masterdata$usr+masterdata$sys, recdata$usr+recdata$sys, cpsdata$usr+cpsdata$sys, origdata$usr+origdata$sys)
ylim <- c(min(all), max(all))

start <- function(fr,col) {
    plot  (fr$n,fr$usr+fr$sys, col=col, pch=20, main=title, xlab = "n", ylab = "Time (s)", ylim=ylim)
    lines (fr$n,fr$usr+fr$sys, col=col, pch=20)
}

draw <- function(fr,col) {
    points (fr$n,fr$usr+fr$sys, col=col, pch=20)
    lines  (fr$n,fr$usr+fr$sys, col=col, pch=20)
}

mastercol="gold"
reccol="darkolivegreen3"
cpscol="blue"
origcol="darkmagenta"

start  (recdata,reccol)
draw (origdata, origcol)
draw  (masterdata,mastercol)
draw  (cpsdata, cpscol)

legend(x="topleft", inset=.05, legend=c("master","recursive", "CPS", "Original CEK machine"), col=c(mastercol, reccol, cpscol, origcol), lty=1, lw=2)

max_stars_repo_name: AriFordsham/plutus

max_issues_repo_name: AriFordsham/plutus

max_forks_repo_name: AriFordsham/plutus


text:
#!/usr/bin/env Rscript

args = commandArgs(trailingOnly=TRUE)

if(length(args) < 2){
  stop("Usage: lefse_prepare_input.r <relASVTable.tsv> <Metadata.tsv>")
}

rel.table <- args[1]
meta.table <- args[2]

OUT="lefse_ready_table.tsv"

#read in relative asv table

tab <-read.table(file = rel.table, sep = '\t', header = FALSE, comment.char = "", skip = 1)

#get column order of ASV table
col_order <- unlist(unname(tab[1,-1]))

#read in metadata file

meta <-read.table(file = meta.table, sep = '\t', header = FALSE)
meta.t <- as.data.frame(t(meta))[c(2,1),]
colnames(meta.t) <- meta.t[2,]
name <- colnames(meta.t)[1]

#re-order metadata columns as in ASV table order
meta.t <- meta.t[,c(name,col_order)]
colnames(meta.t) <- names(tab)

#combine metadata and ASV table
tab.tsv <- rbind(meta.t[1,],tab)

#convert ";" into "|" and remove "|" at the end of taxonomy names

df <- unname(data.frame(sapply(tab.tsv, gsub, pattern = ";", replacement= "|")))
col1 <- strsplit(as.character(df[,1]),split='|',fixed=TRUE)
col1.na_rm <- lapply(col1, function(z){ z[!is.na(z) & z != ""]})
col1.new <- lapply(col1.na_rm, paste, collapse = "|")
col1.df <- unname(data.frame(matrix(unlist(col1.new), nrow=length(col1.new), byrow=TRUE),stringsAsFactors=FALSE))

df[,1] <- col1.df

#write
print (paste("write",OUT))
write.table(df, file = OUT, row.names=FALSE, sep="\t", quote = FALSE)

max_stars_repo_name: PI-Bioinfo/ampliseqPlus

max_issues_repo_name: PI-Bioinfo/ampliseqPlus

max_forks_repo_name: PI-Bioinfo/ampliseqPlus


text:
	loadfunctions = function(
    projectname,
    directorypattern=NULL,
    functionname=NULL,
    RcodeDirectory="R",
    keydirectories=c("r", "\\.r", "\\_r", "rfunctions", "\\.rfunctions", "\\_rfunctions" ),
    toignore = c("retired", "_archive", "archive", "orphan", "request", "example" ),
    filepattern="\\.r$",
    directory=NULL ) {

    # used to load local functions conveniently
    # sequence slightly important ... modify with care
    # the  '\\' used in the above are escape sequences used with regex

    fs = .Platform$file.sep
    keydirectories = paste( "\\", fs, keydirectories, "\\", fs, sep="") # make sure only match directories and not file names
    filestosource = NULL

    for (pn in projectname ) {

      projectdirectory = project.codedirectory( directory, pn )

      for (searchdirectory in c( file.path( projectdirectory, RcodeDirectory ), projectdirectory ) ) {  # first try in RcodeDirectory and then the project if not found in first pass

        projectfiles = NULL
        projectfiles = list.files( path=searchdirectory, pattern=filepattern, full.names=T, recursive=T,  ignore.case=T, include.dirs=F )

        # filter on directory pattern
        if (!is.null(directorypattern)) {
          keep = grep ( directorypattern, projectfiles, ignore.case =T )
          if (length(keep) == 0) {
            # not found, try an approximate match
            keep = agrep( directorypattern, projectfiles, ignore.case =T )
          }
          if (length(keep)>0) projectfiles = projectfiles[ keep ]
        }

        # functionname filter
        if (!is.null(functionname)) {
          keep = grep ( functionname, projectfiles, ignore.case =T )
          if (length(keep) == 0) {
            # if no match try approximate matching
            keep = agrep( functionname, projectfiles, ignore.case =T )
          }
          if (length(keep)>0) {
            for ( nm in keep ) source( file= projectfiles[nm] )
            return( projectfiles[keep] )  # this breaks out of the loops as it is a single file-search and source call
          }
        }

        # remove archived functions, etc. that are to be "ignored"
        rem  = NULL
        for (i in toignore) {
          rem0 = grep ( i, projectfiles,  ignore.case =T )
          if (length( rem0)>0 ) rem = c( rem, rem0)
        }
        if ( length(rem)>0 ) {
          rem = unique(rem)
          projectfiles = projectfiles[-rem]
        }

        # filter on key directories
        # determine correct directory .. first look for a set of key directories, if not then scan all files
        md = NULL
        for ( kd in keydirectories) {
          md = c( md, grep ( kd, projectfiles, ignore.case =T ) )
          md = unique( md )
        }
        if (length(md) > 0 ) {
          projectfiles = projectfiles[ md ]
          filestosource = c( filestosource, projectfiles )
          break()  # no need to continue with the loop
        }
      }

    }

    LoadFiles( filestosource )

    return( filestosource )

	}


text:
model_phyllochron <- function (fixPhyll = 5.0,
         leafNumber = 0.0,
         lincr = 8.0,
         ldecr = 0.0,
         pdecr = 0.4,
         pincr = 1.5,
         ptq = 0.0,
         gAImean = 0.0,
         kl = 0.45,
         pTQhf = 0.0,
         B = 20.0,
         p = 120.0,
         choosePhyllUse = 'Default',
         areaSL = 0.0,
         areaSS = 0.0,
         lARmin = 0.0,
         lARmax = 0.0,
         sowingDensity = 0.0,
         lNeff = 0.0){
    #'- Name: Phyllochron -Version: 1.0, -Time step: 1
    #'- Description:
    #'            * Title: Phyllochron Model
    #'            * Author: Pierre Martre
    #'            * Reference: Modeling development phase in the
    #'                Wheat Simulation Model SiriusQuality.
    #'                See documentation at http://www1.clermont.inra.fr/siriusquality/?page_id=427
    #'            * Institution: INRA Montpellier
    #'            * Abstract: Calculate different types of phyllochron
    #'- inputs:
    #'            * name: fixPhyll
    #'                          ** description : Sowing date corrected Phyllochron
    #'                          ** inputtype : variable
    #'                          ** variablecategory : auxiliary
    #'                          ** datatype : DOUBLE
    #'                          ** default : 5.0
    #'                          ** min : 0.0
    #'                          ** max : 10000.0
    #'                          ** unit : °C d leaf-1
    #'                          ** uri : some url
    #'            * name: leafNumber
    #'                          ** description : Actual number of phytomers
    #'                          ** inputtype : variable
    #'                          ** variablecategory : state
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 25.0
    #'                          ** unit : leaf
    #'                          ** uri : some url
    #'            * name: lincr
    #'                          ** description : Leaf number above which the phyllochron is increased by Pincr
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 8.0
    #'                          ** min : 0.0
    #'                          ** max : 30.0
    #'                          ** unit : leaf
    #'                          ** uri : some url
    #'            * name: ldecr
    #'                          ** description : Leaf number up to which the phyllochron is decreased by Pdecr
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 100.0
    #'                          ** unit : leaf
    #'                          ** uri : some url
    #'            * name: pdecr
    #'                          ** description : Factor decreasing the phyllochron for leaf number less than Ldecr
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.4
    #'                          ** min : 0.0
    #'                          ** max : 10.0
    #'                          ** unit : -
    #'                          ** uri : some url
    #'            * name: pincr
    #'                          ** description : Factor increasing the phyllochron for leaf number higher than Lincr
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 1.5
    #'                          ** min : 0.0
    #'                          ** max : 10.0
    #'                          ** unit : -
    #'                          ** uri : some url
    #'            * name: ptq
    #'                          ** description : Photothermal quotient
    #'                          ** inputtype : variable
    #'                          ** variablecategory : state
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 10000.0
    #'                          ** unit : MJ °C-1 d-1 m-2)
    #'                          ** uri : some url
    #'            * name: gAImean
    #'                          ** description : Green Area Index
    #'                          ** inputtype : variable
    #'                          ** variablecategory : state
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 10000.0
    #'                          ** unit : m2 m-2
    #'                          ** uri : some url
    #'            * name: kl
    #'                          ** description : Exctinction Coefficient
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.45
    #'                          ** min : 0.0
    #'                          ** max : 50.0
    #'                          ** unit : -
    #'                          ** uri : some url
    #'            * name: pTQhf
    #'                          ** description : Slope to intercept ratio for Phyllochron  parametrization with PhotoThermal Quotient
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : °C d leaf-1
    #'                          ** uri : some url
    #'            * name: B
    #'                          ** description : Phyllochron at PTQ equal 1
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 20.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : °C d leaf-1
    #'                          ** uri : some url
    #'            * name: p
    #'                          ** description : Phyllochron (Varietal parameter)
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : DOUBLE
    #'                          ** default : 120.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : °C d leaf-1
    #'                          ** uri : some url
    #'            * name: choosePhyllUse
    #'                          ** description : Switch to choose the type of phyllochron calculation to be used
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : species
    #'                          ** datatype : STRING
    #'                          ** default : Default
    #'                          ** min :
    #'                          ** max :
    #'                          ** unit : -
    #'                          ** uri : some url
    #'            * name: areaSL
    #'                          ** description :  Area Leaf
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : cm2
    #'                          ** uri : some url
    #'            * name: areaSS
    #'                          ** description : Area Sheath
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : cm2
    #'                          ** uri : some url
    #'            * name: lARmin
    #'                          ** description : LAR minimum
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : leaf-1 °C
    #'                          ** uri : some url
    #'            * name: lARmax
    #'                          ** description : LAR maximum
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : leaf-1 °C
    #'                          ** uri : some url
    #'            * name: sowingDensity
    #'                          ** description : Sowing Density
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : plant m-2
    #'                          ** uri : some url
    #'            * name: lNeff
    #'                          ** description : Leaf Number efficace
    #'                          ** inputtype : parameter
    #'                          ** parametercategory : genotypic
    #'                          ** datatype : DOUBLE
    #'                          ** default : 0.0
    #'                          ** min : 0.0
    #'                          ** max : 1000.0
    #'                          ** unit : leaf
    #'                          ** uri : some url
    #'- outputs:
    #'            * name: phyllochron
    #'                          ** description :  the rate of leaf appearance
    #'                          ** variablecategory : state
    #'                          ** datatype : DOUBLE
    #'                          ** min : 0
    #'                          ** max : 1000
    #'                          ** unit :  °C d leaf-1
    #'                          ** uri : some url
    phyllochron <- 0.0
    LAR <- 0.0
    gaiLim <- lNeff * ((areaSL + areaSS) / 10000.0) * sowingDensity
    if (choosePhyllUse == 'Default')
    {
        if (leafNumber < ldecr)
        {
            phyllochron <- fixPhyll * pdecr
        }
        else if ( leafNumber >= ldecr && leafNumber < lincr)
        {
            phyllochron <- fixPhyll
        }
        else
        {
            phyllochron <- fixPhyll * pincr
        }
    }
    if (choosePhyllUse == 'PTQ')
    {
        if (gAImean > gaiLim)
        {
            LAR <- (lARmin + ((lARmax - lARmin) * ptq / (pTQhf + ptq))) / (B * gAImean)
        }
        else
        {
            LAR <- (lARmin + ((lARmax - lARmin) * ptq / (pTQhf + ptq))) / (B * gaiLim)
        }
        phyllochron <- 1.0 / LAR
    }
    if (choosePhyllUse == 'Test')
    {
        if (leafNumber < ldecr)
        {
            phyllochron <- p * pdecr
        }
        else if ( leafNumber >= ldecr && leafNumber < lincr)
        {
            phyllochron <- p
        }
        else
        {
            phyllochron <- p * pincr
        }
    }
    return (list('phyllochron' = phyllochron))
}


text:
JOBID = "276183"
source("Italy/code/utils/make-table.r")

# Prints attack rates to console
scenario_type = "constant-com"
mobility_increase = 0
compliance_decrease = 0
make_table_simulation(paste0('Italy/results/sim-', scenario_type, "-", StanModel, '-', len_forecast, '-', compliance_decrease, '-', JOBID, '-stanfit.Rdata'),
                      date_till_percentage = max(dates[[1]]) + len_forecast)

for (com_dec in c(20, 40)){
  # Prints attack rates to console
  scenario_type = "decrease-com-current"
  mobility_increase = 0
  compliance_decrease = com_dec
  make_table_simulation(paste0('Italy/results/sim-', scenario_type, "-", StanModel, '-', len_forecast, '-', compliance_decrease, '-', JOBID, '-stanfit.Rdata'),
                        date_till_percentage = max(dates[[1]]) + len_forecast)
}


text:
check.dim <- function(data) {
	if (is.na(dim(data)[1])) {
		data <- matrix(data, length(data), 1)
	}
	return(data)
}

#impute using the mean (assuming it's all numeric)
mean.imp <- function(data) {
	data <- check.dim(data)
	#stop('Complete this function.')
	#lapply(data, )
	for (i in 1:ncol(data)){
	  data[is.na(data[,i]),i] = mean(data[,i], na.rm=TRUE)
	}
	#for each variable in data, impute the missing values using the mean of the available values
	return(data)
}

#add missingness to a matrix or dataframe
miss.set.df <- function(dat, prop, type.miss='mcar') {
	dat <- check.dim(dat)
	for (i in 1:dim(dat)[2]) {
		dat[,i] <- miss.set.vec(dat[,i], prop, type.miss)
	}
	return(dat)
}
# what is this ?miss.set.vec

#add misingness to a vector; mcar = at random; mnar = not at random (remove higher values)
miss.set.vec <- function(x, prop, type.miss='mcar') {
	n.miss <- rbinom(1, length(x), prop)
	if (type.miss == 'mcar') {
		miss.idx <- sample(1:length(x), n.miss, replace=F)
	} else if (type.miss == 'mnar') {
		miss.idx <- order(x, decreasing=T)[1:n.miss]
	}
	x[miss.idx] <- NA
	return(x)
}

split.data <- function(data, train.prop, set.seed=NA) {
	if (!is.na(set.seed)) {set.seed(set.seed)}
	train.idx <- sample(1:dim(data)[1], round(dim(data)[1]*train.prop), replace=F)
	test.idx <- setdiff(1:dim(data)[1], train.idx)
	train.set <- data[train.idx,]
	test.set <- data[test.idx,]
	return(list(train=train.set, test=test.set))
}

get.resid.err <- function(train.set, test.set) {
	#Train the models and get the MSE
	#stop('Complete this function.')

	#fit a linear model on the training set (train.set) using only the intercept Y ~ 1
  model1 = lm(Y ~1, data =  train.set)

	#fit a linear model on the training set using all covariates Y ~ .
  modelfull = lm(Y ~., train.set)

	#predict from each model on the test set (test.set)
  predict1 = predict(model1, newdata = test.set[,-91])
  predict.full = predict(modelfull, newdata = test.set[,-91])

	#calculate the mse of the null model
  mse.null = mean((predict1 - test.set$Y)^2)

	#calculate the mse of the full model
  mse.full = mean((predict.full - test.set$Y)^2)

	#calculate the residual error of the ratio of the latter over the former
	pred.res <- mse.full/ mse.null

	return(pred.res)
}

missing_function = function(miss_type_idx, miss_prop_idx, data_idx){

  train.prop <- 0.50
  train.split.seed <- 78901
  outcome.var <- "Y"
  rf.iter <- 12
  mice.iter <- 10

  miss.prop.list <- c(0.10, 0.20, 0.30, 0.40) #proportion of missingness to assign
  miss.type.list <- c('mcar', 'mnar') #mcar = missing completely at random; mnar = missing not at random
  data.idx <- data_idx #the dataset to use
  miss.type.idx <- miss_type_idx #type of missingness to use
  miss.prop.idx <- miss_prop_idx #proportio of missingness to use

  t.st <- Sys.time() #time the run

  data <- data.list[[data.idx]]
  miss.prop <- miss.prop.list[miss.prop.idx]
  miss.type <- miss.type.list[miss.type.idx]

  split <- split.data(data, train.prop, set.seed=train.split.seed)
  train.set <- split$train; test.set <- split$test
  pred.var <- colnames(data)[-which(colnames(data) == outcome.var)]

  ##########################
  ## Set the missing data ##
  train.set.miss <- train.set
  train.set.miss[,pred.var] <- miss.set.df(train.set.miss[,pred.var], miss.prop, miss.type)

  dim(train.set.miss)
  length(pred.var)

  ################
  #MICE Imputation
  mice.data <- train.set.miss
  dim(mice.data)

  mice.data0 = mice.data[,pred.var]
  tempdata = mice(mice.data0,m=mice.iter,maxit=25,meth='pmm',seed=train.split.seed)

  # Will take the exptectation of the imputated data sets to obtain one data set.
  imputed = complete(tempdata, 1)
  for (i in 2:mice.iter){
    print(i)
    print(dim(imputed))
    imputed = imputed + complete(tempdata, i)
    print(dim(imputed))
  }

  imputed = imputed / mice.iter

  #run an imputation on the mice.data dataframe using the mice function from the mice package
  #use 'mice.iter' datasets and cap the iterations at 25
  #DO NOT INCLUDE the outcome 'Y' in the imputation, impute only mice.data[,pred.var]
  #try to parallelize the computation by running several interations in parallel

  mice.data.comp = data.frame(imputed, Y = mice.data$Y)
  #write.csv(mice.data.comp, "mice_data_comp.csv", row.names = FALSE)
  #mice.data.comp = read.csv("mice_data_comp.csv")
  # same, but now as list, mild object
  # dslist <- complete(tempdata, "all")
  # length(dslist)
  # imputed_list_data = mean(dslist)

  # use package miceadds to save the imputed datasets.
  # require(miceadds)
  # write.mice.imputation(mi.res=tempdata, "tempdata1", include.varnames=TRUE,
  #       long=TRUE, mids2spss=TRUE, spss.dec=",", dattype=NULL)

  # # Parallalize mice imputation.
  # total.cores <- detectCores()
  # tempdata_core = parlmice(mice.data0, m = mice.iter, seed = NA, cluster.seed = 500,
  #                          n.core = total.cores,  n.imp.core = NULL, cl.type = "PSOCK")
  mice.data = mice.data.comp
  t.mice = Sys.time() - t.st

  t.st.mean = Sys.time()
  #####################
  #Impute with the mean
  mean.data <- train.set.miss

  #finish the mean.imp function from the functions file and impute mean.data[,pred.var]
  # maybe should only use this mean.data[,pred.var]

  head(is.na(mean.data))
  is.na(dim(mean.data)[1])
  check.dim(mean.data)

  # RUN the mean.imp function.
  mean.data = mean.imp(mean.data)
  all(!is.na(mean.data))

  # Testing if the mean.imp() function works.
  all(sapply(train.set.miss, function(x) mean(x, na.rm  = T)) == sapply(mean.data, mean))
  summary(sapply(train.set.miss, function(x) mean(x, na.rm  = T)))

  t.mean = Sys.time() - t.st.mean

  t.st.rf = Sys.time()
  #########################
  #Random Forest Imputation
  total.cores <- detectCores()
  print(total.cores)
  cl <- makeCluster(total.cores)
  registerDoParallel(cl)

  rf.data <- train.set.miss[,pred.var]
  dim(rf.data)

  #impute using random forest imputation, use 500 trees, and cap the number of iterations at 12 (rf.iter)
  #try to parallelize the forests using the doParallel package to save time
  #?missForest

  class(rf.data)
  set.seed(train.split.seed)
  rf.data.comp = missForest(xmis = rf.data, maxiter = rf.iter, ntree = 500, parallelize = c('forests'))
  rf.data.comp.train = data.frame(rf.data.comp$ximp, Y = train.set.miss$Y)
  #write.csv(rf.data.comp.train, "rf_data_comp_train_final.csv", row.names = F)
  #rf.data.comp.train = read.csv("rf_data_comp_train_final.csv")
  #file.remove("rf.data.comp.train.csv")
  # if (file.exists(fn))
  #   #Delete file if it exists
  #   file.remove(fn)
  rf.data = rf.data.comp.train

  t.rf = Sys.time() - t.st.rf

  ##############################
  #finish the get.resid.err function and calculate the test set errors for each imputed dataset
  mean.imp.err <- get.resid.err(mean.data, test.set)
  mice.imp.err <- get.resid.err(mice.data, test.set)
  rf.imp.err <- get.resid.err(rf.data, test.set)
  no.imp.err <- get.resid.err(train.set, test.set)

  # mean.imp.err
  # mice.imp.err
  # rf.imp.err
  # no.imp.err

  t.end = Sys.time() - t.st
  print(t.end)

  return(data.frame(mean.imp.err, mice.imp.err, rf.imp.err, no.imp.err, t.mice, t.mean, t.rf, t.end))

}


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

#Produce line graph summarising the di-tag size distribution
#Launched by hicup_filter
args <- commandArgs(TRUE)
outdir <- args[1]
file <- args[2]

data <- read.delim(file, header=FALSE, skip=1)

if(length(data) > 0){

	outputfilename <- basename(file)
	outputfilename <- paste(outputfilename, "svg", sep = ".")
	outputfilename <- paste(outdir, outputfilename, sep = "")
	svg(file=outputfilename)

	data <- read.delim(file, header=FALSE, skip=0)

	plot (data, type="l", xlab="Di-tag length (bps)",
	     ylab="Frequency", main="Di-tag frequency Vs Length"
	    )

	garbage <- dev.off()

}else{
  no_data_message = paste("R: No data in ", file)
  print(no_data_message)
}

max_stars_repo_name: PavriLab/hicer-nf

max_issues_repo_name: PavriLab/hicer-nf

max_forks_repo_name: PavriLab/hicer-nf


text:
## De Groot (4th ed.) exemplo 7.6.11

gera_dados <- function(n, theta){
  X <- rexp(n = n, rate = theta)
  return(X)
}

computa_emv <- function(x){
  theta.chapeu <- 1/mean(x)
}

estima <- function(n, theta){
  dados <- gera_dados(n = n, theta = theta)
  est <- computa_emv(dados)
  return(est)
}

dinvgamma <- function(x, alpha, beta, log = FALSE){
  ## densidade de uma gama inversa com parametros alpha e beta
  ## esta é a distribuição **exata** do EMV neste caso
  lconst <- alpha * log(beta) - lgamma(alpha)
  ldens <- -(alpha + 1)*log(x) - beta/x
  ans <- lconst + ldens
  if(!log) ans <- exp(ans)
  return(ans)
}
dinvgamma <- Vectorize(dinvgamma)
#############

M <- 1E4 ## repetições

theta.vdd <- .234

par(mfrow=c(2, 2))

## n = 10
n10 <- 10

ests.n10 <- sapply(1:M, function(i) estima(n = n10, theta = theta.vdd) )

hist(ests.n10, probability = TRUE, xlab = expression(hat(theta)), main = "Tamanho de amostra = 10")
curve(dnorm(x, mean = theta.vdd, sd = sqrt(theta.vdd^2/n10)), lwd = 2, add = TRUE)
curve(dinvgamma(x, alpha = n10, beta = n10*theta.vdd), lwd = 2, lty = 2, col = 2, add = TRUE)
abline(v = theta.vdd, lwd = 2, lty = 3)

## n = 50
n50 <- 50

ests.n50 <- sapply(1:M, function(i) estima(n = n50, theta = theta.vdd) )

hist(ests.n50, probability = TRUE, xlab = expression(hat(theta)), main = "Tamanho de amostra = 50")
curve(dnorm(x, mean = theta.vdd, sd = sqrt(theta.vdd^2/n50)), lwd = 2, add = TRUE)
curve(dinvgamma(x, alpha = n50, beta = n50*theta.vdd), lwd = 2, lty = 2, col = 2, add = TRUE)
abline(v = theta.vdd, lwd = 2, lty = 3)

## n = 100
n100 <- 100

ests.n100 <- sapply(1:M, function(i) estima(n = n100, theta = theta.vdd) )

hist(ests.n100, probability = TRUE, xlab = expression(hat(theta)), main = "Tamanho de amostra = 100")
curve(dnorm(x, mean = theta.vdd, sd = sqrt(theta.vdd^2/n100)), lwd = 2, add = TRUE)
curve(dinvgamma(x, alpha = n100, beta = n100*theta.vdd), lwd = 2, lty = 2, col = 2, add = TRUE)
abline(v = theta.vdd, lwd = 2, lty = 3)

## n = 100
n1000 <- 1000

ests.n1000 <- sapply(1:M, function(i) estima(n = n1000, theta = theta.vdd) )

hist(ests.n1000, probability = TRUE, xlab = expression(hat(theta)), main = "Tamanho de amostra = 1000")
curve(dnorm(x, mean = theta.vdd, sd = sqrt(theta.vdd^2/n1000)), lwd = 2, add = TRUE)
curve(dinvgamma(x, alpha = n1000, beta = n1000*theta.vdd), lwd = 2, lty = 2, col = 2, add = TRUE)
abline(v = theta.vdd, lwd = 2, lty = 3)


text:
#' Generate randomly placed intervals on a genome.
#'
#' @param genome [genome_df]
#' @param length length of intervals
#' @param n number of intervals to generate
#' @param seed seed RNG for reproducible intervals
#' @param sorted return sorted output
#'
#' @details Sorting can be suppressed with `sorted = FALSE`.
#'
#' @return [ivl_df]
#'
#' @family randomizing operations
#'
#' @seealso \url{https://bedtools.readthedocs.io/en/latest/content/tools/random.html}
#'
#' @examples
#' genome <- tibble::tribble(
#'   ~chrom,  ~size,
#'   "chr1",  10000000,
#'   "chr2",  50000000,
#'   "chr3",  60000000,
#'   "chrX",  5000000
#' )
#'
#' bed_random(genome, seed = 10104)
#'
#' # sorting can be suppressed
#' bed_random(genome, sorted = FALSE, seed = 10104)
#'
#' # 500 random intervals of length 500
#' bed_random(genome, length = 500, n = 500, seed = 10104)
#'
#' @export
bed_random <- function(genome, length = 1000, n = 1e6, seed = 0, sorted = TRUE) {
  genome <- check_genome(genome)

  if (!all(genome$size > length)) {
    stop("`length` must be greater than all chrom sizes", call. = FALSE)
  }

  out <- random_impl(genome, length, n, seed)

  if (sorted) {
    out <- bed_sort(out)
  }

  as_tibble(out)
}

max_stars_repo_name: rnabioco/valr

max_issues_repo_name: rnabioco/valr

max_forks_repo_name: rnabioco/valr


text:
 E04UCJ Example Program Results

 Varbl  Istate   Value         Lagr Mult

 V   1   1       1.00000         1.088
 V   2   0       4.74300         0.000
 V   3   0       3.82115         0.000
 V   4   0       1.37941         0.000

 L Con  Istate   Value         Lagr Mult

 L   1   0       10.9436         0.000

 N Con  Istate   Value         Lagr Mult

 N   1   2       40.0000       -0.1615
 N   2   1       25.0000        0.5523

 Final objective value =    17.01402


text:
#' metacell
#'
#' @import dplyr
#' @import ggplot2
#' @import tgconfig
#' @import tgstat
#' @import tgutil
#' @import Matrix
#' @importClassesFrom Matrix Matrix dgCMatrix dgeMatrix
#' @importFrom tidyr pivot_longer pivot_wider spread gather
#' @importFrom dbscan dbscan
#' @importFrom parallel mclapply
#' @importFrom zoo rollmean rollmedian
#' @importFrom pdist pdist
#' @importFrom cluster silhouette
#' @importFrom igraph cluster_louvain graph_from_data_frame
#' @importClassesFrom graph graphNEL
#' @importFrom graph plot addEdge addNode nodeRenderInfo
#' @importFrom Rgraphviz layoutGraph
#' @importFrom umap umap
#' @name metacell
#' @docType package
NULL

# this is important for running this package from scripts
#' @import methods
NULL

max_stars_repo_name: ofirr/metacell

max_issues_repo_name: lbgbox/metacell

max_forks_repo_name: lbgbox/metacell


text:
install.packages("sparklyr")
install.packages("tidyverse")

library(sparklyr)
library(tidyverse)

spark_install(version = "3.0.0")
sc  <- spark_connect(master = "local")

df1 <- as_tibble(iris)
df  <- copy_to(sc, df1)

print(class(df1))
print(class(df))
# df1 is a data.frame, df is a function
# different colnames, one is in spark formatting vs. df

df %>% select(Sepal_Length,Species) %>% head %>% print
df %>% filter(Sepal_Length>5.5) %>% head %>% print
df %>% select(Sepal_Length,Species) %>% filter(Sepal_Length>5.5) %>% head %>% print
df2 <- df %>% group_by(Species) %>% summarize(mean = mean(Sepal_Length), count = n()) %>% head %>% print
df2 <- df2 %>% group_by(Species) %>% summarize(mean = mean(Sepal_Length), count = n()) %>% head %>% print
df2 %>% arrange(Species) %>% head %>% print

max_stars_repo_path: atom code/ps4b.r

max_issues_repo_path: atom code/ps4b.r

max_forks_repo_path: atom code/ps4b.r


text:
# lm within levels of a factor
# result has "final.model" which contains a list "models" of
# lm models indexed by category levels and a string "category" name
"lm.by.category" <- function(x,y,category){
  results <- list()
  final.model <- list()
  final.model$models <- list()
  category <- droplevels(as.factor(as.character(category)))
  yhat <- y
  for(level in levels(category)){
    ix <- category == level
    if(sum(ix) >= 3){
      final.model$models[[level]] <- lm(y[ix] ~ x[ix])
      yhat[ix] <- final.model$models[[level]]$fitted.values
    }
  }
  # make one general model for predicting unrecognized categories
  final.model$general.model <- lm(y ~ x)

  results$final.model <- final.model
  results$rmse <- sqrt(mean((yhat-y)^2))
  results$mape <- 100*mean(abs(yhat-y)/y)
  results$y <- y
  results$yhat <- yhat
  return(results)
}

# lm within levels of a factor
# "final.model" which contains a list "models" of
# lm models indexed by category levels
"predict.lm.by.category" <- function(final.model, newx, category){
  yhat <- numeric(length(newx))
  category <- droplevels(as.factor(as.character(category)))
  for(level in levels(category)){
    ix <- category == level
    if(level %in% names(final.model$models)){
      # this level is in the list of models
      model.i <- final.model$models[[level]] # get model for this category
    }
    else{
      # this level is not in the list of models
      model.i <- final.model$general.model # get general model
    }
    intercept <- model.i$coefficients[1]
    slope <- model.i$coefficients[2]
    yhat[ix] <- slope * newx[ix] + intercept
  }
  return(yhat)
}

"plot.lm.by.category" <- function(x,y,category,x.name='Deadweight Tons',y.name='KG CO2/NM'){
  par(mfrow=c(3,5),mar=c(4,4,5,.4))
  yhat <- y
  category <- droplevels(as.factor(as.character(category)))
  for(level in levels(category)){
    ix <- category == level
    if(sum(ix) >= 5){
      linearmodel <- lm(y[ix] ~ x[ix])
      yhat[ix] <- linearmodel$fitted.values
      rmse <- sqrt(mean((yhat[ix]-y[ix])^2))
      alpha.char <- '22'
      if(sum(ix) < 100 ) alpha.char <- '55'
      plot(x[ix],y[ix],
           pch=21,col=NA,
           bg=sprintf('#000000%s',alpha.char),
           xlab=x.name,
           ylab=y.name,
           main=sprintf('%s \n(RMSE %0.3f)',level, rmse),
           cex.main=1.1,
           cex.lab=0.75)
      abline(linearmodel)
    }
  }
}

max_stars_repo_name: knights-lab/climate-trace-shipping

max_issues_repo_name: knights-lab/climate-trace-shipping

max_forks_repo_name: knights-lab/climate-trace-shipping


text:
digits <- c("-**----*--***--***---*---****--**--****--**---**--",
            "*--*--**-----*----*-*--*-*----*-------*-*--*-*--*-",
            "*--*---*---**---**--****-***--***----*---**---***-",
            "*--*---*--*-------*----*----*-*--*--*---*--*----*-",
            "-**---***-****-***-----*-***---**---*----**---**--",
            "--------------------------------------------------")
numf <- function(s) {
  strsplit(gsub("[^0123456789]", "", s), " ")[[1]]
}

for (line in readLines(tail(commandArgs(), n=1))) {
  cat(sapply(lapply(line, numf), function(s) {
    r <- c("", "", "", "", "", "")
    if (length(s) > 0L) {
      for (i in strsplit(s, "")[[1]]) {
        for (j in 1:6) {
          p <- as.integer(i)*5 + 1
          r[j] <- paste(r[j], substr(digits[j], p, p+4), sep="")
        }
      }
    }
    paste(r, collapse="\n")
  }), sep="\n")
}

max_stars_repo_name: IlkhamGaysin/ce-challenges

max_issues_repo_name: IlkhamGaysin/ce-challenges

max_forks_repo_name: IlkhamGaysin/ce-challenges


text:
pdf_file<-"pdf/maps_germany_cities_piecharts_outline.pdf"
cairo_pdf(bg="grey98", pdf_file,width=13,height=11)
myColour_h<-rgb(240,240,240,maxColorValue=255)

par(omi=c(0.5,0.5,0.5,0),mai=c(0,0,0,0),lend=1,family="Lato Light",las=1)
library(sp)
library(plotrix)
library(gdata)

# Import data and create chart

myData<-"myData/data_maps_germany_cities_piecharts.xlsx"
myCountries<-read.xls(myData,head=T, encoding="latin1")
attach(myCountries)

load("myData/gadm2/DEU_adm1.RData")
plot(gadm,border=rgb(151,151,151,maxColorValue=255),lwd=0.5)

load("myData/gadm2/DEU_adm0.RData")
plot(gadm,border="black",lwd=0.95,add=T)

# Other elements

n<-nrow(myCountries)
for (i in 1:n)
{
myCircle<-c(agriculture[i],forest[i],cultivation[i])
floating.pie(long[i],lat[i],myCircle,radius=0.5,col=c(rgb(215,215,0,150,maxColorValue=255),rgb(34,139,34,150,maxColorValue=255),rgb(178,34,34,150,maxColorValue=255)),border=F)
text(long[i],lat[i],name[i])
}
floating.pie(16,48.25,c(1,1,1),radius=0.5,col=c(rgb(215,215,0,150,maxColorValue=255),rgb(34,139,34,150,maxColorValue=255),rgb(178,34,34,150,maxColorValue=255)),border=F)
text(16.75,47.75,"Settlement area")
text(14.7,48.25,"Forest area")
text(16.75,48.75,"Agriculture")

# Titling

mtext("Shares of agricultural, forest and settlement area",3,line=-0.25,adj=0,cex=2.25,outer=T,family="Lato Black")
mtext("Source: Statistisches Bundesamt Fachserie 3 Reihe 5.1",1,line=0,adj=0.9,cex=1.25,font=3)
dev.off()

max_stars_repo_name: wilsonify/data-visualization

max_issues_repo_name: wilsonify/data-visualization

max_forks_repo_name: wilsonify/data-visualization


text:
# Packages used  : pracma
# To install pracma,type following in command line while connected to internet
# install.packages("pracma")
# package can be included by command " library(pracma) "
# for more information about pracma visit https://cran.r-project.org/web/packages/pracma/index.html

# Example : 4   Chapter : 2.5   Pageno : 85
# Inverse of matrix by guass jordan elimination matrix
library(pracma)
A<-matrix(c(2,4,3,7),ncol=2)
I<-eye(2)
AI<-cbind(A,I)
R<-rref(AI)
X<-R[,c(3:4)]
print("Inverse of A is")
print(X)


text:
#' Function that models a d18O curve through SST and GR sinusoids
#'
#' The core function of the ShellChron growth model. Uses growth
#' rate and SST (Sea Surface Temperature) sinusoids to model d18O
#' data to be matched with the input. In the ShellChron modeling
#' routine, this function is optimized using the SCEUA algorithm
#' and applied on sliding windows through the dataset to estimate
#' the age of each datapoint
#' @param pars List of parameters for temperature and growth rate sinusoids
#' \code{pars <- c(T_amp, T_pha, T_av, G_amp, G_pha, G_av, G_skw)}
#' @param T_per Period of SST sinusoid (in days; default = 365)
#' @param G_per Period of growth rate sinusoid (in days; default = 365)
#' @param years Number of years to be modeled (default = 1)
#' @param t_int Time interval (in days; default = 1)
#' @param transfer_function Transfer function used to convert d18Oc to temperature
#' data.
#' @param d18Ow Either a single value (constant d18Ow) or a vector of length
#' equal to the period in SST data (365 days by default) containing information
#' about seasonality in d18Ow. Defaults to constant d18Ow of 0 permille VSMOW
#' (the modern mean ocean value)
#' @param Dsam Vector of \code{D} values serving as input (keep unit consistent
#' throughout model)
#' @param Osam Vector of \code{d18Oc} values serving as input (in permille VPDB)
#' @param t_maxtemp Timing of the warmest day of the year (in julian day;
#' default = 182.5, or May 26th halfway through the year)
#' @param plot Should results of modeling be plotted? \code{TRUE/FALSE}
#' @param MC Number of Monte Carlo simulations to apply for error propagation
#' Default = 1000
#' @param D_err OPTIONAL: Vector containing errors on \code{Dsam}
#' @param O_err OPTIONAL: Vector containing errors on \code{Osam}
#' @param return String indicating whether to return just the Sum of Squared
#' Residuals ("SSR") or a matrix containing the results of the model and the
#' propagated uncertainties (if applicable)
#' @return Depending on the value of the "return" parameter either a single
#' value representing the Sum of Squared Residuals ("SSR") as a measure for
#' the closeness of the match between modeled d18O and input values, or a
#' matrix containing the full result of the modeling including propagated
#' uncertainties if applicable.
#' @references package dependencies: ggplot2 3.2.1
#' function dependencies: temperature_curve, d18O_model, growth_rate_curve,
#' mc_err_orth
#'
#'   \doi{10.1016/j.palaeo.2017.09.034}
#' @examples
#' # Set parameters
#' G_amp <- 20
#' G_per <- 365
#' G_pha <- 100
#' G_av <- 15
#' G_skw <- 70
#' T_amp <- 20
#' T_per <- 365
#' T_pha <- 150
#' T_av <- 15
#' pars <- c(T_amp, T_pha, T_av, G_amp, G_pha, G_av, G_skw)
#' d18Ow <- 0
#' # Create dummy data
#' Dsam <- seq(1, 40, 1)
#' Osam <- sin((2 * pi * (seq(1, 40, 1) - 8 + 30 / 4)) / 30)
#' # Test returning residual sum of squares for optimization
#' SSR <- growth_model(pars, T_per, G_per, Dsam = Dsam, Osam = Osam,
#'     return = "SSR")
#' # Test returning full model result
#' resmat <- growth_model(pars, T_per, G_per, Dsam = Dsam, Osam = Osam,
#'     return = "result")
#' @export
growth_model <- function(pars, # Growth model function to optimize using sceua
    T_per = 365, # Period of SST sinusoid (default = 365 days)
    G_per = 365, # Period of growth sinusoid (default = 365 days)
    years = 1, # Default year = 1
    t_int = 1, # Time interval, default = 1 day
    transfer_function = "KimONeil97", # Transfer function used to convert d18Oc values of the material into temperature (default = "KimONeil97"; Kim and O'Neil, 1997)
    d18Ow = "default", # d18Ow vector, default = constant at 0 permille VSMOW
    Dsam, # Depth data
    Osam, # d18O data
    t_maxtemp = 182.5, # Time (day) at which maximum temperature is reached. Default is 182.5 (exactly halfway through the year, or 1st of June)
    plot = FALSE, # Plot results?
    MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (error on D and d18O measurements not considered)
    D_err = NULL, # Optional: include standard deviation of uncertainty on depth measurements, default = empty (no error incorporated in the calculation, same effect as MC = 0)
    O_err = NULL, # Optional: include standard deviation of uncertainty on d18O measurements, default = empty (no error incorporated in the calculation, same effect as MC = 0)
    return = "SSR" # What to return? Default = the total sum of squares to be minimized by the SCEUA algorithm
    ){

    # Extract variables
    T_amp <- pars[1]
    T_pha = pars[2]
    T_av = pars[3]

    G_amp = pars[4]
    G_pha = pars[5]
    G_av = pars[6]
    G_skw = pars[7]

    T_par <- c(T_amp, T_per, T_pha, T_av)
    G_par <- c(G_amp, G_per, G_pha, G_av, G_skw)

    # Run functions for SST and GR
    SST <- temperature_curve(T_par, years, t_int)
    d18Oc <- d18O_model(SST, d18Ow, transfer_function)
    if(length(d18Oc) == 1){ # Catch errors in d18O model
        return("ERROR: Supplied transfer function is not recognized")
    }
    GR <- growth_rate_curve(G_par, years, t_int)

    # Create distance (D) vector from GR
    D <- cumsum(GR[,2]) + Dsam[1]

    # Compare model result with data
    Omod <- approx( # Interpolate modeled d18Oc values to positions along the record.
        x = D,
        y = d18Oc[,2],
        xout = Dsam,
        method = "linear",
        rule = 2
    )

    residuals <- Osam - Omod$y # Calculate residuals

    if(return == "SSR"){
        SSR <- sum(residuals^2) # Calculate sum of squared residuals
        return(SSR)
    }else{
        # Compare model result with data
        t <- approx( # Interpolate modeled time values to positions along the record.
            x = D,
            y = d18Oc[,1],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        gr <- approx( # Interpolate modeled growth rate values to positions along the record.
            x = D,
            y = GR[,2],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        Tmod <- approx( # Interpolate modeled temperature values to positions along the record.
            x = D,
            y = SST[,2],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        TY <- (t$y - T_pha + t_maxtemp) %% 365 # Calculate time of year relative to the first T maximum in the record (T_pha) and the assumed day of the year when it occurs (default = 182.5)

        # Prepare export matrix
        resmat <- cbind(Dsam, Osam, Omod$y, residuals, TY, gr$y, Tmod$y) # Calculate matrix of end results
        colnames(resmat) <- c("Dsam", "Osam", "Omod", "Residuals", "Time_of_year", "Instantaneous_growth_rate", "Modeled_temperature")

        # Optional: Incorporate measurement errors on D and d18O
        if(MC > 0){
            print("Propagating measurement uncertainties")
            if(is.null(D_err)){ # If no error on D is given, but error analysis is asked (MC > 0), then D_err is zero for all samples
                D_err <- rep(0, length(Dsam))
            }
            if(is.null(O_err)){ # If no error on d18Oc is given, but error analysis is asked (MC > 0), then O_err is zero for all samples
                O_err = rep(0, length(Osam))
            }
            # Propagate the combined effect of error on D and d18Oc on the model fit
            Proj_err <- mc_err_orth(Dsam, D_err, Osam, O_err, D, d18Oc, MC) # Combine error on D and d18Oc through orthogonal projection on the modeled D-d18Oc curve
            D_err_comb <- Proj_err$X_err # Isolate uncertainty in X-domain for further processing

            Drange <- cbind((Dsam - D_err_comb) %% D[length(D)], (Dsam + D_err_comb) %% D[length(D)]) # Find upper and lower boundaries of D error (1 SD)
            Drange_pos <- apply(abs(outer(Drange, D, FUN = "-")), c(1, 2), which.min) # Find positions of D ranges

            # Approximate SDs of d18Oc, t, gr and Tmod from range of D values within +/- 1 SD
            Omod_SD <- vector(length = length(Drange_pos[,1]))
            t_SD <- vector(length = length(Drange_pos[,1]))
            gr_SD <- vector(length = length(Drange_pos[,1]))
            Tmod_SD <- vector(length = length(Drange_pos[,1]))
            for(j in 1:length(Drange_pos[,1])){ # Loop through all positions in the window and propagate the D error onto the results
                if(Drange_pos[j, 2] < Drange_pos[j, 1]){ # If the D range contains the year edge, incorporate the outer edges of the year.
                # Find the range of modeled values contained in the 1 SD range of D error and approximate SD of modeled values from in and max values
                    Omod_SD[j] <- (max(d18Oc[c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2]), 2]) - min(d18Oc[c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                    t_SD[j] <- length(c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2])) / 2
                    gr_SD[j] <- (max(GR[c(Drange_pos[j, 1]:length(GR[, 1]), 1:Drange_pos[j, 2]), 2]) - min(GR[c(Drange_pos[j, 1]:length(GR[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                    Tmod_SD[j] <- (max(SST[c(Drange_pos[j, 1]:length(SST[, 1]), 1:Drange_pos[j, 2]), 2]) - min(SST[c(Drange_pos[j, 1]:length(SST[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                }else{ # If the D range is contained within one year, range of d18Oc values is located between the lower and upper boundaries of D (+/- 1 SD)
                    Omod_SD[j] <- (max(d18Oc[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(d18Oc[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                    t_SD[j] <- length(Drange_pos[j, 1]:Drange_pos[j, 2]) / 2
                    gr_SD[j] <- (max(GR[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(GR[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                    Tmod_SD[j] <- (max(SST[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(SST[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                }
            }
            resmat <- cbind(resmat, Omod_SD, t_SD, gr_SD, Tmod_SD)
        }else{
            resmat <- cbind(resmat, matrix(rep(NA, 4 * length(resmat[,1])), ncol = 4)) # Add NA-filled columns if no MC error propagation is done
        }
        colnames(resmat) <- c("Dsam",
            "Osam",
            "Omod",
            "Residuals",
            "Time_of_year",
            "Instantaneous_growth_rate",
            "Modeled_temperature",
            "Omod_SD",
            "Time_of_Year_SD",
            "Instantaneous_growth_rate_SD",
            "Modeled_temperature_SD")

        if(plot == TRUE){
            plotdf <- as.data.frame(cbind(resmat, D_err, O_err))
            dev.new()
            ggplot2::ggplot(plotdf, ggplot2::aes(x = Dsam, y = Osam)) +
                ggplot2::geom_point() + # Plot original data
                ggplot2::geom_errorbar(ggplot2::aes(ymin = Osam - O_err, # Add error bars on measurement (1 SD)
                    ymax = Osam + O_err),
                    width = 100,
                    col = "black") +
                ggplot2::geom_errorbarh(ggplot2::aes(xmin = Dsam - D_err,
                    xmax = Dsam + D_err),
                    height = 0.05,
                    col = "black") +
                ggplot2::geom_point(ggplot2::aes(Dsam, Omod), col="red") + # Plot modeled d18Oc on top of data
                ggplot2::geom_errorbar(ggplot2::aes(ymin = Omod - Omod_SD, # Add error bars on model result (1 SD)
                    ymax = Omod + Omod_SD),
                    width = 100,
                    col = "red") +
                ggplot2::geom_errorbarh(ggplot2::aes(Dsam,
                    Omod,
                    xmin = Dsam - D_err,
                    xmax = Dsam + D_err),
                    height = 0.05,
                    col = "red")
        }
        return(resmat)
    }
}


text:
#' gglorenz: Plotting Lorenz Curve with ggplot2
#'
#' The package provides statistical transformations for plotting
#' empirical ordinary Lorenz curve and generalized Lorenz curve.
#'
#' @name gglorenz
#' @docType package
#' @rdname gglorenz-package
#' @author JJ Chen
NULL

max_stars_repo_name: cortinah/gglorenz

max_issues_repo_name: cortinah/gglorenz

max_forks_repo_name: cortinah/gglorenz


text:
# Testing the SVD empirically via Website
# Paul Sztorc

## Load Data ##
print("Loading data..")

tryCatch(expr=setwd("~/GitHub/Truthcoin/lib"), error=function(e) setwd(choose.dir(caption="Failed to set working directory automatically. Choose 'Truthcoin/lib' folder:")) )

DataFile <- "demo/input/Questions_VoteMatrix.csv"
Data <- read.csv(DataFile, stringsAsFactors= FALSE, row.names=1)

print("Load Complete.")
print(" ")

## Get VoteMatrix ##
print("Getting Votes..")

ToMatrix <- function(DF) {
  RowNames <- row.names(DF)
  DFn <- data.frame ( lapply( DF, as.numeric) ) # make all observations numbers
  Matrix <- as.matrix(DFn,ncol = ncol(DF)) # save in matrix format
  row.names(Matrix) <- RowNames  # restore row names (lost during as.numeric)
  return(Matrix)
}

# Remove header info
VoteData <- Data[-1:-4,]
VoteMatrix <- ToMatrix(VoteData)

print(VoteMatrix)
print(" ")

## Rescale ##
print("Rescaling the Scaled Decisions..")

# Scaled claims must become range(0,1)

Scaled <- Data["Qtype",] == "S" # which of these are scaled at all?

# Rescale
ScaleData <- Data[c("Min", "Max"),Scaled] # get the scales from the 'Qtype' row, and nothing else
ScaleMatrix <- ToMatrix(ScaleData)

RescaledVoteMatrix <- VoteMatrix # declare destination data
for(i in 1:ncol(ScaleMatrix)) { # for each scaled decision

  ThisQ <- colnames(ScaleMatrix)[i]
  ThisColumn <- RescaledVoteMatrix[ , colnames(RescaledVoteMatrix)==ThisQ ] # match the right column
  RescaledColumn <- (ThisColumn - ScaleMatrix["Min",i]) / (ScaleMatrix["Max",i] - ScaleMatrix["Min",i])  # "rescale"

  RescaledVoteMatrix[ , colnames(RescaledVoteMatrix)==ThisQ ] <- RescaledColumn # Overwrite
}

print("Rescale Complete.")
print(" ")

## Do Computations ##

# Load what we need
print("Loading SVD Consensus Mechanism..")
source("consensus/ConsensusMechanism.r")

# I've already rescaled, but we still need to pass the boolean - must fix this.
ScaleData <- matrix( c( rep(FALSE,ncol(RescaledVoteMatrix)),
              rep(0,ncol(RescaledVoteMatrix)),
              rep(1,ncol(RescaledVoteMatrix))), 3, byrow=TRUE, dimnames=list(c("Scaled","Min","Max"),colnames(RescaledVoteMatrix)) )

ScaleData[1,] <- Scaled

# Get the Resuls
print("Calculating Results..")
SvdResults <- Factory(RescaledVoteMatrix,Scales = ScaleData)

print("Writing Output .csvs..")
print(" ")

write.csv(SvdResults$Original, file="demo/output/OriginalVoteMatrix.csv")
write.csv(SvdResults$Agents, file="demo/output/agents.csv")
write.csv(SvdResults$Decisions, file="demo/output/decisions.csv")

print("Loading Plot Functions..")
source("consensus/PlotJ.r")

print("Making Plot..")
Plot <- PlotJ(RescaledVoteMatrix, Scales = ScaleData)

svg("demo/output/plot.svg",width = 8.5,height = 11)
Plot
dev.off()

print("Plot Complete.")

print("Done!")

max_stars_repo_name: endolith/Truthcoin

max_issues_repo_name: endolith/Truthcoin

max_forks_repo_name: endolith/Truthcoin


text:
F_statistic=52.01 # computed in example 11.4
F_statistic
t_statistic=-7.212
# both are calculated in 11.4
p=t_statistic^2
# t^2  equals the F value, to within round-off error.
p


text:
# check and install dependent packages in R
pkgs <-c('twitteR','ROAuth','httr','plyr','stringr','ggplot2','plotly')
for(p in pkgs) if(p %in% rownames(installed.packages()) == FALSE) {install.packages(p)}
for(p in pkgs) suppressPackageStartupMessages(library(p, quietly=TRUE, character.only=TRUE))

# Set API Keys.Below keys are sample only.Keys need to be generated using personal twitter account.
api_key <- "dOQFWD0RnmMV4eru4nUp8UG3C"
api_secret <- "0omNIMld9HrrBXLAKVXqaIdpUeeGCzlVr7G6aBCYsCHu3EMEYq"
access_token <- "84611119-3D0eyGI6hCfeto44FKjsT3p67smjYsCOIQiTCYi7v"
access_token_secret <- "DcNi9PzEE4tIkhaRpff4Qa7HUQZDstFSPcDKwIfZu8cRw"
setup_twitter_oauth(api_key, api_secret, access_token, access_token_secret)

# Grab latest tweets from CCI_India timeline
timeline_cciindia <- userTimeline('@CCI_India', n=1000)

#make data frame
df_cciindia <- do.call("rbind", lapply(timeline_cciindia, as.data.frame))
row.names(df_cciindia) <- NULL
View(df_cciindia)
tml_cciindia <- data.frame(df_cciindia$text, df_cciindia$created)
row.names(tml_cciindia) <- NULL
View(tml_cciindia)

# export the results into csv file
write.table(tml_cciindia, "c:/Users/USER/Desktop/iba/tweet_scrap/tml_cciindia.csv", sep=",",row.names = F)


text:
#' Convert json object to a dataframe
#' @param x list object
#' @param ... other parameters passed to \code{\link[tibble]{as_tibble}}:
#' @export
as_tibble.ppmembers <- function(x, ...) {

  # convert list to data frame
  mems <- x$members
  mems <- purrr::transpose(mems)
  mems <- purrr::map(mems, replace_nulls)
  mems <- purrr::simplify_all(mems)
  df <- tibble::as_tibble(mems)

  # fix column types
  fix_column_types(df)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmember <- function(x, ...) {
  # convert list to data frame
  mems <- replace_nulls(x)
  class(mems) <- "list"
  df <- tibble::as_tibble(mems)

  # fix column types
  fix_column_types(df)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmembers_new <- function(x, ...) {
  as_tibble.ppmembers(x, ...)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmembers_location <- function(x, ...) {
  df <- lapply(x, as_tibble.ppmember)
  do.call(rbind, df)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmembers_leaving <- function(x, ...) {
  as_tibble.ppmembers(x, ...)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmember_votes <- function(x, ...) {

  # convert list to data frame
  votes <- x$votes
  votes <- purrr::transpose(votes)
  votes <- purrr::map(votes, replace_nulls)
  votes <- purrr::simplify_all(votes)
  df <- tibble::as_tibble(votes)

  # fix column types
  fix_column_types(df)
}

#' @rdname as_tibble.ppmembers
#' @export
as_tibble.ppmember_bills <- function(x, ...) {

  # convert list to data frame
  bills <- x$bills
  bills <- purrr::transpose(bills)
  bills <- purrr::map(bills, replace_nulls)
  bills <- purrr::simplify_all(bills)
  df <- tibble::as_tibble(bills)

  # fix column types
  fix_column_types(df)
}

max_stars_repo_name: zamorarr/ppcongress

max_issues_repo_name: zamorarr/ppcongress

max_forks_repo_name: zamorarr/ppcongress


text:
#' Create a tidy stats data frame from an anova object
#'
#' \code{tidy_stats.anova} takes an anova object and tidies the output into a
#' data frame.
#'
#' @param model Output of \code{anova()}.
#' @param args Unused.
#'
#' @examples
#' # Regression example
#' ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
#' trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
#' group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
#' weight <- c(ctl, trt)
#'
#' model_0 <- lm(weight ~ 1)
#' model_1 <- lm(weight ~ group)
#'
#' tidy_stats(anova(model_0))
#' tidy_stats(anova(model_1))
#' tidy_stats(anova(model_0, model_1))
#' tidy_stats(anova(model_0, model_1, test = "Chisq"))
#'
#' # Logistic regression example
#' counts <- c(18,17,15,20,10,20,25,13,12)
#' outcome <- gl(3,1,9)
#' treatment <- gl(3,3)
#'
#' glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
#'
#' tidy_stats(anova(glm.D93))
#'
#' @export

tidy_stats.anova <- function(model, args = NULL) {

  # Convert model output to a data frame
  output <- tibble::as_data_frame(model)

  # Rename columns
  output <- rename_columns(output)

  # Set term information
  if (sum(stringr::str_detect(attr(model, "heading"), "Model [0-9]")) == 1) {
    output <- output %>%
      dplyr::mutate(
        term = stringr::str_split(attr(model, "heading")[2], "\n")[[1]],
        term_nr = 1:dplyr::n()
      )
  } else {
    output <- output %>%
      dplyr::mutate(
        term = rownames(model),
        term_nr = 1:dplyr::n()
      )
  }

  # Tidy stats
  output <- output %>%
    tidyr::gather("statistic", "value", -term, -term_nr) %>%
    dplyr::filter(!is.na(value)) %>%
    dplyr::arrange(term_nr)

  # Add method information
  output <- dplyr::mutate(output, method = "ANOVA Table")

  # If it was a chi-square test on an lm model, indicate this in a notes column
  if (sum(stringr::str_detect(names(model), "Pr\\(>Chi\\)")) == 1) {
    output <- dplyr::mutate(output, notes = "Chi-squared test")
  }

  return(output)
}


text:
require("shiny")
require("reshape")
require("ggplot2")
require("plotly")
require("edfReader")
require("eegkit")
require("stringr")

srate=256

eegHead.tabPanel<-function(profile){
	tabPanel("EEG Topoplot",div(id = "eggHead.container"
          , div(class = "plotlybars-wrapper"
            , div( class="plotlybars"
              , div(class="plotlybars-bar b1")
              , div(class="plotlybars-bar b2")
              , div(class="plotlybars-bar b3")
              , div(class="plotlybars-bar b4")
              , div(class="plotlybars-bar b5")
              , div(class="plotlybars-bar b6")
              , div(class="plotlybars-bar b7")
            )
            , div(class="plotlybars-text"
              , p("loading")
            )
          )
          ,plotlyOutput("eegHead.plot",width="640px",height="640px")
		)
      )
#	tabPanel("Box Plot", plotlyOutput("eegHead.plot",width="100%",height ="3000px" ))
}
eegHead.sidebarPanel<-function(profile,input,session){
	conditionalPanel(condition = "input.tabs1 == 'EEG Topoplot'",
					h3("View from top")
        #h4(paste(nrow(profile),"rows")),h5(title),h5(subtitle),
#				     selectInput("eggsignal","Select a EGG signal source",edfHeader.data()$sHeader$label[edfHeader.data()$sHeaders$transducerType=="EEG electrode"])
	)
}

#eegHead.prepareData<-function(input){
	eegHead.coor <-reactive({
		read.csv(fileLocate("eegcoord.csv"))
	})

	eegHead.data <- reactive({
		datafilename=paste(edfHeader.files()[input$edf.file],".Rdata",sep="")
		if(file.exists(datafilename)){
			load(datafilename)
		}
		else{
			snames=rownames(eegHead.coor())
			withProgress(message = 'Making plot', value = 0, {
			ss=lapply(edfSignal.data()[snames], function(x) x$signal[edfSignal.cfdRange()])
			ss2=lapply(ss,function(x) predict(eegsmooth(x,time=(0:(length(x)-1))/256),newdata=(0:(floor(length(x)/16)-1))/16))
			b<-do.call(cbind.data.frame, ss2)
			b=b[seq(1,nrow(b),64),]
			steps <- list()
			for(i in 1:nrow(b)){
				step <- list(args = list('visible', rep(FALSE, nrow(b))),
             	method = 'restyle')
				step$args[[2]][i] = TRUE
				steps[[i]] = step
			}
			coor=eegHead.coor()
			coor$name=rownames(coor)
			p<-plot_ly()
			for(ff in 1:nrow(b)){
				incProgress(1/nrow(b), detail = paste("Doing part", ff))
				coor$signal=unlist(b[ff,rownames(coor)])
				data.loess <- loess(signal ~ x * y, data = coor,control=loess.control(surface="direct"))
				xgrid <-  seq(-10, 10, .1)
				ygrid <-  seq(-10, 10, .1)
				data.fit <-  expand.grid(x = xgrid, y = ygrid)

				mtrx3d <-  predict(data.loess, newdata = data.fit)

				mtrx.melt <- melt(mtrx3d, id.vars = c('x', 'y'), measure.vars = 'signal')
				mtrx.melt$x <- as.numeric(str_sub(mtrx.melt$x, str_locate(mtrx.melt$x, '=')[1,1] + 1))
				mtrx.melt$y <- as.numeric(str_sub(mtrx.melt$y, str_locate(mtrx.melt$y, '=')[1,1] + 1))
				mtrx.melt$inC=(mtrx.melt$x^2+mtrx.melt$y^2)<=100
				mtrx.melt=mtrx.melt[mtrx.melt$inC,]

				vv=F
				if(ff==1){
  					vv=T
				}
#p<-p %>% add_trace(data=mtrx.melt,type = "contour",x=~x,y=~y,z=~value,connectgaps=F,visible=vv,name=paste("v =",ff),showlegend = FALSE,showscale=F,autocontour=F,contours=list(start=bmin,end=bmax))
				p<-p %>% add_trace(data=mtrx.melt,type = "contour",x=~x,y=~y,z=~value,connectgaps=F,visible=vv,name=paste("s =",ff*4-4),showlegend = FALSE,showscale=F)
			}
			})
			p=p %>% add_trace(data=coor,type="scatter",mode="markers+text", x=~x,y=~y,text=~name,textposition = "top right",marker=list(color="white",size=10,line = list(color ="black",width=2)))
			pr=p %>%
				  layout(sliders = list(list(active = 0,
                             steps = steps)))
			save(pr,file=datafilename)
		}
		pr
	})

#}

#eegHead.preparePlot<-function(input,output){
	output$eegHead.plot<-renderPlotly({
		eegHead.data()

	})

#}

max_stars_repo_name: chiajungchang/cfsindividual

max_issues_repo_name: chiajungchang/cfsindividual

max_forks_repo_name: chiajungchang/cfsindividual