Update to Julia v1.0 (#68)

Updates for 0.7 & 1.0

.travis.yml

@@ -3,12 +3,11 @@ os:
    - osx
    - linux
julia:
    - 0.6
    - 0.7
    - 1.0
    - nightly
notifications:
    email: false
# script:
#     - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
#     - julia -e 'Pkg.clone(pwd()); Pkg.build("MultivariateStats"); Pkg.test("MultivariateStats"; coverage=true)';
sudo: false
after_success:
    - julia -e 'cd(Pkg.dir("MultivariateStats")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';
  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

Project.toml

@@ -0,0 +1,22 @@
name = "MultivariateStats"
uuid = "6f286f6a-111f-5878-ab1e-185364afe411"
license = "MIT"
repository = "https://github.com/JuliaStats/MultivariateStats.jl.git"
desc = "A Julia package for multivariate statistics and data analysis"
keywords = ["multivariate statistics", "dimensionality reduction"]

version = "0.6.0"

[deps]
Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]

REQUIRE

@@ -1,3 +1,3 @@
julia 0.6
Compat 0.17.0
StatsBase 0.19.0
julia 0.7
StatsBase
Arpack

src/MultivariateStats.jl

@@ -1,13 +1,9 @@
__precompile__()

module MultivariateStats
    using Compat
    using StatsBase

    using LinearAlgebra
    import Statistics: mean, var, cov, covm
    import Base: length, size, show, dump
    import Base.LinAlg: Cholesky
    import StatsBase: fit, predict
    using  Compat: view
    import StatsBase: fit, predict, ConvergenceException
    import SparseArrays

    export

src/cca.jl

@@ -2,17 +2,17 @@

#### CCA Type

