Add @brm macro

.gitignore

@@ -27,3 +27,7 @@ Manifest*.toml
 # File generated by the Preferences package to store local preferences
 LocalPreferences.toml
 JuliaLocalPreferences.toml
+
+# Confidential client-project examples (local-only, never commit)
+web-macro/examples/bruno-*.jl
+web-macro/src/bruno-*.jl

Project.toml

@@ -24,5 +24,5 @@ FlexiChains = "0.3"
 LinearAlgebra = "1.10"
 PDMats = "0.11"
 Random = "1.10"
-Turing = "0.40"
+Turing = "0.40, 0.41, 0.42"
 julia = "1.10"

scripts/Benchmarking/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"

scripts/Benchmarking/main.jl

@@ -0,0 +1,90 @@
+# Gemini claims you can do this!
+using Pkg
+Pkg.activate(@__DIR__)
+insert!(LOAD_PATH, 2, joinpath(@__DIR__, ".."))
+
+include("../../web-macro/src/macro.jl")
+include("../../web-macro/src/vimpl.jl")
+include("../examples/database.jl")
+
+using Chairmarks, Random
+
+db = Database()
+df = db.dataset(:bambi, :escs)
+fdf = map(Vector{Float64}, (;df.drugs, df.o, df.c, df.e, df.a, df.n))
+# Using fdf is considerably faster than using df
+tmp = @brm fdf """
+    loc ~ o + c + e + a + n
+    log(err_scale) ~ 1
+    drugs ~ Normal(loc, err_scale)
+"""
+display(@brm """
+    loc ~ o + c + e + a + n
+    log(err_scale) ~ 1
+    drugs ~ Normal(loc, err_scale)
+""")
+vtmp = VBRMI(tmp)
+display(vtmp)
+display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity($vtmp, _))
+
+# import Reactant
+
+# rx = Reactant.to_rarray(randn(LogDensityProblems.dimension(vtmp)))
+# errors: NoFieldMatchError(...)
+# _rtmp = Reactant.to_rarray(vtmp)
+# errors: MethodError: no method matching _copyto!(::SubArray{…}, ::Base.Broadcast.Broadcasted{…})
+# rtmp = Reactant.@compile LogDensityProblems.logdensity(vtmp, rx)
+# display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity($rtmp, _))
+# error()
+
+using LogDensityProblemsAD, Mooncake, Enzyme, DifferentiationInterface
+struct ADLogDensity{F, B, E}
+    f::F
+    backend::B
+    extras::E
+end
+
+ADLogDensity(f, backend) = ADLogDensity(
+    f, 
+    backend, 
+    DifferentiationInterface.prepare_gradient(
+        Base.Fix1(LogDensityProblems.logdensity, f), 
+        backend, 
+        zeros(LogDensityProblems.dimension(f))
+    )
+)
+
+LogDensityProblems.capabilities(::ADLogDensity) = LogDensityProblems.LogDensityOrder{1}()
+LogDensityProblems.dimension(p::ADLogDensity) = LogDensityProblems.dimension(p.f)
+LogDensityProblems.logdensity(p::ADLogDensity, x) = LogDensityProblems.logdensity(p.f, x)
+LogDensityProblems.logdensity_and_gradient(
+    p::ADLogDensity, x
+) = DifferentiationInterface.value_and_gradient(
+    Base.Fix1(LogDensityProblems.logdensity, p.f), p.extras, p.backend, x
+)
+
+
+
+mvtmp1 = ADgradient(AutoMooncake(), vtmp)
+mvtmp2 = ADLogDensity(vtmp, AutoMooncake())
+evtmp1 = ADgradient(AutoEnzyme(; mode = Enzyme.set_runtime_activity(Enzyme.Reverse)), vtmp)
+evtmp2 = ADLogDensity(vtmp, AutoEnzyme(; 
+    mode=Enzyme.set_runtime_activity(Enzyme.Reverse),
+    function_annotation=Enzyme.Duplicated
+))
+
+x = randn(LogDensityProblems.dimension(vtmp))
+display(mapreduce(hcat, (mvtmp1, mvtmp2, evtmp1, evtmp2)) do b 
+    LogDensityProblems.logdensity_and_gradient(b, copy(x))[2]
+end)
+
+@info "Primal"
+display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity($vtmp, _))
+@info "Inefficient Mooncake skipped..."
+# display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity_and_gradient($mvtmp1, _))
+@info "Mooncake"
+display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity_and_gradient($mvtmp2, _))
+@warn "wrong Enzyme skipped..."
+# display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity_and_gradient($evtmp1, _))
+@info "correct Enzyme"
+display(@be randn(LogDensityProblems.dimension(vtmp)) LogDensityProblems.logdensity_and_gradient($evtmp2, _))

