PhyloGaussianBeliefProp

AbstractClusterGraphMethod

Abstract type for cluster graph construction algorithms.

BPPosDefException

Exception thrown when a belief message cannot be computed, that is, when the submatrix of the precision J, subsetted to the variables to be integrated out, is not positive definite. It has a message msg field (string), and an info field (integer) inherited from LinearAlgebra.PosDefException, to indicate the location of (one of) the eigenvalue(s) which is (are) less than/equal to 0.

Belief{T<:Real,Vlabel<:AbstractVector,P<:AbstractMatrix{T},V<:AbstractVector{T},M} <: AbstractBelief{T}

A "belief" is an exponential quadratic form, using the canonical parametrization:

C(x | J,h,g) = exp( -(1/2)x'Jx + h'x + g )

It is a normalized distribution density if

g = - (1/2) (log(2πΣ) + μ'Jμ)
  = - entropy of normalized distribution + (1/2) dim(μ) - (1/2) μ'Jμ.

• μ is the mean, of type V (stored but typically not updated)
• h = inv(Σ)μ is the potential, also of type V,
• Σ is the variance matrix (not stored)
• J = inv(Σ) is the precision matrix, of type P
• g is a scalar to get the unnormalized belief: of type MVector{1,T} to be mutable.

See MvNormalCanon{T,P,V} in Distributions.jl

Other fields are used to track which cluster or edge the belief corresponds to, and which traits of which variables are in scope:

• nodelabel of type Vlabel
• ntraits
• inscope
• type: cluster or sepset
• metadata of type M: Symbol for clusters, Tuple{Symbol,Symbol} for sepsets.

Methods for a belief b:

• nodelabels(b): vector of node labels, of type Vlabel, e.g. Int8 if nodes are labelled by their preorder index in the original phylogeny
• ntraits(b): number of traits (dimension of the random variable x above)
• inscope(b): matrix of booleans (trait i in row i and and node j in column j)
• nodedimensions(b): vector of integers, with jth value giving the dimension (number of traits in scope) of node j.
• dimension(b): total dimension of the belief, that is, total number of traits in scope. Without any missing data, that would be ntraits × length of nodelabels.

Belief(nodelabels, numtraits, inscope, belieftype, metadata, T=Float64)

Constructor to allocate memory for one cluster, and initialize objects with 0s to initialize the belief with the constant function exp(0)=1.

Bethe

Subtype of AbstractClusterGraphMethod.

Algorithm

A Bethe cluster graph (also known as factor graph) has:

• a factor-cluster {v, parents(v}} for each node-family in the network, that is, for each non-root node v (a family is a child node and all of its parents)
  • with one exception: if v's family is included in another family, then no factor-cluster is created for v.
• a variable-cluster {v} for each non-leaf node v in the network, or more specifically, for each node v that belongs in more than 1 factor.

Each variable-cluster {v} is joined to the factor-clusters that contain v, by an edge labelled with sepset {v}.

References

D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT Press, 2009. ISBN 9780262013192.

Cliquetree

Subtype of AbstractClusterGraphMethod.

Algorithm

1. moralize the network (connect partners that share a child).
2. triangulate the resulting undirected graph using greedy min-fill, see triangulate_minfill!(graph).
3. extract the maximal cliques of the resulting chordal graph.
4. calculate the edge weight between each pair of maximal cliques as the size of their intersection
5. find a maximum-weight spanning tree
6. label the retained edges (in the spanning tree) by the intersection of the two cliques they connect.

References

D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT Press, 2009. ISBN 9780262013192.

ClusterGraphBelief{B<:Belief, F<:FamilyFactor, M<:MessageResidual}

Structure to hold a vector of beliefs, with cluster beliefs coming first and sepset beliefs coming last. Fields:

• belief: vector of beliefs
• factor: vector of initial cluster beliefs after factor assignment
• nclusters: number of clusters
• cdict: dictionary to get the index of a cluster belief from its node labels
• sdict: dictionary to get the index of a sepset belief from the labels of its two incident clusters
• messageresidual: dictionary to log information about sepset messages, which can be used to track calibration or help adaptive scheduling with residual BP. See MessageResidual. The keys of messageresidual are tuples of cluster labels, similar to a sepset's metadata. Each edge in the cluster graph has 2 messages corresponding to the 2 directions in which a message can be passed, with keys: (label1, label2) and (label2, label1). The cluster receiving the message is the first label in the tuple, and the sending cluster is the second.

Assumptions:

• For a cluster belief, the cluster's nodes are stored in the belief's metadata.
• For a sepset belief, its incident clusters' nodes are in the belief's metadata.

ClusterGraphBelief(belief_vector::Vector{B})

Constructor of a ClusterGraphBelief with belief belief_vector and all other fields constructed accordingly. New memory is allocated for these other fields, e.g. for factors (with data copied from cluster beliefs) and message residuals (with data initialized to 0 but of size matching that from sepset beliefs)

To construct the input vector of beliefs, see init_beliefs_allocate and init_beliefs_assignfactors!

EvolutionaryModel{T}

Evolutionary model type, with T the element type in all parameter vector and matrices. Implemented models include the UnivariateBrownianMotion.

An object of this type must contain at least the following elements:

• μ: the mean of the trait at the root.
• v: the variance of the trait at the root. Can be zero (fixed root) or infinite.

New evolutionary models must implement the following interfaces:

params(obj::EvolutionaryModel)
params_optimize(obj::EvolutionaryModel)
params_original(obj::EvolutionaryModel, transformedparams::AbstractArray)

FamilyFactor(belief::AbstractBelief{T}) where T

Constructor to allocate memory for one family factor, with canonical parameters and metadata initialized to be a copy of those in belief. FamilyFactors metadata are supposed to be symbols, so this constructor should fail if its input is a sepset belief, whose metadata is a Tuple of Symbols.

HeterogeneousBrownianMotion{T,U,V,W} <: HeterogeneousBM{T}

Type for a heterogeneous Brownian motion model, univariate or multivariate. Along each edge, evolution follows a Brownian motion. Each edge can have its own variance rate. This model has no shifts, and no extra hybrid variance. By default, the root is fixed with prior variance 0.

T is the scalar type, U is the type for the root mean (vector of length d, where d is the trait dimension, even if univariate), V is a matrix type for the root variance, and W the matrix type of each variance rate, one per color. For a univariate BM, we may have W=T and V=Vector{T}. For a multivariate BM, we may have W=V<:Matrix{T}. This is such that each field is mutable, so we can update model parameters in place within the model object, itself immutable.

HeterogeneousShiftedBrownianMotion{T,U,V,W} <: HeterogeneousBM{T}

Type for a heterogeneous Brownian motion model like HeterogeneousBrownianMotion but with a possible shift in the mean along each edge.

JoinGraphStructuring

Subtype of AbstractClusterGraphMethod.

Fieldnames

• maxclustersize: upper limit for cluster size. This value must be at least the size of the largest node family in the input phylogenetic network, normally 3 if the network is bicombining (each hybrid node has 2 parents, never more). See nodefamilies(net).

Constructors

• JoinGraphStructuring(maxclustersize, net): checks that the input maxclustersize is valid for net

Algorithm, by Mateescu et al. (2010)

Requires:

• a user-specified maximum cluster size
• an elimination order for the nodes in the HybridNetwork, hopefully yielding a small induced-width, e.g. from a heuristic such as greedy min-fill (see triangulate_minfill!(graph)).

1. Each node in net labels a "bucket", and these buckets are ordered according to the elimination order, e.g. the highest-priority bucket is labelled by the first node in the elimination order.
2. Node families are assigned to buckets based on the highest-priority node they contain. For example, if {1,4,9} forms a node family and if node 4 is higher in the elimination order than nodes 1 or 9, then {1,4,9} gets assigned to bucket 4.
3. Node families are clusters (of nodes), and we refer to the clusters within a bucket as "minibuckets". Minibuckets within a bucket can be merged as long as the size of their union does not exceed the maximum cluster size allowed. Sometimes, this can be done in multiple ways.
4. Starting with the highest-priority bucket, we create new minibuckets by "marginalizing out" the bucket label from each existing minibucket (these are left unchanged in the process).
5. Each new minibucket is joined to its "originator" minibucket by an edge that is labeled by their intersection (i.e. the variables in the new minibucket). Each new minibucket is then reassigned to a new (and necessarily lower-priority) bucket based on its highest priority node. Merging can take place during reassignment as long as the the maximum cluster size is respected. The union of 2 minibuckets retains the edges of each of minibucket.
6. Steps 4 & 5 are carried out for each bucket in order of priority.
7. The resulting minibuckets in each bucket are then joined in a chain (there is some degree of freedom for how this can be done), where each edge is labelled by the bucket label.

References

R. Mateescu, K. Kask, V.Gogate, and R. Dechter. Join-graph propagation algorithms. Journal of Artificial Intelligence Research, 37:279-328, 2010 doi: 10.1613/jair.2842.

LTRIP{T<:Integer}

Subtype of AbstractClusterGraphMethod. A HybridNetwork and a valid LTRIP are passed to clustergraph! to construct a cluster graph from the user-provided clusters based on the Layered Trees Running Intersection Property algorithm of Streicher & du Preez (2017).

Fieldnames

• clusters: vector of clusters, required to be family-preserving with respect to some HybridNetwork – see isfamilypreserving. Within each cluster, nodes (identified by their preorder index in the network) are required to be sorted in decreasing order (for postorder)

Constructors

• LTRIP(net): uses nodefamilies(net) as input clusters, which are guaranteed to be family-preserving
• LTRIP(clusters, net): checks if that clusters provided are family-preserving, then sorts each cluster in decreasing order (modifying them in place!) before creating the LTRIP object.

They assume, with no check, that net already has a preordering.

Algorithm

1. An initial graph G is considered, in which each input cluster is a node. An edge (C1,C2) is added if clusters C1 and C2 share at least 1 node (in net). The weight of edge (C1,C2) is defined as the size of the intersection C1 ∩ C2.
2. For each node n in net, the subgraph of G induced by the clusters containing n, G_n, has its weights adjusted as follows:
   • the edges of maximum weight (within G_n) are identified, then
   • the weight of each edge is increased by the number of max-weight edges that either of its endpoints adjacent to.
   Then, LTRIP finds a maximum-weight spanning tree of G_n. The edges of this tree are all labelled with {n} (or its label or preorder index).
3. The spanning trees for each node are layered on one another to form a cluster graph. In other words, an edge (C1,C2) is added if it is present is any spanning tree. If so, its sepset is the union of its labels across the different spanning trees.

References

S. Streicher and J. du Preez. Graph Coloring: Comparing Cluster Graphs to Factor Graphs. In Proceedings of the ACM Multimedia 2017 Workshop on South African Academic Participation, pages 35-42, 2017. doi: 10.1145/3132711.3132717.

MessageResidual{T<:Real, P<:AbstractMatrix{T}, V<:AbstractVector{T}} <: AbstractResidual{T}

Structure to store the most recent computation history of a message, in the form of the ratio: sent_message / current_sepset_belief, when a message is sent from one cluster to another along a given sepset. At calibration, this ratio is 1. For Gaussian beliefs, this ratio is an exponential quadratic form, stored using its canonical parametrization, excluding the constant.

Fields:

• Δh: canonical parameter vector of the message residual
• ΔJ: canonical parameter matrix of the message residual
• kldiv: kl divergence between the message that was last sent and the sepset belief before the last update
• iscalibrated_resid: true if the last message and prior sepset belief were approximately equal, false otherwise. see iscalibrated_residnorm!
• iscalibrated_kl: same, but in terms of the KL divergence, see iscalibrated_kl!.

MessageResidual(J::AbstractMatrix{T}, h::AbstractVector{T})

Constructor to allocate memory for a MessageResidual with canonical parameters (ΔJ, Δh) of the same dimension and type as J and h, initialized to zeros. kldiv is initalized to [-1.0] and the flags iscalibrated_{resid,kl} are initialized to false if the message is of positive dimension. If the message is empty (ΔJ and Δh of dimension 0) then the message is initialized as being calibrated: kldiv is set to 0 and iscalibrated flags set to true.

(ΔJ, Δh) of zero suggest calibration, but the flags iscalibrated_{resid,kl} being false indicate otherwise.

MvDiagBrownianMotion{T,V} <: HomogeneousBrownianMotion{T}

The multivariate Brownian motion with diagonal variance rate matrix, that is, traits with independent evolution. It is homogeneous across the phylogeny. R is the variance rate (stored as a vector of type V), μ is the prior mean at the root and v the prior variance at the root, 0 by default (and both also of type V)

MvFullBrownianMotion{T,P1,V,P2} <: HomogeneousBrownianMotion{T}

The full multivariate Brownian motion. It is homogeneous across the phylogeny. R is the variance rate (of matrix type P1), μ is the prior mean at the root (of vector type V) and v the prior variance at the root, 0 by default (of matrix type P2).

MvFullBrownianMotion(R::AbstractMatrix, μ, v=nothing)
MvFullBrownianMotion(U::AbstractVector, μ, v=nothing)

Constructor for a full multivariate Brownian motion (homogeneous) with variance rate matrix V and prior mean vector μ at the root. If not provided, the prior variance matrix at the root v is set to 0.

If a vector U is provided, it is used as the Upper cholesky factor of R=U'U, vectorized, so U should be of length p(p+1)/2 where p is the number of traits (also the length of μ).

PaintedParameter{T}

Type with 2 fields:

• parameter: vector whose elements are of type T
• color: DefaultDict dictionary mapping integers to integers, with a default value of 1.

ncolors(pp) returns the number of parameters, that is, the length of pp.parameter This type is meant to store several values for a given evolutionary parameter (say, Brownian motion variance rate), each one being used on some edges or nodes of a phylogenetic network. The default parameter value is the first one.

For an edge number i, color j=pp.color[i] indexes its parameter value, that is, evolution along edge number i should use pp.parameter[j]. This parameter value is obtained with getparameter(pp, j) for the parameter value of color j, or getparameter(pp, edge) for the parameter value of edge.

UnivariateBrownianMotion{T} <: HomogeneousBrownianMotion{T}

The univariate Brownian motion, homogeneous across the phylogeny, that is, with the same variance rate σ2 across all edges. μ is the prior mean at the root. v the prior variance at the root, 0 by default.

UnivariateOrnsteinUhlenbeck{T} <: HomogeneousOrnsteinUhlenbeck{T}

The univariate Ornstein-Uhlenbeck model. It is homogeneous, that is, has the same parameters across all edges in the phylogeny. σ2 is the variance rate, α the selection strength, θ the optimal value, μ the prior mean at the root, and v the prior variance at the root, 0 by default.

absorbevidence!(h,J,g, dataindex, datavalues)

Absorb evidence, at indices dataindex and using datavalues. Warnings:

• a subset of h is modified in place
• traits are assumed to come in the same order in dataindex as in datavalues.

absorbleaf!(h,J,g, rowindex, columntable)

Absorb evidence from a leaf, given in col[rowindex] of each column in the table, then marginalizes out any variable for a missing trait at that leaf. See absorbevidence! and marginalizebelief. Warning: The leaf traits are assumed to correspond to the first variables in h (and J), as is output by factor_treeedge.

addtreenode_belowdegeneratehybrid!(net::HybridNetwork)

If a degenerate hybrid node h1 has 1 child edge of length t>0 to a hybrid child h2: break the edge by adding a tree node at distance t from h1 and 0 from h2. That way, h1 may be removed from scope. This is done iteratively, as h2 may become degenerate after this operation. See shrinkdegenerate_treeedges to remove degenerate internal tree nodes, and hasdegenerate to check if net still has degenerate nodes.

assign!(bucket, new_minibucket, max_minibucket_size)

Merge new_minibucket with one of the minibuckets contained in bucket (in order of decreasing size) subject to the constraint that the resulting minibucket does not exceed max_minibucket_size. If this is not possible, then new_minibucket is added to bucket as a new minibucket. The bucket should be represented as a dictionary: minibucket_size => vector of minibuckets of that size, where each minibucket is itself represented as a vector of indices.

Output:

• (resulting_minibucket, minibucket_merged_into) if a successful merge is found
• (new_minibucket, []) otherwise

average_energy!(ref::Belief, target::AbstractBelief)
average_energy!(ref::Belief, Jₜ, hₜ, gₜ)
average_energy(Jᵣ::Union{LA.Cholesky,PDMat}, μᵣ, Jₜ, hₜ, gₜ)

Average energy (i.e. negative expected log) of a target canonical form with parameters (Jₜ, hₜ, gₜ) with respect to a normalized non-degenerate reference canonical form ref with parameters (Jᵣ, hᵣ). The reference distribution is normalized, so specifying gᵣ is unnecessary. When the target canonical form is also normalized and non-degenerate, this is equal to their cross-entropy:

H(fᵣ, fₜ) = - Eᵣ(log fₜ) = - ∫ fᵣ log fₜ .

ref is assumed to be non-degenerate, that is, Jᵣ should be positive definite.

average_energy! modifies the reference belief by updating ref.μ to J⁻¹h. It calls average_energy after a cholesky decomposition of ref.J, stored in Jᵣ: see getcholesky_μ!.

Calculation:

ref: f(x) = C(x | Jᵣ, hᵣ, _) is the density of 𝒩(μ=Jᵣ⁻¹hᵣ, Σ=Jᵣ⁻¹) target: C(x | Jₜ, hₜ, gₜ) = exp( - (1/2)x'Jₜx + hₜ'x + gₜ )

E[-log C(X | Jₜ, hₜ, gₜ)] where X ∼ C(Jᵣ, hᵣ, _)
= 0.5 (μᵣ'Jₜ μᵣ + tr(Jᵣ⁻¹Jₜ)) - hₜ'μᵣ - gₜ

With empty vectors and matrices (J's of dimension 0×0 and h's of length 0), the result is simply: - gₜ.

average_energy!(ref::Belief, target::AbstractBelief)
average_energy!(ref::Belief, Jₜ, hₜ, gₜ)
average_energy(Jᵣ::Union{LA.Cholesky,PDMat}, μᵣ, Jₜ, hₜ, gₜ)

Average energy (i.e. negative expected log) of a target canonical form with parameters (Jₜ, hₜ, gₜ) with respect to a normalized non-degenerate reference canonical form ref with parameters (Jᵣ, hᵣ). The reference distribution is normalized, so specifying gᵣ is unnecessary. When the target canonical form is also normalized and non-degenerate, this is equal to their cross-entropy:

H(fᵣ, fₜ) = - Eᵣ(log fₜ) = - ∫ fᵣ log fₜ .

ref is assumed to be non-degenerate, that is, Jᵣ should be positive definite.

average_energy! modifies the reference belief by updating ref.μ to J⁻¹h. It calls average_energy after a cholesky decomposition of ref.J, stored in Jᵣ: see getcholesky_μ!.

Calculation:

ref: f(x) = C(x | Jᵣ, hᵣ, _) is the density of 𝒩(μ=Jᵣ⁻¹hᵣ, Σ=Jᵣ⁻¹) target: C(x | Jₜ, hₜ, gₜ) = exp( - (1/2)x'Jₜx + hₜ'x + gₜ )

E[-log C(X | Jₜ, hₜ, gₜ)] where X ∼ C(Jᵣ, hᵣ, _)
= 0.5 (μᵣ'Jₜ μᵣ + tr(Jᵣ⁻¹Jₜ)) - hₜ'μᵣ - gₜ

With empty vectors and matrices (J's of dimension 0×0 and h's of length 0), the result is simply: - gₜ.

betheclustergraph(net)

See Bethe

branch_actualization(obj::EvolutionaryModel, edge::PN.Edge) 
branch_displacement(obj::EvolutionaryModel,  edge::PN.Edge)
branch_precision(obj::EvolutionaryModel,     edge::PN.Edge)
branch_variance(obj::EvolutionaryModel,      edge::PN.Edge)
branch_logdet(obj::EvolutionaryModel,        edge::PN.Edge, precision::AbstractMatrix)
branch_transition_qωjg(obj::EvolutionaryModel,    edge)
branch_transition_qωv!(q, obj::EvolutionaryModel, edge)

Under the most general linear Gaussian model, X₀ given X₁ is Gaussian with conditional mean q X₁ + ω and conditional variance Σ independent of X₁. branch_actualization, branch_displacement and branch_variance return, respectively, q, ω and Σ. branch_precision and branch_variance should return a matrix of symmetric type. branch_variance defaults to the inverse of branch_precision. branch_logdet defaults to g = -0.5*log(|2πΣ|), the log normalizing constant of the Gaussian density in the traditional form. branch_transition_* return or modify in place the corresponding transition matrices.

Under a Brownian motion, we have q=I, ω=0, and conditional variance t*R where R is the model's variance rate and t the branch length.

branch_actualization(obj::EvolutionaryModel, edge::PN.Edge) 
branch_displacement(obj::EvolutionaryModel,  edge::PN.Edge)
branch_precision(obj::EvolutionaryModel,     edge::PN.Edge)
branch_variance(obj::EvolutionaryModel,      edge::PN.Edge)
branch_logdet(obj::EvolutionaryModel,        edge::PN.Edge, precision::AbstractMatrix)
branch_transition_qωjg(obj::EvolutionaryModel,    edge)
branch_transition_qωv!(q, obj::EvolutionaryModel, edge)

Under the most general linear Gaussian model, X₀ given X₁ is Gaussian with conditional mean q X₁ + ω and conditional variance Σ independent of X₁. branch_actualization, branch_displacement and branch_variance return, respectively, q, ω and Σ. branch_precision and branch_variance should return a matrix of symmetric type. branch_variance defaults to the inverse of branch_precision. branch_logdet defaults to g = -0.5*log(|2πΣ|), the log normalizing constant of the Gaussian density in the traditional form. branch_transition_* return or modify in place the corresponding transition matrices.

Under a Brownian motion, we have q=I, ω=0, and conditional variance t*R where R is the model's variance rate and t the branch length.

branch_actualization(obj::EvolutionaryModel, edge::PN.Edge) 
branch_displacement(obj::EvolutionaryModel,  edge::PN.Edge)
branch_precision(obj::EvolutionaryModel,     edge::PN.Edge)
branch_variance(obj::EvolutionaryModel,      edge::PN.Edge)
branch_logdet(obj::EvolutionaryModel,        edge::PN.Edge, precision::AbstractMatrix)
branch_transition_qωjg(obj::EvolutionaryModel,    edge)
branch_transition_qωv!(q, obj::EvolutionaryModel, edge)

Under the most general linear Gaussian model, X₀ given X₁ is Gaussian with conditional mean q X₁ + ω and conditional variance Σ independent of X₁. branch_actualization, branch_displacement and branch_variance return, respectively, q, ω and Σ. branch_precision and branch_variance should return a matrix of symmetric type. branch_variance defaults to the inverse of branch_precision. branch_logdet defaults to g = -0.5*log(|2πΣ|), the log normalizing constant of the Gaussian density in the traditional form. branch_transition_* return or modify in place the corresponding transition matrices.

Under a Brownian motion, we have q=I, ω=0, and conditional variance t*R where R is the model's variance rate and t the branch length.

branch_actualization(obj::EvolutionaryModel, edge::PN.Edge) 
branch_displacement(obj::EvolutionaryModel,  edge::PN.Edge)
branch_precision(obj::EvolutionaryModel,     edge::PN.Edge)
branch_variance(obj::EvolutionaryModel,      edge::PN.Edge)
branch_logdet(obj::EvolutionaryModel,        edge::PN.Edge, precision::AbstractMatrix)
branch_transition_qωjg(obj::EvolutionaryModel,    edge)
branch_transition_qωv!(q, obj::EvolutionaryModel, edge)

Under the most general linear Gaussian model, X₀ given X₁ is Gaussian with conditional mean q X₁ + ω and conditional variance Σ independent of X₁. branch_actualization, branch_displacement and branch_variance return, respectively, q, ω and Σ. branch_precision and branch_variance should return a matrix of symmetric type. branch_variance defaults to the inverse of branch_precision. branch_logdet defaults to g = -0.5*log(|2πΣ|), the log normalizing constant of the Gaussian density in the traditional form. branch_transition_* return or modify in place the corresponding transition matrices.

Under a Brownian motion, we have q=I, ω=0, and conditional variance t*R where R is the model's variance rate and t the branch length.

calibrate!(beliefs::ClusterGraphBelief, schedule, niterations=1;
    auto::Bool=false, info::Bool=false,
    verbose::Bool=true,
    update_residualnorm::Bool=true,
    update_residualkldiv::Bool=false)

Propagate messages in postorder then preorder for each tree in the schedule list, for niterations. Each schedule "tree" should be a tuple of 4 vectors as output by spanningtree_clusterlist, where each vector provides the parent/child label/index of an edge along which to pass a message, and where these edges are listed in preorder. For example, the parent of the first edge is taken to be the root of the schedule tree. Calibration is evaluated after each schedule tree is run, and said to be reached if all message residuals have a small norm, based on iscalibrated_residnorm.

Output: true if calibration is reached, false otherwise.

Optional keyword arguments:

• auto: If true, then belief updates are stopped after calibration is found to be reached. Otherwise belief updates continue for the full number of iterations.
• info: Is true, information is logged with the iteration number and schedule tree index when calibration is reached.
• verbose: log error messages about degenerate messages
• update_residualnorm
• update_residualkldiv

See also: iscalibrated_residnorm and iscalibrated_residnorm! for the tolerance and norm used by default, to declare calibration for a given sepset message (in 1 direction).

calibrate!(beliefs::ClusterGraphBelief, scheduletree::Tuple,
    verbose::Bool=true, up_resnorm::Bool=true, up_reskldiv::Bool=false)

Propage messages along the scheduletree, in postorder then preorder: see propagate_1traversal_postorder!.

Output: true if all message residuals have a small norm, based on iscalibrated_residnorm, false otherwise.

calibrate_exact_cliquetree!(beliefs::ClusterGraphBelief,
    schedule,
    nodevector_preordered,
    tbl::Tables.ColumnTable, taxa::AbstractVector,
    evolutionarymodel_name)

For a Brownian Motion with a fixed root, compute the maximum likelihood estimate of the prior mean at the root and the restricted maximum likelihood (REML) estimate of the variance/covariance rate matrix using analytical formulas relying on belief propagation, using the data in tbl at leaves in the network. These estimates are for the model with a prior variance of 0 at the root, that is, a root state equal to the prior mean.

output: (bestmodel, loglikelihood_score) where bestmodel is an evolutionary model created by evolutionarymodel_name, containing the estimated model parameters.

assumptions:

• taxa should list the taxon names in the same order in which they come in the rows of tbl.
• schedule should provide a schedule to transmit messages between beliefs in beliefs (containing clusters first then sepsets). This schedule is assumed to traverse a clique tree for the input phylogeny, with the root cluster containing the root of the phylogeny in its scope.
• nodevector_preordered should list the nodes in this phylogeny, in preorder.
• beliefs should be of size and scope consistent with evolutionarymodel_name and data in tbl.
• a leaf should either have complete data, or be missing data for all traits.

Steps:

1. Calibrate beliefs in place according to the schedule, under a model with an infinite prior variance at the root.
2. Estimate parameters analytically.
3. Re-calibrate beliefs, to calculate the maximum log-likelihood of the fixed-root model at the estimated optimal parameters, again modifying beliefs in place. (Except that beliefs with the root node in scope are re-assigned to change their scoping dimension.)

Warning: there is no check that the beliefs and schedule are consistent with each other.

calibrate_optimize_cliquetree!(beliefs::ClusterGraphBelief, clustergraph,
    nodevector_preordered, tbl::Tables.ColumnTable, taxa::AbstractVector,
    evolutionarymodel_name, evolutionarymodel_startingparameters)

Optimize model parameters using belief propagation along clustergraph, assumed to be a clique tree for the input network, whose nodes in preorder are nodevector_preordered. Optimization aims to maximize the likelihood of the data in tbl at leaves in the network. The taxon names in taxa should appear in the same order as they come in tbl. The parameters being optimized are the variance rate(s) and prior mean(s) at the root. The prior variance at the root is fixed.

The calibration does a postorder of the clique tree only, to get the likelihood at the root without the conditional distribution at all nodes, modifying beliefs in place. Therefore, if the distribution of ancestral states is sought, an extra preorder calibration would be required. Warning: there is no check that the cluster graph is in fact a clique tree.

calibrate_optimize_clustergraph!(beliefs::ClusterGraphBelief, clustergraph,
    nodevector_preordered, tbl::Tables.ColumnTable, taxa::AbstractVector,
    evolutionarymodel_name, evolutionarymodel_startingparameters,
    max_iterations=100)

Same as calibrate_optimize_cliquetree! above, except that the user can supply an arbitrary clustergraph (including a clique tree) for the input network. Optimization aims to maximize the factored_energy approximation to the ELBO for the log-likelihood of the data (which is also the negative Bethe free_energy). When clustergraph is a clique tree, the factored energy approximation is exactly equal to the ELBO and the log-likelihood.

Cluster beliefs are regularized using regularizebeliefs_bycluster! (other options are likely to be available in future versions) before calibration. The calibration repeatedly loops through a minimal set of spanning trees (see spanningtrees_clusterlist) that covers all edges in the cluster graph, and does a postorder-preorder traversal for each tree. The loop runs till calibration is detected or till max_iterations have passed, whichever occurs first.

check_runningintersection(clustergraph, net)

Vector of tuples. Each tuple is of the form (nodelabel, istree), where nodelabel::Symbol is the label of a node in net and istree is true (false) if the node's cluster subgraph is (is not) a tree. This "cluster subgraph" for a given node is the subgraph of clustergraph induced by the clusters containing the node and by the edges whose sepset contain the node: see nodesubtree.

clustergraph satisfies the generalized running-intersection property if istree is true for all nodes in net.

Warning:

• assumes that net has been preordered, and
• does not check if clustergraph has been correctly constructed.

cliquetree(chordal_graph)

Clique tree U for an input graph g assumed to be chordal (triangulated), e.g. using triangulate_minfill!. Warning: does not check that the graph is already triangulated.

• Each node in U is a maximal clique of g whose data is the tuple of vectors (nodelabels, nodedata) using the labels and data from g, with nodes sorted by decreasing data. If g was originally built from a phylogenetic network using moralize, then the nodes' data are their preorder index, making them sorted in postorder within in each clique. The clique label is the concatenation of the node labels.
• For each edge (clique1, clique2) in U, the edge data hold the sepset (separating set) information as a vector of node data, for nodes shared by both clique1 and clique2. In this sepset, nodes are sorted by decreasing data.

Uses maximal_cliques and kruskal_mst (for min/maximum spanning trees) from Graphs.jl.

clustergraph!(net, method)
clustergraph(net, method)

Cluster graph U for an input network net and a method of cluster graph construction. The following methods are supported:

• Bethe
• LTRIP
• JoinGraphStructuring
• Cliquetree

The first method pre-processes net, which may modify it in place, see preprocessnet!. The second method assumes that net is already pre-processed.

clustergraph!(net, method)
clustergraph(net, method)

Cluster graph U for an input network net and a method of cluster graph construction. The following methods are supported:

• Bethe
• LTRIP
• JoinGraphStructuring
• Cliquetree

The first method pre-processes net, which may modify it in place, see preprocessnet!. The second method assumes that net is already pre-processed.

default_rootcluster(clustergraph, nodevector_preordered)

Index of a cluster that contains the network's root, whose label is assumed to be 1 (preorder index). If multiple clusters contain the network's root, then one is chosen with the smallest number of taxa (leaves in the network).

For cluster with label :lab, its property clustergraph[:lab][2] should list the nodes in the cluster, by the index of each node in nodevector_preordered such that 1 corresponds to the network's root. Typically, this vector is net.nodes_changed after the network is preordered.

default_rootcluster(clustergraph)

Index of a cluster that contains the node with the smallest preorder index. If multiple clusters contain that node, then one is chosen that only has that node. If all clusters containing that node have more than 1 node, then a cluster is chosen containing a node with the second-smallest-preorder index.

For cluster with label :lab, its property clustergraph[:lab][2] should list the nodes in the cluster by their preorder index, sorted in decreasingly order (so the smallest is at the end).

dimension(m::EvolutionaryModel)

Number of traits, e.g. 1 for univariate models.

entropy(J::Cholesky)
entropy(J::AbstractMatrix)
entropy(belief::AbstractBelief)

Entropy of a multivariate Gaussian distribution with precision matrix J, assumed to be square and symmetric (not checked). It is 0 if J is empty (of size 0×0). It may be Inf if J is semi-definite. The second version applies the first to the belief precision belief.J.

entropy is defined for discrete distributions in StatsBase.jl and extended to Gaussian distributions in Distributions.jl around here

factor_hybridnode(evolutionarymodel, ts::AbstractVector, γs)
factor_tree_degeneratehybrid(model,  t0::Real,           γs)

Canonical parameters h,J,g of factor ϕ(X₀, X₁, X₂, ...) from the evolutionary model for a hybrid node: where X₀ is the state at the hybrid node and X₁, X₂, ... the states of the parent nodes. Warning: γs is modified in placed, changed to [1 -γs].

It is assumed that the conditional mean is a simple weighted average:

$E[X_0 | X_1, X_2, ...] = \sum_k \gamma_k X_k = q \mathrm{vec}(X_1,X_2,...) + \omega$

where q has one block for each parent, and each block is diagonal scalar: $\gamma_k I_p$. More complex models could consider adding a shift ω to the conditional mean.

If all the parent hybrid edges edges have length 0, then it is assumed that the model gives a degenerate distribution, with 0 conditional variance. More complex models could consider adding a hybrid conditional variance Σ.

• The first form assumes that at least 1 parent edge length is positive, with conditional variance $\sum_k \gamma_k^2 V_k$ where $V_k$ is the conditional variance from the kth parent edge.
• The second form can be used in case all parent edges have 0 length, to integrate out the hybrid node state and the factor ϕ(X₀, X₁, X₂, ...) when X₀ is its child state, along an edge of length t0 between the hybrid node and its child. This second form is appropriate when this hybrid's child is a tree node, and t0>0.`

In h and J, the first p coordinates are for the hybrid (or its child) and the last coordinates for the parents, in the same order in which the edge lengths and γs are given.

factor_root(m::EvolutionaryModel)

Canonical parameters h,J,g of the prior density at the root, from model m. Assumes that isrootfixed(m) returns false (in which case the root value should be absorbed as evidence and the root removed from scope). More strongly, the root variance is assumed to be invertible, in particular, traits are all non-fixed at the root.

The prior is improper if the prior variance is infinite. In this case this prior is not a distribution (total probability ≠ 1) but is taken as the constant function 1, which corresponds to h,J,g all 0 (and an irrelevant mean).

If the root variance is not invertible (e.g., fixed root), this function fails and should never be called (see isrootfixed)

factor_treeedge(evolutionarymodel, edge)

Canonical parameters h,J,g of factor ϕ(X0,X1) from the given evolutionary model along one edge, where X₀ is the state of the child node and X₁ the state of the parent node. In h and J, the first p coordinates are for the child and the last p for the parent, where p is the number of traits (determined by the model).

Under the most general linear Gaussian model, X₀ given X₁ is Gaussian with conditional mean q X₁ + ω and conditional variance Σ independent of X₁. The generic fallback method uses functions branch_actualization for q, branch_displacement for ω, branch_precision for Σ⁻¹.

Under a Brownian motion, we have q=I, ω=0, and Σ=tR where R is the model's variance rate and t is the length of the branch. In that case, a specific (more efficient) method is implemented, and the default fallback is not used.

factored_energy(beliefs::ClusterGraphBelief)

Factored energy functional for general cluster graphs (Koller & Friedman 2009), which approximates the evidence lower bound (ELBO), a lower bound for the log-likelihood. It is also called the (negative) Bethe free energy in the context of factor graphs It is the sum of the cluster average energies and entropies, minus the sepset entropies. It is assumed but not checked that beliefs are calibrated (neighbor clusters and sepset beliefs are consistent, used as local marginals).

For a calibrated clique tree, the factored energy is equal to the log-likelihood. For a calibrated cluster graph, it can serve as as approximation.

output: tuple of 3 values, the 3rd being the factored energy: (average energy, approximate entropy, factored energy = -energy + entropy).

See also: free_energy, entropy, average_energy!

References

D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT Press, 2009. ISBN 9780262013192.

D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A Review for Statisticians, Journal of the American statistical Association, 112:518, 859-877, 2017, doi: 10.1080/01621459.2017.1285773.

free_energy(beliefs::ClusterGraphBelief)

negative factored_energy to approximate the negative log-likelihood. The approximation is exact on a clique tree after calibration.

getcholesky(J::AbstractMatrix)

Cholesky decomposition of J, assumed to be symmetric positive definite, stored as a PDMat object. Warning: PDMat is not a subtype of Cholesky. PDMats.jl is efficient for structured matrices (e.g diagonal or sparse) and has efficient methods for linear algebra, e.g. \, invquad, X_invA_Xt etc.

getcholesky_μ(J::AbstractMatrix, h)
getcholesky_μ!(belief::Belief)

Tuple (Jchol, μ) where Jchol is a cholesky representation of J or belief.J and μ is J⁻¹h, used to update belief.μ (by the second method).

getcholesky_μ(J::AbstractMatrix, h)
getcholesky_μ!(belief::Belief)

Tuple (Jchol, μ) where Jchol is a cholesky representation of J or belief.J and μ is J⁻¹h, used to update belief.μ (by the second method).

hasdegenerate(net)

true if degenerate nodes remain in scope, that is, if there exists a tree edge of length 0, or if there exists a hybrid node with all parent edges of length 0 and with 2 or more children edges, or with 1 child edge of length 0.

hybridnode_displacement(obj::EvolutionaryModel, parentedges::AbstractVector{PN.Edge})
hybridnode_precision(obj::EvolutionaryModel,    parentedges::AbstractVector{PN.Edge})
hybridnode_variance(obj::EvolutionaryModel,     parentedges::AbstractVector{PN.Edge})

Under the most general weighted average Gaussian model, X₀ given its parents X₁, X₂, ... is Gaussian with conditional mean the weighted average of the parents plus a displacement vector ω and conditional variance Σ independent of X₁, X₂, ... . The weights are given by the inheritance probabilities contained in the PN.Edge objects. hybridnode_displacement and hybridnode_variance return, respectively, ω and Σ.

hybridnode_variance and hybridnode_precision should return a matrix of symmetric type. hybridnode_precision defaults to the inverse of hybridnode_variance. hybridnode_displacement and hybridnode_variance default to a vector or matrix of zeros.

hybridnode_displacement(obj::EvolutionaryModel, parentedges::AbstractVector{PN.Edge})
hybridnode_precision(obj::EvolutionaryModel,    parentedges::AbstractVector{PN.Edge})
hybridnode_variance(obj::EvolutionaryModel,     parentedges::AbstractVector{PN.Edge})

Under the most general weighted average Gaussian model, X₀ given its parents X₁, X₂, ... is Gaussian with conditional mean the weighted average of the parents plus a displacement vector ω and conditional variance Σ independent of X₁, X₂, ... . The weights are given by the inheritance probabilities contained in the PN.Edge objects. hybridnode_displacement and hybridnode_variance return, respectively, ω and Σ.

hybridnode_variance and hybridnode_precision should return a matrix of symmetric type. hybridnode_precision defaults to the inverse of hybridnode_variance. hybridnode_displacement and hybridnode_variance default to a vector or matrix of zeros.

hybridnode_displacement(obj::EvolutionaryModel, parentedges::AbstractVector{PN.Edge})
hybridnode_precision(obj::EvolutionaryModel,    parentedges::AbstractVector{PN.Edge})
hybridnode_variance(obj::EvolutionaryModel,     parentedges::AbstractVector{PN.Edge})

Under the most general weighted average Gaussian model, X₀ given its parents X₁, X₂, ... is Gaussian with conditional mean the weighted average of the parents plus a displacement vector ω and conditional variance Σ independent of X₁, X₂, ... . The weights are given by the inheritance probabilities contained in the PN.Edge objects. hybridnode_displacement and hybridnode_variance return, respectively, ω and Σ.

hybridnode_variance and hybridnode_precision should return a matrix of symmetric type. hybridnode_precision defaults to the inverse of hybridnode_variance. hybridnode_displacement and hybridnode_variance default to a vector or matrix of zeros.

init_beliefs_allocate(tbl::Tables.ColumnTable, taxa, net, clustergraph,
                      evolutionarymodel)

Vector of beliefs, initialized to the constant function exp(0)=1, one for each cluster then one for each sepset in clustergraph. tbl is used to know which leaf in net has data for which trait, so as to remove from the scope each variable without data below it. taxa should be a vector with taxon names in the same order as they come in the table of data tbl. The root is removed from scope if the evolutionary model has a fixed root: so as to use the model's fixed root value as data if the root as zero prior variance. Also removed from scope is any hybrid node that is degenerate and who has a single child edge of positive length.

Warnings: this function might need to be re-run to re-do allocation if

• the data changed: different number of traits, or different pattern of missing data at the tips
• the model changed: with the root changed from fixed to random, see init_beliefs_allocate_atroot! in that case.

init_beliefs_allocate_atroot!(beliefs, factors, messageresiduals, model)

Update the scope and re-allocate memory for cluster & sepset beliefs, factors and messageresiduals to include or exclude the root, depending on whether the root variable is random or fixed in model. To change the dimension of canonical parameters μ,h,J, new memory is allocated and initilized to 0. This function can be used to update beliefs when the root model changes from fixed to non-fixed or vice-versa. It re-allocates less memory than init_beliefs_allocate (which would need to be followed by init_factors_allocate) because clusters and sepsets that do not have the root are not modified.

Assumptions:

• all traits at the root have at least one descendant with non-missing data,
• beliefs were previously initialized with a model that had the same number of traits as the current model.

init_beliefs_assignfactors!(beliefs,
                            evolutionarymodel, columntable, taxa,
                            nodevector_preordered)

Initialize cluster beliefs prior to belief propagation, by assigning each factor to one cluster. Sepset beliefs are reset to 0. There is one factor for each node v in the vector of nodes: the density of X_v conditional on its parent X_pa(v) if v is not the root, or the prior density for X_root.

• for each leaf, the factor is reduced by absorbing the evidence for that leaf, that is, the data found in the columntable, whose rows should be ordered by taxa as they appear in taxa.
• for each leaf, missing trait values are removed from scope.
• for each internal node, any trait not in scope (e.g. if all descendant leaves are missing a value for this trait) is marginalized out of the factor.

Assumptions:

• In vector nodevector_preordered, nodes are assumed to be preordered. Typically, this vector is net.nodes_changed after the network is preordered.
• Belief node labels correspond to the index of each node in nodevector_preordered.
• In beliefs, cluster beliefs come first and sepset beliefs come last, as when created by init_beliefs_allocate

Output: vector node2belief such that, if i is the preorder index of a node in the network, node2belief[i] is the index of the belief that the node family was assigned to.

The beliefs vector is modified in place.

init_beliefs_reset!(beliefs::Vector{<:Belief})

Reset all beliefs (which can be cluster and/or sepset beliefs) to h=0, J=0, g=0 (μ unchanged). They can later be re-initialized for different model parameters and re-calibrated, without re-allocating memory.

init_beliefs_reset_fromfactors!(beliefs::ClusterGraphBelief)

Reset cluster beliefs to existing factors, and sepset beliefs to h=0, J=0, g=0

This is not used so far, as changing model parameters requires a reset of both factors and beliefs, done by init_beliefs_assignfactors!.

init_clustergraph(T::Type{<:Integer}, clustergraph_method::Symbol)

MetaGraph with an empty base graph (0 vertices, 0 edges), meta-graph data clustergraph_method, edge-weight function counting the length of the edge-data vector, and the following types:

• vertex indices in the base graph: T
• vertex labels: Symbol
• vertex data: Tuple{Vector{Symbol}, Vector{T}} to hold information about the variables (nodes in phylogenetic network) in the cluster (vertex in cluster graph): node names as symbols, and node preorder indices
• edge data: Vector{T} to hold information about the sepset: preorder index of nodes in the sepset.

See packages MetaGraphsNext and Graphs.

The empty graph above is of type Graphs.SimpleGraphs.SimpleGraph{T}: undirected, with vertex indices of type T. After addition of n vertices, the vertex indices range from 1 to n, technically in Base.OneTo{T}(n).

init_factors_allocate(beliefs::AbstractVector{<:Belief}, nclusters::Integer)

Vector of nclusters factors of type FamilyFactor, whose canonical parameters and metadata are initialized to be a copy of those in beliefs. Assumption: beliefs[1:nclusters] are cluster beliefs, and beliefs[nclusters+1:end] (if any) are sepset beliefs. This is not checked.

init_factors_frombeliefs!(factors, beliefs, checkmetadata::Bool=false)

Reset all factors by copying h,J,g from beliefs. Assumption: the cluster beliefs match the factors exactly: for a valid factor index i, beliefs[i] is of cluster type and has the same dimension as factors[i].

Set checkmetadata to true to check that beliefs[i] and factors[i] have the same metadata.

init_messagecalibrationflags_reset!(beliefs::ClusterGraphBelief, reset_kl=true)

Reset all non-empty message residuals in beliefs.

init_messagecalibrationflags_reset!(mr::AbstractResidual, reset_kl::Bool)

For a non-empty message residual mr, reset its iscalibrated_* flags to false, and if reset_kl is true, reset its kldiv to -1. Its ΔJ and Δh fields are not reset here, because they are overwritten during a belief propagation step. Nothing is done for empty messages.

init_messageresidual_allocate(beliefs::Vector{B}, nclusters)

Dictionary of 2k residuals of type MessageResidual, whose canonical parameters (Δh,ΔJ) are initialized using MessageResidual, to be of the same size as sepsets in beliefs, where k is length(beliefs) - nclusters. Assumption: the first nclusters beliefs are cluster beliefs, and the next k beliefs are sepset beliefs. This is not checked.

The sepset for edge (label1,label2) is associated with 2 messages, for the 2 directions in which beliefs can be propagated along the edge. The keys for these messages are (label1,label2) and (label2,label1).

inscope_onenode(node_label, b:Belief)

AbstractVector: view of the row vector in b's inscope matrix corresponding to node node_label, indicating whether a trait at that node is in scope or not.

integratebelief!(belief::AbstractBelief)
integratebelief(h,J,g)

(μ,g) from fully integrating the belief, that is: $μ = J^{-1} h$ and $g + (\log|2\pi J^{-1}| + h^{T} J^{-1} h)/2 .$ The first form updates belief.μ.

integratebelief!(obj::ClusterGraphBelief, beliefindex)
integratebelief!(obj::ClusterGraphBelief)
integratebelief!(obj::ClusterGraphBelief, clustergraph, nodevector_preordered)

(μ,g) from fully integrating the object belief indexed beliefindex. This belief is modified, with its belief.μ's values updated to those in μ.

The second method uses the first sepset containing a single node. This is valid if the beliefs are fully calibrated (including a pre-order traversal), but invalid otherwise. The third method uses the default cluster containing the root, see default_rootcluster. This is valid if the same cluster was used as the root of the cluster graph, if this graph is a clique tree, and after a post-order traversal to start the calibration.

iscalibrated_kl!(res::AbstractResidual, atol=1e-5)

True if the KL divergence stored in res.kldiv is within atol of 0; false otherwise. res.iscalibrated_kl is modified accordingly. This KL divergence should have been previously calculated: between a sepset belief, equal to the message that was passed most recently, and its belief just prior to passing that message.

iscalibrated_residnorm!(res::AbstractResidual, atol=1e-5, p::Real=Inf)

True if the canonical parameters res.Δh and res.ΔJ of the message residual have p-norm within atol of 0; false otherwise. res.iscalibrated_resid is updated accordingly.

With p infinite, the max norm is used by default, meaning that res.Δh and res.ΔJ should be close to 0 element-wise.

iscalibrated_residnorm(beliefs::ClusterGraphBelief)
iscalibrated_kl(beliefs::ClusterGraphBelief)

True if all edges in the cluster graph have calibrated messages in both directions, in that their latest message residuals have norm close to 0 (residnorm) or KL divergence close to 0 between the message received and prior sepset belief. False if not all edges have calibrated messages.

This condition is sufficient but not necessary for calibration.

Calibration was determined for each individual message residual by iscalibrated_residnorm! and iscalibrated_kl! using some tolerance value.

isdegenerate(node)

true if all parent edges of node have length 0, false otherwise. Intended for hybrid nodes, as tree edges of length 0 should be suppressed before trait evolution analysis.

isfamilypreserving(clusters, net)

Tuple (ispreserving, isfamily_incluster):

1. ispreserving: true (false) if clusters is (is not) family-preserving with respect to net, that is: if each node family (a node and all of its parents) in net is contained in at least 1 cluster in clusters. clusters should be a vector, where each element describes one cluster, given as a vector of preorder indices. Index i corresponds to node number i in net according the node pre-ordering: net.nodes_changed[i].
2. isfamily_incluster: vector of BitVectors. isfamily_incluster[i][j] is true (false) if the family of node i is (is not) fully contained in cluster [j]. i is taken as the preorder index of a node in net. ispreserving is true if every element (bit vector) in isfamily_incluster contains at least 1 true value.

Warning: assumes that net is preordered, see PhyloNetworks.preorder!).

See also nodefamilies to get node families.

ishybridsinglepositivechild(node)

true if node is a hybrid node with a single child edge of positive length. If it isdegenerate) (all its parent edges have length 0) and if its child is a tree node, then it could be removed from scope: see unscope.

joingraph(net, method::JoinGraphStructuring)

See JoinGraphStructuring

ltripclustergraph(net, method::LTRIP)

See LTRIP

marginalizebelief(belief::AbstractBelief, keep_index)
marginalizebelief(h,J,g, keep_index, beliefmetadata)
marginalizebelief(h,J,g, keep_index, integrate_index, beliefmetadata)

Canonical form (h,J,g) of the input belief, after all variables except those at indices keep_index have been integrated out. If we use I and S subscripts to denote subvectors and submatrices at indices to integrate out (I: integrate_index) and indices to keep (S: save for sepset, keep_index) then the returned belief parameters are:

$h_S - J_{S,I} J_I^{-1} h_I$

$J_S - J_{S,I} J_I^{-1} J_{I,S}$

and

$g + (\log|2\pi J_I^{-1}| + h_I^{T} J_I^{-1} h_I)/2 .$

These operations fail if the Cholesky decomposition of $J_I$ fails. In that case, an error of type BPPosDefException is thrown with a message about the beliefmetadata, which can be handled by downstream functions.

moralize!(net::HybridNetwork, prefix="I")
moralize(net)

Undirected graph g of type MetaGraph with the same nodes as in net, labelled by their names in net, with extra edges to moralize the graph, that is, to connect any nodes with a common child. Node data, accessed as g[:nodelabel], is their index in the network's preordering. Edge data, accessed as g[:label1, :label2] is a type to indicate if the edge was an original tree edge or hybrid edge, or added to moralize the graph. Another type, not used here, if for fill edges that may need to be added to triangulate g (make it chordal).

The first version modifies net to name its internal nodes (used as labels in g) and to create or update its node preordering, then calls the second version.

moralize!(net::HybridNetwork, prefix="I")
moralize(net)

Undirected graph g of type MetaGraph with the same nodes as in net, labelled by their names in net, with extra edges to moralize the graph, that is, to connect any nodes with a common child. Node data, accessed as g[:nodelabel], is their index in the network's preordering. Edge data, accessed as g[:label1, :label2] is a type to indicate if the edge was an original tree edge or hybrid edge, or added to moralize the graph. Another type, not used here, if for fill edges that may need to be added to triangulate g (make it chordal).

The first version modifies net to name its internal nodes (used as labels in g) and to create or update its node preordering, then calls the second version.

nodefamilies(net)

Vector v with elements of type Vector{T}. v[i] first lists i, the preorder index of node net.nodes_changed[i], followed by the preorder index of all of this node's parents in net, sorted in decreasing order. Due to pre-ordering, all of the parents' indices are listed after the node (their child) index. A given node and its parents is called a "node family".

Warning: net is assumed preordered, see preprocessnet! and PhyloNetworks.preorder!).

