ArtPoon / kamphir

Phylogenetic inference using a tree-shape kernel in an Approximate Bayesian Computation framework
BSD 3-Clause "New" or "Revised" License
6 stars 2 forks source link

Bug in simulate.binary.dated.tree.0 [known issue for rcolgem] #6

Closed ArtPoon closed 9 years ago

ArtPoon commented 9 years ago
Called from: FUN(1:10[[5L]], ...)
debug: if (sum(is.nan(A)) > 0) A <- A0
Error in phy$edge.length : $ operator is invalid for atomic vectors
Calls: simulate.binary.dated.tree.fgy -> lapply -> FUN -> write.tree -> .write.tree2
debug: if (nExtant > 1) {
    mstates[isExtant, ] <- t(t(Q) %*% t(mstates[isExtant, ]))
    mstates[isExtant, ] <- abs(mstates[isExtant, ])/rowSums(abs(mstates[isExtant, 
        ]))
    A <- colSums(mstates[isExtant, ])
} else {
    mstates[isExtant, ] <- t(t(Q) %*% mstates[isExtant, ])
    mstates[isExtant, ] <- abs(mstates[isExtant, ])/sum(abs(mstates[isExtant, 
        ]))
    A <- mstates[isExtant, ]
}
debug: mstates[isExtant, ] <- t(t(Q) %*% t(mstates[isExtant, ]))
debug: mstates[isExtant, ] <- abs(mstates[isExtant, ])/rowSums(abs(mstates[isExtant, 
    ]))
debug: A <- colSums(mstates[isExtant, ])
debug: if (isSampleEvent[ih + 1]) {
    sat_h1 <- sampled.at.h(h1)
    isExtant[sat_h1] <- TRUE
    heights[sat_h1] <- h1
} else {
    .F <- fgy$.F
    .G <- fgy$.G
    .Y <- fgy$.Y
    a <- A/.Y
    extantLines <- which(isExtant)
    tryCatch({
        .lambdamat <- (t(t(a)) %*% a) * .F
        kl <- sample.int(m^2, size = 1, prob = as.vector(.lambdamat))
        k <- 1 + ((kl - 1)%%m)
        l <- 1 + floor((kl - 1)/m)
        probstates <- mstates[extantLines, ]
        u_i <- sample.int(nExtant, size = 1, prob = probstates[, 
            k])
        probstates[u_i, ] <- 0
        u <- extantLines[u_i]
        v <- sample(extantLines, size = 1, prob = probstates[, 
            l])
    }, error = function(e) browser())
    ustates[u, ] <- mstates[u, ]
    ustates[v, ] <- mstates[v, ]
    a_u <- pmin(1, mstates[u, ]/.Y)
    a_v <- pmin(1, mstates[v, ]/.Y)
    lambda_uv <- ((a_u) %*% t(a_v)) * .F + ((a_v) %*% t(a_u)) * 
        .F
    palpha <- rowSums(lambda_uv)/sum(lambda_uv)
    alpha <- lineageCounter
    lineageCounter <- lineageCounter + 1
    isExtant[alpha] <- TRUE
    isExtant[u] <- FALSE
    isExtant[v] <- FALSE
    lstates[alpha, ] = (mstates[alpha, ] <- palpha)
    heights[alpha] <- h1
    uv <- c(u, v)
    inEdgeMap[uv] <- alpha
    outEdgeMap[alpha, ] <- uv
    parent[uv] <- alpha
    parentheights[uv] <- h1
    daughters[alpha, ] <- uv
    edge[u, ] <- c(alpha, u)
    edge.length[u] <- h1 - heights[u]
    edge[v, ] <- c(alpha, v)
    edge.length[v] <- h1 - heights[v]
}
debug: .F <- fgy$.F
debug: .G <- fgy$.G
debug: .Y <- fgy$.Y
debug: a <- A/.Y
debug: extantLines <- which(isExtant)
debug: tryCatch({
    .lambdamat <- (t(t(a)) %*% a) * .F
    kl <- sample.int(m^2, size = 1, prob = as.vector(.lambdamat))
    k <- 1 + ((kl - 1)%%m)
    l <- 1 + floor((kl - 1)/m)
    probstates <- mstates[extantLines, ]
    u_i <- sample.int(nExtant, size = 1, prob = probstates[, 
        k])
    probstates[u_i, ] <- 0
    u <- extantLines[u_i]
    v <- sample(extantLines, size = 1, prob = probstates[, l])
}, error = function(e) browser())
debug: .lambdamat <- (t(t(a)) %*% a) * .F
debug: kl <- sample.int(m^2, size = 1, prob = as.vector(.lambdamat))
Called from: value[[3L]](cond)
debug: ustates[u, ] <- mstates[u, ]
debug: ustates[v, ] <- mstates[v, ]
debug: a_u <- pmin(1, mstates[u, ]/.Y)
debug: a_v <- pmin(1, mstates[v, ]/.Y)
debug: lambda_uv <- ((a_u) %*% t(a_v)) * .F + ((a_v) %*% t(a_u)) * .F
debug: palpha <- rowSums(lambda_uv)/sum(lambda_uv)
debug: alpha <- lineageCounter
debug: lineageCounter <- lineageCounter + 1
debug: isExtant[alpha] <- TRUE
debug: isExtant[u] <- FALSE
debug: isExtant[v] <- FALSE
debug: lstates[alpha, ] = (mstates[alpha, ] <- palpha)
debug: heights[alpha] <- h1
debug: uv <- c(u, v)
debug: inEdgeMap[uv] <- alpha
debug: outEdgeMap[alpha, ] <- uv
debug: parent[uv] <- alpha
debug: parentheights[uv] <- h1
debug: daughters[alpha, ] <- uv
debug: edge[u, ] <- c(alpha, u)
debug: edge.length[u] <- h1 - heights[u]
debug: edge[v, ] <- c(alpha, v)
debug: edge.length[v] <- h1 - heights[v]
Error in kids[[parent[i]]] : attempt to select more than one element
Calls: simulate.binary.dated.tree.fgy ... lapply -> FUN -> read.tree -> write.tree -> .write.tree2
ArtPoon commented 9 years ago

Won't fix.