Merge pull request #79 from JuliaRobotics/21Q2/enh/partials

add partialDimMask to product

.github/workflows/ci.yml

@@ -0,0 +1,50 @@
+name: CI
+on:
+  pull_request:
+  push:
+    branches:
+      - master
+
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    env:
+      JULIA_PKG_SERVER: ""
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.6'
+          - 'nightly'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@latest
+      - run: |
+          git config --global user.name Tester
+          git config --global user.email [email protected]
+      - uses: julia-actions/julia-runtest@latest
+        continue-on-error: ${{ matrix.version == 'nightly' }}
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
+

Project.toml

@@ -2,7 +2,7 @@ name = "KernelDensityEstimate"
 uuid = "2472808a-b354-52ea-a80e-1658a3c6056d"
 keywords = ["NBP", "nonparametric", "kernel density", "functional products", "multiscale", "Gibbs sampling", "KDE"]
 desc = "A modified version of the multiscale Gibbs product and KDE"
-version = "0.5.5"
+version = "0.5.6"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/KDE01.jl

@@ -1,8 +1,8 @@
 
 
 function kde!(points::A,
-              addop=(+,),
-              diffop=(-,) ) where {A <: AbstractArray{Float64,2}}
+              addop::Tuple=(+,),
+              diffop::Tuple=(-,) ) where {A <: AbstractArray{Float64,2}}
   #
   dims = size(points,1)
 
@@ -31,11 +31,11 @@ function kde!(points::Array{Float64,1}, addop=(+,), diffop=(-,) )
   return kde!(reshape(points, 1, length(points)), addop, diffop )
 end
 
-function kde!(points::A,
+function kde!(points::AbstractArray{<:Real,2},
               ks::Array{Float64,1},
               weights::Array{Float64,1},
               addop=(+,),
-              diffop=(-,)  ) where {A <: AbstractArray{Float64,2}}
+              diffop=(-,)  )
   #
   Nd, Np = size(points)
   if (length(ks) == 1)
