1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
prepareDataSankoffNew <- function(data){
contrast = attr(data, "contrast")
contrast[contrast == 0] = 1.0e+06
contrast[contrast == 1] <- 0.0
attr(data, "contrast") <- contrast
data
}
sankoffNew <- function (tree, data, cost = NULL, site = 'pscore')
{
if (class(data) != "phyDat")
stop("data must be of class phyDat")
data <- prepareDataSankoffNew(data)
weight <- attr(data, "weight")
levels <- attr(data, "levels")
l = length(levels)
if (is.null(cost)) {
cost <- matrix(1, l, l)
cost <- cost - diag(l)
}
l <- length(data)
nr <- attr(data, "nr")
# on.exit(.C("sankoff_free"))
# .C("sankoff_init", as.integer())
# for (i in 1:length(data)) storage.mode(data[[i]]) = "double"
if(class(tree)=="phylo") return(fit.sankoffNew(tree, data, cost, returnData =site))
if(class(tree)=="multiPhylo"){
if(is.null(tree$TipLabel))tree = unclass(tree)
return(sapply(tree, fit.sankoffNew, data, cost, site))
}
}
fit.sankoffNew <- function (tree, data, cost, returnData = c("pscore", "site", "data"))
{
if (is.null(attr(tree, "order")) || attr(tree, "order") ==
"cladewise")
tree <- reorder(tree, "postorder")
returnData <- match.arg(returnData)
node <- tree$edge[, 1]
edge <- tree$edge[, 2]
weight = attr(data, "weight")
nr = p = attr(data, "nr")
contr = attr(data, "contrast")
q = length(tree$tip.label)
nc = l = attr(data, "nc")
m = length(edge) + 1
dat = vector(mode = "list", length = m)
dat[1:q] = data[tree$tip.label]
node = as.integer(node - 1)
edge = as.integer(edge - 1)
nTips = as.integer(length(tree$tip))
mNodes = as.integer(max(node) + 1)
# tips = as.integer((1:length(tree$tip))-1)
res <- .Call("sankoff3B", dat, as.numeric(cost), as.integer(nr),as.integer(nc),
node, edge, mNodes, nTips, as.double(contr), as.integer(nrow(contr)), PACKAGE="phangorn")
root <- getRoot(tree)
erg <- .Call("C_rowMin", res[[root]], as.integer(nr), as.integer(nc), PACKAGE = "phangorn")
if (returnData=='site') return(erg)
pscore <- sum(weight * erg)
result = pscore
if (returnData=="data"){
result <- list(pscore = pscore, dat = res)
}
result
}
pnodesNew <- function (tree, data, cost)
{
if (is.null(attr(tree, "order")) || attr(tree, "order") ==
"cladewise")
tree <- reorder(tree, "postorder")
node <- tree$edge[, 1]
edge <- tree$edge[, 2]
nr = nrow(data[[1]])
nc = ncol(data[[1]])
node = as.integer(node - 1)
edge = as.integer(edge - 1)
.Call("pNodes", data, as.numeric(cost), as.integer(nr),as.integer(nc),
node, edge, PACKAGE="phangorn")
}
|
