roualdes / bridgestan Goto Github PK

Use stan-dev/math#2743 when it becomes available, after stan-dev/math and stan-dev/cmdstan updates.

This should reduce a copy in BridgeStan's use of stan::math::gradient, e.g.

Also forking off from #90.

Right now print statements in Stan end up getting sent to std::cerr unconditionally. Two things about this

1. If we're going to hardcode anything, we should probably make this std:cout
2. Some languages like Python have ways to redirect ""stdout"". I put this in quotes because Python's sys.stdout and C/C++'s stdout are not necessarily the same, the former does not even need to point to a valid file descriptor but could be a class. Unifying Stan's print statements with this language-specific behavior is tricky, though there are known tricks like what @bob-carpenter links to here: #90 (comment)

This is in many ways trickier than error messaging, since it is much harder to have it be zero-cost when you don't need it (the various redirection tricks you can use in Python all have a cost even if you never print anything. One potential solution is to have it be opt-in behavior, where you could write code like:

import sys
import bridgestan
m = bridgestan.StanModel(...) # initialize model

sys.stdout = # some file we want to redirect to
m.log_density(...) # would print to C's stdout, NOT that file
with bridgestan.redirect_stdout():
    m.log_density(...) # in this scope we WOULD redirect

The Julia tests require STAN_THREADS=true since they test threaded calls to log_density_gradient. We should add a sentence in the doc about this, similar to what we did in the doc for the Rust interface, re PR 130.

The main goal of the refactor is to organize docs into dev-facing and user-facing, make it easy for users and devs to just look at the language they're using or developing, and to put all the documentation source into either top-level or language-specific directories as appropriate. A secondary goal is to collect all the docs into one place, which I think should be the rendered docs on the github.io site.

Lightweight .md pointers from GitHub

Keep the two standard .md files, which are lightweight pointers to the actual doc on the rendered bridgestan.github.io site. This makes it easier to do things like site-specific Google searches and also avoids having to decide for each piece of doc whether it should go in GitHub or in the rendered web pages.

• bridgestan/README.md

  • project overview
  • link to top-level user-facing doc
  • link to top-level dev facing doc

• bridgestan/CONTRIBUTING.md

  • link to top-level dev facing doc

Actual doc rendered for bridgestan.github.io

• bridgestan/docs

  • language-agnostic documentation
  • docs/users/ for user-facing language-agnostic documentation
  • docs/devs/ for dev-facing language-agnostic documentation

• bridgestan/docs/L

  • top-level doc for language L with pointers to dev and user doc
  • bridgestan/docs/L/users for user-facing doc for language L
  • bridgestan/docs/L/devs for dev-facing doc for language L
  • bridgestan/docs/L/api for the API spec for language L

Top-level make target

The doc will be built with a top-level makefile target to generate the bridgestan.github.io doc.

make html

The Julia interface, like the Python interface, allows for make_args and stanc_args to be passed from Julia (or Python) when constructing a Stan model, see line 34 of BridgeStan.jl.

I think it would be helpful to have these arguments more prominently displayed in the Julia doc. For instance, the arguments don't show up in the doc for the function StanModel, and only show up in the doc for compile_model much further down the page.

I'm raising an issue here because it's not obvious to me where in the Julia source code we should include this new documentation.

The two options, as I see it, are

• src/model.jl: right above the definition of mutable struct StanModel, which is where the beginning of the Julia interface doc pulls from. However, this is not where the additional arguments, make_args and stanc_args are defined. So it might be awkward to put documentation here for arguments to a function that are defined else where in the source code.
• src/BridgeStan.jl: this is where the additional arguments make_args and stanc_args are defined, but I don't think this function doc shows up in the Julia interface page.

From our review for the Journal of Open Source Software, openjournals/joss-reviews#5236, it was recommended that we add to our documentation details on how a user could link to a pre-existing Stan directory instead of cloning Stan and Stan Math as submodules from within the BridgeStan repository.

Hi @roualdes

Over the weekend Seth (@sethaxen) and I decided to include bridgestan as a Git subtree and Seth implemented this.

It works fine (and this include Github CI workflows!). But it also always results in:

