Merge pull request #34 from fundamental/compatiability

WIP Fix 0.7/1.0 Incompatiabilities

.travis.yml

@@ -1,12 +1,14 @@
 sudo: false
 language: julia
 julia:
- - 0.6
+ - 0.7
+ - 1.0
+ - nightly
 before_install:
  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
 addons:
   apt:
     packages:
       - hdf5-tools
 script:
- - julia --check-bounds=yes -e 'versioninfo(); Pkg.init(); Pkg.clone(pwd()); Pkg.build("GaussianMixtures"); Pkg.test("GaussianMixtures", coverage=VERSION < v"0.4.0-dev")'
+ - julia --check-bounds=yes -e 'using Pkg; Pkg.clone(pwd()); Pkg.build("GaussianMixtures"); Pkg.test("GaussianMixtures", coverage=true)'

REQUIRE

@@ -1,6 +1,7 @@
-julia 0.6
+julia 0.7
 Clustering 0.6.0
 Distributions
+StatsBase
 PDMats
 Compat
 JLD

src/bayes.jl

@@ -13,9 +13,10 @@
 ## e.g., W₀⁻¹ where others might write W0inv.
 
 using SpecialFunctions: digamma
+using SpecialFunctions: lgamma
 
 ## initialize a prior with minimal knowledge
-function GMMprior{T<:AbstractFloat}(d::Int, alpha::T, beta::T)
+function GMMprior(d::Int, alpha::T, beta::T) where {T<:AbstractFloat}
     m₀ = zeros(T, d)
     W₀ = eye(T, d)
     ν₀ = convert(T,d)
@@ -24,7 +25,7 @@ end
 Base.copy(p::GMMprior) = GMMprior(p.α₀, p.β₀, copy(p.m₀), p.ν₀, copy(p.W₀))
 
 ## initialize from a GMM and nₓ, the number of points used to train the GMM.
-function VGMM{T}(g::GMM{T}, π::GMMprior{T})
+function VGMM(g::GMM{T}, π::GMMprior{T}) where {T}
     nₓ = g.nx
     N = g.w * nₓ
     mx = g.μ
@@ -61,13 +62,13 @@ function GMM(vg::VGMM)
 end
 
 ## m-step given prior and stats
-function mstep{T}(π::GMMprior, N::Vector{T}, mx::Matrix{T}, S::Vector)
+function mstep(π::GMMprior, N::Vector{T}, mx::Matrix{T}, S::Vector) where {T}
     ng = length(N)
-    α = π.α₀ + N                # ng, 10.58
-    ν = π.ν₀ + N + 1            # ng, 10.63
-    β = π.β₀ + N                # ng, 10.60
+    α = π.α₀ .+ N               # ng, 10.58
+    ν = π.ν₀ .+ N .+ 1          # ng, 10.63
+    β = π.β₀ .+ N               # ng, 10.60
     m = similar(mx)             # ng × d
-    W = Array{eltype(FullCov{T})}(ng) # ng × (d*d)
+    W = Array{eltype(FullCov{T})}(undef, ng) # ng × (d*d)
     d = size(mx,2)
     limit = √ eps(eltype(N))
     W₀⁻¹ = inv(π.W₀)
