Skip to content

A few updates #30

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

Merged
merged 7 commits into from
Jan 29, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
*.jl.cov
*.jl.*.cov
*.jl.mem
Manifest.toml
docs/build/
docs/site/
docs/Manifest.toml
15 changes: 8 additions & 7 deletions .travis.yml
Original file line number Diff line number Diff line change
@@ -1,11 +1,12 @@
language: julia
os:
- linux
- osx
- linux
julia:
- 0.7
- 1.0
- nightly
- 0.7
- 1.0
- 1.1
- nightly
notifications:
email: false
sudo: false
email: false
after_success:
- julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())';
13 changes: 13 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
name = "Loess"
uuid = "4345ca2d-374a-55d4-8d30-97f9976e7612"

[deps]
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[extras]
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Random", "Test"]
121 changes: 61 additions & 60 deletions src/Loess.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
__precompile__()


module Loess

import Distances.euclidean
Expand Down Expand Up @@ -38,9 +35,11 @@ Returns:

"""
function loess(xs::AbstractMatrix{T}, ys::AbstractVector{T};
normalize::Bool=true, span::T=0.75, degree::Int=2) where T <: AbstractFloat
normalize::Bool=true,
span::AbstractFloat=0.75,
degree::Integer=2) where T<:AbstractFloat
if size(xs, 1) != size(ys, 1)
error("Predictor and response arrays must of the same length")
throw(DimensionMismatch("Predictor and response arrays must of the same length"))
end

n, m = size(xs)
Expand All @@ -50,7 +49,7 @@ function loess(xs::AbstractMatrix{T}, ys::AbstractVector{T};
# correctly apply predict to unnormalized data. We should have a normalize
# function that just returns a vector of scaling factors.
if normalize && m > 1
xs = tnormalize!(copy(xs))
xs = tnormalize!(copy(xs))
end

kdtree = KDTree(xs, 0.05 * span)
Expand All @@ -59,7 +58,7 @@ function loess(xs::AbstractMatrix{T}, ys::AbstractVector{T};
# map verticies to their index in the bs coefficient matrix
verts = Dict{Vector{T}, Int}()
for (k, vert) in enumerate(kdtree.verts)
verts[vert] = k
verts[vert] = k
end

# Fit each vertex
Expand All @@ -72,44 +71,46 @@ function loess(xs::AbstractMatrix{T}, ys::AbstractVector{T};
vs = Array{T}(undef, q)
for (vert, k) in verts
# reset perm
for i in 1:n
perm[i] = i
end
for i in 1:n
perm[i] = i
end

# distance to each point
for i in 1:n
ds[i] = euclidean(vec(vert), vec(xs[i,:]))
end
for i in 1:n
ds[i] = euclidean(vec(vert), vec(xs[i,:]))
end

# copy the q nearest points to vert into X
partialsort!(perm, q, by=i -> ds[i])
dmax = maximum([ds[perm[i]] for i = 1:q])
# copy the q nearest points to vert into X
partialsort!(perm, q, by=i -> ds[i])
dmax = maximum([ds[perm[i]] for i = 1:q])

for i in 1:q
for i in 1:q
pi = perm[i]
w = tricubic(ds[pi] / dmax)
us[i,1] = w
for j in 1:m
x = xs[pi, j]
wxl = w
for l in 1:degree
w = tricubic(ds[pi] / dmax)
us[i,1] = w
for j in 1:m
x = xs[pi, j]
wxl = w
for l in 1:degree
wxl *= x
us[i, 1 + (j-1)*degree + l] = wxl # w*x^l
end
end
vs[i] = ys[pi] * w
end
bs[k,:] = us \ vs
us[i, 1 + (j-1)*degree + l] = wxl # w*x^l
end
end
vs[i] = ys[pi] * w
end
bs[k,:] = us \ vs
end

LoessModel{T}(xs, ys, bs, verts, kdtree)
end

function loess(xs::AbstractVector{T}, ys::AbstractVector{T};
normalize::Bool=true, span::T=0.75, degree::Int=2) where T <: AbstractFloat
loess(reshape(xs, (length(xs), 1)), ys, normalize=normalize, span=span, degree=degree)
end
loess(xs::AbstractVector{T}, ys::AbstractVector{T}; kwargs...) where {T<:AbstractFloat} =
loess(reshape(xs, (length(xs), 1)), ys; kwargs...)

function loess(xs::AbstractArray{T,N}, ys::AbstractVector{S}; kwargs...) where {T,N,S}
R = float(promote_type(T, S))
loess(convert(AbstractArray{R,N}, xs), convert(AbstractVector{R}, ys); kwargs...)
end


# Predict response values from a trained loess model and predictor observations.
Expand All @@ -126,7 +127,7 @@ end
# A length n' vector of predicted response values.
#
function predict(model::LoessModel{T}, z::T) where T <: AbstractFloat
predict(model, T[z])
predict(model, T[z])
end


Expand All @@ -136,40 +137,40 @@ function predict(model::LoessModel{T}, zs::AbstractVector{T}) where T <: Abstrac
# in the univariate case, interpret a non-singleton zs as vector of
# ponits, not one point
if m == 1 && length(zs) > 1
return predict(model, reshape(zs, (length(zs), 1)))
return predict(model, reshape(zs, (length(zs), 1)))
end

if length(zs) != m
error("$(m)-dimensional model applied to length $(length(zs)) vector")
error("$(m)-dimensional model applied to length $(length(zs)) vector")
end

adjacent_verts = traverse(model.kdtree, zs)

if m == 1
@assert(length(adjacent_verts) == 2)
z = zs[1]
u = (z - adjacent_verts[1][1]) /
(adjacent_verts[2][1] - adjacent_verts[1][1])

y1 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[1][1]]],:])
y2 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[2][1]]],:])
return (1.0 - u) * y1 + u * y2
@assert(length(adjacent_verts) == 2)
z = zs[1]
u = (z - adjacent_verts[1][1]) /
(adjacent_verts[2][1] - adjacent_verts[1][1])

y1 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[1][1]]],:])
y2 = evalpoly(zs, model.bs[model.verts[[adjacent_verts[2][1]]],:])
return (1.0 - u) * y1 + u * y2
else
error("Multivariate blending not yet implemented")
# TODO:
# 1. Univariate linear interpolation between adjacent verticies.
# 2. Blend these estimates. (I'm not sure how this is done.)
error("Multivariate blending not yet implemented")
# TODO:
# 1. Univariate linear interpolation between adjacent verticies.
# 2. Blend these estimates. (I'm not sure how this is done.)
end
end


function predict(model::LoessModel{T}, zs::AbstractMatrix{T}) where T <: AbstractFloat
ys = Array{T}(undef, size(zs, 1))
for i in 1:size(zs, 1)
# the vec() here is not necessary on 0.5 anymore
ys[i] = predict(model, vec(zs[i,:]))
end
ys
ys = Array{T}(undef, size(zs, 1))
for i in 1:size(zs, 1)
# the vec() here is not necessary on 0.5 anymore
ys[i] = predict(model, vec(zs[i,:]))
end
ys
end

"""
Expand Down Expand Up @@ -201,13 +202,13 @@ function evalpoly(xs, bs)
degree = div(length(bs) - 1, m)
y = bs[1]
for i in 1:m
x = xs[i]
x = xs[i]
xx = x
y += xx * bs[1 + (i-1)*degree + 1]
for l in 2:degree
for l in 2:degree
xx *= x
y += xx * bs[1 + (i-1)*degree + l]
end
y += xx * bs[1 + (i-1)*degree + l]
end
end
y
end
Expand All @@ -230,8 +231,8 @@ function tnormalize!(xs::AbstractMatrix{T}, q::T=0.1) where T <: AbstractFloat
n, m = size(xs)
cut = ceil(Int, (q * n))
for j in 1:m
tmp = sort!(xs[:,j])
xs[:,j] ./= mean(tmp[cut+1:n-cut])
tmp = sort!(xs[:,j])
xs[:,j] ./= mean(tmp[cut+1:n-cut])
end
xs
end
Expand Down
Loading