type CCA
struct CCA
    xmean::Vector{Float64}  # sample mean of X: of length dx (can be empty)
    ymean::Vector{Float64}  # sample mean of Y: of length dy (can be empty)  
    ymean::Vector{Float64}  # sample mean of Y: of length dy (can be empty)
    xproj::Matrix{Float64}  # projection matrix for X, of size (dx, p)
    yproj::Matrix{Float64}  # projection matrix for Y, of size (dy, p)
    corrs::Vector{Float64}  # correlations, of length p

    function CCA(xm::Vector{Float64}, 
                 ym::Vector{Float64}, 
    function CCA(xm::Vector{Float64},
                 ym::Vector{Float64},
                 xp::Matrix{Float64},
                 yp::Matrix{Float64}, 
                 yp::Matrix{Float64},
                 crs::Vector{Float64})

        dx, px = size(xp)
@@ -24,7 +24,7 @@ type CCA
        isempty(ym) || length(ym) == dy ||
            throw(DimensionMismatch("Incorrect length of ymean."))

        px == py || 
        px == py ||
            throw(DimensionMismatch("xproj and yproj should have the same number of columns."))

        length(crs) == px ||
@@ -50,8 +50,8 @@ correlations(M::CCA) = M.corrs

## use

xtransform{T<:Real}(M::CCA, X::AbstractVecOrMat{T}) = At_mul_B(M.xproj, centralize(X, M.xmean))
ytransform{T<:Real}(M::CCA, Y::AbstractVecOrMat{T}) = At_mul_B(M.yproj, centralize(Y, M.ymean))
xtransform(M::CCA, X::AbstractVecOrMat{T}) where T<:Real = transpose(M.xproj) * centralize(X, M.xmean)
ytransform(M::CCA, Y::AbstractVecOrMat{T}) where T<:Real = transpose(M.yproj) * centralize(Y, M.ymean)

## show & dump

@@ -79,9 +79,9 @@ end

## ccacov

function ccacov(Cxx::DenseMatrix{Float64}, 
function ccacov(Cxx::DenseMatrix{Float64},
                Cyy::DenseMatrix{Float64},
                Cxy::DenseMatrix{Float64}, 
                Cxy::DenseMatrix{Float64},
                xmean::Vector{Float64},
                ymean::Vector{Float64},
                p::Int)
@@ -91,10 +91,10 @@ function ccacov(Cxx::DenseMatrix{Float64},
    dy, dy2 = size(Cyy)
    dx == dx2 || error("Cxx must be a square matrix.")
    dy == dy2 || error("Cyy must be a square matrix.")
    size(Cxy) == (dx, dy) || 
    size(Cxy) == (dx, dy) ||
        throw(DimensionMismatch("size(Cxy) should be equal to (dx, dy)"))

    isempty(xmean) || length(xmean) == dx || 
    isempty(xmean) || length(xmean) == dx ||
        throw(DimensionMismatch("Incorrect length of xmean."))

    isempty(ymean) || length(ymean) == dy ||
@@ -116,17 +116,17 @@ function _ccacov(Cxx, Cyy, Cxy, xmean, ymean, p::Int)
        # solve Px: (Cxy * inv(Cyy) * Cyx) Px = λ Cxx * Px
        # compute Py: inv(Cyy) * Cyx * Px

        G = cholfact(Cyy) \ Cxy'
        Ex = eigfact(Symmetric(Cxy * G), Symmetric(Cxx))
        G = cholesky(Cyy) \ Cxy'
        Ex = eigen(Symmetric(Cxy * G), Symmetric(Cxx))
        ord = sortperm(Ex.values; rev=true)
        vx, Px = extract_kv(Ex, ord, p)
        Py = qnormalize!(G * Px, Cyy)
    else
        # solve Py: (Cyx * inv(Cxx) * Cxy) Py = λ Cyy Py
        # compute Px: inv(Cx) * Cxy * Py

        H = cholfact(Cxx) \ Cxy
        Ey = eigfact(Symmetric(Cxy'H), Symmetric(Cyy))
        H = cholesky(Cxx) \ Cxy
        Ey = eigen(Symmetric(Cxy'H), Symmetric(Cyy))
        ord = sortperm(Ey.values; rev=true)
        vy, Py = extract_kv(Ey, ord, p)
        Px = qnormalize!(H * Py, Cxx)
@@ -143,18 +143,18 @@ end

## ccasvd

function ccasvd(Zx::DenseMatrix{Float64}, 
                Zy::DenseMatrix{Float64}, 
                xmean::Vector{Float64}, 
                ymean::Vector{Float64}, 
function ccasvd(Zx::DenseMatrix{Float64},
                Zy::DenseMatrix{Float64},
                xmean::Vector{Float64},
                ymean::Vector{Float64},
                p::Int)

    dx, n = size(Zx)
    dy, n2 = size(Zy)
    n == n2 || 
    n == n2 ||
        throw(DimensionMismatch("Zx and Zy must have the same number of columns."))

    isempty(xmean) || length(xmean) == dx || 
    isempty(xmean) || length(xmean) == dx ||
        throw(DimensionMismatch("Incorrect length of xmean."))

    isempty(ymean) || length(ymean) == dy ||
@@ -170,7 +170,7 @@ end
#
#   David Weenink.
#   Canonical Correlation Analysis.
#   Institute of Phonetic Sciences, Univ. of Amsterdam, 
#   Institute of Phonetic Sciences, Univ. of Amsterdam,
#   Proceedings 25 (2003), 81-99.
#
#   Note: in this paper, each row is considered as an observation.
@@ -182,27 +182,27 @@ function _ccasvd(Zx, Zy, xmean, ymean, p::Int)
    n = size(Zx, 2)

    # svd decomposition
    Sx = svdfact(Zx)
    Sy = svdfact(Zy)
    S = svdfact!(A_mul_Bt(Sx.Vt, Sy.Vt)) # svd of Vx * Vy'
    Sx = svd(Zx)
    Sy = svd(Zy)
    S = svd!(Sx.Vt * transpose(Sy.Vt)) # svd of Vx * Vy'

    # compute Px and Py
    ord = sortperm(S.S; rev=true)
    si = ord[1:p]
    Px = scale!(Sx.U, 1.0 ./ Sx.S) * S.U[:, si]
    Py = A_mul_Bt(scale!(Sy.U, 1.0 ./ Sy.S), S.Vt[si, :])
    Px = rmul!(Sx.U, Diagonal(1.0 ./ Sx.S)) * S.U[:, si]
    Py = rmul!(Sy.U, Diagonal(1.0 ./ Sy.S)) * S.V[:, si]

    # scale so that Px' * Cxx * Py == I 
    #           and Py' * Cyy * Py == I, 
    # scale so that Px' * Cxx * Py == I
    #           and Py' * Cyy * Py == I,
    #
    # with Cxx = Zx * Zx' / (n - 1)
    #      Cyy = Zy * Zy' / (n - 1)
    #
    scale!(Px, sqrt(n-1))
    scale!(Py, sqrt(n-1))
    rmul!(Px, sqrt(n-1))
    rmul!(Py, sqrt(n-1))

    # compute correlations
    crs = scale!(coldot(Zx'Px, Zy'Py), inv(n-1))
    crs = rmul!(coldot(Zx'Px, Zy'Py), inv(n-1))

    # construct CCA model
    CCA(xmean, ymean, Px, Py, crs)
@@ -212,17 +212,17 @@ end

function fit(::Type{CCA}, X::DenseMatrix{Float64}, Y::DenseMatrix{Float64};
             outdim::Int=min(min(size(X)...), min(size(Y)...)),
             method::Symbol=:svd, 
             xmean=nothing, 
             method::Symbol=:svd,
             xmean=nothing,
             ymean=nothing)

    dx, n = size(X)
    dy, n2 = size(Y)

    n2 == n || 
    n2 == n ||
        throw(DimensionMismatch("X and Y should have the same number of columns."))

    (n >= dx && n >= dy) || 
    (n >= dx && n >= dy) ||
        warn("CCA would be numerically instable when n < dx or n < dy.")

    xmv = preprocess_mean(X, xmean)
@@ -232,9 +232,9 @@ function fit(::Type{CCA}, X::DenseMatrix{Float64}, Y::DenseMatrix{Float64};
    Zy = centralize(Y, ymv)

    if method == :cov
        Cxx = scale!(A_mul_Bt(Zx, Zx), inv(n - 1))
        Cyy = scale!(A_mul_Bt(Zy, Zy), inv(n - 1))
        Cxy = scale!(A_mul_Bt(Zx, Zy), inv(n - 1))
        Cxx = rmul!(Zx*transpose(Zx), inv(n - 1))
        Cyy = rmul!(Zy*transpose(Zy), inv(n - 1))
        Cxy = rmul!(Zx*transpose(Zy), inv(n - 1))
        M = ccacov(Cxx, Cyy, Cxy, xmv, ymv, outdim)
    elseif method == :svd
        M = ccasvd(Zx, Zy, xmv, ymv, outdim)

src/cmds.jl

@@ -2,7 +2,7 @@

## convert Gram matrix to Distance matrix

function gram2dmat!{DT}(D::AbstractMatrix{DT}, G::AbstractMatrix)
function gram2dmat!(D::AbstractMatrix{DT}, G::AbstractMatrix) where DT
    # argument checking
    m = size(G, 1)
    n = size(G, 2)
@@ -23,13 +23,13 @@ function gram2dmat!{DT}(D::AbstractMatrix{DT}, G::AbstractMatrix)
    return D
end

gram2dmat{T<:Real}(G::AbstractMatrix{T}) = gram2dmat!(similar(G, Base.momenttype(T)), G)
gram2dmat(G::AbstractMatrix{T}) where T<:Real = gram2dmat!(similar(G, T), G)

## convert Distance matrix to Gram matrix

function dmat2gram!{GT}(G::AbstractMatrix{GT}, D::AbstractMatrix)
function dmat2gram!(G::AbstractMatrix{GT}, D::AbstractMatrix) where GT
    # argument checking
    n = Compat.LinAlg.checksquare(D)
    n = LinearAlgebra.checksquare(D)
    size(G) == (n, n) ||
        throw(DimensionMismatch("Sizes of G and D do not match."))

@@ -52,34 +52,35 @@ function dmat2gram!{GT}(G::AbstractMatrix{GT}, D::AbstractMatrix)
    return G
end

dmat2gram{T<:Real}(D::AbstractMatrix{T}) = dmat2gram!(similar(D, Base.momenttype(T)), D)
momenttype(T) = typeof((zero(T) * zero(T) + zero(T) * zero(T))/ 2)
dmat2gram(D::AbstractMatrix{T}) where T<:Real = dmat2gram!(similar(D, momenttype(T)), D)

## classical MDS

function classical_mds{T<:Real}(D::AbstractMatrix{T}, p::Int;
        dowarn::Bool=true)
function classical_mds(D::AbstractMatrix{T}, p::Int;
                       dowarn::Bool=true) where T<:Real

    n = size(D, 1)
    m = min(p, n) #Actual number of eigenpairs wanted

    G = dmat2gram(D)

    #Get m largest eigenpairs
    E = eigfact!(Symmetric(G))
    E = eigen!(Symmetric(G))

    #Sometimes dmat2gram produces a negative definite matrix, and the sign just
    #needs to be flipped. The heuristic to check for this robustly is to check
    #if there is a negative eigenvalue of magnitude larger than the largest
    #positive eigenvalue, and flip the sign of eigenvalues if necessary.

    mineig, maxeig = extrema(E[:values])
    mineig, maxeig = extrema(E.values)
    if mineig < 0 && abs(mineig) > abs(maxeig)
        #do flip
        ord = sortperm(E.values)
        v = -E[:values][ord[1:m]]
        v = -E.values[ord[1:m]]
    else
        ord = sortperm(E.values; rev=true)
        v = E[:values][ord[1:m]]
        v = E.values[ord[1:m]]
    end

    for i = 1:m
@@ -98,14 +99,14 @@ function classical_mds{T<:Real}(D::AbstractMatrix{T}, p::Int;

    #Check if the last considered eigenvalue is degenerate
    if m>0
        nevalsmore = sum(abs.(E[:values][ord[m+1:end]] .- v[m]^2) .< n*eps())
        nevals = sum(abs.(E[:values] .- v[m]^2) .< n*eps())
        nevalsmore = sum(abs.(E.values[ord[m+1:end]] .- v[m]^2) .< n*eps())
        nevals = sum(abs.(E.values .- v[m]^2) .< n*eps())
        if nevalsmore > 1
            dowarn && warn("The last eigenpair is degenerate with $(nevals-1) others; $nevalsmore were ignored. Answer is not unique")
        end
    end
    U = E[:vectors][:, ord[1:m]]
    scale!(U, v)
    U = E.vectors[:, ord[1:m]]
    rmul!(U, Diagonal(v))

    #Add trailing zero coordinates if dimension of embedding space (p) exceeds
    #number of eigenpairs used (m)
@@ -115,4 +116,3 @@ function classical_mds{T<:Real}(D::AbstractMatrix{T}, p::Int;

    U' #Return each coordinate in a column
end