Merge branch 'master' of https://github.com/alan-turing-institute/MLJLinearModels.jl

Author: tlienart

.travis.yml

@@ -4,10 +4,11 @@ os:
   - linux
 julia:
   - 1.0
-  - 1.1
-  - 1.2
-  - 1.3
+  - 1.4
   - nightly
+matrix:
+  allow_failures:
+    - julia: nightly
 notifications:
   email: false
 after_success:

Project.toml

@@ -1,7 +1,7 @@
 name = "MLJLinearModels"
 uuid = "6ee0df7b-362f-4a72-a706-9e79364fb692"
 authors = ["Thibaut Lienart <[email protected]>"]
-version = "0.3.1"
+version = "0.3.3"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -16,8 +16,8 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 DocStringExtensions = "^0.8"
 IterativeSolvers = "^0.8"
 LinearMaps = "^2.6"
-MLJModelInterface = "^0.1"
-Optim = "^0.20"
+MLJModelInterface = "^0.1,^0.2"
+Optim = "^0.20,^0.21"
 Parameters = "^0.12"
 julia = "^1"
 
@@ -28,7 +28,8 @@ PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
 RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DelimitedFiles", "PyCall", "Test", "Random", "RDatasets", "RCall", "MLJBase"]
+test = ["DelimitedFiles", "PyCall", "Test", "Random", "RDatasets", "RCall", "MLJBase", "StableRNGs"]

src/fit/analytical.jl

@@ -11,35 +11,35 @@ Fit a least square regression either with no penalty (OLS) or with a L2 penalty
 Assuming `n` dominates `p`,
 
 * non-iterative (full solve):     O(np²) - dominated by the construction of the
-								  Hessian X'X.
+                                  Hessian X'X.
 * iterative (conjugate gradient): O(κnp) - with κ the number of CG steps
-							 	  (κ ≤ p).
+                                  (κ ≤ p).
 """
 function _fit(glr::GLR{L2Loss,<:L2R}, solver::Analytical, X, y)
-	# full solve
-	if !solver.iterative
-		λ  = getscale(glr.penalty)
-		if iszero(λ)
-			# standard LS solution
-			return augment_X(X, glr.fit_intercept) \ y
-		else
-			# Ridge case -- form the Hat Matrix then solve
-			H = form_XtX(X, glr.fit_intercept, λ)
-			b = X'y
-			glr.fit_intercept && (b = vcat(b, sum(y)))
-			return cholesky!(H) \ b
-		end
-	end
-	# Iterative case, note that there is no augmentation here
-	# it is done implicitly in the application of the Hessian to
-	# avoid copying X
-	# The number of CG steps to convergence is at most `p`
-	p = size(X, 2) + Int(glr.fit_intercept)
-	max_cg_steps = min(solver.max_inner, p)
-	# Form the Hessian map, cost of application H*v is O(np)
-	Hm = LinearMap(Hv!(glr, X, y), p;
-				   ismutating=true, isposdef=true, issymmetric=true)
-	b  = X'y
-	glr.fit_intercept && (b = vcat(b, sum(y)))
-	return cg(Hm, b; maxiter=max_cg_steps)
+    # full solve
+    if !solver.iterative
+        λ  = getscale(glr.penalty)
+        if iszero(λ)
+            # standard LS solution
+            return augment_X(X, glr.fit_intercept) \ y
+        else
+            # Ridge case -- form the Hat Matrix then solve
+            H = form_XtX(X, glr.fit_intercept, λ)
+            b = X'y
+            glr.fit_intercept && (b = vcat(b, sum(y)))
+            return cholesky!(H) \ b
+        end
+    end
+    # Iterative case, note that there is no augmentation here
+    # it is done implicitly in the application of the Hessian to
+    # avoid copying X
+    # The number of CG steps to convergence is at most `p`
+    p = size(X, 2) + Int(glr.fit_intercept)
+    max_cg_steps = min(solver.max_inner, p)
+    # Form the Hessian map, cost of application H*v is O(np)
+    Hm = LinearMap(Hv!(glr, X, y), p;
+                   ismutating=true, isposdef=true, issymmetric=true)
+    b  = X'y
+    glr.fit_intercept && (b = vcat(b, sum(y)))
+    return cg(Hm, b; maxiter=max_cg_steps)
 end

src/fit/default.jl

@@ -10,13 +10,13 @@ export fit
 _solver(::GLR{L2Loss,<:L2R}, np::NTuple{2,Int}) = Analytical()
 
 # Logistic, Multinomial
-_solver(::GLR{LogisticLoss,<:L2R}, 	  np::NTuple{2,Int}) = LBFGS()
+_solver(::GLR{LogisticLoss,<:L2R},    np::NTuple{2,Int}) = LBFGS()
 _solver(::GLR{MultinomialLoss,<:L2R}, np::NTuple{2,Int}) = LBFGS()
 
 # Lasso, ElasticNet, Logistic, Multinomial
 function _solver(glr::GLR{<:SmoothLoss,<:ENR}, np::NTuple{2,Int})
-	(is_l1(glr.penalty) || is_elnet(glr.penalty)) && return FISTA()
-	@error "Not yet implemented."
+    (is_l1(glr.penalty) || is_elnet(glr.penalty)) && return FISTA()
+    @error "Not yet implemented."
 end
 
 # Robust, Quantile
@@ -34,19 +34,19 @@ Fit a generalised linear regression model using an appropriate solver based on
 the loss and penalty of the model. A method can, in some cases, be specified.
 """
 function fit(glr::GLR, X::AbstractMatrix{<:Real}, y::AVR;
-			 solver::Solver=_solver(glr, size(X)))
+             solver::Solver=_solver(glr, size(X)))
     check_nrows(X, y)
