diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index fae822e..637dbd7 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -51,22 +51,15 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v3 - # ProximalSolvers, CovarianceSelection and the ProximalBase that - # CovarianceSelection needs (SliceSum, FiberSum) are not registered yet: all - # from their masters. Julia 1.11 and later take them from [sources] in - # Project.toml; Julia 1.10 ignores [sources], and resolves each Pkg.add with - # the whole project, so they are added in one call. Remove the step and - # [sources] once they are registered. - - name: Add the unregistered dependencies from GitHub on Julia 1.10 - if: matrix.version == 'lts' - shell: julia --color=yes --project=. {0} - run: | - using Pkg - Pkg.add([PackageSpec(url="https://github.com/mlakolar/ProximalBase.jl", rev="master"), - PackageSpec(url="https://github.com/mlakolar/ProximalSolvers.jl", rev="master"), - PackageSpec(url="https://github.com/mlakolar/CovarianceSelection.jl", rev="master")]) + - name: Check naming conventions + run: ${{ runner.os == 'Windows' && 'python' || 'python3' }} scripts/check_naming.py + - name: Resolve dependencies from master checkouts + run: julia --color=yes scripts/instantiate.jl - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 + - name: Test against dependency master checkouts + run: julia --color=yes scripts/test_master.jl + env: + JULIA_TEST_COVERAGE: "true" - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v7 with: diff --git a/.gitignore b/.gitignore index 6732b5d..be00f20 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,9 @@ /exper/sim*/*.pdf /benchmark/compare/*/Manifest*.toml /benchmark/compare/results/ + +# Generated master dependency test environments +.master-test-env/ + +# Local Julia-minor environments are regenerated from master checkouts +/Manifest-v*.toml diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..74edb25 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,21 @@ +# Development conventions + +Types and modules use UpperCamelCase. Descriptive functions, fields and keyword +arguments use snake_case; mutating functions end in `!`. Retain standard numerical +abbreviations, mathematical identifiers and names required by upstream interfaces. +Run `python3 scripts/check_naming.py` to check new declarations. + +Workspace dependencies use sibling master checkouts, including transitive +dependencies. Missing checkouts are cloned from master. Existing checkouts are +never reset or silently updated. Package version numbering and registry changes +are deferred. + +Run the complete tests from the repository root: + +```sh +julia scripts/test_master.jl +``` + +The script prepares an isolated master-based environment with test extras on +Julia 1.10 and newer, then runs the suite without changing package dependencies. +For the simulation environment, use `julia scripts/instantiate.jl exper`. diff --git a/Project.toml b/Project.toml index 1924c69..7636180 100644 --- a/Project.toml +++ b/Project.toml @@ -15,9 +15,9 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [sources] -CovarianceSelection = {url = "https://github.com/mlakolar/CovarianceSelection.jl", rev = "master"} -ProximalBase = {url = "https://github.com/mlakolar/ProximalBase.jl", rev = "master"} -ProximalSolvers = {url = "https://github.com/mlakolar/ProximalSolvers.jl", rev = "master"} +CovarianceSelection = {path = "../CovarianceSelection.jl"} +ProximalBase = {path = "../ProximalBase.jl"} +ProximalSolvers = {path = "../ProximalOPT.jl"} [compat] Aqua = "0.8" diff --git a/README.md b/README.md index 44a2bd5..3052805 100644 --- a/README.md +++ b/README.md @@ -18,8 +18,8 @@ independent Gaussian samples. `ReducedOracleNormal`, and with bootstrap standard errors `SymmetricOracleBoot`, `AsymmetricOracleBoot` * Simultaneous confidence intervals for all entries by the bootstrap: - `bootstrap`, `simulCI`, `simulCIstudentized` -* Support recovery: `supportEstimate` with `ANTSupport`, `BootStdSupport` or + `bootstrap`, `simul_ci`, `simul_ci_studentized` +* Support recovery: `support_estimate` with `ANTSupport`, `BootStdSupport` or `DTraceValidationSupport` DiffPrecTest is one of five Julia packages for high-dimensional statistics, @@ -27,7 +27,7 @@ three for estimation and inference and two for the optimization beneath them: | Package | | |---|---| -| [CoordinateDescent.jl](https://github.com/mlakolar/CoordinateDescent.jl) | the lasso and its paths; the group, square-root, scaled and heteroscedastic lasso; varying-coefficient models | +| [HighDimensionalModels.jl](https://github.com/mlakolar/HighDimensionalModels.jl) | the lasso and its paths; the group, square-root, scaled and heteroscedastic lasso; varying-coefficient models | | [CovarianceSelection.jl](https://github.com/mlakolar/CovarianceSelection.jl) | Gaussian graphical models: the graphical lasso, CLIME and ACLIME, precision matrices by the scaled lasso, the fused graphical lasso, and the difference of two precision matrices | | **DiffPrecTest.jl** (this package) | confidence intervals and tests for the difference of two precision matrices | | [ProximalSolvers.jl](https://github.com/mlakolar/ProximalSolvers.jl) | the solvers: coordinate descent with penalty paths, screening and active sets; proximal gradient methods, proximal Newton and ADMM | @@ -53,8 +53,8 @@ res.p .+ [-1.96, 1.96] .* res.std # simultaneous 95% confidence intervals for the entries of the lower triangle of Δ; # the off-diagonal entries are multiplied by √2 -straps, = bootstrap(X, Y; bootSamples=500) -CI = simulCIstudentized(straps, 0.95) +straps, = bootstrap(X, Y; boot_samples=500) +CI = simul_ci_studentized(straps, 0.95) ``` `Distributions` also exports `estimate`; when both packages are loaded, call @@ -68,7 +68,40 @@ coordinate descent of [ProximalSolvers.jl](https://github.com/mlakolar/ProximalS on losses of [ProximalBase.jl](https://github.com/mlakolar/ProximalBase.jl): the direct difference loss of [CovarianceSelection.jl](https://github.com/mlakolar/CovarianceSelection.jl) and the inverse Hessian losses are its `KronQuadratic`. The functions of the -two stages (`diffEstimation`, `invQSymHessian`, `invAsymHessianEstimation`, -`reducedDiffEstimation`, `invHessianReduced`) take the solver options as a +two stages (`diff_estimation`, `inv_q_sym_hessian`, `inv_asym_hessian_estimation`, +`reduced_diff_estimation`, `inv_hessian_reduced`) take the solver options as a ProximalSolvers `CDSolver`; so are the node-wise lasso regressions of `SeparateNormal`. + +## Naming migration + +Types retain `UpperCamelCase`; descriptive functions, fields and keyword +arguments use `snake_case`. Mathematical identifiers remain available. + +| Previous spelling | Canonical spelling | +|---|---| +| `diffEstimation` | `diff_estimation` | +| `invAsymHessianEstimation`, `invQSymHessian`, `invHessianReduced` | `inv_asym_hessian_estimation`, `inv_q_sym_hessian`, `inv_hessian_reduced` | +| `supportEstimate` | `support_estimate` | +| `simulCI`, `simulCIstudentized` | `simul_ci`, `simul_ci_studentized` | +| `computeSimulationResult` | `compute_simulation_result` | +| `bootSamples`, `estimSupport` | `boot_samples`, `estim_support` | +| `SimulationResult.lenCoverage` | `SimulationResult.len_coverage` | + +Previous public spellings are deprecated aliases for a transition period. +Supplying both spellings of a keyword raises an error. Package versions are +unchanged. Active environments use sibling master checkouts of +CovarianceSelection, ProximalBase and ProximalSolvers through local sources. +The pinned `benchmark/compare/old` environment remains a historical baseline. + +## Development tests + +From the repository root, run: + +```sh +julia scripts/test_master.jl +``` + +The script prepares an isolated test environment using dependency master +checkouts and test extras on Julia 1.10 and newer. See +[CONTRIBUTING.md](CONTRIBUTING.md) for the naming policy and environment setup. diff --git a/benchmark/compare/new/Project.toml b/benchmark/compare/new/Project.toml index 05ef43b..f91d80f 100644 --- a/benchmark/compare/new/Project.toml +++ b/benchmark/compare/new/Project.toml @@ -11,10 +11,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [sources] -CovarianceSelection = {url = "https://github.com/mlakolar/CovarianceSelection.jl", rev = "master"} +CovarianceSelection = {path = "../../../../CovarianceSelection.jl"} DiffPrecTest = {path = "../../.."} -ProximalBase = {url = "https://github.com/mlakolar/ProximalBase.jl", rev = "master"} -ProximalSolvers = {url = "https://github.com/mlakolar/ProximalSolvers.jl", rev = "master"} +ProximalBase = {path = "../../../../ProximalBase.jl"} +ProximalSolvers = {path = "../../../../ProximalOPT.jl"} [compat] CovarianceSelection = "0.1" diff --git a/exper/Project.toml b/exper/Project.toml index e45dd33..9f10392 100644 --- a/exper/Project.toml +++ b/exper/Project.toml @@ -17,10 +17,10 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" # ProximalSolvers is listed for its [sources] entry, until it is registered [sources] -CovarianceSelection = {url = "https://github.com/mlakolar/CovarianceSelection.jl", rev = "master"} +CovarianceSelection = {path = "../../CovarianceSelection.jl"} DiffPrecTest = {path = ".."} -ProximalBase = {url = "https://github.com/mlakolar/ProximalBase.jl", rev = "master"} -ProximalSolvers = {url = "https://github.com/mlakolar/ProximalSolvers.jl", rev = "master"} +ProximalBase = {path = "../../ProximalBase.jl"} +ProximalSolvers = {path = "../../ProximalOPT.jl"} [compat] CSV = "1" diff --git a/exper/README.md b/exper/README.md index eb08423..f09c1b6 100644 --- a/exper/README.md +++ b/exper/README.md @@ -20,17 +20,17 @@ the true support (`O-...`). ## Setup -With Julia 1.11 or later: +From the repository root, prepare the experiment environment on Julia 1.10 +or newer: ```sh -julia --project=. -e 'using Pkg; Pkg.instantiate()' +julia scripts/instantiate.jl exper ``` -`Project.toml` takes DiffPrecTest from this repository, and CovarianceSelection, -ProximalSolvers and ProximalBase from their masters on GitHub, until they are -registered. Julia 1.10 ignores these sources; there, run -`using Pkg; Pkg.develop(path=".."); Pkg.add([PackageSpec(url="https://github.com/mlakolar/ProximalBase.jl", rev="master"), PackageSpec(url="https://github.com/mlakolar/ProximalSolvers.jl", rev="master"), PackageSpec(url="https://github.com/mlakolar/CovarianceSelection.jl", rev="master")])` -once instead. +The bootstrap loads DiffPrecTest from this repository and CovarianceSelection, +ProximalSolvers and ProximalBase from sibling master checkouts, including their +transitive workspace dependencies. Existing checkouts are never reset or +silently updated. ## Running diff --git a/exper/sim10/sim10.jl b/exper/sim10/sim10.jl index 864b453..0b441e7 100644 --- a/exper/sim10/sim10.jl +++ b/exper/sim10/sim10.jl @@ -53,7 +53,7 @@ Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -@time boot_res, eS = bootstrap(X, Y) +@time boot_res, e_s = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res e_s diff --git a/exper/sim10/sim10_o.jl b/exper/sim10/sim10_o.jl index f48ded7..6889d6b 100644 --- a/exper/sim10/sim10_o.jl +++ b/exper/sim10/sim10_o.jl @@ -55,7 +55,7 @@ Sy = Symmetric( Y'Y / n ) # oracle supports: Δ = 0, so every entry is refitted on its own eS_init = [Int64[] for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(ip)_$(rep).jld2" boot_res diff --git a/exper/sim10/sim10_o_power.jl b/exper/sim10/sim10_o_power.jl index 64ff590..afd9147 100644 --- a/exper/sim10/sim10_o_power.jl +++ b/exper/sim10/sim10_o_power.jl @@ -83,6 +83,6 @@ Sy = Symmetric( Y'Y / n ) S = [tril_index(p, r, c) for c = 1:p for r = c:p if abs(Δ[r, c]) > 1e-4] eS_init = [copy(S) for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim11/sim11.jl b/exper/sim11/sim11.jl index 043fb16..1cd576f 100644 --- a/exper/sim11/sim11.jl +++ b/exper/sim11/sim11.jl @@ -55,7 +55,7 @@ Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -@time boot_res, eS = bootstrap(X, Y) +@time boot_res, e_s = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res e_s diff --git a/exper/sim11/sim11_o.jl b/exper/sim11/sim11_o.jl index 17e2091..a628699 100644 --- a/exper/sim11/sim11_o.jl +++ b/exper/sim11/sim11_o.jl @@ -57,7 +57,7 @@ Sy = Symmetric( Y'Y / n ) # oracle supports: Δ = 0, so every entry is refitted on its own eS_init = [Int64[] for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(ip)_$(rep).jld2" boot_res diff --git a/exper/sim11/sim11_o_power.jl b/exper/sim11/sim11_o_power.jl index 8acb7ed..6245eea 100644 --- a/exper/sim11/sim11_o_power.jl +++ b/exper/sim11/sim11_o_power.jl @@ -89,6 +89,6 @@ Sy = Symmetric( Y'Y / n ) S = [tril_index(p, r, c) for c = 1:p for r = c:p if abs(Δ[r, c]) > 1e-4] eS_init = [copy(S) for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim12/sim12_o_power.jl b/exper/sim12/sim12_o_power.jl index 241fd95..135cf7d 100644 --- a/exper/sim12/sim12_o_power.jl +++ b/exper/sim12/sim12_o_power.jl @@ -72,6 +72,6 @@ Sy = Symmetric( Y'Y / n ) S = [tril_index(p, r, c) for c = 1:p for r = c:p if abs(Δ[r, c]) > 1e-4] eS_init = [copy(S) for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim5/sim5a.jl b/exper/sim5/sim5a.jl index 0d7d010..16c58ef 100644 --- a/exper/sim5/sim5a.jl +++ b/exper/sim5/sim5a.jl @@ -48,9 +48,9 @@ Sy = Symmetric( cov(Y) ) Sxtest = Symmetric( cov(Xtest) ) Sytest = Symmetric( cov(Ytest) ) -@time eS = DiffPrecTest.__initSupport(Sx, Sy, X, Y) -@time eΔNormal, eΔNormalP, _ = supportEstimate(ANTSupport(), X, Y, τArr; estimSupport=eS) -@time eΔBoot , eΔBootP, _ = supportEstimate(BootStdSupport(), X, Y, τArr; estimSupport=eS) +@time e_s = DiffPrecTest.__init_support(Sx, Sy, X, Y) +@time eΔNormal, eΔNormalP, _ = support_estimate(ANTSupport(), X, Y, τArr; estim_support=e_s) +@time eΔBoot , eΔBootP, _ = support_estimate(BootStdSupport(), X, Y, τArr; estim_support=e_s) # compete method @@ -58,6 +58,6 @@ maxλ = maximum( abs.(Sx - Sy) ) minλ = maxλ * 0.04 Λarr = exp.(range(log(maxλ), log(minλ), length=50)) -eΔDTr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) +eΔDTr, i2, i_inf, loss2arr, loss_infarr = support_estimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) -@save "$(dir)/res_$(rep).jld2" eΔNormal eΔNormalP eΔBoot eΔBootP eS eΔDTr i2 iInf loss2arr lossInfarr +@save "$(dir)/res_$(rep).jld2" eΔNormal eΔNormalP eΔBoot eΔBootP e_s eΔDTr i2 i_inf loss2arr loss_infarr diff --git a/exper/sim5/sim5b.jl b/exper/sim5/sim5b.jl index 5617e58..f79b190 100644 --- a/exper/sim5/sim5b.jl +++ b/exper/sim5/sim5b.jl @@ -48,9 +48,9 @@ Sy = Symmetric( cov(Y) ) Sxtest = Symmetric( cov(Xtest) ) Sytest = Symmetric( cov(Ytest) ) -@time eS = DiffPrecTest.__initSupport(Sx, Sy, X, Y) -@time eΔNormal, eΔNormalP, _ = supportEstimate(ANTSupport(), X, Y, τArr; estimSupport=eS) -@time eΔBoot , eΔBootP, _ = supportEstimate(BootStdSupport(), X, Y, τArr; estimSupport=eS) +@time e_s = DiffPrecTest.__init_support(Sx, Sy, X, Y) +@time eΔNormal, eΔNormalP, _ = support_estimate(ANTSupport(), X, Y, τArr; estim_support=e_s) +@time eΔBoot , eΔBootP, _ = support_estimate(BootStdSupport(), X, Y, τArr; estim_support=e_s) # compete method @@ -58,6 +58,6 @@ maxλ = maximum( abs.(Sx - Sy) ) minλ = maxλ * 0.04 Λarr = exp.(range(log(maxλ), log(minλ), length=50)) -eΔDTr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) +eΔDTr, i2, i_inf, loss2arr, loss_infarr = support_estimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) -@save "$(dir)/res_$(rep).jld2" eΔNormal eΔNormalP eΔBoot eΔBootP eS eΔDTr i2 iInf loss2arr lossInfarr +@save "$(dir)/res_$(rep).jld2" eΔNormal eΔNormalP eΔBoot eΔBootP e_s eΔDTr i2 i_inf loss2arr loss_infarr diff --git a/exper/sim5/summarize.jl b/exper/sim5/summarize.jl index 6fe6864..77d86c5 100644 --- a/exper/sim5/summarize.jl +++ b/exper/sim5/summarize.jl @@ -38,8 +38,8 @@ TPR_tr = zeros(NUM_REP, numΛ) dtrace = DataFrame(selection=String[], rep=Int[], tpr=Float64[], fpr=Float64[], precision=Float64[]) for (rep, f) in enumerate(files) - eΔNormal, eΔBoot, eΔDTr, i2, iInf = jldopen(joinpath(dir, f), "r") do file - read(file, "eΔNormal"), read(file, "eΔBoot"), read(file, "eΔDTr"), read(file, "i2"), read(file, "iInf") + eΔNormal, eΔBoot, eΔDTr, i2, i_inf = jldopen(joinpath(dir, f), "r") do file + read(file, "eΔNormal"), read(file, "eΔBoot"), read(file, "eΔDTr"), read(file, "i2"), read(file, "i_inf") end for j=1:numτ @@ -58,7 +58,7 @@ for (rep, f) in enumerate(files) TPR_tr[rep, j] = tpr(cm) end - for (selection, i) in (("frobenius", i2), ("max", iInf)) + for (selection, i) in (("frobenius", i2), ("max", i_inf)) cm = getConfusionMatrix(Δ, eΔDTr[i]) push!(dtrace, (selection, rep, tpr(cm), fpr(cm), precision(cm))) end diff --git a/exper/sim8/collect_our.jl b/exper/sim8/collect_our.jl index 4b02deb..b330461 100644 --- a/exper/sim8/collect_our.jl +++ b/exper/sim8/collect_our.jl @@ -18,11 +18,11 @@ NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) covers(CI) = all(CI[:, 1] .<= 0 .<= CI[:, 2]) function coverage(res::Vector{BootstrapEstimates}, α::Float64=0.95) - mean(covers(simulCI(r, α)) for r in res) + mean(covers(simul_ci(r, α)) for r in res) end function coverage_studentized(res::Vector{BootstrapEstimates}, α::Float64=0.95) - mean(covers(simulCIstudentized(r, α)) for r in res) + mean(covers(simul_ci_studentized(r, α)) for r in res) end # loads folderName_rep.jld2 for the available repetitions diff --git a/exper/sim8/compute_coverage.jl b/exper/sim8/compute_coverage.jl index 38ec28e..063ccf2 100644 --- a/exper/sim8/compute_coverage.jl +++ b/exper/sim8/compute_coverage.jl @@ -14,5 +14,5 @@ res = load(fName, "results") # the intervals cover Δ = 0 when every one of them contains 0 covers(CI) = all(CI[:, 1] .<= 0 .<= CI[:, 2]) -@show mean(covers(simulCI(r, alpha)) for r in res) -@show mean(covers(simulCIstudentized(r, alpha)) for r in res) +@show mean(covers(simul_ci(r, alpha)) for r in res) +@show mean(covers(simul_ci_studentized(r, alpha)) for r in res) diff --git a/exper/sim8/sim8.jl b/exper/sim8/sim8.jl index 2af0c3d..d916e11 100644 --- a/exper/sim8/sim8.jl +++ b/exper/sim8/sim8.jl @@ -54,7 +54,7 @@ Sy = Symmetric( Y'Y / n ) -@time boot_res, eS = bootstrap(X, Y) +@time boot_res, e_s = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res e_s diff --git a/exper/sim8/sim8_o.jl b/exper/sim8/sim8_o.jl index 616d3ec..c0304af 100644 --- a/exper/sim8/sim8_o.jl +++ b/exper/sim8/sim8_o.jl @@ -54,6 +54,6 @@ Sy = Symmetric( Y'Y / n ) # oracle supports: Δ = 0, so every entry is refitted on its own eS_init = [Int64[] for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(ip)_$(rep).jld2" boot_res diff --git a/exper/sim8/sim8_o_power.jl b/exper/sim8/sim8_o_power.jl index 15f19ea..16ba18f 100644 --- a/exper/sim8/sim8_o_power.jl +++ b/exper/sim8/sim8_o_power.jl @@ -83,6 +83,6 @@ Sy = Symmetric( Y'Y / n ) S = [tril_index(p, r, c) for c = 1:p for r = c:p if abs(Δ[r, c]) > 1e-4] eS_init = [copy(S) for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim9/sim9.jl b/exper/sim9/sim9.jl index 81d220e..f756694 100644 --- a/exper/sim9/sim9.jl +++ b/exper/sim9/sim9.jl @@ -51,7 +51,7 @@ Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -@time boot_res, eS = bootstrap(X, Y) +@time boot_res, e_s = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res e_s diff --git a/exper/sim9/sim9_o.jl b/exper/sim9/sim9_o.jl index 7048c27..fbc4710 100644 --- a/exper/sim9/sim9_o.jl +++ b/exper/sim9/sim9_o.jl @@ -54,7 +54,7 @@ Sy = Symmetric( Y'Y / n ) # oracle supports: Δ = 0, so every entry is refitted on its own eS_init = [Int64[] for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(ip)_$(rep).jld2" boot_res diff --git a/exper/sim9/sim9_o_power.jl b/exper/sim9/sim9_o_power.jl index 7dcdbf3..2a9e882 100644 --- a/exper/sim9/sim9_o_power.jl +++ b/exper/sim9/sim9_o_power.jl @@ -81,6 +81,6 @@ Sy = Symmetric( Y'Y / n ) S = [tril_index(p, r, c) for c = 1:p for r = c:p if abs(Δ[r, c]) > 1e-4] eS_init = [copy(S) for j=1:div((p + 1)*p, 2)] -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, e_s = bootstrap(X, Y; estim_support=eS_init) @save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/summary.jl b/exper/summary.jl index 0a482d9..98498f5 100644 --- a/exper/summary.jl +++ b/exper/summary.jl @@ -14,8 +14,8 @@ newSummary() = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[ function summarizeSetting!(df, file, p, (row, col), trueValue) res = load(file, "results") for j in eachindex(methodArr) - sr = computeSimulationResult([res[j, i] for i=1:size(res, 2)], trueValue) - push!(df, (p, row, col, trueValue, methodArr[j], sr.bias, sr.coverage, sr.lenCoverage)) + sr = compute_simulation_result([res[j, i] for i=1:size(res, 2)], trueValue) + push!(df, (p, row, col, trueValue, methodArr[j], sr.bias, sr.coverage, sr.len_coverage)) end df end diff --git a/scripts/check_naming.py b/scripts/check_naming.py new file mode 100644 index 0000000..e4cc472 --- /dev/null +++ b/scripts/check_naming.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +"""Check new Julia declarations; compatibility files and mathematical names are explicit exceptions.""" +import argparse +from pathlib import Path +import re + +# Matrix caches and scientific notation retain their conventional identifiers. +MATHEMATICAL = {'A','B','C','D','G','H','I','J','K','L','M','N','P','Q','R','S','T','U','V','W','X','Y','Z', + 'Ax','At','Bt','Xt','Yt','XtY','XtX','sUtY','Ws','SX','SY','Xn','Yn','SXn','SYn', + 'Xa','Ya','SXa','SYa','Xr','Yr','posA','posR','T0','F','Sx','Sy','Sxy','Sxx','Syy','AD','Vt','AU','D_calls','Dt_calls','Dx','Dxp','Dty','Dtyp','Bpath'} +FUNCTION = re.compile(r'^\s*(?:function\s+)([\w.!]+)\s*\(|^\s*([\w.!]+)\([^\n]*\)\s*(?:where\s+[^=]+)?=(?!=)') +FIELD = re.compile(r'^\s*(\w+)\s*::') +LEGACY_KEYWORDS = {'printEvery','normD','maxK','requested_maxK','epsTol','maxIter','checkEvery','bootSamples','estimSupport'} +FACTORIES = {'IndBallLinf','Constrained','_AtomIterate','CDLeastSquaresLoss','CDWeightedLSLoss','CDSqrtLassoLoss','CDSqrtQuadraticLoss','CDQuadraticLoss'} +TYPE = re.compile(r'^\s*(?:mutable\s+)?(?:struct|abstract\s+type|primitive\s+type)\s+(\w+)') + +def check(root): + errors=[] + types={m.group(1) for path in (root/'src').rglob('*.jl') + for line in path.read_text(encoding="utf-8").splitlines() if (m := TYPE.match(line))} + types |= {m.group(1) for path in (root/'src').rglob('*.jl') for m in re.finditer(r'^\s*const\s+(\w+)\s*=',path.read_text(encoding="utf-8"),re.M)} + for path in sorted((root/'src').rglob('*.jl')): + if 'deprecated' in path.stem: continue + in_doc=False; in_struct=False + lines=path.read_text(encoding="utf-8").splitlines() + for number,line in enumerate(lines,1): + if line.count('"""') % 2: in_doc=not in_doc; continue + if in_doc: continue + line=line.split('#',1)[0] + typ=TYPE.match(line) + if typ: + name=typ.group(1) + if not name.lstrip('_')[0].isupper():errors.append((path,number,'type',name)) + in_struct='struct' in line + continue + if re.match(r'^\s*end\b',line):in_struct=False + function=FUNCTION.match(line) + if function: + name=function.group(1) or function.group(2) + # Qualified names implement external interfaces; constructors must name an actual type. + if '.' not in name and name not in types and name not in FACTORIES and name not in MATHEMATICAL and re.search('[A-Z]',name): + errors.append((path,number,'function',name)) + if ';' in line or line.lstrip().startswith('function '): + signature=line[line.find('(')+1:] + depth=1 + signature.count('(')-signature.count(')') + offset=number + while depth>0 and offset ("ProximalBase.jl", "ProximalBase.jl"), + "ProximalSolvers" => ("ProximalOPT.jl", "ProximalSolvers.jl"), + "HighDimensionalModels" => ("HighDimensionalModels.jl", "HighDimensionalModels.jl"), + "CovarianceSelection" => ("CovarianceSelection.jl", "CovarianceSelection.jl"), + "DiffPrecTest" => ("DiffPrecTest.jl", "DiffPrecTest.jl"), + "KernelSmoothing" => ("KernelSmoothing.jl", "KernelSmoothing.jl"), + "FunctionalDataAnalysis" => ("FunctionalDataAnalysis.jl", "FunctionalDataAnalysis.jl"), + "FunctionalGraphicalModels" => ("FunctionalGraphicalModels.jl", "FunctionalGraphicalModels.jl"), +) + +function instantiate_master(environment::AbstractString; + workspace::Union{Nothing,AbstractString}=nothing, + owner::AbstractString=environment) + environment = realpath(environment) + owner = realpath(owner) + workspace = workspace === nothing ? dirname(owner) : realpath(workspace) + owner_project = TOML.parsefile(joinpath(owner, "Project.toml")) + owner_name = get(owner_project, "name", "") + specs = Dict{String,Pkg.PackageSpec}() + visiting = Set{String}() + function visit(project_path) + project = TOML.parsefile(project_path) + deps = merge(get(project, "deps", Dict()), get(project, "extras", Dict())) + for name in sort!(collect(keys(deps))) + haskey(WORKSPACE_PACKAGES, name) || continue + name in visiting && continue + push!(visiting, name) + folder, repository = WORKSPACE_PACKAGES[name] + is_owner = name == owner_name + checkout = is_owner ? owner : joinpath(workspace, folder) + if !isdir(checkout) + url = "https://github.com/mlakolar/$repository" + run(`git clone --branch master --single-branch $url $checkout`) + end + branch = strip(read(`git -C $checkout branch --show-current`, String)) + (is_owner || branch == "master") || error("$name must use its master checkout: $checkout is on $branch") + dependency = TOML.parsefile(joinpath(checkout, "Project.toml")) + dependency["name"] == name || error("wrong package name at $checkout") + dependency["uuid"] == deps[name] || error("wrong UUID for $name at $checkout") + specs[name] = Pkg.PackageSpec(; path=checkout) + visit(joinpath(checkout, "Project.toml")) + end + end + project_path = joinpath(environment, "Project.toml") + original = read(project_path, String) + visit(project_path) + Pkg.activate(environment) + try + if isempty(specs) + Pkg.resolve() + else + Pkg.develop([specs[name] for name in sort!(collect(keys(specs)))]) + end + Pkg.instantiate() + finally + # Transitive master dependencies must not become stale direct dependencies. + write(project_path, original) + end + # Keep extra-only master paths in the manifest for Julia 1.10 test sandboxes. + # Resolving again after restoring Project.toml would prune those entries. + return nothing +end + +if abspath(PROGRAM_FILE) == @__FILE__ + isempty(ARGS) && error("usage: julia instantiate_master.jl PACKAGE_DIRECTORY [WORKSPACE_DIRECTORY]") + instantiate_master(ARGS[1]; workspace=length(ARGS) > 1 ? ARGS[2] : dirname(realpath(ARGS[1]))) +end diff --git a/scripts/test_master.jl b/scripts/test_master.jl new file mode 100644 index 0000000..464cd38 --- /dev/null +++ b/scripts/test_master.jl @@ -0,0 +1,53 @@ +# A dedicated test environment keeps extra-only master dependencies reachable on +# Julia 1.10 and leaves the package Project.toml unchanged for Aqua checks. +using Pkg, TOML +include(joinpath(@__DIR__, "master_dependencies.jl")) +package_root = dirname(@__DIR__) +project_path = joinpath(package_root, "Project.toml") +original_project = read(project_path, String) +project = TOML.parse(original_project) +deps = copy(get(project, "deps", Dict{String,Any}())) +extras = merge(get(project, "weakdeps", Dict()), get(project, "extras", Dict())) +for name in get(get(project, "targets", Dict()), "test", String[]) + deps[name] = extras[name] +end +deps[project["name"]] = project["uuid"] +compat = Dict(name => value for (name, value) in get(project, "compat", Dict()) + if name == "julia" || haskey(deps, name)) +sources = Dict{String,Any}() +for name in keys(deps) + if name == project["name"] + sources[name] = Dict("path" => package_root) + elseif haskey(WORKSPACE_PACKAGES, name) + folder, _ = WORKSPACE_PACKAGES[name] + sources[name] = Dict("path" => joinpath(dirname(package_root), folder)) + elseif haskey(get(project, "sources", Dict()), name) + source = copy(project["sources"][name]) + haskey(source, "path") && (source["path"] = abspath(joinpath(package_root, source["path"]))) + sources[name] = source + end +end +environment = joinpath(package_root, ".master-test-env", "v$(VERSION.major).$(VERSION.minor)") +mkpath(environment) +open(joinpath(environment, "Project.toml"), "w") do io + TOML.print(io, Dict("deps" => deps, "compat" => compat, "sources" => sources)) +end +instantiate_master(environment; workspace=dirname(package_root), owner=package_root) +read(project_path, String) == original_project || error("test setup changed package Project.toml") +installed = Pkg.dependencies() +for (name, uuid) in deps + if name == project["name"] || haskey(WORKSPACE_PACKAGES, name) + expected = name == project["name"] ? package_root : + joinpath(dirname(package_root), WORKSPACE_PACKAGES[name][1]) + realpath(installed[Base.UUID(uuid)].source) == realpath(expected) || + error("test environment did not retain the master checkout of $name") + end +end +println("Master test environment: ", environment) +if !("--setup-only" in ARGS) + depwarn = get(ENV, "JULIA_TEST_DEPWARN", "yes") + depwarn in ("yes", "no", "error") || error("invalid JULIA_TEST_DEPWARN") + coverage = get(ENV, "JULIA_TEST_COVERAGE", "false") == "true" ? "user" : "none" + test_file = joinpath(package_root, "test", "runtests.jl") + run(`$(Base.julia_cmd()) --startup-file=no --check-bounds=yes --depwarn=$depwarn --code-coverage=$coverage --project=$environment $test_file`) +end diff --git a/src/DiffPrecTest.jl b/src/DiffPrecTest.jl index 91effe6..b7c867c 100644 --- a/src/DiffPrecTest.jl +++ b/src/DiffPrecTest.jl @@ -8,8 +8,8 @@ independent samples `X` (`nx × p`) and `Y` (`ny × p`): [`SymmetricNormal`](@ref), [`AsymmetricNormal`](@ref), [`ReducedNormal`](@ref), [`SeparateNormal`](@ref) and the oracle variants); * simultaneous confidence intervals for all entries by the bootstrap - ([`bootstrap`](@ref), [`simulCI`](@ref), [`simulCIstudentized`](@ref)); -* support recovery ([`supportEstimate`](@ref)). + ([`bootstrap`](@ref), [`simul_ci`](@ref), [`simul_ci_studentized`](@ref)); +* support recovery ([`support_estimate`](@ref)). """ module DiffPrecTest @@ -19,6 +19,18 @@ using ProximalSolvers: ProximalSolvers, CDSolver, cdmodel, solve!, solve_path using Random using SpecialFunctions: erfcinv +# Canonical and deprecated keywords share one implementation; reject an +# ambiguous call instead of silently choosing one spelling. +function _deprecated_keyword(value, legacy, canonical::Symbol, old::Symbol, default) + if legacy !== nothing + value === nothing || throw(ArgumentError("pass $canonical or $old, not both")) + Base.depwarn("keyword $old is deprecated; use $canonical.", old) + return legacy + end + value === nothing ? default : value +end + + export DiffPrecResultBoot, DiffPrecResultNormal, @@ -37,22 +49,22 @@ export estimate, compute, variance, - computeSimulationResult, + compute_simulation_result, # first and second stage functions - diffEstimation, invAsymHessianEstimation, invQSymHessian, invHessianReduced, + diff_estimation, inv_asym_hessian_estimation, inv_q_sym_hessian, inv_hessian_reduced, # support estimation ANTSupport, BootStdSupport, BootMaxSupport, DTraceValidationSupport, - supportEstimate, + support_estimate, # bootstrap BootstrapEstimates, bootstrap, - simulCI, simulCIstudentized + simul_ci, simul_ci_studentized @@ -101,7 +113,7 @@ struct SymmetricOracleNormal <: DiffPrecMethod end AsymmetricNormal() Like [`SymmetricNormal`](@ref), but the second stage estimates a row of the -inverse of the asymmetric Hessian `Σy ⊗ Σx` ([`invAsymHessianEstimation`](@ref)) +inverse of the asymmetric Hessian `Σy ⊗ Σx` ([`inv_asym_hessian_estimation`](@ref)) and the refit uses [`AsymmetricOracleNormal`](@ref). """ struct AsymmetricNormal <: DiffPrecMethod end @@ -112,9 +124,9 @@ struct AsymmetricNormal <: DiffPrecMethod end The three stage procedure for one entry of `Δ`: 1. estimate `Δ` by the weighted ℓ1 penalized direct difference estimator - ([`diffEstimation`](@ref)); + ([`diff_estimation`](@ref)); 2. estimate the row of the inverse of the Hessian `(Σx ⊗ Σy + Σy ⊗ Σx) / 2` - that belongs to the entry ([`invQSymHessian`](@ref)); + that belongs to the entry ([`inv_q_sym_hessian`](@ref)); 3. refit on the union of the two supports ([`SymmetricOracleNormal`](@ref)). See [`estimate`](@ref). @@ -151,8 +163,8 @@ struct DecorrelatedScoreOracleNormal <: DiffPrecMethod end ReducedNormal() The procedure of [`SymmetricNormal`](@ref) in the half-vectorized -parametrization `svec(Δ)` of a symmetric `Δ` ([`reducedDiffEstimation`](@ref), -[`invHessianReduced`](@ref), [`ReducedOracleNormal`](@ref)). See +parametrization `svec(Δ)` of a symmetric `Δ` ([`reduced_diff_estimation`](@ref), +[`inv_hessian_reduced`](@ref), [`ReducedOracleNormal`](@ref)). See [`estimate`](@ref). """ struct ReducedNormal <: DiffPrecMethod end @@ -195,14 +207,14 @@ end SimulationResult Summary of repeated estimates of one entry of `Δ` returned by -[`computeSimulationResult`](@ref): `bias`, `coverage` and average length -`lenCoverage` of the confidence intervals, and the mean `p_mean` and variance +[`compute_simulation_result`](@ref): `bias`, `coverage` and average length +`len_coverage` of the confidence intervals, and the mean `p_mean` and variance `p_var` of the estimates. """ struct SimulationResult bias::Float64 coverage::Float64 - lenCoverage::Float64 + len_coverage::Float64 p_mean::Float64 p_var::Float64 end @@ -216,11 +228,11 @@ include("skron.jl") include("ustats.jl") include("support.jl") include("variance.jl") -include("diffEstimation.jl") -include("invHessianEstimation.jl") +include("diff_estimation.jl") +include("inv_hessian_estimation.jl") include("util.jl") include("bootstrap.jl") -include("reducedDiffEstim.jl") +include("reduced_diff_estim.jl") ############################## @@ -234,18 +246,18 @@ include("reducedDiffEstim.jl") Estimate the entry `Δ[row, col]` of `Δ = Ωx - Ωy` from the samples `X` and `Y` with `method` one of [`SymmetricNormal`](@ref) or [`AsymmetricNormal`](@ref). -Returns `(res, Δ, suppΔ, ω, suppω, indS)`: the [`DiffPrecResultNormal`](@ref) +Returns `(res, Δ, suppΔ, ω, suppω, ind_s)`: the [`DiffPrecResultNormal`](@ref) `res`, the first stage estimate `Δ` and its support `suppΔ`, the second stage estimate `ω` and its support `suppω` (both supports as `p × p` `BitMatrix`), -and the linear indices `indS` of the refitted entries. +and the linear indices `ind_s` of the refitted entries. The sample covariances `Sx = X'X/nx`, `Sy = Y'Y/ny` and the first and second stage estimates are computed unless given; when estimating several entries, pass the first stage results of the first call to the later ones. - estimate(method, X, Y, indElem; Sx, Sy, H, Δ, suppΔ, ω, suppω) + estimate(method, X, Y, ind_elem; Sx, Sy, H, Δ, suppΔ, ω, suppω) -The same with `method = ReducedNormal()` for the entry with index `indElem` in +The same with `method = ReducedNormal()` for the entry with index `ind_elem` in the lower triangle (see `ProximalBase.tril_index`). Here `Δ` and `ω` are vectors of length `p(p+1)/2` with supports `suppΔ` and `suppω` given as `BitVector`s. @@ -254,19 +266,19 @@ are vectors of length `p(p+1)/2` with supports `suppΔ` and `suppω` given as Estimator of Xia, Cai and Cai (2015). Returns a [`DiffPrecResultNormal`](@ref). - estimate(method, Sx, Sy, X, Y, row, col, indS) + estimate(method, Sx, Sy, X, Y, row, col, ind_s) -Refit on the support `indS`, given as linear indices into a `p × p` matrix, +Refit on the support `ind_s`, given as linear indices into a `p × p` matrix, with `method` one of [`SymmetricOracleNormal`](@ref), [`AsymmetricOracleNormal`](@ref) (returning a [`DiffPrecResultNormal`](@ref)), [`SymmetricOracleBoot`](@ref) or [`AsymmetricOracleBoot`](@ref) (returning a -[`DiffPrecResultBoot`](@ref); keyword arguments `bootSamples=1000` and +[`DiffPrecResultBoot`](@ref); keyword arguments `boot_samples=1000` and `rng=Random.default_rng()`). The entries `(row, col)` and `(col, row)` are moved -to the front of `indS`, which is modified. +to the front of `ind_s`, which is modified. - estimate(ReducedOracleNormal(), Sx, Sy, X, Y, indElem, indS) + estimate(ReducedOracleNormal(), Sx, Sy, X, Y, ind_elem, ind_s) -Refit on the support `indS`, given as indices into `svec(Δ)`. +Refit on the support `ind_s`, given as indices into `svec(Δ)`. estimate(DecorrelatedScoreOracleNormal(), Sx, Sy, X, Y, row, col, ω, indω, Δn, indΔn) @@ -304,28 +316,28 @@ function estimate(::SymmetricNormal, X, Y, row, col; end # first stage - λ = _penaltyLevel(p) + λ = _penalty_level(p) if Δ === nothing && suppΔ === nothing - Δ = diffEstimation(Sx, Sy, X, Y, λ) - suppΔ = getSupport(Δ) + Δ = diff_estimation(Sx, Sy, X, Y, λ) + suppΔ = get_support(Δ) end if suppΔ === nothing - suppΔ = getSupport(Δ) + suppΔ = get_support(Δ) end # second stage if ω === nothing && suppω === nothing j = tril_index(p, row, col) - ω = invQSymHessian(Sx, Sy, j, X, Y, λ) - suppω = getSupport(ω, p) + ω = inv_q_sym_hessian(Sx, Sy, j, X, Y, λ) + suppω = get_support(ω, p) end if suppω === nothing - suppω = getSupport(ω, p) + suppω = get_support(ω, p) end # refit stage - indS = getLinearSupport(row, col, suppΔ, suppω) - estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS), Δ, suppΔ, ω, suppω, indS + ind_s = get_linear_support(row, col, suppΔ, suppω) + estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, ind_s), Δ, suppΔ, ω, suppω, ind_s end @@ -357,27 +369,27 @@ function estimate(::AsymmetricNormal, X, Y, row, col; end # first stage - λ = _penaltyLevel(p) + λ = _penalty_level(p) if Δ === nothing && suppΔ === nothing - Δ = diffEstimation(Sx, Sy, X, Y, λ) - suppΔ = getSupport(Δ) + Δ = diff_estimation(Sx, Sy, X, Y, λ) + suppΔ = get_support(Δ) end if suppΔ === nothing - suppΔ = getSupport(Δ) + suppΔ = get_support(Δ) end # second stage if ω === nothing && suppω === nothing - ω = invAsymHessianEstimation(Sx, Sy, row, col, X, Y, λ) - suppω = getSupport(ω, p) + ω = inv_asym_hessian_estimation(Sx, Sy, row, col, X, Y, λ) + suppω = get_support(ω, p) end if suppω === nothing - suppω = getSupport(ω, p) + suppω = get_support(ω, p) end # refit stage - indS = getLinearSupport(row, col, suppΔ, suppω) - estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS), Δ, suppΔ, ω, suppω, indS + ind_s = get_linear_support(row, col, suppΔ, suppω) + estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, row, col, ind_s), Δ, suppΔ, ω, suppω, ind_s end @@ -457,7 +469,7 @@ end #################################### -# indS --- coordinates for the nonzero element of the true Δ (LinearIndices) +# ind_s --- coordinates for the nonzero element of the true Δ (LinearIndices) function estimate(::AsymmetricOracleBoot, Sx::Symmetric, Sy::Symmetric, @@ -465,28 +477,29 @@ function estimate(::AsymmetricOracleBoot, Y, row, col, - indS::Vector{Int64}; - bootSamples::Int64=1000, + ind_s::Vector{Int64}; + boot_samples=nothing, bootSamples=nothing, rng::AbstractRNG=Random.default_rng() ) + boot_samples = _deprecated_keyword(boot_samples, bootSamples, :boot_samples, :bootSamples, 1000) nx, p = size(X) ny = size(Y, 1) - makeFirst!(indS, p, row, col) + make_first!(ind_s, p, row, col) - A = kron_sub(Sy, Sx, indS, indS) - Δ = A \ (Sy[indS] - Sx[indS]) + A = kron_sub(Sy, Sx, ind_s, ind_s) + Δ = A \ (Sy[ind_s] - Sx[ind_s]) Δab = row == col ? Δ[1] : ( Δ[1] + Δ[2] ) / 2. - boot_est = zeros(Float64, bootSamples) - for b=1:bootSamples - bX_ind = rand(rng, 1:nx, nx) - bY_ind = rand(rng, 1:ny, ny) - bSx = cov(X[bX_ind, :]) - bSy = cov(Y[bY_ind, :]) - kron_sub!(A, bSy, bSx, indS, indS) - Δb = A \ (bSy[indS] - bSx[indS]) + boot_est = zeros(Float64, boot_samples) + for b=1:boot_samples + b_x_ind = rand(rng, 1:nx, nx) + b_y_ind = rand(rng, 1:ny, ny) + b_sx = cov(X[b_x_ind, :]) + b_sy = cov(Y[b_y_ind, :]) + kron_sub!(A, b_sy, b_sx, ind_s, ind_s) + Δb = A \ (b_sy[ind_s] - b_sx[ind_s]) boot_est[b] = row == col ? Δb[1] : ( Δb[1] + Δb[2] ) / 2. end @@ -494,7 +507,7 @@ function estimate(::AsymmetricOracleBoot, end -# indS --- coordinates for the nonzero element of the true Δ (LinearIndices) +# ind_s --- coordinates for the nonzero element of the true Δ (LinearIndices) function estimate(::SymmetricOracleBoot, Sx::Symmetric, Sy::Symmetric, @@ -502,27 +515,28 @@ function estimate(::SymmetricOracleBoot, Y, row, col, - indS::Vector{Int64}; - bootSamples::Int64=1000, + ind_s::Vector{Int64}; + boot_samples=nothing, bootSamples=nothing, rng::AbstractRNG=Random.default_rng() ) + boot_samples = _deprecated_keyword(boot_samples, bootSamples, :boot_samples, :bootSamples, 1000) nx, p = size(X) ny = size(Y, 1) - makeFirst!(indS, p, row, col) + make_first!(ind_s, p, row, col) - C = skron_sub(Sx, Sy, indS, indS) - Δab = (C \ (Sy[indS] - Sx[indS]))[1] + C = skron_sub(Sx, Sy, ind_s, ind_s) + Δab = (C \ (Sy[ind_s] - Sx[ind_s]))[1] - boot_est = zeros(Float64, bootSamples) - for b=1:bootSamples - bX_ind = rand(rng, 1:nx, nx) - bY_ind = rand(rng, 1:ny, ny) - bSx = cov(X[bX_ind, :]) - bSy = cov(Y[bY_ind, :]) - skron_sub!(C, bSx, bSy, indS, indS) - boot_est[b] = (C \ (bSy[indS] - bSx[indS]))[1] + boot_est = zeros(Float64, boot_samples) + for b=1:boot_samples + b_x_ind = rand(rng, 1:nx, nx) + b_y_ind = rand(rng, 1:ny, ny) + b_sx = cov(X[b_x_ind, :]) + b_sy = cov(Y[b_y_ind, :]) + skron_sub!(C, b_sx, b_sy, ind_s, ind_s) + boot_est[b] = (C \ (b_sy[ind_s] - b_sx[ind_s]))[1] end DiffPrecResultBoot(Δab, boot_est) @@ -540,22 +554,22 @@ function estimate( Y, row, col, - indS::Vector{Int64} + ind_s::Vector{Int64} ) nx, p = size(X) ny = size(Y, 1) - makeFirst!(indS, p, row, col) + make_first!(ind_s, p, row, col) - A = kron_sub(Sy, Sx, indS, indS) - Δ = A \ (Sy[indS] - Sx[indS]) + A = kron_sub(Sy, Sx, ind_s, ind_s) + Δ = A \ (Sy[ind_s] - Sx[ind_s]) - tmp = zeros(length(indS)) + tmp = zeros(length(ind_s)) tmp[1] = 1. ω = A \ tmp - v = variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, indS, Δ, indS) + v = variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind_s, Δ, ind_s) DiffPrecResultNormal(Δ[1], sqrt(v)) end @@ -572,20 +586,20 @@ function estimate( Y, row, col, - indS::Vector{Int64} + ind_s::Vector{Int64} ) nx, p = size(X) ny = size(Y, 1) - makeFirst!(indS, p, row, col) - C = skron_sub(Sx, Sy, indS, indS) + make_first!(ind_s, p, row, col) + C = skron_sub(Sx, Sy, ind_s, ind_s) - tmp = zeros(length(indS)) + tmp = zeros(length(ind_s)) tmp[1] = 1. ω = C \ tmp - Δ = C \ (Sy[indS] - Sx[indS]) - v = variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, indS, Δ, indS) + Δ = C \ (Sy[ind_s] - Sx[ind_s]) + v = variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind_s, Δ, ind_s) DiffPrecResultNormal(Δ[1], sqrt(v)) end @@ -613,7 +627,7 @@ function estimate( v = 0. for ci=1:length(indΔn) for ri=1:length(indω) - @inbounds v += _getElemSKron(Sx, Sy, indω[ri], indΔn[ci]) * ω[ri] * Δn[ci] + @inbounds v += _get_elem_skron(Sx, Sy, indω[ri], indΔn[ci]) * ω[ri] * Δn[ci] end end for ri=1:length(indω) @@ -623,13 +637,13 @@ function estimate( if row == col ind = (col-1)*p + row for ri=1:length(indω) - @inbounds a += ω[ri] * _getElemSKron(Sx, Sy, indω[ri], ind) + @inbounds a += ω[ri] * _get_elem_skron(Sx, Sy, indω[ri], ind) end else ind1 = (col-1)*p + row ind2 = (row-1)*p + col for ri=1:length(indω) - @inbounds a += ω[ri] * (_getElemSKron(Sx, Sy, indω[ri], ind1) + _getElemSKron(Sx, Sy, indω[ri], ind2)) + @inbounds a += ω[ri] * (_get_elem_skron(Sx, Sy, indω[ri], ind1) + _get_elem_skron(Sx, Sy, indω[ri], ind2)) end end Δab = -v / a @@ -647,6 +661,8 @@ function estimate( DiffPrecResultNormal(Δab, sqrt(v)) end +include("deprecated.jl") + include("precompile.jl") end diff --git a/src/bootstrap.jl b/src/bootstrap.jl index 2f06342..5c97817 100644 --- a/src/bootstrap.jl +++ b/src/bootstrap.jl @@ -12,37 +12,39 @@ end """ - bootstrap(X, Y; bootSamples=300, estimSupport=nothing, rng=Random.default_rng()) + bootstrap(X, Y; boot_samples=300, estim_support=nothing, rng=Random.default_rng()) Refit every entry of `svec(Δ)` on its estimated support and compute -`bootSamples` bootstrap replicates of the refitted estimates, resampling the -rows of `X` and `Y`. Returns `(straps, eS)`: the [`BootstrapEstimates`](@ref) +`boot_samples` bootstrap replicates of the refitted estimates, resampling the +rows of `X` and `Y`. Returns `(straps, e_s)`: the [`BootstrapEstimates`](@ref) and the supports, as vectors of indices into `svec(Δ)`, one per entry. The supports are computed by the first and second stage estimators of -[`ReducedNormal`](@ref) unless given in `estimSupport`; the entry itself is +[`ReducedNormal`](@ref) unless given in `estim_support`; the entry itself is moved to the front of each support, which is modified. -Use [`simulCI`](@ref) or [`simulCIstudentized`](@ref) for simultaneous +Use [`simul_ci`](@ref) or [`simul_ci_studentized`](@ref) for simultaneous confidence intervals. """ -function bootstrap(X, Y; bootSamples::Int64=300, estimSupport::Union{Vector{Vector{Int64}},Nothing}=nothing, +function bootstrap(X, Y; boot_samples=nothing, bootSamples=nothing, estim_support=nothing, estimSupport=nothing, rng::AbstractRNG=Random.default_rng()) + boot_samples = _deprecated_keyword(boot_samples, bootSamples, :boot_samples, :bootSamples, 300) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) nx, p = size(X) ny = size(Y, 1) Sx = Symmetric( cov(X) ) Sy = Symmetric( cov(Y) ) - eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport + e_s = estim_support === nothing ? __init_support(Sx, Sy, X, Y) : estim_support rp = div((p+1)*p, 2) - # indS = Vector{Vector{Int64}}(undef, rp) + # ind_s = Vector{Vector{Int64}}(undef, rp) Δhat = Vector{Float64}(undef, rp) - Δb = Matrix{Float64}(undef, rp, bootSamples) + Δb = Matrix{Float64}(undef, rp, boot_samples) - bX = similar(X) - bY = similar(Y) + b_x = similar(X) + b_y = similar(Y) x_ind = Vector{Int64}(undef, nx) y_ind = Vector{Int64}(undef, ny) @@ -51,40 +53,40 @@ function bootstrap(X, Y; bootSamples::Int64=300, estimSupport::Union{Vector{Vect for col=1:p for row=col:p it = it + 1 - indElem = tril_index(p, row, col) - # tmp = estimate(ReducedOracleNormal(), Sx, Sy, X, Y, indElem, eS[it]) + ind_elem = tril_index(p, row, col) + # tmp = estimate(ReducedOracleNormal(), Sx, Sy, X, Y, ind_elem, e_s[it]) - indS = eS[it] - makeFirst!(indS, indElem) - C = skron(Sx, Sy, indS) - b = svec(Sy, indS) - svec(Sx, indS) + ind_s = e_s[it] + make_first!(ind_s, ind_elem) + C = skron(Sx, Sy, ind_s) + b = svec(Sy, ind_s) - svec(Sx, ind_s) Δ = C \ b Δhat[it] = Δ[1] end end - for _b=1:bootSamples + for _b=1:boot_samples rand!(rng, x_ind, 1:nx) rand!(rng, y_ind, 1:ny) - _fill_boot!(bX, X, x_ind) - _fill_boot!(bY, Y, y_ind) + _fill_boot!(b_x, X, x_ind) + _fill_boot!(b_y, Y, y_ind) - bSx = Symmetric( cov(bX) ) - bSy = Symmetric( cov(bY) ) + b_sx = Symmetric( cov(b_x) ) + b_sy = Symmetric( cov(b_y) ) it = 0 for col=1:p for row=col:p it = it + 1 - indElem = tril_index(p, row, col) - # tmp = estimate(ReducedOracleNormal(), bSx, bSy, bX, bY, indexElementLinear, eS[it]) + ind_elem = tril_index(p, row, col) + # tmp = estimate(ReducedOracleNormal(), b_sx, b_sy, b_x, b_y, index_element_linear, e_s[it]) - indS = eS[it] - makeFirst!(indS, indElem) - C = skron(bSx, bSy, indS) - b = svec(bSy, indS) - svec(bSx, indS) + ind_s = e_s[it] + make_first!(ind_s, ind_elem) + C = skron(b_sx, b_sy, ind_s) + b = svec(b_sy, ind_s) - svec(b_sx, ind_s) Δ = C \ b Δb[it, _b] = Δ[1] @@ -92,7 +94,7 @@ function bootstrap(X, Y; bootSamples::Int64=300, estimSupport::Union{Vector{Vect end end - BootstrapEstimates(Δhat, Δb), eS + BootstrapEstimates(Δhat, Δb), e_s end @@ -100,23 +102,23 @@ end """ - simulCI(straps::BootstrapEstimates, α=0.95) + simul_ci(straps::BootstrapEstimates, α=0.95) Simultaneous confidence intervals with coverage `α` for the entries of `svec(Δ)`, `Δhat .± x` with `x` the `α` quantile of the bootstrap distribution of `‖Δhat - Δb‖∞`. Returns a matrix with the lower ends of the intervals in the first column and the upper ends in the second. """ -function simulCI(straps::BootstrapEstimates, α::Float64=0.95) - m, bootSamples = size(straps.Δb) +function simul_ci(straps::BootstrapEstimates, α::Float64=0.95) + m, boot_samples = size(straps.Δb) - infNormDist = Vector{Float64}(undef, bootSamples) + inf_norm_dist = Vector{Float64}(undef, boot_samples) CI = Matrix{Float64}(undef, m, 2) - for b=1:bootSamples - infNormDist[b] = norm_diff(straps.Δhat, view(straps.Δb, :, b), Inf) + for b=1:boot_samples + inf_norm_dist[b] = norm_diff(straps.Δhat, view(straps.Δb, :, b), Inf) end - x = quantile!(infNormDist, α) + x = quantile!(inf_norm_dist, α) CI[:, 1] .= straps.Δhat .- x CI[:, 2] .= straps.Δhat .+ x @@ -124,26 +126,26 @@ function simulCI(straps::BootstrapEstimates, α::Float64=0.95) end """ - simulCIstudentized(straps::BootstrapEstimates, α=0.95) + simul_ci_studentized(straps::BootstrapEstimates, α=0.95) -Like [`simulCI`](@ref), but every entry is standardized by its bootstrap +Like [`simul_ci`](@ref), but every entry is standardized by its bootstrap standard error `w`: the intervals are `Δhat .± x .* w` with `x` the `α` quantile of the bootstrap distribution of `‖(Δhat - Δb) ./ w‖∞`. """ -function simulCIstudentized(straps::BootstrapEstimates, α::Float64=0.95) - m, bootSamples = size(straps.Δb) +function simul_ci_studentized(straps::BootstrapEstimates, α::Float64=0.95) + m, boot_samples = size(straps.Δb) - infNormDist = Vector{Float64}(undef, bootSamples) + inf_norm_dist = Vector{Float64}(undef, boot_samples) w = reshape(std(straps.Δb; dims = 2, corrected = false), :) CI = Matrix{Float64}(undef, m, 2) tmp = Vector{Float64}(undef, m) - for b=1:bootSamples + for b=1:boot_samples tmp .= (straps.Δhat .- straps.Δb[:, b]) ./ w - infNormDist[b] = norm(tmp, Inf) + inf_norm_dist[b] = norm(tmp, Inf) end - x = quantile!(infNormDist, α) + x = quantile!(inf_norm_dist, α) @. CI[:, 1] = straps.Δhat - x * w @. CI[:, 2] = straps.Δhat + x * w @@ -151,11 +153,11 @@ function simulCIstudentized(straps::BootstrapEstimates, α::Float64=0.95) end -function _fill_boot!(bX, X, b_ind) +function _fill_boot!(b_x, X, b_ind) n, p = size(X) for col=1:p for row=1:n - bX[row, col] = X[b_ind[row], col] + b_x[row, col] = X[b_ind[row], col] end end end diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 0000000..aa13b9b --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,28 @@ +# Compatibility spellings for the naming migration. Binding aliases retain +# generic-function identity so downstream method extensions keep dispatching. +Base.@deprecate_binding computeSimulationResult compute_simulation_result +Base.@deprecate_binding diffEstimation diff_estimation +Base.@deprecate_binding invAsymHessianEstimation inv_asym_hessian_estimation +Base.@deprecate_binding invQSymHessian inv_q_sym_hessian +Base.@deprecate_binding invHessianReduced inv_hessian_reduced +Base.@deprecate_binding supportEstimate support_estimate +Base.@deprecate_binding simulCI simul_ci +Base.@deprecate_binding simulCIstudentized simul_ci_studentized +Base.@deprecate_binding getSupport get_support false +Base.@deprecate_binding getReducedSupport get_reduced_support false +Base.@deprecate_binding getLinearSupport get_linear_support false +Base.@deprecate_binding makeFirst! make_first! false +Base.@deprecate_binding createVarianceStatistic create_variance_statistic false +Base.@deprecate_binding createCovarianceStatistic create_covariance_statistic false +Base.@deprecate_binding createSecondMomentStatistic create_second_moment_statistic false +Base.@deprecate_binding reducedDiffEstimation reduced_diff_estimation false + +@inline function Base.getproperty(result::SimulationResult, name::Symbol) + if name === :lenCoverage + Base.depwarn("property lenCoverage is deprecated; use len_coverage.", :lenCoverage) + return getfield(result, :len_coverage) + end + getfield(result, name) +end +Base.propertynames(::SimulationResult, private::Bool=false) = + (fieldnames(SimulationResult)..., :lenCoverage) diff --git a/src/diffEstimation.jl b/src/diff_estimation.jl similarity index 85% rename from src/diffEstimation.jl rename to src/diff_estimation.jl index 79a8f01..6628393 100644 --- a/src/diffEstimation.jl +++ b/src/diff_estimation.jl @@ -2,7 +2,7 @@ # A = Sx⋅Δ⋅Sy -function _computeVarElem( +function _compute_var_elem( q, r, Sx, Sy, A, SxΔ, SyΔ, @@ -100,7 +100,7 @@ end # the weights of the penalty on Δ as a p × p matrix, symmetric -function _computeVar!(ω, Sx, Sy, X, Y, Δ::SymmetricSparseIterate) +function _compute_var!(ω, Sx, Sy, X, Y, Δ::SymmetricSparseIterate) nx, p = size(X) ny = size(Y, 1) @@ -113,7 +113,7 @@ function _computeVar!(ω, Sx, Sy, X, Y, Δ::SymmetricSparseIterate) for col=1:p for row=col:p - ω[row, col] = ω[col, row] = sqrt(_computeVarElem(q, r, Sx, Sy, A, SxΔ, SyΔ, X, Y, row, col)) + ω[row, col] = ω[col, row] = sqrt(_compute_var_elem(q, r, Sx, Sy, A, SxΔ, SyΔ, X, Y, row, col)) end end ω @@ -122,7 +122,7 @@ end # The estimators of both stages solve twice: with the weights ω of the ℓ1 # penalty computed by weights!(ω, x) at the starting point x, then, warm # started, with the weights recomputed at the first solution. -function _twoStage!(x, f, λ, ω, weights!, method::CDSolver) +function _two_stage!(x, f, λ, ω, weights!, method::CDSolver) m = cdmodel(f, x) for stage=1:2 weights!(ω, x) @@ -134,19 +134,19 @@ end """ - diffEstimation(Sx, Sy, X, Y, λ, method=CDSolver()) + diff_estimation(Sx, Sy, X, Y, λ, method=CDSolver()) First stage of [`SymmetricNormal`](@ref) and [`AsymmetricNormal`](@ref): estimate `Δ = Ωx - Ωy` by minimizing the direct difference loss -`tr(Sx⋅Δ⋅Sy⋅Δ)/2 - tr(Δ⋅(Sy - Sx))` (`CovarianceSelection.directDifferenceLoss`) plus a +`tr(Sx⋅Δ⋅Sy⋅Δ)/2 - tr(Δ⋅(Sy - Sx))` (`CovarianceSelection.direct_difference_loss`) plus a weighted ℓ1 penalty with level `λ`, by ProximalSolvers' coordinate descent with the options `method`. The weights are standard errors of the gradient estimated from `X` and `Y`, computed at zero and again at the first solution. Returns a `SymmetricSparseIterate`. """ -function diffEstimation(Sx::Symmetric, Sy::Symmetric, X, Y, λ, method::CDSolver=CDSolver()) +function diff_estimation(Sx::Symmetric, Sy::Symmetric, X, Y, λ, method::CDSolver=CDSolver()) p = size(Sx, 1) x = SymmetricSparseIterate(p) ω = zeros(p, p) - _twoStage!(x, directDifferenceLoss(Sx, Sy), λ, ω, (ω, x) -> _computeVar!(ω, Sx, Sy, X, Y, x), method) + _two_stage!(x, direct_difference_loss(Sx, Sy), λ, ω, (ω, x) -> _compute_var!(ω, Sx, Sy, X, Y, x), method) end diff --git a/src/invHessianEstimation.jl b/src/inv_hessian_estimation.jl similarity index 86% rename from src/invHessianEstimation.jl rename to src/inv_hessian_estimation.jl index 4c4ba44..628aa4d 100644 --- a/src/invHessianEstimation.jl +++ b/src/inv_hessian_estimation.jl @@ -1,6 +1,6 @@ -include("invHessianEstimationASym.jl") -include("invHessianEstimationSymReduced.jl") +include("inv_hessian_estimation_asym.jl") +include("inv_hessian_estimation_sym_reduced.jl") ######################## diff --git a/src/invHessianEstimationASym.jl b/src/inv_hessian_estimation_asym.jl similarity index 77% rename from src/invHessianEstimationASym.jl rename to src/inv_hessian_estimation_asym.jl index 7ee632c..d6c5b7a 100644 --- a/src/invHessianEstimationASym.jl +++ b/src/inv_hessian_estimation_asym.jl @@ -7,7 +7,7 @@ # The loss tr(Σx⋅Θ⋅Σy⋅Θ')/2 - Θ[a, b] of a p × p matrix Θ, held as vec(Θ) in a # SparseIterate, as ProximalBase's KronQuadratic. Its minimizer is the row # (b - 1) p + a of the inverse of Σy ⊗ Σx. -function _inverseKroneckerLoss(Σx::Symmetric, Σy::Symmetric, a::Int, b::Int) +function _inverse_kronecker_loss(Σx::Symmetric, Σy::Symmetric, a::Int, b::Int) E = zeros(size(Σx)) E[a, b] = -1. KronQuadratic(Σx, Σy, E) @@ -16,12 +16,12 @@ end ############################################ # -# asymmetric version of invHessian +# asymmetric version of inv_hessian # ############################################ # A = Sx⋅Θ⋅Sy -function _computeAsymVarElem( +function _compute_asym_var_elem( q, r, A, a, b, SxΘ, ΘSy, @@ -64,7 +64,7 @@ function _computeAsymVarElem( end # the weights of the penalty on vec(Θ); the entry (a, b) is not penalized -function _computeAsymVar!(ω, Sx, Sy, a, b, X, Y, Θ) +function _compute_asym_var!(ω, Sx, Sy, a, b, X, Y, Θ) nx, p = size(X) ny = size(Y, 1) @@ -77,8 +77,8 @@ function _computeAsymVar!(ω, Sx, Sy, a, b, X, Y, Θ) for col=1:p for row=1:p - linInd = (col-1)*p + row - ω[linInd] = sqrt(_computeAsymVarElem(q, r, A, a, b, SxΘ, ΘSy, X, Y, row, col)) + lin_ind = (col-1)*p + row + ω[lin_ind] = sqrt(_compute_asym_var_elem(q, r, A, a, b, SxΘ, ΘSy, X, Y, row, col)) end end ω[(b-1)*p + a] = 0. @@ -88,7 +88,7 @@ end """ - invAsymHessianEstimation(Sx, Sy, ri, ci, X, Y, λ, method=CDSolver()) + inv_asym_hessian_estimation(Sx, Sy, ri, ci, X, Y, λ, method=CDSolver()) Second stage of [`AsymmetricNormal`](@ref): estimate the row `(ci - 1) p + ri` of the inverse of `Σy ⊗ Σx` by minimizing the loss `tr(Σx⋅Θ⋅Σy⋅Θ')/2 - Θ[ri, ci]` @@ -99,11 +99,11 @@ and `Y`, computed at the starting point, the indicator of the entry, and again at the first solution. Returns a `SparseIterate` holding `vec(Θ)`, of length `p²`. """ -function invAsymHessianEstimation(Sx::Symmetric, Sy::Symmetric, ri, ci, X, Y, λ, method::CDSolver=CDSolver()) +function inv_asym_hessian_estimation(Sx::Symmetric, Sy::Symmetric, ri, ci, X, Y, λ, method::CDSolver=CDSolver()) p = size(X, 2) x = SparseIterate(p*p) x[(ci-1)*p + ri] = 1. ω = zeros(p*p) - _twoStage!(x, _inverseKroneckerLoss(Sx, Sy, ri, ci), λ, ω, - (ω, x) -> _computeAsymVar!(ω, Sx, Sy, ri, ci, X, Y, x), method) + _two_stage!(x, _inverse_kronecker_loss(Sx, Sy, ri, ci), λ, ω, + (ω, x) -> _compute_asym_var!(ω, Sx, Sy, ri, ci, X, Y, x), method) end diff --git a/src/invHessianEstimationSymReduced.jl b/src/inv_hessian_estimation_sym_reduced.jl similarity index 86% rename from src/invHessianEstimationSymReduced.jl rename to src/inv_hessian_estimation_sym_reduced.jl index 2c3c67e..37ba0f2 100644 --- a/src/invHessianEstimationSymReduced.jl +++ b/src/inv_hessian_estimation_sym_reduced.jl @@ -10,7 +10,7 @@ # C = -(E_rc + E_cr)/√2 off it, (r, c) the entry of ind. In svec(Θ) the # off-diagonal entries are scaled by √2, so the penalty weight ω of θ[k] is ω on # Θ[r, r], and ω/√2 on each of Θ[r, c] and Θ[c, r]. -function _inverseSymKroneckerLoss(Σx::Symmetric, Σy::Symmetric, ind::Int) +function _inverse_sym_kronecker_loss(Σx::Symmetric, Σy::Symmetric, ind::Int) p = size(Σx, 1) r, c = tril_cartesian(p, ind) C = zeros(p, p) @@ -37,7 +37,7 @@ function _svec(x::SymmetricSparseIterate{T}) where {T} end # the weights W of the penalty on the symmetric Θ from the weights ω of svec(Θ) -function _symWeights!(W, ω, p) +function _sym_weights!(W, ω, p) k = 0 for c=1:p, r=c:p k += 1 @@ -49,11 +49,11 @@ end ############################################ # -# compressed version of invHessian +# compressed version of inv_hessian # ############################################ -function _mul_symX_S(θ::SparseIterate, S::Symmetric) +function _mul_sym_x_s(θ::SparseIterate, S::Symmetric) p = size(S, 1) rp = length(θ) _c = 1. / sqrt(2.) @@ -83,7 +83,7 @@ end # the weights of the penalty on svec(Θ); the entry ind is not penalized -function _computeReducedVar!(ω, Sx, Sy, ind, X, Y, Θ) +function _compute_reduced_var!(ω, Sx, Sy, ind, X, Y, Θ) nx, p = size(X) ny = size(Y, 1) @@ -92,8 +92,8 @@ function _computeReducedVar!(ω, Sx, Sy, ind, X, Y, Θ) _c = 1. / sqrt(2.) - ΘSx = _mul_symX_S(Θ, Sx) - ΘSy = _mul_symX_S(Θ, Sy) + ΘSx = _mul_sym_x_s(Θ, Sx) + ΘSy = _mul_sym_x_s(Θ, Sy) A = Sx * ΘSy # Sx⋅Θ⋅Sy V = X * ΘSy @@ -141,7 +141,7 @@ function _computeReducedVar!(ω, Sx, Sy, ind, X, Y, Θ) end """ - invQSymHessian(Sx, Sy, ind, X, Y, λ, method=CDSolver()) + inv_q_sym_hessian(Sx, Sy, ind, X, Y, λ, method=CDSolver()) Second stage of [`SymmetricNormal`](@ref): estimate the row `ind` of the inverse of `H = Q (Σx ⊗ Σy + Σy ⊗ Σx) Q' / 2`, where `ind` indexes the lower triangle @@ -152,7 +152,7 @@ standard errors of the gradient estimated from `X` and `Y`, computed at the starting point `e_ind / H[ind, ind]` and again at the first solution. Returns a `SparseIterate` of length `p(p+1)/2`. """ -function invQSymHessian(Sx::Symmetric, Sy::Symmetric, ind, X, Y, λ, method::CDSolver=CDSolver()) +function inv_q_sym_hessian(Sx::Symmetric, Sy::Symmetric, ind, X, Y, λ, method::CDSolver=CDSolver()) p = size(X, 2) # coordinate descent on the symmetric Θ = smat(θ), which takes the steps of @@ -164,7 +164,7 @@ function invQSymHessian(Sx::Symmetric, Sy::Symmetric, ind, X, Y, λ, method::CDS ω = zeros(div((p+1)*p, 2)) W = zeros(p, p) - weights!(W, x) = _symWeights!(W, _computeReducedVar!(ω, Sx, Sy, ind, X, Y, _svec(x)), p) - _twoStage!(x, _inverseSymKroneckerLoss(Sx, Sy, ind), λ, W, weights!, method) + weights!(W, x) = _sym_weights!(W, _compute_reduced_var!(ω, Sx, Sy, ind, X, Y, _svec(x)), p) + _two_stage!(x, _inverse_sym_kronecker_loss(Sx, Sy, ind), λ, W, weights!, method) _svec(x) end diff --git a/src/precompile.jl b/src/precompile.jl index a416c35..4f65abd 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -19,7 +19,7 @@ using PrecompileTools: @setup_workload, @compile_workload estimate(AsymmetricNormal(), X, Y, 2, 1) estimate(SeparateNormal(), X, Y, 2, 1) estimate(ReducedNormal(), X, Y, tril_index(p, 2, 1)) - straps, = bootstrap(X, Y; bootSamples = 2) - simulCIstudentized(straps, 0.95) + straps, = bootstrap(X, Y; boot_samples = 2) + simul_ci_studentized(straps, 0.95) end end diff --git a/src/reducedDiffEstim.jl b/src/reduced_diff_estim.jl similarity index 70% rename from src/reducedDiffEstim.jl rename to src/reduced_diff_estim.jl index 1477239..4b43255 100644 --- a/src/reducedDiffEstim.jl +++ b/src/reduced_diff_estim.jl @@ -11,21 +11,21 @@ function estimate( Sy::Symmetric, X, Y, - indElem, - indS::Vector{Int64} + ind_elem, + ind_s::Vector{Int64} ) - makeFirst!(indS, indElem) - C = skron(Sx, Sy, indS) - b = svec(Sy, indS) - svec(Sx, indS) + make_first!(ind_s, ind_elem) + C = skron(Sx, Sy, ind_s) + b = svec(Sy, ind_s) - svec(Sx, ind_s) - tmp = zeros(length(indS)) + tmp = zeros(length(ind_s)) tmp[1] = 1. ω = C \ tmp Δ = C \ b - v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, indS, Δ, indS) + v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, ind_s, Δ, ind_s) - row, col = tril_cartesian(size(Sx, 1), indElem) + row, col = tril_cartesian(size(Sx, 1), ind_elem) if row == col return DiffPrecResultNormal(Δ[1], sqrt(v)) else @@ -39,12 +39,12 @@ end -function _computeVarStep1!( - ω, +function _compute_var_step1!( + ω, Sx::Symmetric, Sy::Symmetric, - X, - Y, + X, + Y, Δ::SparseIterate) nx = size(X, 1) ny = size(Y, 1) @@ -53,7 +53,7 @@ function _computeVarStep1!( r = zeros(ny) # the entries (row, col) of the coordinates, computed once for the inner loops - E = _trilEntries(size(Sx, 1)) + E = _tril_entries(size(Sx, 1)) EΔ = [E[Δ.nzval2ind[inz]] for inz=1:nnz(Δ)] for ri=1:length(ω) @@ -85,13 +85,13 @@ function _computeVarStep1!( σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2 σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2 - + ω[ri] = sqrt( σ1/nx + σ2/ny ) end end """ - reducedDiffEstimation(H, b, Sx, Sy, X, Y, λ, method=CDSolver()) + reduced_diff_estimation(H, b, Sx, Sy, X, Y, λ, method=CDSolver()) First stage of [`ReducedNormal`](@ref): minimize `x'Hx/2 + b'x` plus a weighted ℓ1 penalty with level `λ`, where `H = Q (Sx ⊗ Sy + Sy ⊗ Sx) Q' / 2` and @@ -100,7 +100,7 @@ coordinate descent with the options `method`. The weights are standard errors of the gradient estimated from `X` and `Y`, computed at zero and again at the first solution. Returns a `SparseIterate` of length `p(p+1)/2`. """ -function reducedDiffEstimation( +function reduced_diff_estimation( H::Symmetric, b::Vector, Sx::Symmetric, @@ -108,7 +108,7 @@ function reducedDiffEstimation( X, Y, λ, method::CDSolver=CDSolver()) x = SparseIterate(size(H, 1)) ω = zeros(size(H, 1)) - _twoStage!(x, QuadraticFunction(H, b), λ, ω, (ω, x) -> _computeVarStep1!(ω, Sx, Sy, X, Y, x), method) + _two_stage!(x, QuadraticFunction(H, b), λ, ω, (ω, x) -> _compute_var_step1!(ω, Sx, Sy, X, Y, x), method) end @@ -116,12 +116,12 @@ end -function _computeVarStep2!(ω, - indElem, +function _compute_var_step2!(ω, + ind_elem, Sx::Symmetric, Sy::Symmetric, - X, - Y, + X, + Y, Δ::SparseIterate) nx = size(X, 1) @@ -131,7 +131,7 @@ function _computeVarStep2!(ω, r = zeros(ny) # the entries (row, col) of the coordinates, computed once for the inner loops - E = _trilEntries(size(Sx, 1)) + E = _tril_entries(size(Sx, 1)) EΔ = [E[Δ.nzval2ind[inz]] for inz=1:nnz(Δ)] for ri=1:length(ω) @@ -139,7 +139,7 @@ function _computeVarStep2!(ω, fill!(r, 0.) e = E[ri] - δ = indElem == ri ? - 1. : 0. + δ = ind_elem == ri ? - 1. : 0. @inbounds for k=1:nx vx = view(X, k, :) @@ -165,24 +165,24 @@ function _computeVarStep2!(ω, σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2 σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2 - + ω[ri] = sqrt( σ1/nx + σ2/ny ) end end """ - invHessianReduced(H, indRow, Sx, Sy, X, Y, λ, method=CDSolver()) + inv_hessian_reduced(H, ind_row, Sx, Sy, X, Y, λ, method=CDSolver()) -Second stage of [`ReducedNormal`](@ref): estimate the row `indRow` of the -inverse of `H = Q (Sx ⊗ Sy + Sy ⊗ Sx) Q' / 2` by minimizing `x'Hx/2 - x[indRow]` -plus a weighted ℓ1 penalty with level `λ` that leaves `x[indRow]` unpenalized, +Second stage of [`ReducedNormal`](@ref): estimate the row `ind_row` of the +inverse of `H = Q (Sx ⊗ Sy + Sy ⊗ Sx) Q' / 2` by minimizing `x'Hx/2 - x[ind_row]` +plus a weighted ℓ1 penalty with level `λ` that leaves `x[ind_row]` unpenalized, by ProximalSolvers' coordinate descent with the options `method`. The weights are standard errors of the gradient estimated from `X` and `Y`. Returns a `SparseIterate` of length `p(p+1)/2`. """ -function invHessianReduced(H::Symmetric, - indRow::Int, +function inv_hessian_reduced(H::Symmetric, + ind_row::Int, Sx::Symmetric, Sy::Symmetric, X, @@ -191,19 +191,19 @@ function invHessianReduced(H::Symmetric, method::CDSolver=CDSolver()) b = zeros(size(H, 1)) - b[indRow] = -1. + b[ind_row] = -1. x = SparseIterate(size(H, 1)) - x[indRow] = 1. / H[indRow, indRow] + x[ind_row] = 1. / H[ind_row, ind_row] ω = zeros(size(H, 1)) - weights!(ω, x) = (_computeVarStep2!(ω, indRow, Sx, Sy, X, Y, x); ω[indRow] = 0.; ω) - _twoStage!(x, QuadraticFunction(H, b), λ, ω, weights!, method) + weights!(ω, x) = (_compute_var_step2!(ω, ind_row, Sx, Sy, X, Y, x); ω[ind_row] = 0.; ω) + _two_stage!(x, QuadraticFunction(H, b), λ, ω, weights!, method) end ######################################################################## -# indElem --- index in a lower triangular part of the matrix -# can be obtained by tril_index -function estimate(::ReducedNormal, X, Y, indElem; +# ind_elem --- index in a lower triangular part of the matrix +# can be obtained by tril_index +function estimate(::ReducedNormal, X, Y, ind_elem; Sx::Union{Symmetric, Nothing}=nothing, Sy::Union{Symmetric, Nothing}=nothing, H::Union{Symmetric, Nothing}=nothing, @@ -228,26 +228,26 @@ function estimate(::ReducedNormal, X, Y, indElem; b = svec(Sx) - svec(Sy) # first stage - λ = _penaltyLevel(p) + λ = _penalty_level(p) if Δ === nothing && suppΔ === nothing - Δ = reducedDiffEstimation(H, b, Sx, Sy, X, Y, λ) - suppΔ = getReducedSupport(Δ) + Δ = reduced_diff_estimation(H, b, Sx, Sy, X, Y, λ) + suppΔ = get_reduced_support(Δ) end if suppΔ === nothing - suppΔ = getReducedSupport(Δ) + suppΔ = get_reduced_support(Δ) end # second stage - if ω === nothing && suppω === nothing - ω = invHessianReduced(H, indElem, Sx, Sy, X, Y, λ) - suppω = getReducedSupport(ω) + if ω === nothing && suppω === nothing + ω = inv_hessian_reduced(H, ind_elem, Sx, Sy, X, Y, λ) + suppω = get_reduced_support(ω) end if suppω === nothing - suppω = getReducedSupport(ω) + suppω = get_reduced_support(ω) end # refit stage - indS = (1:size(H,1))[suppΔ .| suppω] - estimate(ReducedOracleNormal(), Sx, Sy, X, Y, indElem, indS), Δ, suppΔ, ω, suppω, indS + ind_s = (1:size(H,1))[suppΔ .| suppω] + estimate(ReducedOracleNormal(), Sx, Sy, X, Y, ind_elem, ind_s), Δ, suppΔ, ω, suppω, ind_s end diff --git a/src/skron.jl b/src/skron.jl index d1517d5..edbfdf3 100644 --- a/src/skron.jl +++ b/src/skron.jl @@ -5,75 +5,75 @@ # Q xx' skron yy' Qt function skron( - x::AbstractVector, - y::AbstractVector, - indRow::AbstractVector, - indCol::AbstractVector) - - nRow = length(indRow) - nCol = length(indCol) - out = zeros( nRow, nCol ) - for ri=1:nRow - for ci=1:nCol - out[ri, ci] = skron(x, y, indRow[ri], indCol[ci]) - end + x::AbstractVector, + y::AbstractVector, + ind_row::AbstractVector, + ind_col::AbstractVector) + + n_row = length(ind_row) + n_col = length(ind_col) + out = zeros( n_row, n_col ) + for ri=1:n_row + for ci=1:n_col + out[ri, ci] = skron(x, y, ind_row[ri], ind_col[ci]) + end end out end function skron( - x::AbstractVector, - y::AbstractVector, - indRow::Int, - indCol::Int + x::AbstractVector, + y::AbstractVector, + ind_row::Int, + ind_col::Int ) p = length(x) - row1, col1 = tril_cartesian(p, indRow) - row2, col2 = tril_cartesian(p, indCol) + row1, col1 = tril_cartesian(p, ind_row) + row2, col2 = tril_cartesian(p, ind_col) vrow = (row1 == col1) ? x[row1] * y[row1] : (x[row1] * y[col1] + y[row1] * x[col1]) / sqrt(2.) vcol = (row2 == col2) ? x[row2] * y[row2] : (x[row2] * y[col2] + y[row2] * x[col2]) / sqrt(2.) - vrow * vcol + vrow * vcol end -# +# # Q Sx skron Sy Qt function skron( - Sx::Symmetric, - Sy::Symmetric) + Sx::Symmetric, + Sy::Symmetric) - p = size(Sx, 1) + p = size(Sx, 1) return skron(Sx, Sy, 1:div(p*(p+1),2)) end function skron( - Sx::Symmetric, + Sx::Symmetric, Sy::Symmetric, - indRow::AbstractVector, - indCol::AbstractVector) - - nRow = length(indRow) - nCol = length(indCol) - out = zeros( nRow, nCol ) - for ri=1:nRow - for ci=1:nCol - out[ri, ci] = skron(Sx, Sy, indRow[ri], indCol[ci]) - end + ind_row::AbstractVector, + ind_col::AbstractVector) + + n_row = length(ind_row) + n_col = length(ind_col) + out = zeros( n_row, n_col ) + for ri=1:n_row + for ci=1:n_col + out[ri, ci] = skron(Sx, Sy, ind_row[ri], ind_col[ci]) + end end out end function skron( - Sx::Symmetric, + Sx::Symmetric, Sy::Symmetric, I::AbstractVector) m = length(I) - + out = zeros( m, m ) for ci = 1:m @@ -88,8 +88,8 @@ function skron( end function skron( - Sx::Symmetric, - Sy::Symmetric, + Sx::Symmetric, + Sy::Symmetric, ind_row::Int, ind_col::Int ) @@ -108,37 +108,37 @@ function skron( if row2 == col2 return ( Sx[row1, row2] * Sy[col1, row2] + Sx[col1, row2] * Sy[row1, row2] ) / sqrt(2.) else - return (Sx[row1, row2] * Sy[col1, col2] + - Sx[row1, col2] * Sy[col1, row2] + - Sx[col1, row2] * Sy[row1, col2] + + return (Sx[row1, row2] * Sy[col1, col2] + + Sx[row1, col2] * Sy[col1, row2] + + Sx[col1, row2] * Sy[row1, col2] + Sx[col1, col2] * Sy[row1, row2] ) / 2. end end end -# +# # Q Sx skron yy' Qt function skron( - Sx::Symmetric, + Sx::Symmetric, y::AbstractVector, - indRow::AbstractVector, - indCol::AbstractVector) - - nRow = length(indRow) - nCol = length(indCol) - out = zeros( nRow, nCol ) - for ri=1:nRow - for ci=1:nCol - out[ri, ci] = skron(Sx, y, indRow[ri], indCol[ci]) - end + ind_row::AbstractVector, + ind_col::AbstractVector) + + n_row = length(ind_row) + n_col = length(ind_col) + out = zeros( n_row, n_col ) + for ri=1:n_row + for ci=1:n_col + out[ri, ci] = skron(Sx, y, ind_row[ri], ind_col[ci]) + end end out end function skron( - Sx::Symmetric, - y::AbstractVector, + Sx::Symmetric, + y::AbstractVector, ind_row::Int, ind_col::Int ) @@ -166,9 +166,9 @@ function skron( if row2 == col2 return ( Sx[row1, row2] * y[col1] * y[row2] + Sx[col1, row2] * y[row1] * y[row2] ) / sqrt(2.) else - return (Sx[row1, row2] * y[col1] * y[col2] + - Sx[row1, col2] * y[col1] * y[row2] + - Sx[col1, row2] * y[row1] * y[col2] + + return (Sx[row1, row2] * y[col1] * y[col2] + + Sx[row1, col2] * y[col1] * y[row2] + + Sx[col1, row2] * y[row1] * y[col2] + Sx[col1, col2] * y[row1] * y[row2] ) / 2. end end @@ -180,12 +180,12 @@ end # the entries (row, col) of the lower triangle of a p × p matrix, in the order # of the coordinates of svec -_trilEntries(p::Int) = [tril_cartesian(p, k) for k=1:div(p * (p + 1), 2)] +_tril_entries(p::Int) = [tril_cartesian(p, k) for k=1:div(p * (p + 1), 2)] # Q(S) == svec(S) function svec(S::Symmetric) - p = size(S, 1) + p = size(S, 1) return svec(S, 1:div(p*(p+1),2)) end @@ -201,14 +201,14 @@ function svec( out end -svec(S::Symmetric, indRow::Int) = svec(S, tril_cartesian(size(S, 1), indRow)) +svec(S::Symmetric, ind_row::Int) = svec(S, tril_cartesian(size(S, 1), ind_row)) svec(S::Symmetric, (row, col)::Tuple{Int, Int}) = (row == col) ? S[row, row] : S[row, col] * sqrt(2.) function svec( - x::AbstractVector, - y::AbstractVector, + x::AbstractVector, + y::AbstractVector, I::AbstractVector) m = length(I) @@ -220,7 +220,7 @@ function svec( end -svec(x::AbstractVector, y::AbstractVector, indRow::Int) = svec(x, y, tril_cartesian(length(x), indRow)) +svec(x::AbstractVector, y::AbstractVector, ind_row::Int) = svec(x, y, tril_cartesian(length(x), ind_row)) function svec( x::AbstractVector, @@ -233,4 +233,4 @@ function svec( else return (x[row] * y[col] + y[row] * x[col]) / sqrt(2.) end -end \ No newline at end of file +end diff --git a/src/support.jl b/src/support.jl index a17e216..6f2624a 100644 --- a/src/support.jl +++ b/src/support.jl @@ -5,7 +5,7 @@ abstract type DiffPrecSupport end Estimate the support of `Δ` by thresholding each entry at `τ` times its standard error from the normal approximation ([`SymmetricOracleNormal`](@ref)). -See [`supportEstimate`](@ref). +See [`support_estimate`](@ref). """ struct ANTSupport <: DiffPrecSupport end @@ -14,7 +14,7 @@ struct ANTSupport <: DiffPrecSupport end Estimate the support of `Δ` by thresholding each entry at `τ` times its bootstrap standard error ([`SymmetricOracleBoot`](@ref)). See -[`supportEstimate`](@ref). +[`support_estimate`](@ref). """ struct BootStdSupport <: DiffPrecSupport end @@ -30,7 +30,7 @@ struct BootMaxSupport <: DiffPrecSupport end # not implemented yet Estimate `Δ` along a path of penalty parameters with the ℓ1 penalized direct difference (D-trace) estimator and choose the penalty on validation data. See -[`supportEstimate`](@ref). +[`support_estimate`](@ref). """ struct DTraceValidationSupport <: DiffPrecSupport end @@ -44,41 +44,41 @@ struct DTraceValidationSupport <: DiffPrecSupport end # For each entry of the lower triangle of Δ, in the order of tril_index, # the indices into svec(Δ) of the union of the supports of the first stage estimate # of Δ and of the second stage estimate of the row of the inverse Hessian. -function __initSupport( +function __init_support( Sx::Symmetric, Sy::Symmetric, X, Y) p = size(X, 2) np = div((p + 1)*p, 2) - estimSupp = Vector{Vector{Int64}}(undef, np) + estim_supp = Vector{Vector{Int64}}(undef, np) # first stage - λ = _penaltyLevel(p) + λ = _penalty_level(p) H = Symmetric( skron(Sx, Sy) ) b = svec(Sx) - svec(Sy) - x1 = reducedDiffEstimation(H, b, Sx, Sy, X, Y, λ) - S1 = getReducedSupport(x1) + x1 = reduced_diff_estimation(H, b, Sx, Sy, X, Y, λ) + S1 = get_reduced_support(x1) # second stage ind = 0 for ci=1:p for ri=ci:p - indElem = tril_index(p, ri, ci) - x2 = invHessianReduced(H, indElem, Sx, Sy, X, Y, λ) - S2 = getReducedSupport(x2) + ind_elem = tril_index(p, ri, ci) + x2 = inv_hessian_reduced(H, ind_elem, Sx, Sy, X, Y, λ) + S2 = get_reduced_support(x2) ind = ind + 1 - estimSupp[ind] = (1:np)[S1 .| S2] + estim_supp[ind] = (1:np)[S1 .| S2] end end - estimSupp + estim_supp end # linear indices into a p × p matrix of a support given as a symmetric p × p # Bool matrix, or as indices into svec(Δ) -_linearSupport(row, col, S::AbstractMatrix{Bool}, p) = getLinearSupport(row, col, S) +_linear_support(row, col, S::AbstractMatrix{Bool}, p) = get_linear_support(row, col, S) -function _linearSupport(row, col, S::AbstractVector{<:Integer}, p) +function _linear_support(row, col, S::AbstractVector{<:Integer}, p) supp = Vector{Int64}() for k in S r, c = tril_cartesian(p, k) @@ -87,7 +87,7 @@ function _linearSupport(row, col, S::AbstractVector{<:Integer}, p) push!(supp, (r-1)*p + c) end end - makeFirst!(supp, p, row, col) + make_first!(supp, p, row, col) end function __threshold(Δpoint::Vector{DiffPrecResultNormal}, p::Int, τ0::Float64) @@ -129,124 +129,130 @@ function __threshold(Δpoint::Vector{DiffPrecResultBoot}, p::Int, τ0::Float64) end -function __getΔ(::ANTSupport, X, Y; estimSupport::Union{AbstractVector,Nothing}=nothing) +function _get_delta(::ANTSupport, X, Y; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) nx, p = size(X) ny = size(Y, 1) Sx = Symmetric( cov(X) ) Sy = Symmetric( cov(Y) ) - eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport + e_s = estim_support === nothing ? __init_support(Sx, Sy, X, Y) : estim_support Δpoint = Array{DiffPrecResultNormal}(undef, div((p+1)*p, 2)) it = 0 for col=1:p for row=col:p it = it + 1 - indS = _linearSupport(row, col, eS[it], p) + ind_s = _linear_support(row, col, e_s[it], p) - Δpoint[it] = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS) + Δpoint[it] = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, ind_s) end end - Δpoint, eS + Δpoint, e_s end """ - supportEstimate(method, X, Y, τ; estimSupport=nothing) + support_estimate(method, X, Y, τ; estim_support=nothing) Estimate the support of `Δ = Ωx - Ωy` with `method` one of [`ANTSupport`](@ref) or [`BootStdSupport`](@ref): every entry in the lower triangle is refitted and kept when it exceeds `τ` times its standard error. -Returns `(Δsupp, Δpoint, eS)`: the thresholded estimate as a sparse `p × p` +Returns `(Δsupp, Δpoint, e_s)`: the thresholded estimate as a sparse `p × p` matrix (a vector of them when `τ` is a vector of thresholds), the refitted entries in the order of `ProximalBase.tril_index`, and the supports used for the refits. -`estimSupport` gives the support used to refit each entry, as computed by the +`estim_support` gives the support used to refit each entry, as computed by the first and second stage estimators unless given: a vector with one element per entry of the lower triangle, either a vector of indices into `svec(Δ)` or a symmetric `p × p` Bool matrix. - supportEstimate(DTraceValidationSupport(), Sx, Sy, SxValid, SyValid, Λarr, method=CDSolver()) + support_estimate(DTraceValidationSupport(), Sx, Sy, SxValid, SyValid, Λarr, method=CDSolver()) Compute the ℓ1 penalized direct difference estimator from the covariances `Sx`, `Sy` for each penalty parameter in `Λarr` (in decreasing order) and evaluate the loss `Sx⋅Δ⋅Sy - Sy + Sx` on the validation covariances `SxValid`, `SyValid`. -Returns `(eΔarr, i2, iInf, loss2arr, lossInfarr)`: the estimates as sparse +Returns `(eΔarr, i2, i_inf, loss2arr, loss_infarr)`: the estimates as sparse matrices, the indices of the penalty parameters that minimize the Frobenius and the maximum norm of the validation loss, and the two losses. """ -function supportEstimate(::ANTSupport, X, Y, τ0::Float64; estimSupport::Union{AbstractVector,Nothing}=nothing) +function support_estimate(::ANTSupport, X, Y, τ0::Float64; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) p = size(X, 2) - Δpoint, eS = __getΔ(ANTSupport(), X, Y; estimSupport=estimSupport) + Δpoint, e_s = _get_delta(ANTSupport(), X, Y; estim_support=estim_support) Δsupp = __threshold(Δpoint, p, τ0) - Δsupp, Δpoint, eS + Δsupp, Δpoint, e_s end -function supportEstimate(::ANTSupport, X, Y, τArr::Vector{Float64}; estimSupport::Union{AbstractVector,Nothing}=nothing) +function support_estimate(::ANTSupport, X, Y, τArr::Vector{Float64}; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) p = size(X, 2) - Δpoint, eS = __getΔ(ANTSupport(), X, Y; estimSupport=estimSupport) + Δpoint, e_s = _get_delta(ANTSupport(), X, Y; estim_support=estim_support) Δsupp = Array{SparseMatrixCSC{Float64,Int64}}(undef, length(τArr)) for i=1:length(τArr) Δsupp[i] = __threshold(Δpoint, p, τArr[i]) end - Δsupp, Δpoint, eS + Δsupp, Δpoint, e_s end -function __getΔ(::BootStdSupport, X, Y; estimSupport::Union{AbstractVector,Nothing}=nothing) +function _get_delta(::BootStdSupport, X, Y; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) nx, p = size(X) ny = size(Y, 1) Sx = Symmetric( cov(X) ) Sy = Symmetric( cov(Y) ) - eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport + e_s = estim_support === nothing ? __init_support(Sx, Sy, X, Y) : estim_support Δpoint = Array{DiffPrecResultBoot}(undef, div((p+1)*p, 2)) it = 0 for col=1:p for row=col:p it = it + 1 - indS = _linearSupport(row, col, eS[it], p) + ind_s = _linear_support(row, col, e_s[it], p) - Δpoint[it] = estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, row, col, indS) + Δpoint[it] = estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, row, col, ind_s) end end - Δpoint, eS + Δpoint, e_s end -function supportEstimate(::BootStdSupport, X, Y, τ0::Float64; estimSupport::Union{AbstractVector,Nothing}=nothing) +function support_estimate(::BootStdSupport, X, Y, τ0::Float64; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) p = size(X, 2) - Δpoint, eS = __getΔ(BootStdSupport(), X, Y; estimSupport=estimSupport) + Δpoint, e_s = _get_delta(BootStdSupport(), X, Y; estim_support=estim_support) Δsupp = __threshold(Δpoint, p, τ0) - Δsupp, Δpoint, eS + Δsupp, Δpoint, e_s end -function supportEstimate(::BootStdSupport, X, Y, τArr::Vector{Float64}; estimSupport::Union{AbstractVector,Nothing}=nothing) +function support_estimate(::BootStdSupport, X, Y, τArr::Vector{Float64}; estim_support=nothing, estimSupport=nothing) + estim_support = _deprecated_keyword(estim_support, estimSupport, :estim_support, :estimSupport, nothing) p = size(X, 2) - Δpoint, eS = __getΔ(BootStdSupport(), X, Y; estimSupport=estimSupport) + Δpoint, e_s = _get_delta(BootStdSupport(), X, Y; estim_support=estim_support) Δsupp = Array{SparseMatrixCSC{Float64,Int64}}(undef, length(τArr)) for i=1:length(τArr) Δsupp[i] = __threshold(Δpoint, p, τArr[i]) end - Δsupp, Δpoint, eS + Δsupp, Δpoint, e_s end #################################### @@ -256,7 +262,7 @@ end #################################### # Λarr should be in descending order -function supportEstimate(::DTraceValidationSupport, +function support_estimate(::DTraceValidationSupport, Sx::Symmetric, Sy::Symmetric, SxValid::Symmetric, SyValid::Symmetric, Λarr::Vector{Float64}, @@ -264,21 +270,21 @@ function supportEstimate(::DTraceValidationSupport, numΛ = length(Λarr) loss2arr = Vector{Float64}(undef, numΛ) - lossInfarr = Vector{Float64}(undef, numΛ) + loss_infarr = Vector{Float64}(undef, numΛ) # each λ warm started from the solution at the previous one function validate(i, λ, x, m) - loss2arr[i] = CovarianceSelection.diffLoss(SxValid, x, SyValid, 2) - lossInfarr[i] = CovarianceSelection.diffLoss(SxValid, x, SyValid, Inf) + loss2arr[i] = CovarianceSelection.diff_loss(SxValid, x, SyValid, 2) + loss_infarr[i] = CovarianceSelection.diff_loss(SxValid, x, SyValid, Inf) false end - path = solve_path(directDifferenceLoss(Sx, Sy), NormL1(1.0), Λarr, method; + path = solve_path(direct_difference_loss(Sx, Sy), NormL1(1.0), Λarr, method; x=SymmetricSparseIterate(size(Sx, 1)), screening=false, callback=validate) eΔarr = path.solutions # find min loss i2 = argmin(loss2arr) - iInf = argmin(lossInfarr) + i_inf = argmin(loss_infarr) - eΔarr, i2, iInf, loss2arr, lossInfarr + eΔarr, i2, i_inf, loss2arr, loss_infarr end diff --git a/src/ustats.jl b/src/ustats.jl index 2a289e1..c95f319 100644 --- a/src/ustats.jl +++ b/src/ustats.jl @@ -67,7 +67,7 @@ function _compute(::OneSampleUstatistics{2}, h, X::AbstractMatrix) v = 0. @inbounds for i1=1:n for i2=i1+1:n - v += h( view(X, i1, :), + v += h( view(X, i1, :), view(X, i2, :) ) end end @@ -135,7 +135,7 @@ function _variance(::TwoSampleUstatistics{1, 1}, h, X::AbstractMatrix, Y::Abstra nx = size(X, 1) ny = size(Y, 1) - q = zeros(nx) + q = zeros(nx) @inbounds for k=1:nx vx = view(X, k, :) for j=1:ny @@ -158,34 +158,34 @@ function _variance(::TwoSampleUstatistics{1, 1}, h, X::AbstractMatrix, Y::Abstra σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2 σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2 - return σ1/nx + σ2/ny + return σ1/nx + σ2/ny end ###################################### -### examples +### examples -function createVarianceStatistic() - @inline function h(x, y) - (x - y)^2 / 2. +function create_variance_statistic() + @inline function h(x, y) + (x - y)^2 / 2. end OneSampleUstatistics{2}(h) end -function createCovarianceStatistic(row::Int, col::Int) - @inline function h(X, Y) - @inbounds (X[row] - Y[row]) * (X[col] - Y[col]) / 2. +function create_covariance_statistic(row::Int, col::Int) + @inline function h(X, Y) + @inbounds (X[row] - Y[row]) * (X[col] - Y[col]) / 2. end OneSampleUstatistics{2}(h) end -function createSecondMomentStatistic() +function create_second_moment_statistic() h(X) = X*X OneSampleUstatistics{1}(h) end -function createSecondMomentStatistic(row::Int, col::Int) +function create_second_moment_statistic(row::Int, col::Int) @inbounds h(X) = X[row] * X[col] OneSampleUstatistics{1}(h) -end \ No newline at end of file +end diff --git a/src/util.jl b/src/util.jl index 9924c46..6e128da 100644 --- a/src/util.jl +++ b/src/util.jl @@ -1,8 +1,8 @@ # quantile function of the standard normal distribution -_normQuantile(α::Real) = -sqrt(2.) * erfcinv(2. * α) +_norm_quantile(α::Real) = -sqrt(2.) * erfcinv(2. * α) # penalty level used by the first and second stage estimators with p variables -_penaltyLevel(p::Integer) = 1.01 * _normQuantile(1. - 0.1 / (p * (p + 1))) +_penalty_level(p::Integer) = 1.01 * _norm_quantile(1. - 0.1 / (p * (p + 1))) # compute kron(A, B)[Irow, Icol] @@ -61,7 +61,7 @@ skron_sub(A, B, Irow, Icol) = skron_sub!(Matrix{Float64}(undef, length(Irow), le -function _getElemKron(A::AbstractMatrix, B::AbstractMatrix, row, col) +function _get_elem_kron(A::AbstractMatrix, B::AbstractMatrix, row, col) p, q = size(B) ac = div(col-1, q) + 1 @@ -72,7 +72,7 @@ function _getElemKron(A::AbstractMatrix, B::AbstractMatrix, row, col) @inbounds A[ar, ac] * B[br, bc] end -function _getElemKron(S::AbstractMatrix, xi::AbstractVector, row, col) +function _get_elem_kron(S::AbstractMatrix, xi::AbstractVector, row, col) p = size(S, 1) ac = div(col-1, p) + 1 @@ -83,7 +83,7 @@ function _getElemKron(S::AbstractMatrix, xi::AbstractVector, row, col) @inbounds S[ar, ac] * xi[br] * xi[bc] end -function _getElemKron(xi::AbstractVector, S::AbstractMatrix, row, col) +function _get_elem_kron(xi::AbstractVector, S::AbstractMatrix, row, col) p = size(S, 1) ac = div(col-1, p) + 1 @@ -94,7 +94,7 @@ function _getElemKron(xi::AbstractVector, S::AbstractMatrix, row, col) @inbounds S[br, bc] * xi[ar] * xi[ac] end -function _getElemSKron(S::AbstractMatrix, xi::AbstractVector, row, col) +function _get_elem_skron(S::AbstractMatrix, xi::AbstractVector, row, col) p = size(S, 1) ac = div(col-1, p) + 1 @@ -105,35 +105,35 @@ function _getElemSKron(S::AbstractMatrix, xi::AbstractVector, row, col) @inbounds ( S[ar, ac] * xi[br] * xi[bc] + S[br, bc] * xi[ar] * xi[ac] ) / 2. end -_getElemSKron(A::AbstractMatrix, B::AbstractMatrix, row, col) = - ( _getElemKron(A, B, row, col) + _getElemKron(B, A, row, col) ) / 2. +_get_elem_skron(A::AbstractMatrix, B::AbstractMatrix, row, col) = + ( _get_elem_kron(A, B, row, col) + _get_elem_kron(B, A, row, col) ) / 2. -# computes variance of Var(Sx - Sy)[indS, indS] -# where indS is a list of elements in [1:p, 1:p] -function _varSxSy(Sx, nx, Sy, ny, indS) +# computes variance of Var(Sx - Sy)[ind_s, ind_s] +# where ind_s is a list of elements in [1:p, 1:p] +function _var_sx_sy(Sx, nx, Sy, ny, ind_s) - varS = zeros(length(indS), length(indS)) # var(Sx) + var(Sy) + var_s = zeros(length(ind_s), length(ind_s)) # var(Sx) + var(Sy) I = CartesianIndices(Sx) - for ia = 1:length(indS) - a = indS[ia] - for ib = 1:length(indS) - b = indS[ib] + for ia = 1:length(ind_s) + a = ind_s[ia] + for ib = 1:length(ind_s) + b = ind_s[ib] i, j = Tuple( I[a] ) k, l = Tuple( I[b] ) - varS[ia, ib] = (Sx[i, k] * Sx[j, l] + Sx[i, l] * Sx[j, k]) / (nx - 1) + (Sy[i, k] * Sy[j, l] + Sy[i, l] * Sy[j, k]) / (ny - 1) + var_s[ia, ib] = (Sx[i, k] * Sx[j, l] + Sx[i, l] * Sx[j, k]) / (nx - 1) + (Sy[i, k] * Sy[j, l] + Sy[i, l] * Sy[j, k]) / (ny - 1) end end - varS + var_s end # entries of a symmetric x larger than τ in absolute value, as a p × p BitMatrix -function getSupport(x::SymmetricSparseIterate, τ::Float64=1e-3) +function get_support(x::SymmetricSparseIterate, τ::Float64=1e-3) p = size(x, 1) S = falses(p, p) @@ -150,7 +150,7 @@ end # x holds either the lower triangle of a symmetric p × p matrix (length p(p+1)/2), # or a p × p matrix stored by columns (length p²) -function getSupport(x::SparseIterate, p::Int, τ::Float64=1e-3) +function get_support(x::SparseIterate, p::Int, τ::Float64=1e-3) S = falses(p, p) if length(x) == p * p for ind=1:p*p @@ -175,9 +175,9 @@ function getSupport(x::SparseIterate, p::Int, τ::Float64=1e-3) S end -getReducedSupport(x::SparseIterate, τ::Float64=1e-3) = abs.(x) .> τ +get_reduced_support(x::SparseIterate, τ::Float64=1e-3) = abs.(x) .> τ -function getLinearSupport(row::Int, col::Int, S1::AbstractMatrix{Bool}, S2::AbstractMatrix{Bool}) +function get_linear_support(row::Int, col::Int, S1::AbstractMatrix{Bool}, S2::AbstractMatrix{Bool}) p = size(S1, 1) supp = Vector{Int64}() @@ -188,10 +188,10 @@ function getLinearSupport(row::Int, col::Int, S1::AbstractMatrix{Bool}, S2::Abst end end end - makeFirst!(supp, p, row, col) + make_first!(supp, p, row, col) end -function getLinearSupport(row::Int, col::Int, S::AbstractMatrix{Bool}) +function get_linear_support(row::Int, col::Int, S::AbstractMatrix{Bool}) p = size(S, 1) supp = Vector{Int64}() @@ -202,11 +202,11 @@ function getLinearSupport(row::Int, col::Int, S::AbstractMatrix{Bool}) end end end - makeFirst!(supp, p, row, col) + make_first!(supp, p, row, col) end # moves (row, col) to the front of supp and (col, row) to the second position -function makeFirst!(supp::Vector{Int64}, p, row, col) +function make_first!(supp::Vector{Int64}, p, row, col) ind = (col - 1) * p + row pos = findfirst(isequal(ind), supp) if pos === nothing @@ -228,10 +228,10 @@ function makeFirst!(supp::Vector{Int64}, p, row, col) end -function makeFirst!(supp::Vector{Int64}, indElem) - pos = findfirst(isequal(indElem), supp) +function make_first!(supp::Vector{Int64}, ind_elem) + pos = findfirst(isequal(ind_elem), supp) if pos === nothing - pushfirst!(supp, indElem) + pushfirst!(supp, ind_elem) else supp[pos], supp[1] = supp[1], supp[pos] end @@ -240,45 +240,45 @@ end """ - computeSimulationResult(res, trueParam, α=0.05) + compute_simulation_result(res, true_param, α=0.05) -Summarize the estimates `res` of one entry of `Δ` with true value `trueParam` +Summarize the estimates `res` of one entry of `Δ` with true value `true_param` over repetitions of a simulation as a [`SimulationResult`](@ref). The confidence intervals at level `1 - α` are `p ± z₁₋α/₂ std` for a vector of [`DiffPrecResultNormal`](@ref), and the `α/2` and `1 - α/2` quantiles of the bootstrap replicates for a vector of [`DiffPrecResultBoot`](@ref). """ -function computeSimulationResult(res::Vector{DiffPrecResultNormal}, trueParam::Float64, α::Float64=0.05) +function compute_simulation_result(res::Vector{DiffPrecResultNormal}, true_param::Float64, α::Float64=0.05) Eestim = 0. bias = 0. coverage = 0 - lenCoverage = 0. + len_coverage = 0. num_rep = length(res) - zα = _normQuantile(1. - α / 2.) + zα = _norm_quantile(1. - α / 2.) for rep=1:num_rep - if res[rep].p - zα * res[rep].std < trueParam && res[rep].p + zα * res[rep].std > trueParam + if res[rep].p - zα * res[rep].std < true_param && res[rep].p + zα * res[rep].std > true_param coverage += 1 end - lenCoverage += 2. * zα * res[rep].std + len_coverage += 2. * zα * res[rep].std Eestim += res[rep].p end - bias = Eestim / num_rep - trueParam + bias = Eestim / num_rep - true_param ep = [res[rep].p for rep=1:num_rep] mean_p = mean(ep) var_p = var(ep) - SimulationResult(bias, coverage / num_rep, lenCoverage / num_rep, mean_p, var_p) + SimulationResult(bias, coverage / num_rep, len_coverage / num_rep, mean_p, var_p) end -function computeSimulationResult(res::Vector{DiffPrecResultBoot}, trueParam::Float64, α::Float64=0.05) +function compute_simulation_result(res::Vector{DiffPrecResultBoot}, true_param::Float64, α::Float64=0.05) Eestim = 0. bias = 0. coverage = 0 - lenCoverage = 0. + len_coverage = 0. num_rep = length(res) z1 = α / 2. @@ -286,19 +286,19 @@ function computeSimulationResult(res::Vector{DiffPrecResultBoot}, trueParam::Flo for rep=1:num_rep lb, ub = quantile(res[rep].boot_p, [z1, z2]) - if lb < trueParam && ub > trueParam + if lb < true_param && ub > true_param coverage += 1 end - lenCoverage += (ub - lb) + len_coverage += (ub - lb) Eestim += res[rep].p end - bias = Eestim / num_rep - trueParam + bias = Eestim / num_rep - true_param ep = [res[rep].p for rep=1:num_rep] mean_p = mean(ep) var_p = var(ep) - SimulationResult(bias, coverage / num_rep, lenCoverage / num_rep, mean_p, var_p) + SimulationResult(bias, coverage / num_rep, len_coverage / num_rep, mean_p, var_p) end diff --git a/src/variance.jl b/src/variance.jl index d8400ac..e9de2eb 100644 --- a/src/variance.jl +++ b/src/variance.jl @@ -27,7 +27,7 @@ function variance( t = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds t += _getElemSKron(Sx, Sy, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds t += _get_elem_skron(Sx, Sy, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) @@ -39,7 +39,7 @@ function variance( v = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds v += _getElemSKron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds v += _get_elem_skron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) @@ -54,7 +54,7 @@ function variance( v = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds v += _getElemSKron(Sx, view(Y, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds v += _get_elem_skron(Sx, view(Y, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) @@ -104,7 +104,7 @@ function variance( end end for ri=1:length(indω) - @inbounds t -= ω[ri] * (svec(Sy, indω[ri]) - svec(Sx, indω[ri])) + @inbounds t -= ω[ri] * (svec(Sy, indω[ri]) - svec(Sx, indω[ri])) end # the entries (row, col) of the coordinates, computed once for the loops below @@ -177,7 +177,7 @@ function variance( t = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds t += _getElemKron(Sy, Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds t += _get_elem_kron(Sy, Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) @@ -189,7 +189,7 @@ function variance( v = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds v += _getElemKron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds v += _get_elem_kron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) @@ -204,7 +204,7 @@ function variance( v = 0. for ci=1:length(indΔ) for ri=1:length(indω) - @inbounds v += _getElemKron(view(Y, k, :), Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] + @inbounds v += _get_elem_kron(view(Y, k, :), Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci] end end for ri=1:length(indω) diff --git a/test/bootstrap.jl b/test/bootstrap.jl index a855c1c..bdfbfcd 100644 --- a/test/bootstrap.jl +++ b/test/bootstrap.jl @@ -12,34 +12,34 @@ include("data.jl") @testset "bootstrap" begin p, n = 5, 500 - X, Y, _ = generateData(p, n) + X, Y, _ = generate_data(p, n) Sx = Symmetric( cov(X) ) Sy = Symmetric( cov(Y) ) np = div(p * (p + 1), 2) # on the full support the estimates are svec(inv(Sx) - inv(Sy)) - eS = [collect(1:np) for _=1:np] - straps, eS1 = bootstrap(X, Y; bootSamples=100, estimSupport=eS, rng=Xoshiro(1)) - @test eS1 === eS + e_s = [collect(1:np) for _=1:np] + straps, eS1 = bootstrap(X, Y; boot_samples=100, estim_support=e_s, rng=Xoshiro(1)) + @test eS1 === e_s @test straps.Δhat ≈ svec(Symmetric(inv(Sx) - inv(Sy))) @test size(straps.Δb) == (np, 100) - @test straps.Δb == bootstrap(X, Y; bootSamples=100, estimSupport=eS, rng=Xoshiro(1))[1].Δb + @test straps.Δb == bootstrap(X, Y; boot_samples=100, estim_support=e_s, rng=Xoshiro(1))[1].Δb # estimated supports - straps, eS = bootstrap(X, Y; bootSamples=100) - @test length(eS) == np - @test all(j -> eS[j][1] == j, 1:np) + straps, e_s = bootstrap(X, Y; boot_samples=100) + @test length(e_s) == np + @test all(j -> e_s[j][1] == j, 1:np) - CI = simulCI(straps, 0.9) + CI = simul_ci(straps, 0.9) @test size(CI) == (np, 2) @test all(CI[:, 1] .< straps.Δhat .< CI[:, 2]) w = CI[:, 2] - CI[:, 1] @test all(w .≈ w[1]) @test CI[:, 1] + CI[:, 2] ≈ 2 * straps.Δhat # wider intervals at a higher level - @test all(simulCI(straps, 0.99)[:, 2] .>= CI[:, 2]) + @test all(simul_ci(straps, 0.99)[:, 2] .>= CI[:, 2]) - CIs = simulCIstudentized(straps, 0.9) + CIs = simul_ci_studentized(straps, 0.9) @test all(CIs[:, 1] .< straps.Δhat .< CIs[:, 2]) sd = vec(std(straps.Δb; dims = 2, corrected = false)) @test (CIs[:, 2] - CIs[:, 1]) ./ sd ≈ fill((CIs[1, 2] - CIs[1, 1]) / sd[1], np) diff --git a/test/data.jl b/test/data.jl index 62dab49..ec0fe80 100644 --- a/test/data.jl +++ b/test/data.jl @@ -1,6 +1,6 @@ # Two Gaussian samples whose precision matrices Ωx and Ωy = Ωx - Δ differ by a # tridiagonal Δ. Returns X, Y and the true difference. -function generateData(p, n) +function generate_data(p, n) # generate model Ωx = Matrix{Float64}(I, p, p) for l=1:p-1 diff --git a/test/deprecated.jl b/test/deprecated.jl new file mode 100644 index 0000000..2d58aca --- /dev/null +++ b/test/deprecated.jl @@ -0,0 +1,30 @@ +using DiffPrecTest, Test, Random, LinearAlgebra + +@testset "naming compatibility" begin + for (old, new) in ((:diffEstimation, :diff_estimation), + (:invAsymHessianEstimation, :inv_asym_hessian_estimation), + (:invQSymHessian, :inv_q_sym_hessian), + (:invHessianReduced, :inv_hessian_reduced), + (:supportEstimate, :support_estimate), + (:simulCI, :simul_ci), (:simulCIstudentized, :simul_ci_studentized), + (:computeSimulationResult, :compute_simulation_result)) + @test Base.isdeprecated(DiffPrecTest, old) + @test getfield(DiffPrecTest, old) === getfield(DiffPrecTest, new) + end + + result = SimulationResult(0., 0.95, 2., 0., 1.) + @test result.len_coverage == 2. + @test_deprecated r"len_coverage" result.lenCoverage + @test :lenCoverage in propertynames(result) + + X = randn(Xoshiro(1), 24, 2) + Y = randn(Xoshiro(2), 24, 2) + supports() = [collect(1:3) for _ in 1:3] + canonical = bootstrap(X, Y; boot_samples=3, estim_support=supports(), rng=Xoshiro(3))[1] + legacy = @test_deprecated r"boot_samples" bootstrap(X, Y; bootSamples=3, estim_support=supports(), rng=Xoshiro(3)) + @test legacy[1].Δb == canonical.Δb + legacy_support = @test_deprecated r"estim_support" bootstrap(X, Y; boot_samples=3, estimSupport=supports(), rng=Xoshiro(3)) + @test legacy_support[1].Δb == canonical.Δb + @test_throws ArgumentError bootstrap(X, Y; boot_samples=3, bootSamples=3) + @test_throws ArgumentError bootstrap(X, Y; estim_support=supports(), estimSupport=supports()) +end diff --git a/test/estimate.jl b/test/estimate.jl index db8f2b2..49da3e4 100644 --- a/test/estimate.jl +++ b/test/estimate.jl @@ -12,7 +12,7 @@ include("data.jl") @testset "oracle estimators" begin p, n = 5, 200 - X, Y, _ = generateData(p, n) + X, Y, _ = generate_data(p, n) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) D = inv(Sx) - inv(Sy) @@ -24,11 +24,11 @@ include("data.jl") for (row, col) in [(1, 1), (3, 2), (2, 3)] @test estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, copy(full)).p ≈ D[row, col] @test estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, row, col, copy(full)).p ≈ D[row, col] - @test estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, row, col, copy(full); bootSamples=2).p ≈ D[row, col] - @test estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, row, col, copy(full); bootSamples=2).p ≈ D[row, col] + @test estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, row, col, copy(full); boot_samples=2).p ≈ D[row, col] + @test estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, row, col, copy(full); boot_samples=2).p ≈ D[row, col] - indElem = tril_index(p, max(row, col), min(row, col)) - @test estimate(ReducedOracleNormal(), Sx, Sy, X, Y, indElem, copy(reducedFull)).p ≈ D[row, col] + ind_elem = tril_index(p, max(row, col), min(row, col)) + @test estimate(ReducedOracleNormal(), Sx, Sy, X, Y, ind_elem, copy(reducedFull)).p ≈ D[row, col] end end @@ -36,11 +36,11 @@ include("data.jl") S = abs.(D) .> 0.1 S[diagind(S)] .= true for (row, col) in [(1, 1), (3, 2), (5, 1)] - res = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, DiffPrecTest.getLinearSupport(row, col, S)) + res = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, DiffPrecTest.get_linear_support(row, col, S)) - indElem = tril_index(p, row, col) - indS = [tril_index(p, r, c) for c=1:p for r=c:p if S[r, c] || (r, c) == (row, col)] - res1 = estimate(ReducedOracleNormal(), Sx, Sy, X, Y, indElem, indS) + ind_elem = tril_index(p, row, col) + ind_s = [tril_index(p, r, c) for c=1:p for r=c:p if S[r, c] || (r, c) == (row, col)] + res1 = estimate(ReducedOracleNormal(), Sx, Sy, X, Y, ind_elem, ind_s) @test res.p ≈ res1.p @test res.std ≈ res1.std @@ -53,16 +53,16 @@ include("data.jl") # is the refitted estimator S = abs.(D) .> 0.1 for (row, col) in [(2, 2), (3, 2)] - indS = DiffPrecTest.getLinearSupport(row, col, S) - res = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS) + ind_s = DiffPrecTest.get_linear_support(row, col, S) + res = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, ind_s) - C = DiffPrecTest.skron_sub(Sx, Sy, indS, indS) - e = zeros(length(indS)) + C = DiffPrecTest.skron_sub(Sx, Sy, ind_s, ind_s) + e = zeros(length(ind_s)) e[1] = 1. ω = C \ e - Δ = C \ (Sy[indS] - Sx[indS]) + Δ = C \ (Sy[ind_s] - Sx[ind_s]) k = row == col ? 1 : 2 - res1 = estimate(DecorrelatedScoreOracleNormal(), Sx, Sy, X, Y, row, col, ω, indS, Δ[k+1:end], indS[k+1:end]) + res1 = estimate(DecorrelatedScoreOracleNormal(), Sx, Sy, X, Y, row, col, ω, ind_s, Δ[k+1:end], ind_s[k+1:end]) @test res.p ≈ res1.p @test res.std ≈ res1.std @@ -70,28 +70,28 @@ include("data.jl") end @testset "bootstrap" begin - indS = DiffPrecTest.getLinearSupport(2, 1, abs.(D) .> 0.1) + ind_s = DiffPrecTest.get_linear_support(2, 1, abs.(D) .> 0.1) for method in (SymmetricOracleBoot(), AsymmetricOracleBoot()) - r1 = estimate(method, Sx, Sy, X, Y, 2, 1, copy(indS); bootSamples=50, rng=Xoshiro(1)) - r2 = estimate(method, Sx, Sy, X, Y, 2, 1, copy(indS); bootSamples=50, rng=Xoshiro(1)) + r1 = estimate(method, Sx, Sy, X, Y, 2, 1, copy(ind_s); boot_samples=50, rng=Xoshiro(1)) + r2 = estimate(method, Sx, Sy, X, Y, 2, 1, copy(ind_s); boot_samples=50, rng=Xoshiro(1)) @test length(r1.boot_p) == 50 @test r1.boot_p == r2.boot_p @test std(r1.boot_p) > 0. end - @test estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, 2, 1, copy(indS); bootSamples=2).p ≈ - estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, 2, 1, copy(indS)).p + @test estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, 2, 1, copy(ind_s); boot_samples=2).p ≈ + estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, 2, 1, copy(ind_s)).p end end @testset "three stage estimators" begin p, n = 6, 2000 - X, Y, Δtrue = generateData(p, n) + X, Y, Δtrue = generate_data(p, n) for (row, col) in [(1, 1), (2, 1), (4, 1)] - res, eΔ, suppΔ, eω, suppω, indS = estimate(SymmetricNormal(), X, Y, row, col) + res, eΔ, suppΔ, eω, suppω, ind_s = estimate(SymmetricNormal(), X, Y, row, col) @test res.std > 0. @test abs(res.p - Δtrue[row, col]) < 4 * res.std - @test indS[1] == (col-1)*p + row + @test ind_s[1] == (col-1)*p + row @test issymmetric(suppΔ) # the order of row and col does not matter @@ -115,12 +115,12 @@ end @test abs(res.p - Δtrue[row, col]) < 4 * res.std @test estimate(SeparateNormal(), X, Y, col, row).p ≈ res.p atol=1e-6 - indElem = tril_index(p, row, col) - res, eΔr, suppΔr, eωr, suppωr, indSr = estimate(ReducedNormal(), X, Y, indElem) + ind_elem = tril_index(p, row, col) + res, eΔr, suppΔr, eωr, suppωr, indSr = estimate(ReducedNormal(), X, Y, ind_elem) @test res.std > 0. @test abs(res.p - Δtrue[row, col]) < 4 * res.std - @test indSr[1] == indElem - @test estimate(ReducedNormal(), X, Y, indElem; Δ=eΔr)[1].p ≈ res.p + @test indSr[1] == ind_elem + @test estimate(ReducedNormal(), X, Y, ind_elem; Δ=eΔr)[1].p ≈ res.p end end @@ -136,7 +136,7 @@ end Y = randn(n, p) / cholesky(Ωy).U' push!(res, estimate(SymmetricNormal(), X, Y, 2, 1)[1]) end - sr = computeSimulationResult(res, Ωx[2, 1] - Ωy[2, 1]) + sr = compute_simulation_result(res, Ωx[2, 1] - Ωy[2, 1]) @test 0.88 <= sr.coverage <= 1. @test abs(sr.bias) < 0.02 end diff --git a/test/invHessian.jl b/test/inv_hessian.jl similarity index 81% rename from test/invHessian.jl rename to test/inv_hessian.jl index ac61935..696f07d 100644 --- a/test/invHessian.jl +++ b/test/inv_hessian.jl @@ -8,7 +8,7 @@ using Random, LinearAlgebra, SparseArrays, Statistics Random.seed!(4) -function genData(p; n = 1000) +function gen_data(p; n = 1000) Sigmay = [0.7^abs(i-j) for i=1:p, j=1:p] X = randn(n,p) Y = randn(n,p) * sqrt(Sigmay) @@ -20,7 +20,7 @@ const tight = CDSolver(; maxpasses=5000, tol=1e-12) @testset "invert Kroneker" begin for rep=1:20 p = 20 - X, Y = genData(p) + X, Y = gen_data(p) hSx, hSy = cov(X), cov(Y) A = kron(hSy, hSx) @@ -33,7 +33,7 @@ const tight = CDSolver(; maxpasses=5000, tol=1e-12) x = SparseIterate(p*p) x1 = SparseIterate(p*p) - f = DiffPrecTest._inverseKroneckerLoss(Symmetric(hSx), Symmetric(hSy), ri, ci) + f = DiffPrecTest._inverse_kronecker_loss(Symmetric(hSx), Symmetric(hSy), ri, ci) b = zeros(Float64, p*p) b[(ci-1)*p+ri] = -1. @@ -49,21 +49,21 @@ end @testset "invert reduced SymKroneker" begin for rep=1:20 p = 10 - X, Y = genData(p) + X, Y = gen_data(p) hSx, hSy = Symmetric(cov(X)), Symmetric(cov(Y)) H = DiffPrecTest.skron(hSx, hSy) m = size(H, 1) λ = 0.05 + 0.25 * rand() # unit weights on svec(Θ) are weights 1/√2 on each off-diagonal entry of Θ - W = DiffPrecTest._symWeights!(zeros(p, p), ones(m), p) + W = DiffPrecTest._sym_weights!(zeros(p, p), ones(m), p) for (ri, ci) in [(1,1), (3,1), (5,4)] ind = tril_index(p, ri, ci) x = SymmetricSparseIterate(p) x1 = SparseIterate(m) - f = DiffPrecTest._inverseSymKroneckerLoss(hSx, hSy, ind) + f = DiffPrecTest._inverse_sym_kronecker_loss(hSx, hSy, ind) b = zeros(m) b[ind] = -1. f1 = QuadraticFunction(H, b) @@ -78,7 +78,7 @@ end @testset "gradient" begin p = 5 - X, Y = genData(p) + X, Y = gen_data(p) hSx, hSy = Symmetric(cov(X)), Symmetric(cov(Y)) # reduced symmetric Hessian: the gradient in svec(Θ) is Hθ - e_ind, which is @@ -93,7 +93,7 @@ end x[r, c] = r == c ? θ[k] : θ[k] / sqrt(2) end @test Vector(DiffPrecTest._svec(x)) ≈ Vector(θ) - G = gradient(DiffPrecTest._inverseSymKroneckerLoss(hSx, hSy, ind), x) + G = gradient(DiffPrecTest._inverse_sym_kronecker_loss(hSx, hSy, ind), x) e = zeros(m) e[ind] = 1. @test [r == c ? G[r, c] : sqrt(2) * G[r, c] for c=1:p for r=c:p] ≈ H * Vector(θ) - e @@ -102,7 +102,7 @@ end A = kron(hSy, hSx) x = SparseIterate(p*p) x .= sprandn(p*p, 0.5) - f = DiffPrecTest._inverseKroneckerLoss(hSx, hSy, 3, 2) + f = DiffPrecTest._inverse_kronecker_loss(hSx, hSy, 3, 2) e = zeros(p*p) e[(2-1)*p + 3] = 1. @test vec(gradient(f, x)) ≈ A * Vector(x) - e @@ -112,7 +112,7 @@ end # with a negligible penalty the second stage estimates are the rows of the # inverse of the sample Hessian p = 5 - X, Y = genData(p; n = 500) + X, Y = gen_data(p; n = 500) n = size(X, 1) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -120,13 +120,13 @@ end H = DiffPrecTest.skron(Sx, Sy) for (ri, ci) in [(1,1), (4,2)] ind = tril_index(p, ri, ci) - ω = invQSymHessian(Sx, Sy, ind, X, Y, 1e-10, tight) + ω = inv_q_sym_hessian(Sx, Sy, ind, X, Y, 1e-10, tight) @test Vector(ω) ≈ inv(H)[:, ind] atol=1e-6 end A = kron(Sy, Sx) for (ri, ci) in [(1,1), (4,2), (2,4)] - ω = invAsymHessianEstimation(Sx, Sy, ri, ci, X, Y, 1e-10, tight) + ω = inv_asym_hessian_estimation(Sx, Sy, ri, ci, X, Y, 1e-10, tight) @test Vector(ω) ≈ inv(A)[:, (ci-1)*p + ri] atol=1e-6 end end @@ -136,11 +136,11 @@ end # Σy by a divides the solution by a² and takes the same number of passes # (powers of 2 keep the scaling exact) p = 5 - X, Y = genData(p) + X, Y = gen_data(p) hSx, hSy = Symmetric(cov(X)), Symmetric(cov(Y)) ind = tril_index(p, 3, 2) - for (loss, x0) in (((Sx, Sy) -> DiffPrecTest._inverseKroneckerLoss(Sx, Sy, 3, 2), () -> SparseIterate(p*p)), - ((Sx, Sy) -> DiffPrecTest._inverseSymKroneckerLoss(Sx, Sy, ind), () -> SymmetricSparseIterate(p))) + for (loss, x0) in (((Sx, Sy) -> DiffPrecTest._inverse_kronecker_loss(Sx, Sy, 3, 2), () -> SparseIterate(p*p)), + ((Sx, Sy) -> DiffPrecTest._inverse_sym_kronecker_loss(Sx, Sy, ind), () -> SymmetricSparseIterate(p))) passes, sols = Int[], Matrix{Float64}[] for a in (0.25, 1.0, 4.0) x = x0() diff --git a/test/reduced.jl b/test/reduced.jl index da3298e..3b1d711 100644 --- a/test/reduced.jl +++ b/test/reduced.jl @@ -11,7 +11,7 @@ using DiffPrecTest: TwoSampleUstatistics, svec, skron Random.seed!(5) include("data.jl") -generateData(p) = generateData(p, 100)[1:2] +generate_data(p) = generate_data(p, 100)[1:2] # minimize x'Ax/2 + b'x + λ ∑ ωᵢ|xᵢ| over x = u - v with u, v ≥ 0 function solve_pen_quad(A, b, λ, ω) @@ -28,23 +28,23 @@ function solve_pen_quad(A, b, λ, ω) end -@testset "_computeVarStep1!" begin - function createVarianceGrad1(Δ::SparseIterate, rowIndex::Int) +@testset "_compute_var_step1!" begin + function create_variance_grad1(Δ::SparseIterate, rowIndex::Int) @inline function h(x, y)::Float64 # skron(vx_j, vy_k)[ri, :] * Δ out = 0. for inz = 1:nnz(Δ) - indColumn = Δ.nzval2ind[inz] + indColumn = Δ.nzval2ind[inz] out += skron(x, y, rowIndex, indColumn) * Δ.nzval[inz] - end + end out + svec(x, x, rowIndex) - svec(y, y, rowIndex) end TwoSampleUstatistics{1, 1}(h) end - function varStep1Alt!(out, X, Y, θ::SparseIterate) + function var_step1_alt!(out, X, Y, θ::SparseIterate) for i=1:length(out) - s = createVarianceGrad1(θ, i) + s = create_variance_grad1(θ, i) out[i] = variance(s, X, Y) end out @@ -52,77 +52,77 @@ end p = 10 m = div(p*(p+1), 2) - X, Y = generateData(p) + X, Y = generate_data(p) n = size(X, 1) Sx = Symmetric( X'X / n ) - Sy = Symmetric( Y'Y / n ) + Sy = Symmetric( Y'Y / n ) out1 = zeros(m) out2 = zeros(m) θ = SparseIterate(m) - DiffPrecTest._computeVarStep1!(out1, Sx, Sy, X, Y, θ); - varStep1Alt!(out2, X, Y, θ); + DiffPrecTest._compute_var_step1!(out1, Sx, Sy, X, Y, θ); + var_step1_alt!(out2, X, Y, θ); @test out1 ≈ sqrt.(out2) θ .= sprandn(m, 0.5) - DiffPrecTest._computeVarStep1!(out1, Sx, Sy, X, Y, θ); - varStep1Alt!(out2, X, Y, θ); + DiffPrecTest._compute_var_step1!(out1, Sx, Sy, X, Y, θ); + var_step1_alt!(out2, X, Y, θ); @test out1 ≈ sqrt.(out2) end -@testset "reducedDiffEstimation" begin +@testset "reduced_diff_estimation" begin - function reducedDiffEstimationAlt(H::Symmetric, b::Vector, Sx, Sy, X, Y, λ) + function reduced_diff_estimation_alt(H::Symmetric, b::Vector, Sx, Sy, X, Y, λ) x = SparseIterate(size(H, 1)) - + ω = Array{eltype(H)}(undef, length(x)) - DiffPrecTest._computeVarStep1!(ω, Sx, Sy, X, Y, x) + DiffPrecTest._compute_var_step1!(ω, Sx, Sy, X, Y, x) xtmp = solve_pen_quad(H, b, λ, ω) for i=1:length(x) x[i] = xtmp[i] end - DiffPrecTest._computeVarStep1!(ω, Sx, Sy, X, Y, x) + DiffPrecTest._compute_var_step1!(ω, Sx, Sy, X, Y, x) solve_pen_quad(H, b, λ, ω) end p = 5 m = div(p*(p+1), 2) - X, Y = generateData(p) + X, Y = generate_data(p) n = size(X, 1) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) H = Symmetric( skron(Sx, Sy) ) b = svec(Sx) - svec(Sy) - λ = DiffPrecTest._penaltyLevel(p) + λ = DiffPrecTest._penalty_level(p) + + out1 = reduced_diff_estimation_alt(H, b, Sx, Sy, X, Y, λ) + out2 = DiffPrecTest.reduced_diff_estimation(H, b, Sx, Sy, X, Y, λ) - out1 = reducedDiffEstimationAlt(H, b, Sx, Sy, X, Y, λ) - out2 = DiffPrecTest.reducedDiffEstimation(H, b, Sx, Sy, X, Y, λ) - @test out1 ≈ Vector(out2) atol=1e-4 end -@testset "_computeVarStep2!" begin - function createVarianceGrad2(Δ::SparseIterate, rowIndex::Int, indElem::Int) +@testset "_compute_var_step2!" begin + function create_variance_grad2(Δ::SparseIterate, rowIndex::Int, ind_elem::Int) @inline function h(x, y)::Float64 # skron(vx_j, vy_k)[ri, :] * Δ out = 0. for inz = 1:nnz(Δ) - indColumn = Δ.nzval2ind[inz] + indColumn = Δ.nzval2ind[inz] out += skron(x, y, rowIndex, indColumn) * Δ.nzval[inz] - end - out + (rowIndex == indElem ? 1. : 0.) + end + out + (rowIndex == ind_elem ? 1. : 0.) end TwoSampleUstatistics{1, 1}(h) end - function varStep2Alt!(out, indElem, X, Y, θ::SparseIterate) + function var_step2_alt!(out, ind_elem, X, Y, θ::SparseIterate) for i=1:length(out) - s = createVarianceGrad2(θ, i, indElem) + s = create_variance_grad2(θ, i, ind_elem) out[i] = variance(s, X, Y) end out @@ -130,63 +130,63 @@ end p = 10 m = div(p*(p+1), 2) - X, Y = generateData(p) + X, Y = generate_data(p) n = size(X, 1) Sx = Symmetric( X'X / n ) - Sy = Symmetric( Y'Y / n ) + Sy = Symmetric( Y'Y / n ) out1 = zeros(m) out2 = zeros(m) for j=1:m θ = SparseIterate(m) - DiffPrecTest._computeVarStep2!(out1, j, Sx, Sy, X, Y, θ); - varStep2Alt!(out2, j, X, Y, θ); + DiffPrecTest._compute_var_step2!(out1, j, Sx, Sy, X, Y, θ); + var_step2_alt!(out2, j, X, Y, θ); @test out1 ≈ sqrt.(out2) θ .= sprandn(m, 0.5) - DiffPrecTest._computeVarStep2!(out1, j, Sx, Sy, X, Y, θ); - varStep2Alt!(out2, j, X, Y, θ); + DiffPrecTest._compute_var_step2!(out1, j, Sx, Sy, X, Y, θ); + var_step2_alt!(out2, j, X, Y, θ); @test out1 ≈ sqrt.(out2) end end -@testset "invHessianReduced" begin - function invHessianReducedAlt(H::Symmetric, indRow, Sx, Sy, X, Y, λ) +@testset "inv_hessian_reduced" begin + function inv_hessian_reduced_alt(H::Symmetric, ind_row, Sx, Sy, X, Y, λ) x = SparseIterate(size(H, 1)) - x[indRow] = 1. / H[indRow, indRow] + x[ind_row] = 1. / H[ind_row, ind_row] local b = zeros(size(H, 1)) - b[indRow] = -1. + b[ind_row] = -1. ω = Array{eltype(H)}(undef, length(x)) - DiffPrecTest._computeVarStep2!(ω, indRow, Sx, Sy, X, Y, x) - ω[indRow] = 0. + DiffPrecTest._compute_var_step2!(ω, ind_row, Sx, Sy, X, Y, x) + ω[ind_row] = 0. xtmp = solve_pen_quad(H, b, λ, ω) for i=1:length(x) x[i] = xtmp[i] end - DiffPrecTest._computeVarStep2!(ω, indRow, Sx, Sy, X, Y, x) - ω[indRow] = 0. + DiffPrecTest._compute_var_step2!(ω, ind_row, Sx, Sy, X, Y, x) + ω[ind_row] = 0. solve_pen_quad(H, b, λ, ω) end p = 5 m = div(p*(p+1), 2) - X, Y = generateData(p) + X, Y = generate_data(p) n = size(X, 1) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) H = Symmetric( skron(Sx, Sy) ) b = svec(Sx) - svec(Sy) - λ = DiffPrecTest._penaltyLevel(p) + λ = DiffPrecTest._penalty_level(p) for j=1:m - out1 = invHessianReducedAlt(H, j, Sx, Sy, X, Y, λ) - out2 = DiffPrecTest.invHessianReduced(H, j, Sx, Sy, X, Y, λ) + out1 = inv_hessian_reduced_alt(H, j, Sx, Sy, X, Y, λ) + out2 = DiffPrecTest.inv_hessian_reduced(H, j, Sx, Sy, X, Y, λ) @test out1 ≈ Vector(out2) atol=1e-4 end @@ -194,17 +194,17 @@ end # with a negligible penalty the estimates are the rows of inv(H) and svec(Δ) opt = CDSolver(; maxpasses=5000, tol=1e-12) for j in [1, 2, 7] - out = DiffPrecTest.invHessianReduced(H, j, Sx, Sy, X, Y, 1e-10, opt) + out = DiffPrecTest.inv_hessian_reduced(H, j, Sx, Sy, X, Y, 1e-10, opt) @test Vector(out) ≈ inv(H)[:, j] atol=1e-6 end - out = DiffPrecTest.reducedDiffEstimation(H, b, Sx, Sy, X, Y, 1e-10, opt) + out = DiffPrecTest.reduced_diff_estimation(H, b, Sx, Sy, X, Y, 1e-10, opt) @test Vector(out) ≈ svec(Symmetric(inv(Sx) - inv(Sy))) atol=1e-6 end @testset "varianceOracle" begin - function createVarianceReducedEstim(ω, I, Δ, J) + function create_variance_reduced_estim(ω, I, Δ, J) function h(x, y)::Float64 dot( (skron(x, y, I, J) * Δ + svec(x, x, I) - svec(y, y, I)), ω ) end @@ -214,22 +214,22 @@ end for rep=1:10 p = 5 m = div(p*(p+1), 2) - X, Y = generateData(p) + X, Y = generate_data(p) n = size(X, 1) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) - indS = 3:9 - C = skron(Sx, Sy, indS) - b = svec(Sy, indS) - svec(Sx, indS) - - tmp = zeros(length(indS)) + ind_s = 3:9 + C = skron(Sx, Sy, ind_s) + b = svec(Sy, ind_s) - svec(Sx, ind_s) + + tmp = zeros(length(ind_s)) tmp[1] = 1. ω = C \ tmp Δ = C \ b - v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, collect(indS), Δ, collect(indS)) + v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, collect(ind_s), Δ, collect(ind_s)) - v2 = variance(createVarianceReducedEstim(ω, indS, Δ, indS), X, Y) + v2 = variance(create_variance_reduced_estim(ω, ind_s, Δ, ind_s), X, Y) @test v ≈ v2 end end diff --git a/test/runtests.jl b/test/runtests.jl index bb294ca..2ebc477 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ tests = [ "ustats", "skron", "variance", - "invHessian", + "inv_hessian", "reduced", "estimate", "support", @@ -21,6 +21,10 @@ tests = [ include(f) end + @testset "deprecated API" begin + run(`$(Base.julia_cmd()) --startup-file=no --depwarn=yes --project=$(dirname(Base.active_project())) $(joinpath(@__DIR__, "deprecated.jl"))`) + end + @testset "Aqua" begin Aqua.test_all(DiffPrecTest) end diff --git a/test/skron.jl b/test/skron.jl index 360ad85..68d1ac3 100644 --- a/test/skron.jl +++ b/test/skron.jl @@ -24,7 +24,7 @@ function Q(p::Int) end Q end - + function Qt(p::Int) _c = 1. / sqrt(2.) rp = div((p+1)*p,2) @@ -61,10 +61,10 @@ end o2 = zeros( num_row_reduced, num_row_reduced ) for ri=1:num_row_reduced for ci=1:num_row_reduced - o2[ri, ci] = skron(x, y, ri, ci) + o2[ri, ci] = skron(x, y, ri, ci) end end - + @test o2 ≈ Q(p) * ( kron(x*x', y*y') + kron(y*y', x*x') ) * Qt(p) / 2. @test skron(x, y, 1:4, 3:6) ≈ (Q(p) * ( kron(x*x', y*y') + kron(y*y', x*x') ) * Qt(p) / 2.)[1:4, 3:6] @@ -84,7 +84,7 @@ end J = 4:7 @test skron(Sx, Sy, I, J) ≈ (Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.)[I, J] - I = [2, 5, 7, 8] + I = [2, 5, 7, 8] @test skron(Sx, Sy, I) ≈ (Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.)[I, I] I = 10:15 @@ -111,10 +111,10 @@ end I = 10:15 @test svec(Sx, I) ≈ (Q(p) * vec( Sx ))[I] - I = [2, 5, 7, 8] - @test svec(Sy, I) ≈ (Q(p) * vec( Sy ))[I] + I = [2, 5, 7, 8] + @test svec(Sy, I) ≈ (Q(p) * vec( Sy ))[I] @test svec(view(X, 1, :), view(Y, 1, :), I) ≈ (Q(p) * kron(view(X, 1, :), view(Y, 1, :)))[I] end -end \ No newline at end of file +end diff --git a/test/support.jl b/test/support.jl index 413b584..4bb8c55 100644 --- a/test/support.jl +++ b/test/support.jl @@ -11,15 +11,15 @@ include("data.jl") @testset "initial support" begin p, n = 5, 500 - X, Y, _ = generateData(p, n) - eS = DiffPrecTest.__initSupport(Symmetric(cov(X)), Symmetric(cov(Y)), X, Y) - @test length(eS) == div(p * (p + 1), 2) - @test all(s -> issorted(s) && all(in(1:length(eS)), s), eS) + X, Y, _ = generate_data(p, n) + e_s = DiffPrecTest.__init_support(Symmetric(cov(X)), Symmetric(cov(Y)), X, Y) + @test length(e_s) == div(p * (p + 1), 2) + @test all(s -> issorted(s) && all(in(1:length(e_s)), s), e_s) end @testset "ANT and BootStd" begin p, n = 5, 1000 - X, Y, Δtrue = generateData(p, n) + X, Y, Δtrue = generate_data(p, n) Sx = Symmetric( cov(X) ) Sy = Symmetric( cov(Y) ) D = inv(Sx) - inv(Sy) @@ -30,18 +30,18 @@ end eSred = [collect(1:np) for _=1:np] eSmat = [trues(p, p) for _=1:np] for method in (ANTSupport(), BootStdSupport()) - Δsupp, Δpoint, eS = supportEstimate(method, X, Y, 0.; estimSupport=eSred) - @test eS === eSred + Δsupp, Δpoint, e_s = support_estimate(method, X, Y, 0.; estim_support=eSred) + @test e_s === eSred @test Matrix(Δsupp) ≈ D @test [Δpoint[tril_index(p, r, c)].p for r=1:p, c=1:p if r >= c] ≈ [D[r, c] for r=1:p, c=1:p if r >= c] - Δsupp1, _, _ = supportEstimate(method, X, Y, 0.; estimSupport=eSmat) + Δsupp1, _, _ = support_estimate(method, X, Y, 0.; estim_support=eSmat) @test Matrix(Δsupp1) ≈ D end # thresholds - Δsupp, Δpoint, eS = supportEstimate(ANTSupport(), X, Y, 3.) - @test eltype(eS) == Vector{Int64} + Δsupp, Δpoint, e_s = support_estimate(ANTSupport(), X, Y, 3.) + @test eltype(e_s) == Vector{Int64} for c=1:p, r=c:p pt = Δpoint[tril_index(p, r, c)] @test Δsupp[r, c] == (abs(pt.p) > 3. * pt.std ? pt.p : 0.) @@ -52,13 +52,13 @@ end @test all(Δsupp[Δtrue .== 0.] .== 0.) τArr = [0., 2., Inf] - Δarr, Δpoint1, _ = supportEstimate(ANTSupport(), X, Y, τArr; estimSupport=eS) + Δarr, Δpoint1, _ = support_estimate(ANTSupport(), X, Y, τArr; estim_support=e_s) @test length(Δarr) == 3 @test all(Δarr[1] .!= 0.) - @test Δarr[2] ≈ supportEstimate(ANTSupport(), X, Y, 2.; estimSupport=eS)[1] + @test Δarr[2] ≈ support_estimate(ANTSupport(), X, Y, 2.; estim_support=e_s)[1] @test nnz(dropzeros(Δarr[3])) == 0 - Δarr, Δpoint, _ = supportEstimate(BootStdSupport(), X, Y, τArr; estimSupport=eS) + Δarr, Δpoint, _ = support_estimate(BootStdSupport(), X, Y, τArr; estim_support=e_s) @test all(Δarr[1] .!= 0.) @test nnz(dropzeros(Δarr[3])) == 0 @test all(r -> length(r.boot_p) == 1000, Δpoint) @@ -66,23 +66,23 @@ end @testset "DTrace" begin p, n = 5, 500 - X, Y, _ = generateData(p, n) - Xv, Yv, _ = generateData(p, n) + X, Y, _ = generate_data(p, n) + Xv, Yv, _ = generate_data(p, n) Sx, Sy = Symmetric(cov(X)), Symmetric(cov(Y)) Sxv, Syv = Symmetric(cov(Xv)), Symmetric(cov(Yv)) Λarr = [0.5, 0.2, 0.1, 0.05, 1e-8] - eΔarr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxv, Syv, Λarr) + eΔarr, i2, i_inf, loss2arr, loss_infarr = support_estimate(DTraceValidationSupport(), Sx, Sy, Sxv, Syv, Λarr) @test length(eΔarr) == length(Λarr) @test i2 == argmin(loss2arr) - @test iInf == argmin(lossInfarr) + @test i_inf == argmin(loss_infarr) @test nnz(eΔarr[1]) <= nnz(eΔarr[end]) @test Matrix(eΔarr[end]) ≈ inv(Sx) - inv(Sy) atol=1e-4 R = Sxv * Matrix(eΔarr[3]) * Syv - Syv + Sxv @test loss2arr[3] ≈ norm(R) - @test lossInfarr[3] ≈ maximum(abs, R) + @test loss_infarr[3] ≈ maximum(abs, R) end end diff --git a/test/ustats.jl b/test/ustats.jl index 3056a03..57de365 100644 --- a/test/ustats.jl +++ b/test/ustats.jl @@ -3,22 +3,22 @@ module UStatsTest using Test using Random, LinearAlgebra, Statistics using DiffPrecTest: compute, variance, OneSampleUstatistics, TwoSampleUstatistics, - createVarianceStatistic, createCovarianceStatistic, createSecondMomentStatistic + create_variance_statistic, create_covariance_statistic, create_second_moment_statistic Random.seed!(2) @testset "one sample" begin x = randn(10) - @test compute(createVarianceStatistic(), x) ≈ var(x) + @test compute(create_variance_statistic(), x) ≈ var(x) - s = createSecondMomentStatistic() + s = create_second_moment_statistic() @test compute(s, x) ≈ sum(abs2, x) / 10 @test variance(s, x) ≈ var( x .* x) X = randn(10, 2) - @test compute(createCovarianceStatistic(1, 2), X) ≈ cov(X[:, 1], X[:, 2]) + @test compute(create_covariance_statistic(1, 2), X) ≈ cov(X[:, 1], X[:, 2]) - s = createSecondMomentStatistic(1, 2) + s = create_second_moment_statistic(1, 2) @test compute(s, X) ≈ mean(X[:, 1] .* X[:, 2]) @test variance(s, X) ≈ var(X[:, 1] .* X[:, 2]) diff --git a/test/util.jl b/test/util.jl index 5d523c6..a07c0fc 100644 --- a/test/util.jl +++ b/test/util.jl @@ -48,7 +48,7 @@ end for ci=1:25 for ri=1:25 - v = DiffPrecTest._getElemSKron(S, view(X, i, :), ri, ci) + v = DiffPrecTest._get_elem_skron(S, view(X, i, :), ri, ci) @test sKron[ri, ci] ≈ v end end @@ -58,7 +58,7 @@ end sKron = ( kron(S, S2) + kron(S2, S) ) / 2. for ci=1:30 for ri=1:30 - v = DiffPrecTest._getElemSKron(S, S2, ri, ci) + v = DiffPrecTest._get_elem_skron(S, S2, ri, ci) @test sKron[ri, ci] ≈ v end end @@ -75,8 +75,8 @@ end for ci=1:25 for ri=1:25 - v = DiffPrecTest._getElemKron(S, view(X, i, :), ri, ci) - v1 = DiffPrecTest._getElemKron(view(X, i, :), S, ri, ci) + v = DiffPrecTest._get_elem_kron(S, view(X, i, :), ri, ci) + v1 = DiffPrecTest._get_elem_kron(view(X, i, :), S, ri, ci) @test k[ri, ci] ≈ v @test k1[ri, ci] ≈ v1 end @@ -87,24 +87,42 @@ end A = randn(3, 4) B = randn(2, 5) k = kron(A, B) - @test [DiffPrecTest._getElemKron(A, B, ri, ci) for ri=1:6, ci=1:20] ≈ k + @test [DiffPrecTest._get_elem_kron(A, B, ri, ci) for ri=1:6, ci=1:20] ≈ k +end + +@testset "covariance difference variance" begin + Sx = [2.0 0.4 -0.2; 0.4 1.5 0.3; -0.2 0.3 1.0] + Sy = [1.2 -0.1 0.2; -0.1 1.8 0.4; 0.2 0.4 1.5] + nx, ny = 13, 29 + p = size(Sx, 1) + + # Gaussian sample-covariance identity: (I + K)(S ⊗ S)/(n - 1), + # where K maps vec(A) to vec(transpose(A)). Independent samples add. + independent = kron(Sx, Sx) / (nx - 1) + kron(Sy, Sy) / (ny - 1) + transpose_order = vec(permutedims(reshape(1:p^2, p, p))) + reference = independent + independent[transpose_order, :] + for indices in (collect(1:p^2), [9, 2, 5, 4, 1, 2], Int[]) + actual = DiffPrecTest._var_sx_sy(Sx, nx, Sy, ny, indices) + @test size(actual) == (length(indices), length(indices)) + @test actual ≈ reference[indices, indices] + end end @testset "normal quantile" begin - @test DiffPrecTest._normQuantile(0.5) == 0. - @test DiffPrecTest._normQuantile(0.975) ≈ 1.9599639845400576 - @test DiffPrecTest._normQuantile(0.025) ≈ -1.9599639845400592 - @test DiffPrecTest._normQuantile(1e-10) ≈ -6.3613409024040575 - @test DiffPrecTest._penaltyLevel(10) ≈ 1.01 * DiffPrecTest._normQuantile(1. - 0.1 / 110) + @test DiffPrecTest._norm_quantile(0.5) == 0. + @test DiffPrecTest._norm_quantile(0.975) ≈ 1.9599639845400576 + @test DiffPrecTest._norm_quantile(0.025) ≈ -1.9599639845400592 + @test DiffPrecTest._norm_quantile(1e-10) ≈ -6.3613409024040575 + @test DiffPrecTest._penalty_level(10) ≈ 1.01 * DiffPrecTest._norm_quantile(1. - 0.1 / 110) end -@testset "getSupport" begin +@testset "get_support" begin p = 4 x = SymmetricSparseIterate(p) x[2, 1] = 0.5 x[4, 3] = -0.5 # negative entries are in the support x[3, 3] = -1e-4 # below the threshold - S = DiffPrecTest.getSupport(x) + S = DiffPrecTest.get_support(x) T = falses(p, p) T[2, 1] = T[1, 2] = T[4, 3] = T[3, 4] = true @test S == T @@ -113,7 +131,7 @@ end θ = SparseIterate(div(p * (p + 1), 2)) θ[tril_index(p, 2, 1)] = 0.5 θ[tril_index(p, 4, 3)] = -0.5 - @test DiffPrecTest.getSupport(θ, p) == T + @test DiffPrecTest.get_support(θ, p) == T # p × p matrix stored by columns, not symmetrized θ = SparseIterate(p * p) @@ -121,11 +139,11 @@ end θ[(3 - 1) * p + 4] = -0.5 T = falses(p, p) T[2, 1] = T[4, 3] = true - @test DiffPrecTest.getSupport(θ, p) == T + @test DiffPrecTest.get_support(θ, p) == T - @test_throws DimensionMismatch DiffPrecTest.getSupport(SparseIterate(5), p) + @test_throws DimensionMismatch DiffPrecTest.get_support(SparseIterate(5), p) - @test DiffPrecTest.getReducedSupport(θ) == (abs.(Vector(θ)) .> 1e-3) + @test DiffPrecTest.get_reduced_support(θ) == (abs.(Vector(θ)) .> 1e-3) end @testset "linear support" begin @@ -133,35 +151,35 @@ end S = falses(p, p) S[3, 1] = S[1, 3] = true S[2, 2] = true - supp = DiffPrecTest.getLinearSupport(2, 1, S) + supp = DiffPrecTest.get_linear_support(2, 1, S) @test supp[1:2] == [2, 5] # (2, 1) and (1, 2) first @test sort(supp[3:end]) == [3, 6, 9] - supp = DiffPrecTest.getLinearSupport(2, 2, S, falses(p, p)) + supp = DiffPrecTest.get_linear_support(2, 2, S, falses(p, p)) @test supp[1] == 6 @test sort(supp) == [3, 6, 9] # (row, col) in the support but not (col, row) - @test DiffPrecTest.makeFirst!([7, 2], p, 2, 1) == [2, 5, 7] - @test DiffPrecTest.makeFirst!([5, 7], p, 2, 1) == [2, 5, 7] - @test DiffPrecTest.makeFirst!([5, 7, 2], p, 2, 1) == [2, 5, 7] + @test DiffPrecTest.make_first!([7, 2], p, 2, 1) == [2, 5, 7] + @test DiffPrecTest.make_first!([5, 7], p, 2, 1) == [2, 5, 7] + @test DiffPrecTest.make_first!([5, 7, 2], p, 2, 1) == [2, 5, 7] # support given by indices into svec(Δ) ind = [tril_index(p, 3, 1), tril_index(p, 2, 2)] - s1 = DiffPrecTest._linearSupport(2, 1, ind, p) - s2 = DiffPrecTest.getLinearSupport(2, 1, S) + s1 = DiffPrecTest._linear_support(2, 1, ind, p) + s2 = DiffPrecTest.get_linear_support(2, 1, S) @test s1[1:2] == s2[1:2] @test sort(s1) == sort(s2) - @test DiffPrecTest._linearSupport(2, 1, S, p) == DiffPrecTest.getLinearSupport(2, 1, S) + @test DiffPrecTest._linear_support(2, 1, S, p) == DiffPrecTest.get_linear_support(2, 1, S) end -@testset "computeSimulationResult" begin - z = DiffPrecTest._normQuantile(0.975) +@testset "compute_simulation_result" begin + z = DiffPrecTest._norm_quantile(0.975) res = [DiffPrecResultNormal(0.1, 0.05), DiffPrecResultNormal(0.3, 0.05)] - sr = computeSimulationResult(res, 0.1) + sr = compute_simulation_result(res, 0.1) @test sr.bias ≈ 0.1 @test sr.coverage == 0.5 - @test sr.lenCoverage ≈ 2 * z * 0.05 + @test sr.len_coverage ≈ 2 * z * 0.05 @test sr.p_mean ≈ 0.2 @test sr.p_var ≈ 0.02 @@ -170,9 +188,9 @@ end shuffle!(b1) res = [DiffPrecResultBoot(0.1, b1), DiffPrecResultBoot(0.4, b2)] b1copy = copy(b1) - sr = computeSimulationResult(res, 0.1) + sr = compute_simulation_result(res, 0.1) @test sr.coverage == 0.5 - @test sr.lenCoverage ≈ 0.2 * 0.95 + @test sr.len_coverage ≈ 0.2 * 0.95 @test sr.bias ≈ 0.15 @test b1 == b1copy # replicates are not reordered end diff --git a/test/variance.jl b/test/variance.jl index 2f6e93d..fba6177 100644 --- a/test/variance.jl +++ b/test/variance.jl @@ -9,7 +9,7 @@ Random.seed!(134) # variance of ω'(H Δ - (Sy - Sx)) computed with dense matrices, where hess(A, B) # is the Hessian H for the covariances A of X and B of Y -function varianceAlt( +function variance_alt( hess, Sx, Sy, @@ -49,10 +49,10 @@ function varianceAlt( return σ1/nx + σ2/ny end -symHessian(A, B) = (kron(A, B) + kron(B, A)) / 2. -asymHessian(A, B) = kron(B, A) +sym_hessian(A, B) = (kron(A, B) + kron(B, A)) / 2. +asym_hessian(A, B) = kron(B, A) -function sampleData(n) +function sample_data(n) Ωx = [1. 0.4; 0.4 1] Ωy = [1. 0.3; 0.3 1] X = randn(n, 2) / cholesky(Ωx).U' @@ -62,45 +62,45 @@ end @testset "sym variance" begin n = 100 - X, Y = sampleData(n) + X, Y = sample_data(n) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) ind = [1, 2, 3, 4] - H = symHessian(Sx, Sy) + H = sym_hessian(Sx, Sy) tmp = zeros(length(ind)) tmp[1] = 1. ω = H \ tmp Δ = H \ (Sy[ind] - Sx[ind]) - @test variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind, Δ, ind) ≈ varianceAlt(symHessian, Sx, Sy, X, Y, ω, ind, Δ, ind) + @test variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind, Δ, ind) ≈ variance_alt(sym_hessian, Sx, Sy, X, Y, ω, ind, Δ, ind) # different supports for ω and Δ ω = randn(3) Δ = randn(2) - @test variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) ≈ varianceAlt(symHessian, Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) + @test variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) ≈ variance_alt(sym_hessian, Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) end @testset "asym variance" begin n = 100 - X, Y = sampleData(n) + X, Y = sample_data(n) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) ind = [1, 2, 3, 4] - H = asymHessian(Sx, Sy) + H = asym_hessian(Sx, Sy) tmp = zeros(length(ind)) tmp[1] = 1. ω = H \ tmp Δ = H \ (Sy[ind] - Sx[ind]) - @test variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind, Δ, ind) ≈ varianceAlt(asymHessian, Sx, Sy, X, Y, ω, ind, Δ, ind) + @test variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind, Δ, ind) ≈ variance_alt(asym_hessian, Sx, Sy, X, Y, ω, ind, Δ, ind) ω = randn(3) Δ = randn(2) - @test variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) ≈ varianceAlt(asymHessian, Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) + @test variance(AsymmetricOracleNormal(), Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) ≈ variance_alt(asym_hessian, Sx, Sy, X, Y, ω, [1, 2, 4], Δ, [2, 3]) end end