diff --git a/.travis.yml b/.travis.yml index d44792f..49586a2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,8 +4,7 @@ os: - linux - osx julia: - - 0.5 - - 0.6 + - 0.7 - nightly matrix: allow_failures: diff --git a/REQUIRE b/REQUIRE index 8be8c67..4cdb5d8 100644 --- a/REQUIRE +++ b/REQUIRE @@ -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 diff --git a/appveyor.yml b/appveyor.yml index e52d22d..eb96cb4 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -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: diff --git a/src/SmoothingSplines.jl b/src/SmoothingSplines.jl index 8d69ba9..a0c97c0 100644 --- a/src/SmoothingSplines.jl +++ b/src/SmoothingSplines.jl @@ -4,6 +4,7 @@ module SmoothingSplines import StatsBase: fit!, fit, RegressionModel, rle, ordinalrank, mean using Reexport +using LinearAlgebra export SmoothingSpline @@ -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} @@ -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 @@ -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] @@ -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 @@ -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) @@ -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 diff --git a/src/matrices.jl b/src/matrices.jl index 500b463..653d11d 100644 --- a/src/matrices.jl +++ b/src/matrices.jl @@ -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 @@ -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] @@ -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()) @@ -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()) @@ -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 @@ -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) @@ -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) @@ -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 .. @@ -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 @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 1d3312d..5718d63 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,8 @@ using SmoothingSplines -using Base.Test +using Test + +using LinearAlgebra +using Random #n=50 #X = rand(n) .* 3 @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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