Write every model as a recipe; all seventeen compile into one shared library

Project.toml

@@ -9,9 +9,11 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 [weakdeps]
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
+ExaModelsCompiler = "3d1e9a26-5b74-4f0c-9a2b-7c8f4e11d3a7"
 
 [extensions]
 COPSBenchmarkExaModels = "ExaModels"
+COPSBenchmarkExaModelsCompiler = "ExaModelsCompiler"
 COPSBenchmarkJuMP = "JuMP"
 
 [compat]

ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl

@@ -5,6 +5,17 @@ import COPSBenchmark
 import COPSBenchmark: ExaModelsBackend
 using ExaModels
 
+# ── Deferred callables in recipes ─────────────────────────────────────────────
+#
+# Every function a recipe defers — start generators, index-set builders, data
+# tables — must be NAMED and owned by THIS PACKAGE: not an anonymous closure,
+# not an extension-owned function, not anything in `Main`.  Two independent
+# mechanisms require it.  The serialized core carries the function's TYPE,
+# which another process can resolve only through the owning package's `PkgId`;
+# and an AOT-compiled library calls argument functions by NAME, which only a
+# package import makes reachable.  `Base.Fix2` over such a function is fine —
+# the wrapper's type parameters stay named.
+
 include("bearing.jl")
 include("camshape.jl")
 include("catmix.jl")

ext/COPSBenchmarkExaModels/bearing.jl

@@ -4,41 +4,49 @@
 # COPS 3.0 - November 2002
 # COPS 3.1 - March 2004
 
-@inline function COPSBenchmark.bearing_model(::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...)
+# Two sizes, so two placeholders.  The grid spacings depend on them and keep
+# their symbolic form; the starting field and the per-column weights are
+# comprehensions over the grid and travel as data.
+@inline function COPSBenchmark.bearing_recipe(
+    ::ExaModelsBackend; T = Float64, backend = nothing,
+)
     b = 10              # grid is (0,2*pi)x(0,2*b)
     e = 0.1             # eccentricity
 
+    core, nx, ny = ExaModels.ExaCore(T; backend = backend, nargs = Val(2))
+    d = ExaModels.ArgNode2(COPSBenchmark.bearing_data, nx, ny)
+
     hx = 2*pi / (nx+1)  # grid spacing
     hy = 2*b / (ny+1)   # grid spacing
-    area = 0.5*hx*hy    # area of triangle
-
-    wq(i) = (1.0 + e*cos((i-1)*hx))^3
-    v0 = [max(sin((i-1)*hx), 0.0) for i in 1:nx+2, j in 1:ny+2]
-
-    core = ExaModels.ExaCore(T; backend=backend, concrete = Val(true))
 
-    ExaModels.@add_var(core, v, 1:nx+2, 1:ny+2; lvar = 0.0, start=v0)
+    ExaModels.@add_var(core, v, 1:nx+2, 1:ny+2; lvar = 0.0, start = d.v0)
 
     ExaModels.@add_obj(
         core,
-        0.5*(hx*hy/6.0) * (wq(i) + 2*wq(i+1))*(((v[i+1,j]-v[i,j])/hx)^2 + ((v[i,j+1]-v[i,j])/hy)^2) for i in 1:nx+1, j in 1:ny+1
+        0.5*(hx*hy/6.0) * (w1 + 2*w2)*(((v[i+1,j]-v[i,j])/hx)^2 + ((v[i,j+1]-v[i,j])/hy)^2) for (i, j, w1, w2) in d.lower
     )
     ExaModels.@add_obj(
         core,
-        0.5*(hx*hy/6.0) * (2*wq(i) + 2*wq(i-1))*(((v[i-1,j]-v[i,j])/hx)^2 + ((v[i,j-1]-v[i,j])/hy)^2) for i in 2:nx+2, j in 2:ny+2
+        0.5*(hx*hy/6.0) * (2*w1 + 2*w2)*(((v[i-1,j]-v[i,j])/hx)^2 + ((v[i,j-1]-v[i,j])/hy)^2) for (i, j, w1, w2) in d.upper
     )
     ExaModels.@add_obj(
         core,
-        -hx*hy*e*sin((i-1)*hx)*v[i, j] for i in 1:nx+2, j in 1:ny+2
+        -hx*hy*e*s*v[i, j] for (i, j, s) in d.lin
     )
 
     ExaModels.@add_con(core, c1, v[i, 1] for i in 1:nx+2)
