From 42003e35c5b11f185e98f212aceb864775d97ccb Mon Sep 17 00:00:00 2001
From: Pariksheet Nanda <pariksheet.nanda@uconn.edu>
Date: Fri, 10 May 2019 18:45:37 -0400
Subject: [PATCH] ENH: Work-in-progress Stan model for #1

---
 README.md                            | 10 +++++++
 examples/03-submit-param-sweep.slurm | 40 ++++++++++++++++++++++++++++
 examples/model.py                    | 27 +++++++++++++++++++
 examples/model_dep.stan              | 27 +++++++++++++++++++
 examples/model_indep.stan            | 16 +++++++++++
 5 files changed, 120 insertions(+)
 create mode 100644 examples/03-submit-param-sweep.slurm
 create mode 100644 examples/model.py
 create mode 100644 examples/model_dep.stan
 create mode 100644 examples/model_indep.stan

diff --git a/README.md b/README.md
index bb902eb..3bc3fd4 100644
--- a/README.md
+++ b/README.md
@@ -139,6 +139,10 @@ task workers complete records and become available.
 
 ```sh
 # From the command-line
+python2.7 -m pip install --user --upgrade pip
+python2.7 -m pip install --user --upgrade pystan
+module purge
+module load gcc/4.9.3		# Needed for pystan 2.19.0
 cd ~/parallel-slurm/examples
 rm -rf joblog submit.out results/ results.csv
 for i in {1..5}; do sbatch 03-submit-param-sweep.slurm; done
@@ -146,6 +150,12 @@ touch submit.out && tail -f submit.out
 # Hit Ctrl+C to exit
 ```
 
+Output:
+
+```
+
+```
+
 ## Next Steps
 
 Hopefully these examples have inspired you to use GNU Parallel to
diff --git a/examples/03-submit-param-sweep.slurm b/examples/03-submit-param-sweep.slurm
new file mode 100644
index 0000000..b32103d
--- /dev/null
+++ b/examples/03-submit-param-sweep.slurm
@@ -0,0 +1,40 @@
+#!/bin/bash
+#SBATCH --ntasks 5
+#SBATCH --output submit.out
+
+#SBATCH --dependency singleton
+#SBATCH --job-name example-03-parameter-sweep
+# Kill job after 5 minutes seconds to show resuming feature.
+#SBATCH --time 5:00
+
+parallel_opts=$(~/parallel-slurm/parallel_opts.sh)
+module load parallel
+
+# Alias our model running program.
+model="python model.py"
+
+# Read the number of simulations to run.
+n_sim=$( $model --sim )
+n_remaining=$( $model --remaining )
+
+# Check if all simulations are completed.
+echo "Started SLURM job $SLURM_JOB_ID"
+if [[ n_remaining -eq 0 ]]
+then
+    echo "Nothing to run; all $n_sim simulations complete."
+    echo "Completed SLURM job $SLURM_JOB_ID in $(sacct -nXj $SLURM_JOB_ID -o elapsed)"
+    exit 0
+fi
+
+# Run an interruption prone program.
+seq $n_sum > param_index
+echo "Running $n_remaining of total $n_sim simulations."
+# It's okay to use
+parallel $parallel_opts \
+	 --joblog joblog \
+	 --resume \
+	 --retries 3 \
+	 --line-buffer \
+	 $model \
+	 ::: param_index
+echo "Completed SLURM job $SLURM_JOB_ID in $(sacct -nXj $SLURM_JOB_ID -o elapsed)"
diff --git a/examples/model.py b/examples/model.py
new file mode 100644
index 0000000..1413db4
--- /dev/null
+++ b/examples/model.py
@@ -0,0 +1,27 @@
+# Example of Surgical institutional ranking from
+# http://www.openbugs.net/Examples/Surgical.html
+
+import os
+
+import pystan
+
+# This example considers mortality rates in 12 hospitals performing
+# cardiac surgery in babies. The data are shown below.
+#
+DATA = {
+    # Number of hospitals.
+    'N': 12,
+    # Number of operations.
+    'n': [ 47, 148, 119, 810, 211, 196, 148, 215, 207,  97, 256, 360],  # noqa: E201, E501
+    # Number of deaths.
+    'r': [  0,  18,   8,  46,   8,  13,   9,  31,  14,   8,  29,  24]   # noqa: E201, E501
+}
+
+if __name__ == "__main__":
+    # Disable threading.
+    os.environ['STAN_NUM_THREADS'] = "1"
+
+    fixed_effects = pystan.StanModel(file="model_indep.stan")
+    random_effects = pystan.StanModel(file="model_dep.stan")
+    burn_in = fixed_effects.sampling(data=DATA, warmup=1000, iter=0)
+    fit = random_effects.sampling(fit=burn_in, data=DATA, warmup=0, iter=10000)
diff --git a/examples/model_dep.stan b/examples/model_dep.stan
new file mode 100644
index 0000000..d953256
--- /dev/null
+++ b/examples/model_dep.stan
@@ -0,0 +1,27 @@
+data {
+  int N;
+  int<lower = 0> n[N];
+  int<lower = 0> r[N];
+}
+
+parameters {
+  real mu;
+  real tau;
+  real b[N];
+  real<lower = 0, upper = 1> p[N];
+  real pop_mean;
+  real sigma;
+}
+
+model {
+  for (i in 1:N) {
+    b[i] ~ normal_lpdf(mu, tau);
+    r[i] ~ binomial_lpmf(n[i], p[i]);
+    p[i] ~ logit(b[i]);
+  }
+
+  pop_mean = exp(mu) / (1 + exp(mu));
+  mu ~ normal(0.0, 1.0E-6);
+  sigma = 1 / sqrt(tau);
+  tau ~ gamma(0.001, 0.001);
+}
diff --git a/examples/model_indep.stan b/examples/model_indep.stan
new file mode 100644
index 0000000..3031202
--- /dev/null
+++ b/examples/model_indep.stan
@@ -0,0 +1,16 @@
+data {
+  int N;
+  int<lower = 0> n[N];
+  int<lower = 0> r[N];
+}
+
+parameters {
+  real<lower = 0, upper = 1> p[N];
+}
+
+model {
+  for (i in 1:N) {
+    p[i] ~ beta_lpdf(1.0, 1.0);
+    r[i] ~ binomial_lpmf(n[i], p[i]);
+  }
+}