scripts/Project.toml

@@ -1,9 +1,16 @@
 [deps]
 BayesianRegressionModels = "cdd3e328-398d-47a6-a87b-6047aaf4b4bc"
-CatalogServer = "a1b2c3d4-0000-0000-0000-000000000001"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DynamicObjects = "23d02862-63fe-4c6e-8fdb-1d52cbbd39d5"
+ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
 FlexiChains = "4a37a8b9-6e57-4b92-8664-298d46e639f7"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
+LogDensityProblemsAD = "996a588d-648d-4e1f-a8f0-a84b347e47b1"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Oxygen = "df9a0d86-3283-4920-82dc-4555fc0d1d8b"
-YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

scripts/examples/database.jl

@@ -0,0 +1,12 @@
+include("all.jl")
+using DynamicObjects
+
+@dynamicstruct "" serial struct Database
+    mod(m::Symbol) = _SOURCE_MODULES[m]
+    @cached dataset(m::Symbol, name::Symbol) = getproperty(mod(m), :load)(Val(name)) 
+    # formula(m::Symbol, name::Symbol) = getproperty(mod(m), :examples)(Val(name))[1]
+    # example(m::Symbol, name::Symbol) = begin
+    #     _formula, datakey = getproperty(mod(m), :examples)(Val(name))
+    #     brm(dataset(m, datakey), _formula)
+    # end
+end 

web-macro/Project.toml

@@ -0,0 +1,35 @@
+name = "BRMMacroWeb"
+uuid = "cd55decf-6be8-4ea5-907a-b2dc35e4cc14"
+version = "0.1.0"
+
+[deps]
+AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67"
+AlgebraOfVega = "a2420894-485a-4f42-a4fc-c4323ec6fb2b"
+BridgeStan = "c88b6f0a-829e-4b0b-94b7-f06ab5908f5a"
+CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DynamicObjects = "23d02862-63fe-4c6e-8fdb-1d52cbbd39d5"
+ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+FlexiChains = "4a37a8b9-6e57-4b92-8664-298d46e639f7"
+HTMX = "27f3e1ef-6ef8-44dc-9e9a-2fb23ed44e83"
+HTMXObjects = "b12ef442-5798-4353-80f3-9562b03a0cb6"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StanBlocks = "2e771a56-c23a-4e0b-9282-20c2e37157e9"
+StanLogDensityProblems = "a545de4d-8dba-46db-9d34-4e41d3f07807"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+TestModules = "63c02187-99fd-4e5c-aaf0-4d6bfebc181c"
+Treebars = "e1e568c4-3a56-40a4-95fa-9b9c6c16fccb"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
+WarmupHMC = "60658175-6863-4866-a322-ab51a11c0cfe"

web-macro/app/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+AlgebraOfVega = "a2420894-485a-4f42-a4fc-c4323ec6fb2b"
+BRMMacroWeb = "cd55decf-6be8-4ea5-907a-b2dc35e4cc14"
+DynamicObjects = "23d02862-63fe-4c6e-8fdb-1d52cbbd39d5"
+HTMX = "27f3e1ef-6ef8-44dc-9e9a-2fb23ed44e83"
+HTMXObjects = "b12ef442-5798-4353-80f3-9562b03a0cb6"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+TestModules = "63c02187-99fd-4e5c-aaf0-4d6bfebc181c"
+Treebars = "e1e568c4-3a56-40a4-95fa-9b9c6c16fccb"
+WarmupHMC = "60658175-6863-4866-a322-ab51a11c0cfe"

