GaussianMixtureAlignment.jl

julia> using GaussianMixtureAlignment; GMA = GaussianMixtureAlignment;

julia> # These are very simple point sets that can be perfectly aligned

julia> xpts = [[0.,0.,0.], [3.,0.,0.,], [0.,4.,0.]];

julia> ypts = [[1.,1.,1.], [1.,-2.,1.], [1.,1.,-3.]];

julia> σ = ϕ = 1.;

julia> gmmx = IsotropicGMM([IsotropicGaussian(x, σ, ϕ) for x in xpts])
IsotropicGMM{3, Float64} with 3 IsotropicGaussian{3, Float64} distributions.


julia> gmmy = IsotropicGMM([IsotropicGaussian(y, σ, ϕ) for y in ypts])
IsotropicGMM{3, Float64} with 3 IsotropicGaussian{3, Float64} distributions.


julia> # Compute the overlap between the two GMMs

julia> overlap(gmmx, gmmy)
1.1908057504684806

julia> res = gogma_align(gmmx, gmmy; nextblockfun=GMA.randomblock, maxsplits=10000);

julia> # upper and lower bounds of alignmnet objective function at search termination

julia> res.upperbound, res.lowerbound
(-3.251290635173653, -5.528700740172087)

julia> # rotation component of the best transformation

julia> res.tform.linear
3×3 AngleAxis{Float64} with indices SOneTo(3)×SOneTo(3)(2.0944, -0.57735, 0.57735, -0.57735):
 -1.38948e-11  -2.5226e-11    1.0
 -1.0           9.12343e-11  -1.38948e-11
 -9.12343e-11  -1.0          -2.52261e-11

julia> # translation component of the best transformation

julia> res.tform.translation
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 1.0000000000115654
 1.000000000557716
 0.9999999997059433

julia> # Compute the overlap between the GMMs after alignment (equal to res.upperbound)

julia> overlap(res.tform(gmmx), gmmy)
3.251290635173653

julia> using GLMakie

julia> # Draw the unaligned GMMs

julia> gmmdisplay(gmmx,gmmy)

julia> # Draw the aligned GMMs

julia> gmmdisplay(res.tform(gmmx), gmmy)

julia> # four points defining a tetrahedron

julia> tetrahedron = [
            [0.,0.,1.],
            [sqrt(8/9), 0., -1/3],
            [-sqrt(2/9),sqrt(2/3),-1/3],
            [-sqrt(2/9),-sqrt(2/3),-1/3]
        ];

julia> # each feature type is given its own IsotropicGMM

julia> chargedGMM = IsotropicGMM([IsotropicGaussian(tetrahedron[1], 1.0, 1.0)]);

julia> stericGMM = IsotropicGMM([IsotropicGaussian(p, 0.5, 1.0) for p in tetrahedron]);

julia> mgmmx = IsotropicMultiGMM(Dict(
           :positive => chargedGMM,
           :steric => stericGMM
       ))
IsotropicMultiGMM{3, Float64, Symbol} with 2 labeled IsotropicGMM{3, Float64} models made up of a total of 5 IsotropicGMM{3, Float64} distributions.

julia> mgmmy = IsotropicMultiGMM(Dict(
           :negative => chargedGMM,
           :steric => stericGMM
       ))
IsotropicMultiGMM{3, Float64, Symbol} with 2 labeled IsotropicGMM{3, Float64} models made up of a total of 5 IsotropicGMM{3, Float64} distributions.

julia> # define interactions between each type of feature

julia> interactions = Dict(
           :positive => Dict(
               :positive => -1.0,
               :negative => 1.0,
           ),
           :negative => Dict(
               :positive => 1.0,
               :negative => -1.0,
           ),
           :steric => Dict(
               :steric => -1.0,
           ),
       );

julia> # apply a random rigid transformation to one of the models

julia> using CoordinateTransformations; randtform = AffineMap((rand(), rand(), rand(), rand(), rand(), rand()));

julia> mgmmx = randtform(mgmmx);

julia> # compute alignment between the two MultiIsotropicGMMs

julia> res = gogma_align(mgmmx, mgmmy; interactions=interactions, nextblockfun=GMA.randomblock, maxsplits=10000);

julia> # Draw the aligned GMMs. Note that the "steric" features (green) repulse one another.

julia> gmmdisplay(res.tform(mgmmx), mgmmy)