Better demos?

docs/make.jl

@@ -72,7 +72,10 @@ makedocs(;
     ],
 )
 
-rm(joinpath(script_dir, "build", "results"), recursive=true) # delete `results` folder before deploying
+try
+    rm(joinpath(script_dir, "build", "results"), recursive=true) # delete `results` folder before deploying
+catch 
+end
 
 deploydocs(;
     repo="github.com/Julia-Tempering/Pigeons.jl",

examples/stan/mRNA.stan

@@ -0,0 +1,55 @@
+functions {
+    // more accurate computation of exp(a)-exp(b), inspired by LSE
+    // uses the fact that expm1(x) does not underflow for x<<<0
+    // (see e.g. https://en.cppreference.com/w/c/numeric/math/expm1),
+    // whereas exp does (https://en.cppreference.com/w/c/numeric/math/exp)
+    // main identity
+    //     exp(a) - exp(b) = exp(max(a,b))[exp(a-max(a,b))-exp(b-max(a,b))] 
+    //     = 1{a>b}exp(a)[1-exp(b-a)] + 1{a<b}exp(b)[exp(a-b)-1]
+    //     = -1{a>b}exp(a)expm1(b-a) + 1{a<b}exp(b)expm1(a-b)
+    real exp_a_minus_exp_b(real a, real b){
+        return a > b ? -exp(a)*expm1(b-a) : exp(b)*expm1(a-b);
+    }
+
+    // compute mean for the Likelihood
+    //     mu(t) = [km0/(delta-beta)][exp(-beta(t-t0)) - exp(-delta(t-t0))]
+    //           = [km0/(delta-beta)] * exp_a_minus_exp_b(-beta(t-t0), -delta(t-t0))
+    // if delta ~ beta,
+    //     exp(-beta(t-t0)) - exp(-delta(t-t0)) ~ (t-t0)(delta-beta)
+    // so
+    //     mu(t) = km0(t-t0)
+    real get_mu(real tmt0, real km0, real beta, real delta){
+        real dmb = delta-beta;
+        return km0 * ( abs(dmb) < machine_precision() ? tmt0 : exp_a_minus_exp_b(-beta*tmt0, -delta*tmt0)/dmb );
+    }
+}
+data {
+    int <lower=0> N; // number of observations
+    array[N] real<lower=0> ts; // time of the observation
+    array[N] real ys; // observed value
+}
+parameters {
+    real<lower=-2,upper=1> lt0;
+    real<lower=-5,upper=5> lkm0;
+    real<lower=-5,upper=5> lbeta;
+    real<lower=-5,upper=5> ldelta;
+    real<lower=-2,upper=2> lsigma;
+}
+transformed parameters{
+    // real t0, km0, beta, delta, sigma;
+    real t0    = pow(10, lt0);
+    real km0   = pow(10, lkm0);
+    real beta  = pow(10, lbeta);
+    real delta = pow(10, ldelta);
+    real sigma = pow(10, lsigma);
+}
+model {
+    // Priors are all uniform so they are implicit
+    // Likelihood:
+    //     y_i|params ~indep N(mu_i, sigma)
+    // with
+    //     mu_i = get_mu(tmt0, km0, beta, delta, sigma)
+    for (i in 1:N) {
+        ys[i] ~ normal(get_mu(ts[i] - t0, km0, beta, delta), sigma);
+    }
+}

ext/PigeonsBridgeStanExt/toy_examples.jl

@@ -33,8 +33,21 @@ Pigeons.stan_banana(dim = 9, scale=1.0) =
         Pigeons.json(; dim, scale)
     )
 
+function Pigeons.stan_mRNA_post_prior_pair()
 
