-
Notifications
You must be signed in to change notification settings - Fork 1
/
RF_FFT.R
126 lines (111 loc) · 3.57 KB
/
RF_FFT.R
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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
library("ape")
library("phangorn")
library("readtext")
#Beta function
Beta=function(m){
if(m<0){return(1)}
return(dfactorial(2*m+1))
}
# the function for computing the number of internal edges
internaledges <- function(tree,ntip){
intedges=array(0,c(1,ntip-1))
edges=tree$edge
for (i in (2*ntip-1):(ntip+1)) {
children=which(edges[,1]==i)
child1=edges[children[1],2]
child2=edges[children[2],2]
if((child1 <= ntip)&(child2 <= ntip)){intedges[i-ntip]=0}
else if((child1<= ntip) & (child2 > ntip)) {intedges[i-ntip]=intedges[child2-ntip]+1}
else if((child2<= ntip )& (child1 > ntip)){intedges[i-ntip]=intedges[child1-ntip]+1}
else {intedges[i-ntip]=intedges[child2-ntip]+intedges[child1-ntip]+2}
}
return(intedges)
}
internalchildren <- function(tree,v,ntip){
edges=tree$edge
children=which(edges[,1]==v)
child1=edges[children[1],2]
child2=edges[children[2],2]
if((child1 > ntip) & (child2 > ntip)){result=c(2,child1,child2)}
else if((child1 > ntip) & (child2 <= ntip)){result=c(1,child1)}
else if((child2 > ntip) & (child1 <= ntip)){result=c(1,child2)}
else {result=0}
return(result)
}
RF_Convolve=function(tree,n){
ntip=n-1
N=tree$Nnode
R=rep(list(matrix(0,(ntip-1),(ntip-1))),N)
edges=internaledges(tree,ntip)
B=c()
for (k in 0:(n-2)) {
B[k+1]=Beta(k)
}
for (v in N:1) {
intchild=internalchildren(tree,v+ntip,ntip)
intedges=edges[v]
if(intchild[1]==0){
R[[v]][1,1]=1
}
else if(intchild[1]==1){
Rchild=R[[intchild[2]-ntip]]
R[[v]][1,intedges+1]=1
R[[v]][2:(ntip-1),1]=rowSums(t(t(Rchild[1:(ntip-2),])*B[1:(ntip-1)]))
R[[v]][2:(ntip-1),2:(ntip-1)]=Rchild[2:(ntip-1),1:((ntip-2))]
}
else {
Rchild1=R[[intchild[2]-ntip]]
Rchild2=R[[intchild[3]-ntip]]
R[[v]][1,intedges+1]=1
R[[v]][3,1]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))
for (s in 4:(ntip-1)) {
R[[v]][s,1]=sum(rowSums(t(t(Rchild1[1:(s-2),])*B[1:(ntip-1)]))*rowSums(t(t(Rchild2[(s-2):1,])*B[1:(ntip-1)])))
}
sum1=matrix(0,(ntip-2),(ntip-2))
sum1[1,1:(ntip-2)]=sum(t(t(Rchild1[1,])*B[1:(ntip-1)]))*Rchild2[1,1:(ntip-2)]
for (s in 3:(ntip-1)) {
temp=colSums(rowSums(t(t(Rchild1[1:(s-1),])*B[1:(ntip-1)]))*Rchild2[(s-1):1,1:(ntip-2)])
sum1[s-1,1:(ntip-2)]=temp
}
sum2=matrix(0,(ntip-2),(ntip-2))
sum2[1,1:(ntip-2)]=sum(t(t(Rchild2[1,])*B[1:(ntip-1)]))*Rchild1[1,1:(ntip-2)]
for (s in 3:(ntip-1)) {
temp=colSums(rowSums(t(t(Rchild2[1:(s-1),])*B[1:(ntip-1)]))*Rchild1[(s-1):1,1:(ntip-2)])
sum2[s-1,1:(ntip-2)]=temp
}
R1=Rchild1[1:(ntip-1),1:(ntip-3)]
R2=Rchild2[1:(ntip-1),1:(ntip-3)]
R1aug=cbind(R1,matrix(0,nrow(R1),ncol(R1)))
R2aug=cbind(R2,matrix(0,nrow(R2),ncol(R2)))
U=round(convolve(t(R1aug), rev(t(R2aug)), type="open"),0)
sum3=t(matrix(c(U,0),2*ncol(R1),2*nrow(R1))[1:ncol(R1),1:nrow(R1)])
sum3=cbind(array(0, dim=c(nrow(R1)-1,1)),sum3[2:nrow(R1),])
R[[v]][2:(ntip-1),2:(ntip-1)]=sum1+sum2+sum3
}
}
return(R)
}
#==========================================
RsT=function(R,n,s){
B=c()
for (k in 0:(n-2)) {
B[k+1]=Beta(k)
}
rst =sum(t(t(R[[1]][s+1,1:(n-2-s)])*B[1:(n-2-s)]))
return(rst)
}
#Compute the value of q_m(T)
qmT=function(R,n,m){
qmt=0
for (s in m:(n-3)) {
rst=RsT(R,n,s)
qmt=qmt+(factorial(s)/(factorial(m)*factorial(s-m)))*rst*(-1)^(s-m)
}
return(qmt)
}
polynomial=function(tree,n){
R=RF_Convolve(tree,n)
for (i in seq(0,2*(n-3),2)) {
print(qmT(R,n,n-3-(i/2)))
}
}