struct ManifoldKernelDensity{M<:AbstractManifold, B<:(Union{var"#s24801", var"#s24800"} where {var"#s24801"<:ApproxManifoldProducts.ManellicTree, var"#s24800"<:BallTreeDensity}), L, P<:AbstractArray}

On-manifold kernel density belief.

Notes

• Allows partials as identified by list of coordinate dimensions e.g. ._partial = [1;3]
  • When building a partial belief, use full points with necessary information in the specified partial coords.

DevNotes

• WIP AMP issue 41, use generic retractions during manifold products.

Likely new bug see KDE #70

ad_lie(M, X)

Construct generic adjoint matrix for Lie group manifolds. Input X is a tangent vector of the Lie algebra of group manifold M.

Notes

• This implementation uses the Lie bracket over affine (or screw) matrices.
  • Ref [Chirikjian 2012, Vol.2, pg.30, eq.10.59b]

antimarginal(newM, u0, mkd, newpartial)

If a marginal (statistics) of a probability reduce the dimensions (i.e. casts a shadow, or projects); then an antimarginal aims to increase dimension of the probability within reason.

Notes

• Marginalization is a integration of for higher dimension to lower dimension, so antimarginal likely involve differentiation (anti-integral) instead.
• In manifold language, this is an embedding into a higher dimensional space.
• See structure from motion in machine vision, or stereo disparity for building depth clouds from 2D images.
• Imagine combining three different partial embedding A=[1, nan, nan, 1.4], B=[nan,2.1,nan,4], C=[nan, nan, 3, nan] which should equal A+B+C = [notnan, notnan, notnan, notnan]

buildHybridManifoldCallbacks(manif)

Lots to do here, see RoME.jl #244 and standardized usage with Manifolds.jl.

Notes

• diffop( test, reference ) <===> ΔX = inverse(test) * reference

DevNotes

• FIXME replace with Manifolds.jl #41, RoME.jl #244

buildTree_Manellic!(M, r_PP; N, weights, kernel, kernel_bw)

Notes:

• Bandwidths for leaves (i.e. kernel_bw) must be passed in as covariances when MvNormalKernel.

DevNotes:

• Design Decision 24Q1, Manellic.MvNormalKernel bandwidth defs should ALWAYS ONLY BE covariances, because
  • Vision state is multiple bandwidth kernels including off diagonals in both tree or leaf kernels
  • Hybrid parametric to leafs convariance continuity
  • https://github.com/JuliaStats/Distributions.jl/blob/a9b0e3c99c8dda367f69b2dbbdfa4530c810e3d7/src/multivariate/mvnormal.jl#L220-L224

calcMean(mkd)
calcMean(mkd, aspartial)

Alias for overloaded Statistics.mean.

calcProductGaussians(
    M,
    kernels;
    μ0,
    weight,
    do_transport_correction
)

EXPERIMENTAL: On-manifold product of Gaussians.

DevNotes

• FIXME do product of concentrated Gaussians on Lie group (approximation):
  • See Section 3.2 and 4 of [Ge, van Goor, Mahony: A Geometric Perspective on using Gaussian Distributions on Lie Groups, 2024].
  • Also see upstream utils, https://juliamanifolds.github.io/Manifolds.jl/stable/features/distributions.html
• FIXME is parallel transport needed when multiplying with covariances from difffent tangent spaces?

calcProductGaussians(
    M,
    μ_,
    Σ_;
    μ0,
    Λ_,
    weight,
    do_transport_correction
)

Calculate covariance weighted mean as product of incoming Gaussian points $μ_$ and coordinate covariances $Σ_$.

Notes

• Return both weighted mean and new covariance (teh congruent product)
• More efficient helper function allows passing keyword inverse covariances Λ_ instead.
• Assume size(Σ_[1],1) == manifold_dimension(M).
• calc lambdas first and use to calculate mean product second.
• https://ccrma.stanford.edu/~jos/sasp/ProductTwoGaussian_PDFs.html
• Pennec, X. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements, HAL Archive, 2011, Inria, France.

DevNotes:

