AOT compilation

[sshin23]

Succeeds #252

[sshin23]

Summary of changes relative to origin/main

1. Core API redesign — functional, immutable model building

ExaCore is now an immutable struct with full type parameters {T,VT,B,S,V,P,O,C,R}. The old mutable struct API (variable(), parameter(), objective(), constraint(), constraint!(), subexpr()) is replaced by a functional API where every add_* call returns (new_core, handle):

c, x  = add_var(c, n; start=0, lvar=-Inf, uvar=Inf)
c, θ  = add_par(c, start_array)
c, o  = add_obj(c, expr for i in ...)
c, g  = add_con(c, expr for i in ...)
c, s  = add_expr(c, expr for i in ...)
c     = add_con!(c, g, contrib for i in ...)   # augment existing constraint

The macro layer (@var, @par, @obj, @con, @con!, @expr) wraps this and updates core in the calling scope.

Renamed types:

• Subexpr / ReducedSubexpr → Expression (always inline, no auxiliary variables)
• ParameterSubexpr → Parameter
• Variable{S,O} → Variable{D,S,O} (added disambiguation index D for variable block ordering)
• Objective{R,F,I} → Objective{F,I}, Constraint{R,F,I,O} → Constraint{F,I,O} (removed nested inner linked-list field)

Removed sentinel types: VariableNull, ParameterNull, ObjectiveNull, ConstraintNull.

Objectives and constraints are now stored as tuples instead of linked-list chains (dispatched via Tuple{} base case + first/Base.tail recursion). This gives fully concrete types at compile time — the key enabler for juliac.

---

2. AOT compilation via juliac --trim=safe

The main goal of this PR. Several changes were needed to make expression trees fully type-stable:

Const{T} <: AbstractNode (src/specialization.jl): a constant wrapper whose value lives in a field (not the type parameter), so Const{Float64} is always the same concrete type regardless of value. Unlike Val{V}, this does not trigger recompilation or produce abstract Val{<:Real} widening that blocks --trim=safe.

Auto-const wrapping in macros (_maybe_const, _auto_const_gen): free numeric symbols and literals inside @obj/@con/@expr generator bodies are automatically wrapped in Const(...), preventing abstract Val{<:Real} when runtime values enter expression trees.

exa_sum / exa_prod with SumNode{I} / ProdNode{I} and Val{N}-based specialization for type-stable sum/product nodes over unit ranges. The @obj/@con/@expr macros rewrite sum(body for j in range) into exa_sum(j -> body, Val(range)), hoisting Val(range) outside the closure boundary in a let block — necessary because Julia's inference cannot propagate constants through closure boundaries.

NaNSource{T} replaces nothing for sparsity-detection passes, giving concrete types to the analysis and eliminating abstract dispatch paths.

Val-based sparsity detection in src/simdfunction.jl: _gr_val / _hr0_val dispatch trees compute gradient/Hessian NTuple sizes purely through types, replacing @generated functions and unique(Any[])/Set{Any} (which juliac cannot handle). Identity-based dedup via _ident_unique replaces Set.

AOT apps and tests: test/LuksanVlcekApp.jl and test/COPSApp.jl are standalone @main executables that compile with juliac --output-exe a --trim=safe --experimental. test/NLPModelsIpoptLite.jl is a minimal Ipopt C-API wrapper (via ccall to libipopt) designed for AOT use (no dynamic dispatch on the solver side, uses @cfunction with concrete model types). test/JuliaCTest/JuliaCTest.jl runs these compilations in CI.

---

3. Multi-generator add_con

A single add_con call can now accept multiple generators to build a constraint block where different generators contribute terms to the same set of constraint rows:

c, g = add_con(c,
    (x[i] for i in 1:N),                  # first contribution
    (i => sin(x[i]) for i in 1:N);        # second contribution (augmentation)
    lcon = 0.0, ucon = 0.0,
)

Internally, the first generator creates the constraint via the standard add_con path (allocating N constraint rows with the given bounds). Each subsequent generator is automatically forwarded to add_con!, which adds its terms to the same rows. The i => expr syntax in augmentation generators maps each expression to the appropriate constraint row via Pair-based indexing.

This is syntactic sugar for the explicit two-step pattern:

c, g = add_con(c, x[i] for i in 1:N; lcon=0.0, ucon=0.0)
c, _ = add_con!(c, g, i => sin(x[i]) for i in 1:N)

The multi-gen form is more concise and keeps all contributions to a constraint block visually grouped together. The keyword arguments (lcon, ucon, start, name) apply to the constraint block created by the first generator; subsequent generators inherit the block dimensions.

---

4. Subexpression changes (@expr / add_expr)

subexpr(core, gen; reduced, parameter_only) on main supported three modes:

1. Lifted (default, reduced=false): introduced auxiliary variables and equality constraints, returning a Subexpr handle indexable like a Variable.
2. Reduced (reduced=true): inlined the expression directly when indexed, returning a ReducedSubexpr — no auxiliary variables.
3. Parameter-only (parameter_only=true): evaluated once when parameters are set, cached values in θ, returning a ParameterSubexpr.

This PR keeps only the reduced (inline) mode, now exposed as add_expr / @expr → Expression:

• ✅ Continued: reduced/inline substitution — s[i] substitutes the expression tree directly into the enclosing objective or constraint. This is now the only behavior of add_expr / @expr. Multi-dimensional indexing is supported.
• ❌ Discontinued: lifted mode (auxiliary variables + equality constraints via Subexpr). This introduced abstract-typed constraint chains that blocked AOT compilation and added unnecessary NLP variables.
• ❌ Discontinued: parameter-only mode (ParameterSubexpr). The recomputation machinery (_recompute_param_subexprs!, param_subexpr_fns, etc.) is removed. Users who need parameter-dependent expressions can use @expr with @par-indexed values instead; recomputation happens naturally through the expression tree.

---

5. ExaModelsLinearAlgebra → ExaModelsStaticArrays (discontinued/replaced)

The ExaModelsLinearAlgebra extension on main (triggered by using LinearAlgebra) provided:

• Null-node pass-throughs for math functions (inv, abs, sqrt, etc.)
• Scalar Null arithmetic with zero/one elimination
• zero/one/conj/adjoint/transpose for AbstractNode

This extension is discontinued. Its functionality is split:

• The Null-node arithmetic and pass-throughs are now handled by the core library (built into specialization.jl and the Const/Null dispatch system).
• The StaticArrays integration (enabling det, cross, norm, dot, etc. on ExaModels expression nodes) is provided by a new ExaModelsStaticArrays extension, triggered by using StaticArrays (which pulls in LinearAlgebra as a co-trigger). This extension defines zero, one, conj, adjoint, transpose, LinearAlgebra.dot, and StaticArrays.arithmetic_closure for AbstractNode.

The LinAlgTest test module (1560 lines) covers AD correctness for all supported linear algebra operations.

---

6. Two-stage stochastic programming

The separate TwoStageExaModel wrapper (552-line src/two_stage.jl) is deleted. Two-stage support is now integrated into the core ExaCore via a tag system:

TwoStageTags{VI} (src/tags.jl) stores per-variable and per-constraint scenario indices. TwoStageExaCore is a type alias for ExaCore{T,VT,B,S} where S <: TwoStageTags. The @var/@con macros accept a scenario= keyword to assign scenario indices.

---

7. Variables constrained to a generator

c, x = add_var(c, f(p) for p in params)

Creates length(params) variables, each tied via an equality constraint to the corresponding generator value. Internally calls add_var for the block, then add_con to bind x[j] == gen.f(orig) for each element.

---

8. Bug fixes

• ExaCore(; backend=CUDABackend()) silently created a CPU model: the no-T constructor did not forward backend to the inner call. Fixed.
• GPU variable bounds: uses convert_array(b, backend) instead of adapt which silently copied GPU arrays to CPU.
• add_var(c; kwargs...) (scalar variable): now returns Var(offset+1) (an AbstractNode) so it's directly usable in expressions.
• obj_weight type safety: casts to eltype(...) to avoid Float64/Float32 mismatch on Metal and other Float32 backends.

---

9. Function registration overhaul

src/functionlist.jl and ext/ExaModelsSpecialFunctions.jl are regenerated by new dev-time scripts (deps/generate_functionlist.jl, deps/generate_specialfunctions.jl) using Symbolics.jl. All Float64 constant literals are replaced with type-generic helpers (_clog2(x), _cpi(x), etc.) for Float32 compatibility. New functions: sinpi, cospi, sinc, expm1, atan(y,x), hypot, and more.

---

10. Backend / extension changes

• Deleted: ExaModelsAMDGPU.jl, ExaModelsCUDA.jl, ExaModelsOneAPI.jl — GPU array conversion now handled via Adapt.jl / KernelAbstractions only.
• Simplified: ExaModelsOpenCL.jl — only sort! override remains.
• New: ExaModelsMetal.jl — Float64→Float32 conversion, default_T(::MetalBackend) = Float32.
• New: ExaModelsStaticArrays.jl — replaces ExaModelsLinearAlgebra (see §5).

---

11. Dependency changes

	Before	After
Adapt	indirect	direct dependency
StaticArrays	—	weak dep (new)
CUDA/AMDGPU/oneAPI	explicit extensions	via KernelAbstractions only
ExaModelsLinearAlgebra	extension	removed (→ ExaModelsStaticArrays)

Package version: 0.9.3 → 0.9.7.

[github-actions]

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/callback_benchmark.jl b/callback_benchmark.jl
index eb5bdfd..a879cb4 100644
--- a/callback_benchmark.jl
+++ b/callback_benchmark.jl
@@ -6,11 +6,13 @@ using Printf
 function btime(f, N)
     f()  # warmup
     GC.gc()
-    return minimum(begin
-        t = time_ns()
-        f()
-        (time_ns() - t) / 1e9
-    end for _ in 1:N)
+    return minimum(
+        begin
+                t = time_ns()
+                f()
+                (time_ns() - t) / 1.0e9
+            end for _ in 1:N
+    )
 end
 
 function benchmark_callbacks(m; N = 20)
@@ -19,39 +21,42 @@ function benchmark_callbacks(m; N = 20)
     nnzj = m.meta.nnzj
     nnzh = m.meta.nnzh
 
-    x    = copy(m.meta.x0)
-    y    = zeros(ncon)
-    c    = zeros(ncon)
-    g    = zeros(nvar)
-    jac  = zeros(nnzj)
+    x = copy(m.meta.x0)
+    y = zeros(ncon)
+    c = zeros(ncon)
+    g = zeros(nvar)
+    jac = zeros(nnzj)
     hess = zeros(nnzh)
 
-    tobj  = btime(() -> ExaModels.obj(m, x), N)
-    tcon  = ncon > 0 ? btime(() -> ExaModels.cons!(m, x, c), N)  : NaN
+    tobj = btime(() -> ExaModels.obj(m, x), N)
+    tcon = ncon > 0 ? btime(() -> ExaModels.cons!(m, x, c), N) : NaN
     tgrad = btime(() -> ExaModels.grad!(m, x, g), N)
-    tjac  = ncon > 0 ? btime(() -> ExaModels.jac_coord!(m, x, jac), N)  : NaN
+    tjac = ncon > 0 ? btime(() -> ExaModels.jac_coord!(m, x, jac), N) : NaN
     thess = btime(() -> ExaModels.hess_coord!(m, x, y, hess), N)
 