-observed_range_squared(x) = (maximum(x) - minimum(x))^2
+    ts = [0.0, 0.1960784314, 0.3921568627, 0.5882352941, 0.7843137255, 0.9803921569, 1.1764705882, 1.3725490196, 1.568627451, 1.7647058824, 1.9607843137, 2.1568627451, 2.3529411765, 2.5490196078, 2.7450980392, 2.9411764706, 3.137254902, 3.3333333333, 3.5294117647, 3.7254901961, 3.9215686275, 4.1176470588, 4.3137254902, 4.5098039216, 4.7058823529, 4.9019607843, 5.0980392157, 5.2941176471, 5.4901960784, 5.6862745098, 5.8823529412, 6.0784313725, 6.2745098039, 6.4705882353, 6.6666666667, 6.862745098, 7.0588235294, 7.2549019608, 7.4509803922, 7.6470588235, 7.8431372549, 8.0392156863, 8.2352941176, 8.431372549, 8.6274509804, 8.8235294118, 9.0196078431, 9.2156862745, 9.4117647059, 9.6078431373, 9.8039215686, 10.0]
+    ys = [-28.7432056846, -23.4516614434, -18.5388766424, -14.8222586806, -11.5382220299, -8.7488350492, -6.3148674646, -4.4116646219, -2.662417487, -1.2595749967, -0.2138655106, 0.7533834515, 1.4456718098, 2.256403532, 2.5768037788, 3.1629072506, 3.4035924126, 3.620654275, 3.6067766064, 3.9944954682, 3.8681716905, 3.9908347507, 3.8867758383, 3.9279954787, 3.9740235009, 3.8318783863, 3.8858488538, 3.6889274942, 3.667765745, 3.5132268101, 3.5563189655, 3.3140202907, 3.1375782101, 3.1215096119, 3.1411191324, 2.8587949725, 3.051174779, 2.7046377879, 2.8663809275, 2.7185400221, 2.4358115411, 2.3923558532, 2.4510098123, 2.3022021852, 2.213355665, 2.0576801516, 1.895682169, 1.8026728933, 1.9212597875, 1.8111628461, 1.8321635329, 1.754736356]
+    N = length(ts)
+
+    model = stan_example_path("mRNA.stan") 
+    data = Pigeons.json(; ts, ys, N)
+    empty_data = Pigeons.json(; ts = [], ys = [], N = 0)
+
+    return (;
+        posterior = StanLogPotential(model, data),
+        prior     = StanLogPotential(model, empty_data)
+    )
+end
 
 # the centered one is the "harder" one, see https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html
 function Pigeons.stan_eight_schools(centered = true)

src/pt/report.jl

@@ -1,7 +1,7 @@
 """ 
 $SIGNATURES
 
-The iterim diagnostics computed and printed to 
+The interim diagnostics computed and printed to 
 standard out at the end of every iteration 
 (this can be disabled using `show_report = false`).
 """

src/submission/MPI.jl

@@ -109,21 +109,22 @@ function mpi_submission_script(exec_folder, mpi_submission::MPI, julia_cmd)
     #    MethodError(f=Core.Compiler.widenconst, args=(Symbol("#342"),), world=0x0000000000001342)
     export JULIA_PKG_PRECOMPILE_AUTO=0
 