nodesubtree(clustergraph::MetaGraph, node_symbol)
nodesubtree(clustergraph::MetaGraph, node_symbol, node_preorderindex)

MetaGraph subgraph of clustergraph induced by the clusters and sepsets containing the node labelled node_symbol (of specified preorder index if known, not checked). If clustergraph satisfies the generalized running-intersection property, then this subgraph should be a tree.

nodesubtree_clusterlist(clustergraph::MetaGraph, nodesymbol)

Spanning tree of the subgraph of clustergraph induced by the clusters and sepsets that contain the node labelled nodesymbol, see nodesubtree. If clustergraph satisfies the generalized running-intersection property, then this subgraph should be a tree anyway, but this is not assumed (via extracting a spanning tree).

Output: tuple of 4 vectors describing a depth-first search traversal of the tree, starting from a cluster containing the node of smallest preorder index, as determined by default_rootcluster. Each element in this tuple is a vector: see spanningtree_clusterlist.

params(m::EvolutionaryModel)

Tuple of parameters, the same that can be used to construct the evolutionary model.

params_optimize(m::EvolutionaryModel)

Tuple of transformed parameters for model m, in an unconstrained space that can be used for numerical optimization.

params_original(m::EvolutionaryModel, transformedparams::AbstractArray)

Tuple of parameters for model m in the original space, corresponding to back-transformed parameters of transformedparams, and that can be used as input to the model constructor.