-	n, p = size(X)
-	p += Int(glr.fit_intercept)
-	# allocate cache for temporary computations of size n/p
-	# which are frequent but otherwise un-important so that
-	# we can reduce the overall number of allocations
-	# these are const Refs defined when the module is loaded
-	c = glr.loss isa MultinomialLoss ? maximum(y) : 0
-	allocate(n, p, c)
-	# effective call to fit routine
+    n, p = size(X)
+    p += Int(glr.fit_intercept)
+    # allocate cache for temporary computations of size n/p
+    # which are frequent but otherwise un-important so that
+    # we can reduce the overall number of allocations
+    # these are const Refs defined when the module is loaded
+    c = glr.loss isa MultinomialLoss ? maximum(y) : 0
+    allocate(n, p, c)
+    # effective call to fit routine
     θ = _fit(glr, solver, X, y)
-	# de-allocate cache
-	deallocate()
-	return θ
+    # de-allocate cache
+    deallocate()
+    return θ
 end

src/glr/constructors.jl

@@ -1,9 +1,9 @@
 export GeneralizedLinearRegression, GLR,
         LinearRegression, RidgeRegression,
         LassoRegression, ElasticNetRegression,
-        LADRegression, MADRegression,
-        LogisticRegression, MultinomialRegression,
-        RobustRegression, HuberRegression, QuantileRegression
+        LADRegression, LogisticRegression,
+        MultinomialRegression, RobustRegression,
+        HuberRegression, QuantileRegression
 
 """
     GeneralizedLinearRegression{L<:Loss, P<:Penalty}

src/mlj/interface.jl

@@ -19,19 +19,20 @@ const ALL_MODELS = (REG_MODELS..., CLF_MODELS...)
 
 function MMI.fit(m::Union{REG_MODELS...}, verb::Int, X, y)
     Xmatrix = MMI.matrix(X)
+    features = (sch = MMI.schema(X)) === nothing ? nothing : sch.names
     reg     = glr(m)
     solver  = m.solver === nothing ? _solver(reg, size(Xmatrix)) : m.solver
     # get the parameters
     θ = fit(reg, Xmatrix, y; solver=solver)
     # return
-    return θ, nothing, NamedTuple{}()
+    return (θ, features), nothing, NamedTuple{}()
 end
 
-MMI.predict(m::Union{REG_MODELS...}, θ, Xnew) = apply_X(MMI.matrix(Xnew), θ)
+MMI.predict(m::Union{REG_MODELS...}, (θ, features), Xnew) = apply_X(MMI.matrix(Xnew), θ)
 
-function MMI.fitted_params(m::Union{REG_MODELS...}, θ)
-    m.fit_intercept && return (coefs = θ[1:end-1], intercept = θ[end])
-    return (coefs = θ, intercept = nothing)
+function MMI.fitted_params(m::Union{REG_MODELS...}, (θ, features))
+    m.fit_intercept && return (coefs = coef_vec(θ[1:end-1], features), intercept = θ[end])
+    return (coefs = coef_vec(θ, features), intercept = nothing)
 end
 
 #= ===========
@@ -40,6 +41,7 @@ end
 
 function MMI.fit(m::Union{CLF_MODELS...}, verb::Int, X, y)
     Xmatrix  = MMI.matrix(X)
+    features = (sch = MMI.schema(X)) === nothing ? nothing : sch.names
     yplain   = convert.(Int, MMI.int(y))
     classes  = MMI.classes(y[1])
     nclasses = length(classes)
@@ -56,10 +58,10 @@ function MMI.fit(m::Union{CLF_MODELS...}, verb::Int, X, y)
     # get the parameters
     θ = fit(clf, Xmatrix, yplain, solver=solver)
     # return
-    return (θ, c, classes), nothing, NamedTuple{}()
+    return (θ, features, c, classes), nothing, NamedTuple{}()
 end
 
-function MMI.predict(m::Union{CLF_MODELS...}, (θ, c, classes), Xnew)
+function MMI.predict(m::Union{CLF_MODELS...}, (θ, features, c, classes), Xnew)
     Xmatrix = MMI.matrix(Xnew)
     preds   = apply_X(Xmatrix, θ, c)
     # binary classification
@@ -73,20 +75,29 @@ function MMI.predict(m::Union{CLF_MODELS...}, (θ, c, classes), Xnew)
     return [MMI.UnivariateFinite(classes, preds[i, :]) for i in 1:size(Xmatrix,1)]
 end
 
-function MMI.fitted_params(m::Union{CLF_MODELS...}, (θ, c, classes))
+function MMI.fitted_params(m::Union{CLF_MODELS...}, (θ, features, c, classes))
+    function _fitted_params(coefs, features, intercept)
+        return (classes = classes, coefs = coef_vec(coefs, features), intercept = intercept)
+    end
     if c > 1
+        W = reshape(θ, :, c)
         if m.fit_intercept
-            W = reshape(θ, div(length(θ), c), c)
-            return (coefs = W, intercept = nothing)
+            return _fitted_params(W, features, W[end, :])
         end
-        W = reshape(θ, p+1, c)
-        return (coefs = W[1:p, :], intercept = W[end, :])
+        return _fitted_params(W[1:end-1, :], features, nothing)
     end
     # single class
-    m.fit_intercept && return (coefs = θ[1:end-1], intercept = θ[end])
-    return (coefs = θ, intercept = nothing)
+    m.fit_intercept && return _fitted_params(θ[1:end-1], features, θ[end])
+    return _fitted_params(θ, features, nothing)
 end
 
+@static VERSION < v"1.1" && (eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1)))
+
+coef_vec(W::AbstractMatrix, features) = [feature => coef for (feature, coef) in zip(features, eachrow(W))]
+coef_vec(θ::AbstractVector, features) = [feature => coef for (feature, coef) in zip(features, θ)]
+coef_vec(W::AbstractMatrix, ::Nothing) = W
+coef_vec(θ::AbstractVector, ::Nothing) = θ
+
 #= =======================
    METADATA FOR ALL MODELS
    ======================= =#