• FIXME is parallel transport needed as products involve covariances from different tangent spaces?
• TODO avoid recomputing covariance matrix inverses all the time

calcProductKernelBTLabels(M, proposals, labels_sampled; ...)
calcProductKernelBTLabels(
    M,
    proposals,
    labels_sampled,
    LOOidx;
    ...
)
calcProductKernelBTLabels(
    M,
    proposals,
    labels_sampled,
    LOOidx,
    gibbsSeq;
    permute,
    weight
)

Calculate one product of proposal kernels, as defined BTLabels.

covTransformNormalized(Σ)

Transform T=RS from unit covariance D to instance covariance Σ = TD.

Notes:

• Geometric interpretation of the covariance matrix, Fig. 10, https://users.cs.utah.edu/~tch/CS6640F2020/resources/A%20geometric%20interpretation%20of%20the%20covariance%20matrix.pdf
  • Eigen decomp: Σ^2 V = VL => Σ^2 = VL(V^-1) = RL(R^-1) = RSS(R^-1) => T=RS

evaluateDensityAtPoints(M, density, eval_at_points)
evaluateDensityAtPoints(
    M,
    density,
    eval_at_points,
    normalize
)

Return vector of weights of evaluated proposal label points against density.

DevNotes:

• TODO should evat points be of equal weights? If multiscale sampling goes down unbalanced trees?
• FIXME how should partials be handled here?
• FIXME, use multipoint evaluation such as NN (not just one point at a time)

evaluateManifoldNaive1D!(ret, idx, pts, bw, x)
evaluateManifoldNaive1D!(ret, idx, pts, bw, x, loo)
evaluateManifoldNaive1D!(ret, idx, pts, bw, x, loo, diffop)

Evalute the KDE naively as equally weighted Gaussian kernels with common bandwidth. This function does, however, allow on-manifold evaluations.

getBandwidthSearchBounds(mtree)

For Manellic tree parent kernels, what is the 'smallest' and 'biggest' covariance.

Notes:

• Thought about det for covariance volume but long access of pancake (smaller volume) is not minimum compared to circular covariance.

getKDEManifoldBandwidths(pts, manif)

Calculate the KDE bandwidths for each dimension independly, as per manifold of each. Return vector of all dimension bandwidths.

getKernelLeaf(mt, i)
getKernelLeaf(mt, i, permuted)

Return leaf kernel associated with input data element i (i.e. permuted=true). Else when set to permuted=false return the sorted leaf_kernel i (different from unsorted input data number).

getKernelLeafAsTreeKer(mtr, idx)
getKernelLeafAsTreeKer(mtr, idx, permuted)

Return leaf kernels as tree kernel types, using regular [1..N] indexing].

Notes:

• use permute=true (default) for sorted index retrieval.

getKernelTree(mtr, currIdx)
getKernelTree(mtr, currIdx, permuted)
getKernelTree(mtr, currIdx, permuted, cov_continuation)

Return kernel from tree by binary tree index, and convert leaf kernels to tree kernel types if necessary.

See also: getKernelLeafAsTreeKer

getManifoldPartial(M, partial; ...)
getManifoldPartial(M, partial, repr; ...)
getManifoldPartial(M, partial, repr, offset; doError)

A so-called full dimension manifold can possibly be reduced to smaller partial manifolds over some of the dimensions, returning a new programatically generated <:AbstractManifold. This funciton can optionally also reduce a point representation for the desired partial dimensions too.

Example

using Manifolds
using ApproxManifoldProducts

# a familiar manifold of translation and rotation in 2D
M = SpecialEuclidean(2)


# make a new partial of only the translation components
M_x, _ = getManifoldPartial(M,[1;])
# returns a new TranslationGroup(1) corresponding to just x dimension

# representation is semidirect product of translation and rotation matrix
u0 = ArrayPartition([0.0;0],[1 0; 0 1.0])

# known coordinates are [x,y,θ], eg
#   vee(M,identity_element(M,u0),log(M,identity_element(M,u0),u0))
#   [0;0;0] in this example

# make another new partial of only the rotation component
M_rot, u_rot = getManifoldPartial(M,[3],u0)
# returns SpecialOrthogonal(2) information