parentinformation(node, net)

Tuple of (edge length, edge γ, index of parent node in net.nodes_changed) for all parent edges of node. Assumes that net has been preordered before.

preprocessnet!(net::HybridNetwork, prefix="I")

Create or update the pre-ordering of nodes in net using PhyloNetworks.preorder!, then name unnamed internal nodes, with names starting with prefix. Nodes in phylogenetic networks need to have names to build cluster graphs, in which a cluster contains network nodes. Pre-ordering is also used to traverse the network for building cluster graphs.

See clustergraph!.

propagate_1traversal_postorder!(beliefs::ClusterGraphBelief, spanningtree...)
propagate_1traversal_preorder!(beliefs::ClusterGraphBelief,  spanningtree...)

Messages are propagated along the spanning tree, from the tips to the root by propagate_1traversal_postorder! and from the root to the tips by propagate_1traversal_preorder!.

The "spanning tree" should be a tuple of 4 vectors as output by spanningtree_clusterlist, meant to list edges in preorder. Its nodes (resp. edges) should correspond to clusters (resp. sepsets) in beliefs: labels and indices in the spanning tree information should correspond to indices in beliefs. This condition holds if beliefs are produced on a given cluster graph and if the tree is produced by spanningtree_clusterlist on the same graph.

Optional positional arguments after spanning tree, in this order (default value):