web-macro/app/main.jl

@@ -0,0 +1,8 @@
+using Revise
+using BRMMacroWeb
+
+begin
+    BRMMacroWeb.terminate()
+    port = length(ARGS) >= 1 ? parse(Int, ARGS[1]) : 8121
+    BRMMacroWeb.serve(; host="0.0.0.0", revise=:lazy, port, async=true)
+end

web-macro/examples/1.1-verify-bernoulli-binomial-done.jl

@@ -0,0 +1,15 @@
+# label: 1.1 verify Bernoulli/Binomial
+# tier: 1
+# status: open
+#=
+**Status: done.** ✓ Confirmed that the existing `FBroadcasted{<:Type{<:Distribution}}` pass-through in `vimpl.jl` handles both Bernoulli and Binomial cleanly.
+
+**Verification.** The form below loads the **cbpp + therapeutic touch** model — a faithful translation of `brms::cbpp_binomial` (categorical predictor + random intercept + Binomial with per-row trial counts) and `kruschke::therapeutic_touch` (hierarchical Bernoulli) into one multi-likelihood model. brms's `incidence | trials(size) ~ ...` sidecar collapses to a plain positional argument: `bin_succ ~ Binomial(bin_n, logistic(η))`. If the gradient sanity check stays green every other `Distribution` family (Beta, Gamma, NegBinomial, …) should be free as well.
+
+=#
+
+log_odds_bin ~ 1 + c2 + (1 | g1)
+bin_succ ~ Binomial(bin_n, logistic(log_odds_bin))
+
+log_odds_b ~ 1 + (1 | g1)
+bin_y ~ Bernoulli(logistic(log_odds_b))

web-macro/examples/1.2-offset-fixed-exposure-already-works-without-a-wrapper.jl

@@ -0,0 +1,21 @@
+# label: 1.2 offset / fixed exposure — already works without a wrapper
+# tier: 1
+# status: open
+#=
+**Status: already works without any new code.** brms needs `offset(z)` because R's formula syntax has no other way to put a "no-coefficient term" into the linear predictor — the only thing on the RHS of `~` is the formula DSL. In our DSL the linear predictor and the likelihood are *separate* `~` lines, and the second one (the likelihood) takes a free-form Julia expression. Anything inside that expression gets evaluated as plain code at materialization time via `vbroadcasted` — function calls dispatch to whatever Julia function the symbol resolves to, and data column references are pulled from the dataframe.
+
+So instead of `count ~ x + offset(log(exposure))`, you write the offset directly inside the likelihood:
+
+```julia
+loc ~ 1 + a
+k1 ~ Poisson(exp(loc + log(exposure)))
+```
+
+The `log(exposure)` here is just `Base.log` applied to the `exposure` data column, broadcasted across rows and added to `loc` (which is the materialized linear predictor). No parameter is allocated for it because `growblock!!` is never called for that branch — there's no `~` on the data side, just an argument to `Poisson(...)`.
+
+**Verification.** Form below loads exactly that model. The VBRMI dim should match the offset-free version (only the population intercept + slope on `a`); the gradient sanity check should stay green; and the materialized `k1` likelihood should incorporate the row-specific exposure shift.
+
+=#
+
+loc ~ 1 + a
+k1 ~ Poisson(exp(loc + log(exposure)))

web-macro/examples/1.3-i-expr-likely-already-works.jl

