从一组集合中查找所有不重叠的集合

4
我的问题:需要从一组集合中找到所有不重叠的集合。
背景:我正在使用比较进化方法研究鸟类特征演化。我有一棵包含大约300个物种的树。这棵树可以被分成子类(即子树)。如果两个子类没有共同的物种,它们是独立的。我正在寻找一种算法(如果可能的话,还需要一个R实现),它将找到所有可能的子类分区,每个子类都有大于10个生物分类单元,并且所有子类都是独立的。每个子类都可以被视为一个集合,当两个子类是独立的(没有共同的物种)时,这些子类就是不相交的集合。
希望这很清楚,希望有人能够提供帮助。
以下代码生成一个示例数据集。其中,子类是一个包含所有可能子类(集合)的列表,我想从中采样X个不相交的集合,而集合长度为Y。
###################################
# Example Dataset
###################################

    library(ape)
    library(phangorn)
    library(TreeSim)
    library(phytools)

    ##simulate a tree

    n.taxa <- 300
    tree <- sim.bd.taxa(n.taxa,1,lambda=.5,mu=0)[[1]][[1]]
    tree$tip.label <- seq(n.taxa)

    ##extract all monophyletic subclades

    get.all.subclades <- function(tree){
    tmp <- vector("list")
    nodes <- sort(unique(tree$edge[,1]))
    i <- 282
    for(i in 1:length(nodes)){
    x <- Descendants(tree,nodes[i],type="tips")[[1]]
    tmp[[i]] <- tree$tip.label[x]
    }
    tmp
    }
    tmp <- get.all.subclades(tree)

    ##set bounds on the maximum and mininum number of tips of the subclades to include

    min.subclade.n.tip <- 10
    max.subclade.n.tip <- 40


    ##function to replace trees of tip length exceeding max and min with NA

    replace.trees <- function(x, min, max){
    if(length(x) >= min & length(x)<= max) x else NA
    }


    #apply testNtip across all the subclades

    tmp2 <- lapply(tmp, replace.trees, min = min.subclade.n.tip, max = max.subclade.n.tip)

    ##remove elements from list with NA, 
    ##all remaining elements are subclades with number of tips between 
##min.subclade.n.tip and max.subclade.n.tip

    subclades <- tmp2[!is.na(tmp2)]

    names(subclades) <- seq(length(subclades))

那么你将每个子树视为一个子集,你想找到一组子集,使得每个元素恰好被包含一次;或者你不介意有些元素被遗漏吗? - Tom
2
如果您指的是“恰好一次”解释,那么这是精确覆盖问题。它是NP完全的,这意味着您不会比简单地回溯和测试所有可能性做得更好。 - ffao
1
我本来想说,这听起来像是一个需要使用 http://en.wikipedia.org/wiki/Knuth%27s_Algorithm_X 的问题。 - Tom
1
如果您的数据集不是非常大,那么首先收集所有 length(clade.x) >=10,然后在该子集的所有组合上运行 intersect - Carl Witthoft
同意@CarlWitthoft的观点,你可以像这样做:i <- 3; j <- 4; length(intersect(A[[i]],A[[j]]))>0 - Stéphane Laurent
显示剩余2条评论
2个回答

2
这是一个测试每对列表元素是否存在零重叠并提取所有非重叠对索引的示例。
findDisjointPairs <- function(X) {
    ## Form a 2-column matrix enumerating all pairwise combos of X's elements
    ij <- t(combn(length(X),2))    
    ## A function that tests for zero overlap between a pair of vectors
    areDisjoint <- function(i, j) length(intersect(X[[i]], X[[j]])) == 0     
    ## Use mapply to test for overlap between each pair and extract indices 
    ## of pairs with no matches
    ij[mapply(areDisjoint, ij[,1], ij[,2]),]
}

## Make some reproducible data and test the function on it
set.seed(1)
A <- replicate(sample(letters, 5), n=5, simplify=FALSE)    
findDisjointPairs(A)
#      [,1] [,2]
# [1,]    1    2
# [2,]    1    4
# [3,]    1    5