• verbose (true): log error messages about degenerate messages that failed to be passed.
• update_residualnorm (true): to update each message residual's iscalibrated_resid
• update_residualkldiv (false): to update each message residual's field kldiv: KL divergence between the new and old sepset beliefs, normalized to be considered as (conditional) distributions.

propagate_belief!(cluster_to, sepset, cluster_from, residual)

Update the canonical parameters of the beliefs in cluster_to and in sepset, by marginalizing the belief in cluster_from to the sepset's variable and passing that message. The change in sepset belief (Δh and ΔJ: new - old) is stored in residual.

Degeneracy

Propagating a belief requires the cluster_from belief to have a non-degenerate J_I: submatrix of J for the indices to be integrated out. Problems arise if this submatrix has one or more 0 eigenvalues, or infinite values (see marginalizebelief). If so, a BPPosDefException is returned but not thrown. Downstream functions should try & catch these failures, and decide how to proceed. See regularizebeliefs_bycluster! to reduce the prevalence of degeneracies.

Output

• nothing if the message was calculated with success
• a BPPosDefException object if marginalization failed. In this case, the error is not thrown: downstream functions should check for failure (and may choose to throw the output error object).

Warnings

• only the h, J and g parameters are updated, not μ.
• Does not check that cluster_from and cluster_to are of cluster type, or that sepset is of sepset type, but does check that the labels and scope of sepset are included in each cluster.