@@ -0,0 +1,24 @@
+# label: 1.3 protect(expr) — brms literal-escape
+# tier: 1
+# status: done (vimpl)
+#=
+**Status: done in vimpl.** Our DSL already supports arbitrary Julia function
+calls on the RHS (e.g. `a^2`, `sqrt(abs(b))`, `log(exposure)`), so `protect()` is
+functionally unnecessary. But brms users who have internalized the `protect()`
+convention can now write it literally — we added a passthrough overload
+(`vbroadcasted(::ExprColumn{typeof(protect)}) = inner`) so `protect(expr)` behaves the same
+as `expr` alone, no parameter penalty.
+
+**Why brms needs protect().** In R, `a:b`, `a*b`, `|` etc. have DSL meaning inside a
+formula — `protect()` escapes them so the inner expression is interpreted as plain
+arithmetic. Our formula DSL is parsed by Julia first, so `a^2` is already a
+plain function call and `protect()` is a no-op.
+
+**Verification.** The form below mixes the "naked" spelling (`a^2`,
+`sqrt(abs(b))`, `log(exposure)`) with protect()-wrapped forms — the VBRMI dim,
+gradient check, and materialized predictor values must be identical.
+
+=#
+
+loc ~ 1 + a + protect(a^2) + protect(sqrt(abs(b))) + log(exposure)
+y1 ~ Normal(loc, 1)

web-macro/examples/1.4-scale-x-standardize-x.jl

@@ -0,0 +1,35 @@
+# label: 1.4 zscale(x) / center(x) / standardize(x)
+# tier: 1
+# status: done (vimpl)
+#=
+**Status: done in vimpl.** Three data-transform wrappers that z-transform (or
+just center) a column at VBRMI-materialization time. brms does `zscale(x)`
+automatically in most vignettes for sampler stability + prior scale-invariance;
+making it available as a formula-level wrapper avoids forcing users to
+pre-transform their DataFrame columns.
+
+**Semantics.**
+- `center(x)`      -> `x - mean(x)`
+- `zscale(x)`       -> `(x - mean(x)) / std(x)`
+- `standardize(x)` -> alias for `zscale(x)`
+
+Each fires inside `vbroadcasted`, so transforms compose with every downstream
+consumer: pop predictor, ranef slope, link functions, interactions, etc.
+`zscale(a):zscale(b)` works as expected (z-transformed both operands before the
+elementwise product).
+
+**Implementation.** Three-line `vbroadcasted` overloads in vimpl.jl using
+inlined `_mean`/`_std` helpers (no new dependency). The transform evaluates
+once at VBRMI construction and caches the resulting vector.
+
+**Verification.** Same model spelled two ways:
+- `loc ~ 1 + zscale(a) + zscale(b)` -- standardized at formula level
+- `loc ~ 1 + a_z + b_z` (where a_z, b_z were pre-standardized in the DataFrame)
+
+VBRMI dim + log-density should be identical; fitted coefficients ≈ the
+unscaled coefficients × `std(x)`.
+
+=#
+
+loc ~ 1 + zscale(a) + center(b)
+y1 ~ Normal(loc, 1)

web-macro/examples/1.5-zerocorr-independent-random-effects.jl

@@ -0,0 +1,33 @@
+# label: 1.5 zerocorr — `(terms || group)` for independent random effects
+# tier: 1
+# status: done (vimpl)
+#=
+**Status: done in vimpl.** brms (via lme4 syntax) lets you opt out of the LKJ
+correlation between multiple random terms on the same grouping factor:
+`(1 + x || g)` says "fit a random intercept and a random slope, but treat them
+as independent — don't estimate a 2x2 Cholesky between them". Useful when
+there's not enough data to estimate the correlation, or when you have prior
+reason to believe the terms are uncorrelated.
+
+**Semantics.** `(1 + a + b || g)` splits into three independent 1-column
+grouped_normal blocks, each keyed distinctly from plain `(... | g)` so the
+uncorrelated and correlated variants never accidentally coalesce. Per-term
+1x1 Cholesky reduces to a single log_scale parameter (same `chol` machinery
+used for the scalar `(1 | g)` case).
+
+**Implementation.** Our `_x` walker wraps `||` as `ExprColumn{typeof(doublepipe)}`.
+We added a `vmeta_sampling_rhs` overload that flattens the LHS on `+`, then
+foldl's over the terms with a synthetic per-term group key
+(`Symbol(name(rhs), :__nocor__, i)`) so each term lands in its own block.
+
+**Verification.** Preset below vs. correlated counterpart
+`loc ~ 1 + (1 + a | g1)`:
+- correlated variant: 1 + `chol(2)` (3 params) + `grouped_normal(8, 2)` (16 params) = 20 params
+- uncorrelated variant: 1 + 2 * (`chol(1)` (1 param) + `grouped_normal(8, 1)` (8 params)) = 19 params
+
+So the only difference is one fewer Cholesky parameter (the off-diagonal).
+
+=#
+
+loc ~ 1 + (1 + a || g1)
+y1 ~ Normal(loc, 1)

