Skip to content

Commit

Permalink
Support Julia 0.7, drop 0.6
Browse files Browse the repository at this point in the history
  • Loading branch information
helgee authored and mauro3 committed Jul 16, 2018
1 parent e83f162 commit 20ee02a
Show file tree
Hide file tree
Showing 6 changed files with 60 additions and 65 deletions.
3 changes: 1 addition & 2 deletions .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,7 @@ os:
- linux
- osx
julia:
- 0.5
- 0.6
- 0.7
- nightly
matrix:
allow_failures:
Expand Down
6 changes: 3 additions & 3 deletions REQUIRE
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
julia 0.5
StatsBase 0.15
Reexport 0.0.3
julia 0.7-beta
StatsBase 0.24.0
Reexport 0.1.0
15 changes: 9 additions & 6 deletions appveyor.yml
Original file line number Diff line number Diff line change
@@ -1,11 +1,14 @@
environment:
matrix:
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
#- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
#- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.7/julia-0.7-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.7/julia-0.7-latest-win64.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"

matrix:
allow_failures:
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"

branches:
only:
Expand Down
27 changes: 14 additions & 13 deletions src/SmoothingSplines.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ module SmoothingSplines

import StatsBase: fit!, fit, RegressionModel, rle, ordinalrank, mean
using Reexport
using LinearAlgebra

export SmoothingSpline

Expand All @@ -13,7 +14,7 @@ LAPACKFloat = Union{Float32,Float64}

include("matrices.jl")

type SmoothingSpline{T<:LAPACKFloat} <: RegressionModel
mutable struct SmoothingSpline{T<:LAPACKFloat} <: RegressionModel
Xorig::Vector{T} # original grid points, but sorted
Yorig::Vector{T} # original values, sorted according to x
Xrank::Vector{Int}
Expand All @@ -27,26 +28,26 @@ type SmoothingSpline{T<:LAPACKFloat} <: RegressionModel
λ::T
end

function fit{T<:LAPACKFloat}(::Type{SmoothingSpline}, X::AbstractVector{T}, Y::AbstractVector{T}, λ::T, wts::AbstractVector{T}=ones(Y))
function fit(::Type{SmoothingSpline}, X::AbstractVector{T}, Y::AbstractVector{T}, λ::T, wts::AbstractVector{T}=fill(1.0, length(Y))) where T <: LAPACKFloat
Xrank = ordinalrank(X) # maybe speed this up when already sorted
Xperm = sortperm(X)
Xorig = X[Xperm]
Yorig = Y[Xperm]

Xdesign, Xcount = rle(Xorig)
ws = zeros(Xdesign)
Ydesign = zeros(Xdesign)
ws = zero(Xdesign)
Ydesign = zero(Xdesign)
running_rle_mean!(Ydesign, ws, Yorig, Xcount, wts[Xperm])

RpαQtQ = QtQpR(diff(Xdesign), λ, ws)
pbtrf!('U', 2, RpαQtQ)

spl = SmoothingSpline{T}(Xorig, Yorig, Xrank, Xdesign, Xcount, Ydesign, ws, RpαQtQ,
zeros(Xdesign), zeros(T,length(Xdesign)-2), λ)
zero(Xdesign), fill(zero(T),length(Xdesign)-2), λ)
fit!(spl)
end

function fit!{T<:LAPACKFloat}(spl::SmoothingSpline{T})
function fit!(spl::SmoothingSpline{T}) where T<:LAPACKFloat
Y = spl.Ydesign
ws = spl.weights
g = spl.g
Expand All @@ -65,17 +66,17 @@ function fit!{T<:LAPACKFloat}(spl::SmoothingSpline{T})
spl
end

function fit!{T<:LAPACKFloat}(spl::SmoothingSpline{T}, Y::AbstractVector{T})
function fit!(spl::SmoothingSpline{T}, Y::AbstractVector{T}) where T<:LAPACKFloat
spl.Y = Y[spl.idx]
fit!(spl)
end


function predict{T<:LAPACKFloat}(spl::SmoothingSpline{T})
function predict(spl::SmoothingSpline{T}) where T<:LAPACKFloat
# need to convert back from RLE encoding
Xcount = spl.Xcount
curridx = 1::Int
g = zeros(length(spl.Yorig))
g = fill(0.0, length(spl.Yorig))
@inbounds for i=1:length(Xcount)
@inbounds for j=1:Xcount[i]
g[curridx] = spl.g[i]
Expand All @@ -85,7 +86,7 @@ function predict{T<:LAPACKFloat}(spl::SmoothingSpline{T})
g[spl.Xrank]
end

function predict{T<:SmoothingSplines.LAPACKFloat}(spl::SmoothingSpline{T}, x::T)
function predict(spl::SmoothingSpline{T}, x::T) where T<:LAPACKFloat
n = length(spl.Xdesign)
idxl = searchsortedlast(spl.Xdesign, x)
idxr = idxl + 1
Expand Down Expand Up @@ -119,8 +120,8 @@ function predict{T<:SmoothingSplines.LAPACKFloat}(spl::SmoothingSpline{T}, x::T)
val
end