julia> include("/Users/rob/.julia/dev/StanSample/example/bernoulli.jl");
┌ Warning: BridgeStan path was not set, compilation will not work until you call `set_bridgestan_path!()`
└ @ StanSample.BridgeStan ~/.julia/dev/StanSample/deps/data/bridgestan/julia/src/BridgeStan.jl:40
[ Info: /Users/rob/.julia/dev/StanSample/example/tmp/bernoulli.stan updated.
4000×1 DataFrame
  Row │ theta     
      │ Float64   
──────┼───────────
    1 │ 0.239221
    2 │ 0.21372

as we discuss here. We define the BRIDGESTAN environment later on as the path to bridgestan now can vary.

I wonder if there is a better solution for this.

A minor issue we bumped into is that cmdstan's makefile is called makefile, BridgeStan's makefile is called Makefile. On MacOS this is not a problem but on Github's CI it does not work it seems if capitalization is not correct.

Thanks for your work!

Rob

See #80 (comment) and discussion:

In both Python and Julia, a user could eventually end up with many downloaded copies due to version bumps over time. For Python, this is easy to detect as they will be in ~/.bridgestan/bridgestan-VERSION. In Julia, it is more difficult, since they will live under ~/.julia/artifacts as git-tree-hashes, not useful names.

When evaluating the density of a model containing a simplex parameter, bridgestan gives a seg fault:

import bridgestan as bs
import numpy as np
model = bs.StanModel.from_stan_file("../test_models/simplex/simplex.stan")
param_length = model.param_num()
theta = np.array([1/param_length] * param_length)
theta_unc = model.param_unconstrain(theta)
log_p = model.log_density(theta_unc)

Output:

[1]    12697 segmentation fault (core dumped)  /usr/bin/python

The same occurs with the R interface.

As a side note, the model in test_models/simplex/simplex.stan had a misspecified prior which I've fixed in #37

Hello,
I am reviewing this package for submission at JOSS. When I tried to run the example.R file from the bridgestan/R folder, I get the following error:

Error in .C("bs_model_construct_R", as.character(data), as.integer(seed),  : 
  "bs_model_construct_R" not available for .C() for package "bernoulli_model"

My guess is that it originates here

We should have a way from within BridgeStan to report which version of CmdStan and which version of Stan it's using.

For CmdStan, we can get it from the makefile in CmdStan from the home directory:

make print-CMDSTAN_VERSION

For Stan, it's available as a three C++ constants.

Here's a list of Julia interface items I'd like to tackle one day.

• adopt more idiomatic package name spelling, BridgeStan instead of Bridgestan
• export K means nothing since K is never defined
• its common Julia practice to use f! when there's an out argument; add log_density, log_density_gradient (no exclamation points) which assign to new memory (instead of using an out argument)
• can drop the ::Vector{Float64} in most use cases of zeros(D); Float64 by default
• not all defined functions are exported
• double check tests for completeness

Experimenting with BridgeStan v1.0 I have no problem on my Mac (M2 currently). But testing on Github CI produces the error ./bin/stanc: 1: Not: not found, hence my question.

StanSample.jl, in subdirectory deps/data/bridgestan, contains a full clone of the BridgeStan repo. (git clone --recurse-submodules https://github.com/roualdes/bridgestan.git).

On the Mac ./bin/stanc is available in debs/data/bridgestan, but not in the repo on Github. .gitignore in debs/data/bridgestan does include an entry bin/*.

I've also tried different ways of "installing" BridgeStan the normal way in Julia, but that runs into path problems, e.g. in example.jl:

julia> BS.set_bridgestan_path!("../")
ERROR: Makefile does not exist at path! Make sure it was installed correctly.
../
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] validate_stan_dir(path::String)
   @ BridgeStan ~/.julia/packages/BridgeStan/iZm5y/src/compile.jl:21
 [3] set_bridgestan_path!(path::String)
   @ BridgeStan ~/.julia/packages/BridgeStan/iZm5y/src/compile.jl:37
 [4] top-level scope
   @ REPL[9]:1

I haven't tried treating BridgeStan similar to cmdstan, i.e. installing it outside of Julia, but that is a major step back from how BridgeStan worked pre-v1.0.

I must be overlooking an option on how to include BridgeStan in another package like StanSample. Any hints would be appreciated.

This is the first part of the CI output (step 8).

[358](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:361)
[ Info: /tmp/jl_zj63SF/bernoulli.stan updated.
[359](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:362)
┌ Warning: BridgeStan compilation of model failed.
[360](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:363)
└ @ StanSample ~/work/StanSample.jl/StanSample.jl/src/stanmodel/SampleModel.jl:139
[361](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:364)
BridgeStan: Error During Test at /home/runner/work/StanSample.jl/StanSample.jl/test/runtests.jl:65
[362](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:365)
  Got exception outside of a @test
[363](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:366)
  LoadError: 
[364](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:367)
  Error when compiling SampleModel bernoulli:
[365](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:368)
    % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
[366](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:369)
                                   Dload  Upload   Total   Spent    Left  Speed
[367](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:370)
  
[368](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:371)
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0
[369](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:372)
100     9    0     9    0     0     37      0 --:--:-- --:--:-- --:--:--    37
[370](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:373)
  ./bin/stanc: 1: Not: not found
[371](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:374)
  make: *** [/home/runner/work/StanSample.jl/StanSample.jl/deps/data/bridgestan/Makefile:51: /tmp/jl_zj63SF/bernoulli.hpp] Error 127
[372](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:375)
  
[373](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:376)
  Stacktrace:
[374](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:377)
    [1] SampleModel(name::String, model::String, tmpdir::String)
[375](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:378)
      @ StanSample ~/work/StanSample.jl/StanSample.jl/src/stanmodel/SampleModel.jl:140
[376](https://github.com/StanJulia/StanSample.jl/actions/runs/3659203720/jobs/6184937979#step:8:379)`

I'm currently trying to access the AD hessian for a very large model I have.

I've been investigating BridgeStan's AD hessian to use in the large model since I was having trouble noticing any big speed-up from compiling the model with BRIDGESTAN_AD_HESSIAN=true (see #52).

I decided to try sketching out a simpler model that has a non-trivial gradient and benchmarking it to see if I could identify whether or not the compiler flag was doing anything. My benchmarking seems to suggest that the AD hessian calculated with BRIDGESTAN_AD_HESSIAN=true is much slower than the default.

Here's the Julia code:

import Pkg; Pkg.activate(@__DIR__)
using BridgeStan
using JSON3
using BenchmarkTools
using Random
using Distributions
using LinearAlgebra

model_path = "gaussian.stan"
data_path = "gaussian-data.json"

# Make a data simulator to change the sample size
K = 10
N = 100
Random.seed!(1)
mu = randn(K)
sigma = rand(K)
draws = map(x -> x, eachrow(Matrix(rand(MvNormal(mu, Diagonal(sigma .^ 2)), N)')))
data = Dict("N" => N, "K" => K, "X" => draws)
open(data_path, "w") do io
    JSON3.write(io, data)
end

# Compile the BridgeStan model
lib_path = BridgeStan.compile_model(
    model_path;
    stanc_args = String["--O1"],
    # make_args = String["BRIDGESTAN_AD_HESSIAN=true"],
)

# Add data to the model
smb = StanModel(lib_path, data_path)

# Find out the number of parameters
num_params = BridgeStan.param_unc_num(smb)

# Run the model
bm = @benchmark lp = BridgeStan.log_density_hessian(smb, zeros(num_params))
display(bm)

and the Stan model, stored as gaussian.stan:

data {
  int<lower=0> N;
  int<lower=0> K;
  matrix[N,K] X;
}
parameters {
  vector[K] mu;
  vector<lower=0>[K]Sigma;
  cholesky_factor_corr[K] Omega;
}
model {
  mu ~ normal(0, 1);
  Omega ~ lkj_corr_cholesky(1);
  Sigma ~ cauchy(0, 2.5);
  for(n in 1:N)
    X[n,:] ~ multi_normal_cholesky(rep_vector(0, K), diag_pre_multiply(Sigma, Omega));
}

The benchmarking is a little confusing, because I would have expected the AD hessian to be much faster than finite differences. Perhaps the log density is simply just very easy to evaluate, so that the finite difference hessian is quicker?

Anyway, here's the results from running it with BRIDGESTAN_AD_HESSIAN=true

BenchmarkTools.Trial: 33 samples with 1 evaluation.
 Range (min … max):  134.849 ms … 164.949 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     153.700 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   151.588 ms ±   8.789 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▃                ▃ ▃              █    ▃▃   ▃            
  ▇▁▇▁▇▁█▁▁▇▁▁▁▇▁▁▁▁▁▁▁▁▁█▇█▁▁▁▇▁▇▁▁▁▁▇▇▇▁█▁▁▁▁██▇▇▁█▁▇▁▇▇▇▁▁▁▇ ▁
  135 ms           Histogram: frequency by time          165 ms <

 Memory estimate: 34.33 KiB, allocs estimate: 6.

and without the flag, in the default mode:

BenchmarkTools.Trial: 194 samples with 1 evaluation.
 Range (min … max):  22.855 ms … 27.557 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     27.042 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   25.794 ms ±  1.804 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                                       █▃▂     
  ▃▃▆▅▆▅▅▄▂▂▃▁▂▁▁▂▁▂▂▂▁▂▁▁▂▂▁▂▁▁▁▁▂▁▁▂▂▁▁▁▁▁▁▃▁▁▁▁▁▁▃▃████▅▃▃ ▂
  22.9 ms         Histogram: frequency by time        27.5 ms <

 Memory estimate: 34.33 KiB, allocs estimate: 6.

Note that the case where BRIDGESTAN_AD_HESSIAN=true takes nearly six times longer on average. Is this expected or am I missing something?

Misc hardware notes:

• Pop OS 22.04
• cmdstan 2.31
• Julia 1.8

Hi, follow up on #68 , I successfully installed the bridge_stan on Google Colab. And set the Python version on Colab to 3.9.16 and used the following code to test indeed the version is correct.

import sys
print(sys.version)

3.9.16 (main, Dec 7 2022, 01:11:51) \n [GCC 9.4.0]

However, I still cannot run the example because I keep getting this error msg now:

ValueError: BridgeStan folder '/usr/local/lib/python3.9' does not contain file 'Makefile', please ensure it is built properly!
If you need to set a different location, call 'set_bridgestan_path()'

Below is the example python code:

import bridgestan as bs
import numpy as np
stan = "/content/bridgestan/test_models/bernoulli/bernoulli.stan"
data = "/content/bridgestan/test_models/bernoulli/bernoulli.data.json"
model = bs.StanModel.from_stan_file(stan, data)

I have also tried some recommendation in #68 , such as changing bridgestan to bridgestan_repo:

!pip install /content/bridgestan_repo/python

Here is the google colab file I'm running

Thank you for this very useful project. I am using BridgeStan in a forthcoming paper and would like to properly cite it. Is there a preferred style or reference to cite? If there is, GitHub supports special CITATION files (https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-citation-files) that would make it easy to document the preferred citation style.

Weirdly specific title, but that's how it appears to be.

Running the following in the folder bridgestan/julia/mwe outputs (-27.5, [0.0, 0.0, -3.0, -4.0, -5.0]) for me:

using Pkg;
Pkg.add(url="https://github.com/roualdes/bridgestan.git", subdir="julia")
using BridgeStan
set_bridgestan_path!(pwd() * "/../..")

model_code = "
parameters {
  vector[5] fvb;
}
model {
  fvb ~ normal(0, 1);
}
"
data_string = "{}"


open("model.stan", "w") do f
    write(f, model_code)
end

open("data.json", "w") do f
    write(f, data_string)
end

model = StanModel(stan_file="model.stan", data="data.json")
println(log_density_gradient(model, collect(1.:5)))

I'm trying to get bridgestan to work on Windows 10 Version 10.0.19045 Build 19045.
I've downloaded both the Pypi version (through pip install) and the git clone version.
I can compile the .stan code to a .so binary and .hpp file with make test_models/bernoulli/bernoulli.stan.
However the compiled binary will not load into Python with the class StanModel:

import bridgestan as bs
import numpy as np

stan = "../test_models/bernoulli/bernoulli.stan"
data = "../test_models/bernoulli/bernoulli.data.json"

model = bs.StanModel.from_stan_file(stan, data)

print(f"This model's name is {model.name()}.")
print(f"It has {model.param_num()} parameters.")

x = np.random.random(model.param_unc_num())
q = np.log(x / (1 - x))  # unconstrained scale
lp, grad = model.log_density_gradient(q, jacobian=False)
print(f"log_density and gradient of Bernoulli model: {(lp, grad)}")

Traceback (most recent call last):
  File "...\bridgestan_test\example.py", line 26, in <module>
    model = bs.StanModel.from_stan_file(stan, data)
  File "...\venv\lib\site-packages\bridgestan\model.py", line 244, in from_stan_file
    return cls(
  File "...\venv\lib\site-packages\bridgestan\model.py", line 71, in __init__
    self.stanlib = ctypes.CDLL(self.lib_path)
  File "...\Python310\lib\ctypes\__init__.py", line 374, in __init__
    self._handle = _dlopen(self._name, mode)
FileNotFoundError: Could not find module '...\.bridgestan\bridgestan-2.0.0\test_models\bernoulli\bernoulli_model.so' (or one of its dependencies). Try using the full path with constructor syntax.

I have tried to trace where it goes wrong and the following also doesn't seem to work:

from pathlib import Path
BS_FILEPATH = r"C:\Users\MyUsername\.bridgestan\bridgestan-2.0.0\python"
compiled = Path(BS_FILEPATH).parent / "test_models/bernoulli/bernoulli_model.so"
print(compiled.exists()) # Returns True

import ctypes
stanlib = ctypes.CDLL(str(compiled))

FileNotFoundError: Could not find module '....bridgestan\bridgestan-2.0.0\test_models\bernoulli\bernoulli_model.so' (or one of its dependencies). Try using the full path with constructor syntax.

I abbreviated my filepaths, but they seem correct nonetheless.

This results the same error. Any ideas on how to solve this problem?

While working on the rust backend we noticed problems that appear when stan modules are loaded and unloaded in separate threads.

On windows we see a deadlock under the following conditions:

Start three threads. Load stan module A in the first thread, while the other threads wait. Load module B in the second thread. Then, unload module A in the first thread and exit the thread. And last, try to load module C in the third thread. This will deadlock. From some debugging it seems that the second thread actually died somewhere while exiting.
This exact behavior seems to depend on the compiler or tbb version, because if we compile the stan model with clang or msvc instead of mingw and a recent TBB, we do not see the same issue (see #105 (comment))

This might be related to TBB, we found some issues that look at bit similar to this on the TBB bugtracker (#105 (comment)).

A reproducer using the rust backend can be found here: https://github.com/aseyboldt/bridgestan/blob/main/rust/tests/create_models.rs#L236

There also seem to be some sporadic issues involving library unloading on linux as well however. With compilation settings O=0 and enabled library unloading I observed a couple of sporadic segfaults when I ran the rust test. In a debugger it looked like the regular expression here was not correctly initialized. I don't know for sure that this is related to unloading the libraries, but I've never seen that issue without it. And I think that since c++11 initialization of static variables should actually be thread safe.

(cc @WardBrian Feel free to add or modify the summary as you like)

As I read our doc, Python and R detail creation of a StanRNG object from methods on their respective StanModel classes, e.g. StanModel.new_rng(seed) and StanModel$new_rng(seed). Neither detail direct calls to StanRNG.

Julia on the other hand details direct calls to StanRNG, but does not have a new_rng(sm::StanModel, seed) function.

I could really go either way on adding new_rng(sm::StanModel, seed) to Julia. On the one hand, it would make things more consistent with the other interfaces, but it would literally just be another name for what StanRNG already does, same arguments and all.

On the other hand, I think it makes sense to doc StanRNG constructors for R and Python.

What do you think?

Expose the version of the Hessian-vector product functional from stan-dev/math.

Here's the functional we need to expose. We want the first one that uses reverse-mode nested in forward mode):

https://github.com/stan-dev/math/blob/develop/stan/math/mix/functor/hessian_times_vector.hpp

We should use the first version in the file (reverse nested in forward), because it's linear. Doing it with finite diffs and multiply will be quadratic (as will be using the second version in the source file).

The motivation is the optimizer (conjugate gradient trust region method) used by the following paper, which uses trust-ncg in scipi.optimize.minimize:

https://arxiv.org/pdf/2304.05527.pdf

@rgiordan said that they used PyMC because it provided Hessian-vector products.

I just tried to build the API doc locally and found out that we need the following pip packages

• pydata-sphinx-theme
• nbsphinx

Easy enough to deal with; no complaints.

What do we think of adding some doc for developers?

Thanks for the package! I'm using BridgeStan in Julia using StanSample.jl. I construct a BridgeStan.StanModel in order to compute log_density_gradient. I then use a sampler implemented in Julia to perform sampling.

Currently when I want to sample multiple chains I need to sample them serially, since if I try run the sampler within a loop marked by Threads.@threads, the result is a segfault at log_density_gradient

Is there a way to use StanModel in a multithreading context?

Since the last release my CI tests for nutpie fail on windows. If my remote debugging is right, this is due to the changes introduced in #118. This adds some code that is executed in __init__, ie right at import time, where it calls os.add_dll_directory with contents of the ~/.bridgestan directory. But if this runs on a fresh install, where stan was not downloaded already, this fails with a FileNotFound error.

Some logs from the CI:

https://github.com/pymc-devs/nutpie/actions/runs/5419095619/jobs/9851843966#step:8:335

If my interpretation of this is correct, this could be fixed by delaying the call to _windows_path_setup until after bridgestan is downloaded, ie by setting a flag to False at import time, and checking that flag each time a model is compiled. (strictly speaking it would be enough to check when a shared lib is loaded, but if, say, someone would want to load the library using a different interface, that interface would then have to do that...)

When I try using BridgeStan I get the following error:

[ Info: Precompiling BridgeStan [c88b6f0a-829e-4b0b-94b7-f06ab5908f5a]
ERROR: LoadError: syntax: expected assignment after "const"
Stacktrace:
  [1] top-level scope
    @ ~/.julia/packages/BridgeStan/YcHCB/src/model.jl:22
  [2] include(mod::Module, _path::String)
    @ Base ./Base.jl:418
  [3] include(x::String)
    @ BridgeStan ~/.julia/packages/BridgeStan/YcHCB/src/BridgeStan.jl:1
  [4] top-level scope
    @ ~/.julia/packages/BridgeStan/YcHCB/src/BridgeStan.jl:24
  [5] include
    @ ./Base.jl:418 [inlined]
  [6] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
    @ Base ./loading.jl:1318
  [7] top-level scope
    @ none:1
  [8] eval
    @ ./boot.jl:373 [inlined]
  [9] eval(x::Expr)
    @ Base.MainInclude ./client.jl:453
 [10] top-level scope
    @ none:1

removing the const qualifiers in model.jl enables me to at least compile/include BridgeStan.

Are you interested in adding BridgeStan.jl to the general Julia registry? This would make it easier for downstream packages like StanSample.jl to make sure they maintain compatibility with BridgeStan, since they could add BridgeStan as a dependency with version bounds. e.g. recently BridgeStan.jl dropped support for Julia versions older than Julia v1.8, which caused StanSample to break when BridgeStan was available and included.

I think the main challenge here is making sure that BridgeStan.jl has access to the rest of the bridgestan repo. Registering the package directly from the BridgeStan subdirectory wouldn't package its parent directory with it.

But it should be fairly straightforward to support this by making the tarball associated with the most recent bridgestan GitHub release an Artifact. Upon BridgeStan.jl installation, the tarball is automatically downloaded and extracted and the path is made available to BridgeStan.jl. The tarball also contains all git submodules. This is, for example, how PosteriorDB.jl interacts with posteriordb.

When a major change is made to bridgestan, the process for updating BridgeStan.jl would be 1) make a GitHub release (generates the tarball) 2) update an Artifacts.toml file in BridgeStan.jl. 3) register a new BridgeStan.jl. Step 2 can be automated.

If you're interested in this, I could open a PR to set up the artifacts.

Initializing a Stan model once, then changing the Stan model in-place (i.e. editing the .stan file) and reinitializing the Stan model does not seem to update the loaded library.

A possible but messy workaround seems to be to give each .so file an additional random string, e.g. in compile.jl

    rstr = randstring(6)
    dummy_file = splitext(absolute_path)[1] * "_$(rstr).stan"
    cp(absolute_path, dummy_file)
    output_file = splitext(absolute_path)[1] * "_$(rstr)_model.so"

Hi @roualdes

Need to make still one minor change to StanSample.jl but I'm wondering if below function can be much simpler:

function sim(smb::BridgeStan.StanModel, x=LinRange(0.1, 0.9, 100))
        q = zeros(length(x))
        ld = zeros(length(x))
        g = zeros(length(x))
        for (i, p) in enumerate(x)
            q[i] = @. log(x[]i / (1 - x[i]))        # unconstrained scale
            lp, grad = BridgeStan.log_density_gradient(smb, [q[i]], jacobian = 0)
            ld[i] = lp[1]
            g[i] = grad[1]
        end
        return DataFrame(x=x, q=q, log_density=ld, gradient=g)
    end

It feels like I'm doing something wrong in the call to log_density_gradient.

I was hoping for (and tried variations on) something like:

function sim(smb::BridgeStan.StanModel, x=LinRange(0.1, 0.9, 100))
    q = @. log(x / (1 - x))        # unconstrained scale
    lp, grad = BridgeStan.log_density_gradient(smb, q, jacobian = 0) # or a broadcasted version
    return DataFrame(x=x, q=q, log_density=ld, gradient=g)
end

Hi, I successfully installed the bridge_stan on Google Colab. However, I cannot really run the example because I keep getting this error msg:

AttributeError                            Traceback (most recent call last)
[<ipython-input-12-af3367af358c>](https://localhost:8080/#) in <module>
      4 stan = "../test_models/bernoulli/bernoulli.stan"
      5 data = "../test_models/bernoulli/bernoulli.data.json"
----> 6 model = bs.StanModel.from_stan_file(stan, data)
      7 
      8 print(f"This model's name is {model.name()}.")

AttributeError: module 'bridgestan' has no attribute 'StanModel'

Below is the example python code:

import bridgestan as bs
import numpy as np

stan = "../test_models/bernoulli/bernoulli.stan"
data = "../test_models/bernoulli/bernoulli.data.json"
model = bs.StanModel.from_stan_file(stan, data)

print(f"This model's name is {model.name()}.")
print(f"It has {model.param_num()} parameters.")

x = np.random.random(model.param_unc_num())
q = np.log(x / (1 - x))  # unconstrained scale
lp, grad = model.log_density_gradient(q, jacobian=False)
print(f"log_density and gradient of Bernoulli model: {(lp, grad)}")

Forking this issue off some discussion in #90 and something @justindomke asked me about when using JAX:

Because of the presence of an RNG in the bs_model class, it is difficult to treat it as fully thread safe.
This is only there because it is required for write_array in the Stan model class, and we furthermore know it is only really used by said function if the generated quantities block calls _rng functions in the Stan code.

Thread-safe options for this which would be relatively small changes:

1. Take in an int seed to each call to param_constrain, which initializes a boost:ecuyer1988 RNG, uses it once, and throws it away. This can get kind of sketchy, see recent changes in Stan: stan-dev/stan#3168. BUT, this plays nicely with the JAX idea of vmap/pmap, etc.
2. Expose the ability to create an RNG object to the interface, and require that param_constrain is passed one of these if include_gq == true. Thread safety then comes from each thread having its own unique RNG object they pass in.

I think option 1 is the easiest to use, since it doesn't require you juggling this other object, but obviously has the potential to be less performant and requires some outside source of entropy to produce fresh seeds.

As above. Trying to compile a model with include directives fails for my (from python, but I don't think it should depend on the interface).

MWE:

import bridgestan as bs
import os
bs.set_bridgestan_path('bridgestan')
os.system("cd bridgestan; git log --name-status HEAD^..HEAD; cd ..")

stan_functions_path = "./functions.stan"
stan_model_path = "./model.stan"
data_path = "./data.json"

with open(stan_functions_path, 'w') as fd:
    fd.write("""
functions {
    int test(){return 1;}
}
""")

with open(stan_model_path, 'w') as fd:
    fd.write("""
#include functions.stan
data {

}
parameters {
    real theta;
}
model {
    theta ~ std_normal();
}
""")

with open(data_path, 'w') as fd:
    fd.write("""{}""")


model = bs.StanModel.from_stan_file(stan_model_path, data_path)

Fails with output

commit f62bf462175df7411893e45823926a3d04eb9263 (HEAD -> main, origin/main, origin/HEAD)
Author: Brian Ward <[email protected]>
Date:   Wed Nov 30 12:54:46 2022 -0500

    Fix broken link

M       docs/install.rst

[...]

--- Translating Stan model to C++ code ---
./bin/stanc  --o=/home/niko/github/mwes/bridgestan/model.hpp /home/niko/github/mwes/bridgestan/model.stan
Makefile:49: recipe for target '/home/niko/github/mwes/bridgestan/model.hpp' failed

stderr:
Syntax error in '/home/niko/github/mwes/bridgestan/model.stan', line 2, column 0, include error:
   -------------------------------------------------
     1:  
     2:  #include functions.stan
         ^
     3:  data {
     4:  
   -------------------------------------------------

Could not find include file functions.stan in specified include paths.
make: *** [/home/niko/github/mwes/bridgestan/model.hpp] Error 1

This would only affect users who directly include bridgestan.h in a project, but it is useful (and fairly standard, I believe) to prefix C function names with a common string to avoid namespace issues that can arise in C. So, our C API would have functions like bs_name, while Python/Julia/R etc would still just use name.

I'm happy to make these changes, the question is just which prefix to use:

• bs short, but commonly used for another meaning
• stan - still short, stan_log_density etc all read naturally, but this isn't an official Stan project
• bridgestan - most accurate and unambiguous, but very long

With #91 and #98 merged, we can start to think about when we want to press release on 2.0

Changes since 1.0.2:

• Fix: Const correctness in C api in #89
• Send Stan outputs to stdout rather than stderr in #94
• Allow Stan outputs to be redirected to higher-level things like Python's sys.stdout in #108
• Expose version numbers in C API in #101
• Breaking change: Propagate error messages through C API in #91
• Allow nullptr as data argument to constructor in #104
• Breaking change: Implement param_constrain in a thread-safe way in #98
• Updated Stan to 2.32.0 in (#109)

Further changes coming at some point (non-breaking, could be 2.1 if we don't want to wait):

• Rust interface #105 (timeline TBD)

I naively expected Sphinx to provide well means to doc other languages than Python. With a bit more research, I'm losing my optimism.

Brute forcing Julia doc into Sphinx looked atrocious.

There's a Julia domain, which looks like it could be tricky to get right in a GitHub Action. I'll try it out locally tomorrow and report back.

I have not found a solution for R.

Here's a hybrid idea. We doc each interface using the common, or at least a preferred, tool within that particular language to produce an HTML file. Include the HTML file in the static directory of the sphinx pages, and then provide hyperlinks to those pages.

Pros: it should work; all interfaces would have doc in a format for which its users would better recognize.

Cons: all interface docs would look different; the seamless "Next >" buttons on the bottom right of the pages would make less sense if we hyperlink out of the languages folder.

Any other ideas?

Some unrelated questions about the README file.

1. Should we list in the README Julia, Python, and R minimum version numbers? It says Python 3, but I feel like we mean Python 3.9+. Julia should be 1.8+. I'd guess R should be 4+, but that's a wild guess. Thoughts?

2. Do we want to advertise/announce/mention the C-example in the README?

3. The sentence "The second command brings in the latest version of Stan." doesn't refer to anything. Should the code chunk above this sentence instead be

cd <cmdstan-dir>
make stan-update
make <bridgestan-dir>/test_models/multi/multi

I was re running some scripts dealing with simplexes, and ran into this.

Basically, for a valid simplex constraint, if I try to use the sampled draws from a stan program in model.param_unconstrain, I end up with

RuntimeError: param_unconstrain() failed with exception: Exception: stan::math::simplex_free:  
Simplex variable is not a valid simplex. sum(Simplex variable) = 0.99999959, but should be 1 

which is of course not a bridgestan issue, and other stan interfaces have also had similar problems I think, anyway this is only because of sig_figs being 6 by default. I'd wager this to be a common issue as bridgestan grows. As for a fix, I don't know if that's easy to do, or even in scope for BridgeStan, but if there'd be a way of adding tolerance for the constraint checks that would be quite nice. For example I have 400 chains with 1000 draws each for 1000 dimensions, so in this case I'd avoid more significant digits solely for the sake of convenience.

Anyway, this is not a bug, but if it's at all helpful for future scope, or if there's an easier fix other than changing sig_figs in Cmdstan itself I'd be very happy to hear about that.

I noticed this is in CmdStanR future milestone list but it sort of bleeds through all Stan Interfaces.

While working on rust langue support I found a few more (somewhat nitpicky) issues with the C API:

File paths

bs_construct takes a file path as one input (as a const *char). It then passes that to std::ifstream. I spend quite some time yesterday trying to figure out how different operating systems handle encodings for filenames and paths, and it seems to be quite a mess. I'm really not sure how I should pass in the path on different operating systems so that it works with non ASCII paths.
Is that really worth the effort to figure that out? I think it would be easier to specify that the input should always be a utf8 encoded json c-string. (Or change it to const *char + strlen?).
Language bindings can then use their languages mechanism to handle file paths, so the C API doesn't have to worry about that at all.

It is also a bit unclear right now what I should do if there is no input data. I think I have to pass in an empty string, but I think a null pointer would be cleaner?

Access to transformed data

I think there is currently no way to access quantities from the transformed data block? It would be nice if there was a way to extract those, so that they could for instance end up in the data group of an InferenceData object.

Integer generated values

If I understand the type system of stan correctly, it has primitive types float64 and int32.
In all functions of the C API that handle parameters, everything is a float64. I think that should be lossless, because every i32 value can be exactly represented as one float64 value.
I am not clear how complex values are encoded in the variable names though, and it would be great to have a way to figure out the actual types of the variables, so that I can convert them later in the traces.

Suggested on the forums here.

I think it is reasonable to accept a string which is either

1. A file path, if it ends with .json (can use stan::io::ends_with), or
2. A JSON object

This can avoid the IO read/write if the JSON file is being created by the calling script anyway

Continuing the conversation from #18 about when we think BridgeStan is ready for an official release and a version number.

Quoting @WardBrian.

After this PR I personally consider this repo to be in a "releasable" state, at least for a 0.1 release if not a 1.0

Seconded.

There are a few things remaining that I have on the list of possible features going forward, such as improved error messages (optional), and utilities in the interfaces for running compilation for you (would be very useful)

The latter would be a great feature, and improved error messages wouldn't hurt. Both though are features, in my mind, beyond the current usable state of BridgeStan.

Hence, I'm in favor of a 0.1 or 1.0 version number.

Any other thoughts, opinions, concerns?

From our review for the Journal of Open Source Software, openjournals/joss-reviews#5236, it was recommended that we address more clearly our software dependencies on Stan and Stan Math, including mentioning these as included git submodules.

From our review for the Journal of Open Source Software, openjournals/joss-reviews#5236, it was recommended that we suggest to users to clone the BridgeStan repository with --shallow-submodules since we don't necessarily need the git histories associated with Stan and Stan Math and this will save much disk space.

On my machine cloning BridgeStan using --shallow-submodules reports

$ du -hs bridgestan 
698M	bridgestan

and without --shallow-submodules

$ du -hs bridgestan 
1.2G	bridgestan

We have this ability in Julia, Python, and R, but not Rust.

Plausibly many users would want to evaluate gradients or hessians at more than one point. At that point being able to give a multidimensional array as theta_unc might be very useful. internally the iteration would happen irregardless so I'm not sure if it helps a lot computationally, but definitely very convenient. As in theta_unc being (N,D) where D is the number of unconstrained parameters. For example if I want to evaluate hessians for all the 1000 draws from one chain of a run, I would do something like this now

for idx, row in enumerate(idata):
        theta = model.param_unconstrain(row)
        model.log_density_hessian(theta, out_grad=grad[idx], out_hess=hessian[idx])

Instead if theta_unc could be a multidimensional array itself I could essential provide the entire sample data for all the unconstrained parameters and then get a 3-dimensional array containing the hessians. Similarly I could get a nice matrix of the unconstrained parameters without a for loop as well, if theta also supported the same structure.

Originally discussed in #21:
We currently have a bunch of folders named after languages (R, julia, python, c-example, stan) and one named src.

@bob-carpenter proposed renaming src to cpp.
I'm of the opinion that src is a special case, since all the other languages need the code contained there to function.

I do on the other hand think that stan sticks out, and we might want to rename it to something like test_models (since that's what they are). If we do that. then we have the convention that src is for the BridgeStan code, and any folder named after a language is for accessing that code from the language it is named after

There are a bunch of autodiff functors that are implemented in Stan but not exposed yet in BridgeStan. The two most basic are already done. Most of them other than directional derivatives require forward-mode autodiff. Please feel free to add more requests to the list.

• gradient
• Hessian
• 3rd order derivatives (called grad_Hessian in stan/math/mix; requires forward mode)
• directional derivatives (called gradient_dot_vector in stan/math/mix; most efficient in forward, can do backward)
• Hessian-vector product (called hessian_times_vector in stan/math/mix; requires forward mode)
• gradient of the trace of Hessian times a matrix, which is required for softabs Riemannian (grad_tr_mat_times_hessian in stan/math/mix; requires forward mode)

There's no inverse Hessian vector product in Stan. I'm not sure the best way to implement that---I think there are a lot of approaches because the direct way is so expensive (vector / Hessian in Stan).

It would be nice to be able to create downstream "user-facing" packages in python that use bridgestan, but allow people to install those packages via a simple pip install mypackage (and have pip dependency management install bridgestan as needed). The challenge with that is that as currently implemented, there are some manual steps to install bridgestan. Would it be possible to have bridgestan try to install stan itself (like pystan does) to avoid those manual steps? (Surely this would require some kind of compiler to exist on the user's system?) Possibly this could be done as some kind of option like pip install bridgestan[batteriesincluded].

There was some discussion about error handling in #88, but I thought a separate issue would be a better place to discuss this.

Currently, the error messages of the C++ exception are printed to stderr, which works fine, but has the disadvantage that many users might not know to look there, and depending on the setup it might not even be easy to find it (for example in a jupyter notebook that is started from jupyterhub or from a system service). Also, if several models are failing at the same time, a user would have to figure out which error belongs to which model manually.

The two (I think) most important kinds of failures would be

1. Errors during model creation (ie bs_construct). Users will I think provide invalid data fairly often, and the error message should contain information about what part of the data was invalid.
2. Errors during bs_log_density_gradient (or for different sampling / VI algorithms the other log_density functions): Error messages often contain information about why there was a divergence during sampling, and sampling algorithms might want to collect this information to help model debugging.

I can see a couple of different options of how this could be achieved, my favorites would probably be those:

1. Store the error information in the bs_model_rng struct as an extra field, and provide an accessor function for it. This unfortunately makes the bs_log_density_gradient function family no longer thread safe. We can work around that issue however by creating new bs_model_rng structs for each thread. With the current API this would require loading the data several times however. This again can be avoided by splitting bs_model_rng into two separate types: One for the dataset and the potential error message from model creating, and a separate context for rng and error messages during sampling.
   This would then look like this:

bs_model *model = bs_model_construct(data)
if (model == 0) {
    // A failure here indicates that we couldn't even allocate the error msg,
    // which would usually not happen...
    return -1;
}

// Callers have to check that there wasn't an error explicitly again!
const char *msg = bs_model_error_msg(model);
if (msg != 0) {
    // Data was invalid...
    // Do something with message
    bs_model_free(model);
}

// In each thread do the following:
{
    bs_ctx *ctx = bs_ctx_construct(model, seed, seed_idx);
    if (ctx == 0)
        stop();
    }

    // Call log_density and other functions in the sampler or...
    for (int i = 0; i < 10; i++) {
        int rc = bs_log_density(ctx, ...);
        if (rc != 0) {
            char *msg = bs_ctx_error_msg(ctx);
            // do something with message
        }
    }

    bs_ctx_free(ctx);  // Frees the msg field if necessary
}

// And finally a single call to free the memory in the
// original thread
bs_model_free(model);  # Frees the msg field if necessary

The additional function declarations:

// Create a new model. A non-null return value still has to
// be checked for errors with `bs_model_error_msg`.
bs_model *bs_model_construct(const char*);

// Check for an error in model creation.
const char *bs_model_error_msg(const bs_model*);

// Create a new ctx object with rng and space for an error message.
// the model has to outlive the ctx.
bs_ctx *ctx bs_ctx_construct(const bs_model*, int, int);

// Most other functions now take a bs_ctx argument.
// The ctx arguments are not const, because we might change the error.
int bs_log_density_gradient(bs_ctx*, ...);

const char *bs_ctx_error_msg(const bs_ctx*);

In this scenario bs_model would be Sync (functions that take a const pointer to bs_model can be called concurrently from different threads), but not Send (ownership of bs_model can not be transferred to a different thread, ie it has to be deallocated from the same thread where it was created. Is that actually a requirement?). bs_ctx doesn't need to be Sync or Send, as every thread can just create a separate context (which should be super cheap).

2. Add char **error_msg arguments to the fallible functions, and add a bs_free_error_msg(char *) function.
   This would then look something like this:

char *msg;
bs_model_rng *ctx = bs_construct(data, *msg);
if (ctx == 0) {
    do_something_with_msg(msg);
    bs_free_error_msg(msg);
}

int rc = bs_log_density(ctx, *msg, ...);
if (rc != 0) {
    do_something_with_msg(msg);
    bs_free_error_msg(msg);
}

I hope this is at least somewhat useful to you, happy to discuss more if there is interest in a change like this.

It seems like the greater BridgeStan community could benefit from setting, and using consistently within our doc, specific nouns to distinguish between BridgeStan users and callers of the C level API.

As I think of it, users refers to people who touch the interface level API. Further, in my opinion, we don't have a relatively settled noun for those who touch exclusively the C level API, distinct from the interface API.

Do you all agree such a distinction would help? If so, what nouns do you suggest?

When I run the R example I got this error

Error in initialize(...) : unused argument (0)

From the function help I think the 0 is an extra argument here

From our review for the Journal of Open Source Software, openjournals/joss-reviews#5236, it was recommended that we add details discussing the similarities and differences across the interfaces.