venelin / PCMBase
Showing 2 of 53 files from the diff.
Other files ignored by Codecov
man/PCMMean.Rd has changed.
man/PCM.Rd has changed.
man/PCMCondVOU.Rd has changed.
man/PCMSim.Rd has changed.
man/PCMLmr.Rd has changed.
man/PCMAbCdEf.Rd has changed.
man/PCMCond.Rd has changed.
NAMESPACE has changed.
man/PCMLik.Rd has changed.
man/PCMTable.Rd has changed.
man/PCMInfo.Rd has changed.
DESCRIPTION has changed.
man/dataFig3.Rd has changed.
NEWS.md has changed.
man/PCMVar.Rd has changed.
man/PCMOptions.Rd has changed.

@@ -79,7 +79,7 @@
Loading
79 79
#' will result in tolerating some parameter values resulting in singular variance covariance
80 80
#' matrix of the transition distribution. Default 1e-4.}
81 81
#' \item{\code{PCMBase.Skip.Singular }}{A logical value indicating whether internal branches with
82 -
#' singular matrix V and shorter than \code{getOption("PCMBase.Threshold.Singular.Skip")}
82 +
#' singular matrix V and shorter than \code{getOption("PCMBase.Threshold.Skip.Singular")}
83 83
#'  should be skipped during likelihood calculation, adding their children
84 84
#' L,m,r values to their parent node. Default TRUE. Note, that setting this option to FALSE
85 85
#' may cause some models to stop working, e.g. the White model. Setting this option to FALSE

@@ -2216,7 +2216,7 @@
Loading
2216 2216
  list(nodes = nodes, positions = positions)
2217 2217
}
2218 2218
2219 -
#' Insert singleton nodes on chosen edges
2219 +
#' Insert tips or singleton nodes on chosen edges
2220 2220
#'
2221 2221
#' @param tree a phylo object
2222 2222
#' @param nodes an integer vector denoting the terminating nodes of the edges
@@ -2226,7 +2226,21 @@
Loading
2226 2226
#'  function several times with the longest position first and so on .
2227 2227
#' @param positions a positive numeric vector of the same length as nodes
2228 2228
#'  denoting the root-ward distances from nodes at which the singleton nodes
2229 -
#'  should be inserted.
2229 +
#'  should be inserted. For PCMTreeInsertTipsOrSingletons this can contains 0's and
2230 +
#'  is set by default to rep(0, length(nodes)).
2231 +
#' @param singleton (PCMTreeInsertTipsOrSingletons only) a logical indicating if a
2232 +
#' singleton node should be inserted and no tip node should be inserted.
2233 +
#' @param tipBranchLengths (PCMTreeInsertTipsOrSingletons only) positive numeric vector of the
2234 +
#' length of \code{nodes}, denoting the lengths of the new edges leading to tips.
2235 +
#' @param nodeLabels (PCMTreeInsertSingletons and PCMTreeInsertTipsOrSingletons) a
2236 +
#' character vector of the same length as \code{nodes}, indicating the names of
2237 +
#' the newly inserted nodes. These names are ignored where \code{positions} is 0. This
2238 +
#' argument is optional and default node labels will be assigned if this is not specified or set
2239 +
#' to NULL. If specified, then it should not contain node-labels already present in the tree.
2240 +
#' @param tipLabels (PCMTreeInsertTipsOrSingletons only) a character vector of the same length as
2241 +
#' \code{nodes} of the new tip-labels. This
2242 +
#' argument is optional and default tip labels will be assigned if this is not specified or set
2243 +
#' to NULL. If specified, then it should not contain tip-labels already present in the tree.
2230 2244
#' @param epoch a numeric indicating a distance from the root at which a
2231 2245
#' singleton node should be inserted in all lineages that are alive at that
2232 2246
#' time.
@@ -2236,8 +2250,8 @@
Loading
2236 2250
#' that this condition is checked only in `PCMTreeInsertSingletonsAtEpoch`.
2237 2251
#'
2238 2252
#' @importFrom ape bind.tree drop.tip
2239 -
#' @return a modified version of tree with inserted singleton nodes at the
2240 -
#'  specified locations
2253 +
#' @return a modified copy of tree.
2254 +
#'
2241 2255
#' @seealso \code{\link{PCMTreeEdgeTimes}} \code{\link{PCMTreeLocateEpochOnBranches}} \code{\link{PCMTreeLocateMidpointsOnBranches}}
2242 2256
#' @examples
2243 2257
#' set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
@@ -2301,6 +2315,50 @@
Loading
2301 2315
#' }
2302 2316
#' @export
2303 2317
PCMTreeInsertSingletons <- function(tree, nodes, positions) {
2318 +
  PCMTreeInsertTipsOrSingletons(tree, nodes, positions, singleton = TRUE)
2319 +
}
2320 +
2321 +
#' @describeIn PCMTreeInsertSingletons
2322 +
#'
2323 +
#' @export
2324 +
PCMTreeInsertSingletonsAtEpoch <- function(tree, epoch, minLength = 0.1) {
2325 +
2326 +
  tree <- PCMTree(tree)
2327 +
2328 +
  nodeTimes <- PCMTreeNodeTimes(tree)
2329 +
2330 +
  points <- PCMTreeLocateEpochOnBranches(tree, epoch)
2331 +
2332 +
  idxPoints <- sapply(seq_along(points$nodes), function(i) {
2333 +
    node <- points$nodes[i]
2334 +
    nodeTime <- nodeTimes[points$nodes[i]]
2335 +
    parentNode <- tree$edge[tree$edge[, 2] == points$nodes[i], 1]
2336 +
    parentTime <- nodeTimes[parentNode]
2337 +
    pos <- points$positions[i]
2338 +
2339 +
    (points$positions[i] > minLength) && (nodeTime - pos - parentTime) > minLength
2340 +
    })
2341 +
2342 +
  if(length(idxPoints) > 0) {
2343 +
    points$nodes <- points$nodes[idxPoints]
2344 +
    points$positions <- points$positions[idxPoints]
2345 +
    if(length(points$nodes) > 0) {
2346 +
      PCMTreeInsertSingletons(tree, points$nodes, points$positions)
2347 +
    } else {
2348 +
      tree
2349 +
    }
2350 +
  } else {
2351 +
    tree
2352 +
  }
2353 +
2354 +
}
2355 +
2356 +
#' @describeIn PCMTreeInsertSingletons
2357 +
#'
2358 +
#' @export
2359 +
PCMTreeInsertTipsOrSingletons <- function(
2360 +
  tree, nodes, positions = rep(0, length(nodes)),
2361 +
  singleton = FALSE, tipBranchLengths = 0.01, nodeLabels = NULL, tipLabels = NULL) {
2304 2362
2305 2363
  tree <- PCMTree(tree)
2306 2364
@@ -2315,10 +2373,46 @@
Loading
2315 2373
     max(nodes, na.rm = TRUE) > PCMTreeNumNodes(tree) ||
2316 2374
     min(nodes, na.rm = TRUE) < 1L) {
2317 2375
    stop(paste0(
2318 -
      "PCMTreeInsertSingletons:: Some of the nodes were NA or they could not",
2376 +
      "PCMTreeInsertTipsOrSingletons:: Some of the nodes were NA or they could not",
2319 2377
      "be matched against nodes in tree."))
2320 2378
  }
