From a2b817d5d8939da231e3d962452cdbb2d70cdd6e Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Wed, 12 Aug 2026 18:43:57 -0400 Subject: [PATCH 1/7] Write ten problems as a recipe and the values that close it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each converted problem gains `_recipe` — its structure, with the size left open — and `_args`, the values that close it; `_model` becomes the two composed, so there is one definition rather than two that can drift. Every model is unchanged: at the sizes the test suite itself uses for callback comparison, all ten agree with the previous constructors on dimensions, bounds, starting point, `minimize`, objective, gradient, constraints, and the Jacobian and Hessian including COO order. Converted: camshape, rocket, steering, channel, gasoil, catmix, methanol, marine, pinene, chain. What the size no longer being a number required: * Scalars derived from it stay in the recipe as deferred expressions -- `h = tf/nh`, `d_theta = 2pi/(5(n+1))` -- including under exponentiation (`h^(k-1)`) and in bounds (`lcon = -alpha*d_theta`). Written `one(T)/nh` rather than `T(1)/T(nh)`: a constructor applied to a placeholder fails where arithmetic defers. * Uniform bounds and starts written as arrays sized by nh become scalars. catmix's `lvar = zeros(T, nh, nc)` is `zero(T)`, which says the same thing without depending on the size at all. * Genuinely patterned data -- comprehensions over the size -- moves to `_args`, bundled as one named tuple rather than many arguments, so the recipe reads `d.con1_itr` where the original read `con1_itr` and the diff stays reviewable through dense collocation arithmetic. * Rows indexing the last point become generators over `nh:nh`; where a fixed set rides along, a two-variable generator; where coefficients ride along, the index travels in the tuple with them (channel's bc3/bc4). methanol keeps scalar starts on purpose: it builds `v0` through three loops and then overwrites it wholesale with 0.001, so nothing size-dependent survives. robot and glider are NOT converted. Both use `@add_expr`, which stores its body as a closure; inside a model-building function that closure captures the `Variable`s, whose types carry the placeholder, and `instantiate` cannot reach inside a closure to resolve them. ExaModels refuses such a core rather than returning a half-instantiated model. Converting them would need either that fixed upstream or the expressions inlined by hand, which would change the sparsity pattern. Co-Authored-By: Claude Opus 5 --- ext/COPSBenchmarkExaModels/camshape.jl | 34 ++++-- ext/COPSBenchmarkExaModels/catmix.jl | 66 +++++++---- ext/COPSBenchmarkExaModels/chain.jl | 51 ++++++--- ext/COPSBenchmarkExaModels/channel.jl | 146 ++++++++++++++----------- ext/COPSBenchmarkExaModels/gasoil.jl | 134 +++++++++++++---------- ext/COPSBenchmarkExaModels/marine.jl | 127 +++++++++++---------- ext/COPSBenchmarkExaModels/methanol.jl | 124 ++++++++++----------- ext/COPSBenchmarkExaModels/pinene.jl | 98 ++++++++++------- ext/COPSBenchmarkExaModels/rocket.jl | 46 ++++++-- ext/COPSBenchmarkExaModels/steering.jl | 54 ++++++--- src/COPSBenchmark.jl | 20 ++++ 11 files changed, 536 insertions(+), 364 deletions(-) diff --git a/ext/COPSBenchmarkExaModels/camshape.jl b/ext/COPSBenchmarkExaModels/camshape.jl index e7484ad..f869d12 100644 --- a/ext/COPSBenchmarkExaModels/camshape.jl +++ b/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..., + ) diff --git a/ext/COPSBenchmarkExaModels/catmix.jl b/ext/COPSBenchmarkExaModels/catmix.jl index 34d419a..bbe1b8a 100644 --- a/ext/COPSBenchmarkExaModels/catmix.jl +++ b/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,76 @@ 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, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + + 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 function COPSBenchmark.catmix_args(::ExaModelsBackend, nh; T = Float64) + ne = 2 + nc = 3 + v_start = T[mod(j, ne) for i in 1:nh, j in 1:ne] + pp_start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne] + return (nh, (; v_start, pp_start)) +end + +@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; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/chain.jl b/ext/COPSBenchmarkExaModels/chain.jl index c7891c8..c5cf307 100644 --- a/ext/COPSBenchmarkExaModels/chain.jl +++ b/ext/COPSBenchmarkExaModels/chain.jl @@ -8,24 +8,28 @@ # 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, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + 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 +46,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 +64,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 +79,25 @@ 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 function COPSBenchmark.chain_args(::ExaModelsBackend, n; T = Float64) + nh = max(2, div(n - 4, 4)) + a = 1 + b = 3 + tmin = b > a ? 1 / 4 : 3 / 4 + u0 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] + x10 = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1] + x20 = [(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] + x30 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] + return (nh, (; u0, x10, x20, x30)) end +@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; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/channel.jl b/ext/COPSBenchmarkExaModels/channel.jl index cb0d639..14799db 100644 --- a/ext/COPSBenchmarkExaModels/channel.jl +++ b/ext/COPSBenchmarkExaModels/channel.jl @@ -5,14 +5,78 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.channel_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# Every table here is a function of the interval length h = tf/nh, so all of it +# is data rather than structure. It travels as a single named tuple rather than +# as ten separate arguments: a placeholder supports field access, so the recipe +# reads `d.con1_itr` exactly as the original read `con1_itr`. +@inline function COPSBenchmark.channel_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) + bc = T[0.0 1.0; 0.0 0.0] + + core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + + ExaModels.@add_var(core, v, nh, 4; start = d.v0) + ExaModels.@add_var(core, w, nh, 4; start = 0.0) + ExaModels.@add_var(core, uc, nh, 4, 4; start = d.uc0) + ExaModels.@add_var(core, Duc, nh, 4, 4; start = 0.0) + + # Constant objective + ExaModels.@add_obj(core, one(T) for _i in 1:1) + + # con1: uc[i,j,s] - v[i,s] - h * sum(w[i,k] * rho[j]^k / k!) + ExaModels.@add_con(core, c1, + uc[i, j, s] - v[i, s] - (a1*w[i,1] + a2*w[i,2] + a3*w[i,3] + a4*w[i,4]) + for (i, j, s, a1, a2, a3, a4) in d.con1_itr + ) + + # con2: Duc[i,j,s] - sum(v[i,k]*Bv[k]) - sum(w[i,k]*Bw[k]) + ExaModels.@add_con(core, c2, + Duc[i, j, s] - (bv1*v[i,1] + bv2*v[i,2] + bv3*v[i,3] + bv4*v[i,4]) + - (bw1*w[i,1] + bw2*w[i,2] + bw3*w[i,3] + bw4*w[i,4]) + for (i, j, s, bv1, bv2, bv3, bv4, bw1, bw2, bw3, bw4) in d.con2_itr + ) + + # Boundary conditions. bc3/bc4 index the last interval, so the index comes + # from a one-element set and the coefficients ride along with it. + ExaModels.@add_con(core, bc1, v[1, 1] - bc[1, 1]) + ExaModels.@add_con(core, bc2, v[1, 2] - bc[2, 1]) + + ExaModels.@add_con(core, bc3, + p1*v[k,1] + p2*v[k,2] + p3*v[k,3] + p4*v[k,4] + + q1*w[k,1] + q2*w[k,2] + q3*w[k,3] + q4*w[k,4] - bc[1, 2] + for (k, p1, p2, p3, p4, q1, q2, q3, q4) in d.bc3_itr + ) + + ExaModels.@add_con(core, bc4, + p1*v[k,2] + p2*v[k,3] + p3*v[k,4] + + q1*w[k,1] + q2*w[k,2] + q3*w[k,3] + q4*w[k,4] - bc[2, 2] + for (k, p1, p2, p3, q1, q2, q3, q4) in d.bc4_itr + ) + + # Continuity + ExaModels.@add_con(core, c3, + cv1*v[i,1] + cv2*v[i,2] + cv3*v[i,3] + cv4*v[i,4] + + cw1*w[i,1] + cw2*w[i,2] + cw3*w[i,3] + cw4*w[i,4] - v[i+1, s] + for (i, s, cv1, cv2, cv3, cv4, cw1, cw2, cw3, cw4) in d.cont_itr + ) + + # Collocation physics + ExaModels.@add_con(core, c4, + e1*w[i,1] + e2*w[i,2] + e3*w[i,3] + e4*w[i,4] - + T(10.0) * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4]) + for (i, j, e1, e2, e3, e4) in d.coll_itr + ) + + return core +end + +@inline function COPSBenchmark.channel_args(::ExaModelsBackend, nh; T = Float64) nc = 4 nd = 4 - R = T(10.0) tf = T(1.0) h = tf / nh - bc = T[0.0 1.0; 0.0 0.0] rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] t = T[(i-1)*h for i in 1:nh+1] @@ -29,22 +93,17 @@ # fac[k+1] = k! fac = T[factorial(k) for k in 0:nc+nd] - # Precompute all coefficients for the collocation constraints - # con1: uc[i,j,s] = v[i,s] + h * sum_{k=1}^{nc}(w[i,k] * rho[j]^k / k!) con1_itr = [ (i, j, s, h*rho[j]^1/fac[2], h*rho[j]^2/fac[3], h*rho[j]^3/fac[4], h*rho[j]^4/fac[5]) for i in 1:nh, j in 1:nc, s in 1:nd ] - # con2: Duc[i,j,s] coefficients depend on (j,s) con2_itr = [ (i, j, s, - # v coefficients: (rho[j]*h)^(k-s) / (k-s)! for k = s:nd (padded to length 4) (s <= 1 ? (rho[j]*h)^(1-s)/fac[1-s+1] : zero(T)), (s <= 2 ? (rho[j]*h)^(2-s)/fac[2-s+1] : zero(T)), (s <= 3 ? (rho[j]*h)^(3-s)/fac[3-s+1] : zero(T)), (s <= 4 ? (rho[j]*h)^(4-s)/fac[4-s+1] : zero(T)), - # w coefficients: h^(nd-s+1) * rho[j]^(k+nd-s) / (k+nd-s)! for k = 1:nc h^(nd-s+1)*rho[j]^(1+nd-s)/fac[1+nd-s+1], h^(nd-s+1)*rho[j]^(2+nd-s)/fac[2+nd-s+1], h^(nd-s+1)*rho[j]^(3+nd-s)/fac[3+nd-s+1], @@ -52,15 +111,12 @@ for i in 1:nh, j in 1:nc, s in 1:nd ] - # continuity coefficients: depend on s cont_itr = [ (i, s, - # v coefficients (s <= 1 ? h^(1-s)/fac[1-s+1] : zero(T)), (s <= 2 ? h^(2-s)/fac[2-s+1] : zero(T)), (s <= 3 ? h^(3-s)/fac[3-s+1] : zero(T)), (s <= 4 ? h^(4-s)/fac[4-s+1] : zero(T)), - # w coefficients h^(nd-s+1)/fac[1+nd-s+1], h^(nd-s+1)/fac[2+nd-s+1], h^(nd-s+1)/fac[3+nd-s+1], @@ -68,69 +124,27 @@ for i in 1:nh-1, s in 1:nd ] - # collocation physics coefficients: depend on j coll_itr = [ - (i, j, - rho[j]^0/fac[1], rho[j]^1/fac[2], rho[j]^2/fac[3], rho[j]^3/fac[4]) + (i, j, rho[j]^0/fac[1], rho[j]^1/fac[2], rho[j]^2/fac[3], rho[j]^3/fac[4]) for i in 1:nh, j in 1:nc ] - # right BC coefficients + # right BC coefficients, carried with the row index they apply to bc3_cv = [h^(k-1)/fac[k] for k in 1:nd] bc3_cw = [h^nd/fac[k+nd] for k in 1:nc] bc4_cv = [h^(k-2)/fac[k-1] for k in 2:nd] bc4_cw = [h^(nd-1)/fac[k+nd-1] for k in 1:nc] + bc3_itr = [(nh, bc3_cv[1], bc3_cv[2], bc3_cv[3], bc3_cv[4], + bc3_cw[1], bc3_cw[2], bc3_cw[3], bc3_cw[4])] + bc4_itr = [(nh, bc4_cv[1], bc4_cv[2], bc4_cv[3], + bc4_cw[1], bc4_cw[2], bc4_cw[3], bc4_cw[4])] - core = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - - ExaModels.@add_var(core, v, nh, nd; start = v0) - ExaModels.@add_var(core, w, nh, nc; start = 0.0) - ExaModels.@add_var(core, uc, nh, nc, nd; start = uc0) - ExaModels.@add_var(core, Duc, nh, nc, nd; start = 0.0) - - # Constant objective - ExaModels.@add_obj(core, one(T) for _i in 1:1) - - # con1: uc[i,j,s] - v[i,s] - h * sum(w[i,k] * rho[j]^k / k!) - ExaModels.@add_con(core, c1, - uc[i, j, s] - v[i, s] - (a1*w[i,1] + a2*w[i,2] + a3*w[i,3] + a4*w[i,4]) - for (i, j, s, a1, a2, a3, a4) in con1_itr - ) - - # con2: Duc[i,j,s] - sum(v[i,k]*Bv[k]) - sum(w[i,k]*Bw[k]) - ExaModels.@add_con(core, c2, - Duc[i, j, s] - (bv1*v[i,1] + bv2*v[i,2] + bv3*v[i,3] + bv4*v[i,4]) - - (bw1*w[i,1] + bw2*w[i,2] + bw3*w[i,3] + bw4*w[i,4]) - for (i, j, s, bv1, bv2, bv3, bv4, bw1, bw2, bw3, bw4) in con2_itr - ) - - # Boundary conditions - ExaModels.@add_con(core, bc1, v[1, 1] - bc[1, 1]) - ExaModels.@add_con(core, bc2, v[1, 2] - bc[2, 1]) - - ExaModels.@add_con(core, bc3, - bc3_cv[1]*v[nh,1] + bc3_cv[2]*v[nh,2] + bc3_cv[3]*v[nh,3] + bc3_cv[4]*v[nh,4] + - bc3_cw[1]*w[nh,1] + bc3_cw[2]*w[nh,2] + bc3_cw[3]*w[nh,3] + bc3_cw[4]*w[nh,4] - bc[1, 2] - ) - - ExaModels.@add_con(core, bc4, - bc4_cv[1]*v[nh,2] + bc4_cv[2]*v[nh,3] + bc4_cv[3]*v[nh,4] + - bc4_cw[1]*w[nh,1] + bc4_cw[2]*w[nh,2] + bc4_cw[3]*w[nh,3] + bc4_cw[4]*w[nh,4] - bc[2, 2] - ) - - # Continuity - ExaModels.@add_con(core, c3, - cv1*v[i,1] + cv2*v[i,2] + cv3*v[i,3] + cv4*v[i,4] + - cw1*w[i,1] + cw2*w[i,2] + cw3*w[i,3] + cw4*w[i,4] - v[i+1, s] - for (i, s, cv1, cv2, cv3, cv4, cw1, cw2, cw3, cw4) in cont_itr - ) + return (nh, (; v0, uc0, con1_itr, con2_itr, cont_itr, coll_itr, bc3_itr, bc4_itr)) +end - # Collocation physics - ExaModels.@add_con(core, c4, - e1*w[i,1] + e2*w[i,2] + e3*w[i,3] + e4*w[i,4] - - R * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4]) - for (i, j, e1, e2, e3, e4) in coll_itr +@inline COPSBenchmark.channel_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.channel_recipe(b; T = T, backend = backend), + COPSBenchmark.channel_args(b, nh; T = T)...; + kwargs..., ) - - return ExaModels.ExaModel(core; kwargs...) -end diff --git a/ext/COPSBenchmarkExaModels/gasoil.jl b/ext/COPSBenchmarkExaModels/gasoil.jl index 8b49476..dc2a3f1 100644 --- a/ext/COPSBenchmarkExaModels/gasoil.jl +++ b/ext/COPSBenchmarkExaModels/gasoil.jl @@ -5,80 +5,38 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.gasoil_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# `h`, the interval length, has a symbolic form and stays in the recipe. The +# tables that depend on it -- the measurement index map `itau`, the collocation +# times, and the starting values built from them -- do not, and travel as data. +@inline function COPSBenchmark.gasoil_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) nc = 4 # number of collocation points ne = 2 # number of differential equations np = 3 # number of ODE parameters nm = 21 # number of measurements - # roots of k-th degree Legendre polynomial rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] - # ODE initial conditions - bc = [1, 1, 2, 0] - # times at which observations made tau = T[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] - # ODEs defined in [0,tf] tf = tau[nm] - # uniform interval length - h = tf / T(nh) - t = T[T(i-1)*h for i in 1:nh+1] zero_T = T(0) + z1 = T[1.0000, 0.0000] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - # Concentrations - z = reshape(T[ - 1.0000, 0.0000, - 0.8105, 0.2000, - 0.6208, 0.2886, - 0.5258, 0.3010, - 0.4345, 0.3215, - 0.3903, 0.3123, - 0.3342, 0.2716, - 0.3034, 0.2551, - 0.2735, 0.2258, - 0.2405, 0.1959, - 0.2283, 0.1789, - 0.2071, 0.1457, - 0.1669, 0.1198, - 0.1530, 0.0909, - 0.1339, 0.0719, - 0.1265, 0.0561, - 0.1200, 0.0460, - 0.0990, 0.0280, - 0.0870, 0.0190, - 0.0770, 0.0140, - 0.0690, 0.0100, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - # Starting-value - for i in 1:itau[1], s in 1:ne - v0[i, s] = T(bc[s]) - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + h = tf / nh # ODE parameters ExaModels.@add_var(core, theta, 1:np; lvar = zero_T, start=zero_T) - # The collocation approximation u is defined by the parameters v and w. - # uc and Duc are, respectively, u and u' evaluated at the collocation points. - ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i =1:nh, s = 1:ne]) + ExaModels.@add_var(core, v, 1:nh, 1:ne; start=d.v_start) ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T) - ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i, s] for i =1:nh, j=1:nc, s = 1:ne]) + ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=d.uc_start) ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T) - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] itr2 = [(j,rho[j]) for j=1:nc] # L2 error - ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj, tauj, tj, zjs) in itr) + ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj, tauj, tj, zjs) in d.itr) # Collocation model ExaModels.@add_con( @@ -95,8 +53,8 @@ ) # Boundary - itr3 = [(s,z[1, s]) for s=1:ne] - ExaModels.@add_con(core, c3, - v[1, s] + z1s for (s,z1s) in itr3) #TODO + itr3 = [(s,z1[s]) for s=1:ne] + ExaModels.@add_con(core, c3, - v[1, s] + z1s for (s,z1s) in itr3) # Continuity ExaModels.@add_con( core, @@ -116,5 +74,67 @@ - Duc[i, j, 2] + theta[1]*uc[i,j,1]^2 - theta[2]*uc[i,j,2] for i=1:nh, j=1:nc ) - return ExaModels.ExaModel(core; kwargs...) + return core +end + +@inline function COPSBenchmark.gasoil_args(::ExaModelsBackend, nh; T = Float64) + nc = 4 + ne = 2 + nm = 21 + rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] + bc = [1, 1, 2, 0] + tau = T[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] + tf = tau[nm] + h = tf / T(nh) + t = T[T(i-1)*h for i in 1:nh+1] + + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(T[ + 1.0000, 0.0000, + 0.8105, 0.2000, + 0.6208, 0.2886, + 0.5258, 0.3010, + 0.4345, 0.3215, + 0.3903, 0.3123, + 0.3342, 0.2716, + 0.3034, 0.2551, + 0.2735, 0.2258, + 0.2405, 0.1959, + 0.2283, 0.1789, + 0.2071, 0.1457, + 0.1669, 0.1198, + 0.1530, 0.0909, + 0.1339, 0.0719, + 0.1265, 0.0561, + 0.1200, 0.0460, + 0.0990, 0.0280, + 0.0870, 0.0190, + 0.0770, 0.0140, + 0.0690, 0.0100, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = T(bc[s]) + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i = 1:nh, s = 1:ne] + uc_start = [v0[i, s] for i = 1:nh, j = 1:nc, s = 1:ne] + itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + + return (nh, (; v_start, uc_start, itr)) end + +@inline COPSBenchmark.gasoil_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.gasoil_recipe(b; T = T, backend = backend), + COPSBenchmark.gasoil_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/marine.jl b/ext/COPSBenchmarkExaModels/marine.jl index 51da740..80a19c0 100644 --- a/ext/COPSBenchmarkExaModels/marine.jl +++ b/ext/COPSBenchmarkExaModels/marine.jl @@ -5,79 +5,36 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.marine_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# Same shape as gasoil: `h` stays symbolic, the measurement table and the +# starting values built from it are data. `v` and `uc` share one start array in +# the original (nc = 1, so the shapes coincide) and share one field here. +@inline function COPSBenchmark.marine_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) nc = 1 # number of collocation points ne = 8 # number of differential equations nm = 21 # number of measurements - rho = T[0.5] # roots of k-th degree Legendre polynomial - tau = collect(T, range(T(0), T(10), 21)) # times at which observations made - tf = tau[nm] # ODEs defined in [0,tf] - h = tf / T(nh) # uniform interval length - t = T[T(i-1)*h for i in 1:nh+1] # partition + rho = T[0.5] + tf = T(10) # tau[nm] zero_T = T(0) - # itau[i] is the largest integer k with t[k] <= tau[i] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - # Observation - z = reshape(T[ - 20000.0, 17000.0, 10000.0, 15000.0, 12000.0, 9000.0, 7000.0, 3000.0, - 12445.0, 15411.0, 13040.0, 13338.0, 13484.0, 8426.0, 6615.0, 4022.0, - 7705.0, 13074.0, 14623.0, 11976.0, 12453.0, 9272.0, 6891.0, 5020.0, - 4664.0, 8579.0, 12434.0, 12603.0, 11738.0, 9710.0, 6821.0, 5722.0, - 2977.0, 7053.0, 11219.0, 11340.0, 13665.0, 8534.0, 6242.0, 5695.0, - 1769.0, 5054.0, 10065.0, 11232.0, 12112.0, 9600.0, 6647.0, 7034.0, - 943.0, 3907.0, 9473.0, 10334.0, 11115.0, 8826.0, 6842.0, 7348.0, - 581.0, 2624.0, 7421.0, 10297.0, 12427.0, 8747.0, 7199.0, 7684.0, - 355.0, 1744.0, 5369.0, 7748.0, 10057.0, 8698.0, 6542.0, 7410.0, - 223.0, 1272.0, 4713.0, 6869.0, 9564.0, 8766.0, 6810.0, 6961.0, - 137.0, 821.0, 3451.0, 6050.0, 8671.0, 8291.0, 6827.0, 7525.0, - 87.0, 577.0, 2649.0, 5454.0, 8430.0, 7411.0, 6423.0, 8388.0, - 49.0, 337.0, 2058.0, 4115.0, 7435.0, 7627.0, 6268.0, 7189.0, - 32.0, 228.0, 1440.0, 3790.0, 6474.0, 6658.0, 5859.0, 7467.0, - 17.0, 168.0, 1178.0, 3087.0, 6524.0, 5880.0, 5562.0, 7144.0, - 11.0, 99.0, 919.0, 2596.0, 5360.0, 5762.0, 4480.0, 7256.0, - 7.0, 65.0, 647.0, 1873.0, 4556.0, 5058.0, 4944.0, 7538.0, - 4.0, 44.0, 509.0, 1571.0, 4009.0, 4527.0, 4233.0, 6649.0, - 2.0, 27.0, 345.0, 1227.0, 3677.0, 4229.0, 3805.0, 6378.0, - 1.0, 20.0, 231.0, 934.0, 3197.0, 3695.0, 3159.0, 6454.0, - 1.0, 12.0, 198.0, 707.0, 2562.0, 3163.0, 3232.0, 5566.0, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - # Starting-value - for i in 1:itau[1], s in 1:ne - v0[i, s] = z[1, s] - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + h = tf / nh # uniform interval length # Growth rates ExaModels.@add_var(core, g, 1:ne-1; lvar = zero_T) # Mortality rates ExaModels.@add_var(core, m, 1:ne; lvar = zero_T) - # The collocation approximation u is defined by the parameters v and w. - # uc and Duc are, respectively, u and u' evaluated at the collocation points. - ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne]) + ExaModels.@add_var(core, v, 1:nh, 1:ne; start=d.v_start) ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T) - ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne]) + ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=d.v_start) ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T) - # error - - - # L2 error - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] itr2 = [(j,rho[j]) for j=1:nc] - ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj,tauj, tj, zjs) in itr) + ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj,tauj, tj, zjs) in d.itr) # Collocation model ExaModels.@add_con( @@ -120,5 +77,61 @@ for i=1:nh, j=1:nc, s=2:ne-1 ) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@inline function COPSBenchmark.marine_args(::ExaModelsBackend, nh; T = Float64) + ne = 8 + nm = 21 + tau = collect(T, range(T(0), T(10), 21)) + tf = tau[nm] + h = tf / T(nh) + t = T[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(T[ + 20000.0, 17000.0, 10000.0, 15000.0, 12000.0, 9000.0, 7000.0, 3000.0, + 12445.0, 15411.0, 13040.0, 13338.0, 13484.0, 8426.0, 6615.0, 4022.0, + 7705.0, 13074.0, 14623.0, 11976.0, 12453.0, 9272.0, 6891.0, 5020.0, + 4664.0, 8579.0, 12434.0, 12603.0, 11738.0, 9710.0, 6821.0, 5722.0, + 2977.0, 7053.0, 11219.0, 11340.0, 13665.0, 8534.0, 6242.0, 5695.0, + 1769.0, 5054.0, 10065.0, 11232.0, 12112.0, 9600.0, 6647.0, 7034.0, + 943.0, 3907.0, 9473.0, 10334.0, 11115.0, 8826.0, 6842.0, 7348.0, + 581.0, 2624.0, 7421.0, 10297.0, 12427.0, 8747.0, 7199.0, 7684.0, + 355.0, 1744.0, 5369.0, 7748.0, 10057.0, 8698.0, 6542.0, 7410.0, + 223.0, 1272.0, 4713.0, 6869.0, 9564.0, 8766.0, 6810.0, 6961.0, + 137.0, 821.0, 3451.0, 6050.0, 8671.0, 8291.0, 6827.0, 7525.0, + 87.0, 577.0, 2649.0, 5454.0, 8430.0, 7411.0, 6423.0, 8388.0, + 49.0, 337.0, 2058.0, 4115.0, 7435.0, 7627.0, 6268.0, 7189.0, + 32.0, 228.0, 1440.0, 3790.0, 6474.0, 6658.0, 5859.0, 7467.0, + 17.0, 168.0, 1178.0, 3087.0, 6524.0, 5880.0, 5562.0, 7144.0, + 11.0, 99.0, 919.0, 2596.0, 5360.0, 5762.0, 4480.0, 7256.0, + 7.0, 65.0, 647.0, 1873.0, 4556.0, 5058.0, 4944.0, 7538.0, + 4.0, 44.0, 509.0, 1571.0, 4009.0, 4527.0, 4233.0, 6649.0, + 2.0, 27.0, 345.0, 1227.0, 3677.0, 4229.0, 3805.0, 6378.0, + 1.0, 20.0, 231.0, 934.0, 3197.0, 3695.0, 3159.0, 6454.0, + 1.0, 12.0, 198.0, 707.0, 2562.0, 3163.0, 3232.0, 5566.0, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = z[1, s] + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i=1:nh, s=1:ne] + itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + return (nh, (; v_start, itr)) +end + +@inline COPSBenchmark.marine_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.marine_recipe(b; T = T, backend = backend), + COPSBenchmark.marine_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/methanol.jl b/ext/COPSBenchmarkExaModels/methanol.jl index 4a9835d..e44cb52 100644 --- a/ext/COPSBenchmarkExaModels/methanol.jl +++ b/ext/COPSBenchmarkExaModels/methanol.jl @@ -5,86 +5,33 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.methanol_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# `v0` is computed elaborately and then overwritten wholesale with 0.001, so +# every start in this model is that scalar -- nothing size-dependent survives. +# Only the measurement table, which depends on the interval length, is data. +@inline function COPSBenchmark.methanol_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) ne = 3 np = 5 nc = 3 - nm = 17 rho = T[0.11270166537926, 0.5, 0.88729833462074] - # times at which observations made - tau = T[ - 0., - 0.050, - 0.065, - 0.080, - 0.123, - 0.233, - 0.273, - 0.354, - 0.397, - 0.418, - 0.502, - 0.553, - 0.681, - 0.750, - 0.916, - 0.937, - 1.122, - ] - tf = tau[nm] # ODEs defined in [0,tf] - h = tf / T(nh) # uniform interval length - t = T[T(i-1)*h for i in 1:nh+1] # partition + tf = T(1.122) # tau[nm]; ODEs defined in [0,tf] zero_T = T(0) two_T = T(2) + bc = T[1, 0, 0] - # itau[i] is the largest integer k with t[k] <= tau[i] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + c, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - # Concentrations - z = reshape(T[ - 1.0000, 0.0000, 0.0000, - 0.7085, 0.1621, 0.0811, - 0.5971, 0.1855, 0.0965, - 0.5537, 0.1989, 0.1198, - 0.3684, 0.2845, 0.1535, - 0.1712, 0.3491, 0.2097, - 0.1198, 0.3098, 0.2628, - 0.0747, 0.3576, 0.2467, - 0.0529, 0.3347, 0.2884, - 0.0415, 0.3388, 0.2757, - 0.0261, 0.3557, 0.3167, - 0.0208, 0.3483, 0.2954, - 0.0085, 0.3836, 0.2950, - 0.0053, 0.3611, 0.2937, - 0.0019, 0.3609, 0.2831, - 0.0018, 0.3485, 0.2846, - 0.0006, 0.3698, 0.2899, - ], ne, nm)' - con1_matrix = [(j, s, itau[j], tau[j], z[j,s], t[itau[j]]) for j in 1:nm, s in 1:ne] - bc = T[1, 0, 0] + h = tf / nh # uniform interval length - # Starting-value - v0 = zeros(T, nh, ne) - for i in 1:itau[1], s in 1:ne - v0[i, s] = bc[s] - end - for j in 2:nm, i in itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - v0 .= T(0.001) - - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) ExaModels.@add_var(c, theta, np; lvar = zero_T, start = fill(T(1), np)) - ExaModels.@add_var(c, v, nh, ne; start = v0) + ExaModels.@add_var(c, v, nh, ne; start = T(0.001)) ExaModels.@add_var(c, w, nh, nc, ne; start = zero_T) - ExaModels.@add_var(c, uc, nh, nc, ne; start = [v0[i,s] for i=1:nh, j=1:nc, s=1:ne]) + ExaModels.@add_var(c, uc, nh, nc, ne; start = T(0.001)) ExaModels.@add_var(c, Duc, nh, nc, ne; start = zero_T) - ExaModels.@add_obj(c, (v[itau,s] + sum(w[itau,k,s]*(tau-t)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - z)^2 for (j,s,itau,tau,z,t) in con1_matrix) + ExaModels.@add_obj(c, (v[itau,s] + sum(w[itau,k,s]*(tau-t)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - z)^2 for (j,s,itau,tau,z,t) in d.con1_matrix) ExaModels.@add_con( c, @@ -132,5 +79,48 @@ theta[4]*uc[i,j,1]) for i=1:nh, j=1:nc ) - return ExaModels.ExaModel(c; kwargs...) + return c end + +@inline function COPSBenchmark.methanol_args(::ExaModelsBackend, nh; T = Float64) + ne = 3 + nm = 17 + tau = T[ + 0., 0.050, 0.065, 0.080, 0.123, 0.233, 0.273, 0.354, 0.397, 0.418, + 0.502, 0.553, 0.681, 0.750, 0.916, 0.937, 1.122, + ] + tf = tau[nm] + h = tf / T(nh) + t = T[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(T[ + 1.0000, 0.0000, 0.0000, + 0.7085, 0.1621, 0.0811, + 0.5971, 0.1855, 0.0965, + 0.5537, 0.1989, 0.1198, + 0.3684, 0.2845, 0.1535, + 0.1712, 0.3491, 0.2097, + 0.1198, 0.3098, 0.2628, + 0.0747, 0.3576, 0.2467, + 0.0529, 0.3347, 0.2884, + 0.0415, 0.3388, 0.2757, + 0.0261, 0.3557, 0.3167, + 0.0208, 0.3483, 0.2954, + 0.0085, 0.3836, 0.2950, + 0.0053, 0.3611, 0.2937, + 0.0019, 0.3609, 0.2831, + 0.0018, 0.3485, 0.2846, + 0.0006, 0.3698, 0.2899, + ], ne, nm)' + + con1_matrix = [(j, s, itau[j], tau[j], z[j,s], t[itau[j]]) for j in 1:nm, s in 1:ne] + return (nh, (; con1_matrix)) +end + +@inline COPSBenchmark.methanol_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.methanol_recipe(b; T = T, backend = backend), + COPSBenchmark.methanol_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/pinene.jl b/ext/COPSBenchmarkExaModels/pinene.jl index 485666d..a4432fc 100644 --- a/ext/COPSBenchmarkExaModels/pinene.jl +++ b/ext/COPSBenchmarkExaModels/pinene.jl @@ -5,63 +5,35 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 # -@inline function COPSBenchmark.pinene_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# Same collocation template as gasoil/marine. `h` keeps its symbolic form; the +# measurement index map and the starting values built from it are data. +@inline function COPSBenchmark.pinene_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) nc = 3 # number of collocation points ne = 5 # number of differential equations np = 5 # number of ODE parameters nm = 8 # number of measurements - # roots of k-th degree Legendre polynomial rho = T[0.11270166537926, 0.5, 0.88729833462074] - # boundary conditions bc = T[100, 0, 0, 0, 0] - # times at which observations made - tau = T[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] - tf = tau[nm] # ODEs defined in [0,tf] - h = tf / T(nh) # uniform interval length - t = T[T(i-1)*h for i in 1:nh+1] # partition + tf = T(36420) # tau[nm] zero_T = T(0) - # itau[i] is the largest integer k with t[k] <= tau[i] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - - # Observations - z = reshape(T[ - 88.35, 7.3, 2.3, 0.4, 1.75, - 76.4, 15.6, 4.5, 0.7, 2.8, - 65.1, 23.1, 5.3, 1.1, 5.8, - 50.4, 32.9, 6.0, 1.5, 9.3, - 37.5, 42.7, 6.0, 1.9, 12.0, - 25.9, 49.1, 5.9, 2.2, 17.0, - 14.0, 57.4, 5.1, 2.6, 21.0, - 4.5, 63.1, 3.8, 2.9, 25.7, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - # Starting-value - for i in 1:itau[1], s in 1:ne - v0[i, s] = bc[s] - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end + core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + h = tf / nh # uniform interval length ExaModels.@add_var(core, theta, 1:np; lvar = zero_T, start=zero_T) # The collocation approximation u is defined by the parameters v and w. # uc and Duc are, respectively, u and u' evaluated at the collocation points. - ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne]) + ExaModels.@add_var(core, v, 1:nh, 1:ne; start=d.v_start) ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T) - ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i,s] for i=1:nh, j=1:nc, s=1:ne]) + ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=d.uc_start) ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T) - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] # l2 error - ExaModels.@add_obj(core, (v[it,s] + sum(w[it,k,s]*(tj-ti)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j, s, it, tj, ti, zjs) in itr) + ExaModels.@add_obj(core, (v[it,s] + sum(w[it,k,s]*(tj-ti)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j, s, it, tj, ti, zjs) in d.itr) # Collocation model itr2 = [(j,rho[j]) for j=1:nc] @@ -92,5 +64,51 @@ ExaModels.@add_con(core, c8, -Duc[i,j,4] + theta[3]*uc[i,j,3] for i=1:nh, j=1:nc) ExaModels.@add_con(core, c9, -Duc[i,j,5] + theta[4]*uc[i,j,3] - theta[5]*uc[i,j,5] for i=1:nh, j=1:nc) - return ExaModels.ExaModel(core; kwargs...) + return core +end + +@inline function COPSBenchmark.pinene_args(::ExaModelsBackend, nh; T = Float64) + nc = 3 + ne = 5 + nm = 8 + bc = T[100, 0, 0, 0, 0] + tau = T[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] + tf = tau[nm] + h = tf / T(nh) + t = T[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(T[ + 88.35, 7.3, 2.3, 0.4, 1.75, + 76.4, 15.6, 4.5, 0.7, 2.8, + 65.1, 23.1, 5.3, 1.1, 5.8, + 50.4, 32.9, 6.0, 1.5, 9.3, + 37.5, 42.7, 6.0, 1.9, 12.0, + 25.9, 49.1, 5.9, 2.2, 17.0, + 14.0, 57.4, 5.1, 2.6, 21.0, + 4.5, 63.1, 3.8, 2.9, 25.7, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = bc[s] + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i=1:nh, s=1:ne] + uc_start = [v0[i, s] for i=1:nh, j=1:nc, s=1:ne] + itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + return (nh, (; v_start, uc_start, itr)) end + +@inline COPSBenchmark.pinene_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.pinene_recipe(b; T = T, backend = backend), + COPSBenchmark.pinene_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/rocket.jl b/ext/COPSBenchmarkExaModels/rocket.jl index 9b3fcaf..3ea3b0a 100644 --- a/ext/COPSBenchmarkExaModels/rocket.jl +++ b/ext/COPSBenchmarkExaModels/rocket.jl @@ -4,7 +4,13 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.rocket_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# The two varying starting points are data, not structure: they are +# comprehensions over the size, which has no symbolic form. They are computed +# in `rocket_args` and supplied as arguments, so the core holds no function of +# this package's. `inv_nh` does have a symbolic form and stays here. +@inline function COPSBenchmark.rocket_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) h_0 = T(1) v_0 = T(0) m_0 = T(1) @@ -14,29 +20,27 @@ v_c = T(620) m_c = T(0.6) half = T(0.5) - inv_nh = T(1) / T(nh) c = half*sqrt(g_0 * h_0) m_f = m_c * m_0 D_c = half * v_c * (m_0 / g_0) T_max = T_c * m_0 * g_0 - core = ExaModels.ExaCore(T; backend= backend, minimize=false, concrete = Val(true)) + core, nh, s_v_start, s_m_start = ExaModels.ExaCore( + T; backend = backend, minimize = false, nargs = Val(3), + ) - # Precompute generator scalars so the broadcast closure stays isbits - # (Metal rejects Type{T} captures). And bind the thrust-variable name as - # Th rather than T to avoid shadowing the Type parameter T inside the - # function — that shadowing caused Julia to Box the captured T value. - s_v_start = T[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] - s_m_start = T[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] + inv_nh = one(T) / nh s_Th_start = T_max*half + ExaModels.@add_var(core, h, 0:nh; start=h_0, lvar = h_0) ExaModels.@add_var(core, v, 0:nh; start=s_v_start, lvar = v_0) ExaModels.@add_var(core, m, 0:nh; start=s_m_start, lvar = m_f, uvar = m_0) ExaModels.@add_var(core, Th, 0:nh; start=s_Th_start, lvar = v_0, uvar = T_max) ExaModels.@add_var(core, step, 1; start=inv_nh, lvar = v_0) - ExaModels.@add_obj(core, h[nh]) + # Indexes the last point, so the index arrives as data. + ExaModels.@add_obj(core, h[k] for k in nh:nh) # Dynamics ExaModels.@add_con(core, c1, - h[i] + h[i-1] + half * step[1] * (v[i] + v[i-1]) for i=1:nh) @@ -47,7 +51,25 @@ ExaModels.@add_con(core, c4, h[0] - h_0) ExaModels.@add_con(core, c5, v[0] - v_0) ExaModels.@add_con(core, c6, m[0] - m_0) - ExaModels.@add_con(core, c7, m[nh] - m_f) + ExaModels.@add_con(core, c7, m[k] - m_f for k in nh:nh) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@inline function COPSBenchmark.rocket_args(::ExaModelsBackend, nh; T = Float64) + m_0 = T(1) + m_f = T(0.6) * m_0 + inv_nh = T(1) / T(nh) + # Precompute generator scalars so the broadcast closure stays isbits + # (Metal rejects Type{T} captures). + s_v_start = T[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] + s_m_start = T[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] + return (nh, s_v_start, s_m_start) +end + +@inline COPSBenchmark.rocket_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.rocket_recipe(b; T = T, backend = backend), + COPSBenchmark.rocket_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/steering.jl b/ext/COPSBenchmarkExaModels/steering.jl index e28f8a2..118a43c 100644 --- a/ext/COPSBenchmarkExaModels/steering.jl +++ b/ext/COPSBenchmarkExaModels/steering.jl @@ -4,31 +4,25 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.steering_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +# The state's starting point is a matrix built by a local helper over the +# discretization, so it is data and moves to `steering_args`. `u_min`/`u_max` +# and the boundary values do not depend on the size and stay here. +@inline function COPSBenchmark.steering_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) a = T(100) # Magnitude of force. # Bounds on the control u_min, u_max = -T(pi)/T(2), T(pi)/T(2) xs = zeros(T, 4) xf = [T(NaN), T(5), T(45), T(0)] half = T(0.5) - inv_nh = T(1) / T(nh) - function gen_x0(k, i) - if i == 1 || i == 4 - return T(0) - elseif i == 2 - return T(5)*T(k)*inv_nh - elseif i == 3 - return T(45)*T(k)*inv_nh - else - return T(0) - end - end + core, nh, x_start = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + inv_nh = one(T) / nh ExaModels.@add_var(core, u, 1:nh+1; lvar = u_min, uvar = u_max, start=T(0)) # control - ExaModels.@add_var(core, x, 1:nh+1, 1:4; start=[gen_x0(i, j) for i=1:nh+1, j=1:4]) # state + ExaModels.@add_var(core, x, 1:nh+1, 1:4; start=x_start) # state ExaModels.@add_var(core, tf, 1; start=T(1)) # final time ExaModels.@add_obj(core, tf[1]) @@ -39,9 +33,33 @@ ExaModels.@add_con(core, c2, -x[i+1,2] + x[i,2] + half*(tf[1]*inv_nh)*(x[i,4] + x[i+1,4]) for i=1:nh) ExaModels.@add_con(core, c3, -x[i+1,3] + x[i,3] + half*(tf[1]*inv_nh)*(a*cos(u[i]) + a*cos(u[i+1])) for i=1:nh) ExaModels.@add_con(core, c4, -x[i+1,4] + x[i,4] + half*(tf[1]*inv_nh)*(a*sin(u[i]) + a*sin(u[i+1])) for i=1:nh) - # Boundary conditions + # Boundary conditions. The final-row index depends on the size, so it comes + # from a one-element set rather than being written into the structure. ExaModels.@add_con(core, c5, -x[1, j] + s for (j,s) in enumerate(xs)) - ExaModels.@add_con(core, c6, -x[nh+1, j] + f for (j,f) in zip(2:4, xf[2:4])) + ExaModels.@add_con(core, c6, -x[k, j] + f for k in (nh+1):(nh+1), (j,f) in zip(2:4, xf[2:4])) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@inline function COPSBenchmark.steering_args(::ExaModelsBackend, nh; T = Float64) + inv_nh = T(1) / T(nh) + function gen_x0(k, i) + if i == 1 || i == 4 + return T(0) + elseif i == 2 + return T(5)*T(k)*inv_nh + elseif i == 3 + return T(45)*T(k)*inv_nh + else + return T(0) + end + end + return (nh, [gen_x0(i, j) for i=1:nh+1, j=1:4]) +end + +@inline COPSBenchmark.steering_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.steering_recipe(b; T = T, backend = backend), + COPSBenchmark.steering_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index de7823d..f51659c 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -9,24 +9,44 @@ struct ExaModelsBackend <: AbstractModelerBackend end # COPS Instances function bearing_model end function camshape_model end +function camshape_recipe end +function camshape_args end function catmix_model end +function catmix_recipe end +function catmix_args end function chain_model end +function chain_recipe end +function chain_args end function channel_model end +function channel_recipe end +function channel_args end function clnlbeam_model end function dirichlet_model end function elec_model end function gasoil_model end +function gasoil_recipe end +function gasoil_args end function glider_model end function henon_model end function lane_emden_model end function marine_model end +function marine_recipe end +function marine_args end function methanol_model end +function methanol_recipe end +function methanol_args end function minsurf_model end function pinene_model end +function pinene_recipe end +function pinene_args end function polygon_model end function robot_model end function rocket_model end +function rocket_recipe end +function rocket_args end function steering_model end +function steering_recipe end +function steering_args end function tetra_duct12_model end function tetra_duct15_model end function tetra_duct20_model end From ac4ef01ffc06afd4fe05a0729b0963d176a8cba7 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Wed, 12 Aug 2026 20:38:28 -0400 Subject: [PATCH 2/7] Write the remaining problems as recipes: robot, glider, elec, polygon, bearing, torsion, minsurf All seventeen problems the suite compares callbacks for are now split into `_recipe` (structure, size open) and `_args` (the values that close it), with `_model` the two composed. Every one is identical to the previous constructor at the sizes the suite itself uses, on dimensions, bounds, starting point, `minimize`, objective, gradient, constraints, and the Jacobian and Hessian including COO order. Notes on the ones that were not mechanical: * elec draws its starting distribution from a seed. That is data: the draw and the seed both live in `elec_args`, so the recipe holds no randomness and one structure serves every seed. * robot and glider used `@add_expr`, whose body is stored as a closure that captures the `Variable`s -- placeholders end up in the closure's type, where `instantiate` cannot reach them. Since `@add_expr` adds no variables and no constraints, only a name, a local helper splices the identical tree at each point of use. glider's eight subexpressions chain through each other exactly as before. * bearing, torsion and minsurf take two sizes, so they carry three placeholders: the two sizes and the data. * minsurf's interior block is indexed by `Int(floor(0.25/hx))`; a constructor applied to a placeholder fails, so those ranges are data -- the iteration set itself, not merely the values in it. Offloading moves values, not arithmetic. bearing first folded `e*sin(x)` into one coefficient, which re-associates `(-hx*hy*e)*sin(x)` and moved twelve gradient entries by one ulp; `sin(x)` alone is passed and `e` stays a factor in the expression. A tolerance-based comparison would not have shown it. Co-Authored-By: Claude Opus 5 --- ext/COPSBenchmarkExaModels/bearing.jl | 50 +++++++++++++------ ext/COPSBenchmarkExaModels/elec.jl | 47 +++++++++++++----- ext/COPSBenchmarkExaModels/glider.jl | 69 +++++++++++++++++---------- ext/COPSBenchmarkExaModels/minsurf.jl | 48 ++++++++++++++----- ext/COPSBenchmarkExaModels/polygon.jl | 37 ++++++++++---- ext/COPSBenchmarkExaModels/robot.jl | 59 ++++++++++++++++------- ext/COPSBenchmarkExaModels/torsion.jl | 46 ++++++++++++------ src/COPSBenchmark.jl | 14 ++++++ 8 files changed, 263 insertions(+), 107 deletions(-) diff --git a/ext/COPSBenchmarkExaModels/bearing.jl b/ext/COPSBenchmarkExaModels/bearing.jl index 5ce70cc..6350001 100644 --- a/ext/COPSBenchmarkExaModels/bearing.jl +++ b/ext/COPSBenchmarkExaModels/bearing.jl @@ -4,41 +4,61 @@ # 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, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) + 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 function COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny; T = Float64) + b = 10 + e = 0.1 + hx = 2*pi / (nx+1) + hy = 2*b / (ny+1) + 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] + # `wq` and the eccentricity term are per-row coefficients; they ride with + # the indices they belong to rather than being recomputed symbolically. + lower = [(i, j, wq(i), wq(i+1)) for i in 1:nx+1, j in 1:ny+1] + upper = [(i, j, wq(i), wq(i-1)) for i in 2:nx+2, j in 2:ny+2] + lin = [(i, j, sin((i-1)*hx)) for i in 1:nx+2, j in 1:ny+2] + return (nx, ny, (; v0, lower, upper, lin)) +end - +@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; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/elec.jl b/ext/COPSBenchmarkExaModels/elec.jl index 31bfddc..251c283 100644 --- a/ext/COPSBenchmarkExaModels/elec.jl +++ b/ext/COPSBenchmarkExaModels/elec.jl @@ -4,7 +4,29 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.elec_model(::ExaModelsBackend, np; seed = 2713, T = Float64, backend = nothing, kwargs...) +# The quasi-uniform starting distribution is a random draw, and the Coulomb +# pair list is a double comprehension over the size. Both are data: they are +# drawn and built in `elec_args`, which owns the seed, so the recipe holds no +# randomness and the structure is the same whatever was drawn. +@inline function COPSBenchmark.elec_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) + core, np, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + + ExaModels.@add_var(core, x, 1:np; start = d.x0) + ExaModels.@add_var(core, y, 1:np; start = d.y0) + ExaModels.@add_var(core, z, 1:np; start = d.z0) + + # Coulomb potential + ExaModels.@add_obj(core, 1.0 / sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2) for (i,j) in d.itr) + + # Unit-ball + ExaModels.@add_con(core, c1, x[i]^2 + y[i]^2 + z[i]^2 - 1 for i=1:np) + + return core +end + +@inline function COPSBenchmark.elec_args(::ExaModelsBackend, np; seed = 2713, T = Float64) Random.seed!(seed) # Set the starting point to a quasi-uniform distribution @@ -12,17 +34,16 @@ theta = (2pi) .* rand(np) phi = pi .* rand(np) - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) - ExaModels.@add_var(core, x, 1:np; start = [cos(theta[i])*sin(phi[i]) for i=1:np]) - ExaModels.@add_var(core, y, 1:np; start = [sin(theta[i])*sin(phi[i]) for i=1:np]) - ExaModels.@add_var(core, z, 1:np; start = [cos(phi[i]) for i=1:np]) - - # Coulomb potential + x0 = [cos(theta[i])*sin(phi[i]) for i=1:np] + y0 = [sin(theta[i])*sin(phi[i]) for i=1:np] + z0 = [cos(phi[i]) for i=1:np] itr = [(i,j) for i in 1:np-1 for j in i+1:np] - ExaModels.@add_obj(core, 1.0 / sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2) for (i,j) in itr) - - # Unit-ball - ExaModels.@add_con(core, c1, x[i]^2 + y[i]^2 + z[i]^2 - 1 for i=1:np) - - return ExaModels.ExaModel(core; kwargs...) + return (np, (; x0, y0, z0, itr)) end + +@inline COPSBenchmark.elec_model(b::ExaModelsBackend, np; seed = 2713, T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.elec_recipe(b; T = T, backend = backend), + COPSBenchmark.elec_args(b, np; seed = seed, T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/glider.jl b/ext/COPSBenchmarkExaModels/glider.jl index 9000e0c..d032ca6 100644 --- a/ext/COPSBenchmarkExaModels/glider.jl +++ b/ext/COPSBenchmarkExaModels/glider.jl @@ -5,7 +5,9 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.glider_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +@inline function COPSBenchmark.glider_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) # Design parameters (T-typed) x_0 = T(0) y_0 = T(1000) @@ -26,45 +28,64 @@ cL_max = T(1.4) cL0 = cL_max / T(2) half = T(0.5) - inv_nh = T(1) / T(nh) r_offset = T(2.5) one_T = T(1) - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) + c, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + + inv_nh = one(T) / nh ExaModels.@add_var(c, t_f, 1; lvar = T(0), start = T(1)) - ExaModels.@add_var(c, x, nh+1; lvar = zeros(T, nh+1), start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh]) - ExaModels.@add_var(c, y, nh+1; start = [y_0 + (T(k)*inv_nh)*(y_f - y_0) for k in 0:nh]) - ExaModels.@add_var(c, vx, nh+1; lvar = zeros(T, nh+1), start = fill(vx_0, nh+1)) - ExaModels.@add_var(c, vy, nh+1; start = fill(vy_0, nh+1)) - ExaModels.@add_var(c, cL, nh+1; lvar = fill(cL_min,nh+1), uvar = fill(cL_max, nh+1), start = fill(cL0, nh+1)) + ExaModels.@add_var(c, x, nh+1; lvar = T(0), start = d.x_start) + ExaModels.@add_var(c, y, nh+1; start = d.y_start) + ExaModels.@add_var(c, vx, nh+1; lvar = T(0), start = vx_0) + ExaModels.@add_var(c, vy, nh+1; start = vy_0) + ExaModels.@add_var(c, cL, nh+1; lvar = cL_min, uvar = cL_max, start = cL0) - # Expressions (matching JuMP @expressions) - ExaModels.@add_expr(c, r, (x[i]/r_0 - r_offset)^2 for i in 1:nh+1) - ExaModels.@add_expr(c, u, u_c*(one_T - r[i])*exp(-r[i]) for i in 1:nh+1) - ExaModels.@add_expr(c, w, vy[i] - u[i] for i in 1:nh+1) - ExaModels.@add_expr(c, v, sqrt(vx[i]^2 + w[i]^2) for i in 1:nh+1) - ExaModels.@add_expr(c, D, half*(cd0+cd1*cL[i]^2)*rho*S*v[i]^2 for i in 1:nh+1) - ExaModels.@add_expr(c, L, half*cL[i]*rho*S*v[i]^2 for i in 1:nh+1) - ExaModels.@add_expr(c, vx_dot, (-L[i]*(w[i]/v[i]) - D[i]*(vx[i]/v[i]))/m for i in 1:nh+1) - ExaModels.@add_expr(c, vy_dot, (L[i]*(vx[i]/v[i]) - D[i]*(w[i]/v[i]))/m - g for i in 1:nh+1) + # These were `@add_expr`s. That macro stores its body as a closure, which a + # recipe cannot instantiate; a local helper is spliced at each point of use + # exactly as the subexpression was, so the trees are unchanged. They chain, + # as the originals did -- u is written in terms of r, v in terms of w. + r(i) = (x[i]/r_0 - r_offset)^2 + u(i) = u_c*(one_T - r(i))*exp(-r(i)) + w(i) = vy[i] - u(i) + v(i) = sqrt(vx[i]^2 + w(i)^2) + D(i) = half*(cd0+cd1*cL[i]^2)*rho*S*v(i)^2 + L(i) = half*cL[i]*rho*S*v(i)^2 + vx_dot(i) = (-L(i)*(w(i)/v(i)) - D(i)*(vx[i]/v(i)))/m + vy_dot(i) = (L(i)*(vx[i]/v(i)) - D(i)*(w(i)/v(i)))/m - g - ExaModels.@add_obj(c, -x[nh+1]) + ExaModels.@add_obj(c, -x[k] for k in (nh+1):(nh+1)) # Dynamics ExaModels.@add_con(c, c1, x[j] - (x[j-1] + half * t_f[1]*inv_nh * (vx[j] + vx[j-1])) for j in 2:nh+1) ExaModels.@add_con(c, c2, y[j] - (y[j-1] + half * t_f[1]*inv_nh * (vy[j] + vy[j-1])) for j in 2:nh+1) - ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + half * t_f[1]*inv_nh * (vx_dot[j] + vx_dot[j-1])) for j in 2:nh+1) - ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + half * t_f[1]*inv_nh * (vy_dot[j] + vy_dot[j-1])) for j in 2:nh+1) + ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + half * t_f[1]*inv_nh * (vx_dot(j) + vx_dot(j-1))) for j in 2:nh+1) + ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + half * t_f[1]*inv_nh * (vy_dot(j) + vy_dot(j-1))) for j in 2:nh+1) # Boundary constraints ExaModels.@add_con(c, c5, x[1] - x_0) ExaModels.@add_con(c, c6, y[1] - y_0) - ExaModels.@add_con(c, c7, y[nh+1] - y_f) + ExaModels.@add_con(c, c7, y[k] - y_f for k in (nh+1):(nh+1)) ExaModels.@add_con(c, c8, vx[1] - vx_0) - ExaModels.@add_con(c, c9, vx[nh+1] - vx_f) + ExaModels.@add_con(c, c9, vx[k] - vx_f for k in (nh+1):(nh+1)) ExaModels.@add_con(c, c10, vy[1] - vy_0) - ExaModels.@add_con(c, c11, vy[nh+1] - vy_f) + ExaModels.@add_con(c, c11, vy[k] - vy_f for k in (nh+1):(nh+1)) + + return c +end - ExaModels.ExaModel(c; kwargs...) +@inline function COPSBenchmark.glider_args(::ExaModelsBackend, nh; T = Float64) + x_0 = T(0); y_0 = T(1000); y_f = T(900); vx_0 = T(13.23) + inv_nh = T(1) / T(nh) + x_start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh] + y_start = [y_0 + (T(k)*inv_nh)*(y_f - y_0) for k in 0:nh] + return (nh, (; x_start, y_start)) end + +@inline COPSBenchmark.glider_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.glider_recipe(b; T = T, backend = backend), + COPSBenchmark.glider_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/minsurf.jl b/ext/COPSBenchmarkExaModels/minsurf.jl index ba4b02a..3a5f88a 100644 --- a/ext/COPSBenchmarkExaModels/minsurf.jl +++ b/ext/COPSBenchmarkExaModels/minsurf.jl @@ -9,20 +9,20 @@ # Argonne National Labs Technical Report ANL/MCS-246 (2004) # classification OBR2-AN-V-V -@inline function COPSBenchmark.minsurf_model(::ExaModelsBackend, nx::Int, ny::Int; T = Float64, backend = nothing, kwargs...) - x_mesh = LinRange(0, 1, nx + 2) # coordinates of the mesh points x - - v0 = zeros(nx + 2, ny + 2) # Surface matrix initialization - for i = 1:(nx + 2), j = 1:(ny + 2) - v0[i, j] = 1 - (2 * x_mesh[i] - 1)^2 - end +# The spacings and the triangle area stay symbolic. The surface initialisation +# is a comprehension over the mesh, and the interior block's index ranges are +# `Int(floor(0.25/hx))`-style constructor calls, which a placeholder cannot take +# -- both are data. +@inline function COPSBenchmark.minsurf_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) + c, nx, ny, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) hx = 1 / (nx + 1) hy = 1 / (ny + 1) area = 1 // 2 * hx * hy - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - ExaModels.@add_var(c, v, nx+2, ny+2; start = v0) + ExaModels.@add_var(c, v, nx+2, ny+2; start = d.v0) ExaModels.@add_obj(c, area * (1 + ((v[i + 1, j] - v[i, j]) / hx)^2 + ((v[i, j + 1] - v[i, j]) / hy)^2)^(1 / 2) for i = 1:(nx + 1), j = 1:(ny + 1)) @@ -35,10 +35,12 @@ v[1, j + 1] for j in 0:ny+1 ) + # The fixed edge indices depend on the size, so they come from one-element + # sets; the enclosing set has a single value, so the row order is unchanged. ExaModels.@add_con( c, c2, - v[nx + 2, j + 1] for j in 0:ny+1 + v[k, j + 1] for k in (nx+2):(nx+2), j in 0:ny+1 ) ExaModels.@add_con( @@ -50,7 +52,7 @@ ExaModels.@add_con( c, c4, - v[i + 1, ny+2] - 1 + (2 * i * hx - 1)^2 for i in 0:nx+1 + v[i + 1, k] - 1 + (2 * i * hx - 1)^2 for i in 0:nx+1, k in (ny+2):(ny+2) ) ExaModels.@add_con( @@ -64,12 +66,32 @@ ExaModels.@add_con( c, c6, - v[i + 1, j + 1] for i in Int(floor(0.25 / hx)):Int(ceil(0.75 / hx)), j in Int(floor(0.25 / hy)):Int(ceil(0.75 / hy)); + v[i + 1, j + 1] for i in d.c6_i, j in d.c6_j; lcon = 1, ucon = Inf ) - return ExaModels.ExaModel(c; kwargs...) + return c end +@inline function COPSBenchmark.minsurf_args(::ExaModelsBackend, nx::Int, ny::Int; T = Float64) + x_mesh = LinRange(0, 1, nx + 2) # coordinates of the mesh points x + v0 = zeros(nx + 2, ny + 2) # Surface matrix initialization + for i = 1:(nx + 2), j = 1:(ny + 2) + v0[i, j] = 1 - (2 * x_mesh[i] - 1)^2 + end + + hx = 1 / (nx + 1) + hy = 1 / (ny + 1) + c6_i = Int(floor(0.25 / hx)):Int(ceil(0.75 / hx)) + c6_j = Int(floor(0.25 / hy)):Int(ceil(0.75 / hy)) + return (nx, ny, (; v0, c6_i, c6_j)) +end + +@inline COPSBenchmark.minsurf_model(b::ExaModelsBackend, nx::Int, ny::Int; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.minsurf_recipe(b; T = T, backend = backend), + COPSBenchmark.minsurf_args(b, nx, ny; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/polygon.jl b/ext/COPSBenchmarkExaModels/polygon.jl index 340ac92..8360ee2 100644 --- a/ext/COPSBenchmarkExaModels/polygon.jl +++ b/ext/COPSBenchmarkExaModels/polygon.jl @@ -6,27 +6,44 @@ # see "Benchmarking Optimization Software with COPS" # Argonne National Labs Technical Report ANL/MCS-246 (2004) -@inline function COPSBenchmark.polygon_model(::ExaModelsBackend, n::Int; T = Float64, backend = nothing, kwargs...) - N = div(n, 2) - - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) +# The public size `n` is halved to get the vertex count, which `polygon_args` +# does along with the angle start and the diameter pair list -- both +# comprehensions over that count. The two rows fixing the last vertex index it, +# so they come from one-element sets. +@inline function COPSBenchmark.polygon_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) + c, N, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) ExaModels.@add_var(c, r, N; lvar = 0.0, uvar = 1.0, start = 1.0) - ExaModels.@add_var(c, θ, N; lvar = 0.0, uvar = T(π), start = [i * π / (N - 1) - π / (N - 1) for i in 1:N]) + ExaModels.@add_var(c, θ, N; lvar = 0.0, uvar = T(π), start = d.θ0) # Objective: maximize area = 0.5 * sum(r[i]*r[i+1]*sin(θ[i+1]-θ[i])) ExaModels.@add_obj(c, -0.5 * r[i] * r[i+1] * sin(θ[i+1] - θ[i]) for i in 1:N-1) # Fix last angle and radius - ExaModels.@add_con(c, c1, θ[N] - T(π)) - ExaModels.@add_con(c, c2, r[N]) + ExaModels.@add_con(c, c1, θ[k] - T(π) for k in N:N) + ExaModels.@add_con(c, c2, r[k] for k in N:N) # Impose ordering on angles: θ[i+1] >= θ[i] ExaModels.@add_con(c, c3, θ[i+1] - θ[i] for i in 1:N-1; lcon = 0.0, ucon = Inf) # Diameter constraint: r[i]^2 + r[j]^2 - 2*r[i]*r[j]*cos(θ[i]-θ[j]) <= 1 - pairs = [(i, j) for i in 1:N-1 for j in i+1:N] - ExaModels.@add_con(c, c4, r[i]^2 + r[j]^2 - 2*r[i]*r[j]*cos(θ[i] - θ[j]) - 1 for (i, j) in pairs; lcon = -Inf, ucon = 0.0) + ExaModels.@add_con(c, c4, r[i]^2 + r[j]^2 - 2*r[i]*r[j]*cos(θ[i] - θ[j]) - 1 for (i, j) in d.pairs; lcon = -Inf, ucon = 0.0) + + return c +end - return ExaModels.ExaModel(c; kwargs...) +@inline function COPSBenchmark.polygon_args(::ExaModelsBackend, n::Int; T = Float64) + N = div(n, 2) + θ0 = [i * π / (N - 1) - π / (N - 1) for i in 1:N] + pairs = [(i, j) for i in 1:N-1 for j in i+1:N] + return (N, (; θ0, pairs)) end + +@inline COPSBenchmark.polygon_model(b::ExaModelsBackend, n::Int; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.polygon_recipe(b; T = T, backend = backend), + COPSBenchmark.polygon_args(b, n; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/robot.jl b/ext/COPSBenchmarkExaModels/robot.jl index dcbf135..8a3ee08 100644 --- a/ext/COPSBenchmarkExaModels/robot.jl +++ b/ext/COPSBenchmarkExaModels/robot.jl @@ -5,7 +5,9 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.robot_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) +@inline function COPSBenchmark.robot_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) # total length of arm L = T(5) @@ -20,19 +22,20 @@ phi0 = pi_T / T(4) half = T(0.5) third = T(1) / T(3) - inv_nh = T(1) / T(nh) two_pi_3 = T(2) * pi_T / T(3) four_pi_3 = T(4) * pi_T / T(3) zero_T = T(0) - core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + + inv_nh = one(T) / nh ExaModels.@add_var(core, rho, nh+1; start=rho0, lvar = zero_T, uvar = L) - ExaModels.@add_var(core, the, nh+1; start=[two_pi_3*(T(k)*inv_nh)^2 for k=1:nh+1], lvar = -pi_T, uvar = pi_T) + ExaModels.@add_var(core, the, nh+1; start=d.the_start, lvar = -pi_T, uvar = pi_T) ExaModels.@add_var(core, phi, nh+1; start=phi0, lvar = zero_T, uvar = pi_T) # Derivatives ExaModels.@add_var(core, rho_dot, nh+1; start=zero_T) - ExaModels.@add_var(core, the_dot, nh+1; start=[four_pi_3*(T(k)*inv_nh) for k=1:nh+1]) + ExaModels.@add_var(core, the_dot, nh+1; start=d.the_dot_start) ExaModels.@add_var(core, phi_dot, nh+1; start=zero_T) # Control ExaModels.@add_var(core, u_rho, nh+1; start=zero_T, lvar = -max_u_rho, uvar = max_u_rho) @@ -41,9 +44,14 @@ # Final time ExaModels.@add_var(core, tf, 1; start = T(1), lvar = zero_T) - # Expressions (matching JuMP @expressions) - ExaModels.@add_expr(core, I_the, ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2*third for i=1:nh+1) - ExaModels.@add_expr(core, I_phi, ((L-rho[i])^3+rho[i]^3)*third for i=1:nh+1) + # The two moments of inertia were `@add_expr`s. That macro stores its body + # as a closure, and in a recipe the closure captures `Variable`s whose types + # carry the placeholder, which `instantiate` cannot reach into. A local + # helper is spliced at the point of use exactly as the subexpression was -- + # `@add_expr` adds no variables and no constraints, only a name -- so the + # expression trees, and the model, are unchanged. + I_the(i) = ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2*third + I_phi(i) = ((L-rho[i])^3+rho[i]^3)*third ExaModels.@add_obj(core, tf[1] for i=1:1) @@ -52,22 +60,39 @@ ExaModels.@add_con(core, c2, - phi[j] + phi[j-1] + half * tf[1]*inv_nh * (phi_dot[j] + phi_dot[j-1]) for j=2:nh+1) ExaModels.@add_con(core, c3, - the[j] + the[j-1] + half * tf[1]*inv_nh * (the_dot[j] + the_dot[j-1]) for j=2:nh+1) ExaModels.@add_con(core, c4, - rho_dot[j] + rho_dot[j-1] + half * tf[1]*inv_nh * (u_rho[j] + u_rho[j-1]) / L for j=2:nh+1) - ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + half * tf[1]*inv_nh * (u_the[j] / I_the[j] + u_the[j-1] / I_the[j-1]) for j=2:nh+1) - ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + half * tf[1]*inv_nh * (u_phi[j] / I_phi[j] + u_phi[j-1] / I_phi[j-1]) for j=2:nh+1) + ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + half * tf[1]*inv_nh * (u_the[j] / I_the(j) + u_the[j-1] / I_the(j-1)) for j=2:nh+1) + ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + half * tf[1]*inv_nh * (u_phi[j] / I_phi(j) + u_phi[j-1] / I_phi(j-1)) for j=2:nh+1) # Boundary conditions ExaModels.@add_con(core, c7, - rho[1] + rho0) ExaModels.@add_con(core, c8, - the[1] + zero_T) ExaModels.@add_con(core, c9, - phi[1] + phi0) - ExaModels.@add_con(core, c10, - rho[nh+1] + rho0) - ExaModels.@add_con(core, c11, - the[nh+1] + two_pi_3) - ExaModels.@add_con(core, c12, - phi[nh+1] + phi0) + ExaModels.@add_con(core, c10, - rho[k] + rho0 for k in (nh+1):(nh+1)) + ExaModels.@add_con(core, c11, - the[k] + two_pi_3 for k in (nh+1):(nh+1)) + ExaModels.@add_con(core, c12, - phi[k] + phi0 for k in (nh+1):(nh+1)) ExaModels.@add_con(core, c13, - rho_dot[1] + zero_T) ExaModels.@add_con(core, c14, - the_dot[1] + zero_T) ExaModels.@add_con(core, c15, - phi_dot[1] + zero_T) - ExaModels.@add_con(core, c16, - rho_dot[nh+1] + zero_T) - ExaModels.@add_con(core, c17, - the_dot[nh+1] + zero_T) - ExaModels.@add_con(core, c18, - phi_dot[nh+1] + zero_T) + ExaModels.@add_con(core, c16, - rho_dot[k] + zero_T for k in (nh+1):(nh+1)) + ExaModels.@add_con(core, c17, - the_dot[k] + zero_T for k in (nh+1):(nh+1)) + ExaModels.@add_con(core, c18, - phi_dot[k] + zero_T for k in (nh+1):(nh+1)) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@inline function COPSBenchmark.robot_args(::ExaModelsBackend, nh; T = Float64) + pi_T = T(pi) + inv_nh = T(1) / T(nh) + two_pi_3 = T(2) * pi_T / T(3) + four_pi_3 = T(4) * pi_T / T(3) + the_start = [two_pi_3*(T(k)*inv_nh)^2 for k=1:nh+1] + the_dot_start = [four_pi_3*(T(k)*inv_nh) for k=1:nh+1] + return (nh, (; the_start, the_dot_start)) +end + +@inline COPSBenchmark.robot_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.robot_recipe(b; T = T, backend = backend), + COPSBenchmark.robot_args(b, nh; T = T)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/torsion.jl b/ext/COPSBenchmarkExaModels/torsion.jl index 56a7bf8..13a671f 100644 --- a/ext/COPSBenchmarkExaModels/torsion.jl +++ b/ext/COPSBenchmarkExaModels/torsion.jl @@ -3,35 +3,51 @@ # Version 2.0 - October 2000 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.torsion_model(::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) +# Two sizes again. The spacings and the triangle area keep their symbolic form; +# the distance-to-boundary field, the row list and the bound vectors built from +# it are comprehensions over the grid and travel as data. +@inline function COPSBenchmark.torsion_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) c_val = T(5.0) - hx = T(1.0 / (nx + 1.0)) - hy = T(1.0 / (ny + 1.0)) - area = T(0.5) * hx * hy - # Distance to boundary: D[k1,k2] for k1 in 1:nx+2, k2 in 1:ny+2 - D = T[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] + core, nx, ny, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) - core = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - ExaModels.@add_var(core, v, nx+2, ny+2; start = D) + hx = one(T) / (nx + 1) + hy = one(T) / (ny + 1) + area = T(0.5) * hx * hy + + ExaModels.@add_var(core, v, nx+2, ny+2; start = d.D) # Objective = area * ((quadLower + quadUpper)/2 - c*(linLower + linUpper)/3) - # quadLower: i in 0:nx, j in 0:ny → k1 in 1:nx+1, k2 in 1:ny+1 ExaModels.@add_obj(core, area / 2 * (((v[k1+1,k2] - v[k1,k2])/hx)^2 + ((v[k1,k2+1] - v[k1,k2])/hy)^2) for k1 in 1:nx+1, k2 in 1:ny+1) - # quadUpper: i in 1:nx+1, j in 1:ny+1 → k1 in 2:nx+2, k2 in 2:ny+2 ExaModels.@add_obj(core, area / 2 * (((v[k1,k2] - v[k1-1,k2])/hx)^2 + ((v[k1,k2] - v[k1,k2-1])/hy)^2) for k1 in 2:nx+2, k2 in 2:ny+2) - # linLower: i in 0:nx, j in 0:ny → k1 in 1:nx+1, k2 in 1:ny+1 ExaModels.@add_obj(core, -area * c_val / 3 * (v[k1+1,k2] + v[k1,k2] + v[k1,k2+1]) for k1 in 1:nx+1, k2 in 1:ny+1) - # linUpper: i in 1:nx+1, j in 1:ny+1 → k1 in 2:nx+2, k2 in 2:ny+2 ExaModels.@add_obj(core, -area * c_val / 3 * (v[k1,k2] + v[k1-1,k2] + v[k1,k2-1]) for k1 in 2:nx+2, k2 in 2:ny+2) # Bound constraints on v: -D <= v <= D (matching JuMP's @constraint formulation) - D_flat = [(k1, k2, D[k1, k2]) for k1 in 1:nx+2, k2 in 1:ny+2] - ExaModels.@add_con(core, c1, v[k1, k2] for (k1, k2, d) in D_flat; lcon = [-d for (_, _, d) in D_flat], ucon = [d for (_, _, d) in D_flat]) + ExaModels.@add_con(core, c1, v[k1, k2] for (k1, k2, dd) in d.D_flat; lcon = d.lcon, ucon = d.ucon) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@inline function COPSBenchmark.torsion_args(::ExaModelsBackend, nx, ny; T = Float64) + hx = T(1.0 / (nx + 1.0)) + hy = T(1.0 / (ny + 1.0)) + D = T[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] + D_flat = [(k1, k2, D[k1, k2]) for k1 in 1:nx+2, k2 in 1:ny+2] + lcon = [-d for (_, _, d) in D_flat] + ucon = [d for (_, _, d) in D_flat] + return (nx, ny, (; D, D_flat, lcon, ucon)) +end + +@inline COPSBenchmark.torsion_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.torsion_recipe(b; T = T, backend = backend), + COPSBenchmark.torsion_args(b, nx, ny; T = T)...; + kwargs..., + ) diff --git a/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index f51659c..dc82378 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -8,6 +8,8 @@ struct ExaModelsBackend <: AbstractModelerBackend end # COPS Instances function bearing_model end +function bearing_recipe end +function bearing_args end function camshape_model end function camshape_recipe end function camshape_args end @@ -23,10 +25,14 @@ function channel_args end function clnlbeam_model end function dirichlet_model end function elec_model end +function elec_recipe end +function elec_args end function gasoil_model end function gasoil_recipe end function gasoil_args end function glider_model end +function glider_recipe end +function glider_args end function henon_model end function lane_emden_model end function marine_model end @@ -36,11 +42,17 @@ function methanol_model end function methanol_recipe end function methanol_args end function minsurf_model end +function minsurf_recipe end +function minsurf_args end function pinene_model end function pinene_recipe end function pinene_args end function polygon_model end +function polygon_recipe end +function polygon_args end function robot_model end +function robot_recipe end +function robot_args end function rocket_model end function rocket_recipe end function rocket_args end @@ -54,6 +66,8 @@ function tetra_foam5_model end function tetra_gear_model end function tetra_hook_model end function torsion_model end +function torsion_recipe end +function torsion_args end function triangle_deer_model end function triangle_pacman_model end function triangle_turtle_model end From c7ec648e0edac0c318d0fd3d82a6a3f713466d70 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Thu, 13 Aug 2026 06:53:23 -0400 Subject: [PATCH 3/7] Make every model AOT-compilable: derive data inside the recipe Moving derived data into `_args` kept the models identical and made all but one of them uncompilable: `P_new` instantiates from integers, so a recipe whose arguments carry tables is refused. The data now stays where the recipe can reach it -- a deferred call to a named package function, `d = ArgNode1(COPSBenchmark.foo_data, nh)` (ArgNode2 for the two-size models), with `_args` returning only the sizes. `_recipe`/`_args`/`_model` keep their meaning; only where the data is computed moves. Because those functions now run inside the compiled library, they inherit `--trim=safe`'s constraints, and two idioms that are fine under the JIT do not survive it: * `T[...]` where `T` is the keyword argument is a runtime-typed comprehension -- `_array_for(T::Type, ...)` cannot resolve. Twenty sites now use a concrete element type; the core converts on append. * elec's seed is baked via `Base.Fix2(elec_data, seed)`, keeping every type in the core named and making the library deterministic per seed. robot and glider return to `@add_expr`, which stores a node rather than a closure as of the ExaModels change this branch builds against. All 20 harness cases remain identical to main's constructors. Compiled: each model individually, then all seventeen into one shared library (16.2 MB, 0 verifier errors), every model exact against its in-Julia counterpart at sizes the library was never compiled at. Co-Authored-By: Claude Opus 5 --- ext/COPSBenchmarkExaModels/bearing.jl | 18 +- ext/COPSBenchmarkExaModels/catmix.jl | 11 +- ext/COPSBenchmarkExaModels/chain.jl | 16 +- ext/COPSBenchmarkExaModels/channel.jl | 74 +---- ext/COPSBenchmarkExaModels/elec.jl | 22 +- ext/COPSBenchmarkExaModels/gasoil.jl | 58 +--- ext/COPSBenchmarkExaModels/glider.jl | 36 +-- ext/COPSBenchmarkExaModels/marine.jl | 52 +--- ext/COPSBenchmarkExaModels/methanol.jl | 39 +-- ext/COPSBenchmarkExaModels/minsurf.jl | 18 +- ext/COPSBenchmarkExaModels/pinene.jl | 42 +-- ext/COPSBenchmarkExaModels/polygon.jl | 10 +- ext/COPSBenchmarkExaModels/robot.jl | 28 +- ext/COPSBenchmarkExaModels/rocket.jl | 20 +- ext/COPSBenchmarkExaModels/steering.jl | 21 +- ext/COPSBenchmarkExaModels/torsion.jl | 13 +- src/COPSBenchmark.jl | 381 +++++++++++++++++++++++++ 17 files changed, 451 insertions(+), 408 deletions(-) diff --git a/ext/COPSBenchmarkExaModels/bearing.jl b/ext/COPSBenchmarkExaModels/bearing.jl index 6350001..9325e57 100644 --- a/ext/COPSBenchmarkExaModels/bearing.jl +++ b/ext/COPSBenchmarkExaModels/bearing.jl @@ -13,7 +13,8 @@ b = 10 # grid is (0,2*pi)x(0,2*b) e = 0.1 # eccentricity - core, nx, ny, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) + 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 @@ -41,20 +42,7 @@ return core end -@inline function COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny; T = Float64) - b = 10 - e = 0.1 - hx = 2*pi / (nx+1) - hy = 2*b / (ny+1) - 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] - # `wq` and the eccentricity term are per-row coefficients; they ride with - # the indices they belong to rather than being recomputed symbolically. - lower = [(i, j, wq(i), wq(i+1)) for i in 1:nx+1, j in 1:ny+1] - upper = [(i, j, wq(i), wq(i-1)) for i in 2:nx+2, j in 2:ny+2] - lin = [(i, j, sin((i-1)*hx)) for i in 1:nx+2, j in 1:ny+2] - return (nx, ny, (; v0, lower, upper, lin)) -end +@inline COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny; T = Float64) = (nx, ny) @inline COPSBenchmark.bearing_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/catmix.jl b/ext/COPSBenchmarkExaModels/catmix.jl index bbe1b8a..ca2d352 100644 --- a/ext/COPSBenchmarkExaModels/catmix.jl +++ b/ext/COPSBenchmarkExaModels/catmix.jl @@ -24,7 +24,8 @@ one_T = T(1) rho_index = [(i, rho[i]) for i in 1:nc] - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.catmix_data, nh) h = one(T) / nh # Final time / nh @@ -83,13 +84,7 @@ return core end -@inline function COPSBenchmark.catmix_args(::ExaModelsBackend, nh; T = Float64) - ne = 2 - nc = 3 - v_start = T[mod(j, ne) for i in 1:nh, j in 1:ne] - pp_start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne] - return (nh, (; v_start, pp_start)) -end +@inline COPSBenchmark.catmix_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.catmix_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/chain.jl b/ext/COPSBenchmarkExaModels/chain.jl index c5cf307..5893816 100644 --- a/ext/COPSBenchmarkExaModels/chain.jl +++ b/ext/COPSBenchmarkExaModels/chain.jl @@ -19,7 +19,8 @@ b = 3 tf = 1.0 - c, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.chain_data, nh) h = tf / nh @@ -82,18 +83,7 @@ return c end -@inline function COPSBenchmark.chain_args(::ExaModelsBackend, n; T = Float64) - nh = max(2, div(n - 4, 4)) - a = 1 - b = 3 - tmin = b > a ? 1 / 4 : 3 / 4 - u0 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] - x10 = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1] - x20 = [(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] - x30 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] - return (nh, (; u0, x10, x20, x30)) -end +@inline COPSBenchmark.chain_args(::ExaModelsBackend, n; T = Float64) = (max(2, div(n - 4, 4)),) @inline COPSBenchmark.chain_model(b::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/channel.jl b/ext/COPSBenchmarkExaModels/channel.jl index 14799db..85d2a28 100644 --- a/ext/COPSBenchmarkExaModels/channel.jl +++ b/ext/COPSBenchmarkExaModels/channel.jl @@ -14,7 +14,8 @@ ) bc = T[0.0 1.0; 0.0 0.0] - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.channel_data, nh) ExaModels.@add_var(core, v, nh, 4; start = d.v0) ExaModels.@add_var(core, w, nh, 4; start = 0.0) @@ -71,76 +72,7 @@ return core end -@inline function COPSBenchmark.channel_args(::ExaModelsBackend, nh; T = Float64) - nc = 4 - nd = 4 - tf = T(1.0) - h = tf / nh - - rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] - t = T[(i-1)*h for i in 1:nh+1] - - # Initial value - v0 = zeros(T, nh, nd) - for i in 1:nh - v0[i, 1] = t[i]^2*(3 - 2*t[i]) - v0[i, 2] = 6*t[i]*(1 - t[i]) - v0[i, 3] = 6*(1 - 2*t[i]) - v0[i, 4] = -12 - end - uc0 = T[v0[i, s] for i in 1:nh, j in 1:nc, s in 1:nd] - - # fac[k+1] = k! - fac = T[factorial(k) for k in 0:nc+nd] - - con1_itr = [ - (i, j, s, h*rho[j]^1/fac[2], h*rho[j]^2/fac[3], h*rho[j]^3/fac[4], h*rho[j]^4/fac[5]) - for i in 1:nh, j in 1:nc, s in 1:nd - ] - - con2_itr = [ - (i, j, s, - (s <= 1 ? (rho[j]*h)^(1-s)/fac[1-s+1] : zero(T)), - (s <= 2 ? (rho[j]*h)^(2-s)/fac[2-s+1] : zero(T)), - (s <= 3 ? (rho[j]*h)^(3-s)/fac[3-s+1] : zero(T)), - (s <= 4 ? (rho[j]*h)^(4-s)/fac[4-s+1] : zero(T)), - h^(nd-s+1)*rho[j]^(1+nd-s)/fac[1+nd-s+1], - h^(nd-s+1)*rho[j]^(2+nd-s)/fac[2+nd-s+1], - h^(nd-s+1)*rho[j]^(3+nd-s)/fac[3+nd-s+1], - h^(nd-s+1)*rho[j]^(4+nd-s)/fac[4+nd-s+1]) - for i in 1:nh, j in 1:nc, s in 1:nd - ] - - cont_itr = [ - (i, s, - (s <= 1 ? h^(1-s)/fac[1-s+1] : zero(T)), - (s <= 2 ? h^(2-s)/fac[2-s+1] : zero(T)), - (s <= 3 ? h^(3-s)/fac[3-s+1] : zero(T)), - (s <= 4 ? h^(4-s)/fac[4-s+1] : zero(T)), - h^(nd-s+1)/fac[1+nd-s+1], - h^(nd-s+1)/fac[2+nd-s+1], - h^(nd-s+1)/fac[3+nd-s+1], - h^(nd-s+1)/fac[4+nd-s+1]) - for i in 1:nh-1, s in 1:nd - ] - - coll_itr = [ - (i, j, rho[j]^0/fac[1], rho[j]^1/fac[2], rho[j]^2/fac[3], rho[j]^3/fac[4]) - for i in 1:nh, j in 1:nc - ] - - # right BC coefficients, carried with the row index they apply to - bc3_cv = [h^(k-1)/fac[k] for k in 1:nd] - bc3_cw = [h^nd/fac[k+nd] for k in 1:nc] - bc4_cv = [h^(k-2)/fac[k-1] for k in 2:nd] - bc4_cw = [h^(nd-1)/fac[k+nd-1] for k in 1:nc] - bc3_itr = [(nh, bc3_cv[1], bc3_cv[2], bc3_cv[3], bc3_cv[4], - bc3_cw[1], bc3_cw[2], bc3_cw[3], bc3_cw[4])] - bc4_itr = [(nh, bc4_cv[1], bc4_cv[2], bc4_cv[3], - bc4_cw[1], bc4_cw[2], bc4_cw[3], bc4_cw[4])] - - return (nh, (; v0, uc0, con1_itr, con2_itr, cont_itr, coll_itr, bc3_itr, bc4_itr)) -end +@inline COPSBenchmark.channel_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.channel_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/elec.jl b/ext/COPSBenchmarkExaModels/elec.jl index 251c283..5d26614 100644 --- a/ext/COPSBenchmarkExaModels/elec.jl +++ b/ext/COPSBenchmarkExaModels/elec.jl @@ -9,9 +9,10 @@ # drawn and built in `elec_args`, which owns the seed, so the recipe holds no # randomness and the structure is the same whatever was drawn. @inline function COPSBenchmark.elec_recipe( - ::ExaModelsBackend; T = Float64, backend = nothing, + ::ExaModelsBackend; seed = 2713, T = Float64, backend = nothing, ) - core, np, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, np = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(Base.Fix2(COPSBenchmark.elec_data, seed), np) ExaModels.@add_var(core, x, 1:np; start = d.x0) ExaModels.@add_var(core, y, 1:np; start = d.y0) @@ -26,24 +27,11 @@ return core end -@inline function COPSBenchmark.elec_args(::ExaModelsBackend, np; seed = 2713, T = Float64) - Random.seed!(seed) - - # Set the starting point to a quasi-uniform distribution - # of electrons on a unit sphere - theta = (2pi) .* rand(np) - phi = pi .* rand(np) - - x0 = [cos(theta[i])*sin(phi[i]) for i=1:np] - y0 = [sin(theta[i])*sin(phi[i]) for i=1:np] - z0 = [cos(phi[i]) for i=1:np] - itr = [(i,j) for i in 1:np-1 for j in i+1:np] - return (np, (; x0, y0, z0, itr)) -end +@inline COPSBenchmark.elec_args(::ExaModelsBackend, np; seed = 2713, T = Float64) = (np,) @inline COPSBenchmark.elec_model(b::ExaModelsBackend, np; seed = 2713, T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( - COPSBenchmark.elec_recipe(b; T = T, backend = backend), + COPSBenchmark.elec_recipe(b; seed = seed, T = T, backend = backend), COPSBenchmark.elec_args(b, np; seed = seed, T = T)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/gasoil.jl b/ext/COPSBenchmarkExaModels/gasoil.jl index dc2a3f1..a6929a5 100644 --- a/ext/COPSBenchmarkExaModels/gasoil.jl +++ b/ext/COPSBenchmarkExaModels/gasoil.jl @@ -22,7 +22,8 @@ zero_T = T(0) z1 = T[1.0000, 0.0000] - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.gasoil_data, nh) h = tf / nh @@ -77,60 +78,7 @@ return core end -@inline function COPSBenchmark.gasoil_args(::ExaModelsBackend, nh; T = Float64) - nc = 4 - ne = 2 - nm = 21 - rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] - bc = [1, 1, 2, 0] - tau = T[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] - tf = tau[nm] - h = tf / T(nh) - t = T[T(i-1)*h for i in 1:nh+1] - - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - - z = reshape(T[ - 1.0000, 0.0000, - 0.8105, 0.2000, - 0.6208, 0.2886, - 0.5258, 0.3010, - 0.4345, 0.3215, - 0.3903, 0.3123, - 0.3342, 0.2716, - 0.3034, 0.2551, - 0.2735, 0.2258, - 0.2405, 0.1959, - 0.2283, 0.1789, - 0.2071, 0.1457, - 0.1669, 0.1198, - 0.1530, 0.0909, - 0.1339, 0.0719, - 0.1265, 0.0561, - 0.1200, 0.0460, - 0.0990, 0.0280, - 0.0870, 0.0190, - 0.0770, 0.0140, - 0.0690, 0.0100, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - for i in 1:itau[1], s in 1:ne - v0[i, s] = T(bc[s]) - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - - v_start = [v0[i, s] for i = 1:nh, s = 1:ne] - uc_start = [v0[i, s] for i = 1:nh, j = 1:nc, s = 1:ne] - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] - - return (nh, (; v_start, uc_start, itr)) -end +@inline COPSBenchmark.gasoil_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.gasoil_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/glider.jl b/ext/COPSBenchmarkExaModels/glider.jl index d032ca6..129e477 100644 --- a/ext/COPSBenchmarkExaModels/glider.jl +++ b/ext/COPSBenchmarkExaModels/glider.jl @@ -31,7 +31,8 @@ r_offset = T(2.5) one_T = T(1) - c, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.glider_data, nh) inv_nh = one(T) / nh @@ -42,26 +43,23 @@ ExaModels.@add_var(c, vy, nh+1; start = vy_0) ExaModels.@add_var(c, cL, nh+1; lvar = cL_min, uvar = cL_max, start = cL0) - # These were `@add_expr`s. That macro stores its body as a closure, which a - # recipe cannot instantiate; a local helper is spliced at each point of use - # exactly as the subexpression was, so the trees are unchanged. They chain, - # as the originals did -- u is written in terms of r, v in terms of w. - r(i) = (x[i]/r_0 - r_offset)^2 - u(i) = u_c*(one_T - r(i))*exp(-r(i)) - w(i) = vy[i] - u(i) - v(i) = sqrt(vx[i]^2 + w(i)^2) - D(i) = half*(cd0+cd1*cL[i]^2)*rho*S*v(i)^2 - L(i) = half*cL[i]*rho*S*v(i)^2 - vx_dot(i) = (-L(i)*(w(i)/v(i)) - D(i)*(vx[i]/v(i)))/m - vy_dot(i) = (L(i)*(vx[i]/v(i)) - D(i)*(w(i)/v(i)))/m - g + # Expressions (matching JuMP @expressions) + ExaModels.@add_expr(c, r, (x[i]/r_0 - r_offset)^2 for i in 1:nh+1) + ExaModels.@add_expr(c, u, u_c*(one_T - r[i])*exp(-r[i]) for i in 1:nh+1) + ExaModels.@add_expr(c, w, vy[i] - u[i] for i in 1:nh+1) + ExaModels.@add_expr(c, v, sqrt(vx[i]^2 + w[i]^2) for i in 1:nh+1) + ExaModels.@add_expr(c, D, half*(cd0+cd1*cL[i]^2)*rho*S*v[i]^2 for i in 1:nh+1) + ExaModels.@add_expr(c, L, half*cL[i]*rho*S*v[i]^2 for i in 1:nh+1) + ExaModels.@add_expr(c, vx_dot, (-L[i]*(w[i]/v[i]) - D[i]*(vx[i]/v[i]))/m for i in 1:nh+1) + ExaModels.@add_expr(c, vy_dot, (L[i]*(vx[i]/v[i]) - D[i]*(w[i]/v[i]))/m - g for i in 1:nh+1) ExaModels.@add_obj(c, -x[k] for k in (nh+1):(nh+1)) # Dynamics ExaModels.@add_con(c, c1, x[j] - (x[j-1] + half * t_f[1]*inv_nh * (vx[j] + vx[j-1])) for j in 2:nh+1) ExaModels.@add_con(c, c2, y[j] - (y[j-1] + half * t_f[1]*inv_nh * (vy[j] + vy[j-1])) for j in 2:nh+1) - ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + half * t_f[1]*inv_nh * (vx_dot(j) + vx_dot(j-1))) for j in 2:nh+1) - ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + half * t_f[1]*inv_nh * (vy_dot(j) + vy_dot(j-1))) for j in 2:nh+1) + ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + half * t_f[1]*inv_nh * (vx_dot[j] + vx_dot[j-1])) for j in 2:nh+1) + ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + half * t_f[1]*inv_nh * (vy_dot[j] + vy_dot[j-1])) for j in 2:nh+1) # Boundary constraints ExaModels.@add_con(c, c5, x[1] - x_0) @@ -75,13 +73,7 @@ return c end -@inline function COPSBenchmark.glider_args(::ExaModelsBackend, nh; T = Float64) - x_0 = T(0); y_0 = T(1000); y_f = T(900); vx_0 = T(13.23) - inv_nh = T(1) / T(nh) - x_start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh] - y_start = [y_0 + (T(k)*inv_nh)*(y_f - y_0) for k in 0:nh] - return (nh, (; x_start, y_start)) -end +@inline COPSBenchmark.glider_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.glider_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/marine.jl b/ext/COPSBenchmarkExaModels/marine.jl index 80a19c0..4f85811 100644 --- a/ext/COPSBenchmarkExaModels/marine.jl +++ b/ext/COPSBenchmarkExaModels/marine.jl @@ -19,7 +19,8 @@ tf = T(10) # tau[nm] zero_T = T(0) - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.marine_data, nh) h = tf / nh # uniform interval length @@ -80,54 +81,7 @@ return core end -@inline function COPSBenchmark.marine_args(::ExaModelsBackend, nh; T = Float64) - ne = 8 - nm = 21 - tau = collect(T, range(T(0), T(10), 21)) - tf = tau[nm] - h = tf / T(nh) - t = T[T(i-1)*h for i in 1:nh+1] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - - z = reshape(T[ - 20000.0, 17000.0, 10000.0, 15000.0, 12000.0, 9000.0, 7000.0, 3000.0, - 12445.0, 15411.0, 13040.0, 13338.0, 13484.0, 8426.0, 6615.0, 4022.0, - 7705.0, 13074.0, 14623.0, 11976.0, 12453.0, 9272.0, 6891.0, 5020.0, - 4664.0, 8579.0, 12434.0, 12603.0, 11738.0, 9710.0, 6821.0, 5722.0, - 2977.0, 7053.0, 11219.0, 11340.0, 13665.0, 8534.0, 6242.0, 5695.0, - 1769.0, 5054.0, 10065.0, 11232.0, 12112.0, 9600.0, 6647.0, 7034.0, - 943.0, 3907.0, 9473.0, 10334.0, 11115.0, 8826.0, 6842.0, 7348.0, - 581.0, 2624.0, 7421.0, 10297.0, 12427.0, 8747.0, 7199.0, 7684.0, - 355.0, 1744.0, 5369.0, 7748.0, 10057.0, 8698.0, 6542.0, 7410.0, - 223.0, 1272.0, 4713.0, 6869.0, 9564.0, 8766.0, 6810.0, 6961.0, - 137.0, 821.0, 3451.0, 6050.0, 8671.0, 8291.0, 6827.0, 7525.0, - 87.0, 577.0, 2649.0, 5454.0, 8430.0, 7411.0, 6423.0, 8388.0, - 49.0, 337.0, 2058.0, 4115.0, 7435.0, 7627.0, 6268.0, 7189.0, - 32.0, 228.0, 1440.0, 3790.0, 6474.0, 6658.0, 5859.0, 7467.0, - 17.0, 168.0, 1178.0, 3087.0, 6524.0, 5880.0, 5562.0, 7144.0, - 11.0, 99.0, 919.0, 2596.0, 5360.0, 5762.0, 4480.0, 7256.0, - 7.0, 65.0, 647.0, 1873.0, 4556.0, 5058.0, 4944.0, 7538.0, - 4.0, 44.0, 509.0, 1571.0, 4009.0, 4527.0, 4233.0, 6649.0, - 2.0, 27.0, 345.0, 1227.0, 3677.0, 4229.0, 3805.0, 6378.0, - 1.0, 20.0, 231.0, 934.0, 3197.0, 3695.0, 3159.0, 6454.0, - 1.0, 12.0, 198.0, 707.0, 2562.0, 3163.0, 3232.0, 5566.0, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - for i in 1:itau[1], s in 1:ne - v0[i, s] = z[1, s] - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - - v_start = [v0[i, s] for i=1:nh, s=1:ne] - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] - return (nh, (; v_start, itr)) -end +@inline COPSBenchmark.marine_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.marine_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/methanol.jl b/ext/COPSBenchmarkExaModels/methanol.jl index e44cb52..d6a7503 100644 --- a/ext/COPSBenchmarkExaModels/methanol.jl +++ b/ext/COPSBenchmarkExaModels/methanol.jl @@ -21,7 +21,8 @@ two_T = T(2) bc = T[1, 0, 0] - c, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.methanol_data, nh) h = tf / nh # uniform interval length @@ -82,41 +83,7 @@ return c end -@inline function COPSBenchmark.methanol_args(::ExaModelsBackend, nh; T = Float64) - ne = 3 - nm = 17 - tau = T[ - 0., 0.050, 0.065, 0.080, 0.123, 0.233, 0.273, 0.354, 0.397, 0.418, - 0.502, 0.553, 0.681, 0.750, 0.916, 0.937, 1.122, - ] - tf = tau[nm] - h = tf / T(nh) - t = T[T(i-1)*h for i in 1:nh+1] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - - z = reshape(T[ - 1.0000, 0.0000, 0.0000, - 0.7085, 0.1621, 0.0811, - 0.5971, 0.1855, 0.0965, - 0.5537, 0.1989, 0.1198, - 0.3684, 0.2845, 0.1535, - 0.1712, 0.3491, 0.2097, - 0.1198, 0.3098, 0.2628, - 0.0747, 0.3576, 0.2467, - 0.0529, 0.3347, 0.2884, - 0.0415, 0.3388, 0.2757, - 0.0261, 0.3557, 0.3167, - 0.0208, 0.3483, 0.2954, - 0.0085, 0.3836, 0.2950, - 0.0053, 0.3611, 0.2937, - 0.0019, 0.3609, 0.2831, - 0.0018, 0.3485, 0.2846, - 0.0006, 0.3698, 0.2899, - ], ne, nm)' - - con1_matrix = [(j, s, itau[j], tau[j], z[j,s], t[itau[j]]) for j in 1:nm, s in 1:ne] - return (nh, (; con1_matrix)) -end +@inline COPSBenchmark.methanol_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.methanol_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/minsurf.jl b/ext/COPSBenchmarkExaModels/minsurf.jl index 3a5f88a..865a0d3 100644 --- a/ext/COPSBenchmarkExaModels/minsurf.jl +++ b/ext/COPSBenchmarkExaModels/minsurf.jl @@ -16,7 +16,8 @@ @inline function COPSBenchmark.minsurf_recipe( ::ExaModelsBackend; T = Float64, backend = nothing, ) - c, nx, ny, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) + c, nx, ny = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + d = ExaModels.ArgNode2(COPSBenchmark.minsurf_data, nx, ny) hx = 1 / (nx + 1) hy = 1 / (ny + 1) @@ -74,20 +75,7 @@ return c end -@inline function COPSBenchmark.minsurf_args(::ExaModelsBackend, nx::Int, ny::Int; T = Float64) - x_mesh = LinRange(0, 1, nx + 2) # coordinates of the mesh points x - - v0 = zeros(nx + 2, ny + 2) # Surface matrix initialization - for i = 1:(nx + 2), j = 1:(ny + 2) - v0[i, j] = 1 - (2 * x_mesh[i] - 1)^2 - end - - hx = 1 / (nx + 1) - hy = 1 / (ny + 1) - c6_i = Int(floor(0.25 / hx)):Int(ceil(0.75 / hx)) - c6_j = Int(floor(0.25 / hy)):Int(ceil(0.75 / hy)) - return (nx, ny, (; v0, c6_i, c6_j)) -end +@inline COPSBenchmark.minsurf_args(::ExaModelsBackend, nx::Int, ny::Int; T = Float64) = (nx, ny) @inline COPSBenchmark.minsurf_model(b::ExaModelsBackend, nx::Int, ny::Int; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/pinene.jl b/ext/COPSBenchmarkExaModels/pinene.jl index a4432fc..3c7e446 100644 --- a/ext/COPSBenchmarkExaModels/pinene.jl +++ b/ext/COPSBenchmarkExaModels/pinene.jl @@ -20,7 +20,8 @@ tf = T(36420) # tau[nm] zero_T = T(0) - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.pinene_data, nh) h = tf / nh # uniform interval length @@ -67,44 +68,7 @@ return core end -@inline function COPSBenchmark.pinene_args(::ExaModelsBackend, nh; T = Float64) - nc = 3 - ne = 5 - nm = 8 - bc = T[100, 0, 0, 0, 0] - tau = T[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] - tf = tau[nm] - h = tf / T(nh) - t = T[T(i-1)*h for i in 1:nh+1] - itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - - z = reshape(T[ - 88.35, 7.3, 2.3, 0.4, 1.75, - 76.4, 15.6, 4.5, 0.7, 2.8, - 65.1, 23.1, 5.3, 1.1, 5.8, - 50.4, 32.9, 6.0, 1.5, 9.3, - 37.5, 42.7, 6.0, 1.9, 12.0, - 25.9, 49.1, 5.9, 2.2, 17.0, - 14.0, 57.4, 5.1, 2.6, 21.0, - 4.5, 63.1, 3.8, 2.9, 25.7, - ], ne, nm)' - - v0 = zeros(T, nh, ne) - for i in 1:itau[1], s in 1:ne - v0[i, s] = bc[s] - end - for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne - v0[i, s] = z[j, s] - end - for i in itau[nm]+1:nh, s in 1:ne - v0[i, s] = z[nm, s] - end - - v_start = [v0[i, s] for i=1:nh, s=1:ne] - uc_start = [v0[i, s] for i=1:nh, j=1:nc, s=1:ne] - itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] - return (nh, (; v_start, uc_start, itr)) -end +@inline COPSBenchmark.pinene_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.pinene_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/polygon.jl b/ext/COPSBenchmarkExaModels/polygon.jl index 8360ee2..5dde638 100644 --- a/ext/COPSBenchmarkExaModels/polygon.jl +++ b/ext/COPSBenchmarkExaModels/polygon.jl @@ -13,7 +13,8 @@ @inline function COPSBenchmark.polygon_recipe( ::ExaModelsBackend; T = Float64, backend = nothing, ) - c, N, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + c, N = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.polygon_data, N) ExaModels.@add_var(c, r, N; lvar = 0.0, uvar = 1.0, start = 1.0) ExaModels.@add_var(c, θ, N; lvar = 0.0, uvar = T(π), start = d.θ0) @@ -34,12 +35,7 @@ return c end -@inline function COPSBenchmark.polygon_args(::ExaModelsBackend, n::Int; T = Float64) - N = div(n, 2) - θ0 = [i * π / (N - 1) - π / (N - 1) for i in 1:N] - pairs = [(i, j) for i in 1:N-1 for j in i+1:N] - return (N, (; θ0, pairs)) -end +@inline COPSBenchmark.polygon_args(::ExaModelsBackend, n::Int; T = Float64) = (div(n, 2),) @inline COPSBenchmark.polygon_model(b::ExaModelsBackend, n::Int; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/robot.jl b/ext/COPSBenchmarkExaModels/robot.jl index 8a3ee08..d1d7b56 100644 --- a/ext/COPSBenchmarkExaModels/robot.jl +++ b/ext/COPSBenchmarkExaModels/robot.jl @@ -26,7 +26,8 @@ four_pi_3 = T(4) * pi_T / T(3) zero_T = T(0) - core, nh, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.robot_data, nh) inv_nh = one(T) / nh @@ -44,14 +45,9 @@ # Final time ExaModels.@add_var(core, tf, 1; start = T(1), lvar = zero_T) - # The two moments of inertia were `@add_expr`s. That macro stores its body - # as a closure, and in a recipe the closure captures `Variable`s whose types - # carry the placeholder, which `instantiate` cannot reach into. A local - # helper is spliced at the point of use exactly as the subexpression was -- - # `@add_expr` adds no variables and no constraints, only a name -- so the - # expression trees, and the model, are unchanged. - I_the(i) = ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2*third - I_phi(i) = ((L-rho[i])^3+rho[i]^3)*third + # Expressions (matching JuMP @expressions) + ExaModels.@add_expr(core, I_the, ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2*third for i=1:nh+1) + ExaModels.@add_expr(core, I_phi, ((L-rho[i])^3+rho[i]^3)*third for i=1:nh+1) ExaModels.@add_obj(core, tf[1] for i=1:1) @@ -60,8 +56,8 @@ ExaModels.@add_con(core, c2, - phi[j] + phi[j-1] + half * tf[1]*inv_nh * (phi_dot[j] + phi_dot[j-1]) for j=2:nh+1) ExaModels.@add_con(core, c3, - the[j] + the[j-1] + half * tf[1]*inv_nh * (the_dot[j] + the_dot[j-1]) for j=2:nh+1) ExaModels.@add_con(core, c4, - rho_dot[j] + rho_dot[j-1] + half * tf[1]*inv_nh * (u_rho[j] + u_rho[j-1]) / L for j=2:nh+1) - ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + half * tf[1]*inv_nh * (u_the[j] / I_the(j) + u_the[j-1] / I_the(j-1)) for j=2:nh+1) - ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + half * tf[1]*inv_nh * (u_phi[j] / I_phi(j) + u_phi[j-1] / I_phi(j-1)) for j=2:nh+1) + ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + half * tf[1]*inv_nh * (u_the[j] / I_the[j] + u_the[j-1] / I_the[j-1]) for j=2:nh+1) + ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + half * tf[1]*inv_nh * (u_phi[j] / I_phi[j] + u_phi[j-1] / I_phi[j-1]) for j=2:nh+1) # Boundary conditions ExaModels.@add_con(core, c7, - rho[1] + rho0) @@ -80,15 +76,7 @@ return core end -@inline function COPSBenchmark.robot_args(::ExaModelsBackend, nh; T = Float64) - pi_T = T(pi) - inv_nh = T(1) / T(nh) - two_pi_3 = T(2) * pi_T / T(3) - four_pi_3 = T(4) * pi_T / T(3) - the_start = [two_pi_3*(T(k)*inv_nh)^2 for k=1:nh+1] - the_dot_start = [four_pi_3*(T(k)*inv_nh) for k=1:nh+1] - return (nh, (; the_start, the_dot_start)) -end +@inline COPSBenchmark.robot_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.robot_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/rocket.jl b/ext/COPSBenchmarkExaModels/rocket.jl index 3ea3b0a..bc39332 100644 --- a/ext/COPSBenchmarkExaModels/rocket.jl +++ b/ext/COPSBenchmarkExaModels/rocket.jl @@ -26,16 +26,17 @@ D_c = half * v_c * (m_0 / g_0) T_max = T_c * m_0 * g_0 - core, nh, s_v_start, s_m_start = ExaModels.ExaCore( - T; backend = backend, minimize = false, nargs = Val(3), + core, nh = ExaModels.ExaCore( + T; backend = backend, minimize = false, nargs = Val(1), ) + d = ExaModels.ArgNode1(COPSBenchmark.rocket_data, nh) inv_nh = one(T) / nh s_Th_start = T_max*half ExaModels.@add_var(core, h, 0:nh; start=h_0, lvar = h_0) - ExaModels.@add_var(core, v, 0:nh; start=s_v_start, lvar = v_0) - ExaModels.@add_var(core, m, 0:nh; start=s_m_start, lvar = m_f, uvar = m_0) + ExaModels.@add_var(core, v, 0:nh; start=d.v_start, lvar = v_0) + ExaModels.@add_var(core, m, 0:nh; start=d.m_start, lvar = m_f, uvar = m_0) ExaModels.@add_var(core, Th, 0:nh; start=s_Th_start, lvar = v_0, uvar = T_max) ExaModels.@add_var(core, step, 1; start=inv_nh, lvar = v_0) @@ -56,16 +57,7 @@ return core end -@inline function COPSBenchmark.rocket_args(::ExaModelsBackend, nh; T = Float64) - m_0 = T(1) - m_f = T(0.6) * m_0 - inv_nh = T(1) / T(nh) - # Precompute generator scalars so the broadcast closure stays isbits - # (Metal rejects Type{T} captures). - s_v_start = T[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] - s_m_start = T[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] - return (nh, s_v_start, s_m_start) -end +@inline COPSBenchmark.rocket_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.rocket_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/steering.jl b/ext/COPSBenchmarkExaModels/steering.jl index 118a43c..c83b642 100644 --- a/ext/COPSBenchmarkExaModels/steering.jl +++ b/ext/COPSBenchmarkExaModels/steering.jl @@ -17,12 +17,13 @@ xf = [T(NaN), T(5), T(45), T(0)] half = T(0.5) - core, nh, x_start = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.steering_data, nh) inv_nh = one(T) / nh ExaModels.@add_var(core, u, 1:nh+1; lvar = u_min, uvar = u_max, start=T(0)) # control - ExaModels.@add_var(core, x, 1:nh+1, 1:4; start=x_start) # state + ExaModels.@add_var(core, x, 1:nh+1, 1:4; start=d.x_start) # state ExaModels.@add_var(core, tf, 1; start=T(1)) # final time ExaModels.@add_obj(core, tf[1]) @@ -41,21 +42,7 @@ return core end -@inline function COPSBenchmark.steering_args(::ExaModelsBackend, nh; T = Float64) - inv_nh = T(1) / T(nh) - function gen_x0(k, i) - if i == 1 || i == 4 - return T(0) - elseif i == 2 - return T(5)*T(k)*inv_nh - elseif i == 3 - return T(45)*T(k)*inv_nh - else - return T(0) - end - end - return (nh, [gen_x0(i, j) for i=1:nh+1, j=1:4]) -end +@inline COPSBenchmark.steering_args(::ExaModelsBackend, nh; T = Float64) = (nh,) @inline COPSBenchmark.steering_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/ext/COPSBenchmarkExaModels/torsion.jl b/ext/COPSBenchmarkExaModels/torsion.jl index 13a671f..a8b2714 100644 --- a/ext/COPSBenchmarkExaModels/torsion.jl +++ b/ext/COPSBenchmarkExaModels/torsion.jl @@ -11,7 +11,8 @@ ) c_val = T(5.0) - core, nx, ny, d = ExaModels.ExaCore(T; backend = backend, nargs = Val(3)) + core, nx, ny = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + d = ExaModels.ArgNode2(COPSBenchmark.torsion_data, nx, ny) hx = one(T) / (nx + 1) hy = one(T) / (ny + 1) @@ -35,15 +36,7 @@ return core end -@inline function COPSBenchmark.torsion_args(::ExaModelsBackend, nx, ny; T = Float64) - hx = T(1.0 / (nx + 1.0)) - hy = T(1.0 / (ny + 1.0)) - D = T[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] - D_flat = [(k1, k2, D[k1, k2]) for k1 in 1:nx+2, k2 in 1:ny+2] - lcon = [-d for (_, _, d) in D_flat] - ucon = [d for (_, _, d) in D_flat] - return (nx, ny, (; D, D_flat, lcon, ucon)) -end +@inline COPSBenchmark.torsion_args(::ExaModelsBackend, nx, ny; T = Float64) = (nx, ny) @inline COPSBenchmark.torsion_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( diff --git a/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index dc82378..bcec00e 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -10,55 +10,425 @@ struct ExaModelsBackend <: AbstractModelerBackend end function bearing_model end function bearing_recipe end function bearing_args end +# Two sizes, so a two-argument deferred call; the recipe then carries only the +# two integers, which is what a builder ABI can accept. +function bearing_data(nx, ny) + b = 10 + e = 0.1 + hx = 2*pi / (nx+1) + hy = 2*b / (ny+1) + 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] + # `wq` and the eccentricity term are per-row coefficients; they ride with + # the indices they belong to rather than being recomputed symbolically. + lower = [(i, j, wq(i), wq(i+1)) for i in 1:nx+1, j in 1:ny+1] + upper = [(i, j, wq(i), wq(i-1)) for i in 2:nx+2, j in 2:ny+2] + lin = [(i, j, sin((i-1)*hx)) for i in 1:nx+2, j in 1:ny+2] + return (; v0, lower, upper, lin) +end function camshape_model end function camshape_recipe end function camshape_args end function catmix_model end function catmix_recipe end function catmix_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function catmix_data(nh; T = Float64) + ne = 2 + nc = 3 + # Concrete element type, not the keyword `T`: inside a compiled library `T` + # is a runtime value and `T[...]` leaves `_array_for(T::Type, ...)` + # unresolved for `--trim=safe`. The core converts on append, so the values + # land in the model's own type either way. + v_start = Float64[mod(j, ne) for i in 1:nh, j in 1:ne] + pp_start = Float64[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne] + return (; v_start, pp_start) +end function chain_model end function chain_recipe end function chain_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function chain_data(nh; T = Float64) + a = 1 + b = 3 + tmin = b > a ? 1 / 4 : 3 / 4 + u0 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] + x10 = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1] + x20 = [(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] + x30 = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1] + return (; u0, x10, x20, x30) +end function channel_model end function channel_recipe end function channel_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function channel_data(nh; T = Float64) + nc = 4 + nd = 4 + tf = T(1.0) + h = tf / nh + + rho = Float64[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] + t = Float64[(i-1)*h for i in 1:nh+1] + + # Initial value + v0 = zeros(T, nh, nd) + for i in 1:nh + v0[i, 1] = t[i]^2*(3 - 2*t[i]) + v0[i, 2] = 6*t[i]*(1 - t[i]) + v0[i, 3] = 6*(1 - 2*t[i]) + v0[i, 4] = -12 + end + uc0 = Float64[v0[i, s] for i in 1:nh, j in 1:nc, s in 1:nd] + + # fac[k+1] = k! + fac = Float64[factorial(k) for k in 0:nc+nd] + + con1_itr = Tuple{Int,Int,Int,T,T,T,T}[ + (i, j, s, h*rho[j]^1/fac[2], h*rho[j]^2/fac[3], h*rho[j]^3/fac[4], h*rho[j]^4/fac[5]) + for i in 1:nh, j in 1:nc, s in 1:nd + ] + + con2_itr = Tuple{Int,Int,Int,T,T,T,T,T,T,T,T}[ + (i, j, s, + (s <= 1 ? (rho[j]*h)^(1-s)/fac[1-s+1] : zero(T)), + (s <= 2 ? (rho[j]*h)^(2-s)/fac[2-s+1] : zero(T)), + (s <= 3 ? (rho[j]*h)^(3-s)/fac[3-s+1] : zero(T)), + (s <= 4 ? (rho[j]*h)^(4-s)/fac[4-s+1] : zero(T)), + h^(nd-s+1)*rho[j]^(1+nd-s)/fac[1+nd-s+1], + h^(nd-s+1)*rho[j]^(2+nd-s)/fac[2+nd-s+1], + h^(nd-s+1)*rho[j]^(3+nd-s)/fac[3+nd-s+1], + h^(nd-s+1)*rho[j]^(4+nd-s)/fac[4+nd-s+1]) + for i in 1:nh, j in 1:nc, s in 1:nd + ] + + cont_itr = Tuple{Int,Int,T,T,T,T,T,T,T,T}[ + (i, s, + (s <= 1 ? h^(1-s)/fac[1-s+1] : zero(T)), + (s <= 2 ? h^(2-s)/fac[2-s+1] : zero(T)), + (s <= 3 ? h^(3-s)/fac[3-s+1] : zero(T)), + (s <= 4 ? h^(4-s)/fac[4-s+1] : zero(T)), + h^(nd-s+1)/fac[1+nd-s+1], + h^(nd-s+1)/fac[2+nd-s+1], + h^(nd-s+1)/fac[3+nd-s+1], + h^(nd-s+1)/fac[4+nd-s+1]) + for i in 1:nh-1, s in 1:nd + ] + + coll_itr = Tuple{Int,Int,T,T,T,T}[ + (i, j, rho[j]^0/fac[1], rho[j]^1/fac[2], rho[j]^2/fac[3], rho[j]^3/fac[4]) + for i in 1:nh, j in 1:nc + ] + + # right BC coefficients, carried with the row index they apply to + bc3_cv = [h^(k-1)/fac[k] for k in 1:nd] + bc3_cw = [h^nd/fac[k+nd] for k in 1:nc] + bc4_cv = [h^(k-2)/fac[k-1] for k in 2:nd] + bc4_cw = [h^(nd-1)/fac[k+nd-1] for k in 1:nc] + bc3_itr = Tuple{Int,T,T,T,T,T,T,T,T}[(nh, bc3_cv[1], bc3_cv[2], bc3_cv[3], bc3_cv[4], + bc3_cw[1], bc3_cw[2], bc3_cw[3], bc3_cw[4])] + bc4_itr = Tuple{Int,T,T,T,T,T,T,T}[(nh, bc4_cv[1], bc4_cv[2], bc4_cv[3], + bc4_cw[1], bc4_cw[2], bc4_cw[3], bc4_cw[4])] + return (; v0, uc0, con1_itr, con2_itr, cont_itr, coll_itr, bc3_itr, bc4_itr) +end function clnlbeam_model end function dirichlet_model end function elec_model end function elec_recipe end function elec_args end +# The draw is data, and the seed is part of it. `Fix2` on a NAMED function keeps +# every type in the core named, and bakes the seed into a compiled library -- +# which is what a library that instantiates from one integer must do. +function elec_data(np, seed) + Random.seed!(seed) + theta = (2pi) .* rand(np) + phi = pi .* rand(np) + x0 = Float64[cos(theta[i])*sin(phi[i]) for i=1:np] + y0 = Float64[sin(theta[i])*sin(phi[i]) for i=1:np] + z0 = Float64[cos(phi[i]) for i=1:np] + itr = Tuple{Int,Int}[(i,j) for i in 1:np-1 for j in i+1:np] + return (; x0, y0, z0, itr) +end function gasoil_model end function gasoil_recipe end function gasoil_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function gasoil_data(nh; T = Float64) + nc = 4 + ne = 2 + nm = 21 + rho = Float64[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] + bc = [1, 1, 2, 0] + tau = Float64[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] + tf = tau[nm] + h = tf / T(nh) + t = Float64[T(i-1)*h for i in 1:nh+1] + + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(Float64[ + 1.0000, 0.0000, + 0.8105, 0.2000, + 0.6208, 0.2886, + 0.5258, 0.3010, + 0.4345, 0.3215, + 0.3903, 0.3123, + 0.3342, 0.2716, + 0.3034, 0.2551, + 0.2735, 0.2258, + 0.2405, 0.1959, + 0.2283, 0.1789, + 0.2071, 0.1457, + 0.1669, 0.1198, + 0.1530, 0.0909, + 0.1339, 0.0719, + 0.1265, 0.0561, + 0.1200, 0.0460, + 0.0990, 0.0280, + 0.0870, 0.0190, + 0.0770, 0.0140, + 0.0690, 0.0100, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = T(bc[s]) + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i = 1:nh, s = 1:ne] + uc_start = [v0[i, s] for i = 1:nh, j = 1:nc, s = 1:ne] + itr = Tuple{Int,Int,Int,T,T,T}[(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + return (; v_start, uc_start, itr) +end function glider_model end function glider_recipe end function glider_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function glider_data(nh; T = Float64) + x_0 = T(0); y_0 = T(1000); y_f = T(900); vx_0 = T(13.23) + inv_nh = T(1) / T(nh) + x_start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh] + y_start = [y_0 + (T(k)*inv_nh)*(y_f - y_0) for k in 0:nh] + return (; x_start, y_start) +end function henon_model end function lane_emden_model end function marine_model end function marine_recipe end function marine_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function marine_data(nh; T = Float64) + ne = 8 + nm = 21 + tau = collect(T, range(T(0), T(10), 21)) + tf = tau[nm] + h = tf / T(nh) + t = Float64[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(Float64[ + 20000.0, 17000.0, 10000.0, 15000.0, 12000.0, 9000.0, 7000.0, 3000.0, + 12445.0, 15411.0, 13040.0, 13338.0, 13484.0, 8426.0, 6615.0, 4022.0, + 7705.0, 13074.0, 14623.0, 11976.0, 12453.0, 9272.0, 6891.0, 5020.0, + 4664.0, 8579.0, 12434.0, 12603.0, 11738.0, 9710.0, 6821.0, 5722.0, + 2977.0, 7053.0, 11219.0, 11340.0, 13665.0, 8534.0, 6242.0, 5695.0, + 1769.0, 5054.0, 10065.0, 11232.0, 12112.0, 9600.0, 6647.0, 7034.0, + 943.0, 3907.0, 9473.0, 10334.0, 11115.0, 8826.0, 6842.0, 7348.0, + 581.0, 2624.0, 7421.0, 10297.0, 12427.0, 8747.0, 7199.0, 7684.0, + 355.0, 1744.0, 5369.0, 7748.0, 10057.0, 8698.0, 6542.0, 7410.0, + 223.0, 1272.0, 4713.0, 6869.0, 9564.0, 8766.0, 6810.0, 6961.0, + 137.0, 821.0, 3451.0, 6050.0, 8671.0, 8291.0, 6827.0, 7525.0, + 87.0, 577.0, 2649.0, 5454.0, 8430.0, 7411.0, 6423.0, 8388.0, + 49.0, 337.0, 2058.0, 4115.0, 7435.0, 7627.0, 6268.0, 7189.0, + 32.0, 228.0, 1440.0, 3790.0, 6474.0, 6658.0, 5859.0, 7467.0, + 17.0, 168.0, 1178.0, 3087.0, 6524.0, 5880.0, 5562.0, 7144.0, + 11.0, 99.0, 919.0, 2596.0, 5360.0, 5762.0, 4480.0, 7256.0, + 7.0, 65.0, 647.0, 1873.0, 4556.0, 5058.0, 4944.0, 7538.0, + 4.0, 44.0, 509.0, 1571.0, 4009.0, 4527.0, 4233.0, 6649.0, + 2.0, 27.0, 345.0, 1227.0, 3677.0, 4229.0, 3805.0, 6378.0, + 1.0, 20.0, 231.0, 934.0, 3197.0, 3695.0, 3159.0, 6454.0, + 1.0, 12.0, 198.0, 707.0, 2562.0, 3163.0, 3232.0, 5566.0, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = z[1, s] + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i=1:nh, s=1:ne] + itr = Tuple{Int,Int,Int,T,T,T}[(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + return (; v_start, itr) +end function methanol_model end function methanol_recipe end function methanol_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function methanol_data(nh; T = Float64) + ne = 3 + nm = 17 + tau = Float64[ + 0., 0.050, 0.065, 0.080, 0.123, 0.233, 0.273, 0.354, 0.397, 0.418, + 0.502, 0.553, 0.681, 0.750, 0.916, 0.937, 1.122, + ] + tf = tau[nm] + h = tf / T(nh) + t = Float64[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(Float64[ + 1.0000, 0.0000, 0.0000, + 0.7085, 0.1621, 0.0811, + 0.5971, 0.1855, 0.0965, + 0.5537, 0.1989, 0.1198, + 0.3684, 0.2845, 0.1535, + 0.1712, 0.3491, 0.2097, + 0.1198, 0.3098, 0.2628, + 0.0747, 0.3576, 0.2467, + 0.0529, 0.3347, 0.2884, + 0.0415, 0.3388, 0.2757, + 0.0261, 0.3557, 0.3167, + 0.0208, 0.3483, 0.2954, + 0.0085, 0.3836, 0.2950, + 0.0053, 0.3611, 0.2937, + 0.0019, 0.3609, 0.2831, + 0.0018, 0.3485, 0.2846, + 0.0006, 0.3698, 0.2899, + ], ne, nm)' + + con1_matrix = Tuple{Int,Int,Int,T,T,T}[(j, s, itau[j], tau[j], z[j,s], t[itau[j]]) for j in 1:nm, s in 1:ne] + return (; con1_matrix) +end function minsurf_model end function minsurf_recipe end function minsurf_args end +# Two sizes, so a two-argument deferred call; the recipe then carries only the +# two integers, which is what a builder ABI can accept. +function minsurf_data(nx, ny) + x_mesh = LinRange(0, 1, nx + 2) # coordinates of the mesh points x + + v0 = zeros(nx + 2, ny + 2) # Surface matrix initialization + for i = 1:(nx + 2), j = 1:(ny + 2) + v0[i, j] = 1 - (2 * x_mesh[i] - 1)^2 + end + + hx = 1 / (nx + 1) + hy = 1 / (ny + 1) + c6_i = Int(floor(0.25 / hx)):Int(ceil(0.75 / hx)) + c6_j = Int(floor(0.25 / hy)):Int(ceil(0.75 / hy)) + return (; v0, c6_i, c6_j) +end function pinene_model end function pinene_recipe end function pinene_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function pinene_data(nh; T = Float64) + nc = 3 + ne = 5 + nm = 8 + bc = Float64[100, 0, 0, 0, 0] + tau = Float64[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] + tf = tau[nm] + h = tf / T(nh) + t = Float64[T(i-1)*h for i in 1:nh+1] + itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] + + z = reshape(Float64[ + 88.35, 7.3, 2.3, 0.4, 1.75, + 76.4, 15.6, 4.5, 0.7, 2.8, + 65.1, 23.1, 5.3, 1.1, 5.8, + 50.4, 32.9, 6.0, 1.5, 9.3, + 37.5, 42.7, 6.0, 1.9, 12.0, + 25.9, 49.1, 5.9, 2.2, 17.0, + 14.0, 57.4, 5.1, 2.6, 21.0, + 4.5, 63.1, 3.8, 2.9, 25.7, + ], ne, nm)' + + v0 = zeros(T, nh, ne) + for i in 1:itau[1], s in 1:ne + v0[i, s] = bc[s] + end + for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne + v0[i, s] = z[j, s] + end + for i in itau[nm]+1:nh, s in 1:ne + v0[i, s] = z[nm, s] + end + + v_start = [v0[i, s] for i=1:nh, s=1:ne] + uc_start = [v0[i, s] for i=1:nh, j=1:nc, s=1:ne] + itr = Tuple{Int,Int,Int,T,T,T}[(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne] + return (; v_start, uc_start, itr) +end function polygon_model end function polygon_recipe end function polygon_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function polygon_data(N; T = Float64) + θ0 = [i * π / (N - 1) - π / (N - 1) for i in 1:N] + pairs = Tuple{Int,Int}[(i, j) for i in 1:N-1 for j in i+1:N] + return (; θ0, pairs) +end function robot_model end function robot_recipe end function robot_args end +# Derived tables, reached from the recipe through a deferred call so the +# recipe keeps a single placeholder -- which is all `P_new` carries. +function robot_data(nh; T = Float64) + pi_T = T(pi) + inv_nh = T(1) / T(nh) + two_pi_3 = T(2) * pi_T / T(3) + four_pi_3 = T(4) * pi_T / T(3) + the_start = [two_pi_3*(T(k)*inv_nh)^2 for k=1:nh+1] + the_dot_start = [four_pi_3*(T(k)*inv_nh) for k=1:nh+1] + return (; the_start, the_dot_start) +end function rocket_model end function rocket_recipe end function rocket_args end +# Derived tables live here, not in `_args`: a recipe that reaches them through +# a deferred call keeps a single placeholder, which is all `P_new` carries. +function rocket_data(nh; T = Float64) + m_0 = T(1) + m_f = T(0.6) * m_0 + inv_nh = T(1) / T(nh) + v_start = Float64[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] + m_start = Float64[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] + return (; v_start, m_start) +end function steering_model end function steering_recipe end function steering_args end +# Derived tables, reached through a deferred call so the recipe keeps a single +# placeholder. Concrete element type: inside a compiled library the keyword `T` +# is a runtime value and `T[...]` will not trim. +function steering_data(nh) + inv_nh = 1.0 / nh + gen_x0(k, i) = i == 2 ? 5.0*k*inv_nh : (i == 3 ? 45.0*k*inv_nh : 0.0) + x_start = Float64[gen_x0(i, j) for i = 1:nh+1, j = 1:4] + return (; x_start) +end function tetra_duct12_model end function tetra_duct15_model end function tetra_duct20_model end @@ -68,6 +438,17 @@ function tetra_hook_model end function torsion_model end function torsion_recipe end function torsion_args end +# Two sizes, so a two-argument deferred call; the recipe then carries only the +# two integers, which is what a builder ABI can accept. +function torsion_data(nx, ny) + hx = Float64(1.0 / (nx + 1.0)) + hy = Float64(1.0 / (ny + 1.0)) + D = Float64[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] + D_flat = [(k1, k2, D[k1, k2]) for k1 in 1:nx+2, k2 in 1:ny+2] + lcon = [-d for (_, _, d) in D_flat] + ucon = [d for (_, _, d) in D_flat] + return (; D, D_flat, lcon, ucon) +end function triangle_deer_model end function triangle_pacman_model end function triangle_turtle_model end From 9f9e52af61625fcb5bff77f7d24f9a4d4723b2bd Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Thu, 13 Aug 2026 07:05:40 -0400 Subject: [PATCH 4/7] test: a recipe testset per problem MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Checks what `*_model` cannot show now that it is defined as the composition: the core declares the arity its `*_args` supplies and its variable count is an expression rather than a number; and it is not consumed by being used — one core instantiates at two sizes and again at the first, each time equal to what `*_model` builds at that size on dimensions, starting point, objective and gradient. Co-Authored-By: Claude Opus 5 --- test/runtests.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index bbc6891..7a1177a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -234,3 +234,48 @@ end @test sort(abs.(eigvals(H_jump))) ≈ sort(abs.(eigvals(H_exa))) atol = 1e-6 end end + +# ── Recipes ─────────────────────────────────────────────────────────────────── +# A recipe is the model's structure with its sizes left open; `*_args` supplies +# the values that close it, and `*_model` is the two composed. Two properties +# `*_model` cannot show are checked here: the core really is open (it declares +# its arity and its variable count is an expression), and it is not consumed by +# being used — one core instantiates at two sizes and again at the first, each +# time matching what `*_model` builds at that size. +RECIPE_INSTANCES = [ + (name, params) for (name, params) in [ + ("bearing", (8, 8)), ("camshape", (60,)), ("catmix", (8,)), + ("chain", (60,)), ("channel", (12,)), ("elec", (15,)), + ("gasoil", (8,)), ("glider", (8,)), ("marine", (8,)), + ("methanol", (8,)), ("minsurf", (8, 8)), ("pinene", (8,)), + ("polygon", (30,)), ("robot", (12,)), ("rocket", (12,)), + ("steering", (12,)), ("torsion", (8, 8)), + ] +] + +@testset "Recipe: $name" for (name, params) in RECIPE_INSTANCES + B = COPSBenchmark.ExaModelsBackend() + recipe = getfield(COPSBenchmark, Symbol(name, "_recipe")) + argsf = getfield(COPSBenchmark, Symbol(name, "_args")) + modelf = getfield(COPSBenchmark, Symbol(name, "_model")) + + core = recipe(B) + @test core.nargs === Val(length(argsf(B, params...))) + @test core.nvar isa ExaModels.AbstractArgNode + + bigger = map(p -> p + 4, params) + first_meta = nothing + for ps in (params, bigger) + m = ExaModels.ExaModel(core, argsf(B, ps...)...) + r = modelf(B, ps...) + @test m.meta.nvar == r.meta.nvar + @test m.meta.ncon == r.meta.ncon + @test m.meta.x0 == r.meta.x0 + x = r.meta.x0 .+ 0.001 .* (1:r.meta.nvar) + @test NLPModels.obj(m, x) == NLPModels.obj(r, x) + @test NLPModels.grad(m, x) == NLPModels.grad(r, x) + ps == params && (first_meta = (m.meta.nvar, m.meta.ncon, copy(m.meta.x0))) + end + again = ExaModels.ExaModel(core, argsf(B, params...)...) + @test (again.meta.nvar, again.meta.ncon, again.meta.x0) == first_meta +end From b3ba5e3b9879eabe96af6d0a866007d4f72c50d8 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Thu, 13 Aug 2026 09:25:55 -0400 Subject: [PATCH 5/7] Thread T through the data functions; align style with LuksanVlcekBenchmark MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sungho's ruling on the Float32 finding: the `*_data` functions regain main's T-generic semantics rather than declaring Float64 canonical. T travels as a leading `::Val{T}` argument carried inside the deferred node — a bare `::Type{T}` stores as an abstract `DataType` field and the deferred call stays dynamic under `--trim=safe`; the `Val` singleton carries the type in its own type, and the where-clause specialization makes the `T[...]` comprehensions static. The `itau` integer index maps are computed in T again, exactly as on main: measured identical to main at Float32 across all eleven T-threaded models, structure and values, and identical at Float64 across all twenty harness cases. Data functions whose main-branch bodies were Float64 arithmetic (bearing, chain, elec, minsurf, polygon) keep size-only signatures — threading T through them would invent a dependency main never had. Style parity with LuksanVlcekBenchmark: `*_args` signatures carry only the sizes they consume (T flows through the recipe; elec's seed is baked there); the deferred-callable rule is stated in the extension entry file in the same wording; the test pin moves to madsuite-org's #308 branch with its repoint condition in the comment. Known regression, not ours: the #308 tip currently fails `--trim=safe` on cores carrying a deferred collect over a product iterator with a symbolic factor — every 2-D COPS model. Verified against an earlier ExaModels lineage with identical model-facing content: gasoil compiles and matches exactly at sizes the library never saw. CI on this PR stays red on the compile-adjacent legs until the branch fix lands. Co-Authored-By: Claude Opus 5 --- .../COPSBenchmarkExaModels.jl | 11 +++ ext/COPSBenchmarkExaModels/bearing.jl | 4 +- ext/COPSBenchmarkExaModels/catmix.jl | 6 +- ext/COPSBenchmarkExaModels/chain.jl | 4 +- ext/COPSBenchmarkExaModels/channel.jl | 6 +- ext/COPSBenchmarkExaModels/elec.jl | 6 +- ext/COPSBenchmarkExaModels/gasoil.jl | 6 +- ext/COPSBenchmarkExaModels/glider.jl | 6 +- ext/COPSBenchmarkExaModels/marine.jl | 6 +- ext/COPSBenchmarkExaModels/methanol.jl | 6 +- ext/COPSBenchmarkExaModels/minsurf.jl | 4 +- ext/COPSBenchmarkExaModels/pinene.jl | 6 +- ext/COPSBenchmarkExaModels/polygon.jl | 4 +- ext/COPSBenchmarkExaModels/robot.jl | 6 +- ext/COPSBenchmarkExaModels/rocket.jl | 6 +- ext/COPSBenchmarkExaModels/steering.jl | 6 +- ext/COPSBenchmarkExaModels/torsion.jl | 6 +- src/COPSBenchmark.jl | 80 +++++++++---------- test/Project.toml | 8 +- 19 files changed, 100 insertions(+), 87 deletions(-) diff --git a/ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl b/ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl index 40e6638..139610e 100644 --- a/ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl +++ b/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") diff --git a/ext/COPSBenchmarkExaModels/bearing.jl b/ext/COPSBenchmarkExaModels/bearing.jl index 9325e57..622b831 100644 --- a/ext/COPSBenchmarkExaModels/bearing.jl +++ b/ext/COPSBenchmarkExaModels/bearing.jl @@ -42,11 +42,11 @@ return core end -@inline COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny; T = Float64) = (nx, ny) +@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; T = T)...; + COPSBenchmark.bearing_args(b, nx, ny)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/catmix.jl b/ext/COPSBenchmarkExaModels/catmix.jl index ca2d352..74472ba 100644 --- a/ext/COPSBenchmarkExaModels/catmix.jl +++ b/ext/COPSBenchmarkExaModels/catmix.jl @@ -25,7 +25,7 @@ rho_index = [(i, rho[i]) for i in 1:nc] core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.catmix_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.catmix_data, (Val(T), nh)) h = one(T) / nh # Final time / nh @@ -84,11 +84,11 @@ return core end -@inline COPSBenchmark.catmix_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@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; T = T)...; + COPSBenchmark.catmix_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/chain.jl b/ext/COPSBenchmarkExaModels/chain.jl index 5893816..efa5720 100644 --- a/ext/COPSBenchmarkExaModels/chain.jl +++ b/ext/COPSBenchmarkExaModels/chain.jl @@ -83,11 +83,11 @@ return c end -@inline COPSBenchmark.chain_args(::ExaModelsBackend, n; T = Float64) = (max(2, div(n - 4, 4)),) +@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; T = T)...; + COPSBenchmark.chain_args(b, n)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/channel.jl b/ext/COPSBenchmarkExaModels/channel.jl index 85d2a28..612f64e 100644 --- a/ext/COPSBenchmarkExaModels/channel.jl +++ b/ext/COPSBenchmarkExaModels/channel.jl @@ -15,7 +15,7 @@ bc = T[0.0 1.0; 0.0 0.0] core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.channel_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.channel_data, (Val(T), nh)) ExaModels.@add_var(core, v, nh, 4; start = d.v0) ExaModels.@add_var(core, w, nh, 4; start = 0.0) @@ -72,11 +72,11 @@ return core end -@inline COPSBenchmark.channel_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.channel_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.channel_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.channel_recipe(b; T = T, backend = backend), - COPSBenchmark.channel_args(b, nh; T = T)...; + COPSBenchmark.channel_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/elec.jl b/ext/COPSBenchmarkExaModels/elec.jl index 5d26614..99c08b5 100644 --- a/ext/COPSBenchmarkExaModels/elec.jl +++ b/ext/COPSBenchmarkExaModels/elec.jl @@ -12,7 +12,7 @@ ::ExaModelsBackend; seed = 2713, T = Float64, backend = nothing, ) core, np = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(Base.Fix2(COPSBenchmark.elec_data, seed), np) + d = ExaModels.ArgCall(COPSBenchmark.elec_data, (np, seed)) ExaModels.@add_var(core, x, 1:np; start = d.x0) ExaModels.@add_var(core, y, 1:np; start = d.y0) @@ -27,11 +27,11 @@ return core end -@inline COPSBenchmark.elec_args(::ExaModelsBackend, np; seed = 2713, T = Float64) = (np,) +@inline COPSBenchmark.elec_args(::ExaModelsBackend, np) = (np,) @inline COPSBenchmark.elec_model(b::ExaModelsBackend, np; seed = 2713, T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.elec_recipe(b; seed = seed, T = T, backend = backend), - COPSBenchmark.elec_args(b, np; seed = seed, T = T)...; + COPSBenchmark.elec_args(b, np)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/gasoil.jl b/ext/COPSBenchmarkExaModels/gasoil.jl index a6929a5..b4a1a2f 100644 --- a/ext/COPSBenchmarkExaModels/gasoil.jl +++ b/ext/COPSBenchmarkExaModels/gasoil.jl @@ -23,7 +23,7 @@ z1 = T[1.0000, 0.0000] core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.gasoil_data, nh) + d = ExaModels.ArgNode2(COPSBenchmark.gasoil_data, Val(T), nh) h = tf / nh @@ -78,11 +78,11 @@ return core end -@inline COPSBenchmark.gasoil_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.gasoil_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.gasoil_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.gasoil_recipe(b; T = T, backend = backend), - COPSBenchmark.gasoil_args(b, nh; T = T)...; + COPSBenchmark.gasoil_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/glider.jl b/ext/COPSBenchmarkExaModels/glider.jl index 129e477..81f53f0 100644 --- a/ext/COPSBenchmarkExaModels/glider.jl +++ b/ext/COPSBenchmarkExaModels/glider.jl @@ -32,7 +32,7 @@ one_T = T(1) c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.glider_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.glider_data, (Val(T), nh)) inv_nh = one(T) / nh @@ -73,11 +73,11 @@ return c end -@inline COPSBenchmark.glider_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.glider_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.glider_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.glider_recipe(b; T = T, backend = backend), - COPSBenchmark.glider_args(b, nh; T = T)...; + COPSBenchmark.glider_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/marine.jl b/ext/COPSBenchmarkExaModels/marine.jl index 4f85811..38ece32 100644 --- a/ext/COPSBenchmarkExaModels/marine.jl +++ b/ext/COPSBenchmarkExaModels/marine.jl @@ -20,7 +20,7 @@ zero_T = T(0) core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.marine_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.marine_data, (Val(T), nh)) h = tf / nh # uniform interval length @@ -81,11 +81,11 @@ return core end -@inline COPSBenchmark.marine_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.marine_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.marine_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.marine_recipe(b; T = T, backend = backend), - COPSBenchmark.marine_args(b, nh; T = T)...; + COPSBenchmark.marine_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/methanol.jl b/ext/COPSBenchmarkExaModels/methanol.jl index d6a7503..dee4cfb 100644 --- a/ext/COPSBenchmarkExaModels/methanol.jl +++ b/ext/COPSBenchmarkExaModels/methanol.jl @@ -22,7 +22,7 @@ bc = T[1, 0, 0] c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.methanol_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.methanol_data, (Val(T), nh)) h = tf / nh # uniform interval length @@ -83,11 +83,11 @@ return c end -@inline COPSBenchmark.methanol_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.methanol_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.methanol_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.methanol_recipe(b; T = T, backend = backend), - COPSBenchmark.methanol_args(b, nh; T = T)...; + COPSBenchmark.methanol_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/minsurf.jl b/ext/COPSBenchmarkExaModels/minsurf.jl index 865a0d3..19d5b6a 100644 --- a/ext/COPSBenchmarkExaModels/minsurf.jl +++ b/ext/COPSBenchmarkExaModels/minsurf.jl @@ -75,11 +75,11 @@ return c end -@inline COPSBenchmark.minsurf_args(::ExaModelsBackend, nx::Int, ny::Int; T = Float64) = (nx, ny) +@inline COPSBenchmark.minsurf_args(::ExaModelsBackend, nx::Int, ny::Int) = (nx, ny) @inline COPSBenchmark.minsurf_model(b::ExaModelsBackend, nx::Int, ny::Int; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.minsurf_recipe(b; T = T, backend = backend), - COPSBenchmark.minsurf_args(b, nx, ny; T = T)...; + COPSBenchmark.minsurf_args(b, nx, ny)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/pinene.jl b/ext/COPSBenchmarkExaModels/pinene.jl index 3c7e446..e31e215 100644 --- a/ext/COPSBenchmarkExaModels/pinene.jl +++ b/ext/COPSBenchmarkExaModels/pinene.jl @@ -21,7 +21,7 @@ zero_T = T(0) core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.pinene_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.pinene_data, (Val(T), nh)) h = tf / nh # uniform interval length @@ -68,11 +68,11 @@ return core end -@inline COPSBenchmark.pinene_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.pinene_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.pinene_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.pinene_recipe(b; T = T, backend = backend), - COPSBenchmark.pinene_args(b, nh; T = T)...; + COPSBenchmark.pinene_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/polygon.jl b/ext/COPSBenchmarkExaModels/polygon.jl index 5dde638..db9a094 100644 --- a/ext/COPSBenchmarkExaModels/polygon.jl +++ b/ext/COPSBenchmarkExaModels/polygon.jl @@ -35,11 +35,11 @@ return c end -@inline COPSBenchmark.polygon_args(::ExaModelsBackend, n::Int; T = Float64) = (div(n, 2),) +@inline COPSBenchmark.polygon_args(::ExaModelsBackend, n::Int) = (div(n, 2),) @inline COPSBenchmark.polygon_model(b::ExaModelsBackend, n::Int; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.polygon_recipe(b; T = T, backend = backend), - COPSBenchmark.polygon_args(b, n; T = T)...; + COPSBenchmark.polygon_args(b, n)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/robot.jl b/ext/COPSBenchmarkExaModels/robot.jl index d1d7b56..5ecdf43 100644 --- a/ext/COPSBenchmarkExaModels/robot.jl +++ b/ext/COPSBenchmarkExaModels/robot.jl @@ -27,7 +27,7 @@ zero_T = T(0) core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.robot_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.robot_data, (Val(T), nh)) inv_nh = one(T) / nh @@ -76,11 +76,11 @@ return core end -@inline COPSBenchmark.robot_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.robot_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.robot_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.robot_recipe(b; T = T, backend = backend), - COPSBenchmark.robot_args(b, nh; T = T)...; + COPSBenchmark.robot_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/rocket.jl b/ext/COPSBenchmarkExaModels/rocket.jl index bc39332..0888a81 100644 --- a/ext/COPSBenchmarkExaModels/rocket.jl +++ b/ext/COPSBenchmarkExaModels/rocket.jl @@ -29,7 +29,7 @@ core, nh = ExaModels.ExaCore( T; backend = backend, minimize = false, nargs = Val(1), ) - d = ExaModels.ArgNode1(COPSBenchmark.rocket_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.rocket_data, (Val(T), nh)) inv_nh = one(T) / nh s_Th_start = T_max*half @@ -57,11 +57,11 @@ return core end -@inline COPSBenchmark.rocket_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.rocket_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.rocket_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.rocket_recipe(b; T = T, backend = backend), - COPSBenchmark.rocket_args(b, nh; T = T)...; + COPSBenchmark.rocket_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/steering.jl b/ext/COPSBenchmarkExaModels/steering.jl index c83b642..9ba4311 100644 --- a/ext/COPSBenchmarkExaModels/steering.jl +++ b/ext/COPSBenchmarkExaModels/steering.jl @@ -18,7 +18,7 @@ half = T(0.5) core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) - d = ExaModels.ArgNode1(COPSBenchmark.steering_data, nh) + d = ExaModels.ArgCall(COPSBenchmark.steering_data, (Val(T), nh)) inv_nh = one(T) / nh @@ -42,11 +42,11 @@ return core end -@inline COPSBenchmark.steering_args(::ExaModelsBackend, nh; T = Float64) = (nh,) +@inline COPSBenchmark.steering_args(::ExaModelsBackend, nh) = (nh,) @inline COPSBenchmark.steering_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.steering_recipe(b; T = T, backend = backend), - COPSBenchmark.steering_args(b, nh; T = T)...; + COPSBenchmark.steering_args(b, nh)...; kwargs..., ) diff --git a/ext/COPSBenchmarkExaModels/torsion.jl b/ext/COPSBenchmarkExaModels/torsion.jl index a8b2714..548a52f 100644 --- a/ext/COPSBenchmarkExaModels/torsion.jl +++ b/ext/COPSBenchmarkExaModels/torsion.jl @@ -12,7 +12,7 @@ c_val = T(5.0) core, nx, ny = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) - d = ExaModels.ArgNode2(COPSBenchmark.torsion_data, nx, ny) + d = ExaModels.ArgCall(COPSBenchmark.torsion_data, (Val(T), nx, ny)) hx = one(T) / (nx + 1) hy = one(T) / (ny + 1) @@ -36,11 +36,11 @@ return core end -@inline COPSBenchmark.torsion_args(::ExaModelsBackend, nx, ny; T = Float64) = (nx, ny) +@inline COPSBenchmark.torsion_args(::ExaModelsBackend, nx, ny) = (nx, ny) @inline COPSBenchmark.torsion_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) = ExaModels.ExaModel( COPSBenchmark.torsion_recipe(b; T = T, backend = backend), - COPSBenchmark.torsion_args(b, nx, ny; T = T)...; + COPSBenchmark.torsion_args(b, nx, ny)...; kwargs..., ) diff --git a/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index bcec00e..46fe218 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -34,15 +34,11 @@ function catmix_recipe end function catmix_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function catmix_data(nh; T = Float64) +function catmix_data(::Val{T}, nh) where {T} ne = 2 nc = 3 - # Concrete element type, not the keyword `T`: inside a compiled library `T` - # is a runtime value and `T[...]` leaves `_array_for(T::Type, ...)` - # unresolved for `--trim=safe`. The core converts on append, so the values - # land in the model's own type either way. - v_start = Float64[mod(j, ne) for i in 1:nh, j in 1:ne] - pp_start = Float64[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne] + v_start = T[mod(j, ne) for i in 1:nh, j in 1:ne] + pp_start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne] return (; v_start, pp_start) end function chain_model end @@ -50,7 +46,7 @@ function chain_recipe end function chain_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function chain_data(nh; T = Float64) +function chain_data(nh) a = 1 b = 3 tmin = b > a ? 1 / 4 : 3 / 4 @@ -66,14 +62,14 @@ function channel_recipe end function channel_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function channel_data(nh; T = Float64) +function channel_data(::Val{T}, nh) where {T} nc = 4 nd = 4 tf = T(1.0) h = tf / nh - rho = Float64[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] - t = Float64[(i-1)*h for i in 1:nh+1] + rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] + t = T[(i-1)*h for i in 1:nh+1] # Initial value v0 = zeros(T, nh, nd) @@ -83,10 +79,10 @@ function channel_data(nh; T = Float64) v0[i, 3] = 6*(1 - 2*t[i]) v0[i, 4] = -12 end - uc0 = Float64[v0[i, s] for i in 1:nh, j in 1:nc, s in 1:nd] + uc0 = T[v0[i, s] for i in 1:nh, j in 1:nc, s in 1:nd] # fac[k+1] = k! - fac = Float64[factorial(k) for k in 0:nc+nd] + fac = T[factorial(k) for k in 0:nc+nd] con1_itr = Tuple{Int,Int,Int,T,T,T,T}[ (i, j, s, h*rho[j]^1/fac[2], h*rho[j]^2/fac[3], h*rho[j]^3/fac[4], h*rho[j]^4/fac[5]) @@ -158,20 +154,20 @@ function gasoil_recipe end function gasoil_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function gasoil_data(nh; T = Float64) +function gasoil_data(::Val{T}, nh) where {T} nc = 4 ne = 2 nm = 21 - rho = Float64[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] + rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703] bc = [1, 1, 2, 0] - tau = Float64[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] + tau = T[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95] tf = tau[nm] h = tf / T(nh) - t = Float64[T(i-1)*h for i in 1:nh+1] + t = T[T(i-1)*h for i in 1:nh+1] itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - z = reshape(Float64[ + z = reshape(T[ 1.0000, 0.0000, 0.8105, 0.2000, 0.6208, 0.2886, @@ -216,7 +212,7 @@ function glider_recipe end function glider_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function glider_data(nh; T = Float64) +function glider_data(::Val{T}, nh) where {T} x_0 = T(0); y_0 = T(1000); y_f = T(900); vx_0 = T(13.23) inv_nh = T(1) / T(nh) x_start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh] @@ -230,16 +226,16 @@ function marine_recipe end function marine_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function marine_data(nh; T = Float64) +function marine_data(::Val{T}, nh) where {T} ne = 8 nm = 21 tau = collect(T, range(T(0), T(10), 21)) tf = tau[nm] h = tf / T(nh) - t = Float64[T(i-1)*h for i in 1:nh+1] + t = T[T(i-1)*h for i in 1:nh+1] itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - z = reshape(Float64[ + z = reshape(T[ 20000.0, 17000.0, 10000.0, 15000.0, 12000.0, 9000.0, 7000.0, 3000.0, 12445.0, 15411.0, 13040.0, 13338.0, 13484.0, 8426.0, 6615.0, 4022.0, 7705.0, 13074.0, 14623.0, 11976.0, 12453.0, 9272.0, 6891.0, 5020.0, @@ -283,19 +279,19 @@ function methanol_recipe end function methanol_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function methanol_data(nh; T = Float64) +function methanol_data(::Val{T}, nh) where {T} ne = 3 nm = 17 - tau = Float64[ + tau = T[ 0., 0.050, 0.065, 0.080, 0.123, 0.233, 0.273, 0.354, 0.397, 0.418, 0.502, 0.553, 0.681, 0.750, 0.916, 0.937, 1.122, ] tf = tau[nm] h = tf / T(nh) - t = Float64[T(i-1)*h for i in 1:nh+1] + t = T[T(i-1)*h for i in 1:nh+1] itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - z = reshape(Float64[ + z = reshape(T[ 1.0000, 0.0000, 0.0000, 0.7085, 0.1621, 0.0811, 0.5971, 0.1855, 0.0965, @@ -342,18 +338,18 @@ function pinene_recipe end function pinene_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function pinene_data(nh; T = Float64) +function pinene_data(::Val{T}, nh) where {T} nc = 3 ne = 5 nm = 8 - bc = Float64[100, 0, 0, 0, 0] - tau = Float64[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] + bc = T[100, 0, 0, 0, 0] + tau = T[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420] tf = tau[nm] h = tf / T(nh) - t = Float64[T(i-1)*h for i in 1:nh+1] + t = T[T(i-1)*h for i in 1:nh+1] itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm] - z = reshape(Float64[ + z = reshape(T[ 88.35, 7.3, 2.3, 0.4, 1.75, 76.4, 15.6, 4.5, 0.7, 2.8, 65.1, 23.1, 5.3, 1.1, 5.8, @@ -385,7 +381,7 @@ function polygon_recipe end function polygon_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function polygon_data(N; T = Float64) +function polygon_data(N) θ0 = [i * π / (N - 1) - π / (N - 1) for i in 1:N] pairs = Tuple{Int,Int}[(i, j) for i in 1:N-1 for j in i+1:N] return (; θ0, pairs) @@ -395,7 +391,7 @@ function robot_recipe end function robot_args end # Derived tables, reached from the recipe through a deferred call so the # recipe keeps a single placeholder -- which is all `P_new` carries. -function robot_data(nh; T = Float64) +function robot_data(::Val{T}, nh) where {T} pi_T = T(pi) inv_nh = T(1) / T(nh) two_pi_3 = T(2) * pi_T / T(3) @@ -409,12 +405,12 @@ function rocket_recipe end function rocket_args end # Derived tables live here, not in `_args`: a recipe that reaches them through # a deferred call keeps a single placeholder, which is all `P_new` carries. -function rocket_data(nh; T = Float64) +function rocket_data(::Val{T}, nh) where {T} m_0 = T(1) m_f = T(0.6) * m_0 inv_nh = T(1) / T(nh) - v_start = Float64[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] - m_start = Float64[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] + v_start = T[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh] + m_start = T[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh] return (; v_start, m_start) end function steering_model end @@ -423,10 +419,10 @@ function steering_args end # Derived tables, reached through a deferred call so the recipe keeps a single # placeholder. Concrete element type: inside a compiled library the keyword `T` # is a runtime value and `T[...]` will not trim. -function steering_data(nh) - inv_nh = 1.0 / nh - gen_x0(k, i) = i == 2 ? 5.0*k*inv_nh : (i == 3 ? 45.0*k*inv_nh : 0.0) - x_start = Float64[gen_x0(i, j) for i = 1:nh+1, j = 1:4] +function steering_data(::Val{T}, nh) where {T} + inv_nh = T(1) / T(nh) + gen_x0(k, i) = i == 2 ? T(5)*T(k)*inv_nh : (i == 3 ? T(45)*T(k)*inv_nh : T(0)) + x_start = T[gen_x0(i, j) for i = 1:nh+1, j = 1:4] return (; x_start) end function tetra_duct12_model end @@ -440,10 +436,10 @@ function torsion_recipe end function torsion_args end # Two sizes, so a two-argument deferred call; the recipe then carries only the # two integers, which is what a builder ABI can accept. -function torsion_data(nx, ny) +function torsion_data(::Val{T}, nx, ny) where {T} hx = Float64(1.0 / (nx + 1.0)) hy = Float64(1.0 / (ny + 1.0)) - D = Float64[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] + D = T[min(min(i, nx-i+1)*hx, min(j, ny-j+1)*hy) for i in 0:nx+1, j in 0:ny+1] D_flat = [(k1, k2, D[k1, k2]) for k1 in 1:nx+2, k2 in 1:ny+2] lcon = [-d for (_, _, d) in D_flat] ucon = [d for (_, _, d) in D_flat] diff --git a/test/Project.toml b/test/Project.toml index d1276b7..19cf9fe 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -11,4 +11,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [sources] -ExaModels = {rev = "main", url = "https://github.com/exanauts/ExaModels.jl"} +# The recipe API (`nargs`, deferred argument nodes) is on ExaModels main and in +# no tagged release yet; without this pin CI resolves the registry's v0.11.x +# and every ExaModels-backend test fails on the `nargs` keyword. This pin +# tracks the #308 branch because robot/glider recipes need `@add_expr`-as-node, +# which is there and not yet on main — REPOINT to main when #308 merges, and +# drop the pin once a release carries the API. +ExaModels = {rev = "abi/multi-model-producer", url = "https://github.com/madsuite-org/ExaModels.jl"} From 97520712d9316bb58c21f57d6be02ad4cf2f15fc Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Fri, 14 Aug 2026 10:40:24 -0400 Subject: [PATCH 6/7] Add the ExaModelsCompiler extension; pin ExaModels to a SHA compile_all is a method on ExaModelsCompiler's generic function, implemented here in an extension triggered by ExaModelsCompiler as a weakdep -- so a benchmark package never acquires a compiler toolchain by being loaded. The seventeen problems are derived from the package's own surface via the *_recipe names rather than hand-listed, so a model added to src/ is compiled without editing the extension. Instantiation spec belongs to the model: the three grid problems take (nx, ny) and the rest a single size. only/exclude go through ExaModelsCompiler.select, which validates names against what the package actually provides. The test pin names a SHA rather than a branch: a branch rev is not reproducible under julia-actions/cache, since the job restores ~/.julia from an earlier run and can reuse that clone. Verified against ExaModels main a9fbf9cc: all 17 models compile into one library (16.07 MB, 0 verifier errors); 16 of 17 read back exact at sizes never compiled for. steering fails with 'steering_block returned nonzero status 2', an upstream metadata fault tracked on the ExaModels side, not a model defect -- the model itself builds and solves. --- Project.toml | 2 ++ ext/COPSBenchmarkExaModelsCompiler.jl | 51 +++++++++++++++++++++++++++ test/Project.toml | 23 ++++++++---- 3 files changed, 69 insertions(+), 7 deletions(-) create mode 100644 ext/COPSBenchmarkExaModelsCompiler.jl diff --git a/Project.toml b/Project.toml index 6eced7f..ea659b7 100644 --- a/Project.toml +++ b/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] diff --git a/ext/COPSBenchmarkExaModelsCompiler.jl b/ext/COPSBenchmarkExaModelsCompiler.jl new file mode 100644 index 0000000..e1519df --- /dev/null +++ b/ext/COPSBenchmarkExaModelsCompiler.jl @@ -0,0 +1,51 @@ +# Compiling every COPS model into one shared library. +# +# The instantiation spec belongs to the model, not to the call: three of these +# problems are posed on a grid and take two sizes where the other fourteen take +# one, and the collocation problems carry their element type into the deferred +# data functions. So the caller says how big and in what precision, and this +# extension knows which model wants that as `(n,)` and which as `(nx, ny)`. +# +# It lives in an extension so the package never depends on the compiler: a +# benchmark set is a modelling package, and a compiler toolchain is not +# something you should acquire by loading one. +module COPSBenchmarkExaModelsCompiler + +import COPSBenchmark as CB +import ExaModelsCompiler +using ExaModelsCompiler: compile_library, select + +# Posed on a 2-D grid; everything else takes a single discretization size. +const GRID_MODELS = (:bearing, :minsurf, :torsion) + +_problems() = sort!([ + Symbol(chopsuffix(string(n), "_recipe")) for n in Base.names(CB; all = true) + if endswith(string(n), "_recipe") && !startswith(string(n), "#") +]) + +function ExaModelsCompiler.compile_all( + ::Val{CB}; + path = "@cops", + sizes = 100, + grid_sizes = (10, 10), + T = Float64, + only = nothing, + exclude = (), + kwargs..., +) + b = CB.ExaModelsBackend() + # Assemble every model, then let `select` filter: it is the shared contract + # and refuses a name this package does not provide, rather than quietly + # compiling a library missing the model the caller asked for. Building a + # recipe is cheap — the expense is the compile, which only the survivors + # reach. + models = map(_problems()) do name + recipe = getfield(CB, Symbol(name, :_recipe)) + args = getfield(CB, Symbol(name, :_args)) + spec = name in GRID_MODELS ? args(b, grid_sizes...) : args(b, sizes) + name => (recipe(b; T = T), spec...) + end + return compile_library(path, select(models; only, exclude)...; kwargs...) +end + +end # module diff --git a/test/Project.toml b/test/Project.toml index 19cf9fe..9ea68e5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -11,10 +11,19 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [sources] -# The recipe API (`nargs`, deferred argument nodes) is on ExaModels main and in -# no tagged release yet; without this pin CI resolves the registry's v0.11.x -# and every ExaModels-backend test fails on the `nargs` keyword. This pin -# tracks the #308 branch because robot/glider recipes need `@add_expr`-as-node, -# which is there and not yet on main — REPOINT to main when #308 merges, and -# drop the pin once a release carries the API. -ExaModels = {rev = "abi/multi-model-producer", url = "https://github.com/madsuite-org/ExaModels.jl"} +# The recipe API (`nargs`, deferred argument nodes, compile_all) is on ExaModels +# main and in no tagged release yet; without this pin CI resolves the registry's +# v0.11.x and every ExaModels-backend test fails on the `nargs` keyword. +# +# Pinned to a SHA rather than `rev = "main"` deliberately. A branch name is not +# reproducible under julia-actions/cache: the job restores ~/.julia from an +# earlier run and reuses that clone without re-fetching, so `main` resolves to +# whatever commit the cache happens to hold. That is how run 31789266458 built +# against a main from before `nargs` existed -- while reporting the pin as +# honoured, since the source URL and rev in the resolved status were both right. +# A SHA either is in the cached clone, in which case it is the code we asked +# for, or it is not and Pkg must fetch it. +# +# BUMP when a needed change lands on ExaModels main; DROP for a version bound at +# the first release carrying the API. +ExaModels = {rev = "a9fbf9ccb5afcfb973bc09ed10d7e6b5f24e07a3", url = "https://github.com/madsuite-org/ExaModels.jl"} From afde0114bf58f3ae7551e60ebb2ff2ba28cba1f5 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Fri, 14 Aug 2026 11:07:56 -0400 Subject: [PATCH 7/7] test: cover compile_all compile_all had no test at all: the provider's whole AOT surface -- which problems it offers, what each is closed with, and the selection contract -- was resting on manual runs. Most of it is testable without invoking a compiler, because select validates names before compile_library is reached, so the error paths exercise compile_all itself rather than a stand-in and cost nothing: - the derived problem list covers the package. The extension builds its models from the *_recipe names rather than a hand-written list precisely so it cannot drift as models are added; this is what makes that a fact rather than an intention. - every (recipe, args...) pair compile_all would hand the compiler closes into a model. This is what breaks when a recipe and its *_args disagree. - an unknown name is refused, rather than silently yielding a library missing the model the caller asked for; so is a selection that excludes everything, and so is a package that implements no method. Then one real compile of a single small model, because a list that assembles is not evidence that anything in it compiles. Bounded deliberately: if it proves too slow for the runner it belongs in a separate job, and CI is what tells us. --- test/Project.toml | 2 ++ test/runtests.jl | 49 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index 9ea68e5..355512a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2" +ExaModelsCompiler = "3d1e9a26-5b74-4f0c-9a2b-7c8f4e11d3a7" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -27,3 +28,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" # BUMP when a needed change lands on ExaModels main; DROP for a version bound at # the first release carrying the API. ExaModels = {rev = "a9fbf9ccb5afcfb973bc09ed10d7e6b5f24e07a3", url = "https://github.com/madsuite-org/ExaModels.jl"} +ExaModelsCompiler = {rev = "a9fbf9ccb5afcfb973bc09ed10d7e6b5f24e07a3", url = "https://github.com/madsuite-org/ExaModels.jl", subdir = "ExaModelsCompiler"} diff --git a/test/runtests.jl b/test/runtests.jl index 7a1177a..d1575ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using JuMP using Ipopt using ExaModels +using ExaModelsCompiler using NLPModels using NLPModelsIpopt using NLPModelsJuMP @@ -279,3 +280,51 @@ RECIPE_INSTANCES = [ again = ExaModels.ExaModel(core, argsf(B, params...)...) @test (again.meta.nvar, again.meta.ncon, again.meta.x0) == first_meta end + +# `compile_all` is the provider's whole AOT surface: which problems it offers, +# the arguments each is closed with, and the selection contract. Most of that +# is testable without invoking a compiler, because `select` validates names +# before `compile_library` is ever reached — so the error paths below exercise +# `compile_all` itself rather than a stand-in. One real compile follows, +# because a list that assembles is not evidence that anything in it compiles. +@testset "compile_all" begin + Bc = COPSBenchmark.ExaModelsBackend() + GRID = (:bearing, :minsurf, :torsion) + + # The list `compile_all` derives from the `*_recipe` names must cover the + # package. Deriving it is what keeps the extension from drifting as models + # are added; this test is what makes that a fact rather than an intention. + recipes = sort([Symbol(chopsuffix(string(n), "_recipe")) + for n in names(COPSBenchmark; all = true) + if endswith(string(n), "_recipe") && !startswith(string(n), "#")]) + @test !isempty(recipes) + @test sort(Symbol.(first.(RECIPE_INSTANCES))) == recipes + + # Every pair `compile_all` would hand the compiler has to close into a + # model. This is what breaks when a recipe and its `*_args` disagree, and + # it costs no compilation to find out. + for (name, params) in RECIPE_INSTANCES + recipe = getfield(COPSBenchmark, Symbol(name, :_recipe)) + argsf = getfield(COPSBenchmark, Symbol(name, :_args)) + m = ExaModels.ExaModel(recipe(Bc; T = Float64), argsf(Bc, params...)...) + @test m.meta.nvar > 0 + end + + # Selection contract: an unknown name is refused rather than silently + # yielding a library missing the model the caller asked for. Refused + # before any compilation, which is what makes this cheap. + @test_throws ArgumentError ExaModelsCompiler.compile_all( + COPSBenchmark; only = [:no_such_problem]) + @test_throws ArgumentError ExaModelsCompiler.compile_all( + COPSBenchmark; exclude = recipes) + @test_throws ArgumentError ExaModelsCompiler.compile_all(Base) + + # One real compile, on a single small non-grid model, exercising the whole + # path: recipe -> library -> load. + mktempdir() do dir + r = ExaModelsCompiler.compile_all(COPSBenchmark; + path = joinpath(dir, "copstest"), + sizes = 12, only = [:chain]) + @test isfile(r.libpath) + end +end