function predict{T<:SmoothingSplines.LAPACKFloat}(spl::SmoothingSpline{T}, xs::AbstractVector{T})
g = zeros(xs)
function predict(spl::SmoothingSpline{T}, xs::AbstractVector{T}) where T<:SmoothingSplines.LAPACKFloat
g = zero(xs)
for (i,x) in enumerate(xs)
# can be made more efficient as in StatsBase ECDF code
g[i] = predict(spl, x)
Expand All @@ -129,7 +130,7 @@ function predict{T<:SmoothingSplines.LAPACKFloat}(spl::SmoothingSpline{T}, xs::A
end

# update g and w in place
function running_rle_mean!{T<:Real}(g::AbstractVector{T}, w::AbstractVector{T}, Y::AbstractVector{T}, rlecount::AbstractVector{Int}, ws::AbstractVector{T})
function running_rle_mean!(g::AbstractVector{T}, w::AbstractVector{T}, Y::AbstractVector{T}, rlecount::AbstractVector{Int}, ws::AbstractVector{T}) where T<:Real
length(g) == length(rlecount) || throw(DimensionMismatch())
length(Y) == length(ws) || throw(DimensionMismatch())
curridx = 1::Int
Expand Down
53 changes: 21 additions & 32 deletions src/matrices.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import Base.LinAlg.LAPACK: chkstride1, chkuplo, BlasInt, liblapack
import LinearAlgebra.LAPACK: chkstride1, chkuplo, BlasInt, liblapack
import LinearAlgebra.BLAS: @blasfunc


immutable ReinschQ{T} <: AbstractMatrix{T}
struct ReinschQ{T} <: AbstractMatrix{T}
h::Vector{T}
end

Expand All @@ -10,7 +10,7 @@ function Base.size(Q::ReinschQ)
(n+1, n-1)
end

function Base.getindex{T}(Q::ReinschQ{T},i::Int,j::Int)
function Base.getindex(Q::ReinschQ{T},i::Int,j::Int) where T
h = Q.h
if (i == j)
1/h[i]
Expand All @@ -23,7 +23,7 @@ function Base.getindex{T}(Q::ReinschQ{T},i::Int,j::Int)
end
end

function Base.LinAlg.At_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVector)
function LinearAlgebra.At_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVector)
n = length(out)
h = Q.h
n == size(Q,2) || throw(DimensionMismatch())
Expand All @@ -37,7 +37,7 @@ function Base.LinAlg.At_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVect
out
end

function Base.LinAlg.A_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVector)
function LinearAlgebra.A_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVector)
n = length(out)
n == size(Q,1) || throw(DimensionMismatch())
length(g) == size(Q,2) || throw(DimensionMismatch())
Expand All @@ -51,7 +51,7 @@ function Base.LinAlg.A_mul_B!(out::AbstractVector, Q::ReinschQ, g::AbstractVecto
end


immutable ReinschR{T} <: AbstractMatrix{T}
struct ReinschR{T} <: AbstractMatrix{T}
h::Vector{T}
end

Expand All @@ -60,7 +60,7 @@ function Base.size(R::ReinschR)
(n-1, n-1)
end

function Base.getindex{T}(R::ReinschR{T}, i::Int, j::Int)
function Base.getindex(R::ReinschR{T}, i::Int, j::Int) where T
# TODO: add check bounds
h = R.h
if (i==j)
Expand All @@ -72,7 +72,7 @@ function Base.getindex{T}(R::ReinschR{T}, i::Int, j::Int)
end
end

function QtQpR{T<:Real}(h::AbstractVector{T}, α::T, w::AbstractVector{T}=ones(T, length(h)+1))
function QtQpR(h::AbstractVector{T}, α::T, w::AbstractVector{T}=ones(T, length(h)+1)) where T<:Real
# maybe has some redundant calculations but should be good enough for now
n = length(h)-1
Q = ReinschQ(h)
Expand All @@ -97,23 +97,12 @@ function QtQpR{T<:Real}(h::AbstractVector{T}, α::T, w::AbstractVector{T}=ones(T
out
end

# temporary hack to use LAPACK wrapper,
# as in
# https://github.com/ApproxFun/BandedMatrices.jl/blob/master/src/blas.jl
if VERSION < v"0.5.0-dev"
macro blasfunc(x)
return :( $(BLAS.blasfunc(x) ))
end
else
import Base.BLAS.@blasfunc
end


for (pbtrf, pbtrs, elty) in
((:dpbtrf_,:dpbtrs_,:Float64),
(:spbtrf_,:spbtrs_,:Float32),
(:zpbtrf_,:zpbtrs_,:Complex128),
(:cpbtrf_,:cpbtrs_,:Complex64))
(:zpbtrf_,:zpbtrs_,:ComplexF64),
(:cpbtrf_,:cpbtrs_,:ComplexF32))
@eval begin
# SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
# * .. Scalar Arguments ..
Expand All @@ -127,10 +116,10 @@ for (pbtrf, pbtrs, elty) in
chkstride1(AB)
n = size(AB, 2)
info = Ref{BlasInt}()
ccall((@blasfunc($pbtrf), liblapack), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
&uplo, &n, &kd, AB, &max(1,stride(AB,2)), info)
ccall((@blasfunc($pbtrf), liblapack), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}),
uplo, n, kd, AB, max(1,stride(AB,2)), info)
AB
end