-    mpiexec $(mpi_submission.mpiexec_args) --merge-stderr-to-stdout --output-filename $exec_folder $julia_cmd_str
+    mpiexec $(cmd_to_string(mpi_submission.mpiexec_args)) --merge-stderr-to-stdout --output-filename $exec_folder $julia_cmd_str
     """
     script_path = "$exec_folder/.submission_script.sh"
     write(script_path, code)
     return script_path
 end
 
+cmd_to_string(cmd::Cmd) = "$cmd"[2:(end-1)]
 
 # Internal: "rosetta stone" of submission commands
 const _rosetta = (;
     queue_concept = [:submit,   :del,     :directive, :job_name,    :output_file,   :error_file,    :submit_dir,            :job_status,    :job_status_all,    :ncpu_info],
 
     # tested:
-    pbs           = [`qsub`,    `qdel`,   "#PBS",     "-N ",        "-o ",          "-e ",          "\$PBS_O_WORKDIR",      `qstat -x`,     `qstat -u`,         `pbsnodes -aSj -F dsv`],
-    slurm         = [`sbatch`,  `scancel`,"#SBATCH",  "--job-name=","-o ",          "-e ",          "\$SLURM_SUBMIT_DIR",   `squeue --job`, `squeue -u`,        `sinfo`],
+    pbs           = [`qsub`,    `qdel`,   "#PBS",     "-N ",        "-o ",          "-e ",          "\$PBS_O_WORKDIR",      `qstat -x`,     `qstat -u`,         `pbsnodes`],
+    slurm         = [`sbatch`,  `scancel`,"#SBATCH",  "--job-name=","-o ",          "-e ",          "\$SLURM_SUBMIT_DIR",   `squeue --job`, `squeue -u`,        `sinfo -o%C`],
 
     # not yet tested:
     lsf           = [`bsub`,    `bkill`,  "#BSUB",    "-J ",        "-o ",          "-e ",          "\$LSB_SUBCWD",         `bjobs`,        `bjobs -u`,         `bhosts`],

src/submission/presets.jl

@@ -12,9 +12,12 @@ UBC Sockeye cluster.
 """
 setup_mpi_sockeye(my_user_allocation_code) =
     setup_mpi(
-        submission_system = :pbs,
-        environment_modules = ["git", "gcc", "intel-mkl", "openmpi"],
-        add_to_submission = ["#PBS -A $my_user_allocation_code"], 
+        submission_system = :slurm,
+        environment_modules = ["gcc", "openmpi", "git"],
+        add_to_submission = [
+            "#SBATCH -A $my_user_allocation_code"
+            "#SBATCH --nodes=1-10000"  # required by cluster
+        ], 
         library_name = "/arc/software/spack-2023/opt/spack/linux-centos7-skylake_avx512/gcc-9.4.0/openmpi-4.1.1-d7o6cdvp67ngi5c5wdcw2qyjyseq3l3o/lib/libmpi"
     )
 

src/submission/submission_utils.jl

@@ -35,20 +35,9 @@ function queue_status()
 end
 
 function queue_ncpus_free()
-    mpi_settings = load_mpi_settings()
-    @assert mpi_settings.submission_system == :pbs "Feature only supported on PBS at the moment"
     r = rosetta()
-    n = 0
-    for line in readlines(`$(r.ncpu_info)`)
-        for item in eachsplit(line, "|")
-            m = match(r"ncpus[(]f[/]t[)][=]([0-9]+)[/].*", item)
-            if m !== nothing
-                suffix = m.captures[1]
-                n += parse(Int, suffix)
-            end
-        end
-    end
-    return n
+    run(`$(r.ncpu_info)`)
+    return nothing
 end
 
 """ 
@@ -59,7 +48,7 @@ Instruct the scheduler to cancel or kill a job.
 function kill_job(result::Result) 
     r = rosetta()
     exec_folder = result.exec_folder 
-    submission_code = readline("$exec_folder/info/submission_output.txt")
+    submission_code = queue_code(result)
     run(`$(r.del) $submission_code`)
     return nothing
 end

src/targets/StanLogPotential.jl

@@ -39,7 +39,7 @@ function stan_funnel end
 function stan_bernoulli end
 function stan_eight_schools end
 function stan_banana end
-
+function stan_mRNA_post_prior_pair end
 
 
 """
@@ -52,6 +52,8 @@ json(; variables...) =
     "{" *
     join(
         map(
-            pair -> "\"$(pair[1])\" : $(pair[2])",
+            pair -> "\"$(pair[1])\" : $(parse_json_item(pair[2]))",
             collect(variables)), ",") *
     "}"
+
+parse_json_item(item) = item == [] ? "[]" : "$item"

test/test_stan.jl

@@ -20,6 +20,11 @@ if !is_windows_in_CI()
 
         # some examples where an error is interpreted as -Inf:
         pigeons(target = Pigeons.stan_funnel(1), record = [online], n_chains = 1, n_rounds = 5, explorer = SliceSampler())
+
+        # check we get reasonable accept on one real example 
+        post_prior = Pigeons.stan_mRNA_post_prior_pair()
+        pt = pigeons(target = post_prior.posterior, reference = post_prior.prior, n_rounds = 5, n_chains = 2)
+        @test minimum(Pigeons.explorer_mh_prs(pt)) > 0.3
     end
 
     @testset "Stan restarts" begin