2321 2379
2380 +
  if(!is.null(nodeLabels)) {
2381 +
    if(!(is.character(nodeLabels) && length(nodeLabels) == length(nodes))) {
2382 +
      stop(paste0(
2383 +
        "PCMTreeInsertTipsOrSingletons:: nodeLabels should be NULL or a character ",
2384 +
        "vector of the same length as nodes"))
2385 +
    }
2386 +
    if(length(intersect(nodeLabels, originalLabels)) > 0) {
2387 +
      stop(paste0(
2388 +
        "PCMTreeInsertTipsOrSingletons:: some labels in nodeLabels are already used as labels in the tree."))
2389 +
    }
2390 +
  }
2391 +
2392 +
  if(!is.null(tipLabels)) {
2393 +
    if(!(is.character(tipLabels) && length(tipLabels) == length(nodes))) {
2394 +
      stop(paste0(
2395 +
        "PCMTreeInsertTipsOrSingletons:: tipLabels should be NULL or a character ",
2396 +
        "vector of the same length as nodes"))
2397 +
    }
2398 +
    if(length(intersect(tipLabels, originalLabels)) > 0) {
2399 +
      stop(paste0(
2400 +
        "PCMTreeInsertTipsOrSingletons:: some labels in tipLabels are already used as labels in the tree."))
2401 +
    }
2402 +
  }
2403 +
2404 +
  if(!is.numeric(tipBranchLengths)) {
2405 +
    stop("PCMTreeInsertTipsOrSingletons:: tipBranchLengths should be positive numeric of length 1 or equal to the length of nodes.")
2406 +
  }
2407 +
2408 +
  if(length(tipBranchLengths) == 1) {
2409 +
    tipBranchLengths <- rep(tipBranchLengths, length(nodes))
2410 +
  } else if(length(tipBranchLengths) != length(nodes)) {
2411 +
    stop("PCMTreeInsertTipsOrSingletons:: tipBranchLengths should be positive numeric of length 1 or equal to the length of nodes.")
2412 +
  } else if(isTRUE(any(tipBranchLengths) < 0)) {
2413 +
    stop("PCMTreeInsertTipsOrSingletons:: all entries in tipBranchLengths should be non-negative.")
2414 +
  }