# make another new partial of only  y and θ
M_yθ, u_yθ = getManifoldPartial(M,[2;3],u0)
# returns new manifold information as ProductArray(TranslationGroup(1),SpecialOrthogonal(2))

Notes

• Partial dimensions of interest are defined via a AbstractVector{Int} of coordinate dimensions.
• assumed to follow coordinates as transition step towards more general solution
• assume ProductManifold is the only way to stitch multiple manifolds together
• Still experimental.

DevNotes

• FIXME any semidirect product or action information is lost and naive Product manifold assumptions are currently made.

Related

getManifold, [AbstractManifold], [ProductManifold], [GroupManifold], [ArrayPartition]

isPartial(mkd)

Return true if this ManifoldKernelDensity is a partial.

ker(M, p, q)
ker(M, p, q, sigma)
ker(M, p, q, sigma, distFnc)

Normal kernel used for Hilbert space embeddings.

makePointFromCoords(M, coords)
makePointFromCoords(M, coords, u0)
makePointFromCoords(M, coords, u0, ϵ)
makePointFromCoords(M, coords, u0, ϵ, retraction_method)

Helper function to convert coordinates to a desired on-manifold point.

DevNotes

• FIXME need much better consolidation or even removal of this function entirely.
  • This function makes too strong an assumption of groups

Notes

• u0 is used to identify the data type for a point
• Pass in a different exp if needed.

manifoldLooCrossValidation(pts, bw; own, diffop)

Calculate negative entropy with leave one out (j'th element) cross validation.

Background

From: Silverman, B.: Density Estimation for Statistics and Data Analysis, 1986, p.52

Probability density function p(x), as estimated by kernels

\[hatp_{-j}(x) = 1/(N-1) Σ_{i != j}^N frac{1}{sqrt{2pi}σ } exp{ -frac{(x-μ)^2}{2 σ^2} }\]

and has Cross Validation number as the average log evaluations of leave one out hatp_{-j}(x):

\[CV(p) = 1/N Σ_i^N log hat{p}_{-j}(x_i)\]

This quantity CV is related to an entropy H(p) estimate via:

\[H(p) = -CV(p)\]

manifoldProduct(ff; ...)
manifoldProduct(
    ff,
    mani;
    makeCopy,
    Niter,
    ndims,
    N,
    u0,
    oldPoints,
    addEntropy,
    recordLabels,
    selectedLabels,
    _randU,
    _randN,
    logger,
    bws
)

Approximate the pointwise the product of functionals on manifolds using KernelDensityEstimate.

Notes:

• Always pass full beliefs, for partials use for e.g. partialDimsWorkaround=[1;3;6]
• Can also multiply different partials together

Example

# setup
M = TranslationGroup(3)
N = 75
p = manikde!(M, [randn(3) for _ in 1:N])
q = manikde!(M, [randn(3) .+ 1 for _ in 1:N])

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q])

# direct histogram plot
using Gadfly
plot( x=getPoints(pq)[1,:], y=getPoints(pq)[2,:], Geom.histogram2d )

# TODO, convenient plotting (work in progress...)

mmd(MF, a, b; ...)
mmd(MF, a, b, N; ...)
mmd(MF, a, b, N, M; ...)
mmd(MF, a, b, N, M, threads; bw)

MMD disparity (i.e. 'distance') measure based on Kernel Hilbert Embeddings between two beliefs.

Notes:

• This is a wrapper to the in-place mmd! function.

Related

mmd!, ker

mmd!(MF, val, a, b; ...)
mmd!(MF, val, a, b, N; ...)
mmd!(MF, val, a, b, N, M; ...)
mmd!(MF, val, a, b, N, M, threads; bw)

MMD disparity (i.e. 'distance') measure based on Kernel Hilbert Embeddings between two beliefs.

Notes:

• This is the in-place version (well more in-place than mmd)

DevNotes:

• TODO dont assume equally weighted particles
• TODO profile SIMD vs SLEEF
• TODO optimize memory
• TODO make multithreaded

See also: mmd, ker