-    return (nvar=nvar, ncon=ncon, tobj=tobj, tcon=tcon, tgrad=tgrad, tjac=tjac, thess=thess)
+    return (nvar = nvar, ncon = ncon, tobj = tobj, tcon = tcon, tgrad = tgrad, tjac = tjac, thess = thess)
 end
 
 function print_header(title)
     println()
-    println("=" ^ 80)
-    @printf("  %-20s  %6s %6s | %8s %8s %8s %8s %8s\n",
-        title, "nvar", "ncon", "obj", "cons", "grad", "jac", "hess")
-    println("=" ^ 80)
+    println("="^80)
+    @printf(
+        "  %-20s  %6s %6s | %8s %8s %8s %8s %8s\n",
+        title, "nvar", "ncon", "obj", "cons", "grad", "jac", "hess"
+    )
+    return println("="^80)
 end
 
 function print_row(name, r)
-    @printf("  %-20s  %6s %6s | %8.2e %8s %8.2e %8s %8.2e\n",
+    return @printf(
+        "  %-20s  %6s %6s | %8.2e %8s %8.2e %8s %8.2e\n",
         name,
-        r.nvar < 1_000_000 ? @sprintf("%6d", r.nvar) : @sprintf("%5.0fk", r.nvar/1000),
-        r.ncon < 1_000_000 ? @sprintf("%6d", r.ncon) : @sprintf("%5.0fk", r.ncon/1000),
+        r.nvar < 1_000_000 ? @sprintf("%6d", r.nvar) : @sprintf("%5.0fk", r.nvar / 1000),
+        r.ncon < 1_000_000 ? @sprintf("%6d", r.ncon) : @sprintf("%5.0fk", r.ncon / 1000),
         r.tobj,
-        isnan(r.tcon)  ? "  N/A   " : @sprintf("%8.2e", r.tcon),
+        isnan(r.tcon) ? "  N/A   " : @sprintf("%8.2e", r.tcon),
         r.tgrad,
-        isnan(r.tjac)  ? "  N/A   " : @sprintf("%8.2e", r.tjac),
+        isnan(r.tjac) ? "  N/A   " : @sprintf("%8.2e", r.tjac),
         r.thess,
     )
 end
@@ -71,22 +76,22 @@ end
 
 # ── COPS ───────────────────────────────────────────────────────────────────────
 const COPS_INSTANCES = [
-    (:bearing_model,          (50, 50)),
-    (:chain_model,            (800,)),
-    (:camshape_model,         (1000,)),
-    (:catmix_model,           (500,)),
-    (:channel_model,          (1000,)),
-    (:gasoil_model,           (500,)),
-    (:glider_model,           (500,)),
-    (:marine_model,           (500,)),
-    (:methanol_model,         (500,)),
-    (:minsurf_model,          (100, 100)),
-    (:pinene_model,           (500,)),
-    (:robot_model,            (500,)),
-    (:rocket_model,           (2000,)),
-    (:steering_model,         (1000,)),
-    (:torsion_model,          (100, 100)),
-    (:channel_model,          (1000,)),
+    (:bearing_model, (50, 50)),
+    (:chain_model, (800,)),
+    (:camshape_model, (1000,)),
+    (:catmix_model, (500,)),
+    (:channel_model, (1000,)),
+    (:gasoil_model, (500,)),
+    (:glider_model, (500,)),
+    (:marine_model, (500,)),
+    (:methanol_model, (500,)),
+    (:minsurf_model, (100, 100)),
+    (:pinene_model, (500,)),
+    (:robot_model, (500,)),
+    (:rocket_model, (2000,)),
+    (:steering_model, (1000,)),
+    (:torsion_model, (100, 100)),
+    (:channel_model, (1000,)),
 ]
 
 print_header("COPS")
@@ -95,7 +100,7 @@ for (sym, params) in COPS_INSTANCES
     try
         m = model_func(COPSBenchmark.ExaModelsBackend(), params...)
         r = benchmark_callbacks(m)
-        print_row("$sym($(join(params,',')))", r)
+        print_row("$sym($(join(params, ',')))", r)
     catch e
         @printf("  %-20s  ERROR: %s\n", string(sym), e)
     end
diff --git a/deps/generate_specialfunctions.jl b/deps/generate_specialfunctions.jl
index bb5222e..ec54e2c 100644
--- a/deps/generate_specialfunctions.jl
+++ b/deps/generate_specialfunctions.jl
@@ -152,7 +152,7 @@ end
 function main()
     content = generate()
 
-    if "--write" in ARGS
+    return if "--write" in ARGS
         outfile = joinpath(dirname(@__DIR__), "ext", "functionlist.jl")
         write(outfile, content)
         println("Wrote $(length(content)) bytes to $outfile")
diff --git a/docs/src/distillation.jl b/docs/src/distillation.jl
index 52ed949..232f2cd 100644
--- a/docs/src/distillation.jl
+++ b/docs/src/distillation.jl
@@ -46,14 +46,14 @@ function distillation_column_model(T = 3; backend = nothing)
 
     c = ExaCore(; backend)
 
-    @var(c, xA, 0:T, 0:(NT+1); start = 0.5)
-    @var(c, yA, 0:T, 0:(NT+1); start = 0.5)
+    @var(c, xA, 0:T, 0:(NT + 1); start = 0.5)
+    @var(c, yA, 0:T, 0:(NT + 1); start = 0.5)
     @var(c, u, 0:T; start = 1.0)
     @var(c, V, 0:T; start = 1.0)
     @var(c, L2, 0:T; start = 1.0)
 
-    @obj(c, (yA[t, 1] - ybar)^2 for t = 0:T)
-    @obj(c, (u[t] - ubar)^2 for t = 0:T)
+    @obj(c, (yA[t, 1] - ybar)^2 for t in 0:T)
+    @obj(c, (u[t] - ubar)^2 for t in 0:T)
 
     @con(c, xA[0, i] - xA0 for (i, xA0) in xA0s)
     @con(
@@ -86,8 +86,8 @@ function distillation_column_model(T = 3; backend = nothing)
         (1 / Ar) * (L2[t] * xA[t, NT] - (F - D) * xA[t, NT+1] - V[t] * yA[t, NT+1]) for
         t = 1:T
     )
-    @con(c, V[t] - u[t] * D - D for t = 0:T)
-    @con(c, L2[t] - u[t] * D - F for t = 0:T)
+    @con(c, V[t] - u[t] * D - D for t in 0:T)
+    @con(c, L2[t] - u[t] * D - F for t in 0:T)
     @con(
         c,
         yA[t, i] * (1 - xA[t, i]) - alpha * xA[t, i] * (1 - yA[t, i]) for (t, i) in itr2
diff --git a/docs/src/gpu.jl b/docs/src/gpu.jl
index 93f88a6..fc49707 100644
--- a/docs/src/gpu.jl
+++ b/docs/src/gpu.jl
@@ -23,9 +23,9 @@ end
 function luksan_vlcek_model(N)
 
     c = ExaCore()
-    @var(c, x, N; start = (luksan_vlcek_x0(i) for i = 1:N))
-    @con(c, luksan_vlcek_con(x, i) for i = 1:(N-2))
-    @obj(c, luksan_vlcek_obj(x, i) for i = 2:N)
+    @var(c, x, N; start = (luksan_vlcek_x0(i) for i in 1:N))
+    @con(c, luksan_vlcek_con(x, i) for i in 1:(N - 2))
+    @obj(c, luksan_vlcek_obj(x, i) for i in 2:N)
 
     return ExaModel(c)
 end
@@ -34,9 +34,9 @@ end
 function luksan_vlcek_model(N, backend = nothing)
 
     c = ExaCore(; backend = backend) # specify the backend
-    @var(c, x, N; start = (luksan_vlcek_x0(i) for i = 1:N))
-    @con(c, luksan_vlcek_con(x, i) for i = 1:(N-2))
-    @obj(c, luksan_vlcek_obj(x, i) for i = 2:N)
+    @var(c, x, N; start = (luksan_vlcek_x0(i) for i in 1:N))
+    @con(c, luksan_vlcek_con(x, i) for i in 1:(N - 2))
+    @obj(c, luksan_vlcek_obj(x, i) for i in 2:N)
 
     return ExaModel(c)
 end
diff --git a/docs/src/guide.jl b/docs/src/guide.jl
index 4c48a4f..b0e5700 100644
--- a/docs/src/guide.jl
+++ b/docs/src/guide.jl
@@ -29,7 +29,7 @@ c = ExaCore()
 # ## Variables
 # Now, let's create the optimization variables. From the problem definition, we can see that we will need $N$ scalar variables. We will choose $N=10$, and create the variable $x\in\mathbb{R}^{N}$ with the following command:
 N = 10
-@var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
+@var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
 # This creates the variable `x`, which we will be able to refer to when we create constraints/objective constraints. Also, this modifies the information in the `ExaCore` object properly so that later an optimization model can be properly created with the necessary information. Observe that we have used the keyword argument `start` to specify the initial guess for the solution. The variable upper and lower bounds can be specified in a similar manner. For example, if we wanted to set the lower bound of the variable `x` to 0.0 and the upper bound to 10.0, we could do it as follows:
 # ```julia
 # @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N), lvar = 0.0, uvar = 10.0)
@@ -37,7 +37,7 @@ N = 10
 
 # ## Objective
 # The objective can be set as follows:
-@obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)
+@obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in 2:N)
 # !!! note
 #     Note that the terms here are summed, without explicitly using `sum( ... )` syntax.
 
diff --git a/docs/src/opf.jl b/docs/src/opf.jl
index 2a4199d..c1a44cc 100644
--- a/docs/src/opf.jl
+++ b/docs/src/opf.jl
@@ -101,7 +101,7 @@ function ac_power_model(filename; backend = nothing, T = Float64)
 
     w = ExaCore(T; backend = backend)
 
-    @var(w, va, length(data.bus);)
+    @var(w, va, length(data.bus))
 
     @var(
         w,
diff --git a/docs/src/parameters.jl b/docs/src/parameters.jl
index 6c86169..081e81e 100644
--- a/docs/src/parameters.jl
+++ b/docs/src/parameters.jl
@@ -14,10 +14,10 @@ c_param = ExaCore()
 
 # Define the variables as before:
 N = 10
-@var(c_param, x_p, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
+@var(c_param, x_p, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
 
 # Now we can use the parameters in our objective function just like variables:
-@obj(c_param, θ[1] * (x_p[i-1]^2 - x_p[i])^2 + (x_p[i-1] - θ[2])^2 for i = 2:N)
+@obj(c_param, θ[1] * (x_p[i - 1]^2 - x_p[i])^2 + (x_p[i - 1] - θ[2])^2 for i in 2:N)
 
 # Add the same constraints as before:
 @con(
diff --git a/docs/src/performance.jl b/docs/src/performance.jl
index a83fedb..f6ce053 100644
--- a/docs/src/performance.jl
+++ b/docs/src/performance.jl
@@ -9,8 +9,8 @@ using ExaModels
 t = @elapsed begin
     c = ExaCore()
     N = 10
-    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
-    @obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)
+    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
+    @obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in 2:N)
     @con(
         c,
         3x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2])sin(x[i+1] + x[i+2]) + 4x[i+1] - x[i]exp(x[i] - x[i+1]) - 3 for i = 1:(N-2)
@@ -24,8 +24,8 @@ println("$t seconds elapsed")
 t = @elapsed begin
     c = ExaCore()
     N = 10
-    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
-    @obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)
+    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
+    @obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in 2:N)
     @con(
         c,
         3x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2])sin(x[i+1] + x[i+2]) + 4x[i+1] - x[i]exp(x[i] - x[i+1]) - 3 for i = 1:(N-2)
