Nothing
sim.BipartiteEvol=function(nx,ny=nx,NG,dSpace=Inf,D=1,muP,muH,alphaP=0,alphaH=0,iniP=0,iniH=0,nP=1,nH=1,rP=1,rH=1,effect=1,
verbose=100,thin=1,P=NULL,H=NULL){
N=nx*ny # number of individuals
timeStep=floor(sqrt(NG/(nH+nP)))
oneByOne=F
parentsP=1:N
parentsH=1:N
pmut_act=rep(0,N)
hmut_act=rep(0,N)
changeP=c()
changeH=c()
kernel=sapply(1:nx,function(x){sapply(1:ny,function(y){exp(-(1/2)*((min(abs(1-x),abs(nx+1-x))^2+min(abs(1-y),abs(ny+1-y))^2)/(dSpace^2)))})})
Ys=rep(1:ny,each=nx,times=1)
Xs=rep(1:nx,times=ny)
# choice of the fitness function
f=function(X,x,r,alpha,pos){
.Call("fitnessFunction", X=as.numeric(t(X)), x=as.numeric(x),r=as.numeric(r),alpha=as.numeric(alpha), Ncol=as.integer(N), D=as.integer(D),
dSpace=as.numeric(dSpace), Xs=as.numeric(Xs),Ys=as.numeric(Ys),I=as.integer(pos),PACKAGE = "RPANDA")
}
# trait initial values
if(is.null(P)) P=matrix(iniP,nrow=D,ncol=N)
if(is.null(H)) H=matrix(iniH,nrow=D,ncol=N)
# creation of lists to help with the execution time
Phist=list(a=lapply(1:D,function(i){Matrix::Matrix(0,nrow=1,ncol=N,sparse=T)}))
Hhist=list(a=lapply(1:D,function(i){Matrix::Matrix(0,nrow=1,ncol=N,sparse=T)}))
for(j in 1:D){
Phist[[1]][[j]][1,]=P[j,]
Hhist[[1]][[j]][1,]=H[j,]
}
Pmut=list(a=Matrix::Matrix(0,nrow=1,ncol=N,sparse=T))
Hmut=Pmut
Pgen=Pmut
Hgen=Pmut
# fitness=sapply(1:N,function(i){p=P[,i];h=H[,i];c(f(p,h,alphaP,rP),f(p,h,alphaH,rH))}) # vectors of actual fitness
listInd=1
time=0
t=0 # number of time steps executed in the loop
# beginning of the main loop
for(u in 1:max(1,ceiling(NG/(timeStep*thin)))){
NTime=ceiling((min(u*timeStep*thin,NG)-time)/thin) # number of time steps in one loop run
t=0
# creation of the matrices in which we record what happens
phist=lapply(1:D,function(e){Matrix::Matrix(0,nrow=NTime,ncol=N,sparse=T)}) # trait of P when there is a birth/death
hhist=phist # trait of H when there is a birth/death
pmut=phist[[1]] # number of mutations for P
hmut=phist[[1]] # number of mutations for H
pgen=phist[[1]] # parents for P
hgen=phist[[1]] # parents for H
# begining of the internal loop _ execution between two timeStep
while(time<min(u*timeStep*thin,NG)){
# initialisation
if((t) %% thin ==0){
parentsP=1:N
parentsH=1:N
pmut_act=rep(0,N)
hmut_act=rep(0,N)
changeP=c()
changeH=c()
}
t=t+1
time=time+1
if(time %% verbose == 0 ) {cat("\r",time,"\r")}
for(k in 1:nH){
deadH=sample(1:N,1) # which H individual are dying
newH=0 # which H individual giving birth
# space=t(kernel[c(Y1,Y2),c(X1,X2)])
p=P[,deadH]
prob=f(H,p,rH,alphaH,deadH)
if(any(prob==Inf)|any(is.na(prob))){
is.inf=(prob==Inf|is.na(prob))
prob[!is.inf]=0.0000001
prob[is.inf]=1}
if(all(prob==0)) {
y=floor((deadH-1)/nx)+1
x=deadH-nx*(y-1)
X1=x:1
X2=2:(nx+1-x)
if(x==nx) X2=c()
Y1=y:1
Y2=2:(ny+1-y)
if(y==ny) Y2=c()
prob=t(kernel[c(Y1,Y2),c(X1,X2)])}
prob[prob==0]=min(prob[prob>0])/10
newH=sample(1:N,1,prob = prob)
mutH=(runif(1) < muH) # mutation ?
H[,deadH]=H[,newH] # new H trait values
changeH=unique(c(changeH,deadH,newH)) # which H changed since the last recording
parentsH[deadH]=parentsH[newH] # parents of new Hs
hmut_act[deadH]=hmut_act[parentsH[newH]]+mutH # number of mutations for H since last recording
# modification of the trait by mutation
# for H
if(sum(mutH)>0){
if(oneByOne){
where=sample(D,sum(mutH),replace = T)
H[D*(deadH[mutH]-1)+where]=H[D*(deadH[mutH]-1)+where]+rnorm(sum(mutH),mean=0,sd=effect)
}else{
H[,deadH[mutH]]=H[,deadH[mutH]]+rnorm(D*sum(mutH),mean = 0,sd=effect)
}
}
}
for(k in 1:nP){
deadP=sample(1:N,1) # which H individual are dying
newP=0 # which H individual giving birth
# space=t(kernel[c(Y1,Y2),c(X1,X2)])
h=H[,deadP]
prob=f(P,h,rP,alphaP,deadP)
prob[prob<0]
if(any(prob==Inf)|any(is.na(prob))){
is.inf=(prob==Inf|is.na(prob))
prob[!is.inf]=0.0000001
prob[is.inf]=1}
if(all(prob==0)) {
y=floor((deadP-1)/nx)+1
x=deadP-nx*(y-1)
X1=x:1
X2=2:(nx+1-x)
if(x==nx) X2=c()
Y1=y:1
Y2=2:(ny+1-y)
if(y==ny) Y2=c()
prob=t(kernel[c(Y1,Y2),c(X1,X2)])
}
prob[prob==0]=min(prob[prob>0])/10
newP=sample(1:N,1,prob = prob)
mutP=(runif(1) < muP) # mutation ?
P[,deadP]=P[,newP] # new H trait values
changeP=unique(c(changeP,deadP,newP)) # which H changed since the last recording
parentsP[deadP]=parentsP[newP] # parents of new Hs
pmut_act[deadP]=pmut_act[parentsP[newP]]+mutP # number of mutations for H since last recording
# modification of the trait by mutation
# for H
if(sum(mutP)>0){
if(oneByOne){
where=sample(D,sum(mutP),replace = T)
P[D*(deadP[mutP]-1)+where]=P[D*(deadP[mutP]-1)+where]+rnorm(sum(mutP),mean=0,sd=effect)
}else{
P[,deadP[mutP]]=P[,deadP[mutP]]+rnorm(D*sum(mutP),mean = 0,sd=effect)
}
}
}
# recording
if((t) %% thin ==0){
for(i in 1:D){
phist[[i]][t/thin,changeP]=P[i,changeP]
hhist[[i]][t/thin,changeH]=H[i,changeH]
}
pgen[t/thin,changeP]=parentsP[changeP]
hgen[t/thin,changeH]=parentsH[changeH]
pmut[t/thin,]=pmut_act
hmut[t/thin,]=hmut_act
}
}
Phist[[listInd]]=phist
Hhist[[listInd]]=hhist
Pmut[[listInd]]=pmut
Hmut[[listInd]]=hmut
Pgen[[listInd]]=pgen
Hgen[[listInd]]=hgen
listInd=listInd+1
}
return(list(Pgenealogy=Pgen,Hgenealogy=Hgen,xP=Phist,xH=Hhist,P=P,H=H,Pmut=Pmut,Hmut=Hmut,iniP=iniP,iniH=iniH,thin.factor=thin))
}
######## Construct the genealogy ###################################
# out an object created by the simulation function
# treeP, treeH former genealogy the tree have to be graphted on
make_gen.BipartiteEvol=function(out, treeP=NULL, treeH=NULL, verbose=T){
N=ncol(out$Pmut[[1]]) # number of individuals
NG=sum(sapply(1:length(out$Pmut),function(i){nrow(out$Pmut[[i]])})) # total time steps number
D=length(out$xP[[1]]) # dimention of trait space
time.step=nrow(out$Pmut)
# auxiliary function creating a genalogy for one type
aux=function(P.trait,Gen,X,Mut,ini,thin,tree){
edgeP=matrix(ncol=2) # the edges of the genealogy tree
xP=lapply(1:D,function(i){c()}) # trait values at the beginning of the edges
currentP=1:N # individuals at current time that have descendants in the present
node=N+1 # current node name (1:N are the tips)
nodes=1:N # current nodes
nodes_age=rep(NG+1,N) # current nodes age
edge.length=c() # length of the edges of the genealogy
nMut=rep(0,N) # number of mutation between a node and its parent node
t=NG # actual time
listInd=length(Mut)+1 # indice of the current matrix in the list
time=NG # time at the first row of the current matrix
newInd=NULL # who was born at time t
# begining of main loop
while (t>0 & length(currentP)>1){
# ie we are not yet at the root and there is more than one individual with descendant in the present
listInd=listInd-1
time.step=nrow(Mut[[listInd]])
time=time-time.step
if(time<0){
time.step=time.step+time
time=0
}
x=lapply(1:D,function(i){X[[listInd]][[i]]}) # trait values
gen=Gen[[listInd]] # parents
mut=Mut[[listInd]] # mutation number
# internal loop
while(t>time){
# ie we are not at the first row of the matrix yet
prev.t=t # ???
if(verbose) cat(paste("\r",length(currentP),t,"\r"))
if(is.null(newInd)){
newInd=which(gen[t-time,currentP]>0)
}
# if there was a birth at time t
if (sum(newInd)>0){
newCurrentP=currentP
newNodes=nodes
fathers=gen[t-time,currentP[newInd]] # parents of the new individuals
for(i in 1:length(newInd)){
father=fathers[i]
newCurrentP[newInd[i]]=father
nMut[nodes[newInd[i]]]=nMut[nodes[newInd[i]]]+(mut[t-time,currentP])[newInd[i]]
if(sum(fathers[1:i]==father)==1){
# ie the individual do not have the same parent than another born at this time already treated
if(!(father %in% currentP[newInd]) & father %in% currentP){
# ie the parent is in currentP (else there is a pb !) and did not die at time t
# we add a new node
newNodes[newInd[i]]=node
newNodes[currentP==father]=node
nodes_age=c(nodes_age,t)
nMut=c(nMut,0)
# and the corresponding offspring edges
# edge 1
edgeP=rbind(edgeP,c(node,nodes[newInd[i]]))
for(j in 1:D){
xP[[j]]=c(xP[[j]],(x[[j]][t-time,currentP])[newInd[i]])
}
edge.length=c(edge.length,nodes_age[nodes[newInd[i]]]-t)
# edge 2
edgeP=rbind(edgeP,c(node,nodes[currentP==father]))
for(j in 1:D){
xP[[j]]=c(xP[[j]],x[[j]][t-time,father])
}
edge.length=c(edge.length,nodes_age[nodes[currentP==father]]-t)
node=node+1
}else{
# ie the parent died at time t
if(sum(fathers==father)>1){
# ie there is another current individual born from the same partent at time t
# we add a new node
newNodes[newInd[i]]=node
nodes_age=c(nodes_age,t)
nMut=c(nMut,0)
# and the first of its offspring edges
edgeP=rbind(edgeP,c(node,nodes[newInd[i]]))
for(j in 1:D){
xP[[j]]=c(xP[[j]],(x[[j]][t-time,currentP])[newInd[i]])
}
edge.length=c(edge.length,nodes_age[nodes[newInd[i]]]-t)
node=node+1
}
# else we don't do anything, the current individual identity is only changed in current P (already done)
}
}else{
# ie there is an individual born at time t (and already seen) with the same parent
if(!(father %in% currentP[newInd]) & father %in% currentP){
# ie the parent did not die at time t and is in currentP (always, else there is a pb)
# we select the parent node
Node=newNodes[currentP==father]
}else{
# ie the parent died at time t
# we select the parent node
Node=(newNodes[newCurrentP==father])[1]
}
#and add the corresponding edge
edgeP=rbind(edgeP,c(Node,nodes[newInd[i]]))
for(j in 1:D){
xP[[j]]=c(xP[[j]],(x[[j]][t-time,currentP])[newInd[i]])
}
newNodes[newInd[i]]=Node
edge.length=c(edge.length,nodes_age[nodes[newInd[i]]]-t)
}
}
currentP=unique(newCurrentP)
nodes=unique(newNodes)
}
# select next time
if(t==(time+1)){
# ie the matrix has no line left
t=t-1
newInd=NULL
}else{
newInd=which(gen[t-time-1,currentP]>0)
if(sum(newInd)>0 ){
# ie there is a change in the next time step
t=t-1
}else{
# we look for the next change
newInd=NULL
if(t>1 & length(currentP)>1){
mat=Matrix::which(gen[1:(t-time-1),currentP]>0,arr.ind = T)
if(length(mat)==0) {
t=time
}else{
t=max(mat[,1])+time
ind=mat[mat[,1]==(t-time),2]
}
}else{
t=time
}}}
} # end internal loop
} # end main loop
edgeP=edgeP[-1,] # edgeP was initialized with a row of NA, we suppress it
if(length(currentP)>1){
# ie all the individual did not coalesce
if(is.null(tree)){
for(i in 1:length(currentP)){
edgeP=rbind(edgeP,c(node,nodes[i]))
for(j in 1:D){xP[[j]]=c(xP[[j]],ini)}
edge.length=c(edge.length,nodes_age[nodes[i]])
}
}else{
# we grapht the new tree on the former tree
M=max(max(edgeP),N)
newmut=rep(0,max(tree$edge))
for(j in 1:nrow(tree$edge)){ newmut[tree$edge[j,2]]=tree$nMut[j]}
nMut=c(nMut,newmut)
tree$edge=tree$edge+M
edgeP=rbind(edgeP,tree$edge)
edge.length=c(edge.length,tree$edge.length/thin)
for(j in 1:D) xP[[j]]=c(xP[[j]],tree$x[[j]])
for (i in 1:length(currentP)){
ind=which(tree$tip.label==currentP[i])+M
ind.edge=which(edgeP[,2]==ind)
nMut[nodes[i]]=nMut[nodes[i]]+nMut[edgeP[ind.edge,2]]
edgeP[ind.edge,2]=nodes[i]
edge.length[ind.edge]=edge.length[ind.edge]+nodes_age[nodes[i]]-1
}
}
}
# create a phylo object
if(is.null(tree) | length(currentP)==1){
treeP=list(Nnode=max(edgeP)-N,edge=edgeP,tip.label=1:N,edge.length=(edge.length)*thin)
}else{
treeP=list(Nnode=length(unique(as.vector(edgeP)))-(2*N-length(currentP)),edge=edgeP,tip.label=c(1:(max(edgeP))),edge.length=(edge.length)*thin)
}
class(treeP)="phylo"
X=rep(0,nrow(edgeP))
nMut=nMut[treeP$edge[,2]]
traits=xP
traits[[D+1]]=nMut
if(verbose) cat(".")
# make the tree well conformed and supress extinct lineages
treeP=rigth.order.BE(treeP,traits)
if(!is.null(tree) & (2*N-length(currentP))>N){
extinct=treeP$tree$tip.label[treeP$tree$tip.label>N]
treeP=prune.with.traits(treeP$tree,extinct,treeP$trait,extensiveTrait = D+1)}
if(verbose) cat(".")
# add the traits in the phylo object
xP=treeP$trait[1:D]
nMut=treeP$trait[[D+1]]
x.tip=matrix(0,ncol=N,nrow=D)
for(j in 1:D){
x.tip[j,]=P.trait[j,treeP$tree$tip.label]
}
if(verbose) cat(".")
treeP=treeP$tree
treeP$x=xP
treeP$ini=max(node.depth.edgelength(treeP))
treeP$nMut=nMut
treeP$x.tip=x.tip
return(treeP)
}
if(verbose) print("P")
P=aux(out$P,out$Pgenealogy,out$xP,out$Pmut,out$iniP,out$thin.factor,treeP)
if(verbose) print("H")
H=aux(out$H,out$Hgenealogy,out$xH,out$Hmut,out$iniH,out$thin.factor,treeH)
return(list(P=P,H=H))
}
######## Create a species tree from the genealogy #################
# genealogy an object created with make.gen
# threshold number of mutation to be in different species
# distanceP, distanceH distance (ie nb of mutations) matrix between the individual
# if NULL it is computed within the function
define_species.BipartiteEvol=function(genealogy,threshold=1,distanceH=NULL,distanceP=NULL, verbose=T, monophyly=TRUE, seed=NULL){
if (monophyly==TRUE){
D=length(genealogy$P$x) # dimension of trait space
# auxiliary function to define the species for one type
aux = function(gen,distance){
N=length(gen$tip.label) # number of individuals
# first we define the genetic types of the individuals
if(threshold==1){
type=rep(1,N)
new.type=2
for(i in 1:nrow(gen$edge)){
if(gen$nMut[i]>0){
if(gen$edge[i,2]<=N){
type[gen$tip.label[gen$edge[i,2]]]=new.type
new.type=new.type+1
}else{
type[as.integer(extract.clade(gen,gen$edge[i,2])$tip.label)]=new.type}
new.type=new.type+1
}
}
}else{
if(is.null(distance)){
# we compute the distance matrix if it is not given as an argument
distance=matrix(0,nrow=N,ncol=N)
for(i in 1:nrow(gen$edge)){
if(gen$nMut[i]>0){
#add the number of mutation between the individuals separated by this edge
if (verbose) cat("\r",i,"\r")
if(gen$edge[i,2]<=N){
tip=as.integer(gen$tip.label[gen$edge[i,2]])
distance[tip,(1:N)[-tip]]=distance[tip,(1:N)[-tip]]+gen$nMut[i]
distance[(1:N)[-tip],tip]=distance[(1:N)[-tip],tip]+gen$nMut[i]
}else{
tip=as.integer(extract.clade(gen,gen$edge[i,2])$tip.label)
distance[tip,(1:N)[-tip]]=distance[tip,(1:N)[-tip]]+gen$nMut[i]
distance[(1:N)[-tip],tip]=distance[(1:N)[-tip],tip]+gen$nMut[i]
}
}
}
distance=distance[as.integer(gen$tip.label),as.integer(gen$tip.label)]
}
# we use the distance matrix to define the genetic type of individuals
phylo=gen
phylo$edge.length=phylo$nMut
type=1:N
for(i in 1:N){
ind=((1:i)[distance[i,1:i]<threshold])[1]
type[phylo$tip.label[i]]=type[phylo$tip.label[ind]]
if (verbose) cat("\r",i,"\r")
}
}
# species are then defined as monophyletic types
species=rep(1,N) # the species of each tip
newSpecies=2
nextNode=N+1
# test of the tips descendig from a particular node _ main loop
while(length(nextNode)>0){
i=nextNode[1] # candidate node
nextNode=nextNode[-1] # update candidate nodes vector
a=extract.clade(gen,i) # subtree with root i
ty=as.integer(gen$tip.label)
ty=ty[which(!(ty %in% as.integer(a$tip.label)))]
ty=type[ty] # type of tips not in a
if(length(intersect(type[as.integer(a$tip.label)],ty))==0){
# seeing how nextNode is built, this should always be the case... Test needed
n=length(a$tip.label) # number of tips in the subtree
offspring=gen$edge[gen$edge[,1]==i,2] # offspring nodes of i in the main tree
offspringSubTree=a$edge[a$edge[,1]==(n+1),2] # offspring node of i in the subtree
if(length(offspringSubTree)>1){
for(j in 1:(length(offspringSubTree))){
if(offspringSubTree[j]<=n){
# ie offspringSubTree[j] is a tip
ty=as.integer(a$tip.label)
ty=ty[which(!(ty %in% as.integer(a$tip.label[offspringSubTree[j]])))]
ty=type[ty] # type of tips in a that are not offspring[j]
if(!(type[as.integer(a$tip.label[offspringSubTree[j]])] %in% ty)){
# offspringSubTree[j] is a species
species[as.integer(a$tip.label[offspringSubTree[j]])]=newSpecies
newSpecies=newSpecies+1
}
}else{
# ie offspringSubTree[j] is an internal node
a1=extract.clade(a,offspringSubTree[j])
ty=as.integer(a$tip.label)
ty=ty[which(!(ty %in% as.integer(a1$tip.label)))]
ty=type[ty] # type of tips in a but not in a1
if(length(intersect(type[as.integer(a1$tip.label)],ty))==0){
# the tips in a1 are in a different species than the other tips in a
nextNode=c(nextNode,offspring[j]) # new candidate node
species[as.integer(a1$tip.label)]=newSpecies
newSpecies=newSpecies+1
}
}
}
}else{print("error : one node has only one offspring... ???")}
}else{print(paste("error : this node should not be in next node :",i))}
} # end of main loop
# rename the species so that there is no gap in the names
M=rep(0,max(species))
name=1
for(i in 1:length(species)){
if(M[species[i]]==0){
M[species[i]]=name
name=name+1
}
species[i]=M[species[i]]
}
return(species)
}
# auxiliary function to build the species tree from the genealogy and the species of each tip
make.phylo=function(gen,spec){
N=length(gen$tip.label) # number of individuals
abundance=sapply(unique(spec),function(i){sum(spec==i)}) # abundance of each species
mean.trait=sapply(unique(spec),function(i){sapply(1:D,function(e){mean(gen$x.tip[e,spec[gen$tip.label]==i])})}) # mean trait of each species
if(D == 1){mean.trait = matrix(mean.trait,nrow = D)}
if(max(spec)==1){return(list(abundance=abundance,mean.trait=mean.trait))}
keep.tip=c() # which tips will be in the species tree
keep.tip=sapply(unique(spec),function(i){vect=(1:N)[spec[gen$tip.label]==i];
if(length(vect)==1){return(vect)};
sample(vect,1)})
tree=prune.with.traits(gen,gen$tip.label[-keep.tip],gen$x) # suppress the other tips
tree$abundance=abundance[spec[as.integer(tree$tree$tip.label)]]
tree$mean.trait=matrix(mean.trait[,spec[as.integer(tree$tree$tip.label)]], nrow = D, byrow = F)
tree$tree$tip.label=spec[as.integer(tree$tree$tip.label)]
return(tree)
}
P=aux(genealogy$P,distanceP)
H=aux(genealogy$H,distanceH)
Pphylo=make.phylo(genealogy$P,P)
Hphylo=make.phylo(genealogy$H,H)
} else { # Monophyly == False
D=length(genealogy$P$x) # dimention of trait space
# auxiliary function to define the species for one type
aux = function(gen,distance){
N=length(gen$tip.label) # number of individuals
# first we define the genetic types of the individuals
if(threshold!=1){
print("Threshold set to 1 when monophyly==FALSE")
threshold <- 1}
if(threshold==1){
type=rep(1,N)
new.type=2
for(i in 1:nrow(gen$edge)){
if(gen$nMut[i]>0){
if(gen$edge[i,2]<=N){
type[gen$tip.label[gen$edge[i,2]]]=new.type
new.type=new.type+1
}else{
type[as.integer(extract.clade(gen,gen$edge[i,2])$tip.label)]=new.type}
new.type=new.type+1
}
}
}
# rename the species (=type) so that there is no gap in the names
M=rep(0,max(type))
name=1
for(i in 1:length(type)){
if(M[type[i]]==0){
M[type[i]]=name
name=name+1
}
type[i]=M[type[i]]
}
return(type)
}
# auxiliary function to build the species tree from the genealogy and the species of each tip
make.phylo=function(gen,spec){
N=length(gen$tip.label)
list_species <- unique(spec)
unique_type <- c()
pruned_gen <- gen
pruned_gen$tip.label <- as.character(pruned_gen$tip.label)
pruned_gen <- multi2di(pruned_gen)
if (!is.null(seed)) {
set.seed(seed)
sampled_tips <- sample(1:N)} else {sampled_tips <- 1:N}
for (i in sampled_tips){
if (spec[i] %in% unique_type){
pruned_gen <- drop.tip(pruned_gen,tip=as.character(i), trim.internal = TRUE)
}else{unique_type <- c(unique_type, spec[i])}
}
pruned_gen$tip.label <- spec[as.numeric(pruned_gen$tip.label)]
abundance=sapply(unique(spec),function(i){sum(spec==i)}) # abundance of each species
mean.trait=sapply(unique(spec),function(i){sapply(1:D,function(e){mean(gen$x.tip[e,spec[gen$tip.label]==i])})}) # mean trait of each species
if(D == 1){mean.trait = matrix(mean.trait,nrow = D)}
if(max(spec)==1){return(list(abundance=abundance,mean.trait=mean.trait))}
tree=list(tree=pruned_gen)
tree$abundance=abundance[tree$tree$tip.label]
tree$mean.trait=matrix(mean.trait[,tree$tree$tip.label], nrow = D, byrow = F)
# tree$mean.trait and tree$abundance order like tree$tree$tip.label
return(tree)
}
P=aux(genealogy$P,distanceP)
H=aux(genealogy$H,distanceH)
Pphylo=make.phylo(genealogy$P,P)
Hphylo=make.phylo(genealogy$H,H)
}
return(list(P=P,H=H,Pphylo=Pphylo,Hphylo=Hphylo))
}
######## Build the network ##################
build_network.BipartiteEvol=function(gen, spec){
P=spec$P # what Pspecies is each P individual in ?
H=spec$H # what Hspecies is each H individual in ?
N=length(P) # number of individual in each clade
nSP=max(P) # number of Pspecies
nSH=max(H)
link=Matrix::Matrix(0,nrow = nSP,ncol = nSH,sparse = T) # initialize the link matrix
for(i in 1:N){
link[P[i],H[i]]=link[P[i],H[i]]+1
}
return(link)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.