@@ -61,7 +61,7 @@ end
 
 Construct a BallTreeDensity object using `points` for centers and bandwidth `ks`.
 """
-function kde!(points::A, ks::Array{Float64,1}, addop=(+,), diffop=(-,)) where {A <: AbstractArray{Float64,2}}
+function kde!(points::A, ks::Array{Float64,1}, addop::Tuple=(+,), diffop::Tuple=(-,)) where {A <: AbstractArray{Float64,2}}
 
   Nd, Np = size(points)
   weights = ones(Np)
@@ -75,7 +75,7 @@ function kde!(points::A, ks::Array{Float64,1}, addop=(+,), diffop=(-,)) where {A
   kde!(points, ks, weights, addop, diffop)
 end
 
-function kde!(points::Array{Float64,1}, ks::Array{Float64,1}, addop=(+,), diffop=(-,))
+function kde!(points::Array{Float64,1}, ks::Array{Float64,1}, addop::Tuple=(+,), diffop::Tuple=(-,))
   Np = length(points)
   pts = zeros(1,Np)
   pts[1,:] = points

src/MSGibbs01.jl

@@ -29,6 +29,7 @@ mutable struct GbGlb
     # labelsChoosen[sample][densityId][level]
     labelsChoosen::Dict{Int,Dict{Int,Dict{Int,Int}}}
     recordChoosen::Bool
+    partialDimMask::Vector{BitArray{1}}
 end
 
 function makeEmptyGbGlb(;recordChoosen::Bool=false)
@@ -48,11 +49,14 @@ function makeEmptyGbGlb(;recordChoosen::Bool=false)
                 ones(Int,1,0),
                 ones(Int,1,0),
                 zeros(Int,0),
-                0, 0,
+                0, 
+                0,
                 Float64[0.0;], Float64[0.0;],
-                zeros(0), zeros(0),
+                zeros(0), 
+                zeros(0),
                 Dict{Int,Dict{Int,Dict{Int,Int}}}(),
-                recordChoosen  )
+                recordChoosen,
+                Vector{BitVector}()  )
 end
 
 
@@ -70,6 +74,8 @@ function printGlbs(g::GbGlb, tag::String="")
   @show g.levelListNew
   @show round.(g.newPoints, digits=4)
   @show g.newIndices
+  @show g.partialDimMask
+  nothing
 end
 
 
@@ -185,10 +191,10 @@ function gaussianProductMeanCov!( glb::GbGlb,
   else
     # on manifold development
     @inbounds @fastmath @simd for j in 1:glb.Ndens
-      if (j!=skip)
+      if (j!=skip) && glb.partialDimMask[j][dim]
         # who populated glb.particles?
         glb.calclambdas[j] = 1.0/glb.variance[dim,j]
-        glb.calcmu[j] = glb.particles[dim,j] #[dim+glb.Ndim*(j-1)]
+        glb.calcmu[j] = glb.particles[dim,j]
       else
         # adding zeros does not influence the result because `lambda_i = 0`
         glb.calclambdas[j] = 0.0
@@ -495,7 +501,9 @@ function gibbs1(Ndens::Int, trees::Array{BallTreeDensity,1},
                 getMu=(getEuclidMu,),
                 getLambda=(getEuclidLambda,),
                 glbs = makeEmptyGbGlb(),
-                addEntropy::Bool=true  )
+                addEntropy::Bool=true,
+                ndims::Int=maximum(Ndim.(trees)),
+                partialDimMask::AbstractVector{<:BitVector} = [ones(Int,ndims) .== 1 for i in 1:Ndens]  )
   #
 
   # number densities in the product
@@ -510,7 +518,8 @@ function gibbs1(Ndens::Int, trees::Array{BallTreeDensity,1},
   glbs.randU = randU
   glbs.randN = randN
   # the dimension of the incoming densities
-  glbs.Ndim = trees[1].bt.dims
+  glbs.Ndim = maximum(Ndim.(trees))
+  glbs.partialDimMask = partialDimMask
 
   maxNp = 0                         # largest # of particles we deal with
   for tree in trees
@@ -612,27 +621,29 @@ function prodAppxMSGibbsS(npd0::BallTreeDensity,
                           getMu::T3=(getEuclidMu,),
                           getLambda::T4=(getEuclidLambda,),
                           glbs = makeEmptyGbGlb(),
-                          addEntropy::Bool=true  ) where {T1<:Tuple,T2<:Tuple,T3<:Tuple,T4<:Tuple}
+                          addEntropy::Bool=true,
+                          ndims::Int=maximum(Ndim.(trees)),
+                          partialDimMask::AbstractVector{BitVector} = [ones(Int,ndims) .== 1 for i in 1:length(trees)]  
+                        ) where {T1<:Tuple,T2<:Tuple,T3<:Tuple,T4<:Tuple}
   #
   #
 
-
   # See  Ihler,Sudderth,Freeman,&Willsky, "Efficient multiscale sampling from products
   #         of Gaussian mixtures", in Proc. Neural Information Processing Systems 2003
   Ndens = length(trees)              # of densities
-  Ndim  = trees[1].bt.dims           # of dimensions
+  # ndims  = trees[1].bt.dims           # of dimensions
   Np    = Npts(npd0)                # of points to sample
 
   # prepare stack manifold add and diff operations functions (manifolds must match dimension)
-  addopT = length(addop)!=Ndim ? ([ (addop[1]) for i in 1:Ndim]...,) : addop
-  diffopT = length(diffop)!=Ndim ? ([ (diffop[1]) for i in 1:Ndim]...,) : diffop
-  getMuT = length(getMu)!=Ndim ? ([ getMu[1] for i in 1:Ndim]...,) : getMu
-  getLambdaT = length(getLambda)!=Ndim ? ([ getLambda[1] for i in 1:Ndim]...,) : getLambda
+  addopT = length(addop)!=ndims ? ([ (addop[1]) for i in 1:ndims]...,) : addop
+  diffopT = length(diffop)!=ndims ? ([ (diffop[1]) for i in 1:ndims]...,) : diffop
+  getMuT = length(getMu)!=ndims ? ([ getMu[1] for i in 1:ndims]...,) : getMu
+  getLambdaT = length(getLambda)!=ndims ? ([ getLambda[1] for i in 1:ndims]...,) : getLambda
 
   # skipping analytic functions for now TODO ??
   UseAn = false
-  #??pointsM = zeros(Ndim, Np)
-  points = zeros(Ndim*Np)
+  #??pointsM = zeros(ndims, Np)
+  points = zeros(ndims*Np)
   #??plhs[1] = mxCreateNumericMatrix(Ndens, Np, mxUINT32_CLASS, mxREAL);
   indices=ones(Int,Ndens, Np)
   maxNp = Np
@@ -648,15 +659,19 @@ function prodAppxMSGibbsS(npd0::BallTreeDensity,
   # Generate enough random numbers to get us through the rest of this
   if true
     randU = rand(Int(Np*Ndens*(Niter+2)*Nlevels))
-    randN = randn(Int(Ndim*Np*(Nlevels+1)))
+    randN = randn(Int(ndims*Np*(Nlevels+1)))
   else
     # FIXME using DelimitedFiles
     randU = vec(readdlm("randU.csv"))
     randN = vec(readdlm("randN.csv"))
   end
 
-  gibbs1(Ndens, trees, Np, Niter, points, indices, randU, randN, addop=addopT, diffop=diffopT, getMu=getMuT, getLambda=getLambdaT, glbs=glbs, addEntropy=addEntropy );
-  return reshape(points, Ndim, Np), indices
+  gibbs1( Ndens, trees, Np, Niter, points, indices, 
+          randU, randN, addop=addopT, diffop=diffopT, ndims=ndims,
+          getMu=getMuT, getLambda=getLambdaT, glbs=glbs, 
+          addEntropy=addEntropy, partialDimMask=partialDimMask );
+  #
+  return reshape(points, ndims, Np), indices
 end
 
 
@@ -673,7 +688,7 @@ function *( trees::Vector{BallTreeDensity};
   end
 
   numpts = round(Int, Statistics.mean(Npts.(trees)))
-  d = Ndim(trees[1])
+  d = maximum(Ndim.(trees))
   for p in trees
     d != Ndim(p) ? error("kdes must have same dimension") : nothing
   end

test/runtests.jl

@@ -253,3 +253,6 @@ pp = convert(BallTreeDensity, ps)
 @test norm(getBW(pp)-getBW(p)) < 1e-4
 
 end
+
+
+include("testPartialProd.jl")

test/testPartialProd.jl

@@ -0,0 +1,61 @@
+
+
+using KernelDensityEstimate
+
+
+##
+
+@testset "test partial product" begin
+
+##
+
+pts1 = rand(2,100) .+ 10.0
+mask1 = 1 .== ones(2)
+mask1[2] = false
+
+pts2 = rand(2,100)
+mask2 = 1 .== ones(2)
+
+pts3 = rand(2,100) .- 10.0
+mask3 = 1 .== ones(2)
+mask3[1] = false
+
+
+##
+
+P1 = kde!(pts1)
+P2 = kde!(pts2)
+P3 = kde!(pts3)
+
+bw1 = getBW(P1)[:,1]
+bw2 = getBW(P2)[:,1]
+bw3 = getBW(P3)[:,1]
+
+pts1[2,:] .= 9999999.0
+pts3[1,:] .= 9999999.0
+
+P1 = kde!(pts1, bw1)
+P3 = kde!(pts3, bw3)
+
+
+mask = [mask1, mask2, mask3]
+
+##
+
+dummy = kde!(rand(2,100));
+
+pGM, = prodAppxMSGibbsS(dummy, [P1;P2;P3], nothing, nothing, partialDimMask=mask )
+
+# P = kde!(pGM)
+
+@test 80 < sum(0 .< pGM[1,:] .< 10 )
+
+@test 80 < sum(-10 .< pGM[2,:] .< 0)
+
+
+##
+
+end
+
+
+#