Skip to content

Commit b6cc040

Browse files
tlienarttlienart
authored
removing all global pointers (#67)
* removing all global pointers * patch release
1 parent d3e9994 commit b6cc040

File tree

11 files changed

+118
-205
lines changed

11 files changed

+118
-205
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "MLJLinearModels"
22
uuid = "6ee0df7b-362f-4a72-a706-9e79364fb692"
33
authors = ["Thibaut Lienart <[email protected]>"]
4-
version = "0.3.4"
4+
version = "0.3.5"
55

66
[deps]
77
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

‎src/MLJLinearModels.jl‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,6 @@ const MMI = MLJModelInterface
1414
const AVR = AbstractVector{<:Real}
1515
const Option{T} = Union{Nothing,T}
1616

17-
include("scratchspace.jl")
18-
1917
include("utils.jl")
2018

2119
# > Loss / penalty definitions <

‎src/fit/default.jl‎

Lines changed: 1 addition & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -36,17 +36,5 @@ the loss and penalty of the model. A method can, in some cases, be specified.
3636
function fit(glr::GLR, X::AbstractMatrix{<:Real}, y::AVR;
3737
solver::Solver=_solver(glr, size(X)))
3838
check_nrows(X, y)
39-
n, p = size(X)
40-
p += Int(glr.fit_intercept)
41-
# allocate cache for temporary computations of size n/p
42-
# which are frequent but otherwise un-important so that
43-
# we can reduce the overall number of allocations
44-
# these are const Refs defined when the module is loaded
45-
c = glr.loss isa MultinomialLoss ? maximum(y) : 0
46-
allocate(n, p, c)
47-
# effective call to fit routine
48-
θ = _fit(glr, solver, X, y)
49-
# de-allocate cache
50-
deallocate()
51-
return θ
39+
return _fit(glr, solver, X, y)
5240
end

‎src/glr/d_l2loss.jl‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,10 @@
1515
function Hv!(glr::GLR{L2Loss,<:L2R}, X, y)
1616
n, p = size(X)
1717
λ = getscale(glr.penalty)
18+
# scratch allocation
19+
SCRATCH_N = zeros(n)
1820
if glr.fit_intercept
21+
SCRATCH_P = zeros(p)
1922
# H = [X 1]'[X 1] + λ I
2023
# rows a 1:p = [X'X + λI | X'1]
2124
# row e end = [1'X | n+λι] where ι is 1 if glr.penalize_intercept
@@ -25,11 +28,11 @@ function Hv!(glr::GLR{L2Loss,<:L2R}, X, y)
2528
a = 1:p
2629
Hvₐ = view(Hv, a)
2730
vₐ = view(v, a)
28-
Xt1 = view(SCRATCH_P[], a)
31+
Xt1 = SCRATCH_P
2932
copyto!(Xt1, sum(X, dims=1)) # -- X'1 (note: sum will allocate)
3033
vₑ = v[end]
3134
# update for the first p rows -- (X'X + λI)v[1:p] + (X'1)v[end]
32-
Xvₐ = SCRATCH_N[]
35+
Xvₐ = SCRATCH_N
3336
mul!(Xvₐ, X, vₐ)
3437
mul!(Hvₐ, X', Xvₐ)
3538
Hvₐ .+= λ .* vₐ .+ Xt1 .* vₑ
@@ -38,7 +41,7 @@ function Hv!(glr::GLR{L2Loss,<:L2R}, X, y)
3841
end
3942
else
4043
(Hv, v) -> begin
41-
Xv = SCRATCH_N[]
44+
Xv = SCRATCH_N
4245
mul!(Xv, X, v) # -- Xv
4346
mul!(Hv, X', Xv) # -- X'Xv
4447
Hv .+= λ .* v # -- X'Xv + λv
@@ -59,9 +62,11 @@ end
5962

6063
function smooth_fg!(glr::GLR{L2Loss,<:ENR}, X, y)
6164
λ = getscale_l2(glr.penalty)
65+
# scratch allocation
66+
SCRATCH_N = zeros(size(X, 1))
6267
(g, θ) -> begin
6368
# cache contains the residuals (Xθ-y)
64-
r = SCRATCH_N[]
69+
r = SCRATCH_N
6570
get_residuals!(r, X, θ, y) # -- r = Xθ-y
6671
apply_Xt!(g, X, r)
6772
g .+= λ .* θ

‎src/glr/d_logistic.jl‎

Lines changed: 50 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -10,18 +10,23 @@
1010
# ---------------------------------------------------------
1111

1212
function fgh!(glr::GLR{LogisticLoss,<:L2R}, X, y)
13-
J = objective(glr) # GLR objective (loss+penalty)
14-
p = size(X, 2)
15-
λ = getscale(glr.penalty)
13+
J = objective(glr) # GLR objective (loss+penalty)
14+
n, p = size(X)
15+
λ = getscale(glr.penalty)
16+
# scratch allocation
17+
SCRATCH_N = zeros(n)
18+
SCRATCH_N2 = zeros(n)
19+
SCRATCH_N3 = zeros(n)
1620
if glr.fit_intercept
21+
SCRATCH_P = zeros(p)
1722
(f, g, H, θ) -> begin
18-
Xθ = SCRATCH_N[]
23+
Xθ = SCRATCH_N
1924
apply_X!(Xθ, X, θ) # -- Xθ = apply_X(X, θ)
2025
# precompute σ(yXθ) use -σ(-x) = (σ(x)-1)
21-
w = SCRATCH_N2[]
26+
w = SCRATCH_N2
2227
w .= σ.(Xθ .* y) # -- w = σ.(Xθ .* y)
2328
g === nothing || begin
24-
t = SCRATCH_N3[]
29+
t = SCRATCH_N3
2530
t .= y .* (w .- 1.0) # -- t = y .* (w .- 1.0)
2631
apply_Xt!(g, X, t) # -- g = X't
2732
g .+= λ .* θ
@@ -33,7 +38,7 @@ function fgh!(glr::GLR{LogisticLoss,<:L2R}, X, y)
3338
# probably not really worth it
3439
ΛX = w .* X # !! big allocs
3540
mul!(view(H, 1:p, 1:p), X', ΛX) # -- H[1:p,1:p] = X'ΛX
36-
ΛXt1 = view(SCRATCH_P[], 1:p)
41+
ΛXt1 = SCRATCH_P
3742
copyto!(ΛXt1, sum(ΛX, dims=1)) # -- (ΛX)'1
3843
@inbounds for i = 1:p
3944
H[i, end] = H[end, i] = ΛXt1[i] # -- H[:,p+1] = (ΛX)'1
@@ -47,12 +52,12 @@ function fgh!(glr::GLR{LogisticLoss,<:L2R}, X, y)
4752
# see comments above, same computations just no additional things for
4853
# fit_intercept
4954
(f, g, H, θ) -> begin
50-
Xθ = SCRATCH_N[]
55+
Xθ = SCRATCH_N
5156
apply_X!(Xθ, X, θ)
52-
w = SCRATCH_N2[]
57+
w = SCRATCH_N2
5358
w .= σ.(y .* Xθ)
5459
g === nothing || begin
55-
t = SCRATCH_N3[]
60+
t = SCRATCH_N3
5661
t .= y .* (w .- 1.0)
5762
apply_Xt!(g, X, t)
5863
g .+= λ .* θ
@@ -67,26 +72,30 @@ function fgh!(glr::GLR{LogisticLoss,<:L2R}, X, y)
6772
end
6873

6974
function Hv!(glr::GLR{LogisticLoss,<:L2R}, X, y)
70-
p = size(X, 2)
71-
λ = getscale(glr.penalty)
75+
n, p = size(X)
76+
λ = getscale(glr.penalty)
77+
# scratch allocation
78+
SCRATCH_N = zeros(n)
79+
SCRATCH_N2 = zeros(n)
7280
if glr.fit_intercept
81+
SCRATCH_P = zeros(p)
7382
# H = [X 1]'Λ[X 1] + λ I
7483
# rows a 1:p = [X'ΛX + λI | X'Λ1]
7584
# row e end = [1'ΛX | sum(a)+λ]
7685
(Hv, θ, v) -> begin
77-
Xθ = SCRATCH_N[]
86+
Xθ = SCRATCH_N
7887
apply_X!(Xθ, X, θ) # -- Xθ = apply_X(X, θ)
79-
w = SCRATCH_N2[]
88+
w = SCRATCH_N2
8089
w .= σ.(Xθ .* y) # -- w = σ.(Xθ .* y)
8190
# view on the first p rows
8291
a = 1:p
8392
Hvₐ = view(Hv, a)
8493
vₐ = view(v, a)
85-
XtΛ1 = view(SCRATCH_P[], 1:p)
94+
XtΛ1 = view(SCRATCH_P, 1:p)
8695
mul!(XtΛ1, X', w) # -- X'Λ1; O(np)
8796
vₑ = v[end]
8897
# update for the first p rows -- (X'X + λI)v[1:p] + (X'1)v[end]
89-
Xvₐ = SCRATCH_N[]
98+
Xvₐ = SCRATCH_N
9099
mul!(Xvₐ, X, vₐ)
91100
Xvₐ .*= w # -- ΛXvₐ
92101
mul!(Hvₐ, X', Xvₐ) # -- (X'ΛX)vₐ
@@ -97,13 +106,13 @@ function Hv!(glr::GLR{LogisticLoss,<:L2R}, X, y)
97106
end
98107
else
99108
(Hv, θ, v) -> begin
100-
Xθ = SCRATCH_N[]
109+
Xθ = SCRATCH_N
101110
apply_X!(Xθ, X, θ)
102-
w = SCRATCH_N2[]
111+
w = SCRATCH_N2
103112
w .= σ.(Xθ .* y) # -- σ(yXθ)
104-
Xv = SCRATCH_N3[]
113+
Xv = SCRATCH_N
105114
mul!(Xv, X, v)
106-
Xv .*= SCRATCH_N2[] # -- ΛXv
115+
Xv .*= SCRATCH_N2 # -- ΛXv
107116
mul!(Hv, X', Xv) # -- X'ΛXv
108117
Hv .+= λ .* v
109118
end
@@ -144,23 +153,29 @@ function fg!(glr::GLR{MultinomialLoss,<:L2R}, X, y)
144153
n, p = size(X)
145154
c = length(unique(y))
146155
λ = getscale(glr.penalty)
156+
SCRATCH_N = zeros(n)
157+
SCRATCH_NC = zeros(n, c)
158+
SCRATCH_NC2 = zeros(n, c)
159+
SCRATCH_NC3 = zeros(n, c)
160+
SCRATCH_NC4 = zeros(n, c)
161+
SCRATCH_PC = zeros(p+Int(glr.fit_intercept), c)
147162
(f, g, θ) -> begin
148-
P = SCRATCH_NC[]
149-
apply_X!(P, X, θ, c) # O(npc) store n * c
150-
M = SCRATCH_NC2[]
163+
P = SCRATCH_NC
164+
apply_X!(P, X, θ, c, SCRATCH_PC) # O(npc) store n * c
165+
M = SCRATCH_NC2
151166
M .= exp.(P) # O(npc) store n * c
152167
g === nothing || begin
153-
ΛM = SCRATCH_NC3[]
168+
ΛM = SCRATCH_NC3
154169
ΛM .= M ./ sum(M, dims=2) # O(nc) store n * c
155-
Q = SCRATCH_NC4[]
170+
Q = SCRATCH_NC4
156171
@inbounds for i = 1:n, j=1:c
157172
Q[i, j] = ifelse(y[i] == j, 1.0, 0.0)
158173
end
159174
∑ΛM = sum(ΛM, dims=1)
160175
∑Q = sum(Q, dims=1)
161176
R = ΛM
162177
R .-= Q
163-
G = SCRATCH_PC[]
178+
G = SCRATCH_PC
164179
if glr.fit_intercept
165180
mul!(view(G, 1:p, :), X', R)
166181
@inbounds for k in 1:c
@@ -179,11 +194,11 @@ function fg!(glr::GLR{MultinomialLoss,<:L2R}, X, y)
179194
# ms = maximum(P, dims=2)
180195
# ss = sum(M ./ exp.(ms), dims=2)
181196
ms = maximum(P, dims=2)
182-
ems = SCRATCH_N[]
197+
ems = SCRATCH_N
183198
@inbounds for i in 1:n
184199
ems[i] = exp(ms[i])
185200
end
186-
ΛM = SCRATCH_NC2[] # note that _NC is already linked to P
201+
ΛM = SCRATCH_NC2 # note that _NC is already linked to P
187202
ΛM .= M ./ ems
188203
ss = sum(ΛM, dims=2)
189204
t = 0.0
@@ -207,13 +222,13 @@ function Hv!(glr::GLR{MultinomialLoss,<:L2R}, X, y)
207222
# allocate less but is likely slower; maybe in the future we could have a
208223
# keyword indicating which one the user wants to use.
209224
(Hv, θ, v) -> begin
210-
P = apply_X(X, θ, c) # P_ik = <x_i, θ_k> // dims n * c; O(npc)
211-
Q = apply_X(X, v, c) # Q_ik = <x_i, v_k> // dims n * c; O(npc)
212-
M = exp.(P) # M_ik = exp<x_i, w_k> // dims n * c;
213-
MQ = M .* Q # // dims n * c; O(nc)
214-
ρ = 1 ./ sum(M, dims=2) # ρ_i = 1/Z_i = 1/∑_k exp<x_i, w_k>
215-
κ = sum(MQ, dims=2) # κ_i = ∑_k exp<x_i, w_k><x_i, v_k>
216-
γ = κ .* ρ.^2 # γ_i = κ_i / Z_i^2
225+
P = apply_X(X, θ, c) # P_ik = <x_i, θ_k> // dims n * c; O(npc)
226+
Q = apply_X(X, v, c) # Q_ik = <x_i, v_k> // dims n * c; O(npc)
227+
M = exp.(P) # M_ik = exp<x_i, w_k> // dims n * c;
228+
MQ = M .* Q # // dims n * c; O(nc)
229+
ρ = 1 ./ sum(M, dims=2) # ρ_i = 1/Z_i = 1/∑_k exp<x_i, w_k>
230+
κ = sum(MQ, dims=2) # κ_i = ∑_k exp<x_i, w_k><x_i, v_k>
231+
γ = κ .* ρ.^2 # γ_i = κ_i / Z_i^2
217232
# computation of Hv
218233
U = (ρ .* MQ) .- (γ .* M) # // dims n * c; O(nc)
219234
Hv_mat = X' * U # // dims n * c; O(npc)

0 commit comments

Comments
 (0)