@@ -39,8 +39,8 @@ println("$t seconds elapsed")
 # But instead, if you create a function, we can significantly reduce the model creation time.
 function luksan_vlcek_model(N)
     c = ExaCore()
-    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
-    @obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i = 2:N)
+    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
+    @obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in 2:N)
     @con(
         c,
         3x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2])sin(x[i+1] + x[i+2]) + 4x[i+1] -
@@ -74,8 +74,8 @@ function luksan_vlcek_model_concrete(N)
     arr1 = Array(2:N)
     arr2 = Array(1:(N-2))
 
-    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
-    @obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i in arr1)
+    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
+    @obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in arr1)
     @con(
         c,
         3x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2])sin(x[i+1] + x[i+2]) + 4x[i+1] -
@@ -90,8 +90,8 @@ function luksan_vlcek_model_non_concrete(N)
     arr1 = Array{Any}(2:N)
     arr2 = Array{Any}(1:(N-2))
 
-    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i = 1:N))
-    @obj(c, 100 * (x[i-1]^2 - x[i])^2 + (x[i-1] - 1)^2 for i in arr1)
+    @var(c, x, N; start = (mod(i, 2) == 1 ? -1.2 : 1.0 for i in 1:N))
+    @obj(c, 100 * (x[i - 1]^2 - x[i])^2 + (x[i - 1] - 1)^2 for i in arr1)
     @con(
         c,
         3x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2])sin(x[i+1] + x[i+2]) + 4x[i+1] -
diff --git a/docs/src/quad.jl b/docs/src/quad.jl
index 49a5ad4..7c5692e 100644
--- a/docs/src/quad.jl
+++ b/docs/src/quad.jl
@@ -22,11 +22,11 @@ function quadrotor_model(N = 3; backend = nothing)
 
     c = ExaCore(; backend = backend)
 
-    @var(c, x, 1:(N+1), 1:n)
+    @var(c, x, 1:(N + 1), 1:n)
     @var(c, u, 1:N, 1:p)
 
     @con(c, x[1, i] - x0 for (i, x0) in x0s)
-    @con(c, -x[i+1, 1] + x[i, 1] + (x[i, 2]) * dt for i = 1:N)
+    @con(c, -x[i + 1, 1] + x[i, 1] + (x[i, 2]) * dt for i in 1:N)
     @con(
         c,
         -x[i+1, 2] +
@@ -36,7 +36,7 @@ function quadrotor_model(N = 3; backend = nothing)
             u[i, 1] * sin(x[i, 7]) * sin(x[i, 9])
         ) * dt for i = 1:N
     )
-    @con(c, -x[i+1, 3] + x[i, 3] + (x[i, 4]) * dt for i = 1:N)
+    @con(c, -x[i + 1, 3] + x[i, 3] + (x[i, 4]) * dt for i in 1:N)
     @con(
         c,
         -x[i+1, 4] +
@@ -46,7 +46,7 @@ function quadrotor_model(N = 3; backend = nothing)
             u[i, 1] * sin(x[i, 7]) * cos(x[i, 9])
         ) * dt for i = 1:N
     )
-    @con(c, -x[i+1, 5] + x[i, 5] + (x[i, 6]) * dt for i = 1:N)
+    @con(c, -x[i + 1, 5] + x[i, 5] + (x[i, 6]) * dt for i in 1:N)
     @con(
         c,
         -x[i+1, 6] + x[i, 6] + (u[i, 1] * cos(x[i, 7]) * cos(x[i, 8]) - 9.8) * dt for
@@ -77,7 +77,7 @@ function quadrotor_model(N = 3; backend = nothing)
 
     @obj(c, 0.5 * R * (u[i, j]^2) for (i, j, R) in itr0)
     @obj(c, 0.5 * Q * (x[i, j] - d)^2 for (i, j, Q, d) in itr1)
-    @obj(c, 0.5 * Qf * (x[N+1, j] - d)^2 for (j, Qf, d) in itr2)
+    @obj(c, 0.5 * Qf * (x[N + 1, j] - d)^2 for (j, Qf, d) in itr2)
 
     m = ExaModel(c)
 
diff --git a/docs/src/two_stage.jl b/docs/src/two_stage.jl
index 44aa5ed..d5e836e 100644
--- a/docs/src/two_stage.jl
+++ b/docs/src/two_stage.jl
@@ -20,13 +20,13 @@ core = TwoStageExaCore()
 @var(core, d; start = 1.0, lvar = 0.0, uvar = Inf, scenario = 0)  ## design variable d
 
 # For the recourse variables `v`, we specify `scenario = [i for i=1:ns, j=1:nv]` to indicate that each variable `v[s,i]` belongs to scenario `s`. This allows us to define scenario-specific constraints and objectives that involve these recourse variables.
-@var(core, v, ns, nv; start = 1.0, lvar = 0.0, uvar = Inf, scenario = [i for i=1:ns, j=1:nv])  ## recourse variables v
+@var(core, v, ns, nv; start = 1.0, lvar = 0.0, uvar = Inf, scenario = [i for i in 1:ns, j in 1:nv])  ## recourse variables v
 
 # Now we can define the constraints and objective function. The `scenario` keyword argument in the `constraint` and `objective` functions allows us to specify which scenario(s) each constraint or objective term belongs to. 
-@con(core, v[s,1] - v[s,2]^2 for s in 1:ns; lcon = 0.0, scenario = 1:ns)
+@con(core, v[s, 1] - v[s, 2]^2 for s in 1:ns; lcon = 0.0, scenario = 1:ns)
 
 @obj(core, d^2)
-@obj(core, weight * (v[s,i] - d)^2 for s in 1:ns, i in 1:nv)
+@obj(core, weight * (v[s, i] - d)^2 for s in 1:ns, i in 1:nv)
 
 m = ExaModel(core)
 
diff --git a/ext/ExaModelsKernelAbstractions.jl b/ext/ExaModelsKernelAbstractions.jl
index d899059..95bde22 100644
--- a/ext/ExaModelsKernelAbstractions.jl
+++ b/ext/ExaModelsKernelAbstractions.jl
@@ -121,7 +121,7 @@ end
 _conaug_structure!(T, backend, ::Tuple{}, sparsity) = nothing
 function _conaug_structure!(T, backend, (con, cons...), sparsity)
     _conaug_structure!(T, backend, cons, sparsity)
-    con isa ExaModels.ConstraintAug && !isempty(con.itr) &&
+    return con isa ExaModels.ConstraintAug && !isempty(con.itr) &&
         kers(backend)(sparsity, con.f, con.itr, con.oa; ndrange = length(con.itr))
 end
 @kernel function kers(sparsity, @Const(f), @Const(itr), @Const(oa))
@@ -134,7 +134,7 @@ end
 _grad_structure!(T, backend, ::Tuple{}, gsparsity) = nothing
 function _grad_structure!(T, backend, (obj, objs...), gsparsity)
     _grad_structure!(T, backend, objs, gsparsity)
-    ExaModels.sgradient!(backend, gsparsity, obj, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN))
+    return ExaModels.sgradient!(backend, gsparsity, obj, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN))
 end
 
 function ExaModels.jac_structure!(
@@ -150,7 +150,7 @@ end
 _jac_structure!(T, backend, ::Tuple{}, rows, cols) = nothing
 function _jac_structure!(T, backend, (con, cons...), rows, cols)
     _jac_structure!(T, backend, cons, rows, cols)
-    ExaModels.sjacobian!(backend, rows, cols, con, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN))
+    return ExaModels.sjacobian!(backend, rows, cols, con, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN))
 end
 
 
@@ -169,12 +169,12 @@ end
 _obj_hess_structure!(T, backend, ::Tuple{}, rows, cols) = nothing
 function _obj_hess_structure!(T, backend, (obj, objs...), rows, cols)
     _obj_hess_structure!(T, backend, objs, rows, cols)
-    ExaModels.shessian!(backend, rows, cols, obj, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN), T(NaN))
+    return ExaModels.shessian!(backend, rows, cols, obj, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN), T(NaN))
 end
 _con_hess_structure!(T, backend, ::Tuple{}, rows, cols) = nothing
 function _con_hess_structure!(T, backend, (con, cons...), rows, cols)
     _con_hess_structure!(T, backend, cons, rows, cols)
-    ExaModels.shessian!(backend, rows, cols, con, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN), T(NaN))
+    return ExaModels.shessian!(backend, rows, cols, con, ExaModels.NaNSource{T}(), ExaModels.NaNSource{T}(), T(NaN), T(NaN))
 end
 
 
@@ -220,7 +220,7 @@ end
 _cons_nln!(backend, y, ::Tuple{}, x, θ) = nothing
 function _cons_nln!(backend, y, (con, cons...), x, θ)
     _cons_nln!(backend, y, cons, x, θ)
-    if con isa ExaModels.Constraint && !isempty(con.itr)
+    return if con isa ExaModels.Constraint && !isempty(con.itr)
         kerf(backend)(y, con.f, con.itr, x, θ; ndrange = length(con.itr))
     end
 end
@@ -230,7 +230,7 @@ end
 _conaugs!(backend, y, ::Tuple{}, x, θ) = nothing
 function _conaugs!(backend, y, (con, cons...), x, θ)
     _conaugs!(backend, y, cons, x, θ)
-    if con isa ExaModels.ConstraintAug && !isempty(con.itr)
+    return if con isa ExaModels.ConstraintAug && !isempty(con.itr)
         kerf2(backend)(y, con.f, con.itr, x, θ, con.oa; ndrange = length(con.itr))
     end
 end
@@ -261,7 +261,7 @@ end
 _grad!(backend, y, ::Tuple{}, x, θ) = nothing
 function _grad!(backend, y, (obj, objs...), x, θ)
     _grad!(backend, y, objs, x, θ)