Expand All @@ -150,12 +139,12 @@ for (pbtrf, pbtrs, elty) in
if n != size(B,1)
throw(DimensionMismatch("Matrix AB has dimensions $(size(AB)), but right hand side matrix B has dimensions $(size(B))"))
end
ccall((@blasfunc($pbtrs), liblapack), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}),
&uplo, &n, &kd, &size(B,2), AB, &max(1,stride(AB,2)),
B, &max(1,stride(B,2)), info)
ccall((@blasfunc($pbtrs), liblapack), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
Ref{BlasInt}),
uplo, n, kd, size(B,2), AB, max(1,stride(AB,2)),
B, max(1,stride(B,2)), info)
B
end
end
Expand Down
21 changes: 12 additions & 9 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
using SmoothingSplines
using Base.Test
using Test

using LinearAlgebra
using Random

#n=50
#X = rand(n) .* 3
Expand All @@ -17,7 +20,7 @@ X = [0.0; cumsum(h)]
λ = 1.0

ws1 = ones(Float64, length(Y))
ws2 = 1.0 + rand(n)
ws2 = 1.0 .+ rand(n)
ws_dict = Dict("uniform weights"=> ws1, "random weights"=>ws2)
for (k,ws) in ws_dict
println(k)
Expand All @@ -26,12 +29,12 @@ for (k,ws) in ws_dict
Qfull = reshape(Float64[x for x in Q], size(Q))
Rfull = reshape(Float64[x for x in R], size(R))
QtYfull = Qfull'*Y
QtY = At_mul_B!(zeros(n-2), Q, Y)
QtY = At_mul_B!(fill(0.0, n-2), Q, Y)

@test QtY QtYfull


tmpQtQpRfull = Rfull + λ*Qfull'*diagm(1./ws)*Qfull
tmpQtQpRfull = Rfull + λ*Qfull'*Matrix(Diagonal(1.0 ./ ws))*Qfull
tmpQtQpR =SmoothingSplines. QtQpR(h, λ, ws)
@test vec(tmpQtQpR[3,:]) diag(tmpQtQpRfull)
@test vec(tmpQtQpR[2,2:end]) diag(tmpQtQpRfull,1)
Expand All @@ -43,8 +46,8 @@ for (k,ws) in ws_dict
SmoothingSplines.pbtrs!('U', 2, tmpQtQpR, γ)
@test γ γfull

gfull = Y - λ*diagm(1./ws)*Qfull*γfull
g = Y - λ*A_mul_B!(zeros(Y), Q, γ)./ws
gfull = Y - λ*Matrix(Diagonal(1.0./ws))*Qfull*γfull
g = Y - λ*A_mul_B!(zero(Y), Q, γ)./ws
@test g gfull

fullfit = fit(SmoothingSpline, X, Y, λ, ws)
Expand All @@ -54,7 +57,7 @@ println("testing predict function")
srand(1)
n=50
X = rand(n) .* 3
Y = 2 .* X.^2 - X .- randn(n)
Y = 2 .* X.^2 .- X .- randn(n)

spl = fit(SmoothingSpline, X, Y, 1.0)

Expand All @@ -76,7 +79,7 @@ Y = [2.0,10.0,4.0,22.0,16.0,10.0,18.0,26.0,34.0,17.0,28.0,14.0,20.0,24.0,
32.0,40.0,50.0,42.0,56.0,76.0,84.0,36.0,46.0,68.0,32.0,48.0,52.0,56.0,
64.0,66.0,54.0,70.0,92.0,93.0,120.0,85.0]
# so that λ for this package is comparable to λ in smooth.spline
X = (X-minimum(X))/(maximum(X)-minimum(X))
X = (X .- minimum(X))/(maximum(X)-minimum(X))

# compare to R.stats:
# attach(cars)
Expand All @@ -87,4 +90,4 @@ Rpred = [ 1.657809, 11.682913, 15.064409, 18.482339, 21.947078, 25.465623, 29.0
32.707552, 36.423468, 40.195168, 44.054188, 48.025219, 52.115627, 56.323874,
60.673887, 69.808295, 74.533703, 79.313597, 84.103688]
cars_fit = fit(SmoothingSpline, X, Y, λ)
@test_approx_eq_eps cars_fit.g Rpred 1e-4
@test cars_fit.g Rpred rtol=1e-4

0 comments on commit 20ee02a

Please sign in to comment.