2415 +
2322 2416
  PCMTreeSetLabels(tree, paste0("x_", PCMTreeGetLabels(tree)))
2323 2417
2324 2418
  names(originalLabels) <- PCMTreeGetLabels(tree)
@@ -2329,47 +2423,67 @@
Loading
2329 2423
  edgeTimes <- PCMTreeEdgeTimes(tree)
2330 2424
  rownames(edgeTimes) <- edgeNames
2331 2425
2332 -
  # which edges should be processed
2333 2426
  names(nodes) <- names(positions) <- PCMTreeGetLabels(tree)[nodes]
2334 -
  edgesToCut <- tree$edge[, 2] %in% nodes
2335 -
  edgesToCutNames <- edgeNames[edgesToCut]
2336 -
2337 -
  tipTree <- list(edge = matrix(c(2, 1), nrow = 1, ncol = 2),
2338 -
                  tip.label = "y_1",
2339 -
                  node.label = "y_2",
2340 -
                  edge.length = c(1.0),
2341 -
                  Nnode = 1)
2342 -
  class(tipTree) <- "phylo"
2343 -
2344 -
  for(edgeName in edgesToCutNames) {
2345 -
    # find the number of the ending node for the edge. DON'T USE CASHED LABELS
2346 -
    # HERE, but use PCMTreeGetLabels(tree)!
2347 -
    node <- match(edgeName, PCMTreeGetLabels(tree))
2348 -
2349 -
    # find the position (root-ward offset from node) where to cut
2350 -
2351 -
    # add, then drop the tipTree using ape's functions
2352 -
    tree <- bind.tree(tree, tipTree, node, positions[edgeName])
2353 -
    tree <- drop.tip(tree, "y_1", collapse.singles = FALSE)
2354 -
2355 -
    # inserted edge
2356 -
    edgeNameNew <- paste0(
2357 -
      "i_",
2358 -
      round(edgeTimes[edgeName, 2] - positions[edgeName], 2), "_", edgeName)
2359 -
2360 -
    tree$node.label[is.na(tree$node.label)] <- edgeNameNew
2361 -
2362 -
    # The edge leading to the new internal node take the same part as ITS
2363 -
    # DAUGHTER edge.
2364 -
    # If the edge we just cut was leading to a partition node, then this
2365 -
    # node is no more a partition node, because the newly inserted node is
2366 -
    # its parent and it has to becomes the partition node.
2367 -
    # Also the part should be renamed. The regime for the (now renamed) part is
2368 -
    # preserved.
2369 -
    idxPartNodeInPartRegimes <- match(edgeName, names(partRegimeWithXLabels))
2370 -
    if(!is.na(idxPartNodeInPartRegimes)) {
2371 -
      names(partRegimeWithXLabels)[idxPartNodeInPartRegimes] <- edgeNameNew
2427 +
2428 +
  # which edges should be processed in order of nodes
2429 +
  edgesToCutNames <- edgeNames[match(nodes, tree$edge[, 2])]
2430 +
2431 +
2432 +
  for(i in seq_along(edgesToCutNames)) {
2433 +
    edgeName <- edgesToCutNames[i]
2434 +
2435 +
    if(singleton && positions[edgeName] == 0) {
2436 +
      warning("PCMTreeInsertTipsOrSingletons:: Skipping ", edgeName, " since position is 0 and singleton is set to TRUE.")
2437 +
      next
2438 +
    } else {
2439 +
      # find the number of the ending node for the edge. DON'T USE CASHED LABELS
2440 +
      # HERE, but use PCMTreeGetLabels(tree)!
2441 +
      node <- match(edgeName, PCMTreeGetLabels(tree))
2442 +
2443 +
      # inserted edge
2444 +
      edgeNameNew <- if(!is.null(nodeLabels)) {
2445 +
        nodeLabels[i]
2446 +
      } else {
2447 +
        paste0(
2448 +
          "i_",
2449 +
          round(edgeTimes[edgeName, 2] - positions[edgeName], 2), "_", edgeName)
2450 +
      }
2451 +
2452 +
      tipLabelNew <- if(!is.null(tipLabels)) {
2453 +
        tipLabels[i]
2454 +
      } else {
2455 +
        paste0("t", "_", edgeNameNew)
2456 +
      }
2457 +
2458 +
      tipTree <- structure(
2459 +
        list(edge = matrix(c(2, 1), nrow = 1, ncol = 2),
2460 +
             tip.label = tipLabelNew,
2461 +
             edge.length = tipBranchLengths[i],
2462 +
             Nnode = 1),
2463 +
        class = "phylo")
2464 +
2465 +
      # bind the tipTree
2466 +
      tree <- bind.tree(tree, tipTree, node, positions[edgeName])
2467 +
      if(singleton) {
2468 +
        # drop the tip, without removing its singleton parent node
2469 +
        tree <- drop.tip(tree, tipLabelNew, collapse.singles = FALSE)
2470 +
      }
2471 +
2472 +
      tree$node.label[is.na(tree$node.label)] <- edgeNameNew
2473 +
2474 +
      # The edge leading to the new internal node take the same part as ITS
2475 +
      # DAUGHTER edge.
2476 +
      # If the edge we just cut was leading to a partition node, then this
2477 +
      # node is no more a partition node, because the newly inserted node is
2478 +
      # its parent and it has to becomes the partition node.
2479 +
      # Also the part should be renamed. The regime for the (now renamed) part is
2480 +
      # preserved.
2481 +
      idxPartNodeInPartRegimes <- match(edgeName, names(partRegimeWithXLabels))
2482 +
      if(!is.na(idxPartNodeInPartRegimes)) {
2483 +
        names(partRegimeWithXLabels)[idxPartNodeInPartRegimes] <- edgeNameNew
2484 +
      }
2372 2485
    }
2486 +
2373 2487
  }
2374 2488
2375 2489
  PCMTreeSetPartRegimes(
@@ -2386,40 +2500,6 @@
Loading
2386 2500
  tree
2387 2501
}
2388 2502
2389 -
#' @describeIn PCMTreeInsertSingletons
2390 -
#'
2391 -
#' @export
2392 -
PCMTreeInsertSingletonsAtEpoch <- function(tree, epoch, minLength = 0.1) {
2393 -
2394 -
  tree <- PCMTree(tree)
2395 -
2396 -
  nodeTimes <- PCMTreeNodeTimes(tree)
2397 -
2398 -
  points <- PCMTreeLocateEpochOnBranches(tree, epoch)
2399 -
2400 -
  idxPoints <- sapply(seq_along(points$nodes), function(i) {
2401 -
    node <- points$nodes[i]
2402 -
    nodeTime <- nodeTimes[points$nodes[i]]
2403 -
    parentNode <- tree$edge[tree$edge[, 2] == points$nodes[i], 1]
2404 -
    parentTime <- nodeTimes[parentNode]
2405 -
    pos <- points$positions[i]
2406 -
2407 -
    (points$positions[i] > minLength) && (nodeTime - pos - parentTime) > minLength
2408 -
    })
2409 -
2410 -
  if(length(idxPoints) > 0) {
2411 -
    points$nodes <- points$nodes[idxPoints]
2412 -
    points$positions <- points$positions[idxPoints]
2413 -
    if(length(points$nodes) > 0) {
2414 -
      PCMTreeInsertSingletons(tree, points$nodes, points$positions)
2415 -
    } else {
2416 -
      tree
2417 -
    }
2418 -
  } else {
2419 -
    tree
2420 -
  }
2421 -
2422 -
}
2423 2503
2424 2504
#' Find the nearest node to a given time from the root (epoch) on each lineage crossing this epoch
2425 2505
#' @param tree a phylo
@@ -2532,7 +2612,17 @@
Loading
2532 2612
        node = N+1L,
2533 2613
        regime = tree$part.regime[PCMTreeGetPartsForNodes(tree, N+1L)]))
2534 2614
2535 -
    plotTree <- ggtree::`%<+%`(ggtree::ggtree(tree, ...), data)
2615 +
    # Try to prevent plotting errors due to bugs in the function ladderize in ape package
2616 +
    # only changing default values here.
2617 +
    listParams <- c(list(tr = tree), list(...))
2618 +
    if(!"ladderize" %in% names(listParams)) {
2619 +
      listParams[["ladderize"]] <- FALSE
2620 +
    }
2621 +
    if(!"right" %in% names(listParams)) {
2622 +
      listParams[["right"]] <- TRUE
2623 +
    }
2624 +
2625 +
    plotTree <- ggtree::`%<+%`(do.call(ggtree::ggtree, listParams), data)
2536 2626
2537 2627
    plotTree + aes(color = regime) +
2538 2628
      scale_color_manual(name = "regime", values = palette)
Files Coverage
R 35.71%
Project Totals (16 files) 35.71%
1
comment: false
2

3
ignore:
4
  - R/AutoGeneratedModelTypes.R
5

6
coverage:
7
  status:
8
    project:
9
      default:
10
        target: auto
11
        threshold: 1%
12
    patch:
13
      default:
14
        target: auto
15
        threshold: 1%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading