Skip to content

Commit

Permalink
simplify returns for speedup (szcf-weiya/paperMonoDecomp#26)
Browse files Browse the repository at this point in the history
  • Loading branch information
@szcf-weiya
szcf-weiya committed Oct 8, 2023
1 parent b1b5bbd commit c7b44ec
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 20 deletions.
33 changes: 17 additions & 16 deletions src/mono_decomp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,21 +232,22 @@ function mono_decomp_cs(x::AbstractVector{T}, y::AbstractVector{T};
workspace.H = construct_H(J)
end
if s_is_μ
γhat = _optim(y, workspace, [s])[:]
_γhat = _optim(y, workspace, [s])[:]
else
γhat = _optim(y, workspace.J, workspace.B, s, workspace.H)
_γhat = _optim(y, workspace.J, workspace.B, s, workspace.H)
end
γup = γhat[1:J]
γdown = γhat[J+1:2J]
yhat = B * (γup .+ γdown)
γup = _γhat[1:J]
γdown = _γhat[J+1:2J]
γhat = γup .+ γdown
yhat = B * γhat
return MonoDecomp(γup, γdown, γhat, yhat, 0.0, s, workspace)
end

MDCS = mono_decomp_cs

function cubic_spline(x::AbstractVector{T}, y::AbstractVector{T}, xnew::AbstractVector{T}; J = 4) where T <: AbstractFloat
spl = MonotoneSplines.bspline(x, y, J)
return MonotoneSplines.predict(spl, x), MonotoneSplines.predict(spl, xnew), spl
return MonotoneSplines.predict(spl, x), MonotoneSplines.predict(spl, xnew)#, spl
end

mono_decomp_cs(J) = (x, y; kw...) -> mono_decomp_cs(x, y; J = J, kw...)
Expand Down Expand Up @@ -312,9 +313,9 @@ function cv_mono_decomp_cs(x::AbstractVector{T}, y::AbstractVector{T}, xnew::Abs
# workspace = nothing
@assert s_is_μ
γhats = _optim(y[train_idx], workspace, ss)
ynews = predict(workspace, x[test_idx], γhats[1:J, :] + γhats[J+1:2J, :])
ynews = predict(workspace, x[test_idx], γhats[1:J, :] .+ γhats[J+1:2J, :])
for (i, s) in enumerate(ss)
errs[k, j, i] = norm(ynews[:, i] - y[test_idx])^2 / length(test_idx)
errs[k, j, i] = norm(ynews[:, i] .- y[test_idx])^2 / length(test_idx)
end
end
end
Expand All @@ -334,7 +335,7 @@ function cv_mono_decomp_cs(x::AbstractVector{T}, y::AbstractVector{T}, xnew::Abs
savefig(cvplot(μerr, σerr, 1.0 .* Js, ss, lbl = ["J", "μ"], nfold = nfold, ind0 = ind), figname)
end
D = mono_decomp_cs(Jopt)(x, y; s = sopt)
return D, ss[argmin(μerr)[2]], ss, μerr, σerr # use for recording without 1se
return D, ss[argmin(μerr)[2]], μerr, σerr # use for recording without 1se
end

"""
Expand All @@ -350,7 +351,7 @@ function cv_mono_decomp_cs(x::AbstractVector{T}, y::AbstractVector{T}; ss = 10.0
x0 = x,
Js = 4:50,
one_se_rule = false,
one_se_rule_pre = false) where T <: AbstractFloat
one_se_rule_pre = false)::Tuple{MonoDecomp{T}, T, Array{T}, Array{T}} where T <: AbstractFloat
if fixJ
if length(Js) == 1
J = Js[1]
Expand All @@ -359,7 +360,7 @@ function cv_mono_decomp_cs(x::AbstractVector{T}, y::AbstractVector{T}; ss = 10.0
J, yhat, yhatnew = cv_cubic_spline(x, y, x0, one_se_rule = one_se_rule_pre, nfold = nfold_pre, Js = Js,
figname = isnothing(figname) ? figname : figname[1:end-4] * "_bspl.png")
end
return cv_mono_decomp_cs(x, y, x0, Js = J:J, ss = ss, figname = figname, nfold = nfold, one_se_rule = one_se_rule)..., yhat, yhatnew
return cv_mono_decomp_cs(x, y, x0, Js = J:J, ss = ss, figname = figname, nfold = nfold, one_se_rule = one_se_rule)
else
return cv_mono_decomp_cs(x, y, x0, ss = ss, Js = Js, figname = figname, nfold = nfold, one_se_rule = one_se_rule)
end
Expand Down Expand Up @@ -1066,15 +1067,15 @@ Monotone decomposition with smoothing splines.
function mono_decomp_ss(workspace::WorkSpaceSS, x::AbstractVector{T}, y::AbstractVector{T}, λ::AbstractFloat, s::AbstractFloat; s_is_μ = true, prop_nknots = 1.0, strict = false)::MonoDecomp{T} where T <: AbstractFloat
build_model!(workspace, x, prop_nknots = prop_nknots)
if s_is_μ
γhat = _optim(y, workspace.J, workspace.B, nothing, workspace.H, L = workspace.L, t = nothing, λ = λ, μ = s, strict = strict)
_γhat = _optim(y, workspace.J, workspace.B, nothing, workspace.H, L = workspace.L, t = nothing, λ = λ, μ = s, strict = strict)
else
γhat = _optim(y, workspace.J, workspace.B, s, workspace.H, L = workspace.L, t = nothing, λ = λ, strict = strict)
_γhat = _optim(y, workspace.J, workspace.B, s, workspace.H, L = workspace.L, t = nothing, λ = λ, strict = strict)
end
# calculate properties of monotone decomposition
J = workspace.J
γup = γhat[1:J]
γdown = γhat[J+1:2J]
γhat = γup + γdown
γup = _γhat[1:J]
γdown = _γhat[J+1:2J]
γhat = γup .+ γdown
yhat = workspace.B * γhat
return MonoDecomp(γup, γdown, γhat, yhat, λ, s, workspace)
end
Expand Down
9 changes: 5 additions & 4 deletions src/mono_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,13 +268,13 @@ function mono_test_bootstrap_cs(x::AbstractVector{T}, y::AbstractVector{T}; nrep
figname = nothing,
nfold = 10, nfold_pre = 10,
kw...)::Tuple{T, MonoDecomp{T}} where T <: Real
D, μ, μs = cv_mono_decomp_cs(x, y, ss = μs, one_se_rule = one_se_rule, fixJ = fixJ, Js = Js, one_se_rule_pre = one_se_rule_pre, figname = figname, nfold = nfold, nfold_pre = nfold_pre)
D, μ = cv_mono_decomp_cs(x, y, ss = μs, one_se_rule = one_se_rule, fixJ = fixJ, Js = Js, one_se_rule_pre = one_se_rule_pre, figname = figname, nfold = nfold, nfold_pre = nfold_pre)
@debug D.γdown
J = D.workspace.J
# scatter(x, y)
# plot!(x, res.yhat)
c = mean(D.yhat) / 2
error = y - D.yhat
error = y .- D.yhat
@debug maximum(max.(error))
# σ = std(err)
tobs = var(D.γdown)
Expand All @@ -301,7 +301,8 @@ function mono_test_bootstrap_sup(x::AbstractVector{T}, y::AbstractVector{T};
nblock = 10,
kw...) where T <: AbstractFloat
n = length(y)
D1, μ0, μs0 = cv_mono_decomp_cs(x, y, ss = 10.0 .^ (-6:0.5:6), one_se_rule = true, fixJ = fixJ, nfold = nfold)
μs0 = 10.0 .^ (-6:0.5:6)
D1, μ0 = cv_mono_decomp_cs(x, y, ss = μs0, one_se_rule = true, fixJ = fixJ, nfold = nfold)
μ1 = D1.μ
J = D1.workspace.J
# μ0 < μ1
Expand Down Expand Up @@ -383,7 +384,7 @@ function construct_bootstrap_y(y::AbstractVector{T}, e::AbstractVector{T}, B::Ab
idx = sample(1:n, n)
ei = e[idx] .- mean(e[idx]) .+ mean(e)
end
yi = B * γ .+ c + ei
yi = B * γ .+ c .+ ei
if debias_mean_yi
yi = yi .- mean(yi) .+ mean(y)
end
Expand Down

0 comments on commit c7b44ec

Please sign in to comment.