normDistAccAt!(ret, idx, x, sigma)
normDistAccAt!(ret, idx, x, sigma, w)

Probability density function p(x), as estimated by kernels

\[hatp_{-j}(x) = 1/(N-1) Σ_{i != j}^N frac{1}{sqrt{2pi}σ } exp{ -frac{(x-μ)^2}{2 σ^2} }\]

parallel_transport_best(M, p, X, d)

EXPERIMENTAL: ApproxManifoldProducts hard coded versions for best parallel transport available (as per development cycle).

Inputs:

• M is a Manifold (must be a Lie group)
• p is the expansion point on manifold
• X is a tangent vector which is to be transported
• d is a tangent vector from a starting point in the direction and distance to transport

Notes

• Default is transport without curvature estimate provided by upstream Manifolds.jl

See also: `paralleltransportcurvature2ndlie', Jr, Manifolds.parallel_transport_direction, Manifolds.parallel_transport_to, Manifolds.parallel_transport_along

productbelief(
    denspts,
    manifold,
    dens,
    N;
    asPartial,
    dbg,
    logger
)

Take product of dens (including optional partials beliefs) as proposals to be multiplied together.

Notes

• Return points of full dimension, even if only partial dimensions in proposals.
  • 'Other' dimensions left unchanged from incoming denspts
• d dimensional product approximation
• Incorporate ApproxManifoldProducts to process variables in individual batches.

DevNotes

• TODO Consolidate with AMP.manifoldProduct, especially concerning partials.

sampleProductSeqGibbsBTLabel(M, proposals; ...)
sampleProductSeqGibbsBTLabel(M, proposals, MC; ...)
sampleProductSeqGibbsBTLabel(
    M,
    proposals,
    MC,
    label_pools;
    ...
)
sampleProductSeqGibbsBTLabel(
    M,
    proposals,
    MC,
    label_pools,
    labels_sampled;
    MAX_RECURSE_DEPTH
)

Notes:

• Advise, 2<=MC to ensure multiscale works during decent transitions (TBD obsolete requirement)
• To force sequential Gibbs on leaves only, use: label_pools = [[(length(getPoints(prop))+1):(2*length(getPoints(prop)));] for prop in proposals]
• Taken from: Sudderth, E.B., Ihler, A.T., Isard, M., Freeman, W.T. and Willsky, A.S., 2010. Nonparametric belief propagation. Communications of the ACM, 53(10), pp.95-103.

splitPointsEigen(M, r_PP; ...)
splitPointsEigen(M, r_PP, weights; kernel, kernel_bw)

Give vector of manifold points and split along largest covariance (i.e. major direction)

DeVNotes:

• FIXME: upgrade to Manopt version
  • https://github.com/JuliaRobotics/ApproxManifoldProducts.jl/issues/277
• FIXME use manifold mean and cov calculation instead

updateProductSample(
    dest,
    proposals,
    manifolds,
    smplIdx,
    labels
)

A clunky repeat calculation of one product kernel.

evaluate(mt, pt)
evaluate(mt, pt, LOO)
evaluate(mt, pt, LOO, force_kbw)

Evaluate the belief density for a given Manellic tree.

DevNotes:

• Computational Geometry
  • use geometric computing for faster evaluation
• Dual tree evaluations
  • Holmes, M.P., Gray, A.G. and Isbell Jr, C.L., 2010. Fast kernel conditional density estimation: A dual-tree Monte Carlo approach. Computational statistics & data analysis, 54(7), pp.1707-1718.
  • Curtin, R., March, W., Ram, P., Anderson, D., Gray, A. and Isbell, C., 2013, May. Tree-independent dual-tree algorithms. In International Conference on Machine Learning (pp. 1435-1443). PMLR.
• Fast kernels
• Parallel transport shortcuts?

getPoints(x)
getPoints(x, )

Return underlying points used to construct the ManifoldKernelDensity.

Notes

• Return type is ::Vector{P} where P represents a Manifold point type (e.g. group element or coordinates).
• Second argument controls whether partial dimensions only should be returned (=true default).

DevNotes

• Currently converts down to manifold from matrix of coordinates (legacy), to be deprecated TODO