https://github.com/cran/ape
Raw File
Tip revision: 24545d71f754d68573d0c28b11bd23acc96e322b authored by Emmanuel Paradis on 20 August 2004, 00:00:00 UTC
version 1.2-5
Tip revision: 24545d7
collapsed.intervals.R
### collapsed.intervals.R  (2002-09-12)
###
###     Collapsed coalescent intervals (e.g. for the skyline plot)
###
### Copyright 2002 Korbinian Strimmer <strimmer@stat.uni-muenchen.de>
###
### This file is part of the `ape' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program 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 this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

# construct collapsed intervals from coalescent intervals
collapsed.intervals <- function(ci, epsilon=0.0)
{
  if (class(ci) != "coalescentIntervals")
    stop("object \"ci\" is not of class \"coalescentIntervals\"")

  sz <- ci$interval.length
  lsz <- length(sz)
  idx <- c <- 1:lsz

  p <- 1
  w <- 0
  
  # starting from tips collapes intervals
  # until total size is >= epsilon
  for (i in 1:lsz)
  {
    idx[[i]] <- p
    w <- w + sz[[i]]
    if (w >= epsilon)
    {
      p <- p+1
      w <- 0
    }
  }

  # if last interval is smaller than epsilon merge
  # with second last interval
  lastInterval <- idx==p
  if ( sum(sz[lastInterval]) < epsilon )
  {
    p <- p-1
    idx[lastInterval] <- p
  }

  obj <- list(
     lineages=ci$lineages,
     interval.length=ci$interval.length,
     collapsed.interval=idx, # collapsed intervals (via reference)
     interval.count=ci$interval.count,
     collapsed.interval.count = idx[[ci$interval.count]],
     total.depth =ci$total.depth,
     epsilon = epsilon
    )
  class(obj) <- "collapsedIntervals"
  
  return(obj)
}
back to top