regularizebeliefs_1clustersepset!(cluster::AbstractBelief, sepset::AbstractBelief, ϵ)

Modify beliefs of 1 cluster and 1 of sepset, assumed to be neighbors in a cluster graph (such that the sepset's scope is a subset of the cluster's scope) by adding ϵ to all diagonal elements of the sepset precision matrice J and to the corresponding diagonal elements of the cluster precision, so as to preserve the full graphical model.

Used by regularizebeliefs_bycluster! and regularizebeliefs_onschedule!.

regularizebeliefs_bycluster!(beliefs::ClusterGraphBelief, clustergraph)
regularizebeliefs_bycluster!(beliefs::ClusterGraphBelief, clustergraph, cluster_label)

Modify beliefs of cluster graph by adding positive values to some diagonal elements of precision matrices J, while preserving the full graphical model (product of cluster beliefs over product of sepset beliefs, invariant during belief propagation) so that all beliefs are non-degenerate.

This regularization could be done after initialization with init_beliefs_assignfactors! for example.

The goal is that at each later step of belief propagation, the sending cluster has a non-degenerate (positive definite) precision matrix for the variables to be integrated, so that the message to be sent is well-defined (i.e. can be computed) and positive semidefinite.

Algorithm

For each cluster Ci (or for only for 1 cluster, labeled cluster_label):

