From 744acfbdd92655a88b38b180b84a1410759ba5e2 Mon Sep 17 00:00:00 2001 From: Mladen Kolar Date: Thu, 22 Oct 2026 15:48:23 +0400 Subject: [PATCH] Update the simulation scripts in exper/ for DiffPrecTest 0.2 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Environment and cluster jobs - exper/Project.toml with DiffPrecTest from the repository and CovSel from GitHub, and exper/README.md describing the simulations - JLD2 replaces JLD, Plots replaces PyPlot, and each script loads only the packages it uses; data are generated as matrices - The sbatch files take JULIA and SCRATCH from the environment instead of Midway2 paths, run Julia with --project, precompile once before the workers start, and use one BLAS thread per worker; the start scripts pass their arguments to sbatch (for the partition) and give the estimated and oracle runs of sims 8-11 different job names, so that their GNU parallel job logs do not collide - collect.sh loops over the coverage simulations and collect_results.jl skips missing repetitions Fixes - The oracle scripts of sims 8-12 passed supports as Bool matrices, which bootstrap no longer accepts; they now pass svec indices, one copy per entry - sim8.jl and sim10.jl were switched to oracle supports when the _o copies were made; they use estimated supports again, like sim9.jl and sim11.jl, and sim8.jl uses the dimensions its start script ran - sim8_YinXia.jl looped over four dimensions with three defined, and sim12_YinXia_power.jl skipped α = 1 - The coverage tables threw an InexactError when coverage * 1000 was not an integer; they are written with MultiRow cells, which LaTeXTabulars 1.0 no longer accepts as raw strings - sim5/summarize.jl rebuilt a p = 20 model that matched neither sim5a nor sim5b; it now reads the results of either and plots ROC curves, with the model shared in sim5/model.jl - The summaries of sims 1-4, 6 and 7 share summary.jl - Remove sim1/test.jl, sim2/sim2_test.jl and sim3/test.jl, scratch scripts that called functions that no longer exist Co-Authored-By: Claude Opus 5.5 --- .gitignore | 6 + README.md | 4 +- exper/Project.toml | 33 +++++ exper/README.md | 67 +++++++++ exper/cleanup.sh | 24 +--- exper/collect.sh | 120 +++++----------- exper/collect_results.jl | 35 ++--- exper/qq_plot.jl | 25 ++-- exper/sim1/plot_results.jl | 78 ++++------ exper/sim1/sbatch_sim1 | 31 ++-- exper/sim1/sim1.jl | 19 +-- exper/sim1/start_sbatch.sh | 9 +- exper/sim1/summarize.jl | 84 ++--------- exper/sim1/test.jl | 64 --------- exper/sim10/sbatch_sim10 | 34 +++-- exper/sim10/sbatch_sim10_power | 33 +++-- exper/sim10/sim10.jl | 26 ++-- exper/sim10/sim10_YinXia.jl | 24 ++-- exper/sim10/sim10_YinXia_power.jl | 32 +++-- exper/sim10/sim10_o.jl | 26 ++-- exper/sim10/sim10_o_power.jl | 32 +++-- exper/sim10/start_sbatch.sh | 13 +- exper/sim10/start_sbatch_power.sh | 18 ++- exper/sim11/sbatch_sim11 | 34 +++-- exper/sim11/sbatch_sim11_power | 36 ++--- exper/sim11/sim11.jl | 20 +-- exper/sim11/sim11_YinXia.jl | 24 ++-- exper/sim11/sim11_YinXia_power.jl | 32 +++-- exper/sim11/sim11_o.jl | 26 ++-- exper/sim11/sim11_o_power.jl | 32 +++-- exper/sim11/start_sbatch.sh | 13 +- exper/sim11/start_sbatch_power.sh | 18 ++- exper/sim12/sbatch_sim12_power | 29 ++-- exper/sim12/sim12_YinXia_power.jl | 31 ++-- exper/sim12/sim12_o_power.jl | 34 ++--- exper/sim12/start_sbatch_power.sh | 18 ++- exper/sim2/sbatch_sim2 | 31 ++-- exper/sim2/sim2.jl | 20 +-- exper/sim2/sim2_test.jl | 73 ---------- exper/sim2/start_sbatch.sh | 9 +- exper/sim2/summarize.jl | 72 ++-------- exper/sim3/sbatch_sim3 | 29 ++-- exper/sim3/sim3.jl | 20 +-- exper/sim3/start_sbatch.sh | 9 +- exper/sim3/summarize.jl | 154 +++++++------------- exper/sim3/test.jl | 151 -------------------- exper/sim4/sbatch_sim4 | 31 ++-- exper/sim4/sim4.jl | 20 +-- exper/sim4/start_sbatch.sh | 9 +- exper/sim4/summarize.jl | 150 +++++++------------- exper/sim5/model.jl | 18 +++ exper/sim5/sbatch_sim5 | 31 ++-- exper/sim5/sbatch_sim5b | 31 ++-- exper/sim5/sim5a.jl | 64 +++------ exper/sim5/sim5b.jl | 64 +++------ exper/sim5/summarize.jl | 227 +++++++++--------------------- exper/sim6/sbatch_sim6 | 30 ++-- exper/sim6/sim6.jl | 20 +-- exper/sim6/start_sbatch.sh | 9 +- exper/sim6/summarize.jl | 145 ++++++------------- exper/sim7/sbatch_sim7 | 31 ++-- exper/sim7/sim7.jl | 20 +-- exper/sim7/start_sbatch.sh | 9 +- exper/sim7/summarize.jl | 97 +++++++------ exper/sim8/aggregate_results.jl | 32 ++--- exper/sim8/collect_our.jl | 89 +++++------- exper/sim8/compute_coverage.jl | 28 ++-- exper/sim8/sbatch_sim8 | 34 +++-- exper/sim8/sbatch_sim8_power | 32 +++-- exper/sim8/sim8.jl | 28 ++-- exper/sim8/sim8_YinXia.jl | 28 ++-- exper/sim8/sim8_YinXia_power.jl | 32 +++-- exper/sim8/sim8_o.jl | 26 ++-- exper/sim8/sim8_o_power.jl | 32 +++-- exper/sim8/start_sbatch.sh | 13 +- exper/sim8/start_sbatch_power.sh | 18 ++- exper/sim9/sbatch_sim9 | 35 +++-- exper/sim9/sbatch_sim9_power | 33 +++-- exper/sim9/sim9.jl | 20 +-- exper/sim9/sim9_YinXia.jl | 24 ++-- exper/sim9/sim9_YinXia_power.jl | 32 +++-- exper/sim9/sim9_o.jl | 26 ++-- exper/sim9/sim9_o_power.jl | 32 +++-- exper/sim9/start_sbatch.sh | 13 +- exper/sim9/start_sbatch_power.sh | 18 ++- exper/start_all.sh | 21 +-- exper/summary.jl | 46 ++++++ 87 files changed, 1523 insertions(+), 1927 deletions(-) create mode 100644 exper/Project.toml create mode 100644 exper/README.md mode change 100644 => 100755 exper/cleanup.sh mode change 100644 => 100755 exper/collect.sh mode change 100644 => 100755 exper/sim1/start_sbatch.sh delete mode 100644 exper/sim1/test.jl mode change 100644 => 100755 exper/sim10/start_sbatch.sh mode change 100644 => 100755 exper/sim10/start_sbatch_power.sh mode change 100644 => 100755 exper/sim11/start_sbatch.sh mode change 100644 => 100755 exper/sim11/start_sbatch_power.sh mode change 100644 => 100755 exper/sim12/start_sbatch_power.sh delete mode 100644 exper/sim2/sim2_test.jl mode change 100644 => 100755 exper/sim2/start_sbatch.sh mode change 100644 => 100755 exper/sim3/start_sbatch.sh delete mode 100644 exper/sim3/test.jl mode change 100644 => 100755 exper/sim4/start_sbatch.sh create mode 100644 exper/sim5/model.jl mode change 100644 => 100755 exper/sim6/start_sbatch.sh mode change 100644 => 100755 exper/sim7/start_sbatch.sh mode change 100644 => 100755 exper/sim8/start_sbatch.sh mode change 100644 => 100755 exper/sim8/start_sbatch_power.sh mode change 100644 => 100755 exper/sim9/start_sbatch.sh mode change 100644 => 100755 exper/sim9/start_sbatch_power.sh mode change 100644 => 100755 exper/start_all.sh create mode 100644 exper/summary.jl diff --git a/.gitignore b/.gitignore index 18e98c3..d61aaa6 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,9 @@ .vscode/ *~ *jld +*.jld2 +/exper/Manifest*.toml +/exper/sim*/results/ +/exper/sim*/*.csv +/exper/sim*/*.tex +/exper/sim*/*.pdf diff --git a/README.md b/README.md index d9fb356..f9c2c20 100644 --- a/README.md +++ b/README.md @@ -44,8 +44,8 @@ CI = simulCIstudentized(straps, 0.95) `Distributions` also exports `estimate`; when both packages are loaded, call `DiffPrecTest.estimate`. -The `exper/` folder holds simulation scripts written for Julia 1.5 and version -0.1 of the package; they are kept as a record. +The `exper/` folder holds the simulation scripts, with their own environment; +see [`exper/README.md`](exper/README.md). Package depends on [ProximalBase.jl](https://github.com/mlakolar/ProximalBase.jl), [CoordinateDescent.jl](https://github.com/mlakolar/CoordinateDescent.jl) and diff --git a/exper/Project.toml b/exper/Project.toml new file mode 100644 index 0000000..812c4d1 --- /dev/null +++ b/exper/Project.toml @@ -0,0 +1,33 @@ +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +CovSel = "f50c8ead-2a5d-5d7a-afef-a6d87c1f6489" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DiffPrecTest = "bb4bfe13-bfd7-5c6f-a2d2-8072376c6806" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +LaTeXTabulars = "266f59ce-6e72-579c-98bb-27b39b5c037e" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +ProximalBase = "0a755509-c5b1-5d6f-9b24-59bcda3c37a3" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" + +[sources] +CovSel = {url = "https://github.com/mlakolar/CovSel.jl", rev = "master"} +DiffPrecTest = {path = ".."} + +[compat] +CSV = "1" +CovSel = "0.2" +DataFrames = "1" +DiffPrecTest = "0.2" +Distributions = "0.25" +JLD2 = "0.6" +LaTeXStrings = "1" +LaTeXTabulars = "1" +Plots = "1" +ProximalBase = "0.3, 0.4" +StatsBase = "0.34" +julia = "1.10" diff --git a/exper/README.md b/exper/README.md new file mode 100644 index 0000000..4d93867 --- /dev/null +++ b/exper/README.md @@ -0,0 +1,67 @@ +# Simulations + +Simulation scripts for DiffPrecTest. Each `simN` folder holds the scripts of one +simulation; one run of a script computes one repetition, and the Slurm jobs run +1000 repetitions (100 for sim5). + +| Simulation | Question | Model | Scripts | +|---|---|---|---| +| sim1, sim2 | coverage of the confidence interval for one entry, `Δ = 0` | `Ω` banded (sim1), `Ω[i, j] = 0.9^\|i-j\|` (sim2) | `simN.jl`, `summarize.jl`, `plot_results.jl` (sim1) | +| sim3, sim4 | coverage for one entry, sparse `Δ` | `Ωx` banded (sim3), `Ωx[i, j] = 0.7^\|i-j\|` (sim4) | `simN.jl`, `summarize.jl` | +| sim6, sim7 | coverage for one entry, banded `Δ` | banded `Ωx` and `Ωy` | `simN.jl`, `summarize.jl` | +| sim5 | support recovery | block diagonal `Ωx`, diagonal `Ωy` (`model.jl`), `p = 60` (a) and `100` (b) | `sim5a.jl`, `sim5b.jl`, `summarize.jl` | +| sim8–sim11 | simultaneous bootstrap intervals for all entries, `Δ = 0` | banded, hub, random sparse and block `Ωx` | `simN.jl` (estimated supports), `simN_o.jl` (oracle supports), `simN_YinXia.jl` | +| sim8–sim12 | power of the bootstrap test along `Ωy = (1 - α) Ωx + α Ωy₁` | as above; banded `Ωx`, `Ωy` (sim12) | `simN_o_power.jl`, `simN_YinXia_power.jl` | + +The coverage simulations compare nine methods: the three stage procedure +(`Sym-N`, `Asym-N`), Xia, Cai and Cai (2015) (`YinXia`), bootstrap standard +errors on the estimated support (`Sym-B`, `Asym-B`), and the same procedures on +the true support (`O-...`). + +## Setup + +With Julia 1.11 or later: + +```sh +julia --project=. -e 'using Pkg; Pkg.instantiate()' +``` + +`Project.toml` takes DiffPrecTest from this repository and CovSel from GitHub. +Julia 1.10 ignores these sources; there, run +`using Pkg; Pkg.develop(path=".."); Pkg.add(url="https://github.com/mlakolar/CovSel.jl")` +once instead. + +## Running + +One repetition, from a simulation folder: + +```sh +cd sim1 +julia --project=.. sim1.jl 1 1 1 /tmp/sim1 # rep ip iElem dir +``` + +All repetitions on a Slurm cluster, with GNU parallel running one repetition per +task (`./start_all.sh` submits sims 1-7): + +```sh +export SCRATCH=/path/to/scratch JULIA=/path/to/julia +cd sim1 +./start_sbatch.sh -p # arguments go to sbatch +``` + +Results go to `$SCRATCH/diffTest/`. The `simN_YinXia*.jl` scripts +run all repetitions in one process and write to the current folder; set +`NUM_REP` for fewer. + +## Summaries + +```sh +./collect.sh # sims 1-4, 6, 7 into sim*/results +julia --project=. sim1/summarize.jl # sim1.csv, sim1_table.tex +julia --project=. sim1/plot_results.jl # Q-Q plots and histograms +julia --project=. sim5/summarize.jl $SCRATCH/diffTest/sim5a a +julia --project=. sim8/collect_our.jl # bootstrap coverage, sims 8-11 +``` + +The scripts were first run with Julia 1.5 in 2019–2020. Julia's default random +number generator changed in Julia 1.7, so the same seeds give different data. diff --git a/exper/cleanup.sh b/exper/cleanup.sh old mode 100644 new mode 100755 index 9650b3e..a92fcac --- a/exper/cleanup.sh +++ b/exper/cleanup.sh @@ -1,21 +1,7 @@ #!/bin/bash -rm *~ -rm sim1/*log -rm sim1/*out -rm sim1/*~ -rm sim2/*log -rm sim2/*out -rm sim2/*~ -rm sim3/*log -rm sim3/*out -rm sim3/*~ -rm sim4/*log -rm sim4/*out -rm sim4/*~ -# rm sim5/*log -# rm sim5/*out -# rm sim5/*~ -rm sim6/*log -rm sim6/*out -rm sim6/*~ +## removes the GNU parallel job logs and editor backups written next to the +## scripts + +cd "$(dirname "$0")" +rm -f *~ sim*/*~ sim*/__*.runtask.log diff --git a/exper/collect.sh b/exper/collect.sh old mode 100644 new mode 100755 index 9940f5d..5bddf5a --- a/exper/collect.sh +++ b/exper/collect.sh @@ -1,87 +1,33 @@ -#sim 6 - -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_1_1 sim6_res_1_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_1_2 sim6_res_1_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_1_3 sim6_res_1_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_1_4 sim6_res_1_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_1_5 sim6_res_1_5 - -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_2_1 sim6_res_2_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_2_2 sim6_res_2_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_2_3 sim6_res_2_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_2_4 sim6_res_2_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim6/res_2_5 sim6_res_2_5 - - - - - - -# sim1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_1_1 sim1_res_1_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_1_2 sim1_res_1_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_1_3 sim1_res_1_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_2_1 sim1_res_2_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_2_2 sim1_res_2_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim1/res_2_3 sim1_res_2_3 - - -# sim2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_1_1 sim2_res_1_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_1_2 sim2_res_1_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_1_3 sim2_res_1_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_2_1 sim2_res_2_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_2_2 sim2_res_2_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim2/res_2_3 sim2_res_2_3 - -# sim3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_1_1 sim3_res_1_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_1_2 sim3_res_1_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_1_3 sim3_res_1_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_1_4 sim3_res_1_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_1_5 sim3_res_1_5 - -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_2_1 sim3_res_2_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_2_2 sim3_res_2_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_2_3 sim3_res_2_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_2_4 sim3_res_2_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_2_5 sim3_res_2_5 - -# # jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_3_1 sim3_res_3_1 -# # jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_3_2 sim3_res_3_2 -# # jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_3_3 sim3_res_3_3 -# # jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_3_4 sim3_res_3_4 -# # jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim3/res_3_5 sim3_res_3_5 - - -# sim4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_1_1 sim4_res_1_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_1_2 sim4_res_1_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_1_3 sim4_res_1_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_1_4 sim4_res_1_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_1_5 sim4_res_1_5 - -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_2_1 sim4_res_2_1 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_2_2 sim4_res_2_2 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_2_3 sim4_res_2_3 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_2_4 sim4_res_2_4 -# jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim4/res_2_5 sim4_res_2_5 - - - - - -# sim7 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_1_1 sim7_res_1_1 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_1_2 sim7_res_1_2 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_1_3 sim7_res_1_3 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_1_4 sim7_res_1_4 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_1_5 sim7_res_1_5 - -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_2_1 sim7_res_2_1 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_2_2 sim7_res_2_2 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_2_3 sim7_res_2_3 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_2_4 sim7_res_2_4 -jl1 collect_results.jl /scratch/midway2/mkolar/diffTest/sim7/res_2_5 sim7_res_2_5 - - +#!/bin/bash + +## collects the results of the coverage simulations from $SCRATCH/diffTest into +## sim*/results, where the summarize.jl scripts read them +## +## ./collect.sh [sim ...] +## +## collects sims 1-4, 6 and 7 unless given. JULIA (default: julia) and SCRATCH +## (default: $HOME/scratch) can be set in the environment. + +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} + +cd "$(dirname "$0")" + +sims=${@:-1 2 3 4 6 7} + +for sim in $sims +do + case $sim in + 1|2) numElem=3 ;; + *) numElem=5 ;; + esac + mkdir -p sim${sim}/results + for ip in 1 2 + do + for iElem in $(seq 1 $numElem) + do + echo "sim${sim} ${ip} ${iElem} ..." + $JULIA --project=. collect_results.jl $SCRATCH/diffTest/sim${sim}/res_${ip}_${iElem} sim${sim}/results/sim${sim}_res_${ip}_${iElem} + done + done +done diff --git a/exper/collect_results.jl b/exper/collect_results.jl index fde2621..a2b1d2d 100644 --- a/exper/collect_results.jl +++ b/exper/collect_results.jl @@ -1,26 +1,27 @@ +# Collect the estimates of one setting of a coverage simulation (sims 1-4, 6 and +# 7) into a matrix with one column per repetition: +# +# julia --project=. collect_results.jl resPath outPrefix +# +# reads resPath_rep.jld2 for rep = 1, ..., NUM_REP (default 1000) and writes +# outPrefix.jld2. Missing repetitions are skipped with a warning. using DiffPrecTest -using JLD +using JLD2 resPath = ARGS[1] outPrefix = ARGS[2] -# resPath = /scratch/mkolar/diffTest/sim1/res_2_2_$(rep).jld +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) -NUM_REP = 1000 -results = Array{Any}(undef, 9, NUM_REP) +estimates = Vector{Any}[] for rep=1:NUM_REP - global results - fname = "$(resPath)_$(rep).jld" - - # try - file = jldopen(fname, "r") - res = read(file, "est") - close(file) - - results[:, rep] = deepcopy(res) - # catch - # @show fname - # end + fname = "$(resPath)_$(rep).jld2" + if isfile(fname) + push!(estimates, load(fname, "est")) + else + @warn "missing $fname" + end end +results = reduce(hcat, estimates) -@save "$(outPrefix).jld" results +@save "$(outPrefix).jld2" results diff --git a/exper/qq_plot.jl b/exper/qq_plot.jl index c6a1475..eeec5ce 100644 --- a/exper/qq_plot.jl +++ b/exper/qq_plot.jl @@ -1,23 +1,22 @@ -using Distributions, StatsBase -using PyPlot +using Distributions +using Plots - -function qqplot(z; color="black") +# normal Q-Q plot of the sample z added to plt +function qqplot!(plt, z; color=:black) n = length(z) - grid = [(1 / (n + 1)):(1 / (n + 1)):(1.0 - (1 / (n + 1)));] + grid = (1:n) ./ (n + 1) qz = quantile(z, grid) - qd = quantile.(Ref(Distributions.Normal()), grid) + qd = quantile.(Normal(), grid) lims = 3.290 - x = range(-lims, stop=lims) - plot(x, x, color="grey", linestyle=":", linewidth=.25) - scatter(qz, qd, s=.75, color=color) + plot!(plt, [-lims, lims], [-lims, lims]; color=:grey, linestyle=:dot, linewidth=0.5, label="") + scatter!(plt, qz, qd; markersize=1, markerstrokewidth=0, color=color, label="") - xlim([-lims, lims]) - ylim([-lims, lims]) + xlims!(plt, -lims, lims) + ylims!(plt, -lims, lims) - nothing -end \ No newline at end of file + plt +end diff --git a/exper/sim1/plot_results.jl b/exper/sim1/plot_results.jl index 95e5e39..19bcfd7 100644 --- a/exper/sim1/plot_results.jl +++ b/exper/sim1/plot_results.jl @@ -1,71 +1,53 @@ -using JLD +# Normal Q-Q plots and histograms of the standardized estimates p / std of Sym-N +# (and of O-Sym-N in the Q-Q plots) in sim1, from the files +# results/sim1_res_ip_iElem.jld2 written by collect.sh. Writes qq_plot_sim1.pdf +# and hist_plot_sim1.pdf. +# +# julia --project=.. plot_results.jl using DiffPrecTest +using JLD2 using LaTeXStrings include("../qq_plot.jl") +cd(@__DIR__) + pArr = [100, 200] elemArr = [(5,5), (8, 7), (50, 25)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] -## qq-plots +standardized(res, j) = [res[j, i].p / res[j, i].std for i=1:size(res, 2)] +panelTitle(ip, iElem) = latexstring("p = $(pArr[ip]), \\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0") -fig_num = 0 +## qq-plots -fig = figure(figsize=(7.08661, 4.72441), dpi=1000) +panels = [] for ip=1:2 - for iElem=1:3 - global fig_num - - res = load("results/sim1_res_$(ip)_$(iElem).jld", "results") - fig_num += 1 + for iElem=1:3 + res = load("results/sim1_res_$(ip)_$(iElem).jld2", "results") - qq_my = [res[1, i].p / res[1, i].std for i=1:1000] - qq_oracle = [res[6, i].p / res[6, i].std for i=1:1000] - - ax = subplot(2, 3, fig_num) - qqplot(qq_oracle, color="grey") - qqplot(qq_my) - title( - latexstring("p = $(pArr[ip]), \\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"), - size="small" - ) - ax[:tick_params]("both", labelsize="xx-small", length=2, pad=2) + plt = plot(title=panelTitle(ip, iElem), titlefontsize=8, tickfontsize=5) + qqplot!(plt, standardized(res, 6), color=:grey) + qqplot!(plt, standardized(res, 1)) + push!(panels, plt) end end -tight_layout() -savefig("qq_plot_sim1.pdf") -close(fig) +savefig(plot(panels...; layout=(2, 3), size=(709, 472)), "qq_plot_sim1.pdf") ## hist-plots -fig_num = 0 - x = range(-3.290, stop=3.290, length=100) y = pdf.(Normal(), x) -fig = figure(figsize=(7.08661, 4.72441), dpi=1000) +panels = [] for ip=1:2 - for iElem=1:3 - global fig_num - - res = load("results/sim1_res_$(ip)_$(iElem).jld", "results") - fig_num += 1 - - ax = subplot(2, 3, fig_num) - plt[:hist]([res[1, i].p / res[1, i].std for i=1:1000], 100, density=true) - plot(x,y) - xlim(-3.290,3.290) - title( - latexstring("p = $(pArr[ip]), \\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"), - size="small" - ) - ax[:tick_params]("both", labelsize="xx-small", length=2, pad=2) + for iElem=1:3 + res = load("results/sim1_res_$(ip)_$(iElem).jld2", "results") + + plt = histogram(standardized(res, 1); bins=100, normalize=:pdf, label="", + title=panelTitle(ip, iElem), titlefontsize=8, tickfontsize=5) + plot!(plt, x, y; label="") + xlims!(plt, -3.290, 3.290) + push!(panels, plt) end end -tight_layout() -savefig("hist_plot_sim1.pdf") -close(fig) - - -## irina plots \ No newline at end of file +savefig(plot(panels...; layout=(2, 3), size=(709, 472)), "hist_plot_sim1.pdf") diff --git a/exper/sim1/sbatch_sim1 b/exper/sim1/sbatch_sim1 index d6c8466..353b8e7 100644 --- a/exper/sim1/sbatch_sim1 +++ b/exper/sim1/sbatch_sim1 @@ -1,20 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim1_${ip}_${iElem} sbatch_sim1 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim1, logs to +## $SCRATCH/output/sim1. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -##SBATCH --job-name=diff_sim1 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive -SCRATCH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim1 +LOG=$SCRATCH/output/sim1 -mkdir -p $SCRATCH/diffTest/sim1 -mkdir -p $SCRATCH/output/sim1 +mkdir -p $OUT +mkdir -p $LOG -module load parallel +command -v parallel > /dev/null || module load parallel -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH/output/sim1/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH/output/sim1/$SLURM_JOB_NAME.out.{1} $JULIA sim1.jl {1} ${1} ${2} $SCRATCH/diffTest/sim1 ::: {1..1000} +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' + +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim1.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim1/sim1.jl b/exper/sim1/sim1.jl index a88ddd4..000097f 100644 --- a/exper/sim1/sim1.jl +++ b/exper/sim1/sim1.jl @@ -1,9 +1,12 @@ +# Coverage of confidence intervals for one entry of Δ when Ωx = Ωy (Δ = 0). +# +# julia --project=.. sim1.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -43,8 +46,8 @@ dist_Y = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(134 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -65,4 +68,4 @@ ri, ci = elemArr[iElem] @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim1/start_sbatch.sh b/exper/sim1/start_sbatch.sh old mode 100644 new mode 100755 index 8c14b67..a68e660 --- a/exper/sim1/start_sbatch.sh +++ b/exper/sim1/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" + for iElem in 1 2 3 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim1_${ip}_${iElem} sbatch_sim1 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim1_${ip}_${iElem} sbatch_sim1 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim1/summarize.jl b/exper/sim1/summarize.jl index d5a85e0..42bb41b 100644 --- a/exper/sim1/summarize.jl +++ b/exper/sim1/summarize.jl @@ -1,83 +1,23 @@ -using JLD -using DiffPrecTest -using DataFrames, DataFramesMeta, CSV -using LaTeXStrings, LaTeXTabulars +# Coverage, length and bias of the confidence intervals in sim1 (Δ = 0), from +# the files results/sim1_res_ip_iElem.jld2 written by collect.sh. Writes +# sim1.csv and sim1_table.tex. +# +# julia --project=.. summarize.jl +using CSV -cd("exper/sim1") +include("../summary.jl") + +cd(@__DIR__) pArr = [100, 200, 500] elemArr = [(5,5), (8, 7), (50, 25)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) +df = newSummary() for ip=1:2 for iElem=1:3 - res = load("results/sim1_res_$(ip)_$(iElem).jld", "results") - - for j in 1:9 - global pArr - global elemArr - global methodArr - global df - - sr = computeSimulationResult([res[j, i] for i=1:1000], 0.) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end + summarizeSetting!(df, "results/sim1_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], 0.) end end CSV.write("sim1.csv", df) - - -### create a table - -latex_tabular("sim1_table.tex", - Tabular("lllccc"), - [ - Rule(:top), - ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"], - Rule(:mid), - vcat([ - hcat( - (a = fill("", 9, 1); - b = latexstring("p = $(pArr[ip])"); - a[1] = "\\multirow{9}{*}{$(b)}"; - a - ), - vcat([ - hcat( - (a = fill("", 3, 1); - b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"); - a[1] = "\\multirow{3}{*}{$(b)}"; - a - ), - Matrix{Any}( - @linq df |> - where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"), - :p .== pArr[ip], - :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |> select(:method, - coverage = convert.(Int, :coverage * 1000), - lenCI = round.(:lenCI, digits=3), - bias = round.(:bias * 1000, digits=1)) - ) - ) - for iElem=1:3 - ]... - ) - ) - for ip=1:2 - ]...), - Rule(:bottom) - ] -) - - - - - - - - - - +coverageTable("sim1_table.tex", df, pArr[1:2], elemArr, zeros(2, 3)) diff --git a/exper/sim1/test.jl b/exper/sim1/test.jl deleted file mode 100644 index 83fac1b..0000000 --- a/exper/sim1/test.jl +++ /dev/null @@ -1,64 +0,0 @@ -using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD - -pArr = [100, 200, 500, 300] -elemArr = [(5,5), (8, 7), (50, 25)] -n = 100 -est = Array{Any}(undef, 5) # number of methods - - -# rep = parse(Int,ARGS[1]) -# ip = parse(Int,ARGS[2]) -# iElem = parse(Int,ARGS[3]) -# dir = ARGS[4] - -rep = 1 -ip = 2 -iElem = 1 - -Random.seed!(134) - -p = pArr[ip] -# generate model -Ω = Matrix{Float64}(I, p, p) -for l=1:p-1 - Ω[l , l+1] = 0.6 - Ω[l+1, l ] = 0.6 -end -for l=1:p-2 - Ω[l , l+2] = 0.3 - Ω[l+2, l ] = 0.3 -end -d = Vector{Float64}(undef, p) -rand!(Uniform(0.5, 2.5), d) -d .= sqrt.(d) -D = Diagonal(d) -Σ = inv(Symmetric(D * Ω * D)) - -dist_X = MvNormal(convert(Matrix, Σ)) -dist_Y = MvNormal(convert(Matrix, Σ)) - -# generate data -Random.seed!(134 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' - -@show (ip, iElem) - -ri, ci = elemArr[iElem] -indE = (ci - 1) * p + ri -indOracle = [indE] - -@time est[1], _, _, indS = DiffPrecTest.estimate(SymmetricNormal(), X, Y, indE) -@show indS -@time est[2] = DiffPrecTest.estimate(SeparateNormal(), X, Y, indE) -@time est[3] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indS) -@time est[4] = DiffPrecTest.estimate(SymmetricOracleNormal(), Symmetric(cov(X)), n, Symmetric(cov(Y)), n, indOracle) -@time est[5] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indOracle) - - -# @save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est diff --git a/exper/sim10/sbatch_sim10 b/exper/sim10/sbatch_sim10 index ee5489e..0dd47b9 100644 --- a/exper/sim10/sbatch_sim10 +++ b/exper/sim10/sbatch_sim10 @@ -1,23 +1,31 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim10_${ip} sbatch_sim10 ip script +## with script sim10 (estimated supports) or sim10_o (oracle supports) +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/$SCRIPT, logs to +## $SCRATCH/output/$SCRIPT. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -module load parallel -module load hdf5 - -SCRATCH_PATH=/scratch/midway2/mkolar -TMPOUT=$SCRATCH_PATH/output/sim10_o -OUT=$SCRATCH_PATH/diffTest/sim10_o -JULIA_PATH=/home/mkolar/project-midway/julia-1.5.0/bin/julia -FNAME=sim10_o.jl +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +SCRIPT=${2:-sim10} +OUT=$SCRATCH/diffTest/$SCRIPT +LOG=$SCRATCH/output/$SCRIPT mkdir -p $OUT -mkdir -p $TMPOUT +mkdir -p $LOG + +command -v parallel > /dev/null || module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$TMPOUT/$SLURM_JOB_NAME.err.{1} --output=$TMPOUT/$SLURM_JOB_NAME.out.{1} $JULIA_PATH $FNAME {1} ${1} $OUT ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT $SCRIPT.jl {1} ${1} $OUT ::: {1..1000} diff --git a/exper/sim10/sbatch_sim10_power b/exper/sim10/sbatch_sim10_power index 117d0ce..d1298ce 100644 --- a/exper/sim10/sbatch_sim10_power +++ b/exper/sim10/sbatch_sim10_power @@ -1,24 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch_power.sh as +## sbatch --job-name=diff_sim10_${ip}_${ialpha}_power sbatch_sim10_power ip ialpha +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim10_o_power_${1}_${2}, logs to +## $SCRATCH/output/sim10_o_power_${1}_${2}. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -module load parallel -module load hdf5 - - -SCRATCH_PATH=/scratch/midway2/mkolar -TMPOUT=$SCRATCH_PATH/output/sim10_o_power -OUT=$SCRATCH_PATH/diffTest/sim10_o_power_${1}_${2} -JULIA_PATH=/home/mkolar/project-midway/julia-1.5.0/bin/julia -FNAME=sim10_o_power.jl +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim10_o_power_${1}_${2} +LOG=$SCRATCH/output/sim10_o_power_${1}_${2} mkdir -p $OUT -mkdir -p $TMPOUT +mkdir -p $LOG + +command -v parallel > /dev/null || module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$TMPOUT/$SLURM_JOB_NAME.err.{1} --output=$TMPOUT/$SLURM_JOB_NAME.out.{1} $JULIA_PATH $FNAME {1} ${1} ${2} $OUT ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim10_o_power.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim10/sim10.jl b/exper/sim10/sim10.jl index 7f5bdb0..864b453 100644 --- a/exper/sim10/sim10.jl +++ b/exper/sim10/sim10.jl @@ -1,9 +1,13 @@ +# Simultaneous confidence intervals for all entries of Δ by the bootstrap, with +# estimated supports, when Ωx = Ωy (random sparse Ωx). +# +# julia --project=.. sim10.jl rep ip dir +# +# saves the bootstrap estimates and supports to dir/res_ip_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -42,18 +46,14 @@ dist_X = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(12346 + rep) -X = rand(dist_X, n)' -Y = rand(dist_X, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = falses(p, p) -end -@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init) +@time boot_res, eS = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld" boot_res +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS diff --git a/exper/sim10/sim10_YinXia.jl b/exper/sim10/sim10_YinXia.jl index 486109e..bda949b 100644 --- a/exper/sim10/sim10_YinXia.jl +++ b/exper/sim10/sim10_YinXia.jl @@ -1,18 +1,22 @@ +# The estimator of Xia, Cai and Cai (2015) for all entries of Δ when Ωx = Ωy +# (random sparse Ωx). Runs all repetitions: +# +# julia --project=.. sim10_YinXia.jl +# +# saves the estimates to YinXia_res_ip.jld2; set NUM_REP for fewer repetitions. using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) bnm = Binomial(1, 0.05) -for ip=1:3 +for ip=eachindex(pArr) p = pArr[ip] Random.seed!(54298) @@ -45,8 +49,8 @@ for ip=1:3 # generate data Random.seed!(12346 + rep) - X = rand(dist_X, n)' - Y = rand(dist_X, n)' + X = permutedims(rand(dist_X, n)) + Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -61,5 +65,5 @@ for ip=1:3 end - @save "YinXia_res_$(ip).jld" res + @save "YinXia_res_$(ip).jld2" res end diff --git a/exper/sim10/sim10_YinXia_power.jl b/exper/sim10/sim10_YinXia_power.jl index 2cbfcd5..dfbd1a1 100644 --- a/exper/sim10/sim10_YinXia_power.jl +++ b/exper/sim10/sim10_YinXia_power.jl @@ -1,24 +1,26 @@ +# The estimator of Xia, Cai and Cai (2015) for all entries of Δ along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (random sparse Ωx). Runs all repetitions: +# +# julia --project=.. sim10_YinXia_power.jl +# +# saves the estimates to YinXia_power_res_ip_ialpha.jld2; set NUM_REP for fewer +# repetitions. using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using StatsBase: sample +using Distributions +using JLD2 -# rep = parse(Int,ARGS[1]) -# ip = parse(Int,ARGS[2]) -# dir = ARGS[3] - alpha = 0.1:0.1:1 pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) bnm = Binomial(1, 0.05) -for ip=1:3 - for ialpha=1:10 +for ip=eachindex(pArr) + for ialpha=eachindex(alpha) p = pArr[ip] Random.seed!(54298) @@ -73,8 +75,8 @@ for ip=1:3 # generate data Random.seed!(12346 + rep*100 + ialpha) - X = rand(dist_X, n)' - Y = rand(dist_Y, n)' + X = permutedims(rand(dist_X, n)) + Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -89,6 +91,6 @@ for ip=1:3 end - @save "YinXia_power_res_$(ip)_$(ialpha).jld" res + @save "YinXia_power_res_$(ip)_$(ialpha).jld2" res end end diff --git a/exper/sim10/sim10_o.jl b/exper/sim10/sim10_o.jl index db58321..f48ded7 100644 --- a/exper/sim10/sim10_o.jl +++ b/exper/sim10/sim10_o.jl @@ -1,9 +1,13 @@ +# Simultaneous confidence intervals for all entries of Δ by the bootstrap, with +# oracle supports, when Ωx = Ωy (random sparse Ωx). +# +# julia --project=.. sim10_o.jl rep ip dir +# +# saves the bootstrap estimates to dir/res_ip_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -42,18 +46,16 @@ dist_X = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(12346 + rep) -X = rand(dist_X, n)' -Y = rand(dist_X, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = falses(p, p) -end +# 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) -@save "$(dir)/res_$(ip)_$(rep).jld" boot_res +@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 7807e5b..48f6731 100644 --- a/exper/sim10/sim10_o_power.jl +++ b/exper/sim10/sim10_o_power.jl @@ -1,9 +1,15 @@ +# Power of the bootstrap test of Ωx = Ωy, with oracle supports, along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (random sparse Ωx). +# +# julia --project=.. sim10_o_power.jl rep ip ialpha dir +# +# saves the bootstrap estimates to dir/res_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using ProximalBase: sub2indLowerTriangular +using StatsBase: sample +using Distributions +using JLD2 @show gethostname() @@ -66,19 +72,17 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(12346 + rep*100 + ialpha) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -S = BitArray( abs(Δ[r, c]) > 1e-4 for r = 1:p, c = 1:p ) - -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = S -end +# oracle supports: the nonzero entries of Δ as indices into svec(Δ); bootstrap +# reorders each support, so every entry gets its own copy +S = [sub2indLowerTriangular(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) -@save "$(dir)/res_$(rep).jld" boot_res +@save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim10/start_sbatch.sh b/exper/sim10/start_sbatch.sh old mode 100644 new mode 100755 index 09a9a8e..468efd5 --- a/exper/sim10/start_sbatch.sh +++ b/exper/sim10/start_sbatch.sh @@ -1,8 +1,15 @@ #!/bin/bash -for ip in "1" "2" "3" +## submits one job per dimension for estimated (sim10) and oracle (sim10_o) +## supports; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for script in sim10 sim10_o do - echo "sbatch ${ip} ..." - sbatch --job-name=diff_sim10_${ip} sbatch_sim10 ${ip} + for ip in 1 2 3 + do + echo "sbatch ${script} ${ip} ..." + sbatch "$@" --job-name=diff_${script}_${ip} sbatch_sim10 ${ip} ${script} echo "... done" + done done diff --git a/exper/sim10/start_sbatch_power.sh b/exper/sim10/start_sbatch_power.sh old mode 100644 new mode 100755 index 52bdbfe..d5ce6b0 --- a/exper/sim10/start_sbatch_power.sh +++ b/exper/sim10/start_sbatch_power.sh @@ -1,11 +1,15 @@ #!/bin/bash -for ip in "1" "2" "3" +## submits one job per dimension and signal strength; arguments are passed on to +## sbatch, e.g. +## ./start_sbatch_power.sh -p + +for ip in 1 2 3 do - for ialpha in {1..10} - do - echo "sbatch ${ip} ${ialpha} ..." - sbatch --job-name=diff_sim10_${ip}_${ialpha}_power sbatch_sim10_power ${ip} ${ialpha} - echo "... done" - done + for ialpha in {1..10} + do + echo "sbatch ${ip} ${ialpha} ..." + sbatch "$@" --job-name=diff_sim10_${ip}_${ialpha}_power sbatch_sim10_power ${ip} ${ialpha} + echo "... done" + done done diff --git a/exper/sim11/sbatch_sim11 b/exper/sim11/sbatch_sim11 index c0d0641..d63946f 100644 --- a/exper/sim11/sbatch_sim11 +++ b/exper/sim11/sbatch_sim11 @@ -1,23 +1,31 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim11_${ip} sbatch_sim11 ip script +## with script sim11 (estimated supports) or sim11_o (oracle supports) +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/$SCRIPT, logs to +## $SCRATCH/output/$SCRIPT. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=200 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -module load parallel -module load hdf5 - -SCRATCH_PATH=/scratch/midway2/mkolar -TMPOUT=$SCRATCH_PATH/output/sim11_o -OUT=$SCRATCH_PATH/diffTest/sim11_o -JULIA_PATH=/home/mkolar/project-midway/julia-1.1.1/bin/julia -FNAME=sim11_o.jl +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +SCRIPT=${2:-sim11} +OUT=$SCRATCH/diffTest/$SCRIPT +LOG=$SCRATCH/output/$SCRIPT mkdir -p $OUT -mkdir -p $TMPOUT +mkdir -p $LOG + +command -v parallel > /dev/null || module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$TMPOUT/$SLURM_JOB_NAME.err.{1} --output=$TMPOUT/$SLURM_JOB_NAME.out.{1} $JULIA_PATH $FNAME {1} ${1} $OUT ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT $SCRIPT.jl {1} ${1} $OUT ::: {1..1000} diff --git a/exper/sim11/sbatch_sim11_power b/exper/sim11/sbatch_sim11_power index 0c9413a..d758547 100644 --- a/exper/sim11/sbatch_sim11_power +++ b/exper/sim11/sbatch_sim11_power @@ -1,27 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir - -##SBATCH --ntasks=200 -##SBATCH --partition=broadwl +## submitted by start_sbatch_power.sh as +## sbatch --job-name=diff_sim11_${ip}_${ialpha}_power sbatch_sim11_power ip ialpha +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim11_o_power_${1}_${2}, logs to +## $SCRATCH/output/sim11_o_power_${1}_${2}. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=28 -#SBATCH --partition=standard #SBATCH --time=7-00:00:00 -# module load parallel -# module load hdf5 - -SCRATCH_PATH=/scratch/mkolar -TMPOUT=$SCRATCH_PATH/output/sim11_o_power -OUT=$SCRATCH_PATH/diffTest/sim11_o_power_${1}_${2} -#JULIA_PATH=/home/mkolar/project-midway/julia-1.1.1/bin/julia -JULIA_PATH=/project/mkolar_space/julia-1.1.1/bin/julia -FNAME=sim11_o_power.jl +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim11_o_power_${1}_${2} +LOG=$SCRATCH/output/sim11_o_power_${1}_${2} mkdir -p $OUT -mkdir -p $TMPOUT +mkdir -p $LOG + +command -v parallel > /dev/null || module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$TMPOUT/$SLURM_JOB_NAME.err.{1} --output=$TMPOUT/$SLURM_JOB_NAME.out.{1} $JULIA_PATH $FNAME {1} ${1} ${2} $OUT ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim11_o_power.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim11/sim11.jl b/exper/sim11/sim11.jl index 9624b84..043fb16 100644 --- a/exper/sim11/sim11.jl +++ b/exper/sim11/sim11.jl @@ -1,9 +1,13 @@ +# Simultaneous confidence intervals for all entries of Δ by the bootstrap, with +# estimated supports, when Ωx = Ωy (Ωx the inverse of a block diagonal Σx). +# +# julia --project=.. sim11.jl rep ip dir +# +# saves the bootstrap estimates and supports to dir/res_ip_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -44,8 +48,8 @@ dist_X = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(12346 + rep) -X = rand(dist_X, n)' -Y = rand(dist_X, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -54,4 +58,4 @@ Sy = Symmetric( Y'Y / n ) @time boot_res, eS = bootstrap(X, Y) -@save "$(dir)/res_$(ip)_$(rep).jld" boot_res eS +@save "$(dir)/res_$(ip)_$(rep).jld2" boot_res eS diff --git a/exper/sim11/sim11_YinXia.jl b/exper/sim11/sim11_YinXia.jl index 12d7ae7..5f5bd3a 100644 --- a/exper/sim11/sim11_YinXia.jl +++ b/exper/sim11/sim11_YinXia.jl @@ -1,16 +1,20 @@ +# The estimator of Xia, Cai and Cai (2015) for all entries of Δ when Ωx = Ωy +# (Ωx the inverse of a block diagonal Σx). Runs all repetitions: +# +# julia --project=.. sim11_YinXia.jl +# +# saves the estimates to YinXia_res_ip.jld2; set NUM_REP for fewer repetitions. using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) -for ip=1:3 +for ip=eachindex(pArr) p = pArr[ip] Random.seed!(54298) @@ -47,8 +51,8 @@ for ip=1:3 # generate data Random.seed!(12346 + rep) - X = rand(dist_X, n)' - Y = rand(dist_X, n)' + X = permutedims(rand(dist_X, n)) + Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -63,5 +67,5 @@ for ip=1:3 end - @save "YinXia_res_$(ip).jld" res + @save "YinXia_res_$(ip).jld2" res end diff --git a/exper/sim11/sim11_YinXia_power.jl b/exper/sim11/sim11_YinXia_power.jl index a2bd097..97a05b7 100644 --- a/exper/sim11/sim11_YinXia_power.jl +++ b/exper/sim11/sim11_YinXia_power.jl @@ -1,22 +1,24 @@ +# The estimator of Xia, Cai and Cai (2015) for all entries of Δ along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (Ωx the inverse of a block diagonal Σx). Runs all repetitions: +# +# julia --project=.. sim11_YinXia_power.jl +# +# saves the estimates to YinXia_power_res_ip_ialpha.jld2; set NUM_REP for fewer +# repetitions. using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using StatsBase: sample +using Distributions +using JLD2 -# rep = parse(Int,ARGS[1]) -# ip = parse(Int,ARGS[2]) -# dir = ARGS[3] - alpha = 0.1:0.1:1 pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) -for ip=1:3 - for ialpha=1:10 +for ip=eachindex(pArr) + for ialpha=eachindex(alpha) p = pArr[ip] Random.seed!(54298) @@ -77,8 +79,8 @@ for ip=1:3 # generate data Random.seed!(12346 + rep*100 + ialpha) - X = rand(dist_X, n)' - Y = rand(dist_Y, n)' + X = permutedims(rand(dist_X, n)) + Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -93,6 +95,6 @@ for ip=1:3 end - @save "YinXia_power_res_$(ip)_$(ialpha).jld" res + @save "YinXia_power_res_$(ip)_$(ialpha).jld2" res end end diff --git a/exper/sim11/sim11_o.jl b/exper/sim11/sim11_o.jl index 0988fbd..17e2091 100644 --- a/exper/sim11/sim11_o.jl +++ b/exper/sim11/sim11_o.jl @@ -1,9 +1,13 @@ +# Simultaneous confidence intervals for all entries of Δ by the bootstrap, with +# oracle supports, when Ωx = Ωy (Ωx the inverse of a block diagonal Σx). +# +# julia --project=.. sim11_o.jl rep ip dir +# +# saves the bootstrap estimates to dir/res_ip_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -44,18 +48,16 @@ dist_X = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(12346 + rep) -X = rand(dist_X, n)' -Y = rand(dist_X, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_X, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = falses(p, p) -end +# 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) -@save "$(dir)/res_$(ip)_$(rep).jld" boot_res +@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 70b36bf..b3a3cbb 100644 --- a/exper/sim11/sim11_o_power.jl +++ b/exper/sim11/sim11_o_power.jl @@ -1,9 +1,15 @@ +# Power of the bootstrap test of Ωx = Ωy, with oracle supports, along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (Ωx the inverse of a block diagonal Σx). +# +# julia --project=.. sim11_o_power.jl rep ip ialpha dir +# +# saves the bootstrap estimates to dir/res_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using ProximalBase: sub2indLowerTriangular +using StatsBase: sample +using Distributions +using JLD2 @show gethostname() @@ -72,19 +78,17 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(12346 + rep*100 + ialpha) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -S = BitArray( abs(Δ[r, c]) > 1e-4 for r = 1:p, c = 1:p ) - -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = S -end +# oracle supports: the nonzero entries of Δ as indices into svec(Δ); bootstrap +# reorders each support, so every entry gets its own copy +S = [sub2indLowerTriangular(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) -@save "$(dir)/res_$(rep).jld" boot_res +@save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim11/start_sbatch.sh b/exper/sim11/start_sbatch.sh old mode 100644 new mode 100755 index f1142be..82c8891 --- a/exper/sim11/start_sbatch.sh +++ b/exper/sim11/start_sbatch.sh @@ -1,8 +1,15 @@ #!/bin/bash -for ip in "1" "2" "3" +## submits one job per dimension for estimated (sim11) and oracle (sim11_o) +## supports; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for script in sim11 sim11_o do - echo "sbatch ${ip} ..." - sbatch --job-name=diff_sim11_${ip} sbatch_sim11 ${ip} + for ip in 1 2 3 + do + echo "sbatch ${script} ${ip} ..." + sbatch "$@" --job-name=diff_${script}_${ip} sbatch_sim11 ${ip} ${script} echo "... done" + done done diff --git a/exper/sim11/start_sbatch_power.sh b/exper/sim11/start_sbatch_power.sh old mode 100644 new mode 100755 index fb69cdc..b9fda2d --- a/exper/sim11/start_sbatch_power.sh +++ b/exper/sim11/start_sbatch_power.sh @@ -1,11 +1,15 @@ #!/bin/bash -for ip in "1" "2" "3" +## submits one job per dimension and signal strength; arguments are passed on to +## sbatch, e.g. +## ./start_sbatch_power.sh -p + +for ip in 1 2 3 do - for ialpha in {1..10} - do - echo "sbatch ${ip} ${ialpha} ..." - sbatch --job-name=diff_sim11_${ip}_${ialpha}_power sbatch_sim11_power ${ip} ${ialpha} - echo "... done" - done + for ialpha in {1..10} + do + echo "sbatch ${ip} ${ialpha} ..." + sbatch "$@" --job-name=diff_sim11_${ip}_${ialpha}_power sbatch_sim11_power ${ip} ${ialpha} + echo "... done" + done done diff --git a/exper/sim12/sbatch_sim12_power b/exper/sim12/sbatch_sim12_power index 0ff452c..23a9514 100644 --- a/exper/sim12/sbatch_sim12_power +++ b/exper/sim12/sbatch_sim12_power @@ -1,20 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch_power.sh as +## sbatch --job-name=diff_sim12_${ip}_${ialpha}_power sbatch_sim12_power ip ialpha +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim12_o_power_${1}_${2}, logs to +## $SCRATCH/output/sim12_o_power_${1}_${2}. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=28 -#SBATCH --partition=standard #SBATCH --time=7-00:00:00 -SCRATCH_PATH=/scratch/mkolar -TMPOUT=$SCRATCH_PATH/output/sim12_o_power_${1}_${2} -OUT=$SCRATCH_PATH/diffTest/sim12_o_power_${1}_${2} -JULIA_PATH=/project/mkolar_space/julia-1.1.1/bin/julia -FNAME=sim12_o_power.jl +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim12_o_power_${1}_${2} +LOG=$SCRATCH/output/sim12_o_power_${1}_${2} mkdir -p $OUT -mkdir -p $TMPOUT +mkdir -p $LOG +command -v parallel > /dev/null || module load parallel -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$TMPOUT/$SLURM_JOB_NAME.err.{1} --output=$TMPOUT/$SLURM_JOB_NAME.out.{1} $JULIA_PATH $FNAME {1} ${1} ${2} $OUT ::: {1..1000} +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' + +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim12_o_power.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim12/sim12_YinXia_power.jl b/exper/sim12/sim12_YinXia_power.jl index 4c6bc81..04d99ad 100644 --- a/exper/sim12/sim12_YinXia_power.jl +++ b/exper/sim12/sim12_YinXia_power.jl @@ -1,22 +1,23 @@ +# The estimator of Xia, Cai and Cai (2015) for all entries of Δ along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (banded Ωx and Ωy). Runs all repetitions: +# +# julia --project=.. sim12_YinXia_power.jl +# +# saves the estimates to YinXia_power_res_ip_ialpha.jld2; set NUM_REP for fewer +# repetitions. using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 -# rep = parse(Int,ARGS[1]) -# ip = parse(Int,ARGS[2]) -# dir = ARGS[3] - alpha = 0.:0.1:1. pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) -for ip=1:3 - for ialpha=1:10 +for ip=eachindex(pArr) + for ialpha=eachindex(alpha) p = pArr[ip] Random.seed!(7689) @@ -63,8 +64,8 @@ for ip=1:3 # generate data Random.seed!(12346 + rep*100 + ialpha) - X = rand(dist_X, n)' - Y = rand(dist_Y, n)' + X = permutedims(rand(dist_X, n)) + Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -79,6 +80,6 @@ for ip=1:3 end - @save "YinXia_power_res_$(ip)_$(ialpha).jld" res + @save "YinXia_power_res_$(ip)_$(ialpha).jld2" res end end diff --git a/exper/sim12/sim12_o_power.jl b/exper/sim12/sim12_o_power.jl index 4f64840..7bb1d2c 100644 --- a/exper/sim12/sim12_o_power.jl +++ b/exper/sim12/sim12_o_power.jl @@ -1,9 +1,14 @@ +# Power of the bootstrap test of Ωx = Ωy, with oracle supports, along the path +# Ωy = (1 - α) Ωx + α Ωy₁ (banded Ωx and Ωy). +# +# julia --project=.. sim12_o_power.jl rep ip ialpha dir +# +# saves the bootstrap estimates to dir/res_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using ProximalBase: sub2indLowerTriangular +using Distributions +using JLD2 @show gethostname() @@ -15,7 +20,6 @@ dir = ARGS[4] alpha = 0.:0.1:1. pArr = [50, 100, 150] n = 300 -NUM_REP = 1000 p = pArr[ip] @@ -55,23 +59,19 @@ D = Diagonal(d) dist_X = MvNormal(convert(Matrix, Σx)) dist_Y = MvNormal(convert(Matrix, Σy)) -rp = div(p * (p+1), 2) -res = Matrix{DiffPrecResultNormal}(undef, rp, NUM_REP) - # generate data Random.seed!(12346 + rep*100 + ialpha) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) -S = BitArray( abs(Δ[r, c]) > 1e-4 for r = 1:p, c = 1:p ) -eS_init = Array{BitArray}(undef, div((p + 1)*p, 2)) -for j=1:div((p + 1)*p, 2) - eS_init[j] = S -end +# oracle supports: the nonzero entries of Δ as indices into svec(Δ); bootstrap +# reorders each support, so every entry gets its own copy +S = [sub2indLowerTriangular(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) -@save "$(dir)/res_$(rep).jld" boot_res +@save "$(dir)/res_$(rep).jld2" boot_res diff --git a/exper/sim12/start_sbatch_power.sh b/exper/sim12/start_sbatch_power.sh old mode 100644 new mode 100755 index 44e315b..eef2873 --- a/exper/sim12/start_sbatch_power.sh +++ b/exper/sim12/start_sbatch_power.sh @@ -1,11 +1,15 @@ #!/bin/bash -for ip in "1" "2" "3" +## submits one job per dimension and signal strength; arguments are passed on to +## sbatch, e.g. +## ./start_sbatch_power.sh -p + +for ip in 1 2 3 do - for ialpha in {1..11} - do - echo "sbatch ${ip} ${ialpha} ..." - sbatch --job-name=diff_sim12_${ip}_${ialpha}_power sbatch_sim12_power ${ip} ${ialpha} - echo "... done" - done + for ialpha in {1..11} + do + echo "sbatch ${ip} ${ialpha} ..." + sbatch "$@" --job-name=diff_sim12_${ip}_${ialpha}_power sbatch_sim12_power ${ip} ${ialpha} + echo "... done" + done done diff --git a/exper/sim2/sbatch_sim2 b/exper/sim2/sbatch_sim2 index 7f61b89..a2c45f6 100644 --- a/exper/sim2/sbatch_sim2 +++ b/exper/sim2/sbatch_sim2 @@ -1,20 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim2_${1}_${2} sbatch_sim2 +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim2_${ip}_${iElem} sbatch_sim2 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim2, logs to +## $SCRATCH/output/sim2. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive -module load parallel -module load hdf5 +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim2 +LOG=$SCRATCH/output/sim2 -SCRATCH_PATH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +mkdir -p $OUT +mkdir -p $LOG -mkdir -p $SCRATCH_PATH/diffTest/sim2 -mkdir -p $SCRATCH_PATH/output/sim2 +command -v parallel > /dev/null || module load parallel -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH_PATH/output/sim2/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH_PATH/output/sim2/$SLURM_JOB_NAME.out.{1} $JULIA sim2.jl {1} ${1} ${2} $SCRATCH_PATH/diffTest/sim2 ::: {1..1000} +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' + +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim2.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim2/sim2.jl b/exper/sim2/sim2.jl index ffb803f..aecdffe 100644 --- a/exper/sim2/sim2.jl +++ b/exper/sim2/sim2.jl @@ -1,9 +1,13 @@ +# Coverage of confidence intervals for one entry of Δ when Ωx = Ωy (Δ = 0) with +# Ω[i, j] = 0.9^|i - j|. +# +# julia --project=.. sim2.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -42,8 +46,8 @@ dist_Y = MvNormal(convert(Matrix, Σ)) # generate data Random.seed!(134 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -64,4 +68,4 @@ ri, ci = elemArr[iElem] @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim2/sim2_test.jl b/exper/sim2/sim2_test.jl deleted file mode 100644 index 5d433d8..0000000 --- a/exper/sim2/sim2_test.jl +++ /dev/null @@ -1,73 +0,0 @@ -using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD - - -@show gethostname() - -pArr = [25, 100, 200] -elemArr = [(5,5), (8, 7), (50, 25)] -n = 300 -est = Array{Any}(undef, 9) # number of methods - - -# rep = parse(Int,ARGS[1]) -# ip = parse(Int,ARGS[2]) -# iElem = parse(Int,ARGS[3]) -# dir = ARGS[4] - -rep = 1 -ip = 1 -iElem = 1 -dir = "." - - -Random.seed!(134) - -p = pArr[ip] -# generate model -Ω = Matrix{Float64}(I, p, p) -ρ = 0.9 -for k=1:p-1 - for l=1:p-k - Ω[l , l+k] = ρ^k - Ω[l+k, l ] = ρ^k - end -end -d = Vector{Float64}(undef, p) -rand!(Uniform(0.5, 2.5), d) -d .= sqrt.(d) -D = Diagonal(d) -Σ = inv(Symmetric(D * Ω * D)) - -dist_X = MvNormal(convert(Matrix, Σ)) -dist_Y = MvNormal(convert(Matrix, Σ)) - -# generate data -Random.seed!(134 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' - -Sx = Symmetric( X'X / n ) -Sy = Symmetric( Y'Y / n ) - -@show (ip, iElem) -supportOracle = Vector{Int64}() -ri, ci = elemArr[iElem] - -@time est[1], eΔ, esuppΔ, eω, esuppω, eSupport = DiffPrecTest.estimate(SymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy) -@time est[2], _, _, _, _, _ = DiffPrecTest.estimate(AsymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy, Δ=eΔ, suppΔ=esuppΔ) -@time est[3] = DiffPrecTest.estimate(SeparateNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy) -@time est[4] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport) -@time est[5] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport) - -@time est[6] = DiffPrecTest.estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle) -@time est[7] = DiffPrecTest.estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle) -@time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) - - -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est diff --git a/exper/sim2/start_sbatch.sh b/exper/sim2/start_sbatch.sh old mode 100644 new mode 100755 index 9f959cf..de47687 --- a/exper/sim2/start_sbatch.sh +++ b/exper/sim2/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" + for iElem in 1 2 3 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim2_${ip}_${iElem} sbatch_sim2 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim2_${ip}_${iElem} sbatch_sim2 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim2/summarize.jl b/exper/sim2/summarize.jl index 920f511..2b7fc63 100644 --- a/exper/sim2/summarize.jl +++ b/exper/sim2/summarize.jl @@ -1,73 +1,23 @@ -using JLD -using DiffPrecTest -using DataFrames, DataFramesMeta, CSV -using LaTeXStrings, LaTeXTabulars +# Coverage, length and bias of the confidence intervals in sim2 (Δ = 0), from +# the files results/sim2_res_ip_iElem.jld2 written by collect.sh. Writes +# sim2.csv and sim2_table.tex. +# +# julia --project=.. summarize.jl +using CSV +include("../summary.jl") -cd("exper/sim2") +cd(@__DIR__) pArr = [100, 200] elemArr = [(5,5), (8, 7), (50, 25)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) +df = newSummary() for ip=1:2 for iElem=1:3 - res = load("results/sim2_res_$(ip)_$(iElem).jld", "results") - - for j in 1:9 - global pArr - global elemArr - global methodArr - global df - - sr = computeSimulationResult([res[j, i] for i=1:1000], 0.) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end + summarizeSetting!(df, "results/sim2_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], 0.) end end CSV.write("sim2.csv", df) - -### create a table - -latex_tabular("sim2_table.tex", - Tabular("lllccc"), - [ - Rule(:top), - ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"], - Rule(:mid), - vcat([ - hcat( - (a = fill("", 9, 1); - b = latexstring("p = $(pArr[ip])"); - a[1] = "\\multirow{9}{*}{$(b)}"; - a - ), - vcat([ - hcat( - (a = fill("", 3, 1); - b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"); - a[1] = "\\multirow{3}{*}{$(b)}"; - a - ), - Matrix{Any}( - @linq df |> - where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"), - :p .== pArr[ip], - :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |> select(:method, - coverage = convert.(Int, :coverage * 1000), - lenCI = round.(:lenCI, digits=3), - bias = round.(:bias * 1000, digits=1)) - ) - ) - for iElem=1:3 - ]... - ) - ) - for ip=1:2 - ]...), - Rule(:bottom) - ] -) \ No newline at end of file +coverageTable("sim2_table.tex", df, pArr[1:2], elemArr, zeros(2, 3)) diff --git a/exper/sim3/sbatch_sim3 b/exper/sim3/sbatch_sim3 index eb7adf2..cea3dfa 100644 --- a/exper/sim3/sbatch_sim3 +++ b/exper/sim3/sbatch_sim3 @@ -1,20 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim3_${ip}_${iElem} sbatch_sim3 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim3, logs to +## $SCRATCH/output/sim3. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim3 +LOG=$SCRATCH/output/sim3 -SCRATCH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +mkdir -p $OUT +mkdir -p $LOG -mkdir -p $SCRATCH/diffTest/sim3 -mkdir -p $SCRATCH/output/sim3 +command -v parallel > /dev/null || module load parallel -module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH/output/sim3/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH/output/sim3/$SLURM_JOB_NAME.out.{1} $JULIA sim3.jl {1} ${1} ${2} $SCRATCH/diffTest/sim3 ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim3.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim3/sim3.jl b/exper/sim3/sim3.jl index ba5ef3c..6eba668 100644 --- a/exper/sim3/sim3.jl +++ b/exper/sim3/sim3.jl @@ -1,9 +1,13 @@ +# Coverage of confidence intervals for one entry of Δ, with Ωx banded and Δ +# nonzero on five diagonal and five off-diagonal entries. +# +# julia --project=.. sim3.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -58,8 +62,8 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(1234 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -79,4 +83,4 @@ ri, ci = elemArr[iElem] @time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim3/start_sbatch.sh b/exper/sim3/start_sbatch.sh old mode 100644 new mode 100755 index 23e5df5..5de1ad8 --- a/exper/sim3/start_sbatch.sh +++ b/exper/sim3/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" "4" "5" + for iElem in 1 2 3 4 5 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim3_${ip}_${iElem} sbatch_sim3 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim3_${ip}_${iElem} sbatch_sim3 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim3/summarize.jl b/exper/sim3/summarize.jl index 78afe36..baf1bd6 100644 --- a/exper/sim3/summarize.jl +++ b/exper/sim3/summarize.jl @@ -1,118 +1,58 @@ -using JLD -using DiffPrecTest -using DataFrames, DataFramesMeta, CSV -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using Random, Distributions -using LaTeXStrings, LaTeXTabulars +# Coverage, length and bias of the confidence intervals in sim3, from the files +# results/sim3_res_ip_iElem.jld2 written by collect.sh. Writes sim3.csv and +# sim3_table.tex. +# +# julia --project=.. summarize.jl +using CSV +using LinearAlgebra, Random +using Distributions -cd("exper/sim3") +include("../summary.jl") + +cd(@__DIR__) pArr = [100, 200] elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) - -trueValue = zeros(Float64, 2, 5) +# the true difference, generated as in sim3.jl +function trueDelta(p) + Random.seed!(1234) + Ωx = Matrix{Float64}(I, p, p) + for l=1:p-1 + Ωx[l , l+1] = 0.6 + Ωx[l+1, l ] = 0.6 + end + for l=1:p-2 + Ωx[l , l+2] = 0.3 + Ωx[l+2, l ] = 0.3 + end + Δ = zeros(p, p) + for j=1:5 + i = 5 + (j-1) + Δ[i, i] = rand(Uniform(0.1, 0.2)) + end + for j=1:5 + i = 5 + (j-1) + v = rand(Uniform(0.2, 0.5)) + Δ[i , i+1] = v + Δ[i+1, i ] = v + end + d = Vector{Float64}(undef, p) + rand!(Uniform(0.5, 2.5), d) + d .= sqrt.(d) + D = Diagonal(d) + D * Δ * D +end +df = newSummary() +trueValue = zeros(Float64, 2, length(elemArr)) for ip=1:2 - for iElem=1:5 - global trueValue - - p = pArr[ip] - Random.seed!(1234) - - # generate model - Ωx = Matrix{Float64}(I, p, p) - for l=1:p-1 - Ωx[l , l+1] = 0.6 - Ωx[l+1, l ] = 0.6 - end - for l=1:p-2 - Ωx[l , l+2] = 0.3 - Ωx[l+2, l ] = 0.3 - end - Δ = zeros(p, p) - # generate Delta - for j=1:5 - i = 5 + (j-1) - Δ[i, i] = rand(Uniform(0.1, 0.2)) - end - for j=1:5 - i = 5 + (j-1) - v = rand(Uniform(0.2, 0.5)) - Δ[i , i+1] = v - Δ[i+1, i ] = v - end - Ωy = Ωx - Δ - d = Vector{Float64}(undef, p) - rand!(Uniform(0.5, 2.5), d) - d .= sqrt.(d) - D = Diagonal(d) - Σx = inv(Symmetric(D * Ωx * D)) - Σy = inv(Symmetric(D * Ωy * D)) - - tΔ = D * Δ * D - - ri, ci = elemArr[iElem] - indE = (ci - 1) * p + ri - - trueValue[ip, iElem] = tΔ[ri, ci] - - res = load("results/sim3_res_$(ip)_$(iElem).jld", "results") - - for j in 1:9 - sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE]) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end + tΔ = trueDelta(pArr[ip]) + for iElem in eachindex(elemArr) + trueValue[ip, iElem] = tΔ[elemArr[iElem]...] + summarizeSetting!(df, "results/sim3_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], trueValue[ip, iElem]) end end - CSV.write("sim3.csv", df) - - - -### create a table - -latex_tabular("sim3_table.tex", - Tabular("lllccc"), - [ - Rule(:top), - ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"], - Rule(:mid), - vcat([ - hcat( - (a = fill("", 15, 1); - b = latexstring("p = $(pArr[ip])"); - a[1] = "\\multirow{15}{*}{$(b)}"; - a - ), - vcat([ - hcat( - (a = fill("", 3, 1); - b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = $(trueValue[ip, iElem])"); - a[1] = "\\multirow{3}{*}{$(b)}"; - a - ), - Matrix{Any}( - @linq df |> - where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"), - :p .== pArr[ip], - :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |> select(:method, - coverage = convert.(Int, :coverage * 1000), - lenCI = round.(:lenCI, digits=3), - bias = round.(:bias * 1000, digits=1)) - ) - ) - for iElem=1:5 - ]... - ) - ) - for ip=1:2 - ]...), - Rule(:bottom) - ] -) +coverageTable("sim3_table.tex", df, pArr[1:2], elemArr, trueValue) diff --git a/exper/sim3/test.jl b/exper/sim3/test.jl deleted file mode 100644 index 30883a1..0000000 --- a/exper/sim3/test.jl +++ /dev/null @@ -1,151 +0,0 @@ -using Revise -using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD - -#using Plots - -pArr = [13, 200] -elemArr = [(5,5), (3, 2), (50, 25), (21, 20), (30, 30)] -n = 300 -est = Array{Any}(undef, 5) # number of methods - -out = Vector{DiffPrecResultNormal}(undef, 1000) - -ip = 1 -iElem = 2 - -p = pArr[ip] -Random.seed!(1234) - -# generate model -Ωx = Matrix{Float64}(I, p, p) -for l=1:p-1 - Ωx[l , l+1] = 0.6 - Ωx[l+1, l ] = 0.6 -end -for l=1:p-2 - Ωx[l , l+2] = 0.3 - Ωx[l+2, l ] = 0.3 -end -Δ = zeros(p, p) -# generate Delta -for j=1:p - #i = 5 + (j-1) - i = j - Δ[i, i] = rand(Uniform(0.1, 0.2)) -end -for j=1:p-1 - #i = 5 + (j-1) - i = j - v = rand(Uniform(0.2, 0.5)) - Δ[i , i+1] = v - Δ[i+1, i ] = v -end -Ωy = Ωx - Δ -d = Vector{Float64}(undef, p) -rand!(Uniform(0.5, 2.5), d) -d .= sqrt.(d) -D = Diagonal(d) -Σx = inv(Symmetric(D * Ωx * D)) -Σy = inv(Symmetric(D * Ωy * D)) - -tΔ = D * Δ * D - -dist_X = MvNormal(convert(Matrix, Σx)) -dist_Y = MvNormal(convert(Matrix, Σy)) - -ri, ci = elemArr[iElem] -indE = sub2indLowerTriangular(p, ri, ci) - -indOracle = Vector{Int64}() -for ci=1:p - for ri=ci:p - global indOracle - - if !iszero(tΔ[ri, ci]) - push!(indOracle, sub2indLowerTriangular(p, ri, ci)) - end - end -end - - -Random.seed!(1234) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' - -# Sx = Symmetric(cov(X)) -# Sy = Symmetric(cov(Y)) - -# C = DiffPrecTest.skron(Sx, Sy) -# b = DiffPrecTest.svec(Sy) - DiffPrecTest.svec(Sx) - -# f = CDQuadraticLoss(C, b) - - -@time res1, _ = DiffPrecTest.estimate(ReducedNormal(), X, Y, indE); -@time res2, _ = DiffPrecTest.estimate(SymmetricNormal(), X, Y, 3, 2); - -for rep = 1:1000 - global out - global indOracle - global indE - - # generate data - # Random.seed!(1234 + rep) - X = rand(dist_X, n)' - Y = rand(dist_Y, n)' - - #@time out[rep], _, _, _, _, indS = DiffPrecTest.estimate(ReducedNormal(), X, Y, indE) - #@time out[rep], _, _, _, _, indS = DiffPrecTest.estimate(SymmetricNormal(), X, Y, 3, 2) - - @time out[rep] = DiffPrecTest.estimate(ReducedOracleNormal(), Symmetric(cov(X)), Symmetric(cov(Y)), X, Y, indE, indOracle) -end - -computeSimulationResult(out, tΔ[elemArr[iElem]...]) - -true_v = tΔ[elemArr[iElem]...] -histogram( [(out[i].p - true_v) / out[i].std for i = 1:1000] ) - - - -# @time est[2] = DiffPrecTest.estimate(SeparateNormal(), X, Y, indE) -# @time est[3] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indS) -# @time est[4] = DiffPrecTest.estimate(SymmetricOracleNormal(), Symmetric(cov(X)), n, Symmetric(cov(Y)), n, indOracle) -# @time est[5] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indOracle) - - - - -###################### -# test variance - -Random.seed!(1234+5) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' - -Sx = Symmetric(cov(X)) -Sy = Symmetric(cov(Y)) - -indS = indOracle - -C = DiffPrecTest.skron(Sx, Sy, indS) -b = DiffPrecTest.svec(Sy, indS) - DiffPrecTest.svec(Sx, indS) - -tmp = zeros(length(indS)) -tmp[1] = 1. -ω = C \ tmp -Δ = C \ b - -@time v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, indS, Δ, indS) - -@time variance(DiffPrecTest.createVarianceReducedEstim(ω, indS, Δ, indS), X, Y) - - - - - - diff --git a/exper/sim4/sbatch_sim4 b/exper/sim4/sbatch_sim4 index 7e02af5..ac3ed62 100644 --- a/exper/sim4/sbatch_sim4 +++ b/exper/sim4/sbatch_sim4 @@ -1,22 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir - +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim4_${ip}_${iElem} sbatch_sim4 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim4, logs to +## $SCRATCH/output/sim4. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive - -SCRATCH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim4 +LOG=$SCRATCH/output/sim4 -mkdir -p $SCRATCH/diffTest/sim4 -mkdir -p $SCRATCH/output/sim4 +mkdir -p $OUT +mkdir -p $LOG -module load parallel +command -v parallel > /dev/null || module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH/output/sim4/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH/output/sim4/$SLURM_JOB_NAME.out.{1} $JULIA sim4.jl {1} ${1} ${2} $SCRATCH/diffTest/sim4 ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim4.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim4/sim4.jl b/exper/sim4/sim4.jl index 7d583ca..356c0ec 100644 --- a/exper/sim4/sim4.jl +++ b/exper/sim4/sim4.jl @@ -1,9 +1,13 @@ +# Coverage of confidence intervals for one entry of Δ, with Ωx[i, j] = 0.7^|i - j| +# and Δ nonzero on five diagonal and five off-diagonal entries. +# +# julia --project=.. sim4.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -57,8 +61,8 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(1234 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -78,4 +82,4 @@ ri, ci = elemArr[iElem] @time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim4/start_sbatch.sh b/exper/sim4/start_sbatch.sh old mode 100644 new mode 100755 index b499dd4..7aaf136 --- a/exper/sim4/start_sbatch.sh +++ b/exper/sim4/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" "4" "5" + for iElem in 1 2 3 4 5 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim4_${ip}_${iElem} sbatch_sim4 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim4_${ip}_${iElem} sbatch_sim4 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim4/summarize.jl b/exper/sim4/summarize.jl index 5f0fd1c..639677c 100644 --- a/exper/sim4/summarize.jl +++ b/exper/sim4/summarize.jl @@ -1,113 +1,57 @@ -using JLD -using DiffPrecTest -using DataFrames, DataFramesMeta, CSV -using Random, LinearAlgebra, Distributions -using LaTeXStrings, LaTeXTabulars +# Coverage, length and bias of the confidence intervals in sim4, from the files +# results/sim4_res_ip_iElem.jld2 written by collect.sh. Writes sim4.csv and +# sim4_table.tex. +# +# julia --project=.. summarize.jl +using CSV +using LinearAlgebra, Random +using Distributions -cd("exper/sim4") +include("../summary.jl") -pArr = [100, 200] -elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) - -trueValue = zeros(Float64, 2, 5) +cd(@__DIR__) -for ip=1:2 - for iElem=1:5 - p = pArr[ip] - Random.seed!(1234) +pArr = [100, 200, 500] +elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)] - # generate model - Ωx = Matrix{Float64}(I, p, p) - ρ = 0.7 - for k=1:p-1 - for l=1:p-k - Ωx[l , l+k] = ρ^k - Ωx[l+k, l ] = ρ^k - end - end - Δ = zeros(p, p) - # generate Delta - for j=1:5 - i = 5 + (j-1) - Δ[i, i] = rand(Uniform(0.1, 0.2)) +# the true difference, generated as in sim4.jl +function trueDelta(p) + Random.seed!(1234) + Ωx = Matrix{Float64}(I, p, p) + ρ = 0.7 + for k=1:p-1 + for l=1:p-k + Ωx[l , l+k] = ρ^k + Ωx[l+k, l ] = ρ^k end - for j=1:5 - i = 5 + (j-1) - v = rand(Uniform(0.2, 0.5)) - Δ[i , i+1] = v - Δ[i+1, i ] = v - end - Ωy = Ωx - Δ - d = Vector{Float64}(undef, p) - rand!(Uniform(0.5, 2.5), d) - d .= sqrt.(d) - D = Diagonal(d) - Σx = inv(Symmetric(D * Ωx * D)) - Σy = inv(Symmetric(D * Ωy * D)) - - tΔ = D * Δ * D - - res = load("results/sim4_res_$(ip)_$(iElem).jld", "results") - - ri, ci = elemArr[iElem] - indE = (ci - 1) * p + ri - - trueValue[ip, iElem] = tΔ[ri, ci] + end + Δ = zeros(p, p) + for j=1:5 + i = 5 + (j-1) + Δ[i, i] = rand(Uniform(0.1, 0.2)) + end + for j=1:5 + i = 5 + (j-1) + v = rand(Uniform(0.2, 0.5)) + Δ[i , i+1] = v + Δ[i+1, i ] = v + end + d = Vector{Float64}(undef, p) + rand!(Uniform(0.5, 2.5), d) + d .= sqrt.(d) + D = Diagonal(d) + D * Δ * D +end - for j in 1:9 - sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE]) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end +df = newSummary() +trueValue = zeros(Float64, 2, length(elemArr)) +for ip=1:2 + tΔ = trueDelta(pArr[ip]) + for iElem in eachindex(elemArr) + trueValue[ip, iElem] = tΔ[elemArr[iElem]...] + summarizeSetting!(df, "results/sim4_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], trueValue[ip, iElem]) end end CSV.write("sim4.csv", df) - - - -### create a table - -latex_tabular("sim4_table.tex", - Tabular("lllccc"), - [ - Rule(:top), - ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"], - Rule(:mid), - vcat([ - hcat( - (a = fill("", 15, 1); - b = latexstring("p = $(pArr[ip])"); - a[1] = "\\multirow{15}{*}{$(b)}"; - a - ), - vcat([ - hcat( - (a = fill("", 3, 1); - b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = $(trueValue[ip, iElem])"); - a[1] = "\\multirow{3}{*}{$(b)}"; - a - ), - Matrix{Any}( - @linq df |> - where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"), - :p .== pArr[ip], - :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |> select(:method, - coverage = convert.(Int, :coverage * 1000), - lenCI = round.(:lenCI, digits=3), - bias = round.(:bias * 1000, digits=1)) - ) - ) - for iElem=1:5 - ]... - ) - ) - for ip=1:2 - ]...), - Rule(:bottom) - ] -) - - +coverageTable("sim4_table.tex", df, pArr[1:2], elemArr, trueValue) diff --git a/exper/sim5/model.jl b/exper/sim5/model.jl new file mode 100644 index 0000000..75f6fe3 --- /dev/null +++ b/exper/sim5/model.jl @@ -0,0 +1,18 @@ +# Model of sim5: Ωx is block diagonal with banded blocks of sizes p/2, p/4 and +# p/4 scaled by 1, 2 and 4, and Ωy is the diagonal of Ωx, so that Δ = Ωx - Ωy +# is the off-diagonal part of Ωx. p must be divisible by 4. +function sim5model(p) + Ωx = Matrix{Float64}(I, p, p) + mf = [1., 0.5, 0.4] + for (α, bp, ep) in ((1., 1, div(p, 2)), (2., div(p, 2) + 1, div(3p, 4)), (4., div(3p, 4) + 1, p)) + for k=0:2 + v = α * mf[k+1] + for l=bp:ep-k + Ωx[l , l+k] = v + Ωx[l+k, l ] = v + end + end + end + Ωy = Matrix(Diagonal(Ωx)) + Ωx, Ωy +end diff --git a/exper/sim5/sbatch_sim5 b/exper/sim5/sbatch_sim5 index 87904a5..d1dfedf 100644 --- a/exper/sim5/sbatch_sim5 +++ b/exper/sim5/sbatch_sim5 @@ -1,23 +1,30 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir - +## submit as +## sbatch sbatch_sim5 +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim5a, logs to +## $SCRATCH/output/sim5a. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 #SBATCH --job-name=diff_sim5a -##SBATCH --exclusive - -SCRATCH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim5a +LOG=$SCRATCH/output/sim5a +mkdir -p $OUT +mkdir -p $LOG -mkdir -p $SCRATCH/diffTest/sim5a -mkdir -p $SCRATCH/output/sim5a +command -v parallel > /dev/null || module load parallel -module load parallel +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH/output/sim5a/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH/output/sim5a/$SLURM_JOB_NAME.out.{1} $JULIA sim5a.jl {1} $SCRATCH/diffTest/sim5a ::: {1..100} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim5a.jl {1} $OUT ::: {1..100} diff --git a/exper/sim5/sbatch_sim5b b/exper/sim5/sbatch_sim5b index 6c4a2d9..2eca64f 100644 --- a/exper/sim5/sbatch_sim5b +++ b/exper/sim5/sbatch_sim5b @@ -1,23 +1,30 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir - +## submit as +## sbatch sbatch_sim5b +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim5b, logs to +## $SCRATCH/output/sim5b. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 #SBATCH --job-name=diff_sim5b -##SBATCH --exclusive +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim5b +LOG=$SCRATCH/output/sim5b -module load parallel -module load hdf5 +mkdir -p $OUT +mkdir -p $LOG -SCRATCH_PATH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +command -v parallel > /dev/null || module load parallel -mkdir -p $SCRATCH_PATH/diffTest/sim5b -mkdir -p $SCRATCH_PATH/output/sim5b +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH_PATH/output/sim5b/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH_PATH/output/sim5b/$SLURM_JOB_NAME.out.{1} $JULIA sim5b.jl {1} $SCRATCH_PATH/diffTest/sim5b ::: {1..100} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim5b.jl {1} $OUT ::: {1..100} diff --git a/exper/sim5/sim5a.jl b/exper/sim5/sim5a.jl index 51e57f1..0d7d010 100644 --- a/exper/sim5/sim5a.jl +++ b/exper/sim5/sim5a.jl @@ -1,9 +1,16 @@ +# Support recovery with p = 60 (see model.jl): thresholded refits with normal +# (ANTSupport) and bootstrap (BootStdSupport) standard errors, and the D-trace +# estimator path with the penalty chosen on validation data. +# +# julia --project=.. sim5a.jl rep dir +# +# saves the estimates to dir/res_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random, Statistics +using Distributions +using JLD2 + +include("model.jl") rep = parse(Int,ARGS[1]) dir = ARGS[2] @@ -12,42 +19,7 @@ p = 60 n = 300 # generate model -Ωx = Matrix{Float64}(I, p, p) -mf = [1., 0.5, 0.4] -# block 1 -α = 1. -bp = 1 -ep = 30 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 2 -α = 2. -bp = 31 -ep = 45 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 3 -α = 4. -bp = 46 -ep = 60 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -Ωy = Matrix(Diagonal(Ωx)) +Ωx, Ωy = sim5model(p) Δ = Ωx - Ωy Σx = inv(Symmetric(Ωx)) @@ -58,10 +30,10 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(2345 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' -Xtest = rand(dist_X, n)' # validation data for -Ytest = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) +Xtest = permutedims(rand(dist_X, n)) # validation data for the D-trace estimator +Ytest = permutedims(rand(dist_Y, n)) ################# # @@ -88,4 +60,4 @@ minλ = maxλ * 0.04 eΔDTr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) -@save "$(dir)/res_$(rep).jld" 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 eS eΔDTr i2 iInf loss2arr lossInfarr diff --git a/exper/sim5/sim5b.jl b/exper/sim5/sim5b.jl index 7eabed2..5617e58 100644 --- a/exper/sim5/sim5b.jl +++ b/exper/sim5/sim5b.jl @@ -1,9 +1,16 @@ +# Support recovery with p = 100 (see model.jl): thresholded refits with normal +# (ANTSupport) and bootstrap (BootStdSupport) standard errors, and the D-trace +# estimator path with the penalty chosen on validation data. +# +# julia --project=.. sim5b.jl rep dir +# +# saves the estimates to dir/res_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random, Statistics +using Distributions +using JLD2 + +include("model.jl") rep = parse(Int,ARGS[1]) dir = ARGS[2] @@ -12,42 +19,7 @@ p = 100 n = 300 # generate model -Ωx = Matrix{Float64}(I, p, p) -mf = [1., 0.5, 0.4] -# block 1 -α = 1. -bp = 1 -ep = 50 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 2 -α = 2. -bp = 51 -ep = 75 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 3 -α = 4. -bp = 76 -ep = 100 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -Ωy = Matrix(Diagonal(Ωx)) +Ωx, Ωy = sim5model(p) Δ = Ωx - Ωy Σx = inv(Symmetric(Ωx)) @@ -58,10 +30,10 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(2345 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' -Xtest = rand(dist_X, n)' # validation data for -Ytest = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) +Xtest = permutedims(rand(dist_X, n)) # validation data for the D-trace estimator +Ytest = permutedims(rand(dist_Y, n)) ################# # @@ -88,4 +60,4 @@ minλ = maxλ * 0.04 eΔDTr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr) -@save "$(dir)/res_$(rep).jld" 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 eS eΔDTr i2 iInf loss2arr lossInfarr diff --git a/exper/sim5/summarize.jl b/exper/sim5/summarize.jl index 6ae0daa..046ff6e 100644 --- a/exper/sim5/summarize.jl +++ b/exper/sim5/summarize.jl @@ -1,178 +1,79 @@ -using Revise - -using JLD -using DiffPrecTest +# Support recovery in sim5: false and true positive rates, averaged over the +# repetitions, of the thresholded refits along the thresholds τ and of the +# D-trace estimator along its penalty path, and the rates of the D-trace +# estimator at the penalties chosen on validation data. +# +# julia --project=.. summarize.jl dir variant +# +# where dir holds the files res_rep.jld2 written by sim5a.jl (variant a) or +# sim5b.jl (variant b). Writes sim5_roc.csv, sim5_dtrace.csv +# and sim5_roc.pdf. +using LinearAlgebra, Statistics +using CovSel: getConfusionMatrix, tpr, fpr, precision using DataFrames, CSV -using Random, LinearAlgebra, Distributions, SparseArrays, Statistics -using CovSel +using JLD2 using Plots -dir = "/home/mkolar/.julia/dev/DiffPrecTest/exper/sim5/res" +include("model.jl") -p = 20 -n = 300 +dir = ARGS[1] +variant = ARGS[2] +p = Dict("a" => 60, "b" => 100)[variant] -# generate model -Ωx = Matrix{Float64}(I, p, p) -mf = [1., 0.5, 0.4] -# block 1 -α = 1. -bp = 1 -ep = 10 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 2 -α = 2. -bp = 11 -ep = 15 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -# block 3 -α = 4. -bp = 16 -ep = 20 -for k=0:2 - v = α * mf[k+1] - for l=bp:ep-k - Ωx[l , l+k] = v - Ωx[l+k, l ] = v - end -end -Ωy = Matrix(Diagonal(Ωx)) +Ωx, Ωy = sim5model(p) Δ = Ωx - Ωy -Σx = inv(Symmetric(Ωx)) -Σy = inv(Symmetric(Ωy)) - -dist_X = MvNormal(convert(Matrix, Σx)) -dist_Y = MvNormal(convert(Matrix, Σy)) - - -τArr = collect( range(10,0,length=500) ) - - - -NUM_REP=200 -FPR_mine = zeros(NUM_REP, length(τArr)) -TPR_mine = zeros(NUM_REP, length(τArr)) -FPR_tr = zeros(NUM_REP, 50) -TPR_tr = zeros(NUM_REP, 50) +files = filter(f -> occursin(r"^res_\d+\.jld2$", f), readdir(dir)) +NUM_REP = length(files) +@show NUM_REP + +numτ = 500 +numΛ = 50 +FPR_normal = zeros(NUM_REP, numτ) +TPR_normal = zeros(NUM_REP, numτ) +FPR_boot = zeros(NUM_REP, numτ) +TPR_boot = zeros(NUM_REP, numτ) +FPR_tr = zeros(NUM_REP, numΛ) +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") + end -for rep=1:NUM_REP - global FPR_mine, TPR_mine, FPR_tr, TPR_tr, Δ, τArr - fname = "$(dir)/res_$(rep).jld" + for j=1:numτ + cm = getConfusionMatrix(Δ, eΔNormal[j]) + FPR_normal[rep, j] = fpr(cm) + TPR_normal[rep, j] = tpr(cm) - file = jldopen(fname, "r") - eΔNormal = read(file, "eΔNormal") - eΔDTr = read(file, "eΔDTr") - close(file) + cm = getConfusionMatrix(Δ, eΔBoot[j]) + FPR_boot[rep, j] = fpr(cm) + TPR_boot[rep, j] = tpr(cm) + end - for j=1:length(τArr) - confusionMatrix = CovSel.getConfusionMatrix(Δ, eΔNormal[j]) - FPR_mine[rep, j] = CovSel.fpr(confusionMatrix) - TPR_mine[rep, j] = CovSel.tpr(confusionMatrix) + for j=1:numΛ + cm = getConfusionMatrix(Δ, eΔDTr[j]) + FPR_tr[rep, j] = fpr(cm) + TPR_tr[rep, j] = tpr(cm) end - for j=1:50 - confusionMatrix = CovSel.getConfusionMatrix(Δ, eΔDTr[j]) - FPR_tr[rep, j] = CovSel.fpr(confusionMatrix) - TPR_tr[rep, j] = CovSel.tpr(confusionMatrix) + for (selection, i) in (("frobenius", i2), ("max", iInf)) + cm = getConfusionMatrix(Δ, eΔDTr[i]) + push!(dtrace, (selection, rep, tpr(cm), fpr(cm), precision(cm))) end end -pyplot() -plot(mean(FPR_mine, dims=1)', mean(TPR_mine, dims=1)') -plot!(mean(FPR_tr, dims=1)', mean(TPR_tr, dims=1)') - - -plot(FPR_mine[1, :], TPR_mine[1, :]) -FPR_mine[1, :] -rep = 1 -fname = "$(dir)/res_$(rep).jld" -file = jldopen(fname, "r") -eΔNormal = read(file, "eΔNormal") - - -# using JLD -# using DiffPrecTest -# using DataFrames, CSV -# using Random, LinearAlgebra, Distributions, SparseArrays -# using CovSel -# -# dir = ARGS[1] -# -# p = 100 -# n = 300 -# -# # generate model -# Ωx = Matrix{Float64}(I, p, p) -# ρ = 0.5 -# for k=1:p-1 -# for l=1:p-k -# Ωx[l , l+k] = ρ^k -# Ωx[l+k, l ] = ρ^k -# end -# end -# -# Ωy = copy(Ωx) -# k = div(p, 4) -# for l=1:p-k -# Ωy[l , l+k] = 0.2 -# Ωy[l+k, l ] = 0.2 -# end -# Δ = Ωx - Ωy -# -# NUM_REP = 100 -# resN = Array{CovSel.ConfusionMatrix}(undef, NUM_REP) -# resBoot = Array{CovSel.ConfusionMatrix}(undef, NUM_REP) -# resDTr2 = Array{CovSel.ConfusionMatrix}(undef, NUM_REP) -# resDTrInf = Array{CovSel.ConfusionMatrix}(undef, NUM_REP) -# i = 0 -# for rep=1:NUM_REP -# global i, resN, resBoot, resDTr2, resDTrInf -# fname = "$(dir)/res_$(rep).jld" -# -# try -# i += 1 -# file = jldopen(fname, "r") -# eΔNormal = read(file, "eΔNormal") -# eΔBoot = read(file, "eΔBoot") -# eΔDTr = read(file, "eΔDTr") -# i2 = read(file, "i2") -# iInf = read(file, "iInf") -# close(file) -# -# resN[i] = CovSel.getConfusionMatrix(Δ, eΔNormal) -# resBoot[i] = CovSel.getConfusionMatrix(Δ, eΔBoot) -# resDTr2[i] = CovSel.getConfusionMatrix(Δ, eΔDTr[i2]) -# resDTrInf[i] = CovSel.getConfusionMatrix(Δ, eΔDTr[iInf]) -# catch -# @show fname -# end -# end -# -# @show mean([CovSel.tpr(resN[j]) for j=1:i]) -# @show mean([CovSel.fpr(resN[j]) for j=1:i]) -# @show mean([CovSel.precision(resN[j]) for j=1:i]) -# -# @show mean([CovSel.tpr(resBoot[j]) for j=1:i]) -# @show mean([CovSel.fpr(resBoot[j]) for j=1:i]) -# @show mean([CovSel.precision(resBoot[j]) for j=1:i]) -# -# @show mean([CovSel.tpr(resDTr2[j]) for j=1:i]) -# @show mean([CovSel.fpr(resDTr2[j]) for j=1:i]) -# @show mean([CovSel.precision(resDTr2[j]) for j=1:i]) -# -# @show mean([CovSel.tpr(resDTrInf[j]) for j=1:i]) -# @show mean([CovSel.fpr(resDTrInf[j]) for j=1:i]) -# @show mean([CovSel.precision(resDTrInf[j]) for j=1:i]) +roc = vcat( + DataFrame(method="Normal", fpr=vec(mean(FPR_normal; dims=1)), tpr=vec(mean(TPR_normal; dims=1))), + DataFrame(method="Boot", fpr=vec(mean(FPR_boot; dims=1)), tpr=vec(mean(TPR_boot; dims=1))), + DataFrame(method="DTrace", fpr=vec(mean(FPR_tr; dims=1)), tpr=vec(mean(TPR_tr; dims=1))), +) +CSV.write("sim5$(variant)_roc.csv", roc) +CSV.write("sim5$(variant)_dtrace.csv", combine(groupby(dtrace, :selection), [:tpr, :fpr, :precision] .=> mean)) + +plt = plot(xlabel="false positive rate", ylabel="true positive rate", legend=:bottomright, title="p = $(p)") +for g in groupby(roc, :method; sort=false) + plot!(plt, g.fpr, g.tpr; label=first(g.method)) +end +savefig(plt, "sim5$(variant)_roc.pdf") diff --git a/exper/sim6/sbatch_sim6 b/exper/sim6/sbatch_sim6 index 588c041..93e3b76 100644 --- a/exper/sim6/sbatch_sim6 +++ b/exper/sim6/sbatch_sim6 @@ -1,19 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim6_${ip}_${iElem} sbatch_sim6 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim6, logs to +## $SCRATCH/output/sim6. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive -SCRATCH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim6 +LOG=$SCRATCH/output/sim6 -mkdir -p $SCRATCH/diffTest/sim6 -mkdir -p $SCRATCH/output/sim6 +mkdir -p $OUT +mkdir -p $LOG -module load parallel +command -v parallel > /dev/null || module load parallel -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH/output/sim6/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH/output/sim6/$SLURM_JOB_NAME.out.{1} $JULIA sim6.jl {1} ${1} ${2} $SCRATCH/diffTest/sim6 ::: {1..1000} +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' + +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim6.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim6/sim6.jl b/exper/sim6/sim6.jl index dc84145..7062e25 100644 --- a/exper/sim6/sim6.jl +++ b/exper/sim6/sim6.jl @@ -1,9 +1,13 @@ +# Coverage of confidence intervals for one entry of Δ, with banded Ωx and Ωy that +# differ in the second off-diagonal. +# +# julia --project=.. sim6.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -56,8 +60,8 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(1234 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -77,4 +81,4 @@ ri, ci = elemArr[iElem] @time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim6/start_sbatch.sh b/exper/sim6/start_sbatch.sh old mode 100644 new mode 100755 index d64b637..9cf6dbe --- a/exper/sim6/start_sbatch.sh +++ b/exper/sim6/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" "4" "5" + for iElem in 1 2 3 4 5 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim6_${ip}_${iElem} sbatch_sim6 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim6_${ip}_${iElem} sbatch_sim6 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim6/summarize.jl b/exper/sim6/summarize.jl index d3549c5..3123afc 100644 --- a/exper/sim6/summarize.jl +++ b/exper/sim6/summarize.jl @@ -1,110 +1,57 @@ -using JLD -using DiffPrecTest -using DataFrames, DataFramesMeta, CSV -using Random, LinearAlgebra, Distributions -using LaTeXStrings, LaTeXTabulars +# Coverage, length and bias of the confidence intervals in sim6, from the files +# results/sim6_res_ip_iElem.jld2 written by collect.sh. Writes sim6.csv and +# sim6_table.tex. +# +# julia --project=.. summarize.jl +using CSV +using LinearAlgebra, Random +using Distributions -cd("exper/sim6") +include("../summary.jl") + +cd(@__DIR__) pArr = [100, 200] elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) -trueValue = zeros(Float64, 2, 5) +# the true difference, generated as in sim6.jl +function trueDelta(p) + Random.seed!(7658) + Ωx = Matrix{Float64}(I, p, p) + for l=1:p-1 + Ωx[l , l+1] = 0.3 + Ωx[l+1, l ] = 0.3 + end + for l=1:p-2 + Ωx[l , l+2] = 0.2 + Ωx[l+2, l ] = 0.2 + end + Ωy = Matrix{Float64}(I, p, p) + for l=1:p-1 + Ωy[l , l+1] = 0.3 + Ωy[l+1, l ] = 0.3 + end + for l=1:p-2 + Ωy[l , l+2] = -0.1 + Ωy[l+2, l ] = -0.1 + end + Δ = Ωx - Ωy + d = Vector{Float64}(undef, p) + rand!(Uniform(0.5, 2.5), d) + d .= sqrt.(d) + D = Diagonal(d) + D * Δ * D +end +df = newSummary() +trueValue = zeros(Float64, 2, length(elemArr)) for ip=1:2 - for iElem=1:5 - p = pArr[ip] - Random.seed!(7658) - - # generate model - Ωx = Matrix{Float64}(I, p, p) - for l=1:p-1 - Ωx[l , l+1] = 0.3 - Ωx[l+1, l ] = 0.3 - end - for l=1:p-2 - Ωx[l , l+2] = 0.2 - Ωx[l+2, l ] = 0.2 - end - Ωy = Matrix{Float64}(I, p, p) - for l=1:p-1 - Ωy[l , l+1] = 0.3 - Ωy[l+1, l ] = 0.3 - end - for l=1:p-2 - Ωy[l , l+2] = -0.1 - Ωy[l+2, l ] = -0.1 - end - Δ = Ωx - Ωy - # generate Delta - d = Vector{Float64}(undef, p) - rand!(Uniform(0.5, 2.5), d) - d .= sqrt.(d) - D = Diagonal(d) - Σx = inv(Symmetric(D * Ωx * D)) - Σy = inv(Symmetric(D * Ωy * D)) - - tΔ = D * Δ * D - - res = load("results/sim6_res_$(ip)_$(iElem).jld", "results") - - ri, ci = elemArr[iElem] - indE = (ci - 1) * p + ri - - trueValue[ip, iElem] = tΔ[ri, ci] - - for j in 1:9 - sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE]) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end + tΔ = trueDelta(pArr[ip]) + for iElem in eachindex(elemArr) + trueValue[ip, iElem] = tΔ[elemArr[iElem]...] + summarizeSetting!(df, "results/sim6_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], trueValue[ip, iElem]) end end CSV.write("sim6.csv", df) - - - -### create a table - -latex_tabular("sim6_table.tex", - Tabular("lllccc"), - [ - Rule(:top), - ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"], - Rule(:mid), - vcat([ - hcat( - (a = fill("", 15, 1); - b = latexstring("p = $(pArr[ip])"); - a[1] = "\\multirow{15}{*}{$(b)}"; - a - ), - vcat([ - hcat( - (a = fill("", 3, 1); - b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = $(trueValue[ip, iElem])"); - a[1] = "\\multirow{3}{*}{$(b)}"; - a - ), - Matrix{Any}( - @linq df |> - where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"), - :p .== pArr[ip], - :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |> select(:method, - coverage = convert.(Int, :coverage * 1000), - lenCI = round.(:lenCI, digits=3), - bias = round.(:bias * 1000, digits=1)) - ) - ) - for iElem=1:5 - ]... - ) - ) - for ip=1:2 - ]...), - Rule(:bottom) - ] -) +coverageTable("sim6_table.tex", df, pArr[1:2], elemArr, trueValue) diff --git a/exper/sim7/sbatch_sim7 b/exper/sim7/sbatch_sim7 index 2e60ae1..62312a0 100644 --- a/exper/sim7/sbatch_sim7 +++ b/exper/sim7/sbatch_sim7 @@ -1,22 +1,29 @@ #!/bin/bash -## call this file as -## sbatch --job-name=diff_sim1_${1}_${2} sbatch_sim1 dir +## submitted by start_sbatch.sh as +## sbatch --job-name=diff_sim7_${ip}_${iElem} sbatch_sim7 ip iElem +## +## JULIA (default: julia) and SCRATCH (default: $HOME/scratch) can be set in the +## environment. Results go to $SCRATCH/diffTest/sim7, logs to +## $SCRATCH/output/sim7. Pass the partition for your cluster to sbatch +## (sbatch -p ...). #SBATCH --ntasks=50 -#SBATCH --partition=broadwl #SBATCH --time=1-12:00:00 -##SBATCH --exclusive -module load parallel -module load hdf5 +JULIA=${JULIA:-julia} +SCRATCH=${SCRATCH:-$HOME/scratch} +PROJECT=$SLURM_SUBMIT_DIR/.. +OUT=$SCRATCH/diffTest/sim7 +LOG=$SCRATCH/output/sim7 -SCRATCH_PATH=/scratch/midway2/mkolar -JULIA=/project2/mkolar/julia-1.5.0/bin/julia +mkdir -p $OUT +mkdir -p $LOG +command -v parallel > /dev/null || module load parallel -mkdir -p $SCRATCH_PATH/diffTest/sim7 -mkdir -p $SCRATCH_PATH/output/sim7 +# one BLAS thread per worker; install and precompile once before the workers start +export OPENBLAS_NUM_THREADS=1 +$JULIA --project=$PROJECT -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' - -parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$SCRATCH_PATH/output/sim7/$SLURM_JOB_NAME.err.{1} --output=$SCRATCH_PATH/output/sim7/$SLURM_JOB_NAME.out.{1} $JULIA sim7.jl {1} ${1} ${2} $SCRATCH_PATH/diffTest/sim7 ::: {1..1000} +parallel --delay .2 -j $SLURM_NTASKS --joblog __$SLURM_JOB_NAME.runtask.log --resume srun --exclusive -N1 -n1 --error=$LOG/$SLURM_JOB_NAME.err.{1} --output=$LOG/$SLURM_JOB_NAME.out.{1} $JULIA --project=$PROJECT sim7.jl {1} ${1} ${2} $OUT ::: {1..1000} diff --git a/exper/sim7/sim7.jl b/exper/sim7/sim7.jl index 70e1296..b6d519d 100644 --- a/exper/sim7/sim7.jl +++ b/exper/sim7/sim7.jl @@ -1,9 +1,13 @@ +# Coverage of confidence intervals for one entry of Δ, with banded Ωx and Ωy that +# differ in the diagonal and the second off-diagonal. +# +# julia --project=.. sim7.jl rep ip iElem dir +# +# saves the estimates of all methods to dir/res_ip_iElem_rep.jld2 using DiffPrecTest -using Statistics, StatsBase, LinearAlgebra -using SparseArrays -using ProximalBase, CoordinateDescent, CovSel -using Random, Distributions -using JLD +using LinearAlgebra, Random +using Distributions +using JLD2 @show gethostname() @@ -56,8 +60,8 @@ dist_Y = MvNormal(convert(Matrix, Σy)) # generate data Random.seed!(1234 + rep) -X = rand(dist_X, n)' -Y = rand(dist_Y, n)' +X = permutedims(rand(dist_X, n)) +Y = permutedims(rand(dist_Y, n)) Sx = Symmetric( X'X / n ) Sy = Symmetric( Y'Y / n ) @@ -77,4 +81,4 @@ ri, ci = elemArr[iElem] @time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) @time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle) -@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport +@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld2" est eΔ esuppΔ eω esuppω eSupport diff --git a/exper/sim7/start_sbatch.sh b/exper/sim7/start_sbatch.sh old mode 100644 new mode 100755 index b4184c1..30ee055 --- a/exper/sim7/start_sbatch.sh +++ b/exper/sim7/start_sbatch.sh @@ -1,11 +1,14 @@ #!/bin/bash -for ip in "1" "2" +## submits one job per setting; arguments are passed on to sbatch, e.g. +## ./start_sbatch.sh -p + +for ip in 1 2 do - for iElem in "1" "2" "3" "4" "5" + for iElem in 1 2 3 4 5 do echo "sbatch ${ip} ${iElem} ..." - sbatch --job-name=diff_sim7_${ip}_${iElem} sbatch_sim7 ${ip} ${iElem} + sbatch "$@" --job-name=diff_sim7_${ip}_${iElem} sbatch_sim7 ${ip} ${iElem} echo "... done" done done diff --git a/exper/sim7/summarize.jl b/exper/sim7/summarize.jl index d245f20..0568ee4 100644 --- a/exper/sim7/summarize.jl +++ b/exper/sim7/summarize.jl @@ -1,60 +1,57 @@ -using JLD -using DiffPrecTest -using DataFrames, CSV -using Random, LinearAlgebra, Distributions +# Coverage, length and bias of the confidence intervals in sim7, from the files +# results/sim7_res_ip_iElem.jld2 written by collect.sh. Writes sim7.csv and +# sim7_table.tex. +# +# julia --project=.. summarize.jl +using CSV +using LinearAlgebra, Random +using Distributions + +include("../summary.jl") + +cd(@__DIR__) pArr = [100, 200] elemArr = [(5,5), (8, 7), (50, 25), (22, 20), (32, 30)] -methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"] - -df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[]) +# the true difference, generated as in sim7.jl +function trueDelta(p) + Random.seed!(7689) + Ωx = Matrix{Float64}(I, p, p) + for l=1:p-1 + Ωx[l , l+1] = 0.3 + Ωx[l+1, l ] = 0.3 + end + for l=1:p-2 + Ωx[l , l+2] = 0.2 + Ωx[l+2, l ] = 0.2 + end + Ωy = Matrix{Float64}(I, p, p) * 1.1 + for l=1:p-1 + Ωy[l , l+1] = 0.3 + Ωy[l+1, l ] = 0.3 + end + for l=1:p-2 + Ωy[l , l+2] = -0.1 + Ωy[l+2, l ] = -0.1 + end + Δ = Ωx - Ωy + d = Vector{Float64}(undef, p) + rand!(Uniform(0.5, 2.5), d) + d .= sqrt.(d) + D = Diagonal(d) + D * Δ * D +end +df = newSummary() +trueValue = zeros(Float64, 2, length(elemArr)) for ip=1:2 - for iElem=1:5 - p = pArr[ip] - Random.seed!(7689) - - # generate model - Ωx = Matrix{Float64}(I, p, p) - for l=1:p-1 - Ωx[l , l+1] = 0.3 - Ωx[l+1, l ] = 0.3 - end - for l=1:p-2 - Ωx[l , l+2] = 0.2 - Ωx[l+2, l ] = 0.2 - end - Ωy = Matrix{Float64}(I, p, p) * 1.1 - for l=1:p-1 - Ωy[l , l+1] = 0.3 - Ωy[l+1, l ] = 0.3 - end - for l=1:p-2 - Ωy[l , l+2] = -0.1 - Ωy[l+2, l ] = -0.1 - end - Δ = Ωx - Ωy - # generate Delta - d = Vector{Float64}(undef, p) - rand!(Uniform(0.5, 2.5), d) - d .= sqrt.(d) - D = Diagonal(d) - Σx = inv(Symmetric(D * Ωx * D)) - Σy = inv(Symmetric(D * Ωy * D)) - - tΔ = D * Δ * D - - res = load("../sim7_res_$(ip)_$(iElem).jld", "results") - - ri, ci = elemArr[iElem] - indE = (ci - 1) * p + ri - - for j in 1:9 - sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE]) - push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage]) - end + tΔ = trueDelta(pArr[ip]) + for iElem in eachindex(elemArr) + trueValue[ip, iElem] = tΔ[elemArr[iElem]...] + summarizeSetting!(df, "results/sim7_res_$(ip)_$(iElem).jld2", pArr[ip], elemArr[iElem], trueValue[ip, iElem]) end end CSV.write("sim7.csv", df) +coverageTable("sim7_table.tex", df, pArr[1:2], elemArr, trueValue) diff --git a/exper/sim8/aggregate_results.jl b/exper/sim8/aggregate_results.jl index e3130f5..6c76ed6 100644 --- a/exper/sim8/aggregate_results.jl +++ b/exper/sim8/aggregate_results.jl @@ -1,27 +1,25 @@ -# jl1 aggregate_results.jl /scratch/midway2/mkolar/diffTest/sim8/res_1 our_res_1 - +# Collect the bootstrap estimates of one setting of sims 8-11 into one file: +# +# julia --project=.. aggregate_results.jl $SCRATCH/diffTest/sim8/res_1 our_res_1 +# +# reads resPath_rep.jld2 for rep = 1, ..., NUM_REP (default 1000) and writes +# outPrefix.jld2. Missing repetitions are skipped with a warning. using DiffPrecTest -using JLD +using JLD2 resPath = ARGS[1] outPrefix = ARGS[2] -NUM_REP = 1000 +NUM_REP = parse(Int, get(ENV, "NUM_REP", "1000")) -results = Vector{BootstrapEstimates}(undef, NUM_REP) +results = BootstrapEstimates[] for rep=1:NUM_REP - global results - fname = "$(resPath)_$(rep).jld" - - try - file = jldopen(fname, "r") - res = read(file, "boot_res") - close(file) - - results[rep] = deepcopy(res) - catch - @show fname + fname = "$(resPath)_$(rep).jld2" + if isfile(fname) + push!(results, load(fname, "boot_res")) + else + @warn "missing $fname" end end -@save "$(outPrefix).jld" results \ No newline at end of file +@save "$(outPrefix).jld2" results diff --git a/exper/sim8/collect_our.jl b/exper/sim8/collect_our.jl index b72289d..4b02deb 100644 --- a/exper/sim8/collect_our.jl +++ b/exper/sim8/collect_our.jl @@ -1,69 +1,48 @@ -using JLD +# Coverage of the simultaneous bootstrap confidence intervals in sims 8-11 +# (Δ = 0), with estimated and oracle supports, from the files written by the +# Slurm jobs to $SCRATCH/diffTest/