diff --git a/src/services/ManellicTree.jl b/src/services/ManellicTree.jl index 528dd4f..aa2ba8e 100644 --- a/src/services/ManellicTree.jl +++ b/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) diff --git a/src/services/ManifoldKernelDensity.jl b/src/services/ManifoldKernelDensity.jl index 49d0010..e692be6 100644 --- a/src/services/ManifoldKernelDensity.jl +++ b/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 diff --git a/test/manellic/testManellicTree.jl b/test/manellic/testManellicTree.jl index 9c9c1e7..175f54d 100755 --- a/test/manellic/testManellicTree.jl +++ b/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 diff --git a/test/testMarginalProducts.jl b/test/testMarginalProducts.jl index 9cddfbe..0a42dcf 100644 --- a/test/testMarginalProducts.jl +++ b/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