Commit 9ad2fc1c authored by Bastian Schiffthaler's avatar Bastian Schiffthaler

Reversed inf LOD change as that is currently bugged

parent b227bb2e
# This file was generated by Rcpp::compileAttributes
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
GET_RF_MAT_NO_LOD <- function(seqnum, nmrk, CC, CR, RC, RR, minLOD, maxRF) {
......
......@@ -126,7 +126,8 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
fun.order = NULL, phase.cores = 4,
ripple.cores = 1, verbosity = NULL, max.dist = Inf,
ws = 4, increase.every = 4, max.tries = 10,
min.tries = 0, ...)
min.tries = 0, seeds = NULL, optimize = "likelihood",
...)
{
#TODO: error checks...
#Create initial set of batches
......@@ -141,9 +142,16 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
LGs <- list()
#The first batch is run in full again to get all necessary data (phases etc.)
tryCatch({
LG <- map(make.seq(get(input.seq$twopt), batches[[1]],
twopt = input.seq$twopt), phase.cores = phase.cores,
verbosity = verbosity)
if(is.null(seeds))
{
LG <- map(make.seq(get(input.seq$twopt), batches[[1]],
twopt = input.seq$twopt), phase.cores = phase.cores,
verbosity = verbosity)
} else {
LG <- seeded.map(make.seq(get(input.seq$twopt), batches[[1]],
twopt = input.seq$twopt), phase.cores = phase.cores,
verbosity = verbosity, seeds = seeds)
}
}, error = function(e) {
warning("Error during initial map calculation.",
"Trying to fix by reordering...")
......@@ -152,7 +160,7 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
LG <- ripple_ord(input.seq = s,ws = ws, phase.cores = phase.cores,
ripple.cores = ripple.cores, method = "one",
no_reverse = TRUE, verbosity = verbosity, start = 1,
batches = batches)
batches = batches, optimize = optimize)
if(LG$seq.like == -Inf)
stop("Could not fix issue. You need to reorder or provide more",
" informative markers")
......@@ -172,7 +180,7 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
if(round %% increase.every == 0) increment <- increment + 1
LG <- fun.order(LG, ripple.cores = ripple.cores, start = 1,
verbosity = verbosity, batches = batches,
ws = ws + increment, ...)
ws = ws + increment, optimize = optimize, ...)
round <- round + 1
}
LGs[[1]] <- LG
......@@ -202,7 +210,8 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
LG <- ripple_ord(input.seq = s,ws = ws, phase.cores = phase.cores,
ripple.cores = ripple.cores, method = "one",
no_reverse = TRUE, verbosity = verbosity,
start = overlap + 2, batches = batches)
start = overlap + 2, batches = batches,
optimize = optimize)
if(LG$seq.like == -Inf)
stop("Could not fix issue. You need to reorder or provide more",
" informative markers")
......@@ -222,7 +231,7 @@ map_overlapping_batches <- function(input.seq, size = 50, overlap = 15,
if(round %% increase.every == 0) increment <- increment + 1
LG <- fun.order(LG, ripple.cores = ripple.cores, start=overlap+2,
verbosity = verbosity, batches = batches,
ws = ws + increment, ...)
ws = ws + increment, optimize = optimize, ...)
round <- round + 1
}
}
......
......@@ -129,7 +129,7 @@ generate_one <- function(input.seq, p, ws, no_reverse)
ripple_window <- function(input.seq, ws=4, tol=10E-4, phase.cores = 4,
ripple.cores = 4, start = 1, verbosity = NULL,
type = "one", n = NULL, pref = NULL,
no_reverse = TRUE) {
no_reverse = TRUE, optimize = "likelihood") {
## checking for correct objects
if(!any(class(input.seq)=="sequence")) {
......@@ -147,10 +147,10 @@ ripple_window <- function(input.seq, ws=4, tol=10E-4, phase.cores = 4,
## allocate variables
rf.init <- rep(NA,len-1)
phase <- rep(NA,len-1)
tot <- prod(1:ws)
best.ord.phase <- matrix(NA,tot,len-1)
best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
all.data <- list()
# tot <- prod(1:ws)
# best.ord.phase <- matrix(NA,tot,len-1)
# best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
# all.data <- list()
## gather two-point information
list.init <- phases(input.seq)
......@@ -167,45 +167,106 @@ ripple_window <- function(input.seq, ws=4, tol=10E-4, phase.cores = 4,
if(type == "rand") all.ord <- generate_rand(input.seq, p, ws, n, pref)
if(type == "one") all.ord <- generate_one(input.seq, p, ws, no_reverse)
poss <- mclapply(1:nrow(all.ord), mc.allow.recursive = TRUE,
mc.cores = ripple.cores, function(i){
if("position" %in% verbosity)
{
message("Trying order ",i," of ",nrow(all.ord),
" for start position ",p)
}
mp <- list(seq.like = -Inf)
tryCatch({
if(start > 1)
{
seeds <- input.seq$seq.phases[1:(start - 1)]
mp <- seeded.map(make.seq(get(input.seq$twopt),
all.ord[i,],
twopt = input.seq$twopt),
phase.cores = phase.cores,
verbosity = verbosity,
seeds = seeds)
}
else
if(optimize == "likelihood" || optimize == "size")
{
poss <- mclapply(1:nrow(all.ord), mc.allow.recursive = TRUE,
mc.cores = ripple.cores, function(i){
if("position" %in% verbosity)
{
mp <- map(make.seq(get(input.seq$twopt), all.ord[i,],
twopt = input.seq$twopt),
phase.cores = phase.cores,
verbosity = verbosity)
message("Trying order ",i," of ",nrow(all.ord),
" for start position ",p)
}
}, error = function(e){},
finally = {
return(mp)
mp <- list(seq.like = -Inf, seq.rf = -1)
tryCatch({
if(start > 1)
{
seeds <- input.seq$seq.phases[1:(start - 1)]
mp <- seeded.map(make.seq(get(input.seq$twopt),
all.ord[i,],
twopt = input.seq$twopt),
phase.cores = phase.cores,
verbosity = verbosity,
seeds = seeds)
}
else
{
mp <- map(make.seq(get(input.seq$twopt), all.ord[i,],
twopt = input.seq$twopt),
phase.cores = phase.cores,
verbosity = verbosity)
}
}, error = function(e){},
finally = {
return(mp)
})
})
})
} else if(optimize == "count") {
all.ord <- rbind(all.ord, input.seq$seq.num)
COUNT<-function(X, sequence){ ## See eq. 1 on the paper (Van Os et al., 2005)
return(sum(diag(X[sequence[-length(sequence)],sequence[-1]]),na.rm=TRUE))
}
r<-get_mat_rf_out(input.seq, LOD=FALSE, max.rf=0.5, min.LOD=0)
r[is.na(r)]<-0.5
diag(r)<-0
X<-r*get(input.seq$data.name, pos=1)$n.ind
count <- mclapply(1:nrow(all.ord), mc.allow.recursive = TRUE,
mc.cores = ripple.cores, function(i){
if("position" %in% verbosity)
{
message("Trying order ",i," of ",nrow(all.ord),
" for start position ",p)
}
se <- match(colnames(get(input.seq$data.name)$geno)[all.ord[i,]],
colnames(X))
return(COUNT(X, se))
})
}
if(optimize == "likelihood"){
best <- which.max(sapply(poss,"[[","seq.like"))
if(poss[[best]]$seq.like > input.seq$seq.like)
{
return(poss[[best]])
} else {
return(input.seq)
}
}
if(optimize == "size"){
# Check which map is shortest. If there are errors in the model, return NULL
tryCatch({
best <- which.min(sapply(poss,function(f){
if(length(f$seq.rf) < 2 | any(f$seq.rf < 0 | f$seq.rf > 0.5))
return(NULL)
sum(kosambi(f$seq.rf))
}))
}, error = function(e){
save(input.seq, poss, file = "~/bad_poss.RData")
}, finally = {})
# In case all models have errors, just return the input.
if(length(poss[[best]]$seq.rf) < 2 | any(poss[[best]]$seq.rf < 0 | poss[[best]]$seq.rf > 0.5))
return(input.seq)
if(sum(kosambi(poss[[best]]$seq.rf)) < sum(kosambi(input.seq$seq.rf)))
{
return(poss[[best]])
} else {
return(input.seq)
}
}
best <- which.max(sapply(poss,"[[","seq.like"))
if(poss[[best]]$seq.like > input.seq$seq.like)
if(optimize == "count")
{
return(poss[[best]])
} else {
return(input.seq)
best.ord <- which.min(unlist(count))
print(all.ord[best.ord,])
print(all.ord[nrow(all.ord),])
#print(unlist(count))
mp <- map(make.seq(get(input.seq$twopt), all.ord[best.ord, ],
twopt = input.seq$twopt))
return(mp)
}
}
......@@ -256,7 +317,7 @@ ripple_window <- function(input.seq, ws=4, tol=10E-4, phase.cores = 4,
ripple_ord <- function(input.seq, ws = 4, tol = 10E-4, phase.cores = 4,
ripple.cores = 4, method = "one", n = NULL,
pref = "neutral", start = 1, verbosity = NULL,
batches = NULL, no_reverse = TRUE)
batches = NULL, no_reverse = TRUE, optimize = "likelihood")
{
LG <- input.seq
if(start + ws > length(input.seq$seq.num)) return(LG)
......@@ -269,7 +330,7 @@ ripple_ord <- function(input.seq, ws = 4, tol = 10E-4, phase.cores = 4,
phase.cores = phase.cores, ripple.cores = ripple.cores,
start = i, verbosity = verbosity,
no_reverse = no_reverse, type = method, n = n,
pref = pref)
pref = pref, optimize = optimize)
toc <- Sys.time()
timings <- c(timings, as.numeric(difftime(toc, tic, units = "secs")))
if("time" %in% verbosity)
......
// This file was generated by Rcpp::compileAttributes
// Generated by using Rcpp::compileAttributes() -> do not edit by hand
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#include <RcppArmadillo.h>
......@@ -10,8 +10,8 @@ using namespace Rcpp;
SEXP GET_RF_MAT_NO_LOD(SEXP seqnum, SEXP nmrk, SEXP CC, SEXP CR, SEXP RC, SEXP RR, SEXP minLOD, SEXP maxRF);
RcppExport SEXP onemap_GET_RF_MAT_NO_LOD(SEXP seqnumSEXP, SEXP nmrkSEXP, SEXP CCSEXP, SEXP CRSEXP, SEXP RCSEXP, SEXP RRSEXP, SEXP minLODSEXP, SEXP maxRFSEXP) {
BEGIN_RCPP
Rcpp::RObject __result;
Rcpp::RNGScope __rngScope;
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< SEXP >::type seqnum(seqnumSEXP);
Rcpp::traits::input_parameter< SEXP >::type nmrk(nmrkSEXP);
Rcpp::traits::input_parameter< SEXP >::type CC(CCSEXP);
......@@ -20,30 +20,30 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< SEXP >::type RR(RRSEXP);
Rcpp::traits::input_parameter< SEXP >::type minLOD(minLODSEXP);
Rcpp::traits::input_parameter< SEXP >::type maxRF(maxRFSEXP);
__result = Rcpp::wrap(GET_RF_MAT_NO_LOD(seqnum, nmrk, CC, CR, RC, RR, minLOD, maxRF));
return __result;
rcpp_result_gen = Rcpp::wrap(GET_RF_MAT_NO_LOD(seqnum, nmrk, CC, CR, RC, RR, minLOD, maxRF));
return rcpp_result_gen;
END_RCPP
}
// READ_OUTCROSS
SEXP READ_OUTCROSS(SEXP file);
RcppExport SEXP onemap_READ_OUTCROSS(SEXP fileSEXP) {
BEGIN_RCPP
Rcpp::RObject __result;
Rcpp::RNGScope __rngScope;
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< SEXP >::type file(fileSEXP);
__result = Rcpp::wrap(READ_OUTCROSS(file));
return __result;
rcpp_result_gen = Rcpp::wrap(READ_OUTCROSS(file));
return rcpp_result_gen;
END_RCPP
}
// CCOUNT
SEXP CCOUNT(SEXP X, SEXP sequence);
RcppExport SEXP onemap_CCOUNT(SEXP XSEXP, SEXP sequenceSEXP) {
BEGIN_RCPP
Rcpp::RObject __result;
Rcpp::RNGScope __rngScope;
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< SEXP >::type X(XSEXP);
Rcpp::traits::input_parameter< SEXP >::type sequence(sequenceSEXP);
__result = Rcpp::wrap(CCOUNT(X, sequence));
return __result;
rcpp_result_gen = Rcpp::wrap(CCOUNT(X, sequence));
return rcpp_result_gen;
END_RCPP
}
......@@ -356,6 +356,9 @@ Rcpp::List est_rf_out(Rcpp::NumericVector geno,
}
for(int k=0; k < 4; k++)
{
if(r[k] < rf_TOL_min) r[k]=rf_TOL_min;
if(r[k] > rf_TOL_max) r[k]=rf_TOL_max;
/*
if(r[k] <= arma::datum::eps) {
r[k] = 0;
r[k+4] = R_PosInf;
......@@ -363,7 +366,7 @@ Rcpp::List est_rf_out(Rcpp::NumericVector geno,
if(r[k] >= (1 - arma::datum::eps)) {
r[k]= 1;
r[k+4]=R_NegInf;
}
} */
}
if(mrk < 0)
{
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment