multidim manellic on SE2

src/services/ManellicTree.jl

@@ -248,7 +248,7 @@ function Base.convert(
 ) where {P,T,M,iM}
   #
   _matType(::Type{Distributions.PDMats.PDMat{_F,_M}}) where {_F,_M} = _M
-  μ = P(src.μ)
+  μ = convert(P,src.μ) # P(src.μ)
   p = MvNormal(_matType(M)(cov(src.p)))
   sqrt_iΣ = iM(src.sqrt_iΣ)
   MvNormalKernel{P,T,M,iM}(μ, p, sqrt_iΣ, src.weight)

src/services/ManifoldKernelDensity.jl

@@ -108,6 +108,7 @@ function manikde!_manellic(
   M::AbstractManifold,
   pts::AbstractVector;
   bw=diagm(ones(manifold_dimension(M))),
+  algo = Optim.NelderMead()
 )
   #
 
@@ -138,9 +139,9 @@ function manikde!_manellic(
     res = Optim.optimize(
       _cost, 
       diag(bw), # FIXME Optim API issue, if using bw::matrix then steps not PDMat (NelderMead) 
-      Optim.NelderMead()
+      algo
     )
-    diagm(Optim.minimizer(res))
+    diagm(abs.(Optim.minimizer(res)))
   end
 
   # reuse (heavy lift parts of) earlier tree build

test/manellic/testManellicTree.jl

@@ -7,6 +7,7 @@ using LinearAlgebra
 using StaticArrays
 using TensorCast
 using Manifolds
+import Rotations as Rot_
 using Distributions
 import ApproxManifoldProducts: ManellicTree, eigenCoords, splitPointsEigen
 
@@ -36,8 +37,8 @@ function testEigenCoords(
 
   # spot check
   @show _ax_ERR = log_lie(SpecialOrthogonal(2), (r_R_ax_')*r_R_ax)[1,2]
-  @show testval = isapprox(0, _ax_ERR; atol = 5/length(ax_CC))
-  @assert testval "Spot check failed on eigen split of manifold points, the estimated point rotation matrix did not match construction."
+  @show testval = isapprox(0, _ax_ERR; atol = 6/length(ax_CC))
+  @assert testval "Spot check failed on eigen split of manifold points, the estimated point rotation matrix did not match construction. length(ax_CC)=$(length(ax_CC))"
 
   r_CC, r_R_ax_, pidx, r_CV
 end
@@ -807,7 +808,7 @@ end
 end
 
 
-@testset "Multidimensional LOOCV bandwidth optimization" begin
+@testset "Multidimensional LOOCV bandwidth optimization, TranslationGroup(2)" begin
 ##
 
 M = TranslationGroup(2)
@@ -829,7 +830,7 @@ end
 
 @test res.ls_success
 
-@show best_cov = Optim.minimizer(res)
+@show best_cov = abs.(Optim.minimizer(res))
 
 @test isapprox([0.5; 0.5], best_cov; atol=0.3)
 
@@ -843,6 +844,44 @@ mkd = ApproxManifoldProducts.manikde!_manellic(M,pts)
 end
 
 
+
+@testset "Multidimensional LOOCV bandwidth optimization, SpecialEuclidean(2)" begin
+##
+
+M = SpecialEuclidean(2)
+pts = [ArrayPartition(randn(2),Rot_.RotMatrix{2}(0.1*randn()).mat) for _ in 1:64]
+
+bw = [1.0; 1.0; 0.3]
+mtree = ApproxManifoldProducts.buildTree_Manellic!(M, pts; kernel_bw=bw,kernel=AMP.MvNormalKernel)
+
+cost4(σ) = begin
+  AMP.entropy(mtree, diagm(σ.^2))
+end
+
+# and optimize with "update" kernel bandwith cost
+@time res = Optim.optimize(
+  cost4, 
+  bw, 
+  Optim.NelderMead()
+);
+
+@test res.ls_success
+
+@show best_cov = abs.(Optim.minimizer(res))
+
+@test isapprox([0.6; 0.6; 0.06], best_cov; atol=0.35)
+
+
+mkd = ApproxManifoldProducts.manikde!_manellic(M,pts)
+
+@test isapprox([0.6 0 0; 0 0.6 0; 0 0 0.06], getBW(mkd)[1]; atol=0.35)
+
+
+##
+end
+
+
+
 ##
 # # using GLMakie
 

test/testMarginalProducts.jl

@@ -47,8 +47,8 @@ P12 = manifoldProduct([P1;P2], recordLabels=true, selectedLabels=sl, addEntropy=
 # @test isapprox( mean(P12)[1], 0, atol=1 )
 # @test isapprox( mean(P12)[2], 0, atol=1 )
 
-(x->println()).(1:5)
-@show sl;
+
+# @show sl;
 
 P12
 
@@ -99,7 +99,7 @@ P12_ = manifoldProduct([P1;P2_], recordLabels=true, selectedLabels=sl, addEntrop
 @test isapprox( mean(P12_)[1], 0, atol=1 )
 @test isapprox( mean(P12_)[2], 0, atol=1 )
 
-(x->println()).(1:5)
+
 # @show sl
 
 P12_
@@ -156,7 +156,7 @@ P123_ = manifoldProduct([P1;P2_;P3_], recordLabels=true, selectedLabels=sl, addE
 @test isapprox( mean(P123_)[1], 0, atol=1 )
 @test isapprox( mean(P123_)[2], 0, atol=1 )
 
-(x->println()).(1:5)
+
 # @show sl
 
 P123_
@@ -207,7 +207,7 @@ P = manifoldProduct([P2;P1;P3], recordLabels=true, selectedLabels=sl, addEntropy
 
 @test !isPartial(P)
 
-(x->println()).(1:5)
+
 # @show sl;
 P
 
@@ -267,7 +267,7 @@ P_ = manifoldProduct([P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=f
 
 @test !isPartial(P_)
 
-(x->println()).(1:5)
+
 # @show sl
 
 ##
@@ -329,7 +329,7 @@ P45__ = manifoldProduct([P4;P4_;P5;P5_], recordLabels=true, selectedLabels=sl, a
 
 @test !isPartial(P45__)
 
-(x->println()).(1:5)
+
 # @show sl;
 
 P45__
@@ -393,7 +393,7 @@ P = manifoldProduct([P2;P1;P3], recordLabels=true, selectedLabels=sl, addEntropy
 
 @test !isPartial(P)
 
-(x->println()).(1:5)
+
 # @show sl;
 P
 
@@ -456,7 +456,7 @@ P = manifoldProduct([P1;P2;P3], recordLabels=true, selectedLabels=sl, addEntropy
 
 @test !isPartial(P)
 
-(x->println()).(1:5)
+
 # @show sl;
 P
 
@@ -517,7 +517,7 @@ P = manifoldProduct([P1;P3], recordLabels=true, selectedLabels=sl, addEntropy=fa
 @test isPartial(P)
 @test P._partial == [1;3]
 
-(x->println()).(1:5)
+
 # @show sl;
 P