-    ExaModels.sgradient!(backend, y, obj, x, θ, one(eltype(y)))
+    return ExaModels.sgradient!(backend, y, obj, x, θ, one(eltype(y)))
 end
 
 function ExaModels.jac_coord!(
@@ -276,7 +276,7 @@ end
 _jac_coord!(backend, y, ::Tuple{}, x, θ) = nothing
 function _jac_coord!(backend, y, (con, cons...), x, θ)
     _jac_coord!(backend, y, cons, x, θ)
-    ExaModels.sjacobian!(backend, y, nothing, con, x, θ, one(eltype(y)))
+    return ExaModels.sjacobian!(backend, y, nothing, con, x, θ, one(eltype(y)))
 end
 
 function ExaModels.jprod_nln!(
@@ -456,12 +456,12 @@ end
 _obj_hess_coord!(backend, hess, ::Tuple{}, x, θ, obj_weight) = nothing
 function _obj_hess_coord!(backend, hess, (obj, objs...), x, θ, obj_weight)
     _obj_hess_coord!(backend, hess, objs, x, θ, obj_weight)
-    ExaModels.shessian!(backend, hess, nothing, obj, x, θ, obj_weight, zero(eltype(hess)))
+    return ExaModels.shessian!(backend, hess, nothing, obj, x, θ, obj_weight, zero(eltype(hess)))
 end
 _con_hess_coord!(backend, hess, ::Tuple{}, x, θ, y) = nothing
 function _con_hess_coord!(backend, hess, (con, cons...), x, θ, y)
     _con_hess_coord!(backend, hess, cons, x, θ, y)
-    ExaModels.shessian!(backend, hess, nothing, con, x, θ, y, zero(eltype(hess)))
+    return ExaModels.shessian!(backend, hess, nothing, con, x, θ, y, zero(eltype(hess)))
 end
 
 
diff --git a/src/gradient.jl b/src/gradient.jl
index bce5d7a..2906550 100644
--- a/src/gradient.jl
+++ b/src/gradient.jl
@@ -99,44 +99,44 @@ end
 end
 
 @inline function grpass(
-    d::D,
-    comp::Nothing,
-    y,
-    o1,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-) where {D<:Union{AdjointNull,ParIndexed,Real}}
+        d::D,
+        comp::Nothing,
+        y,
+        o1,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+    ) where {D <: Union{AdjointNull, ParIndexed, Real}}
     return cnt
 end
 @inline function grpass(
-    d::D,
-    comp::Nothing,
-    y,
-    o1,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-) where {D<:AdjointNode1}
+        d::D,
+        comp::Nothing,
+        y,
+        o1,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+    ) where {D <: AdjointNode1}
     return grpass(d.inner, nothing, y, o1, cnt, adj * d.y)
 end
 @inline function grpass(
-    d::D,
-    comp::Nothing,
-    y,
-    o1,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-) where {D<:AdjointNode2}
+        d::D,
+        comp::Nothing,
+        y,
+        o1,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+    ) where {D <: AdjointNode2}
     cnt = grpass(d.inner1, nothing, y, o1, cnt, adj * d.y1)
     return grpass(d.inner2, nothing, y, o1, cnt, adj * d.y2)
 end
 @inline function grpass(
-    d::D,
-    comp::Nothing,
-    y,
-    o1,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-) where {D<:AdjointNodeVar}
+        d::D,
+        comp::Nothing,
+        y,
+        o1,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+    ) where {D <: AdjointNodeVar}
     mapping, uniques = cnt
     idx = _grpass_find_ident(d.i, uniques, 1)
     if idx === 0
diff --git a/src/graph.jl b/src/graph.jl
index 2ae3928..8c0af9c 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -351,12 +351,12 @@ end
 
 # ── Primal evaluation (x::AbstractVector → scalar) ───────────────────────────
 
-@inline (n::SumNode{Tuple{}})(i, x::V, θ) where {T, V<:AbstractVector{T}} = zero(T)
-@inline (n::SumNode)(i, x::V, θ) where {T, V<:AbstractVector{T}} =
+@inline (n::SumNode{Tuple{}})(i, x::V, θ) where {T, V <: AbstractVector{T}} = zero(T)
+@inline (n::SumNode)(i, x::V, θ) where {T, V <: AbstractVector{T}} =
     mapreduce(inner -> inner(i, x, θ), +, n.inners)
 
-@inline (n::ProdNode{Tuple{}})(i, x::V, θ) where {T, V<:AbstractVector{T}} = one(T)
-@inline (n::ProdNode)(i, x::V, θ) where {T, V<:AbstractVector{T}} =
+@inline (n::ProdNode{Tuple{}})(i, x::V, θ) where {T, V <: AbstractVector{T}} = one(T)
+@inline (n::ProdNode)(i, x::V, θ) where {T, V <: AbstractVector{T}} =
     mapreduce(inner -> inner(i, x, θ), *, n.inners)
 
 # ── Adjoint tree (gradient) ───────────────────────────────────────────────────
diff --git a/src/hessian.jl b/src/hessian.jl
index 64f92a0..5b9223b 100644
--- a/src/hessian.jl
+++ b/src/hessian.jl
@@ -268,15 +268,15 @@ end
 end
 
 @inline function hdrpass(
-    t1::T1,
-    t2::T2,
-    comp::Nothing,
-    y1,
-    y2,
-    o2,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-) where {T1<:SecondAdjointNodeVar,T2<:SecondAdjointNodeVar}
+        t1::T1,
+        t2::T2,
+        comp::Nothing,
+        y1,
+        y2,
+        o2,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+    ) where {T1 <: SecondAdjointNodeVar, T2 <: SecondAdjointNodeVar}
     pair = (t1.i, t2.i)
     mapping, uniques = cnt
     idx = _hpass_find_pair(pair, uniques, 1)
@@ -524,7 +524,7 @@ end
     y1,
     y2,
     o2,
-    cnt::Vector,
+        cnt::Vector,
     adj,
 )
     push!(cnt, (t1.i, t2.i))
@@ -544,15 +544,15 @@ end
 end
 
 @inline function hrpass(
-    t::T,
-    comp::Nothing,
-    y1,
-    y2,
-    o2,
-    cnt::Tuple{<:Tuple,<:Tuple},
-    adj,
-    adj2,
-) where {T<:SecondAdjointNodeVar}
+        t::T,
+        comp::Nothing,
+        y1,
+        y2,
+        o2,
+        cnt::Tuple{<:Tuple, <:Tuple},
+        adj,
+        adj2,
+    ) where {T <: SecondAdjointNodeVar}
     pair = (t.i, t.i)
     mapping, uniques = cnt
     idx = _hpass_find_pair(pair, uniques, 1)
diff --git a/src/nlp.jl b/src/nlp.jl
index c99900e..663dcca 100644
--- a/src/nlp.jl
+++ b/src/nlp.jl
@@ -16,11 +16,11 @@ A handle to a block of optimization variables added to an [`ExaCore`](@ref) via
 individual entries in objective and constraint expressions. Retrieve solution
 values with [`solution`](@ref).
 """
-struct Variable{D,S,O} <: AbstractVariable
+struct Variable{D, S, O} <: AbstractVariable
     size::S
     length::O
     offset::O
-    Variable(size::S, length::O, offset::O, D) where {S, O} = new{D,S,O}(size, length, offset)
+    Variable(size::S, length::O, offset::O, D) where {S, O} = new{D, S, O}(size, length, offset)
 end
 Base.show(io::IO, v::Variable) = print(
     io,
@@ -49,10 +49,10 @@ end
 Base.show(io::IO, s::Expression) = print(
     io,
     """
-Subexpression (reduced)
+    Subexpression (reduced)
 
-  s ∈ R^{$(join(size(s.size), " × "))}
-""",
+      s ∈ R^{$(join(size(s.size), " × "))}
+    """,
 )
 
 """
@@ -84,7 +84,7 @@ An objective term group added to an [`ExaCore`](@ref) via [`add_obj`](@ref) /
 [`@obj`](@ref). All `Objective` objects in a core are summed at evaluation time
 to form the total objective value.
 """
-struct Objective{F,I} <: AbstractObjective
+struct Objective{F, I} <: AbstractObjective
     f::F
     itr::I
 end
@@ -108,7 +108,7 @@ A block of constraints added to an [`ExaCore`](@ref) via [`add_con`](@ref) /
 Row `k` of this block maps to global constraint index `offset + k`. Dual
 solution values can be retrieved with [`multipliers`](@ref).
 """
-struct Constraint{F,I,O} <: AbstractConstraint
+struct Constraint{F, I, O} <: AbstractConstraint
     f::F
     itr::I
     offset::O
@@ -136,7 +136,7 @@ by `idx` at evaluation time. Multiple `ConstraintAug` objects can be stacked on
 the same base constraint to aggregate contributions from several data sources
 (e.g. summing arc flows into nodal balance constraints).
 """
-struct ConstraintAug{F,I,D} <: AbstractConstraint
+struct ConstraintAug{F, I, D} <: AbstractConstraint
     f::F
     itr::I
     oa::Int
@@ -155,7 +155,7 @@ Constraint Augmentation
 """,
 )
 
-abstract type AbstractExaCore{T,VT,B,S} end
+abstract type AbstractExaCore{T, VT, B, S} end
 
 """
     ExaCore([array_eltype::Type; backend = nothing, minimize = true, name = :Generic])
@@ -199,7 +199,7 @@ An ExaCore
   number of constraint patterns: ... 0

"""
-struct ExaCore{T,VT<:AbstractVector{T}, B, S, V, P, O, C, R} <: AbstractExaCore{T,VT,B,S}
+struct ExaCore{T, VT <: AbstractVector{T}, B, S, V, P, O, C, R} <: AbstractExaCore{T, VT, B, S}
name::Symbol
backend::B
var::V
@@ -228,32 +228,32 @@ struct ExaCore{T,VT<:AbstractVector{T}, B, S, V, P, O, C, R} <: AbstractExaCore{
end

@inline function _exa_core(

• ;
• name = :Generic,
• backend = nothing,
• var = (),
• par = (),
• obj = (),
• cons = (),
• nvar = 0,
• npar = 0,
• ncon = 0,
• nconaug = 0,
• nobj = 0,
• nnzc = 0,
• nnzg = 0,
• nnzj = 0,
• nnzh = 0,
• x0 = convert_array(zeros(default_T(backend), 0), backend),
• θ = similar(x0, 0),
• lvar = similar(x0),
• uvar = similar(x0),
• y0 = similar(x0),
• lcon = similar(x0),
• ucon = similar(x0),
• minimize = true,
• tags = nothing,
• refs = (;)

•    ;
  

•    name = :Generic,
  

•    backend = nothing,
  

•    var = (),
  

•    par = (),
  

•    obj = (),
  

•    cons = (),
  

•    nvar = 0,
  

•    npar = 0,
  

•    ncon = 0,
  

•    nconaug = 0,
  

•    nobj = 0,
  

•    nnzc = 0,
  

•    nnzg = 0,
  

•    nnzj = 0,
  

•    nnzh = 0,
  

•    x0 = convert_array(zeros(default_T(backend), 0), backend),
  

•    θ = similar(x0, 0),
  

•    lvar = similar(x0),
  

•    uvar = similar(x0),
  

•    y0 = similar(x0),
  

•    lcon = similar(x0),
  

•    ucon = similar(x0),
  

•    minimize = true,
  

•    tags = nothing,
  

•    refs = (;)
  

  )

  return ExaCore(
  @@ -285,9 +285,9 @@ end
  )
  end

-@inline ExaCore(::Type{T}; backend = nothing, kwargs...) where {T<:AbstractFloat} = _exa_core(; x0 = convert_array(zeros(T, 0), backend), backend, kwargs...)
+@inline ExaCore(::Type{T}; backend = nothing, kwargs...) where {T <: AbstractFloat} = _exa_core(; x0 = convert_array(zeros(T, 0), backend), backend, kwargs...)
@inline ExaCore(; backend = nothing, kwargs...) = ExaCore(default_T(backend); backend, kwargs...)
-@inline ExaCore(c::C; kwargs...) where C <: ExaCore = _exa_core(
+@inline ExaCore(c::C; kwargs...) where {C <: ExaCore} = _exa_core(
;
zip(fieldnames(C), ntuple(i -> getfield(c, i), Val(fieldcount(C))))...,
kwargs...,
@@ -309,7 +309,7 @@ An ExaCore
Backend: ......................... $B

number of objective patterns: .... $(depth(c.obj))

• number of constraint patterns: ... $(depth(c.cons))

•  number of constraint patterns: ... $(depth(c.cons))
  

""",
)