-    ExaModels.@add_con(core, c2, v[i, ny+2] for i in 1:nx+2)
+    ExaModels.@add_con(core, c2, v[i, k] for i in 1:nx+2, k in (ny+2):(ny+2))
     ExaModels.@add_con(core, c3, v[1, i] for i in 1:ny+2)
-    ExaModels.@add_con(core, c4, v[nx+2, i] for i in 1:ny+2)
+    ExaModels.@add_con(core, c4, v[k, i] for k in (nx+2):(nx+2), i in 1:ny+2)
 
-    return ExaModels.ExaModel(core; kwargs...)
+    return core
 end
 
+@inline COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny) = (nx, ny)
 
-
+@inline COPSBenchmark.bearing_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) =
+    ExaModels.ExaModel(
+        COPSBenchmark.bearing_recipe(b; T = T, backend = backend),
+        COPSBenchmark.bearing_args(b, nx, ny)...;
+        kwargs...,
+    )

ext/COPSBenchmarkExaModels/camshape.jl

@@ -4,19 +4,27 @@
 # COPS 3.0 - November 2002
 # COPS 3.1 - March 2004
 
-@inline function COPSBenchmark.camshape_model(::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...)
-
+# The structure, with the number of discretization points left open.  `d_theta`
+# depends on that number, so it cannot be a constant here — it is a deferred
+# expression, and resolves to an ordinary Float64 when the model is built.  The
+# rows that index the last point are written as generators over `n:n` so the
+# index arrives as data rather than as a number the structure has to know.
+@inline function COPSBenchmark.camshape_recipe(
+    ::ExaModelsBackend; T = Float64, backend = nothing,
+)
     R_v = 1.0         # design parameter related to the valve shape
     R_max = 2.0       # maximum allowed radius of the cam
     R_min = 1.0       # minimum allowed radius of the cam
     alpha = 1.5       # curvature limit parameter
 
-    d_theta = 2*pi/(5*(n+1))   # angle between discretization points
+    core, n = ExaModels.ExaCore(
+        T; backend = backend, minimize = false, nargs = Val(1),
+    )
 
-    core = ExaModels.ExaCore(T; backend= backend, minimize=false, concrete = Val(true))
+    d_theta = 2*pi/(5*(n+1))   # angle between discretization points
 
     # radius of the cam at discretization points
-    ExaModels.@add_var(core, r, 1:n; lvar =R_min, uvar = R_max, start=(R_min+R_max)/2.0)
+    ExaModels.@add_var(core, r, 1:n; lvar = R_min, uvar = R_max, start = (R_min+R_max)/2.0)
 
     ExaModels.@add_obj(core, (pi*R_v)/n * r[i] for i in 1:n)
 
@@ -39,12 +47,12 @@
     ExaModels.@add_con(
         core,
         c4,
-        - r[n-1]*r[n] - r[n]*R_max + 2*r[n-1]*R_max*cos(d_theta); lcon = -Inf, ucon = 0.0
+        - r[k-1]*r[k] - r[k]*R_max + 2*r[k-1]*R_max*cos(d_theta) for k in n:n; lcon = -Inf, ucon = 0.0
     )
     ExaModels.@add_con(
         core,
         c5,
-        - 2*R_max*r[n] + 2*r[n]^2*cos(d_theta); lcon = -Inf, ucon = 0.0
+        - 2*R_max*r[k] + 2*r[k]^2*cos(d_theta) for k in n:n; lcon = -Inf, ucon = 0.0
     )
     # Curvature
     ExaModels.@add_con(
@@ -60,9 +68,17 @@
     ExaModels.@add_con(
         core,
         c8,
-        (R_max - r[n]); lcon = -alpha*d_theta, ucon = alpha*d_theta
+        (R_max - r[k]) for k in n:n; lcon = -alpha*d_theta, ucon = alpha*d_theta
     )
 
-    return ExaModels.ExaModel(core; kwargs...)
+    return core
 end
 
+@inline COPSBenchmark.camshape_args(::ExaModelsBackend, n) = (n,)
+
+@inline COPSBenchmark.camshape_model(b::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...) =
+    ExaModels.ExaModel(
+        COPSBenchmark.camshape_recipe(b; T = T, backend = backend),
+        COPSBenchmark.camshape_args(b, n)...;
+        kwargs...,
+    )

ext/COPSBenchmarkExaModels/catmix.jl

@@ -3,12 +3,15 @@
 # COPS 3.0 - November 2002
 # COPS 3.1 - March 2004
 
-@inline function COPSBenchmark.catmix_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...)
+# The uniform bounds and zero starts were written as full arrays sized by nh;
+# as scalars they say the same thing without depending on the size.  The two
+# patterned starts do depend on it and travel as data.
+@inline function COPSBenchmark.catmix_recipe(
+    ::ExaModelsBackend; T = Float64, backend = nothing,
+)
     ne = 2
     nc = 3
 