1. Find a regularization parameter adaptively for that cluster: ϵ = maximum absolute value of all entries in Ci's precision matrix J, and of the machine epsilon. Then loop through its incident edges:
2. For each neighbor cluster Cj and associated sepset Sij, add ϵ > 0 to the diagonal entries of Ci's precision matrix J corresponding to the traits in Sij.
3. To preserve the graphical model's joint distribution for the full set of variables (invariant during BP), the same ϵ is added to each diagonal entry of Sij's precision matrix.

regularizebeliefs_bynodesubtree!(beliefs::ClusterGraphBelief, clustergraph)

Modify beliefs of cluster graph by adding positive values to some diagonal elements of precision matrices J, while preserving the full graphical model (product of cluster beliefs over product of sepset beliefs, invariant during belief propagation) so that all beliefs are non-degenerate.

The goal is the same as regularizebeliefs_bycluster! and regularizebeliefs_onschedule!, but the algorithm is different.

Algorithm

For each node (or variable) v:

1. Consider the subgraph T of clusters & edges that have v. If clustergraph has the generalized running-intersection property, this subgraph is a tree.
2. Root T at a cluster containing a node with the largest postorder index.
3. Find a regularization parameter adaptively for that node: ϵ = maximum absolute value of all entries in Ci's precision matrix J, and of the machine epsilon, over clusters Ci in T.
4. For each trait j, find the subtree Tj of clusters and sepsets that have trait j of node v in their scope.
5. For each cluster and sepset in Tj, except at its cluster root: add ϵ on the diagonal of their belief precision matrix J corresponding to trait j of node v.
6. Check that graphical model invariant is preserved, that is: for each trait j, ϵ was added to the same number of clusters as number of sepsets.