他仍然需要根据所需的元素大小进行过滤,但除此之外,这看起来非常不错。 - Carl Witthoft
1
@CarlWitthoft -- 确实。不过这很简单,我会让OP在最合适的地方将其添加到他/她的工作流程中。 (可以在areDisjoint函数的主体中仅添加一行length(X [[i]]> = 10&length(X [[j]])> = 10&,但我不知道是否需要优化此函数以提高速度,如果需要,则会做更复杂的事情。) - Josh O'Brien
嗨@JoshO'Brien。感谢您的评论和函数,它运行得很好。现在的问题是从集合的宇宙中找到所有可能的组合。当要计算的可能组合数量变大时,例如choose(50,20),我的计算机没有足够的内存。有没有办法从可能组合的非常大的列表中随机抽样组合?类似于combn(50,20)。如果我能做到这一点,那么我只需要随机选择集合的组合,直到循环累积了足够数量的独立亚支系的随机集合,以供我的分析使用。谢谢! - user1322491
如果我理解你的问题,以下是一个能帮到你的想法。replicate(5, sample(1:50, size=20)) (这个例子将给你五个包含从1到50之间选中的20个整数的集合,然后你可以使用它们作为索引来访问你的子集列表。) - Josh O'Brien

1
以下是一些可能有用的函数。
第一个函数计算列表集合的所有可能不相交的集合组合。
我使用“collection”而不是“partition”,因为集合的集合不一定覆盖宇宙(即所有集合的并集)。
该算法是递归的,并且仅适用于少量可能的集合。这并不一定意味着它无法处理大量的集合列表,因为该函数在每次迭代中都会删除相交的集合。
如果代码不清楚,请询问,我将添加注释。
输入必须是一个命名列表,结果将是一个包含集合名称的字符向量的集合列表。
DisjointCollectionsNotContainingX <- function(L, branch=character(0), x=numeric(0))
{
    filter <- vapply(L, function(y) length(intersect(x, y))==0, logical(1))

    L <- L[filter]

    result <- list(branch)

    for( i in seq_along(L) )
    {
        result <- c(result, Recall(L=L[-(1:i)], branch=c(branch, names(L)[i]), x=union(x, L[[i]])))
    }

    result
}

这只是一个隐藏辅助参数的包装器:
DisjointCollections <- function(L) DisjointCollectionsNotContainingX(L=L)

下一个函数可以用来验证给定的集合列表是否不重叠且“最大”。对于每个集合,它将测试:
1. 所有集合实际上是不相交的 2. 添加另一个集合会导致非不相交集合或现有集合。
ValidateDC <- function(L, DC)
{
    for( collection in DC )
    {
        for( i in seq_along(collection) )
        {
            others <- Reduce(f=union, x=L[collection[-i]])

            if( length(intersect(L[collection[i]], others)) > 0 ) return(FALSE)
        }

        elements <- Reduce(f=union, x=L[collection])

        for( k in seq_along(L) ) if( ! (names(L)[k] %in% collection) )
        {
            if( length(intersect(elements, L[[k]])) == 0 )
            {
                check <- vapply(DC, function(z) setequal(c(collection, names(L)[k]), z), logical(1))

                if( ! any(check) ) return(FALSE)
            }
        }
    }

    TRUE
}

例子:

L <- list(A=c(1,2,3), B=c(3,4), C=c(5,6), D=c(6,7,8))

> ValidateDC(L,DisjointCollections(L))
[1] TRUE

> DisjointCollections(L)
[[1]]
character(0)

[[2]]
[1] "A"

[[3]]
[1] "A" "C"

[[4]]
[1] "A" "D"

[[5]]
[1] "B"

[[6]]
[1] "B" "C"

[[7]]
[1] "B" "D"

[[8]]
[1] "C"

[[9]]
[1] "D"

请注意,同时包含AB的集合不会显示,因为它们存在非空交集。同时包含CD的集合也不会出现。其他情况都可以。
注意:空集合character(0)始终是有效的组合。
创建所有可能的不相交集合后,您可以应用任何过滤器以继续进行。

编辑:

  1. 我已经从第一个函数中删除了代码行 if( length(L)==0 ) return(list(branch)),这行代码不是必要的。

  2. 性能:如果集合之间存在较大重叠,则该函数运行速度较快。例如:

    set.seed(1)

    L <- lapply(1:50, function(.)sample(x=1200, size=20))

    names(L) <- c(LETTERS, letters)[1:50]

    system.time(DC <- DisjointCollections(L))

结果:

#   user  system elapsed 
#   9.91    0.00    9.92

发现的收藏总数:
> length(DC)
[1] 121791

感谢提供这些函数。我试用了它们,对于小于15的一小组集合,它们运行得非常好。但是当我尝试在长度大于40的大型列表上使用时,R崩溃了。我可能会扩大我的分析范围,并需要一个能够处理长度大于50的集合列表的函数。我应该澄清一下,我最感兴趣的是从长度大于X的不相交集合列表中进行随机抽样。X可能会是8或更多。我承认我不太理解您的Recall()函数是如何工作的,但是否有任何方法修改它以仅对长度大于X的不相交集合进行抽样?谢谢! Glenn - user1322491
嗨,Glenn @user1322491。 Recall 是被调用函数的简称。它是递归函数的习语。这些函数可以进一步优化。但是我认为它们对具有相当重叠的集合列表效果不错。不幸的是,事实并非如此。一些问题:1.您的集合中重叠的情况是否罕见? 2. R在大于40个列表时会卡死或者仅是内存不足?3.通过随机采样,您是否意味着在可能组合的集合上进行等概率采样? - Ferdinand.kraft
嗨,Ferdinand @Ferdinand.kraft
  1. 在任何给定的集合组合中,重叠的可能性比不重叠的要大,但我很难给出实际的比例。
  2. 当我将您的函数应用于一组包含45个变量长度和重叠的集合时,R停止响应。我甚至无法点击应用程序以强制退出(我使用Mac)。最终我不得不手动重新启动计算机。
- user1322491
我有点好奇。你能上传你的列表并发布链接吗?你确定你有一个命名的数字向量列表吗?顺便说一句,与此同时我正在尝试不同的方法 :-) - Ferdinand.kraft
嗨,费迪南德,我修改了我的原始帖子。现在它包括创建一个虚拟数据集的代码,可以创建一个命名的数字向量列表。为了使其工作,您将不得不加载开头列出的系统发育学软件包。每个数字向量都是来自单系亚支的物种名称向量,列表是所有大于10个物种且小于40个物种的单系亚支。感谢您的帮助! - user1322491
我尝试了你的列表,几个小时后出现了“节点堆栈溢出”错误 :-)。看起来你的集合之间有很少的重叠,并且可能的集合数量很大。但是如果我们稍微减少一些集合数量,就可以让它正常工作:使用L=subclades[1:30],我在不到5分钟内得到了约100万个集合。 - Ferdinand.kraft

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接