-    h = T(1) / T(nh)   # Final time / nh (was Int/Int → Float64)
-
     rho = T[
         0.11270166537926,
         0.50000000000000,
@@ -21,61 +24,71 @@
     one_T = T(1)
     rho_index = [(i, rho[i]) for i in 1:nc]
 
-    c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true))
-    ExaModels.@add_var(c, u, nh, nc; lvar = zeros(T, nh, nc), uvar = ones(T, nh, nc), start = zeros(T, nh, nc))
-    ExaModels.@add_var(c, v, nh, ne; start = T[mod(j, ne) for i in 1:nh, j in 1:ne])
-    ExaModels.@add_var(c, w, nh, nc, ne; start = zeros(T, nh, nc, ne))
-    ExaModels.@add_var(c, pp, nh, nc, ne; start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne])
-    ExaModels.@add_var(c, Dpp, nh, nc, ne; start = zeros(T, nh, nc, ne))
-    ExaModels.@add_var(c, ppf, ne; start = T[mod(i,ne) for i in 1:ne])
+    core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1))
+    d = ExaModels.ArgCall(COPSBenchmark.catmix_data, (Val(T), nh))
+
+    h = one(T) / nh   # Final time / nh
+
+    ExaModels.@add_var(core, u, nh, nc; lvar = zero(T), uvar = one(T), start = zero(T))
+    ExaModels.@add_var(core, v, nh, ne; start = d.v_start)
+    ExaModels.@add_var(core, w, nh, nc, ne; start = zero(T))
+    ExaModels.@add_var(core, pp, nh, nc, ne; start = d.pp_start)
+    ExaModels.@add_var(core, Dpp, nh, nc, ne; start = zero(T))
+    ExaModels.@add_var(core, ppf, ne; start = T[mod(i,ne) for i in 1:ne])
 
-    ExaModels.@add_obj(c, neg_one + ppf[1] + ppf[2])
-    ExaModels.@add_obj(c, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc)
+    ExaModels.@add_obj(core, neg_one + ppf[1] + ppf[2])
+    ExaModels.@add_obj(core, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc)
 
     ExaModels.@add_con(
-        c,
+        core,
         c1,
         pp[i, k, s] - v[i, s] - h*sum(w[i, j, s]*(rho^j/T(factorial(j))) for j in 1:nc) for i=1:nh,  (k, rho) in rho_index, s=1:ne
     )
 
     ExaModels.@add_con(
-        c,
+        core,
         c2,
         Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/T(factorial(j-1))) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne
     )
 
     ExaModels.@add_con(
-        c,
+        core,
         c3,
-        ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / T(factorial(j)) for j in 1:nc) for s in 1:ne
+        ppf[s] - v[k, s] - h * sum(w[k, j, s] / T(factorial(j)) for j in 1:nc) for k in nh:nh, s in 1:ne
     )
 
     ExaModels.@add_con(
-        c,
+        core,
         c4,
         v[i, s] + sum(w[i, j, s] * h / T(factorial(j)) for j in 1:nc) - v[i+1, s] for i in 1:nh-1, s in 1:ne
     )
 