@@ -320,7 +320,7 @@ An abstract type for ExaModel, which is a subtype of `NLPModels.AbstractNLPModel
"""
abstract type AbstractExaModel{T,VT,E} <: NLPModels.AbstractNLPModel{T,VT} end

-struct ExaModel{T,VT,E,V,P,O,C,S,R} <: AbstractExaModel{T,VT,E}
+struct ExaModel{T, VT, E, V, P, O, C, S, R} <: AbstractExaModel{T, VT, E}
name::Symbol
vars::V
pars::P
@@ -380,7 +380,7 @@ julia> result = ipopt(m; print_level=0) # solve the problem

"""
-function ExaModel(c::C; prod = false, kwargs...) where {C<:ExaCore}
+function ExaModel(c::C; prod = false, kwargs...) where {C <: ExaCore}
    return ExaModel(
        c.name,
        c.var,
@@ -411,17 +411,17 @@ end

build_extension(c::ExaCore; kwargs...) = nothing

-@inline function Base.getindex(v::V, i) where {D, V<:Variable{D}}
+@inline function Base.getindex(v::V, i) where {D, V <: Variable{D}}
    _bound_check(v.size, i)
-    _indexed_var(i, v.offset - _start(v.size[1]) + 1)
+    return _indexed_var(i, v.offset - _start(v.size[1]) + 1)
end
# For symbolic (AbstractNode) indices: use Node2 directly so the offset (runtime Int)
# is stored as a plain Int64 child — giving concrete type Node2{+, I, Int64}.
# Going through _add_node_real would wrap the runtime offset in Val(d2::Int64),
# which is type-unstable and breaks juliac --trim=safe.
-@inline _indexed_var(i::I, o::Int) where {I<:AbstractNode} = Var(Node2(+, i, o))
+@inline _indexed_var(i::I, o::Int) where {I <: AbstractNode} = Var(Node2(+, i, o))
@inline _indexed_var(i, o) = Var(i + o)
-@inline function Base.getindex(v::V, is...) where {D, V<:Variable{D}}
+@inline function Base.getindex(v::V, is...) where {D, V <: Variable{D}}
    @assert(length(is) == length(v.size), "Variable index dimension error")
    _bound_check(v.size, is)
    Var(v.offset + idxx(is .- (_start.(v.size) .- 1), _length.(v.size)))
@@ -568,7 +568,7 @@ Variable
function add_var(
    c::C,
    ns...;
-    name = nothing,
+        name = nothing,
    start = zero(T),
    lvar = T(-Inf),
    uvar = T(Inf),
@@ -586,20 +586,20 @@ function add_var(
    append_var_tags(c.tags, c.backend, total(ns); kwargs...)
    v = Variable(ns, len, o, length(c.var) + 1)

-    (ExaCore(c; var = (v, c.var...), nvar=nvar, x0=x0, lvar=lvar, uvar=uvar, refs = add_refs(c.refs, name, v)), v)
+    return (ExaCore(c; var = (v, c.var...), nvar = nvar, x0 = x0, lvar = lvar, uvar = uvar, refs = add_refs(c.refs, name, v)), v)
end

# Explicit 2D overload — prevents juliac from widening ns to Tuple{T,Vararg{T}}
function add_var(
-    c::C,
-    n1::N1,
-    n2::N2;
-    name = nothing,
-    start = zero(T),
-    lvar = T(-Inf),
-    uvar = T(Inf),
-    kwargs...
-) where {T, C<:ExaCore{T}, N1<:Union{Integer,AbstractUnitRange}, N2<:Union{Integer,AbstractUnitRange}}
+        c::C,
+        n1::N1,
+        n2::N2;
+        name = nothing,
+        start = zero(T),
+        lvar = T(-Inf),
+        uvar = T(Inf),
+        kwargs...
+    ) where {T, C <: ExaCore{T}, N1 <: Union{Integer, AbstractUnitRange}, N2 <: Union{Integer, AbstractUnitRange}}
    ns = (n1, n2)
    o = c.nvar
    len = total(ns)
@@ -609,21 +609,21 @@ function add_var(
    uvar = append!(c.backend, c.uvar, uvar, len)
    append_var_tags(c.tags, c.backend, len; kwargs...)
    v = Variable(ns, len, o, length(c.var) + 1)
-    (ExaCore(c; var = (v, c.var...), nvar=nvar, x0=x0, lvar=lvar, uvar=uvar, refs = add_refs(c.refs, name, v)), v)
+    return (ExaCore(c; var = (v, c.var...), nvar = nvar, x0 = x0, lvar = lvar, uvar = uvar, refs = add_refs(c.refs, name, v)), v)
end

# Explicit 3D overload — prevents juliac from widening ns to Tuple{T,Vararg{T}}
function add_var(
-    c::C,
-    n1::N1,
-    n2::N2,
-    n3::N3;
-    name = nothing,
-    start = zero(T),
-    lvar = T(-Inf),
-    uvar = T(Inf),
-    kwargs...
-) where {T, C<:ExaCore{T}, N1<:Union{Integer,AbstractUnitRange}, N2<:Union{Integer,AbstractUnitRange}, N3<:Union{Integer,AbstractUnitRange}}
+        c::C,
+        n1::N1,
+        n2::N2,
+        n3::N3;
+        name = nothing,
+        start = zero(T),
+        lvar = T(-Inf),
+        uvar = T(Inf),
+        kwargs...
+    ) where {T, C <: ExaCore{T}, N1 <: Union{Integer, AbstractUnitRange}, N2 <: Union{Integer, AbstractUnitRange}, N3 <: Union{Integer, AbstractUnitRange}}
    ns = (n1, n2, n3)
    o = c.nvar
    len = total(ns)
@@ -633,7 +633,7 @@ function add_var(
    uvar = append!(c.backend, c.uvar, uvar, len)
    append_var_tags(c.tags, c.backend, len; kwargs...)
    v = Variable(ns, len, o, length(c.var) + 1)
-    (ExaCore(c; var = (v, c.var...), nvar=nvar, x0=x0, lvar=lvar, uvar=uvar, refs = add_refs(c.refs, name, v)), v)
+    return (ExaCore(c; var = (v, c.var...), nvar = nvar, x0 = x0, lvar = lvar, uvar = uvar, refs = add_refs(c.refs, name, v)), v)
end

@inline add_refs(refs, ::Nothing, var) = refs
@@ -645,29 +645,29 @@ end
# The @var / @obj / @con / @par / @expr macros emit calls to these methods
# so that `name::Val` is a positional argument, not a keyword argument.

-@inline function add_var(c::C, ns_and_name::Val; kwargs...) where {T, C<:ExaCore{T}}
-    add_var(c; name = ns_and_name, kwargs...)
+@inline function add_var(c::C, ns_and_name::Val; kwargs...) where {T, C <: ExaCore{T}}
+    return add_var(c; name = ns_and_name, kwargs...)
end
-@inline function add_var(c::C, n1, name::Val; kwargs...) where {T, C<:ExaCore{T}}
-    add_var(c, n1; name = name, kwargs...)
+@inline function add_var(c::C, n1, name::Val; kwargs...) where {T, C <: ExaCore{T}}
+    return add_var(c, n1; name = name, kwargs...)
end
-@inline function add_var(c::C, n1, n2, name::Val; kwargs...) where {T, C<:ExaCore{T}}
-    add_var(c, n1, n2; name = name, kwargs...)
+@inline function add_var(c::C, n1, n2, name::Val; kwargs...) where {T, C <: ExaCore{T}}
+    return add_var(c, n1, n2; name = name, kwargs...)
end
-@inline function add_var(c::C, n1, n2, n3, name::Val; kwargs...) where {T, C<:ExaCore{T}}
-    add_var(c, n1, n2, n3; name = name, kwargs...)
+@inline function add_var(c::C, n1, n2, n3, name::Val; kwargs...) where {T, C <: ExaCore{T}}
+    return add_var(c, n1, n2, n3; name = name, kwargs...)
end
-@inline function add_var(c::C, gen::Base.Generator, name::Val; kwargs...) where {T, C<:ExaCore{T}}
-    add_var(c, gen; name = name, kwargs...)
+@inline function add_var(c::C, gen::Base.Generator, name::Val; kwargs...) where {T, C <: ExaCore{T}}
+    return add_var(c, gen; name = name, kwargs...)
end
-@inline function add_par(c::C, start::AbstractArray, name::Val) where {T, C<:ExaCore{T}}
-    add_par(c, start; name = name)
+@inline function add_par(c::C, start::AbstractArray, name::Val) where {T, C <: ExaCore{T}}
+    return add_par(c, start; name = name)
end
-@inline function add_obj(c::C, gen, name::Val) where {T, C<:ExaCore{T}}
-    add_obj(c, gen; name = name)
+@inline function add_obj(c::C, gen, name::Val) where {T, C <: ExaCore{T}}
+    return add_obj(c, gen; name = name)
end
-@inline function add_obj(c::C, expr::N, pars, name::Val) where {T, C<:ExaCore{T}, N<:AbstractNode}
-    add_obj(c, expr, pars; name = name)
+@inline function add_obj(c::C, expr::N, pars, name::Val) where {T, C <: ExaCore{T}, N <: AbstractNode}
+    return add_obj(c, expr, pars; name = name)
end


@@ -693,7 +693,7 @@ Parameter
  θ ∈ R^{10}

"""
-function add_par(c::C, start::AbstractArray; name = nothing) where {T,C<:ExaCore{T}}
+function add_par(c::C, start::AbstractArray; name = nothing) where {T, C <: ExaCore{T}}

 ns = Base.size(start)
 o = c.npar

@@ -701,7 +701,7 @@ function add_par(c::C, start::AbstractArray; name = nothing) where {T,C<:ExaCore
npar = c.npar + len
θ = append!(c.backend, c.θ, start, len)
p = Parameter(ns, len, o)

• (ExaCore(c; par = (p, c.par...), θ=θ, npar=npar, refs = add_refs(c.refs, name, p)), p)

• return (ExaCore(c; par = (p, c.par...), θ = θ, npar = npar, refs = add_refs(c.refs, name, p)), p)
  end

"""
@@ -737,7 +737,7 @@ function set_parameter!(c::ExaCore, param::Parameter, values::AbstractArray)
return nothing
end

-function add_var(c::C; kwargs...) where {T,C<:ExaCore{T}}
+function add_var(c::C; kwargs...) where {T, C <: ExaCore{T}}
c, v = add_var(c, 1; kwargs...)
return (c, v[1])
end
@@ -751,14 +751,14 @@ tying each to the corresponding generator expression.
Returns (core, Variable).
"""
function add_var(

• c::C,
• gen::Base.Generator;
• name = nothing,
• start = zero(T),
• lvar = T(-Inf),
• uvar = T(Inf),
• kwargs...
  -) where {T,C<:ExaCore{T}}

•    c::C,
  

•    gen::Base.Generator;
  

•    name = nothing,
  

•    start = zero(T),
  

•    lvar = T(-Inf),
  

•    uvar = T(Inf),
  

•    kwargs...
  

• ) where {T, C <: ExaCore{T}}
  gen = _adapt_gen(gen)
  n = length(gen.iter)
  c, x = add_var(c, n; name = name, start = start, lvar = lvar, uvar = uvar, kwargs...)
  @@ -769,8 +769,8 @@ function add_var(
  return (c, x)
  end

-function add_var(c::C, name::Symbol, args...; kwargs...) where {T,C<:ExaCore{T}}

• c, v= add_var(c, args...; name, kwargs...)
  +function add_var(c::C, name::Symbol, args...; kwargs...) where {T, C <: ExaCore{T}}

• c, v = add_var(c, args...; name, kwargs...)
  return c
  end

@@ -800,12 +800,12 @@ Objective
where |P| = 10

"""
-@inline function add_obj(c::C, gen; name = nothing) where {T, C<:ExaCore{T}}
+@inline function add_obj(c::C, gen; name = nothing) where {T, C <: ExaCore{T}}
    gen = _adapt_gen(gen)
    f = SIMDFunction(T, gen, c.nobj, c.nnzg, c.nnzh)
    pars = gen.iter

-    _add_objective(c, f, pars, name)
+    return _add_objective(c, f, pars, name)
end

"""
@@ -820,10 +820,10 @@ When `name` is given as `Val(:name)`, the objective is also accessible as
Prefer the generator form (`add_obj(core, gen)`) for typical use; this form
is intended for code that builds expression trees programmatically.
"""
-@inline function add_obj(c::C, expr::N, pars = 1:1; name = nothing) where {T,C<:ExaCore{T},N<:AbstractNode}
+@inline function add_obj(c::C, expr::N, pars = 1:1; name = nothing) where {T, C <: ExaCore{T}, N <: AbstractNode}
    f = _simdfunction(T, expr, c.nobj, c.nnzg, c.nnzh)

-      _add_objective(c, f, pars, name)
+    return _add_objective(c, f, pars, name)
end

@inline function _add_objective(c, f, pars, name = nothing)
@@ -833,7 +833,7 @@ end
    nnzh = c.nnzh + nitr * f.o2step

    obj = Objective(f, convert_array(pars, c.backend))
-    (ExaCore(c; nobj=nobj, nnzg=nnzg, nnzh=nnzh, obj=(obj, c.obj...), refs = add_refs(c.refs, name, obj)), obj)
+    return (ExaCore(c; nobj = nobj, nnzg = nnzg, nnzh = nnzh, obj = (obj, c.obj...), refs = add_refs(c.refs, name, obj)), obj)
end


@@ -871,7 +871,7 @@ Constraint
function add_con(
    c::C,
    gen::Base.Generator;
-    name = nothing,
+        name = nothing,
    start = zero(T),
    lcon = zero(T),
    ucon = zero(T),
@@ -882,7 +882,7 @@ function add_con(
    f = SIMDFunction(T, gen, c.ncon, c.nnzj, c.nnzh)
    pars = gen.iter

-    _add_constraint(c, f, pars, start, lcon, ucon, name; kwargs...)
+    return _add_constraint(c, f, pars, start, lcon, ucon, name; kwargs...)
end

"""
@@ -908,11 +908,11 @@ c, g = add_con(c,

"""
function add_con(

• c::C,
• gen::Base.Generator,
• gens::Base.Generator...;
• kwargs...
  -) where {T,C<:ExaCore{T}}

•    c::C,
  

•    gen::Base.Generator,
  

•    gens::Base.Generator...;
  

•    kwargs...
  

• ) where {T, C <: ExaCore{T}}
  c, con = add_con(c, gen; kwargs...)
  for g in gens
  c, _ = add_con!(c, con, g)
  @@ -936,7 +936,7 @@ function add_con(
  c::C,
  expr::N,
  pars = 1:1;

• name = nothing,

•    name = nothing,
  

  start = zero(T),
  lcon = zero(T),
  ucon = zero(T),
  @@ -945,7 +945,7 @@ function add_con(

  f = _simdfunction(T,expr, c.ncon, c.nnzj, c.nnzh)

• _add_constraint(c, f, pars, start, lcon, ucon, name; kwargs...)

• return _add_constraint(c, f, pars, start, lcon, ucon, name; kwargs...)
  end

"""
@@ -982,7 +982,7 @@ Constraint
function add_con(
c::C,
n;

• name = nothing,

•    name = nothing,
  

  start = zero(T),
  lcon = zero(T),
  ucon = zero(T),
  @@ -991,22 +991,22 @@ function add_con(

  f = _simdfunction(T, Null(nothing), c.ncon, c.nnzj, c.nnzh)

• _add_constraint(c, f, 1:n, start, lcon, ucon, name; kwargs...)

• return _add_constraint(c, f, 1:n, start, lcon, ucon, name; kwargs...)
  end

Positional-name forwarding for add_con — avoids Core.kwcall for name.

-@inline function add_con(c::C, gen::Base.Generator, name::Val; kwargs...) where {T, C<:ExaCore{T}}

• add_con(c, gen; name = name, kwargs...)
  +@inline function add_con(c::C, gen::Base.Generator, name::Val; kwargs...) where {T, C <: ExaCore{T}}

• return add_con(c, gen; name = name, kwargs...)
  end
  -@inline function add_con(c::C, expr::N, pars, name::Val; kwargs...) where {T, C<:ExaCore{T}, N<:AbstractNode}

• add_con(c, expr, pars; name = name, kwargs...)
  +@inline function add_con(c::C, expr::N, pars, name::Val; kwargs...) where {T, C <: ExaCore{T}, N <: AbstractNode}

• return add_con(c, expr, pars; name = name, kwargs...)
  end
  -@inline function add_con(c::C, n, name::Val; kwargs...) where {T, C<:ExaCore{T}}

• add_con(c, n; name = name, kwargs...)
  +@inline function add_con(c::C, n, name::Val; kwargs...) where {T, C <: ExaCore{T}}

• return add_con(c, n; name = name, kwargs...)
  end

-function _add_constraint(c::C, f, pars, start, lcon, ucon, name = nothing; kwargs...) where {C<:ExaCore}
+function _add_constraint(c::C, f, pars, start, lcon, ucon, name = nothing; kwargs...) where {C <: ExaCore}
nitr = length(pars)
o = c.ncon
ncon = c.ncon + nitr
@@ -1020,7 +1020,7 @@ function _add_constraint(c::C, f, pars, start, lcon, ucon, name = nothing; kwarg

 con = Constraint(f, convert_array(pars, c.backend), o)

• (ExaCore(c; ncon=ncon, nnzj=nnzj, nnzh=nnzh, y0=y0, lcon=lcon, ucon=ucon, cons=(con, c.cons...), refs = add_refs(c.refs, name, con)), con)

• return (ExaCore(c; ncon = ncon, nnzj = nnzj, nnzh = nnzh, y0 = y0, lcon = lcon, ucon = ucon, cons = (con, c.cons...), refs = add_refs(c.refs, name, con)), con)
  end

@@ -1088,13 +1088,13 @@ add_con!(c, bus, arc.bus => p[arc.i] for arc in data.arc) # add arc
add_con!(c, bus, gen.bus => -pg[gen.i] for gen in data.gen) # subtract generation

"""
-function add_con!(c::C, c1, gen::Base.Generator) where {T, C<:ExaCore{T}}
+function add_con!(c::C, c1, gen::Base.Generator) where {T, C <: ExaCore{T}}

    gen = _adapt_gen(gen)
    f = SIMDFunction(T, gen, offset0(c1, 0), c.nnzj, c.nnzh)
    pars = gen.iter

-    _add_constraint!(c, f, pars, _constraint_dims(c1))
+    return _add_constraint!(c, f, pars, _constraint_dims(c1))
end

# Extract the dimensions of the original constraint's iterator
@@ -1108,7 +1108,7 @@ function _add_constraint!(c, f, pars, dims)
    nnzj = c.nnzj + nitr * f.o1step
    nnzh = c.nnzh + nitr * f.o2step
    con = ConstraintAug(f, convert_array(pars, c.backend), oa, dims)
-    (ExaCore(c; nconaug=nconaug, nnzj=nnzj, nnzh=nnzh, cons=(con, c.cons...)), con)
+    return (ExaCore(c; nconaug = nconaug, nnzj = nnzj, nnzh = nnzh, cons = (con, c.cons...)), con)
end

# Helper to infer dimensions from iterator
@@ -1157,7 +1157,7 @@ c, s = add_expr(c, x[i, k]^2 for (i, k) in itr)
"""
# Positional-name forwarding for add_expr — avoids Core.kwcall for `name`.
@inline function add_expr(c::C, gen::Base.Generator, name::Val) where {T, C <: ExaCore{T}}
-    add_expr(c, gen; name = name)
+    return add_expr(c, gen; name = name)
end

function add_expr(c::C, gen::Base.Generator; name = nothing) where {T, C <: ExaCore{T}}
@@ -1179,7 +1179,7 @@ end
_jac_structure!(T, cons::Tuple{}, rows, cols) = nothing
@inline function _jac_structure!(T, cons::Tuple, rows, cols)
    _jac_structure!(T, Base.tail(cons), rows, cols)
-    sjacobian!(rows, cols, first(cons), NaNSource{T}(), NaNSource{T}(), T(NaN))
+    return sjacobian!(rows, cols, first(cons), NaNSource{T}(), NaNSource{T}(), T(NaN))
end

function hess_structure!(m::AbstractExaModel{T}, rows::AbstractVector, cols::AbstractVector) where T
@@ -1191,13 +1191,13 @@ end
_obj_hess_structure!(T, objs::Tuple{}, rows, cols) = nothing
@inline function _obj_hess_structure!(T, objs::Tuple, rows, cols)
    _obj_hess_structure!(T, Base.tail(objs), rows, cols)
-    shessian!(rows, cols, first(objs), NaNSource{T}(), NaNSource{T}(), T(NaN), T(NaN))
+    return shessian!(rows, cols, first(objs), NaNSource{T}(), NaNSource{T}(), T(NaN), T(NaN))
end

_con_hess_structure!(T, cons::Tuple{}, rows, cols) = nothing
@inline function _con_hess_structure!(T, cons::Tuple, rows, cols)
    _con_hess_structure!(T, Base.tail(cons), rows, cols)
-    shessian!(rows, cols, first(cons), NaNSource{T}(), NaNSource{T}(), T(NaN), T(NaN))
+    return shessian!(rows, cols, first(cons), NaNSource{T}(), NaNSource{T}(), T(NaN), T(NaN))
end

function obj(m::AbstractExaModel, x::AbstractVector)
@@ -1211,7 +1211,7 @@ end
    @simd for k in eachindex(obj.itr)
        @inbounds s += obj.f(obj.itr[k], x, θ)
    end
-    s
+    return s
end
_obj(obj::Tuple{}, x, θ) = zero(eltype(x))

@@ -1224,7 +1224,7 @@ end
@inline function _cons_nln!(cons::Tuple, x, θ, g)
    con = first(cons)
    _cons_nln!(Base.tail(cons), x, θ, g)
-    @simd for i in eachindex(con.itr)
+    return @simd for i in eachindex(con.itr)
        g[offset0(con, i)] += con.f(con.itr[i], x, θ)
    end
end
@@ -1240,7 +1240,7 @@ end

@inline function _grad!(objs::Tuple, x, θ, f)
    _grad!(Base.tail(objs), x, θ, f)
-    gradient!(f, first(objs), x, θ, one(eltype(f)))
+    return gradient!(f, first(objs), x, θ, one(eltype(f)))
end
_grad!(objs::Tuple{}, x, θ, f) = nothing

@@ -1253,7 +1253,7 @@ end
_jac_coord!(cons::Tuple{}, x, θ, jac) = nothing
@inline function _jac_coord!(cons::Tuple, x, θ, jac)
    _jac_coord!(Base.tail(cons), x, θ, jac)
-    sjacobian!(jac, nothing, first(cons), x, θ, one(eltype(jac)))
+    return sjacobian!(jac, nothing, first(cons), x, θ, one(eltype(jac)))
end

function jprod_nln!(m::AbstractExaModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
@@ -1265,7 +1265,7 @@ end
_jprod_nln!(cons::Tuple{}, x, θ, v, Jv) = nothing
@inline function _jprod_nln!(cons::Tuple, x, θ, v, Jv)
    _jprod_nln!(Base.tail(cons), x, θ, v, Jv)
-    sjacobian!((Jv, v), nothing, first(cons), x, θ, one(eltype(Jv)))
+    return sjacobian!((Jv, v), nothing, first(cons), x, θ, one(eltype(Jv)))
end

function jtprod_nln!(m::AbstractExaModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector)
@@ -1277,7 +1277,7 @@ end
_jtprod_nln!(cons::Tuple{}, x, θ, v, Jtv) = nothing
@inline function _jtprod_nln!(cons::Tuple, x, θ, v, Jtv)
    _jtprod_nln!(Base.tail(cons), x, θ, v, Jtv)
-    sjacobian!(nothing, (Jtv, v), first(cons), x, θ, one(eltype(Jtv)))
+    return sjacobian!(nothing, (Jtv, v), first(cons), x, θ, one(eltype(Jtv)))
end

function hess_coord!(
@@ -1307,13 +1307,13 @@ end
_obj_hess_coord!(objs::Tuple{}, x, θ, hess, obj_weight) = nothing
@inline function _obj_hess_coord!(objs::Tuple, x, θ, hess, obj_weight)
    _obj_hess_coord!(Base.tail(objs), x, θ, hess, obj_weight)
-    shessian!(hess, nothing, first(objs), x, θ, obj_weight, zero(eltype(hess)))
+    return shessian!(hess, nothing, first(objs), x, θ, obj_weight, zero(eltype(hess)))
end

_con_hess_coord!(cons::Tuple{}, x, θ, y, hess, obj_weight) = nothing
@inline function _con_hess_coord!(cons::Tuple, x, θ, y, hess, obj_weight)
    _con_hess_coord!(Base.tail(cons), x, θ, y, hess, obj_weight)
-    shessian!(hess, nothing, first(cons), x, θ, y, zero(eltype(hess)))
+    return shessian!(hess, nothing, first(cons), x, θ, y, zero(eltype(hess)))
end

function hprod!(
@@ -1345,13 +1345,13 @@ end
_obj_hprod!(objs::Tuple{}, x, θ, v, Hv, obj_weight) = nothing
@inline function _obj_hprod!(objs::Tuple, x, θ, v, Hv, obj_weight)
    _obj_hprod!(Base.tail(objs), x, θ, v, Hv, obj_weight)
-    shessian!((Hv, v), nothing, first(objs), x, θ, obj_weight, zero(eltype(Hv)))
+    return shessian!((Hv, v), nothing, first(objs), x, θ, obj_weight, zero(eltype(Hv)))
end

_con_hprod!(cons::Tuple{}, x, θ, y, v, Hv, obj_weight) = nothing
@inline function _con_hprod!(cons::Tuple, x, θ, y, v, Hv, obj_weight)
    _con_hprod!(Base.tail(cons), x, θ, y, v, Hv, obj_weight)
-    shessian!((Hv, v), nothing, first(cons), x, θ, y, zero(eltype(Hv)))
+    return shessian!((Hv, v), nothing, first(cons), x, θ, y, zero(eltype(Hv)))
end

@inbounds @inline offset0(a, i) = offset0(a.f, i)
@@ -1424,7 +1424,7 @@ function solution(result::SolverCore.AbstractExecutionStats, x)
    o = x.offset
    len = total(x.size)
    s = size(x.size)
-    return reshape(view(result.solution, (o+1):(o+len)), s...)
+    return reshape(view(result.solution, (o + 1):(o + len)), s...)
end

"""
@@ -1461,7 +1461,7 @@ function multipliers_L(result::SolverCore.AbstractExecutionStats, x)
    o = x.offset
    len = total(x.size)
    s = size(x.size)
-    return reshape(view(result.multipliers_L, (o+1):(o+len)), s...)
+    return reshape(view(result.multipliers_L, (o + 1):(o + len)), s...)
end

"""
@@ -1498,7 +1498,7 @@ function multipliers_U(result::SolverCore.AbstractExecutionStats, x)
    o = x.offset
    len = total(x.size)
    s = size(x.size)
-    return reshape(view(result.multipliers_U, (o+1):(o+len)), s...)
+    return reshape(view(result.multipliers_U, (o + 1):(o + len)), s...)
end

solution(result::SolverCore.AbstractExecutionStats, x::Var) = result.solution[x.i]
@@ -1542,8 +1542,8 @@ _adapt_gen(gen) = Base.Generator(gen.f, collect(gen.iter))
_adapt_gen(gen::Base.Generator{P}) where {P<:Union{AbstractArray,AbstractRange}} = gen

function Base.getproperty(core::E, name::Symbol) where {E <: Union{ExaCore, ExaModel}}
-    if hasfield(E, name)
-        getfield(core,name)
+    return if hasfield(E, name)
+        getfield(core, name)
    elseif hasfield(typeof(core.refs), name)
        getfield(core.refs, name)
    else
@@ -1792,9 +1792,9 @@ aug = @con!(c, g, i => sin(x[i+1]) for i in 4:6)
macro con!(exs...)
    isempty(exs) && error("@constraint! requires core, existing constraint, and generator arguments")
    parts, kwargs = _split_macro_args(exs)
-    core  = parts[1]
-    c1    = parts[2]
-    args  = [_auto_const_gen(a) for a in parts[3:end]]
+    core = parts[1]
+    c1 = parts[2]
+    args = [_auto_const_gen(a) for a in parts[3:end]]

    con = gensym(:con)
    return quote
diff --git a/src/register.jl b/src/register.jl
index a5f7022..7f16039 100644
--- a/src/register.jl
+++ b/src/register.jl
@@ -100,7 +100,7 @@ julia> @register_bivariate(relu23, drelu231, drelu232, ddrelu2311, ddrelu2312, d
macro register_bivariate(f, df1, df2, ddf11, ddf12, ddf22)
    return esc(
        quote
-            if ExaModels._needs_overload($f, Tuple{ExaModels.AbstractNode,ExaModels.AbstractNode})
+            if ExaModels._needs_overload($f, Tuple{ExaModels.AbstractNode, ExaModels.AbstractNode})
                @inline function $f(
                    d1::D1,
                    d2::D2,
@@ -109,7 +109,7 @@ macro register_bivariate(f, df1, df2, ddf11, ddf12, ddf22)
                end
            end

-            if ExaModels._needs_overload($f, Tuple{ExaModels.AbstractNode,Real})
+            if ExaModels._needs_overload($f, Tuple{ExaModels.AbstractNode, Real})
                @inline function $f(
                    d1::D1,
                    d2::D2,
@@ -118,7 +118,7 @@ macro register_bivariate(f, df1, df2, ddf11, ddf12, ddf22)
                end
            end

-            if ExaModels._needs_overload($f, Tuple{Real,ExaModels.AbstractNode})
+            if ExaModels._needs_overload($f, Tuple{Real, ExaModels.AbstractNode})
                @inline function $f(
                    d1::D1,
                    d2::D2,
diff --git a/src/simdfunction.jl b/src/simdfunction.jl
index 5c788df..1d19b99 100644
--- a/src/simdfunction.jl
+++ b/src/simdfunction.jl
@@ -42,7 +42,7 @@ Returns a `SIMDFunction` using the `gen`.
    _simdfunction(T, gen.f(ParSource()), o0, o1, o2)
end

-@inline function _simdfunction(T, f::F, o0, o1, o2) where {F<:Real}
+@inline function _simdfunction(T, f::F, o0, o1, o2) where {F <: Real}
    f = replace_T(T, f)
    SIMDFunction(
        f,
@@ -98,15 +98,15 @@ end
# Each function returns Val{N}() where N is encoded in the return type via dispatch.
# This lets ntuple(f, Val{N}()) produce a concrete NTuple{N,...}.

-@inline _add_vals(::Val{A}, ::Val{B}) where {A,B} = Val(A + B)
+@inline _add_vals(::Val{A}, ::Val{B}) where {A, B} = Val(A + B)

# --- Gradient count: mirrors grpass leaf dispatch ---
_gr_val(::Type{<:AdjointNodeVar}) = Val(1)
_gr_val(::Type{<:AdjointNull}) = Val(0)
_gr_val(::Type{<:Real}) = Val(0)
_gr_val(::Type{<:ParIndexed}) = Val(0)
-_gr_val(::Type{AdjointNode1{F,T,I}}) where {F,T,I} = _gr_val(I)
-_gr_val(::Type{AdjointNode2{F,T,I1,I2}}) where {F,T,I1,I2} =
+_gr_val(::Type{AdjointNode1{F, T, I}}) where {F, T, I} = _gr_val(I)
+_gr_val(::Type{AdjointNode2{F, T, I1, I2}}) where {F, T, I1, I2} =
    _add_vals(_gr_val(I1), _gr_val(I2))

# --- Hessian count: mirrors hrpass0 → hrpass → hdrpass dispatch ---
@@ -120,45 +120,47 @@ const _LinearHr1F = Union{
_hr0_val(::Type{<:SecondAdjointNodeVar}) = Val(0)
_hr0_val(::Type{<:SecondAdjointNull}) = Val(0)
_hr0_val(::Type{<:Real}) = Val(0)
-_hr0_val(::Type{SecondAdjointNode1{F,T,I}}) where {F<:_LinearHr1F,T,I} = _hr0_val(I)
-_hr0_val(::Type{SecondAdjointNode2{typeof(+),T,I1,I2}}) where {T,I1,I2} =
+_hr0_val(::Type{SecondAdjointNode1{F, T, I}}) where {F <: _LinearHr1F, T, I} = _hr0_val(I)
+_hr0_val(::Type{SecondAdjointNode2{typeof(+), T, I1, I2}}) where {T, I1, I2} =
    _add_vals(_hr0_val(I1), _hr0_val(I2))
-_hr0_val(::Type{SecondAdjointNode2{typeof(-),T,I1,I2}}) where {T,I1,I2} =
+_hr0_val(::Type{SecondAdjointNode2{typeof(-), T, I1, I2}}) where {T, I1, I2} =
    _add_vals(_hr0_val(I1), _hr0_val(I2))
_hr0_val(T::Type) = _hrpass_val(T)  # fallthrough to hrpass

_hrpass_val(::Type{<:SecondAdjointNull}) = Val(0)
_hrpass_val(::Type{<:Real}) = Val(0)
_hrpass_val(::Type{<:SecondAdjointNodeVar}) = Val(1)
-_hrpass_val(::Type{SecondAdjointNode1{F,T,I}}) where {F,T,I} = _hrpass_val(I)
-_hrpass_val(::Type{SecondAdjointNode2{F,T,I1,I2}}) where {F,T,I1,I2} =
+_hrpass_val(::Type{SecondAdjointNode1{F, T, I}}) where {F, T, I} = _hrpass_val(I)
+_hrpass_val(::Type{SecondAdjointNode2{F, T, I1, I2}}) where {F, T, I1, I2} =
    _add_vals(_add_vals(_hrpass_val(I1), _hrpass_val(I2)), _hdrpass_val(I1, I2))

_hdrpass_val(::Type{<:SecondAdjointNull}, ::Type) = Val(0)
_hdrpass_val(::Type, ::Type{<:SecondAdjointNull}) = Val(0)
_hdrpass_val(::Type{<:SecondAdjointNodeVar}, ::Type{<:SecondAdjointNodeVar}) = Val(1)
-_hdrpass_val(::Type{<:SecondAdjointNodeVar}, ::Type{SecondAdjointNode1{F,T,I}}) where {F,T,I} =
+_hdrpass_val(::Type{<:SecondAdjointNodeVar}, ::Type{SecondAdjointNode1{F, T, I}}) where {F, T, I} =
    _hdrpass_fixedvar_val(I)
-_hdrpass_val(::Type{SecondAdjointNode1{F,T,I}}, ::Type{<:SecondAdjointNodeVar}) where {F,T,I} =
+_hdrpass_val(::Type{SecondAdjointNode1{F, T, I}}, ::Type{<:SecondAdjointNodeVar}) where {F, T, I} =
    _hdrpass_fixedvar_val(I)
-_hdrpass_val(::Type{SecondAdjointNode1{F1,T1,I1}}, ::Type{SecondAdjointNode1{F2,T2,I2}}) where {F1,T1,I1,F2,T2,I2} =
+_hdrpass_val(::Type{SecondAdjointNode1{F1, T1, I1}}, ::Type{SecondAdjointNode1{F2, T2, I2}}) where {F1, T1, I1, F2, T2, I2} =
    _hdrpass_val(I1, I2)
-_hdrpass_val(::Type{SecondAdjointNode2{F,T,I1,I2}}, ::Type{SecondAdjointNode2{G,U,J1,J2}}) where {F,T,I1,I2,G,U,J1,J2} =
-    _add_vals(_add_vals(_hdrpass_val(I1,J1), _hdrpass_val(I1,J2)),
-              _add_vals(_hdrpass_val(I2,J1), _hdrpass_val(I2,J2)))
-_hdrpass_val(::Type{SecondAdjointNode2{F,T,I1,I2}}, ::Type{SecondAdjointNode1{G,U,J}}) where {F,T,I1,I2,G,U,J} =
+_hdrpass_val(::Type{SecondAdjointNode2{F, T, I1, I2}}, ::Type{SecondAdjointNode2{G, U, J1, J2}}) where {F, T, I1, I2, G, U, J1, J2} =
+    _add_vals(
+    _add_vals(_hdrpass_val(I1, J1), _hdrpass_val(I1, J2)),
+    _add_vals(_hdrpass_val(I2, J1), _hdrpass_val(I2, J2))
+)
+_hdrpass_val(::Type{SecondAdjointNode2{F, T, I1, I2}}, ::Type{SecondAdjointNode1{G, U, J}}) where {F, T, I1, I2, G, U, J} =
    _add_vals(_hdrpass_val(I1, J), _hdrpass_val(I2, J))
-_hdrpass_val(::Type{SecondAdjointNode1{F,T,I}}, ::Type{SecondAdjointNode2{G,U,J1,J2}}) where {F,T,I,G,U,J1,J2} =
+_hdrpass_val(::Type{SecondAdjointNode1{F, T, I}}, ::Type{SecondAdjointNode2{G, U, J1, J2}}) where {F, T, I, G, U, J1, J2} =
    _add_vals(_hdrpass_val(I, J1), _hdrpass_val(I, J2))
-_hdrpass_val(::Type{<:SecondAdjointNodeVar}, ::Type{SecondAdjointNode2{F,T,I1,I2}}) where {F,T,I1,I2} =
+_hdrpass_val(::Type{<:SecondAdjointNodeVar}, ::Type{SecondAdjointNode2{F, T, I1, I2}}) where {F, T, I1, I2} =
    _add_vals(_hdrpass_fixedvar_val(I1), _hdrpass_fixedvar_val(I2))
-_hdrpass_val(::Type{SecondAdjointNode2{F,T,I1,I2}}, ::Type{<:SecondAdjointNodeVar}) where {F,T,I1,I2} =
+_hdrpass_val(::Type{SecondAdjointNode2{F, T, I1, I2}}, ::Type{<:SecondAdjointNodeVar}) where {F, T, I1, I2} =
    _add_vals(_hdrpass_fixedvar_val(I1), _hdrpass_fixedvar_val(I2))

_hdrpass_fixedvar_val(::Type{<:SecondAdjointNodeVar}) = Val(1)
-_hdrpass_fixedvar_val(::Type{<:Union{SecondAdjointNull,Real}}) = Val(0)
-_hdrpass_fixedvar_val(::Type{SecondAdjointNode1{F,T,I}}) where {F,T,I} = _hdrpass_fixedvar_val(I)
-_hdrpass_fixedvar_val(::Type{SecondAdjointNode2{F,T,I1,I2}}) where {F,T,I1,I2} =
+_hdrpass_fixedvar_val(::Type{<:Union{SecondAdjointNull, Real}}) = Val(0)
+_hdrpass_fixedvar_val(::Type{SecondAdjointNode1{F, T, I}}) where {F, T, I} = _hdrpass_fixedvar_val(I)
+_hdrpass_fixedvar_val(::Type{SecondAdjointNode2{F, T, I1, I2}}) where {F, T, I1, I2} =
    _add_vals(_hdrpass_fixedvar_val(I1), _hdrpass_fixedvar_val(I2))


@@ -175,14 +177,14 @@ end
    i2 = replace_T(t, n.inner2)
    return Node2{F,typeof(i1),typeof(i2)}(i1, i2)
end
-@inline replace_T(t, n::Null{T}) where T <: Real = Null{t}(t(n.value))
+@inline replace_T(t, n::Null{T}) where {T <: Real} = Null{t}(t(n.value))
@inline replace_T(::Type{T1}, n::T2) where {T1, T2 <: Real} = T1(n)
@inline replace_T(::Type{T1}, ::Val{V}) where {T1, V} = Val(T1(V))
@inline function replace_T(t, n::SumNode{I}) where {I}
    inners = map(x -> replace_T(t, x), n.inners)
-    SumNode(inners)
+    return SumNode(inners)
end
@inline function replace_T(t, n::ProdNode{I}) where {I}
    inners = map(x -> replace_T(t, x), n.inners)
-    ProdNode(inners)
+    return ProdNode(inners)
end
diff --git a/src/specialization.jl b/src/specialization.jl
index 346f982..9a36b63 100644
--- a/src/specialization.jl
+++ b/src/specialization.jl
@@ -44,12 +44,12 @@ AbstractNodes, etc.) passes through unchanged.
# (not iterator vars, not function names, not array/dot bases) via _maybe_const().

const _MAYBE_CONST_REF = GlobalRef(@__MODULE__, :_maybe_const)
-const _EXA_SUM_REF    = GlobalRef(@__MODULE__, :exa_sum)
-const _EXA_PROD_REF   = GlobalRef(@__MODULE__, :exa_prod)
+const _EXA_SUM_REF = GlobalRef(@__MODULE__, :exa_sum)
+const _EXA_PROD_REF = GlobalRef(@__MODULE__, :exa_prod)

# Collect iterator variable names from the LHS of a for-clause
function _collect_iter_vars!(vars::Set{Symbol}, lhs)
-    if lhs isa Symbol
+    return if lhs isa Symbol
        push!(vars, lhs)
    elseif lhs isa Expr && lhs.head == :tuple
        for arg in lhs.args
@@ -85,7 +85,7 @@ function _process_generator(gen, outer_vars::Set{Symbol})

    # hoisted: accumulates (sym => val_expr) pairs that must be computed OUTSIDE
    # the generator closure for juliac constant propagation to resolve Val{N}.
-    hoisted = Pair{Symbol,Any}[]
+    hoisted = Pair{Symbol, Any}[]

    body = gen.args[1]
    if body isa Expr && body.head == :generator
@@ -113,8 +113,10 @@ end
# not an array base (first arg of ref[]), not a dot-access target.
# `hoisted` accumulates (sym => expr) pairs for Val bindings that must be
# evaluated outside the generator closure (for juliac constant propagation).
-function _wrap_free_symbols(expr, iter_vars::Set{Symbol},
-                            hoisted::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[])
+function _wrap_free_symbols(
+        expr, iter_vars::Set{Symbol},
+        hoisted::Vector{Pair{Symbol, Any}} = Pair{Symbol, Any}[]
+    )
    if expr isa Symbol
        return expr in iter_vars ? expr : Expr(:call, _MAYBE_CONST_REF, expr)
    elseif !(expr isa Expr)
@@ -129,20 +131,22 @@ function _wrap_free_symbols(expr, iter_vars::Set{Symbol},
        # hoisted Val(length(range)) binding.  Hoisting outside the generator closure
        # lets juliac --trim=safe resolve Val{N} concretely via constant propagation.
        if length(expr.args) == 2 &&
-           expr.args[2] isa Expr && expr.args[2].head == :generator &&
-           length(expr.args[2].args) == 2 &&   # body + exactly one for-spec
-           (fn === :sum || fn === :prod ||
-            fn isa GlobalRef && (fn.name === :sum || fn.name === :prod))
-            gen        = expr.args[2]
-            body       = gen.args[1]
-            spec       = gen.args[2]            # Expr(:(=), iter_var_lhs, range_expr)
-            iter_lhs   = spec.args[1]
+                expr.args[2] isa Expr && expr.args[2].head == :generator &&
+                length(expr.args[2].args) == 2 &&   # body + exactly one for-spec
+                (
+                fn === :sum || fn === :prod ||
+                    fn isa GlobalRef && (fn.name === :sum || fn.name === :prod)
+            )
+            gen = expr.args[2]
+            body = gen.args[1]
+            spec = gen.args[2]            # Expr(:(=), iter_var_lhs, range_expr)
+            iter_lhs = spec.args[1]
            range_expr = spec.args[2]
            nested_vars = copy(iter_vars)
            _collect_iter_vars!(nested_vars, iter_lhs)
            wrapped_body = _wrap_free_symbols(body, nested_vars, hoisted)
            ref = (fn === :sum || fn isa GlobalRef && fn.name === :sum) ?
-                  _EXA_SUM_REF : _EXA_PROD_REF
+                _EXA_SUM_REF : _EXA_PROD_REF
            # Hoist Val(range) outside the generator closure so juliac can
            # resolve the concrete type (e.g. Val{1:3}).
            val_sym = gensym(:exa_sum_val)
@@ -298,10 +302,10 @@ See [`exa_sum`](@ref) for supported iterators and juliac usage notes.
# resolve the concrete Val{UnitRange{Int64}(1,3)} type.
# Users can also call this directly: `exa_sum(j -> x[j], Val(1:nc))`
@inline function exa_sum(f, ::Val{r}) where {r}
-    _exa_sum_range(f, first(r), Val(length(r)))
+    return _exa_sum_range(f, first(r), Val(length(r)))
end
@inline function exa_prod(f, ::Val{r}) where {r}
-    _exa_prod_range(f, first(r), Val(length(r)))
+    return _exa_prod_range(f, first(r), Val(length(r)))
end

@inline _exa_sum_range(f, lo::Int, ::Val{N}) where {N} =
diff --git a/test/ADTest/ADTest.jl b/test/ADTest/ADTest.jl
index 66205fd..4f22afe 100644
--- a/test/ADTest/ADTest.jl
+++ b/test/ADTest/ADTest.jl
@@ -139,7 +139,7 @@ function sgradient(f, x)

    ff = f(ExaModels.VarSource())
    d = ff(ExaModels.Identity(), ExaModels.AdjointNodeSource(ExaModels.NaNSource{T}()), ExaModels.NaNSource{T}())
-    y1, a1 = ExaModels.grpass(d, nothing, nothing, ExaModels.NaNSource{T}(), ((),()), T(NaN))
+    y1, a1 = ExaModels.grpass(d, nothing, nothing, ExaModels.NaNSource{T}(), ((), ()), T(NaN))

    comp = ExaModels.Compressor(y1)

@@ -161,7 +161,7 @@ function sgradient_with_params(f, x, θ)

    ff =...*[Comment body truncated]*