regularizebeliefs_onschedule!(beliefs::ClusterGraphBelief, clustergraph)

Modify beliefs of the cluster graph so that all beliefs are non-degenerate, by (1) adding positive values to some diagonal elements of precision matrices J while preserving the full graphical model, and (2) propagating messages by marginalizing cluster beliefs.

The goal is the same as regularizebeliefs_bycluster! and regularizebeliefs_bynodesubtree!, but the algorithm is different.

Algorithm

The outcomes depends on a "schedule", that is, on an order in which clusters are considered. Here, the method simply takes the order given by labels(clustergraph). Mark all messages Ci → Cj as un-sent. For each cluster Ci in order:

1. For each neighbor cluster Cj of Ci, if the message Cj → Ci has not been sent, then:
   • add ϵ_i to the diagonal entries of sepset Sij's precision matrix Jij
   • add ϵ_i to the diagonal entries of Ci's precision matrix Ji that correspond to the inscope variables of Sij
   • mark Cj → Ci as sent.
   Notes:
   • ϵ_i is taken as the maximum of all absolute values of entries in Ji, and of the square root of machine epsilon.
   • This is equivalent to "sending" a default message, that is not computed from Cj's belief, with diagonal precision matrix ϵ⋅I and other canonical parameters h=0, g=0.