web-macro/examples/1.6-cache-levels-level_map-dense-gc_idx.jl

@@ -0,0 +1,17 @@
+# label: 1.6 cache levels / level_map / dense / gc_idx
+# tier: 1
+# status: open
+#=
+**What it is.** Stop rebuilding the dense level mapping (`Dict(level => row_index)`) and the gc_idx vector on every `VBRMI(brmi)` call. Cache them once per source data column.
+
+**Why it matters.** Today every `VBRMI` build re-traces the categorical / grouping columns, sorts unique values, builds a Dict, and walks the column to dense-encode it. For models with many categorical columns or many `VBRMI` rebuilds (e.g. during AD), this adds up.
+
+**Implementation.** Pick a storage layout for the per-column metadata. Two candidates:
+- A new `meta.factor` NamedTuple keyed by source column name, holding `(; levels, level_map, dense, gc_idx)` per column. Built lazily on first reference, indexed via `name(column)`.
+- Attach the metadata to `meta.materialized[column_name]` directly. More tightly coupled but avoids a parallel NamedTuple.
+
+The TODO already lives at `vimpl.jl:78–86`. Once a layout is picked, refactor `_gc_idx` and the inline dense map in the categorical path to read from the cache, falling back to a build-on-miss helper.
+
+**Verification.** No behavioral change — the gradient sanity check should stay green. Benchmark `VBRMI(brmi)` with Chairmarks before/after and confirm a measurable speedup on a model with multiple categorical/grouping columns.
+
+=#

web-macro/examples/1.7-categoricalarrays-pooledarrays-integration.jl

@@ -0,0 +1,19 @@
+# label: 1.7 CategoricalArrays / PooledArrays integration
+# tier: 1
+# status: open
+#=
+**What it is.** When the input column is already a `CategoricalVector` or a `PooledArray`, the dense level mapping is already computed and stored in the column's `.refs` field. Use it directly instead of rebuilding via `Dict`.
+
+**Why it matters.** Most real-world DataFrames use `CategoricalArrays.jl` for factor columns. Skipping the rebuild eliminates allocation entirely for the common case and gets us "for free" interop with the standard categorical-data ecosystem.
+
+**Implementation.** Two design choices:
+- **Hard dep**: add `CategoricalArrays` to vimpl.jl's deps, dispatch on `CategoricalVector`, read `levelcode.(col)` and `levels(col)` directly.
+- **Duck-typed**: sniff for the `.refs` field and `levels` method without importing the package, falling back to the generic Dict path.
+
+Recommend hard dep — it's the standard for tabular Julia code, and the duck-type path is more code with no real win. Same for `PooledArrays`.
+
+The actual integration is small once the design is picked: a method specialization in `_gc_idx` and in the categorical-predictor path. Composes with the caching TODO above.
+
+**Verification.** Preset (or test) that builds a DataFrame with a `CategoricalVector` column and uses it as a grouping factor / categorical predictor. Confirm the gradient sanity check stays green and the per-VBRMI allocation count drops.
+
+=#