Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

differentiable kdes #77

Closed
wants to merge 6 commits into from
Closed

differentiable kdes #77

wants to merge 6 commits into from

Conversation

gszep
Copy link

@gszep gszep commented Oct 11, 2019

I would like the kde method to be used as part of a machine learning pipline. Therefore we need all methods to be extended to TrackedArray types. The plan would be to extend univariate.jl first. Which can be done by by-passing the use of rfft,irfft from FFTW.jl (which is not compatible with TrackedArrays) and use the conv layer from Flux.jl

closes #76

@gszep
Copy link
Author

gszep commented Oct 11, 2019

the following now works:

using Flux,KernelDensity
data = param(randn(100))
xrange = -3:0.1:3
kde(data,xrange,bandwidth=0.1)

return dropdims(result, dims = (1+N,2+N))
end

# patches for TrackedReal and Vector{TrackedReal}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

all of these methods below are type piracy, so probably better to make a PR to Tracker.jl

Copy link
Author

@gszep gszep Oct 19, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see.. suppose these methods are added to Tracker.jl how would you rewrite the conv methods in univariate.jl and bivariate.jl so that they are compatible?

Copy link
Author

@gszep gszep left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the main compatibility issue here is the conv method which calls FFTW which does not accept tracked variables.. open to suggestions on how to patch this

Comment on lines +22 to +40
"""n-dimensional convolution"""
function conv(x::AbstractArray{T,N}, w::AbstractArray{T,N}) where {T,N}
wdim = Int.(ceil.((size(w).-1)./2))
padding = Iterators.flatten([ (wdim[i],wdim[i]) for i=1:length(wdim) ]) |> collect

dims = DenseConvDims((size(x)...,1,1),(size(w)...,1,1); padding=padding )
result = Tracker.conv( reshape(x,(size(x)...,1,1)), reshape(w,(size(w)...,1,1)), dims)
return dropdims(result, dims = (1+N,2+N))
end

# patches for TrackedReal and Vector{TrackedReal}
conv(x::AbstractArray{TrackedReal{T},N}, w::AbstractArray) where {T,N} = conv(Tracker.collect(x),w)
conv(x::AbstractArray, w::AbstractArray{TrackedReal{T},N}) where {T,N} = conv(x,Tracker.collect(w))
conv(x::AbstractArray{TrackedReal{T},N}, w::AbstractArray{TrackedReal{T},N}) where {T,N} = conv(Tracker.collect(x),Tracker.collect(w))

round(::Type{R}, t::TrackedReal) where {R<:Real} = round(R, t.data)
round(t::TrackedReal, mode::RoundingMode) = round(t.data, mode)
Integer(x::TrackedReal) = Integer(x.data)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

bypasses the need for FFTW which is not compatible with tracked types (JuliaMath/AbstractFFTs.jl#35). Tracker.conv throws outofMemoryError() for large kernels w which is why tests are failing

density::Matrix{Float64}
density::AbstractMatrix{}
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

relax type casting for future compatibility.. will make a PR to Tracker.jl soon

grid = zeros(Float64, nx, ny)
grid = zeros(eltype(xdata),nx,ny)
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

relaxing types

Comment on lines 122 to 127
function conv(k::UnivariateKDE, dist::UnivariateDistribution)
# Transform to Fourier basis
K = length(k.density)
ft = rfft(k.density)

# Convolve fft with characteristic function of kernel
# empirical cf
# = \sum_{n=1}^N e^{i*t*X_n} / N
# = \sum_{k=0}^K e^{i*t*(a+k*s)} N_k / N
# = e^{i*t*a} \sum_{k=0}^K e^{-2pi*i*k*(-t*s*K/2pi)/K} N_k / N
# = A * fft(N_k/N)[-t*s*K/2pi + 1]
c = -twoπ/(step(k.x)*K)
for j = 0:length(ft)-1
ft[j+1] *= cf(dist,j*c)
end

dens = irfft(ft, K)
# fix rounding error.
for i = 1:K
dens[i] = max(0.0,dens[i])
end

# Invert the Fourier transform to get the KDE
UnivariateKDE(k.x, dens)
half_grid = range(step(k.x),5*std(dist),step=step(k.x))
grid = [-reverse(half_grid);0;half_grid]
density = conv(k.density, pdf.(dist,grid)) * step(k.x)
UnivariateKDE(k.x, density)
end
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what would your recommendation be for this? given that the methods in KernelDensity.jl are type piracy?

@gszep
Copy link
Author

gszep commented Oct 20, 2019

closing this request since all changes have been moved into patches https://github.com/gszep/FluxContinuation/blob/master/patches/KernelDensity.jl

@tpapp
Copy link
Collaborator

tpapp commented Oct 21, 2019

How are those patches meant to be applied?

@gszep
Copy link
Author

gszep commented Oct 21, 2019

download that file and instead of using KernelDensity write include("pathtopatch/KernelDensity.jl")

@gszep
Copy link
Author

gszep commented Oct 21, 2019

note that currently the Fourier transform for Tracked types is redirected to a naivedft which unfortunately is N^2. this makes the patch unusable for more than ~4000 tracked data inputs. I'm working on a solution to use the NlogN fft algorithms

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

extend methods for TrackedArray and TrackedReal
3 participants