2. For each neighbor cluster Cj of Ci, if Ci → Cj has not been sent then
   • send a message from Ci to Cj by marginalizing Ci's belief (that is, using belief propagation).
   • mark the message Ci → Cj as sent.

residual_kldiv!(residual::AbstractResidual, sepset::AbstractBelief,
    canonicalparams::Tuple)

Update residual.kldiv with the Kullback-Leibler divergence between a message sent through a sepset (normalized to a probability distribution), and the sepset belief before the belief update (also normalized). sepset should contain the updated belief, and residual the difference in the J and h parameters due to the belief update (after - before), such that the previous belief is: sepset belief - residual. As a side product, sepset.μ is updated.

Output: true if the KL divergence is close to 0, false otherwise. See iscalibrated_kl! for the tolerance.

If the current or previous sepset belief is degenerate, in the sense that its precision matrix is not positive definite and the belief cannot be normalized to a proper distribution, then residual and sepset are not updated, and false is returned. No warning and no error is sent, because sepset beliefs are initialized at 0 and this case is expected to be frequent before enough messages are sent.

Calculation:

sepset after belief-update (i.e. message sent): C(x | Jₘ, hₘ, gₘ) ∝ density for X ~ 𝒩(μ=Jₘ⁻¹hₘ, Σ=Jₘ⁻¹) sepset before belief-update: C(x | Jₛ, hₛ, gₛ) residual: ΔJ = Jₘ - Jₛ, Δh = hₘ - hₛ p: dimension of X (number of variables: number of nodes * number of traits). Below, we use the notation Δg for the change in constants to normalize each message, which is not gₘ-gₛ because the stored beliefs are not normalized.

KL(C(Jₘ, hₘ, _) || C(Jₛ, hₛ, _))
= Eₘ[log C(x | Jₘ,hₘ,_)/C(x | Jₛ,hₛ,_)] where x ∼ C(Jₘ,hₘ,_)
= Eₘ[-(1/2) x'ΔJx + Δh'x + Δg)]
= (  tr(JₛJₘ⁻¹) - p + (μₛ-μₘ)'Jₛ(μₛ-μₘ) + log(det(Jₘ)/det(Jₛ))  ) /2

See also: average_energy!, which only requires the sepset belief to be positive definite.

scopeindex(node_label, sepset::AbstractBelief, cluster::AbstractBelief)

Tuple of 2 index vectors (ind_in_sepset, ind_in_cluster) in the sepset and in the cluster in-scope variables (that is, in the cluster's μ,h,J vectors and matrices) of the shared in-scope traits for node node_label, such that sepset.μ[ind_in_sepset] correspond to all the node's traits in the sepset scope and cluster.μ[ind_in_cluster] correspond to the same traits in the cluster scope, which may be a subset of all the node's traits in scope for that cluster. If not, an error is thrown.

scopeindex(j::Integer, belief::AbstractBelief)

Indices in the belief's μ,h,J vectors and matrices of the traits in scope for node indexed j in nodelabels(belief).

scopeindex(sepset::AbstractBelief, cluster::AbstractBelief)

Indices ind, in the cluster in-scope variables (that is, in the cluster's μ,h,J vectors and matrices) of the sepset in-scope variables, such that cluster.μ[ind] correspond to the same variables as sepset.μ, for example. These sepset in-scope variables can be a subset of traits for each node in the sepset, as indicated by inscope(sepset).

Warning: the labels in the sepset are assumed to be ordered as in the cluster.

An error is returned if sepset contains labels not in the cluster, or if a variable in the sepset's scope is not in scope in the cluster.

scopeindex(node_labels::Union{Tuple,AbstractVector}, belief::AbstractBelief)

Indices in the belief's μ,h,J vectors and matrices of the variables for nodes labeled node_labels. The belief's inscope matrix of booleans says which node (column) and trait (row) is in the belief's scope. These variables are vectorized by stacking up columns, that is, listing all in-scope traits of the first node, then all in-scope traits of the second node etc.

shrinkdegenerate_treeedges(net::HybridNetwork)

Network obtained from net with any non-external tree edge of length 0 suppressed. Returns an error if any edge length is missing or negative, or if any γ is missing or non positive. It is assumed that γs sum to 1 across partner hybrid edges.

spanningtree_clusterlist(clustergraph, root_index)
spanningtree_clusterlist(clustergraph, nodevector_preordered)

Build the depth-first search spanning tree of the cluster graph, starting from the cluster indexed root_index in the underlying simple graph; find the associated topological ordering of the clusters (preorder); then return a tuple of these four vectors:

1. parent_labels: labels of the parents' child clusters. The first one is the root.
2. child_labels: labels of clusters in pre-order, except for the cluster choosen to be the root.
3. parent_indices: indices of the parent clusters
4. child_indices: indices of the child clusters, listed in preorder as before.

In the second version in which root_index is not provided, the root of the spanning tree is chosen to be a cluster that contains the network's root. If multiple clusters contain the network's root, then one is chosen containing the smallest number of taxa: see default_rootcluster.

spanningtrees_clusterlist(clustergraph, nodevector_preordered)

Vector of spanning trees for clustergraph, that together cover all edges. Each spanning tree is specified as a tuple of 4 vectors describing a depth-first search traversal of the tree, starting from a cluster that contains the network's root, as in spanningtree_clusterlist.

The spanning trees are iteratively obtained using Kruskal's minimum-weight spanning tree algorithm, with edge weights defined as the number of previous trees covering the edge.

triangulate_minfill!(graph)

Ordering for node elimination, chosen to greedily minimize the number of fill edges necessary to eliminate the node (to connect all its neighbors with each other). Ties are broken by favoring the post-ordering of nodes. graph is modified with these extra fill edges, making it chordal.

unscope(node)

true if node is a hybrid node with a single child edge of positive length, and if its child node is a tree node. If it isdegenerate) (all its parent edges have length 0) then it could be removed from scope: see addtreenode_belowdegeneratehybrid!.