@@ -77,10 +78,11 @@ function mstep{T}(π::GMMprior, N::Vector{T}, mx::Matrix{T}, S::Vector)
             Δ = vec(mx[k,:]) - π.m₀ ## v0.5 arraymageddon
             ## do some effort to keep the matrix positive definite
             third = π.β₀ * N[k] / (π.β₀ + N[k]) * (Δ * Δ') # guarantee symmety in Δ Δ'
-            W[k] = chol(inv(cholfact(Symmetric(W₀⁻¹ + N[k]*S[k] + third)))) # 10.62
+            W[k] = cholesky(inv(cholesky(Symmetric(W₀⁻¹ + N[k]*S[k] +
+                                                   third)))).U # 10.62
         else
             m[k,:] = zeros(d)
-            W[k] = chol(eye(d))
+            W[k] = cholesky(eye(d))
         end
     end
     return α, β, m, ν, W
@@ -92,7 +94,7 @@ function expectations(vg::VGMM)
     Elogπ = digamma.(vg.α) .- digamma(sum(vg.α)) # 10.66, size ng
     ElogdetΛ = similar(Elogπ) # size ng
     for k in 1:ng
-        ElogdetΛ[k] = sum(digamma.(0.5(vg.ν[k] .+ 1 - collect(1:d)))) + d*log(2) + mylogdet(vg.W[k]) # 10.65
+        ElogdetΛ[k] = sum(digamma.(0.5(vg.ν[k] .+ 1 .- collect(1:d)))) + d*log(2) + mylogdet(vg.W[k]) # 10.65
     end
     return Elogπ, ElogdetΛ
 end
@@ -114,7 +116,7 @@ function logρ(vg::VGMM, x::Matrix, ex::Tuple)
         ### d/vg.β[k] + vg.ν[k] * (x_i - m_k)' W_k (x_i = m_k) forall i
         ## Δ = (x_i - m_k)' W_k (x_i = m_k)
         xμTΛxμ!(Δ, x, vec(vg.m[k,:]), vg.W[k])
-        EμΛ[:,k] = d/vg.β[k] .+ vg.ν[k] * sum(abs2, Δ, 2)
+        EμΛ[:,k] = d/vg.β[k] .+ vg.ν[k] * sum(abs2, Δ, dims=2)
     end
     Elogπ, ElogdetΛ = ex
     (Elogπ + 0.5ElogdetΛ .- 0.5d*log(2π))' .- 0.5EμΛ
@@ -147,15 +149,15 @@ end
 ## Like for the GMM, we return nₓ, (some value), zeroth, first, second order stats
 ## All return values can be accumulated, except r, which we need for
 ## lowerbound ElogpZπ and ElogqZ
-function stats{T}(vg::VGMM, x::Matrix{T}, ex::Tuple)
+function stats(vg::VGMM, x::Matrix{T}, ex::Tuple) where {T}
     ng = vg.n
     (nₓ, d) = size(x)
     if nₓ == 0
         return 0, zero(RT), zeros(RT, ng), zeros(RT, ng, d), [zeros(RT, d,d) for k=1:ng]
     end
     r, ElogpZπqZ = rnk(vg, x, ex) # nₓ × ng
-    N = vec(sum(r, 1))          # ng
-    F = r' * x                  # ng × d
+    N = vec(sum(r, dims=1))       # ng
+    F = r' * x                    # ng × d
     ## S_k = sum_i r_ik x_i x_i'
     ## Sm = x' * hcat([broadcast(*, r[:,k], x) for k=1:ng]...)
     SS = similar(x, nₓ, d*ng)   # nₓ*nd*ng mem, nₓ*nd*ng multiplications
@@ -183,7 +185,7 @@ function stats(vg::VGMM, d::Data, ex::Tuple; parallel=false)
 end
 
 ## trace(A*B) = sum(A' .* B)
-function trAB{T1,T2}(A::Matrix{T1}, B::Matrix{T2})
+function trAB(A::Matrix{T1}, B::Matrix{T2}) where {T1,T2}
     RT = promote_type(T1,T2)
     nr, nc = size(A)
     size(B) == (nc, nr) || error("Inconsistent matrix size")
@@ -205,7 +207,8 @@ function lowerbound(vg::VGMM, N::Vector, mx::Matrix, S::Vector,
     α, β, m, ν, W = vg.α, vg.β, vg.m, vg.ν, vg.W # VGMM vars
     gaussians = 1:ng
     ## B.79
-    logB(W, ν) = -0.5ν * (mylogdet(W) + d*log(2)) - d*(d-1)/4*log(π) - sum(lgamma.(0.5(ν + 1 - collect(1:d))))
+    logB(W, ν) = -0.5ν * (mylogdet(W) + d*log(2)) .- d*(d-1)/4*log(π) -
+    sum(lgamma.(0.5(ν + 1 .- collect(1:d))))
     ## E[log p(x|Ζ,μ,Λ)] 10.71
     Elogll = 0.
     for k in gaussians
@@ -215,7 +218,7 @@ function lowerbound(vg::VGMM, N::Vector, mx::Matrix, S::Vector,
                              - ν[k] * (trAB(S[k], Wk) + Wk * Δ ⋅ Δ)) # check chol efficiency
     end                        # 10.71
     ## E[log p(Z|π)] from rnk() 10.72
-    Elogpπ = lgamma(ng*α₀) - ng*lgamma(α₀) - (α₀-1)sum(Elogπ) # E[log p(π)] 10.73
+    Elogpπ = lgamma(ng*α₀) .- ng*lgamma(α₀) .- (α₀-1)sum(Elogπ) # E[log p(π)] 10.73
     ElogpμΛ = ng*logB(W₀, ν₀)   # E[log p(μ, Λ)] B.79
     for k in gaussians
         Δ = vec(m[k,:]) - m₀        # d ## v0.5 arraymageddon

src/compat.jl

@@ -1,21 +1,23 @@
+using LinearAlgebra
+
 UTriangular(a::Matrix) = UpperTriangular(a)
-if VERSION < v"0.5.0-dev"
-    Base.cholfact(s::Symmetric) = cholfact(full(s))
-end
 
 ## NumericExtensions is no longer supported, underoptimized implementation:
-function logsumexp{T<:AbstractFloat}(x::AbstractVector{T})
+function logsumexp(x::AbstractVector{T}) where {T<:AbstractFloat}
     m = maximum(x)
     log(sum(exp.(x .- m))) + m
 end
-logsumexp{T<:AbstractFloat}(x::Matrix{T}, dim::Integer) = mapslices(logsumexp, x, dim)
+logsumexp(x::Matrix{T}, dim::Integer) where {T<:AbstractFloat} = mapslices(logsumexp, x, dims=dim)
+
+eye(n::Int) = Matrix(1.0I,n,n)
+eye(::Type{Float64}, n::Int) = Matrix(1.0I,n,n)
 
 ## Also NumericExtensions' semantics of dot() is no longer supported.
-function Base.dot{T<:AbstractFloat}(x::Matrix{T}, y::Matrix{T})
+function LinearAlgebra.dot(x::Matrix{T}, y::Matrix{T}) where {T<:AbstractFloat}
     size(x) == size(y) || error("Matrix sizes must match")
     dot(vec(x), vec(y))
 end
-function Base.dot{T<:AbstractFloat}(x::Matrix{T}, y::Matrix{T}, dim::Integer)
+function LinearAlgebra.dot(x::Matrix{T}, y::Matrix{T}, dim::Integer) where {T<:AbstractFloat}
     size(x) == size(y) || error("Matrix sizes must match")
     if dim==1
         r = zeros(T, 1, size(x,2))
@@ -34,13 +36,3 @@ function Base.dot{T<:AbstractFloat}(x::Matrix{T}, y::Matrix{T}, dim::Integer)
     end
     r
 end
-
-if VERSION < v"0.5.0-dev+2023"
-    displaysize(io::IO) = Base.tty_size()
-end
-
-## this we need for xμTΛxμ!
-#Base.A_mul_Bc!(A::StridedMatrix{Float64}, B::AbstractTriangular{Float32}) = A_mul_Bc!(A, convert(AbstractMatrix{Float64}, B))
-#Base.A_mul_Bc!(A::Matrix{Float32}, B::AbstractTriangular{Float64}) = A_mul_Bc!(A, convert(AbstractMatrix{Float32}, B))
-## this for diagstats
-#Base.BLAS.gemm!(a::Char, b::Char, alpha::Float64, A::Matrix{Float32}, B::Matrix{Float64}, beta::Float64, C::Matrix{Float64}) = Base.BLAS.gemm!(a, b, alpha, float64(A), B, beta, C)

src/data.jl

@@ -1,16 +1,18 @@
 ## data.jl Julia code to handle matrix-type data on disc
+using Statistics
+using Distributed
 
 ## report the kind of Data structure from the type instance
-kind{T,S<:AbstractString}(d::Data{T,S}) = :file
-kind{T}(d::Data{T,Matrix{T}}) = :matrix
-Base.eltype{T}(d::Data{T}) = T
+kind(d::Data{T,S}) where {T,S<:AbstractString} = :file
+kind(d::Data{T,Matrix{T}}) where {T} = :matrix
+Base.eltype(d::Data{T}) where {T} = T
 
 ## constructor for a plain matrix.  rowvectors: data points x represented as rowvectors
-Data{T}(x::Matrix{T}) = Data(Matrix{T}[x])
+Data(x::Matrix{T}) where {T} = Data(Matrix{T}[x])
 
 ## constructor for a vector of files
 ## Data([strings], type, loadfunction)
-function Data{S<:AbstractString}(files::Vector{S}, datatype::DataType, load::Function)
+function Data(files::Vector{S}, datatype::DataType, load::Function) where {S<:AbstractString}
     Data(files, datatype, @compat Dict(:load => load))
 end
 
@@ -32,7 +34,7 @@ function JLD.save(file::AbstractString, x::Matrix)
 end
 
 ## Data([strings], type; load=loadfunction, size=sizefunction)
-function Data{S<:AbstractString}(files::Vector{S}, datatype::DataType; kwargs...)
+function Data(files::Vector{S}, datatype::DataType; kwargs...) where {S<:AbstractString}
     all([isa((k,v), (Symbol,Function)) for (k,v) in kwargs]) || error("Wrong type of argument", args)
     d = Dict{Symbol,Function}([kwargs...])
     if !haskey(d, :load)
@@ -60,15 +62,21 @@ function Base.getindex(x::Data, i::Int)
 end
 
 ## A range as index (including [1:1]) returns a Data object, not the data.
-function Base.getindex{T,VT}(x::Data{T,VT}, r::Range)
+function Base.getindex(x::Data{T,VT}, r::UnitRange{Int}) where {T,VT}
     Data{T, VT}(x.list[r], x.API)
 end
 
 ## define an iterator for Data
 Base.length(x::Data) = length(x.list)
-Base.start(x::Data) = 0
-Base.next(x::Data, state::Int) = x[state+1], state+1
-Base.done(x::Data, state::Int) = state == length(x)
+function Base.iterate(x::Data, state=1)
+    count = state
+
+    if(count > length(x.list))
+        return nothing
+    else
+        return (x[count], count+1)
+    end
+end
 
 ## This function is like pmap(), but executes each element of Data on a predestined
 ## worker, so that file caching at the local machine is beneficial.
@@ -105,7 +113,7 @@ end
 function dmapreduce(f::Function, op::Function, x::Data)
     nₓ = length(x)
     nw = nworkers()
-    results = Array{Any}(nw) ## will contain pointers for parallel return value.
+    results = Array{Any}(undef, nw) ## will contain pointers for parallel return value.
     valid = Any[false for i=1:nw] # can't use bitarray as parallel return value, must be pointers
     id=0
     nextid() = (id += 1)
@@ -132,21 +140,21 @@ function dmapreduce(f::Function, op::Function, x::Data)
             end
         end
     end
-    reduce(op, results[find(valid)])
+    reduce(op, results[findall(valid)])
 end
 
 ## stats: compute nth order stats for array (this belongs in stats.jl)
-function stats{T<:AbstractFloat}(x::Matrix{T}, order::Int=2; kind=:diag, dim=1)
+function stats(x::Matrix{T}, order::Int=2; kind=:diag, dim=1) where {T<:AbstractFloat}
     if dim==1
         n, d = size(x)
     else
         d, n = size(x)
     end
     if kind == :diag
         if order == 2
-            return n, vec(sum(x, dim)), vec(sum(abs2, x, dim))
+            return n, vec(sum(x, dims=dim)), vec(sum(abs2, x, dims=dim))
         elseif order == 1
-            return n, vec(sum(x, dim))
+            return n, vec(sum(x, dims=dim))
         else
             sx = zeros(T, order, d)
             for j=1:d
@@ -168,7 +176,7 @@ function stats{T<:AbstractFloat}(x::Matrix{T}, order::Int=2; kind=:diag, dim=1)
     elseif kind == :full
         order == 2 || error("Can only do covar starts for order=2")
         ## lazy implementation
-        sx = vec(sum(x, dim))
+        sx = vec(sum(x, dims=dim))
         sxx = x' * x
         return n, sx, sxx
     else
@@ -218,17 +226,17 @@ end
 
 Base.sum(d::Data, dim::Int) = retranspose(stats(d,1, dim=dim)[2], dim)
 
-function Base.mean(d::Data)
+function Statistics.mean(d::Data)
     n, sx = stats(d, 1)
     sum(sx) / (n*length(sx))
 end
 
- function Base.mean(d::Data, dim::Int)
+ function Statistics.mean(d::Data, dim::Int)
      n, sx = stats(d, 1, dim=dim)
      return retranspose(sx ./ n, dim)
 end
 
-function Base.var(d::Data)
+function Statistics.var(d::Data)
     n, sx, sxx = stats(d, 2)
     n *= length(sx)
     ssx = sum(sx)                       # keep type stability...
@@ -237,13 +245,13 @@ function Base.var(d::Data)
     return (ssxx - n*μ^2) / (n - 1)
 end
 
-function Base.var(d::Data, dim::Int)
+function Statistics.var(d::Data, dim::Int)
     n, sx, sxx = stats(d, 2, dim=dim)
     μ = sx ./ n
     return retranspose((sxx - n*μ.^2) ./ (n-1), dim)
 end
 
-function Base.cov(d::Data)
+function Statistics.cov(d::Data)
     n, sx, sxx = stats(d, 2, kind=:full)
     μ = sx ./ n
     (sxx - n*μ*μ') ./ (n-1)
@@ -278,7 +286,7 @@ end
 
 Base.collect(d::Data) = vcat([x for x in d]...)
 
-function Base.convert{Td,Ts}(::Type{Data{Td}}, d::Data{Ts})
+function Base.convert(::Type{Data{Td}}, d::Data{Ts}) where {Td,Ts}
     Td == Ts && return d
     if kind(d) == :file
         api = copy(d.API)

src/distributions.jl

@@ -20,7 +20,7 @@ function GMM(m::MixtureModel{Multivariate,Continuous,DiagNormal})
 end
 
 ## conversion to MixtureModel
-function Distributions.MixtureModel{T<:AbstractFloat}(gmm::GMM{T})
+function Distributions.MixtureModel(gmm::GMM{T}) where {T<:AbstractFloat}
     if gmm.d == 1
         mixtures = [Normal(gmm.μ[i,1], gmm.Σ[i,1]) for i=1:gmm.n]
     elseif kind(gmm) == :full