-
-
     ExaModels.@add_con(
-        c,
+        core,
         c5,
         Dpp[i,j,1] - u[i,j] * (ten_T*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc
     )
 
     ExaModels.@add_con(
-        c,
+        core,
         c6,
         Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - ten_T*pp[i,j,2]) + (one_T - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc
     )
 
-
     ExaModels.@add_con(
-        c,
+        core,
         c7,
         v[1, s] - bc for (s, bc) in [(i, bc[i]) for i in 1:ne]
     )
 
-    return ExaModels.ExaModel(c; kwargs...)
+    return core
 end
+
+@inline COPSBenchmark.catmix_args(::ExaModelsBackend, nh) = (nh,)
+
+@inline COPSBenchmark.catmix_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) =
+    ExaModels.ExaModel(
+        COPSBenchmark.catmix_recipe(b; T = T, backend = backend),
+        COPSBenchmark.catmix_args(b, nh)...;
+        kwargs...,
+    )

ext/COPSBenchmarkExaModels/chain.jl

@@ -8,24 +8,29 @@
 #   see "Benchmarking Optimization Software with COPS"
 #   Argonne National Labs Technical Report ANL/MCS-246 (2004)
 
-@inline function COPSBenchmark.chain_model(::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...)
-    nh = max(2, div(n - 4, 4))
-
+# The public size `n` is not the discretization: `nh = max(2, div(n-4, 4))` is.
+# `chain_args` does that reduction and builds the four starting curves, which
+# are comprehensions over nh; the recipe is written against nh directly.
+@inline function COPSBenchmark.chain_recipe(
+    ::ExaModelsBackend; T = Float64, backend = nothing,
+)
     L = 4
     a = 1
     b = 3
-    tmin = b > a ? 1 / 4 : 3 / 4
     tf = 1.0
+
+    c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1))
+    d = ExaModels.ArgNode1(COPSBenchmark.chain_data, nh)
+
     h = tf / nh
 
-    c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true))
-    ExaModels.@add_var(c, u, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
-    ExaModels.@add_var(c, x1, nh + 1; start = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1])
-    ExaModels.@add_var(c, x2, nh + 1; start = [(4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a) *
-        (4 * abs(b - a) * (k / nh - tmin)) for k in 1:nh+1])
-    ExaModels.@add_var(c, x3, nh + 1;  start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+    ExaModels.@add_var(c, u, nh + 1; start = d.u0)
+    ExaModels.@add_var(c, x1, nh + 1; start = d.x10)
+    ExaModels.@add_var(c, x2, nh + 1; start = d.x20)
+    ExaModels.@add_var(c, x3, nh + 1; start = d.x30)
 
-    ExaModels.@add_obj(c, x2[nh + 1])
+    # Indexes the last point, so the index arrives as data.
+    ExaModels.@add_obj(c, x2[k] for k in (nh+1):(nh+1))
 
     ExaModels.@add_con(
         c,
@@ -42,7 +47,7 @@
     ExaModels.@add_con(
         c,
         c3,
-        x1[nh + 1] - b
+        x1[k] - b for k in (nh+1):(nh+1)
     )
 
     ExaModels.@add_con(
@@ -60,7 +65,7 @@
     ExaModels.@add_con(
         c,
         c6,
-        x3[nh+1] - L
+        x3[k] - L for k in (nh+1):(nh+1)
     )
 
     ExaModels.@add_con(
@@ -75,6 +80,14 @@
         x3[j + 1] - x3[j] - 1 / 2 * h * (sqrt(1 + u[j]^2) + sqrt(1 + u[j + 1]^2)) for j in 1:nh
     )
 
-    return ExaModels.ExaModel(c; kwargs...)
+    return c
 end
 
+@inline COPSBenchmark.chain_args(::ExaModelsBackend, n) = (max(2, div(n - 4, 4)),)
+
+@inline COPSBenchmark.chain_model(b::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...) =
+    ExaModels.ExaModel(
+        COPSBenchmark.chain_recipe(b; T = T, backend = backend),
+        COPSBenchmark.chain_args(b, n)...;
+        kwargs...,
+    )