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/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 5ce70cc..622b831 100644 --- a/ext/COPSBenchmarkExaModels/bearing.jl +++ b/ext/COPSBenchmarkExaModels/bearing.jl @@ -4,41 +4,49 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.bearing_model(::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) +# Two sizes, so two placeholders. The grid spacings depend on them and keep +# their symbolic form; the starting field and the per-column weights are +# comprehensions over the grid and travel as data. +@inline function COPSBenchmark.bearing_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) b = 10 # grid is (0,2*pi)x(0,2*b) e = 0.1 # eccentricity + core, nx, ny = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + d = ExaModels.ArgNode2(COPSBenchmark.bearing_data, nx, ny) + hx = 2*pi / (nx+1) # grid spacing hy = 2*b / (ny+1) # grid spacing - area = 0.5*hx*hy # area of triangle - - wq(i) = (1.0 + e*cos((i-1)*hx))^3 - v0 = [max(sin((i-1)*hx), 0.0) for i in 1:nx+2, j in 1:ny+2] - - core = ExaModels.ExaCore(T; backend=backend, concrete = Val(true)) - ExaModels.@add_var(core, v, 1:nx+2, 1:ny+2; lvar = 0.0, start=v0) + ExaModels.@add_var(core, v, 1:nx+2, 1:ny+2; lvar = 0.0, start = d.v0) ExaModels.@add_obj( core, - 0.5*(hx*hy/6.0) * (wq(i) + 2*wq(i+1))*(((v[i+1,j]-v[i,j])/hx)^2 + ((v[i,j+1]-v[i,j])/hy)^2) for i in 1:nx+1, j in 1:ny+1 + 0.5*(hx*hy/6.0) * (w1 + 2*w2)*(((v[i+1,j]-v[i,j])/hx)^2 + ((v[i,j+1]-v[i,j])/hy)^2) for (i, j, w1, w2) in d.lower ) ExaModels.@add_obj( core, - 0.5*(hx*hy/6.0) * (2*wq(i) + 2*wq(i-1))*(((v[i-1,j]-v[i,j])/hx)^2 + ((v[i,j-1]-v[i,j])/hy)^2) for i in 2:nx+2, j in 2:ny+2 + 0.5*(hx*hy/6.0) * (2*w1 + 2*w2)*(((v[i-1,j]-v[i,j])/hx)^2 + ((v[i,j-1]-v[i,j])/hy)^2) for (i, j, w1, w2) in d.upper ) ExaModels.@add_obj( core, - -hx*hy*e*sin((i-1)*hx)*v[i, j] for i in 1:nx+2, j in 1:ny+2 + -hx*hy*e*s*v[i, j] for (i, j, s) in d.lin ) ExaModels.@add_con(core, c1, v[i, 1] for i in 1:nx+2) - ExaModels.@add_con(core, c2, v[i, ny+2] for i in 1:nx+2) + ExaModels.@add_con(core, c2, v[i, k] for i in 1:nx+2, k in (ny+2):(ny+2)) ExaModels.@add_con(core, c3, v[1, i] for i in 1:ny+2) - ExaModels.@add_con(core, c4, v[nx+2, i] for i in 1:ny+2) + ExaModels.@add_con(core, c4, v[k, i] for k in (nx+2):(nx+2), i in 1:ny+2) - return ExaModels.ExaModel(core; kwargs...) + return core end +@inline COPSBenchmark.bearing_args(::ExaModelsBackend, nx, ny) = (nx, ny) - +@inline COPSBenchmark.bearing_model(b::ExaModelsBackend, nx, ny; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.bearing_recipe(b; T = T, backend = backend), + COPSBenchmark.bearing_args(b, nx, ny)...; + kwargs..., + ) 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..74472ba 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,71 @@ one_T = T(1) rho_index = [(i, rho[i]) for i in 1:nc] - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - ExaModels.@add_var(c, u, nh, nc; lvar = zeros(T, nh, nc), uvar = ones(T, nh, nc), start = zeros(T, nh, nc)) - ExaModels.@add_var(c, v, nh, ne; start = T[mod(j, ne) for i in 1:nh, j in 1:ne]) - ExaModels.@add_var(c, w, nh, nc, ne; start = zeros(T, nh, nc, ne)) - ExaModels.@add_var(c, pp, nh, nc, ne; start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne]) - ExaModels.@add_var(c, Dpp, nh, nc, ne; start = zeros(T, nh, nc, ne)) - ExaModels.@add_var(c, ppf, ne; start = T[mod(i,ne) for i in 1:ne]) + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.catmix_data, (Val(T), nh)) + + h = one(T) / nh # Final time / nh + + ExaModels.@add_var(core, u, nh, nc; lvar = zero(T), uvar = one(T), start = zero(T)) + ExaModels.@add_var(core, v, nh, ne; start = d.v_start) + ExaModels.@add_var(core, w, nh, nc, ne; start = zero(T)) + ExaModels.@add_var(core, pp, nh, nc, ne; start = d.pp_start) + ExaModels.@add_var(core, Dpp, nh, nc, ne; start = zero(T)) + ExaModels.@add_var(core, ppf, ne; start = T[mod(i,ne) for i in 1:ne]) - ExaModels.@add_obj(c, neg_one + ppf[1] + ppf[2]) - ExaModels.@add_obj(c, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc) + ExaModels.@add_obj(core, neg_one + ppf[1] + ppf[2]) + ExaModels.@add_obj(core, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc) ExaModels.@add_con( - c, + core, c1, pp[i, k, s] - v[i, s] - h*sum(w[i, j, s]*(rho^j/T(factorial(j))) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne ) ExaModels.@add_con( - c, + core, c2, Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/T(factorial(j-1))) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne ) ExaModels.@add_con( - c, + core, c3, - ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / T(factorial(j)) for j in 1:nc) for s in 1:ne + ppf[s] - v[k, s] - h * sum(w[k, j, s] / T(factorial(j)) for j in 1:nc) for k in nh:nh, s in 1:ne ) ExaModels.@add_con( - c, + core, c4, v[i, s] + sum(w[i, j, s] * h / T(factorial(j)) for j in 1:nc) - v[i+1, s] for i in 1:nh-1, s in 1:ne ) - - ExaModels.@add_con( - c, + core, c5, Dpp[i,j,1] - u[i,j] * (ten_T*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc ) ExaModels.@add_con( - c, + core, c6, Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - ten_T*pp[i,j,2]) + (one_T - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc ) - ExaModels.@add_con( - c, + core, c7, v[1, s] - bc for (s, bc) in [(i, bc[i]) for i in 1:ne] ) - return ExaModels.ExaModel(c; kwargs...) + return core end + +@inline COPSBenchmark.catmix_args(::ExaModelsBackend, nh) = (nh,) + +@inline COPSBenchmark.catmix_model(b::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.catmix_recipe(b; T = T, backend = backend), + COPSBenchmark.catmix_args(b, nh)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/chain.jl b/ext/COPSBenchmarkExaModels/chain.jl index c7891c8..efa5720 100644 --- a/ext/COPSBenchmarkExaModels/chain.jl +++ b/ext/COPSBenchmarkExaModels/chain.jl @@ -8,24 +8,29 @@ # see "Benchmarking Optimization Software with COPS" # Argonne National Labs Technical Report ANL/MCS-246 (2004) -@inline function COPSBenchmark.chain_model(::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...) - nh = max(2, div(n - 4, 4)) - +# The public size `n` is not the discretization: `nh = max(2, div(n-4, 4))` is. +# `chain_args` does that reduction and builds the four starting curves, which +# are comprehensions over nh; the recipe is written against nh directly. +@inline function COPSBenchmark.chain_recipe( + ::ExaModelsBackend; T = Float64, backend = nothing, +) L = 4 a = 1 b = 3 - tmin = b > a ? 1 / 4 : 3 / 4 tf = 1.0 + + c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode1(COPSBenchmark.chain_data, nh) + h = tf / nh - c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - ExaModels.@add_var(c, u, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1]) - ExaModels.@add_var(c, x1, nh + 1; start = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1]) - ExaModels.@add_var(c, x2, nh + 1; start = [(4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a) * - (4 * abs(b - a) * (k / nh - tmin)) for k in 1:nh+1]) - ExaModels.@add_var(c, x3, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1]) + ExaModels.@add_var(c, u, nh + 1; start = d.u0) + ExaModels.@add_var(c, x1, nh + 1; start = d.x10) + ExaModels.@add_var(c, x2, nh + 1; start = d.x20) + ExaModels.@add_var(c, x3, nh + 1; start = d.x30) - ExaModels.@add_obj(c, x2[nh + 1]) + # Indexes the last point, so the index arrives as data. + ExaModels.@add_obj(c, x2[k] for k in (nh+1):(nh+1)) ExaModels.@add_con( c, @@ -42,7 +47,7 @@ ExaModels.@add_con( c, c3, - x1[nh + 1] - b + x1[k] - b for k in (nh+1):(nh+1) ) ExaModels.@add_con( @@ -60,7 +65,7 @@ ExaModels.@add_con( c, c6, - x3[nh+1] - L + x3[k] - L for k in (nh+1):(nh+1) ) ExaModels.@add_con( @@ -75,6 +80,14 @@ x3[j + 1] - x3[j] - 1 / 2 * h * (sqrt(1 + u[j]^2) + sqrt(1 + u[j + 1]^2)) for j in 1:nh ) - return ExaModels.ExaModel(c; kwargs...) + return c end +@inline COPSBenchmark.chain_args(::ExaModelsBackend, n) = (max(2, div(n - 4, 4)),) + +@inline COPSBenchmark.chain_model(b::ExaModelsBackend, n; T = Float64, backend = nothing, kwargs...) = + ExaModels.ExaModel( + COPSBenchmark.chain_recipe(b; T = T, backend = backend), + COPSBenchmark.chain_args(b, n)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/channel.jl b/ext/COPSBenchmarkExaModels/channel.jl index cb0d639..612f64e 100644 --- a/ext/COPSBenchmarkExaModels/channel.jl +++ b/ext/COPSBenchmarkExaModels/channel.jl @@ -5,88 +5,22 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.channel_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...) - nc = 4 - nd = 4 - R = T(10.0) - tf = T(1.0) - h = tf / nh - +# 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] - 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] - - # 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], - 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 - ] - - # 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], - h^(nd-s+1)/fac[4+nd-s+1]) - 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]) - for i in 1:nh, j in 1:nc - ] + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.channel_data, (Val(T), nh)) - # right BC coefficients - 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] - - 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) + 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) @@ -94,43 +28,55 @@ # 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 + 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 con2_itr + for (i, j, s, bv1, bv2, bv3, bv4, bw1, bw2, bw3, bw4) in d.con2_itr ) - # Boundary conditions + # 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, - 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] + 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, - 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] + 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 cont_itr + 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] - - 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 + 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 ExaModels.ExaModel(core; kwargs...) + return core end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/elec.jl b/ext/COPSBenchmarkExaModels/elec.jl index 31bfddc..99c08b5 100644 --- a/ext/COPSBenchmarkExaModels/elec.jl +++ b/ext/COPSBenchmarkExaModels/elec.jl @@ -4,25 +4,34 @@ # COPS 3.0 - November 2002 # COPS 3.1 - March 2004 -@inline function COPSBenchmark.elec_model(::ExaModelsBackend, np; seed = 2713, T = Float64, backend = nothing, kwargs...) - Random.seed!(seed) +# 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; seed = 2713, T = Float64, backend = nothing, +) + core, np = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.elec_data, (np, seed)) - # Set the starting point to a quasi-uniform distribution - # of electrons on a unit sphere - 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]) + 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 - 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) + 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 ExaModels.ExaModel(core; kwargs...) + return core end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/gasoil.jl b/ext/COPSBenchmarkExaModels/gasoil.jl index 8b49476..b4a1a2f 100644 --- a/ext/COPSBenchmarkExaModels/gasoil.jl +++ b/ext/COPSBenchmarkExaModels/gasoil.jl @@ -5,80 +5,39 @@ # 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgNode2(COPSBenchmark.gasoil_data, Val(T), nh) - # 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 +54,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 +75,14 @@ - 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/glider.jl b/ext/COPSBenchmarkExaModels/glider.jl index 9000e0c..81f53f0 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,18 +28,20 @@ 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.glider_data, (Val(T), nh)) + + 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) @@ -49,7 +53,7 @@ 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[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) @@ -60,11 +64,20 @@ # 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)) - ExaModels.ExaModel(c; kwargs...) + return c end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/marine.jl b/ext/COPSBenchmarkExaModels/marine.jl index 51da740..38ece32 100644 --- a/ext/COPSBenchmarkExaModels/marine.jl +++ b/ext/COPSBenchmarkExaModels/marine.jl @@ -5,79 +5,37 @@ # 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.marine_data, (Val(T), nh)) - # 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 +78,14 @@ for i=1:nh, j=1:nc, s=2:ne-1 ) - return ExaModels.ExaModel(core; kwargs...) + return core end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/methanol.jl b/ext/COPSBenchmarkExaModels/methanol.jl index 4a9835d..dee4cfb 100644 --- a/ext/COPSBenchmarkExaModels/methanol.jl +++ b/ext/COPSBenchmarkExaModels/methanol.jl @@ -5,86 +5,34 @@ # 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) - - # 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] - - # 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] - # 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)) + c, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.methanol_data, (Val(T), nh)) + + h = tf / nh # uniform interval length + 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 +80,14 @@ theta[4]*uc[i,j,1]) for i=1:nh, j=1:nc ) - return ExaModels.ExaModel(c; kwargs...) + return c end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/minsurf.jl b/ext/COPSBenchmarkExaModels/minsurf.jl index ba4b02a..19d5b6a 100644 --- a/ext/COPSBenchmarkExaModels/minsurf.jl +++ b/ext/COPSBenchmarkExaModels/minsurf.jl @@ -9,20 +9,21 @@ # 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + d = ExaModels.ArgNode2(COPSBenchmark.minsurf_data, nx, ny) 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 +36,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 +53,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 +67,19 @@ 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/pinene.jl b/ext/COPSBenchmarkExaModels/pinene.jl index 485666d..e31e215 100644 --- a/ext/COPSBenchmarkExaModels/pinene.jl +++ b/ext/COPSBenchmarkExaModels/pinene.jl @@ -5,63 +5,36 @@ # 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] + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.pinene_data, (Val(T), nh)) - # 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 = 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 +65,14 @@ 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/polygon.jl b/ext/COPSBenchmarkExaModels/polygon.jl index 340ac92..db9a094 100644 --- a/ext/COPSBenchmarkExaModels/polygon.jl +++ b/ext/COPSBenchmarkExaModels/polygon.jl @@ -6,27 +6,40 @@ # 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 = 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 = [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 ExaModels.ExaModel(c; kwargs...) + return c end + +@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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/robot.jl b/ext/COPSBenchmarkExaModels/robot.jl index dcbf135..5ecdf43 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,21 @@ 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.robot_data, (Val(T), nh)) + + 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) @@ -59,15 +63,24 @@ 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/rocket.jl b/ext/COPSBenchmarkExaModels/rocket.jl index 9b3fcaf..0888a81 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,28 @@ 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 = ExaModels.ExaCore( + T; backend = backend, minimize = false, nargs = Val(1), + ) + d = ExaModels.ArgCall(COPSBenchmark.rocket_data, (Val(T), nh)) - # 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, 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) - 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 +52,16 @@ 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/steering.jl b/ext/COPSBenchmarkExaModels/steering.jl index e28f8a2..9ba4311 100644 --- a/ext/COPSBenchmarkExaModels/steering.jl +++ b/ext/COPSBenchmarkExaModels/steering.jl @@ -4,31 +4,26 @@ # 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 = ExaModels.ExaCore(T; backend= backend, concrete = Val(true)) + + core, nh = ExaModels.ExaCore(T; backend = backend, nargs = Val(1)) + d = ExaModels.ArgCall(COPSBenchmark.steering_data, (Val(T), 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=[gen_x0(i, j) for i=1:nh+1, j=1:4]) # 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]) @@ -39,9 +34,19 @@ 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 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)...; + kwargs..., + ) diff --git a/ext/COPSBenchmarkExaModels/torsion.jl b/ext/COPSBenchmarkExaModels/torsion.jl index 56a7bf8..548a52f 100644 --- a/ext/COPSBenchmarkExaModels/torsion.jl +++ b/ext/COPSBenchmarkExaModels/torsion.jl @@ -3,35 +3,44 @@ # 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 = ExaModels.ExaCore(T; backend = backend, nargs = Val(2)) + d = ExaModels.ArgCall(COPSBenchmark.torsion_data, (Val(T), nx, ny)) + + hx = one(T) / (nx + 1) + hy = one(T) / (ny + 1) + area = T(0.5) * hx * hy - core = ExaModels.ExaCore(T; backend = backend, concrete = Val(true)) - ExaModels.@add_var(core, v, nx+2, ny+2; start = D) + 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 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)...; + kwargs..., + ) 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/src/COPSBenchmark.jl b/src/COPSBenchmark.jl index de7823d..46fe218 100644 --- a/src/COPSBenchmark.jl +++ b/src/COPSBenchmark.jl @@ -8,25 +8,423 @@ struct ExaModelsBackend <: AbstractModelerBackend end # COPS Instances 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(::Val{T}, nh) where {T} + 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 (; 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) + 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(::Val{T}, nh) where {T} + 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 = 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(::Val{T}, nh) where {T} + 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 = 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(::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] + 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(::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 = 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 = 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(::Val{T}, nh) where {T} + 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 = 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(::Val{T}, nh) where {T} + 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 = 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) + θ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(::Val{T}, nh) where {T} + 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(::Val{T}, nh) where {T} + m_0 = T(1) + m_f = T(0.6) * m_0 + inv_nh = T(1) / T(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 +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(::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 function tetra_duct15_model end function tetra_duct20_model end @@ -34,6 +432,19 @@ 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 +# 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(::Val{T}, nx, ny) where {T} + hx = Float64(1.0 / (nx + 1.0)) + hy = Float64(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 (; D, D_flat, lcon, ucon) +end function triangle_deer_model end function triangle_pacman_model end function triangle_turtle_model end diff --git a/test/Project.toml b/test/Project.toml index d1276b7..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" @@ -11,4 +12,20 @@ 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, 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"} +ExaModelsCompiler = {rev = "a9fbf9ccb5afcfb973bc09ed10d7e6b5f24e07a3", url = "https://github.com/madsuite-org/ExaModels.jl", subdir = "ExaModelsCompiler"} diff --git a/test/runtests.jl b/test/runtests.jl index bbc6891..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 @@ -234,3 +235,96 @@ 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 + +# `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