diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ad37434..c1298b0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,7 +29,6 @@ jobs: strategy: matrix: python-version: ["3.11", "3.12"] - julia-version: ["1.11", "1.12"] steps: - uses: actions/checkout@v4 @@ -39,21 +38,14 @@ jobs: - uses: julia-actions/setup-julia@v2 with: - version: ${{ matrix.julia-version }} + version: "1.12" - uses: julia-actions/cache@v2 - name: Install juliacall run: pip install juliacall - - name: Preinstall Julia packages - run: | - julia -e ' - using Pkg - Pkg.add(["MathOptInterface", "HiGHS"]) - using MathOptInterface - using HiGHS - ' - - name: Run solve tests run: python tests/test_solve.py + env: + JUMPY_BACKEND: juliacall diff --git a/.github/workflows/juliac.yml b/.github/workflows/juliac.yml new file mode 100644 index 0000000..52ead0f --- /dev/null +++ b/.github/workflows/juliac.yml @@ -0,0 +1,72 @@ +name: juliac backend + +on: + push: + branches: [main] + pull_request: + branches: [main] + workflow_dispatch: + +jobs: + test-entrypoints: + name: Julia entry-point tests (in-process) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - uses: julia-actions/setup-julia@v2 + with: + version: "1.12" + + - uses: julia-actions/cache@v2 + + - name: Instantiate + run: julia --project=julia -e 'using Pkg; Pkg.instantiate()' + + - name: Run entry-point tests + run: julia --project=julia julia/test/runtests.jl + + build-lib: + name: Build and test compiled library (${{ matrix.os }}) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + steps: + - uses: actions/checkout@v4 + + - uses: julia-actions/setup-julia@v2 + with: + version: "1.12" + + - uses: julia-actions/cache@v2 + + - name: Install JuliaC + run: julia --project=@juliac -e 'using Pkg; Pkg.add("JuliaC")' + + - name: Build shared library + working-directory: julia + run: > + julia --project=@juliac -m JuliaC + --output-lib build/libjumpy_highs + --compile-ccallable + --jl-option handle-signals=no + --bundle build + --verbose + . + + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Run end-to-end tests against the compiled library + run: python3 tests/test_solve.py + env: + JUMPY_BACKEND: juliac + + - name: Upload bundle + uses: actions/upload-artifact@v4 + with: + name: jumpy-backend-${{ matrix.os }} + path: julia/build diff --git a/README.md b/README.md index 7bdf5cf..fd3b2fd 100644 --- a/README.md +++ b/README.md @@ -49,8 +49,8 @@ m = jp.Model() x = m.variables(100, lower=0, name="x") # A constraint group: 99 constraints from one template -i = jp.Iterator(range(99)) -m.constraint_group([i], x[i] + x[i + 1] <= 10) +i = m.iterator(range(99)) +m.constraint_group(x[i] + x[i + 1] <= 10) m.objective = jp.minimize(x[0] + x[1]) m.optimize() @@ -65,35 +65,35 @@ Constraint groups are the core feature. Instead of building constraints one by o ### Basic group ```python -i = jp.Iterator(range(1000000)) -m.constraint_group([i], x[i] <= 10) +i = m.iterator(range(1000000)) +m.constraint_group(x[i] <= 10) # One template in Python → 1,000,000 constraints in Julia ``` ### Multi-dimensional ```python -i = jp.Iterator(range(100)) -j = jp.Iterator(range(100)) -m.constraint_group([i, j], x[100 * i + j] >= 0) +i = m.iterator(range(100)) +j = m.iterator(range(100)) +m.constraint_group(x[100 * i + j] >= 0) # 10,000 constraints from one template ``` ### With data ```python -costs = jp.Parameter([...], name="costs") -demand = jp.Parameter([...], name="demand") +costs = m.parameter([...], name="costs") +demand = m.parameter([...], name="demand") -i = jp.Iterator(range(n)) -m.constraint_group([i], costs[i] * x[i] >= demand[i]) +i = m.iterator(range(n)) +m.constraint_group(costs[i] * x[i] >= demand[i]) ``` ### Nonlinear ```python -i = jp.Iterator(range(n)) -m.constraint_group([i], jp.sin(x[i]) + jp.exp(x[i]) <= 1.0) +i = m.iterator(range(n)) +m.constraint_group(jp.sin(x[i]) + jp.exp(x[i]) <= 1.0) ``` ### Individual constraints @@ -111,12 +111,14 @@ m.constraint(x[0] + x[1] == 5) | Method | Description | |---|---| | `m = jp.Model()` | Create a new model | -| `m.variables(n, lower=, upper=, name=)` | Add `n` variables, returns a `VariableVector` | -| `m.variable(lower=, upper=, name=)` | Add a single variable | -| `m.constraint_group(iterators, template)` | Add a constraint group | +| `m.variables(n, lower=, upper=, name=, binary=, integer=)` | Add `n` variables, returns a `VariableVector` | +| `m.variable(lower=, upper=, name=, binary=, integer=)` | Add a single variable | +| `m.constraint_group(template)` | Add a constraint group (iterators are discovered from the template) | | `m.constraint(con)` | Add an individual constraint | | `m.objective = jp.minimize(expr)` | Set a minimization objective | | `m.objective = jp.maximize(expr)` | Set a maximization objective | +| `m.iterator(range(n))` | An index set for constraint groups | +| `m.parameter(values, name=)` | A data vector, symbolically indexable | | `m.optimize()` | Solve the model | | `m.value(var)` | Get the solved value of a variable | @@ -134,39 +136,39 @@ jp.log(x) jp.sqrt(x) jp.jp_abs(x) `VariableVector` and `Parameter` support both concrete and symbolic indexing: ```python -x[0] # concrete: returns Variable -x[i] # symbolic: returns IndexedVariable (template node) +x[0] # concrete: returns a Variable +x[i] # symbolic: a getindex template node over the block x[10*i + j] # symbolic arithmetic on the index -costs[0] # concrete: returns Constant -costs[i] # symbolic: returns IndexedParameter (template node) +costs[0] # concrete: returns a float +costs[i] # symbolic: a getindex template node over the data ``` ## Architecture JuMPy has two layers: -1. **Python package** (`jumpy`): Builds expression graphs using operator overloading. The graph maps directly to `MOI.ScalarNonlinearFunction`. Serializes models to a flat array format for FFI. +1. **Python package** (`jumpy`): Operator overloading builds MOI functions *eagerly* — every operation is one MOI call through the backend's `ops` object. There is no Python-side expression tree and no conversion step. -2. **Compiled Julia library** (built with juliac): Bundles MathOptInterface + GenOpt + Bridges + HiGHS into a shared library with a C ABI. Reconstructs expression trees, expands constraint groups, and solves. +2. **Compiled Julia library** (built with juliac): exposes the MOI API as C entry points, one per MOI call — `jumpy_scalar_nonlinear` is the compiled `MOI.ScalarNonlinearFunction` constructor, and so on. GenOpt is compiled in: templates reference iterators by identity (`GenOpt.IteratorRef`) and groups are expanded in Julia. ``` src/jumpy/ -├── expressions.py # Expression tree nodes with operator overloading -├── iterators.py # Iterator — an Expr node usable in index arithmetic -├── serialize.py # Flatten expression trees for the C ABI -└── model.py # Model class, constraint groups, solver interface +├── expressions.py # Node handles with operator overloading (eager MOI calls) +├── bridge_juliacall.py # MOI ops via juliacall +├── backend.py # Backend selection; MOI ops via ctypes (juliac) +└── model.py # Model class: variables, groups, objective, solve ``` ## How it maps to Julia | Python | Julia | |---|---| -| `Iterator(range(n))` | `GenOpt.Iterator(n, values)` | -| `x[i]` (symbolic) | `MOI.VariableIndex` resolved during expansion | +| `m.iterator(range(n))` | `GenOpt.IteratorRef(GenOpt.Iterator(values))` | +| `x[i]` (symbolic) | `getindex` node over `GenOpt.ContiguousArrayOfVariables` | | `x[i] + x[i+1] <= 10` | `MOI.ScalarNonlinearFunction` template | -| `m.constraint_group([i], ...)` | `GenOpt.IteratedFunction` | -| `Parameter([...])` | Data vector passed alongside iterators | +| `m.constraint_group(...)` | `GenOpt.FunctionGenerator` (iterators discovered by identity) | +| `m.parameter([...])` | Data vector, `getindex` resolved during expansion | ## Development @@ -176,8 +178,17 @@ python3 tests/test_expressions.py # Run example python3 examples/basic.py + +# Build the compiled backend (see julia/README.md), then: +JUMPY_BACKEND=juliac python3 tests/test_solve.py ``` +The compiled backend lives in [`julia/`](julia/): a small Julia package +(`JuMPyHiGHS`) exposing C entry points that mirror the MOI API around a raw +`HiGHS.Optimizer`, compiled with +[JuliaC](https://github.com/JuliaLang/JuliaC.jl). See +[`julia/README.md`](julia/README.md) for the C ABI and build instructions. + ## Related projects - [JuMP](https://github.com/jump-dev/JuMP.jl) — the Julia optimization modeling language diff --git a/examples/basic.py b/examples/basic.py index b1cacdf..ddd8049 100644 --- a/examples/basic.py +++ b/examples/basic.py @@ -1,8 +1,8 @@ """ Basic JuMPy usage example. -Shows the complete user-facing API. All expression graphs are built once -in Python; GeneratorOptInterface expands them in compiled Julia. +Every operation is one MOI call into the backend; constraint groups are +expanded by GenOpt in compiled Julia, not in Python. """ import sys @@ -19,44 +19,25 @@ # Instead of 99 individual constraints built in Python (slow!), # we define ONE template. GenOpt expands it in compiled Julia. -i = jp.Iterator(range(99)) -m.constraint_group([i], x[i] + x[i + 1] <= 10) - -# ── Nonlinear constraint group ──────────────────────────────────────────────── - -j = jp.Iterator(range(100)) -m.constraint_group([j], jp.sin(x[j]) + jp.exp(x[j]) <= 1.0) +i = m.iterator(range(99)) +m.constraint_group(x[i] + x[i + 1] <= 10) # ── Multi-dimensional constraint group ──────────────────────────────────────── -p = jp.Iterator(range(10)) -q = jp.Iterator(range(10)) -m.constraint_group([p, q], x[10 * p + q] >= 0) +p = m.iterator(range(10)) +q = m.iterator(range(10)) +m.constraint_group(x[10 * p + q] >= 0) # ── Constraint group with data parameters ───────────────────────────────────── -costs = jp.Parameter([float(k) * 0.5 for k in range(100)], name="costs") -k = jp.Iterator(range(100)) -m.constraint_group([k], costs[k] * x[k] <= 50) - -# ── Objective ───────────────────────────────────────────────────────────────── - -m.objective = jp.minimize(x[0] + x[1] + x[2]) - -# ── Inspect ─────────────────────────────────────────────────────────────────── +costs = m.parameter([float(k) * 0.5 + 1.0 for k in range(100)], name="costs") +k = m.iterator(range(100)) +m.constraint_group(costs[k] * x[k] <= 50) -print(f"Variables: {m._num_vars}") -print(f"Constraint groups: {len(m._constraint_groups)}") -for idx, group in enumerate(m._constraint_groups): - print(f" Group {idx}: {group}") -print(f"Objective: {m.objective}") +# ── Objective and solve ─────────────────────────────────────────────────────── -# Serialize to see the flat representation -data = m._serialize() -print(f"\nSerialized:") -print(f" {data['num_vars']} variables") -print(f" {len(data['constraint_groups'])} constraint group(s)") -print(f" {len(data['parameters'])} parameter array(s)") -print(f" objective sense: {data['objective']['sense']}") +m.objective = jp.minimize(x[0] + x[1]) +m.optimize() -# m.optimize() # ← calls compiled Julia library (not yet available) +print("x[0] =", m.value(x[0])) +print("x[1] =", m.value(x[1])) diff --git a/julia/Manifest.toml b/julia/Manifest.toml new file mode 100644 index 0000000..f1895f7 --- /dev/null +++ b/julia/Manifest.toml @@ -0,0 +1,360 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.12.6" +manifest_format = "2.0" +project_hash = "ed3c5c14e386088c6d95d3839c11b662f2ce94fb" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +version = "1.11.0" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1b96ea4a01afe0ea4090c5c8039690672dd13f2e" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.9+0" + +[[deps.CodecBzip2]] +deps = ["Bzip2_jll", "TranscodingStreams"] +git-tree-sha1 = "84990fa864b7f2b4901901ca12736e45ee79068c" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.8.5" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.8" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools"] +git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.1" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.3.0+1" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +version = "1.11.0" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "79a2aca180a85c690c58a020d47b426954b590f8" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.16.0" + +[[deps.DocStringExtensions]] +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.5" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] +git-tree-sha1 = "2c5d0b0e12088cde2cf84afb2784415b1ea3dfee" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "1.4.1" + + [deps.ForwardDiff.extensions] + ForwardDiffStaticArraysExt = "StaticArrays" + + [deps.ForwardDiff.weakdeps] + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.GenOpt]] +deps = ["JuMP", "MathOptInterface"] +git-tree-sha1 = "27eae5583786e923eecf133c8661679934392f65" +uuid = "f2c049d8-7489-4223-990c-4f1c121a4cde" +version = "0.2.1" + +[[deps.HiGHS]] +deps = ["HiGHS_jll", "LinearAlgebra", "MathOptIIS", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools", "SparseArrays"] +git-tree-sha1 = "01a5241985559c08a5baadbcebd6d87daaf84a84" +uuid = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +version = "1.24.1" + +[[deps.HiGHS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Zlib_jll", "libblastrampoline_jll"] +git-tree-sha1 = "5814a4409f49e8430c184cbe4bc19fa2957bbf0a" +uuid = "8fd58aa0-07eb-5a78-9b36-339c94fd15ea" +version = "1.15.1+0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +version = "1.11.0" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "b2d91fe939cae05960e760110b328288867b5758" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.6" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.8.0" + +[[deps.JSON]] +deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"] +git-tree-sha1 = "c89d196f5ffb64bfbf80985b699ea913b0d2c211" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "1.6.1" + + [deps.JSON.extensions] + JSONArrowExt = ["ArrowTypes"] + + [deps.JSON.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" + +[[deps.JuMP]] +deps = ["LinearAlgebra", "MacroTools", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays"] +git-tree-sha1 = "6941586d9cf3c0af718bc6e6250dcf24057d412e" +uuid = "4076af6c-e467-56ae-b986-b466b2749572" +version = "1.30.1" + + [deps.JuMP.extensions] + JuMPDimensionalDataExt = "DimensionalData" + + [deps.JuMP.weakdeps] + DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" + +[[deps.JuMPyHiGHS]] +deps = ["GenOpt", "HiGHS", "MathOptInterface"] +path = "." +uuid = "502c632e-507f-494b-8876-b49ee56148ca" +version = "0.1.0" + +[[deps.JuliaSyntaxHighlighting]] +deps = ["StyledStrings"] +uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" +version = "1.12.0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +version = "1.12.0" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "bba2d9aa057d8f126415de240573e86a8f39d2a1" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "1.0.1" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +version = "1.11.0" + +[[deps.MacroTools]] +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.16" + +[[deps.Markdown]] +deps = ["Base64", "JuliaSyntaxHighlighting", "StyledStrings"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +version = "1.11.0" + +[[deps.MathOptIIS]] +deps = ["MathOptInterface"] +git-tree-sha1 = "3b3d69130d8ab8c39d5fa4d30e20a8e6428c9d37" +uuid = "8c4f8055-bd93-4160-a86b-a0c04941dbff" +version = "0.2.0" + +[[deps.MathOptInterface]] +deps = ["CodecBzip2", "CodecZlib", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] +git-tree-sha1 = "9f23c8c1667bd0b0e611110aaf80aa91c1bdf274" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "1.51.1" + + [deps.MathOptInterface.extensions] + MathOptInterfaceBenchmarkToolsExt = "BenchmarkTools" + MathOptInterfaceCliqueTreesExt = "CliqueTrees" + + [deps.MathOptInterface.weakdeps] + BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" + CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" + +[[deps.MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "dc5b2c4c111c46bc79ac4405eeb563523b39c004" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "1.8.0" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "dbd2e8cd2c1c27f0b584f6661b4309609c5a685e" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.1.4" + +[[deps.OpenBLAS32_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "libblastrampoline_jll"] +git-tree-sha1 = "8b492aefdd20fb9dc1ebc377ff7e5fa1591c9acc" +uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" +version = "0.3.33+2" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.29+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.7+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.6+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "94ba93778373a53bfd5a0caaf7d809c445292ff4" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.2" + +[[deps.Parsers]] +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "32a4e09c5f29402573d673901778a0e03b0807b9" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.8.6" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "edbeefc7a4889f528644251bdb5fc9ab5348bc2c" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.3.4" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "8b770b60760d4451834fe79dd483e318eee709c4" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.5.2" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +version = "1.11.0" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.12.0" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "6547cbdd8ce32efba0d21c5a40fa96d1a3548f9f" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.8.0" + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + + [deps.SpecialFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.4" + +[[deps.StructUtils]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "82bee338d650aa515f31866c460cb7e3bcef90b8" +uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42" +version = "2.8.2" + + [deps.StructUtils.extensions] + StructUtilsMeasurementsExt = ["Measurements"] + StructUtilsStaticArraysCoreExt = ["StaticArraysCore"] + StructUtilsTablesExt = ["Tables"] + + [deps.StructUtils.weakdeps] + Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" + +[[deps.StyledStrings]] +uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" +version = "1.11.0" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.8.3+2" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +version = "1.11.0" + +[[deps.TranscodingStreams]] +git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.11.3" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +version = "1.11.0" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.3.1+2" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.15.0+0" diff --git a/julia/Project.toml b/julia/Project.toml new file mode 100644 index 0000000..554aedf --- /dev/null +++ b/julia/Project.toml @@ -0,0 +1,14 @@ +name = "JuMPyHiGHS" +uuid = "502c632e-507f-494b-8876-b49ee56148ca" +version = "0.1.0" + +[deps] +GenOpt = "f2c049d8-7489-4223-990c-4f1c121a4cde" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + +[compat] +GenOpt = "0.2.1" +HiGHS = "1" +MathOptInterface = "1" +julia = "1.12" diff --git a/julia/README.md b/julia/README.md new file mode 100644 index 0000000..442624d --- /dev/null +++ b/julia/README.md @@ -0,0 +1,107 @@ +# JuMPyHiGHS + +The compiled Julia backend for JuMPy: C entry points that mirror the +[MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl) API, +one MOI call per entry point. Python builds models by calling these +eagerly (`src/jumpy/expressions.py` + `model.py`), identically for the +juliacall and juliac backends; `jumpy_scalar_nonlinear` is simply the +compiled counterpart of `jl.MOI.ScalarNonlinearFunction(...)`. Nothing is +HiGHS-specific except the `Optimizer` constant in `src/JuMPyHiGHS.jl`: any +MOI optimizer can be compiled behind the same entry points. + +The optimizer is used raw — **no** `MOI.Bridges`, **no** +`CachingOptimizer`. Whatever it does not support (for HiGHS: nonlinear +functions) is reported as an error. GenOpt is compiled in: templates +reference iterators by identity (`GenOpt.IteratorRef`) and +`jumpy_add_group_constraint` expands them here with the same loop as +`GenOpt.FunctionGeneratorBridge`, one scalar constraint per combination of +iterator values. + +## C ABI + +See `src/JuMPyHiGHS.jl` for the full conventions. + +A model is an opaque `void*` pointing at the Julia-side model object; MOI +functions are opaque `void*` nodes built with the constructor entry points. +Both are rooted on the Julia side (the GC cannot see references held by C): +nodes belong to the model that built them, and everything is freed by +`jumpy_free_model`, after which no pointer from that model may be used. + +| Function | MOI equivalent | +|---|---| +| `jumpy_new_model() -> void*` | `Optimizer()` (NULL on error) | +| `jumpy_free_model(m) -> int32` | release the model and its nodes | +| `jumpy_add_variables(m, count) -> int64` | `MOI.add_variables`; returns the 0-based start index | +| `jumpy_constant(m, value) -> void*` | a `Float64` node | +| `jumpy_variable(m, index) -> void*` | `MOI.VariableIndex` of the 0-based column `index` | +| `jumpy_scalar_nonlinear(m, head, args**, nargs) -> void*` | `MOI.ScalarNonlinearFunction(Symbol(head), Any[args...])` | +| `jumpy_iterator(m, values*, len) -> void*` | `GenOpt.IteratorRef(GenOpt.Iterator(values))`, usable in template expressions | +| `jumpy_contiguous_variables(m, offset, count) -> void*` | `GenOpt.ContiguousArrayOfVariables`, 1-based-indexable block of variables | +| `jumpy_float_array(m, values*, len) -> void*` | a data vector, 1-based-indexable in templates | +| `jumpy_add_constraint(m, f, sense, rhs) -> int64` | `MOI.add_constraint(f, set)` with set `{0: LessThan, 1: GreaterThan, 2: EqualTo}(rhs)` or `{3: ZeroOne, 4: Integer}`; function constants are normalized into the set; variable bounds are just variable nodes | +| `jumpy_add_group_constraint(m, f, sense) -> int64` | expand the template over its iterators (GenOpt), one scalar constraint each; returns the count | +| `jumpy_set_objective_sense(m, sense) -> int32` | `MOI.set(MOI.ObjectiveSense())`; 0 = min, 1 = max | +| `jumpy_set_objective_function(m, f) -> int32` | `MOI.set(MOI.ObjectiveFunction{F}(), f)` | +| `jumpy_optimize(m) -> int32` | `MOI.optimize!`; returns `Int(MOI.TerminationStatusCode)` (`OPTIMAL == 1`) | +| `jumpy_primal_status(m) -> int32` | `Int(MOI.ResultStatusCode)` (`FEASIBLE_POINT == 1`) | +| `jumpy_get_values(m, out*, len) -> int64` | `MOI.VariablePrimal`; copies into `out` | +| `jumpy_objective_value(m) -> float64` | `MOI.ObjectiveValue` | + +Affine expressions built as `ScalarNonlinearFunction` trees are narrowed to +`ScalarAffineFunction` with `MOI.Nonlinear.SymbolicAD.simplify` before being +passed to the optimizer, so HiGHS accepts them. + +The consumer must initialize the Julia runtime once after loading the +library, by calling `jl_init_with_image_handle(dlopen_handle)` (see +`_load_lib` in `src/jumpy/backend.py`). + +## Building + +Requires Julia 1.12+, a C compiler, and the +[JuliaC.jl](https://github.com/JuliaLang/JuliaC.jl) frontend: + +```bash +julia --project=@juliac -e 'using Pkg; Pkg.add("JuliaC")' +``` + +Then, from this directory: + +```bash +julia --project=@juliac -m JuliaC \ + --output-lib build/libjumpy_highs \ + --compile-ccallable \ + --jl-option handle-signals=no \ + --bundle build \ + . +``` + +Notes: + +- `--jl-option handle-signals=no` is required because the library is loaded + into a Python process; Julia's signal handlers would conflict with Python's. +- `--bundle` makes the output relocatable: `build/lib/` contains + `libjumpy_highs.so` next to the Julia runtime libraries + (`build/lib/julia/`), and `build/share/julia/artifacts/` contains the + HiGHS_jll artifact with `libhighs.so`. This is the two-shared-library + layout: `libjumpy_highs.so` (our entry points + Julia runtime image) loads + `libhighs.so` (the solver distributed by HiGHS_jll) dynamically. +- No `--trim` for now: the MOI wrapper relies on dynamic dispatch that + `--trim=safe` cannot verify yet. The untrimmed library is large but + correct; trimming is an optimization to revisit. + +## Testing + +In-process tests of the entry points (fast, no compilation): + +```bash +julia --project=. test/runtests.jl +``` + +End-to-end through the compiled library and Python ctypes: + +```bash +cd .. && JUMPY_BACKEND=juliac python3 tests/test_solve.py +``` + +The Python loader searches `$JUMPY_LIB`, the installed package's `lib/` +directory, then `julia/build/lib/` (this development layout). diff --git a/julia/src/JuMPyHiGHS.jl b/julia/src/JuMPyHiGHS.jl new file mode 100644 index 0000000..c317331 --- /dev/null +++ b/julia/src/JuMPyHiGHS.jl @@ -0,0 +1,353 @@ +module JuMPyHiGHS + +# C entry points for the JuMPy juliac backend. +# +# The ABI mirrors the MOI API one function per entry point, so that the +# Python code building the model is the same for the juliacall and juliac +# backends: `jumpy_scalar_nonlinear` is the compiled counterpart of +# `jl.MOI.ScalarNonlinearFunction(...)`, `jumpy_add_constraint` of +# `MOI.Utilities.normalize_and_add_constraint(...)`, and so on. Nothing here +# is HiGHS-specific except the `Optimizer` constant below: compiling any +# other MOI optimizer behind the same entry points is a one-line change. +# +# The optimizer is used raw: no `MOI.Bridges`, no +# `MOI.Utilities.CachingOptimizer`. GenOpt is compiled in and +# jumpy_add_group_constraint expands templates here. Whatever functions and +# sets the optimizer does not support are reported as errors. +# +# Conventions across the C ABI: +# - a model is an opaque pointer returned by jumpy_new_model; it stays +# valid until jumpy_free_model, after which it must not be used +# - MOI functions are opaque pointers built with jumpy_constant / +# jumpy_variable / jumpy_scalar_nonlinear; they belong to the model +# that built them and are freed with it +# - variables are 0-based column indices in the order they were added +# - constraint sense: 0 = <=, 1 = >=, 2 = ==, 3 = binary, 4 = integer +# - objective sense: 0 = min, 1 = max +# - entry points return -1 (NULL, NaN) on error, after printing to stderr +# (a Julia exception must never propagate across the C boundary) + +import GenOpt +import HiGHS +import MathOptInterface as MOI + +# The only solver-specific line in this package. +const Optimizer = HiGHS.Optimizer + +mutable struct ModelHandle + optimizer::Optimizer + variables::Vector{MOI.VariableIndex} + # Roots the expression nodes handed out as pointers: the Julia GC + # cannot see references held by the C caller. + nodes::Vector{Base.RefValue{Any}} +end + +# Same rooting, for the models themselves: alive until jumpy_free_model. +const KEEP_ALIVE = IdDict{ModelHandle,Nothing}() +const LOCK = ReentrantLock() + +function _get(model::Ptr{Cvoid}) + model == C_NULL && error("Model pointer is NULL") + return unsafe_pointer_to_objref(model)::ModelHandle +end + +# Expression nodes (Float64, MOI.VariableIndex, MOI functions) are boxed in +# a Ref so that immutable values also get a stable pointer. +function _box(handle::ModelHandle, value)::Ptr{Cvoid} + node = Base.RefValue{Any}(value) + push!(handle.nodes, node) + return pointer_from_objref(node) +end + +function _unbox(node::Ptr{Cvoid}) + node == C_NULL && error("Expression pointer is NULL") + return (unsafe_pointer_to_objref(node)::Base.RefValue{Any})[] +end + +macro _catch(default, expr) + quote + try + $(esc(expr)) + catch err + print(stderr, "JuMPyHiGHS error: ") + showerror(stderr, err) + println(stderr) + $(esc(default)) + end + end +end + +# -- Model lifecycle ---------------------------------------------------------- + +Base.@ccallable function jumpy_new_model()::Ptr{Cvoid} + @_catch C_NULL begin + optimizer = Optimizer() + MOI.set(optimizer, MOI.Silent(), true) + handle = ModelHandle(optimizer, MOI.VariableIndex[], Base.RefValue{Any}[]) + Base.@lock LOCK KEEP_ALIVE[handle] = nothing + pointer_from_objref(handle) + end +end + +Base.@ccallable function jumpy_free_model(model::Ptr{Cvoid})::Cint + @_catch Cint(-1) begin + Base.@lock LOCK delete!(KEEP_ALIVE, _get(model)) + Cint(0) + end +end + +# -- Variables ---------------------------------------------------------------- + +# MOI.add_variables. Returns the 0-based index of the first added variable. +# Bounds are constraints: pass a variable node to jumpy_add_constraint. +Base.@ccallable function jumpy_add_variables( + model::Ptr{Cvoid}, + count::Clonglong, +)::Clonglong + @_catch Clonglong(-1) begin + handle = _get(model) + start = length(handle.variables) + append!(handle.variables, MOI.add_variables(handle.optimizer, count)) + Clonglong(start) + end +end + +# -- MOI function constructors -------------------------------------------------- + +Base.@ccallable function jumpy_constant( + model::Ptr{Cvoid}, + value::Cdouble, +)::Ptr{Cvoid} + @_catch C_NULL _box(_get(model), value) +end + +# MOI.VariableIndex of the 0-based column `index`. +Base.@ccallable function jumpy_variable( + model::Ptr{Cvoid}, + index::Clonglong, +)::Ptr{Cvoid} + @_catch C_NULL begin + handle = _get(model) + _box(handle, handle.variables[index+1]) + end +end + +# MOI.ScalarNonlinearFunction(Symbol(head), Any[args...]). +Base.@ccallable function jumpy_scalar_nonlinear( + model::Ptr{Cvoid}, + head::Cstring, + args::Ptr{Ptr{Cvoid}}, + nargs::Clonglong, +)::Ptr{Cvoid} + @_catch C_NULL begin + handle = _get(model) + func = MOI.ScalarNonlinearFunction( + Symbol(unsafe_string(head)), + Any[_unbox(unsafe_load(args, k)) for k in 1:nargs], + ) + _box(handle, func) + end +end + +# GenOpt.IteratorRef over the given values: a template node usable in +# jumpy_scalar_nonlinear args, expanded by jumpy_add_group_constraint. +Base.@ccallable function jumpy_iterator( + model::Ptr{Cvoid}, + values::Ptr{Cdouble}, + len::Clonglong, +)::Ptr{Cvoid} + @_catch C_NULL begin + handle = _get(model) + iterator = GenOpt.Iterator([unsafe_load(values, k) for k in 1:len]) + _box(handle, GenOpt.IteratorRef(iterator)) + end +end + +# GenOpt.ContiguousArrayOfVariables: the block of `count` variables starting +# at 0-based column `offset`, indexable (1-based) inside a template. +Base.@ccallable function jumpy_contiguous_variables( + model::Ptr{Cvoid}, + offset::Clonglong, + count::Clonglong, +)::Ptr{Cvoid} + @_catch C_NULL begin + handle = _get(model) + _box(handle, GenOpt.ContiguousArrayOfVariables(offset, (Int64(count),))) + end +end + +# A data vector, indexable (1-based) inside a template. +Base.@ccallable function jumpy_float_array( + model::Ptr{Cvoid}, + values::Ptr{Cdouble}, + len::Clonglong, +)::Ptr{Cvoid} + @_catch C_NULL begin + handle = _get(model) + _box(handle, [unsafe_load(values, k) for k in 1:len]) + end +end + +# -- Constraints -------------------------------------------------------------- + +# Simplify returns a ScalarAffineFunction when the expression is affine, so +# optimizers without nonlinear support (like HiGHS) accept it. Never narrow +# below ScalarAffineFunction: `x >= 0` as a constraint must stay a row, not +# become a VariableIndex bound (same semantics as JuMP's @constraint). Only +# a function that already is a VariableIndex — the bounds path — is a bound. +function _simplify(func::MOI.ScalarNonlinearFunction) + f = MOI.Nonlinear.SymbolicAD.simplify(func) + if f isa MOI.VariableIndex || f isa Float64 + return convert(MOI.ScalarAffineFunction{Float64}, f) + end + return f +end +_simplify(func) = func + +function _scalar_set(sense::Cint, rhs::Float64) + if sense == 0 + return MOI.LessThan(rhs) + elseif sense == 1 + return MOI.GreaterThan(rhs) + elseif sense == 2 + return MOI.EqualTo(rhs) + elseif sense == 3 + return MOI.ZeroOne() + elseif sense == 4 + return MOI.Integer() + end + return error("Invalid constraint sense: $sense") +end + +# MOI.add_constraint(func, set) where set is +# {0: LessThan, 1: GreaterThan, 2: EqualTo}(rhs) or {3: ZeroOne, 4: Integer} +# (rhs ignored). Function constants are +# normalized into the set. Returns the raw MOI constraint index value. +Base.@ccallable function jumpy_add_constraint( + model::Ptr{Cvoid}, + func::Ptr{Cvoid}, + sense::Cint, + rhs::Cdouble, +)::Clonglong + @_catch Clonglong(-1) begin + handle = _get(model) + ci = MOI.Utilities.normalize_and_add_constraint( + handle.optimizer, + _simplify(_unbox(func)), + _scalar_set(sense, rhs), + ) + Clonglong(ci.value) + end +end + +# Expand a template containing GenOpt.IteratorRef nodes into one scalar +# constraint `expanded(func) sense 0` per combination of iterator values — +# the same expansion loop as GenOpt.FunctionGeneratorBridge, on the raw +# optimizer. Returns the number of constraints added. +Base.@ccallable function jumpy_add_group_constraint( + model::Ptr{Cvoid}, + func::Ptr{Cvoid}, + sense::Cint, +)::Clonglong + @_catch Clonglong(-1) begin + handle = _get(model) + set = _scalar_set(sense, 0.0) + template, iterators = GenOpt.collect_iterator_refs( + _unbox(func)::MOI.ScalarNonlinearFunction, + ) + sizes = Tuple(length.(iterators)) + for idx in CartesianIndices(sizes) + values = [iterators[k].values[idx[k]] for k in eachindex(iterators)] + expanded = GenOpt._expand(template, values) + MOI.Utilities.normalize_and_add_constraint( + handle.optimizer, + _simplify(expanded), + set, + ) + end + Clonglong(prod(sizes)) + end +end + +# -- Objective ---------------------------------------------------------------- + +# MOI.set(MOI.ObjectiveSense()); 0 = min, 1 = max. +Base.@ccallable function jumpy_set_objective_sense( + model::Ptr{Cvoid}, + sense::Cint, +)::Cint + @_catch Cint(-1) begin + handle = _get(model) + moi_sense = sense == 0 ? MOI.MIN_SENSE : MOI.MAX_SENSE + MOI.set(handle.optimizer, MOI.ObjectiveSense(), moi_sense) + Cint(0) + end +end + +# MOI.set(MOI.ObjectiveFunction{F}(), func). +Base.@ccallable function jumpy_set_objective_function( + model::Ptr{Cvoid}, + func::Ptr{Cvoid}, +)::Cint + @_catch Cint(-1) begin + handle = _get(model) + f = _simplify(_unbox(func)) + if f isa MOI.VariableIndex || f isa Float64 + # HiGHS has no VariableIndex or constant objective; the affine + # one is equivalent. + f = convert(MOI.ScalarAffineFunction{Float64}, f) + end + MOI.set(handle.optimizer, MOI.ObjectiveFunction{typeof(f)}(), f) + Cint(0) + end +end + +# -- Solve and solution retrieval ---------------------------------------------- + +# Returns Int(MOI.TerminationStatusCode); MOI.OPTIMAL is 1. +Base.@ccallable function jumpy_optimize(model::Ptr{Cvoid})::Cint + @_catch Cint(-1) begin + handle = _get(model) + MOI.optimize!(handle.optimizer) + Cint(Integer(MOI.get(handle.optimizer, MOI.TerminationStatus()))) + end +end + +# Returns Int(MOI.ResultStatusCode) of the primal; MOI.FEASIBLE_POINT is 1. +Base.@ccallable function jumpy_primal_status(model::Ptr{Cvoid})::Cint + @_catch Cint(-1) begin + handle = _get(model) + Cint(Integer(MOI.get(handle.optimizer, MOI.PrimalStatus()))) + end +end + +# Writes the primal values of the first min(len, num variables) variables +# into `out`. Returns the number of values written. +Base.@ccallable function jumpy_get_values( + model::Ptr{Cvoid}, + out::Ptr{Cdouble}, + len::Clonglong, +)::Clonglong + @_catch Clonglong(-1) begin + handle = _get(model) + n = min(len, length(handle.variables)) + values = MOI.get( + handle.optimizer, + MOI.VariablePrimal(), + handle.variables[1:n], + ) + for k in 1:n + unsafe_store!(out, values[k], k) + end + Clonglong(n) + end +end + +Base.@ccallable function jumpy_objective_value(model::Ptr{Cvoid})::Cdouble + @_catch Cdouble(NaN) begin + handle = _get(model) + Cdouble(MOI.get(handle.optimizer, MOI.ObjectiveValue())) + end +end + +end # module JuMPyHiGHS diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl new file mode 100644 index 0000000..86a219b --- /dev/null +++ b/julia/test/runtests.jl @@ -0,0 +1,149 @@ +# Tests the C entry points in-process (no compilation needed). + +import JuMPyHiGHS +using Test + +variable(m, i) = JuMPyHiGHS.jumpy_variable(m, Clonglong(i)) +constant(m, v) = JuMPyHiGHS.jumpy_constant(m, Cdouble(v)) + +function snf(m, head::String, args...) + argv = collect(Ptr{Cvoid}, args) + GC.@preserve head argv JuMPyHiGHS.jumpy_scalar_nonlinear( + m, + Cstring(pointer(head)), + pointer(argv), + Clonglong(length(argv)), + ) +end + +add_constraint(m, f, sense, rhs) = + JuMPyHiGHS.jumpy_add_constraint(m, f, Cint(sense), Cdouble(rhs)) + +function set_objective(m, sense, f) + JuMPyHiGHS.jumpy_set_objective_sense(m, Cint(sense)) == 0 && + JuMPyHiGHS.jumpy_set_objective_function(m, f) == 0 +end + +function get_values(m, n::Int) + out = zeros(Cdouble, n) + written = GC.@preserve out JuMPyHiGHS.jumpy_get_values( + m, pointer(out), Clonglong(n), + ) + @test written == n + return out +end + +# min x + y s.t. x + y >= 10, x, y >= 0 +@testset "simple LP" begin + m = JuMPyHiGHS.jumpy_new_model() + @test m != C_NULL + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(2)) == 0 + # bounds are VariableIndex-in-GreaterThan constraints, as in MOI + @test add_constraint(m, variable(m, 0), 1, 0.0) >= 0 + @test add_constraint(m, variable(m, 1), 1, 0.0) >= 0 + f = snf(m, "+", variable(m, 0), variable(m, 1)) + @test f != C_NULL + @test add_constraint(m, f, 1, 10.0) >= 0 + @test set_objective(m, 0, f) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 # MOI.OPTIMAL + @test JuMPyHiGHS.jumpy_primal_status(m) == 1 # MOI.FEASIBLE_POINT + out = get_values(m, 2) + @test out[1] + out[2] ≈ 10.0 atol = 1e-6 + @test JuMPyHiGHS.jumpy_objective_value(m) ≈ 10.0 atol = 1e-6 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +# max x s.t. 0 <= x <= 42 +@testset "maximize with bounds" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(1)) == 0 + @test add_constraint(m, variable(m, 0), 1, 0.0) >= 0 + @test add_constraint(m, variable(m, 0), 0, 42.0) >= 0 + @test set_objective(m, 1, variable(m, 0)) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 + @test get_values(m, 1)[1] ≈ 42.0 atol = 1e-6 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +# equality with a function constant, normalized into the set: +# x + 1 == 5 => x == 4 +# also exercises simplification: 2 * (x + 3) <= 14 => x <= 4 +@testset "constant normalization and simplification" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(1)) == 0 + x = variable(m, 0) + @test add_constraint(m, snf(m, "+", x, constant(m, 1.0)), 2, 5.0) >= 0 + scaled = snf(m, "*", constant(m, 2.0), snf(m, "+", x, constant(m, 3.0))) + @test add_constraint(m, scaled, 0, 14.0) >= 0 + @test set_objective(m, 0, x) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 + @test get_values(m, 1)[1] ≈ 4.0 atol = 1e-6 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +# constraint group: x[i] >= demand[i] for i in 0..2, demand = (1, 2, 3) +# min sum(x) => x = (1, 2, 3) +@testset "constraint group" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(3)) == 0 + values = Cdouble[0.0, 1.0, 2.0] # Python-style 0-based iterator values + demand = Cdouble[1.0, 2.0, 3.0] + i = GC.@preserve values JuMPyHiGHS.jumpy_iterator(m, pointer(values), Clonglong(3)) + x = JuMPyHiGHS.jumpy_contiguous_variables(m, Clonglong(0), Clonglong(3)) + d = GC.@preserve demand JuMPyHiGHS.jumpy_float_array(m, pointer(demand), Clonglong(3)) + @test i != C_NULL && x != C_NULL && d != C_NULL + i1 = snf(m, "+", i, constant(m, 1.0)) # 0-based -> 1-based + template = snf(m, "-", snf(m, "getindex", x, i1), snf(m, "getindex", d, i1)) + n = JuMPyHiGHS.jumpy_add_group_constraint(m, template, Cint(1)) + @test n == 3 + obj = snf(m, "+", variable(m, 0), snf(m, "+", variable(m, 1), variable(m, 2))) + @test set_objective(m, 0, obj) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 + @test get_values(m, 3) ≈ [1.0, 2.0, 3.0] atol = 1e-6 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +# binary knapsack: max 2 x0 + x1 s.t. x0 + x1 <= 1, x binary => (1, 0) +@testset "binary variables" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(2)) == 0 + @test add_constraint(m, variable(m, 0), 3, 0.0) >= 0 # ZeroOne + @test add_constraint(m, variable(m, 1), 3, 0.0) >= 0 + @test add_constraint(m, snf(m, "+", variable(m, 0), variable(m, 1)), 0, 1.0) >= 0 + obj = snf(m, "+", snf(m, "*", constant(m, 2.0), variable(m, 0)), variable(m, 1)) + @test set_objective(m, 1, obj) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 + @test get_values(m, 2) ≈ [1.0, 0.0] atol = 1e-6 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +# `x >= 0` as a *constraint* is a ScalarAffineFunction row (like JuMP's +# @constraint), so it must not clash with an existing variable bound. +@testset "constraint on bounded variable is a row, not a bound" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(1)) == 0 + x = variable(m, 0) + @test add_constraint(m, x, 1, 0.0) >= 0 # bound: raw VariableIndex + # constraint: x - 0 >= 0, simplifies to an affine row, no bound clash + @test add_constraint(m, snf(m, "-", x, constant(m, 0.0)), 1, 0.0) >= 0 + @test set_objective(m, 0, x) + @test JuMPyHiGHS.jumpy_optimize(m) == 1 + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +@testset "unsupported inputs return -1" begin + m = JuMPyHiGHS.jumpy_new_model() + @test JuMPyHiGHS.jumpy_add_variables(m, Clonglong(1)) == 0 + redirect_stderr(devnull) do + # sin(x) <= 1 stays a ScalarNonlinearFunction, which the raw + # HiGHS.Optimizer does not support + @test add_constraint(m, snf(m, "sin", variable(m, 0)), 0, 1.0) == -1 + # NULL pointers (any other invalid pointer is undefined behavior, + # as in any C API) + @test JuMPyHiGHS.jumpy_optimize(C_NULL) == -1 + @test isnan(JuMPyHiGHS.jumpy_objective_value(C_NULL)) + end + @test JuMPyHiGHS.jumpy_free_model(m) == 0 +end + +println("All JuMPyHiGHS tests passed.") diff --git a/src/jumpy/__init__.py b/src/jumpy/__init__.py index 17ffff1..06f3e21 100644 --- a/src/jumpy/__init__.py +++ b/src/jumpy/__init__.py @@ -1,38 +1,32 @@ """ -JuMPy: A Python interface to MathOptInterface via GeneratorOptInterface. +JuMPy: A Python interface to MathOptInterface via GenOpt. -Builds expression graphs in Python, hands them off to a compiled Julia library -(MOI + GeneratorOptInterface + Bridges + HiGHS) for constraint expansion and solving. +Models are built eagerly: every operation performs the corresponding MOI +call, either in the compiled Julia library (juliac backend, no Julia +installation needed) or through juliacall. """ from jumpy.expressions import ( - Variable, - VariableVector, - Constant, - Parameter, - Expr, - Func, Constraint, + Node, Objective, + Parameter, + Variable, + VariableVector, ) from jumpy.expressions import sin, cos, exp, log, sqrt, abs as jp_abs -from jumpy.iterators import Iterator -from jumpy.model import Model, minimize, maximize, sum_over +from jumpy.model import Model, minimize, maximize __all__ = [ "Model", + "Node", "Variable", "VariableVector", - "Constant", "Parameter", - "Expr", - "Func", "Constraint", "Objective", - "Iterator", "minimize", "maximize", - "sum_over", "sin", "cos", "exp", diff --git a/src/jumpy/backend.py b/src/jumpy/backend.py index dfff732..0ec2eca 100644 --- a/src/jumpy/backend.py +++ b/src/jumpy/backend.py @@ -1,140 +1,202 @@ """ -Solver backend abstraction. +Backend selection and the compiled-library (juliac) MOI ops. -Two backends: - - "juliac" (default): calls a precompiled shared library via ctypes. - No Julia installation required. - - "juliacall": calls MOI + GenOpt + HiGHS through juliacall. - Requires Julia (installed lazily by juliacall on first use). +An "ops" object maps each MOI call either to the compiled shared library +(JuliacOps below, via ctypes) or to Julia through juliacall +(jumpy.bridge_juliacall.JuliaCallOps). Models call the ops directly; the +two implementations expose the same methods. """ from __future__ import annotations -from abc import ABC, abstractmethod -from typing import TYPE_CHECKING +import ctypes -if TYPE_CHECKING: - from jumpy.model import Model +def get_ops(backend: str): + if backend == "juliac": + return JuliacOps(_load_lib()) + if backend == "juliacall": + from jumpy.bridge_juliacall import JuliaCallOps -class Backend(ABC): - """Abstract solver backend.""" + return JuliaCallOps() + raise ValueError(f"Unknown backend '{backend}'. Choose from: juliac, juliacall") - @abstractmethod - def optimize(self, model: Model) -> list[float]: - """Solve the model and return the solution vector.""" - ... +# The Julia runtime can only be initialized once per process. +_LIB = None -class JuliacBackend(Backend): - """ - Default backend: calls a precompiled Julia shared library via ctypes. - The library is built with juliac from: - MOI + GenOpt + Bridges + HiGHS +def _load_lib(): + global _LIB + if _LIB is not None: + return _LIB + import os + import platform - No Julia installation required. - """ + soext = {"Linux": ".so", "Darwin": ".dylib", "Windows": ".dll"}[platform.system()] + lib_name = "libjumpy_highs" + soext + + candidates = [] + if "JUMPY_LIB" in os.environ: + candidates.append(os.environ["JUMPY_LIB"]) + pkg_dir = os.path.dirname(__file__) + # Wheel layout: shipped inside the package. + candidates.append(os.path.join(pkg_dir, "lib", lib_name)) + # Development layout: JuliaC bundle in /julia/build. + repo = os.path.dirname(os.path.dirname(pkg_dir)) + candidates.append(os.path.join(repo, "julia", "build", "lib", lib_name)) - def __init__(self): - self._lib = None - - def _load_lib(self): - if self._lib is not None: - return - import ctypes - import importlib.resources - # TODO: resolve platform-specific library path - # For now, search standard locations - import os - lib_names = [ - "libjumpy_backend.so", - "libjumpy_backend.dylib", - "jumpy_backend.dll", - ] - for name in lib_names: - for search_dir in [os.path.dirname(__file__), os.getcwd(), "/usr/local/lib"]: - path = os.path.join(search_dir, name) - if os.path.exists(path): - self._lib = ctypes.CDLL(path) - return + path = next((p for p in candidates if os.path.exists(p)), None) + if path is None: raise FileNotFoundError( - "Could not find the compiled JuMPy backend library.\n" - "The juliac-compiled shared library (libjumpy_backend.so) is not installed.\n" - "Either:\n" + f"Could not find the compiled JuMPy backend library ({lib_name}). " + "Searched:\n " + "\n ".join(candidates) + "\nEither:\n" " 1. Install the pre-built wheel: pip install jumpy\n" - " 2. Use the juliacall backend: jp.Model(backend='juliacall')\n" + " 2. Build it locally: see julia/README.md\n" + " 3. Use the juliacall backend: jp.Model(backend='juliacall')\n" ) - def optimize(self, model: Model) -> list[float]: - self._load_lib() - data = model._serialize() - # TODO: implement ctypes calls to the compiled library - raise NotImplementedError( - "juliac backend not yet compiled. " - "Use jp.Model(backend='juliacall') for now." - ) + # RTLD_GLOBAL so that libjulia symbols are visible process-wide, + # which the Julia runtime requires. + lib = ctypes.CDLL(path, mode=ctypes.RTLD_GLOBAL) + + # Initialize the Julia runtime from the image embedded in the library. + init = lib.jl_init_with_image_handle + init.argtypes = [ctypes.c_void_p] + init.restype = None + init(lib._handle) + + c_longlong = ctypes.c_longlong + c_int = ctypes.c_int + c_double = ctypes.c_double + c_void_p = ctypes.c_void_p + p_double = ctypes.POINTER(c_double) + p_void = ctypes.POINTER(c_void_p) + + # A model is an opaque pointer to a Julia object, valid until + # jumpy_free_model. MOI functions are opaque pointers built with the + # constructor entry points; they belong to the model and are freed + # with it. + lib.jumpy_new_model.argtypes = [] + lib.jumpy_new_model.restype = c_void_p + lib.jumpy_free_model.argtypes = [c_void_p] + lib.jumpy_free_model.restype = c_int + lib.jumpy_add_variables.argtypes = [c_void_p, c_longlong] + lib.jumpy_add_variables.restype = c_longlong + lib.jumpy_constant.argtypes = [c_void_p, c_double] + lib.jumpy_constant.restype = c_void_p + lib.jumpy_variable.argtypes = [c_void_p, c_longlong] + lib.jumpy_variable.restype = c_void_p + lib.jumpy_scalar_nonlinear.argtypes = [c_void_p, ctypes.c_char_p, p_void, c_longlong] + lib.jumpy_scalar_nonlinear.restype = c_void_p + lib.jumpy_iterator.argtypes = [c_void_p, p_double, c_longlong] + lib.jumpy_iterator.restype = c_void_p + lib.jumpy_contiguous_variables.argtypes = [c_void_p, c_longlong, c_longlong] + lib.jumpy_contiguous_variables.restype = c_void_p + lib.jumpy_float_array.argtypes = [c_void_p, p_double, c_longlong] + lib.jumpy_float_array.restype = c_void_p + lib.jumpy_add_constraint.argtypes = [c_void_p, c_void_p, c_int, c_double] + lib.jumpy_add_constraint.restype = c_longlong + lib.jumpy_add_group_constraint.argtypes = [c_void_p, c_void_p, c_int] + lib.jumpy_add_group_constraint.restype = c_longlong + lib.jumpy_set_objective_sense.argtypes = [c_void_p, c_int] + lib.jumpy_set_objective_sense.restype = c_int + lib.jumpy_set_objective_function.argtypes = [c_void_p, c_void_p] + lib.jumpy_set_objective_function.restype = c_int + lib.jumpy_optimize.argtypes = [c_void_p] + lib.jumpy_optimize.restype = c_int + lib.jumpy_primal_status.argtypes = [c_void_p] + lib.jumpy_primal_status.restype = c_int + lib.jumpy_get_values.argtypes = [c_void_p, p_double, c_longlong] + lib.jumpy_get_values.restype = c_longlong + lib.jumpy_objective_value.argtypes = [c_void_p] + lib.jumpy_objective_value.restype = c_double + + _LIB = lib + return lib + + +# Set codes of jumpy_add_constraint: {0: LessThan, 1: GreaterThan, +# 2: EqualTo}(rhs), {3: ZeroOne, 4: Integer} (rhs ignored). +_SENSE_CODES = {"<=": 0, ">=": 1, "==": 2, "binary": 3, "integer": 4} + + +class JuliacOps: + """ + The compiled-library implementation of the MOI ops. Each method is one + C call into the entry point wrapping the same MOI function the + juliacall ops call. + """ + def __init__(self, lib): + self._lib = lib + self._m = lib.jumpy_new_model() + if not self._m: # NULL + raise RuntimeError("Failed to create model") -class JuliaCallBackend(Backend): - """ - Optional backend: calls Julia directly through juliacall. + def free(self): + self._lib.jumpy_free_model(self._m) - Requires `pip install jumpy[juliacall]`. Julia is installed lazily - by juliacall on first use if not already present. + def _node(self, node): + if not node: # NULL + raise RuntimeError("Failed to build MOI function") + return node - This backend has full flexibility — it can use any solver or MOI - feature, not just what's compiled into the juliac library. - """ + # -- MOI functions --------------------------------------------------------- + + def constant(self, value): + return self._node(self._lib.jumpy_constant(self._m, value)) - def __init__(self): - self._jl = None - - def _init_julia(self): - if self._jl is not None: - return - try: - from juliacall import Main as jl - except ImportError: - raise ImportError( - "juliacall is not installed.\n" - "Install it with: pip install jumpy[juliacall]\n" - "This will also install Julia automatically if needed." - ) from None - # Install and load Julia packages on first use - jl.seval("using Pkg") - for pkg in ["MathOptInterface", "HiGHS", "GenOpt"]: - jl.seval(f""" - if !haskey(Pkg.project().dependencies, "{pkg}") - Pkg.add("{pkg}") - end - """) - jl.seval("import MathOptInterface as MOI") - jl.seval("import GenOpt") - jl.seval("import HiGHS") - # TODO: load GenOpt once it's registered / available - self._jl = jl - - def optimize(self, model: Model) -> list[float]: - self._init_julia() - jl = self._jl - return self._build_and_solve(jl, model) - - def _build_and_solve(self, jl, model: Model) -> list[float]: - from jumpy.bridge_juliacall import build_moi_model - return build_moi_model(jl, model) - - -_BACKENDS = { - "juliac": JuliacBackend, - "juliacall": JuliaCallBackend, -} - - -def get_backend(name: str) -> Backend: - cls = _BACKENDS.get(name) - if cls is None: - raise ValueError( - f"Unknown backend '{name}'. Choose from: {list(_BACKENDS.keys())}" + def variable(self, index): + return self._node(self._lib.jumpy_variable(self._m, index)) + + def scalar_nonlinear(self, head, args): + argv = (ctypes.c_void_p * len(args))(*args) + return self._node( + self._lib.jumpy_scalar_nonlinear(self._m, head.encode(), argv, len(args)) ) - return cls() + + def iterator(self, values): + data = (ctypes.c_double * len(values))(*values) + return self._node(self._lib.jumpy_iterator(self._m, data, len(values))) + + def contiguous_variables(self, start, count): + return self._node(self._lib.jumpy_contiguous_variables(self._m, start, count)) + + def float_array(self, values): + data = (ctypes.c_double * len(values))(*values) + return self._node(self._lib.jumpy_float_array(self._m, data, len(values))) + + # -- Model building ---------------------------------------------------------- + + def add_variables(self, count): + start = self._lib.jumpy_add_variables(self._m, count) + if start < 0: + raise RuntimeError("Failed to add variables") + return start + + def add_constraint(self, func, sense, rhs): + ci = self._lib.jumpy_add_constraint(self._m, func, _SENSE_CODES[sense], rhs) + if ci < 0: + raise RuntimeError("Failed to add constraint") + + def add_constraint_group(self, func, sense, linear): + n = self._lib.jumpy_add_group_constraint(self._m, func, _SENSE_CODES[sense]) + if n < 0: + raise RuntimeError("Failed to add constraint group") + + def set_objective(self, sense, func): + if self._lib.jumpy_set_objective_sense(self._m, 0 if sense == "min" else 1) != 0: + raise RuntimeError("Failed to set objective sense") + if self._lib.jumpy_set_objective_function(self._m, func) != 0: + raise RuntimeError("Failed to set objective function") + + def optimize(self): + return self._lib.jumpy_optimize(self._m) + + def get_values(self, count): + out = (ctypes.c_double * count)() + written = self._lib.jumpy_get_values(self._m, out, count) + if written != count: + raise RuntimeError(f"Expected {count} solution values, got {written}") + return list(out) diff --git a/src/jumpy/bridge_juliacall.py b/src/jumpy/bridge_juliacall.py index 16992f7..83b3033 100644 --- a/src/jumpy/bridge_juliacall.py +++ b/src/jumpy/bridge_juliacall.py @@ -1,390 +1,165 @@ """ -Bridge between JuMPy's Python expression graph and MathOptInterface via juliacall. +The juliacall implementation of the MOI ops. -Translates Python Expr nodes into Julia MOI + GenOpt types. -This module is only imported when backend="juliacall" is used. - -The Python side does NO iteration over constraints or variables. -It builds one expression template per constraint group, hands it to GenOpt -as a FunctionGenerator, and lets Julia handle all expansion. +Each method is one MOI call through juliacall. The compiled (juliac) +backend implements the same ops against the C entry points of the shared +library (jumpy.backend.JuliacOps). """ from __future__ import annotations -from typing import TYPE_CHECKING - -if TYPE_CHECKING: - from jumpy.model import Model - -from jumpy.expressions import ( - BinaryOp, - Constant, - Func, - IndexedParameter, - IndexedVariable, - UnaryOp, - Variable, -) -from jumpy.iterators import Iterator - - -_HELPERS_DEFINED = False -_any_vec = None - - -def _define_helpers(jl): - """Define Julia helper functions once.""" - global _HELPERS_DEFINED, _any_vec - if _HELPERS_DEFINED: - return - # jl.Any[...] is broken in PythonCall with Julia 1.12+ - _any_vec = jl.seval("(args...) -> Any[args...]") - jl.seval(""" - function _jumpy_add_variables!(optimizer, count, lower, upper, binary, integer) - vars = MOI.add_variables(optimizer, count) - if !isnothing(lower) - for v in vars - MOI.add_constraint(optimizer, v, MOI.GreaterThan(lower)) - end - end - if !isnothing(upper) - for v in vars - MOI.add_constraint(optimizer, v, MOI.LessThan(upper)) - end - end - if binary - for v in vars - MOI.add_constraint(optimizer, v, MOI.ZeroOne()) - end - elseif integer - for v in vars - MOI.add_constraint(optimizer, v, MOI.Integer()) - end +_JL = None + + +def _julia(): + global _JL + if _JL is not None: + return _JL + try: + from juliacall import Main as jl + except ImportError: + raise ImportError( + "juliacall is not installed.\n" + "Install it with: pip install jumpy[juliacall]\n" + "This will also install Julia automatically if needed." + ) from None + # Install and load Julia packages on first use + jl.seval("using Pkg") + for pkg in ["MathOptInterface", "HiGHS", "GenOpt"]: + jl.seval(f""" + if !haskey(Pkg.project().dependencies, "{pkg}") + Pkg.add("{pkg}") end - return vars - end - - function _jumpy_add_constraint_group!(optimizer, func, sense) - n = prod(length.(func.iterators)) - if sense == "<=" - set = MOI.Nonpositives(n) - elseif sense == ">=" - set = MOI.Nonnegatives(n) - elseif sense == "==" - set = MOI.Zeros(n) - else - error("Unknown sense: $sense") + """) + jl.seval("import MathOptInterface as MOI") + jl.seval("import GenOpt") + jl.seval("import HiGHS") + _JL = jl + return jl + + +class JuliaCallOps: + def __init__(self): + jl = _julia() + self._jl = jl + # jl.Any[...] is broken in PythonCall with Julia 1.12+ + self._any_vec = jl.seval("(args...) -> Any[args...]") + # {} type application is not expressible in Python syntax + self._objective_attr = jl.seval("f -> MOI.ObjectiveFunction{typeof(f)}()") + self._to_affine = jl.seval("f -> convert(MOI.ScalarAffineFunction{Float64}, f)") + self._optimizer = jl.seval(""" + let optimizer = MOI.instantiate( + MOI.OptimizerWithAttributes(HiGHS.Optimizer, "output_flag" => false), + with_bridge_type = Float64, + ) + MOI.Bridges.add_bridge(optimizer, GenOpt.FunctionGeneratorBridge{Float64}) + optimizer end - MOI.add_constraint(optimizer, func, set) - end - - function _jumpy_make_generator(func, iterators, target_type) - return GenOpt.FunctionGenerator{target_type}(func, iterators) - end - - function _jumpy_create_optimizer() - optimizer = MOI.instantiate( - MOI.OptimizerWithAttributes(HiGHS.Optimizer, "output_flag" => false), - with_bridge_type = Float64, - ) - MOI.Bridges.add_bridge(optimizer, GenOpt.FunctionGeneratorBridge{Float64}) - return optimizer - end + """) + self._variables = [] - function _jumpy_get_solution(optimizer, vars) - return [MOI.get(optimizer, MOI.VariablePrimal(), v) for v in vars] - end - """) - _HELPERS_DEFINED = True + def free(self): + pass # the optimizer is garbage-collected with this object + # -- MOI functions --------------------------------------------------------- -def build_moi_model(jl, model: Model) -> list[float]: - """ - Build an MOI model in Julia from a JuMPy Model and solve it. + def constant(self, value): + return value - Constraint groups are passed as GenOpt.FunctionGenerator objects - so that expansion happens entirely in Julia. - """ + def variable(self, index): + return self._variables[index] - _define_helpers(jl) + def scalar_nonlinear(self, head, args): + jl = self._jl + return jl.MOI.ScalarNonlinearFunction(jl.Symbol(head), self._any_vec(*args)) - optimizer = jl._jumpy_create_optimizer() + def iterator(self, values): + jl = self._jl + return jl.GenOpt.IteratorRef(jl.GenOpt.Iterator(jl.seval("collect")(values))) - # Add variables — one bulk call per block - all_jl_vars = [] - for block in model._var_blocks: - lower = float(block.lower) if block.lower is not None else jl.nothing - upper = float(block.upper) if block.upper is not None else jl.nothing - block_vars = jl._jumpy_add_variables_b( - optimizer, block.count, lower, upper, - block.binary, - block.integer, + def contiguous_variables(self, start, count): + return self._jl.seval( + f"GenOpt.ContiguousArrayOfVariables({start}, ({count},))" ) - all_jl_vars.append(block_vars) - # Add constraint groups via GenOpt - for group in model._constraint_groups: - _add_constraint_group(jl, optimizer, all_jl_vars, model, group) - - # Add individual constraints - for con in model._individual_constraints: - _add_individual_constraint(jl, optimizer, all_jl_vars, con) - - # Set objective - if model._objective is not None: - sense = ( - jl.MOI.MIN_SENSE - if model._objective.sense == "min" - else jl.MOI.MAX_SENSE + def float_array(self, values): + return self._jl.seval("collect")(values) + + def _simplify(self, func): + """ + Narrow an affine ScalarNonlinearFunction to ScalarAffineFunction — + but never below: `x >= 0` as a constraint must stay a row, not + become a VariableIndex bound (same semantics as JuMP's @constraint). + Only a function that already is a VariableIndex — the bounds path — + is a bound. + """ + jl = self._jl + if jl.isa(func, jl.MOI.ScalarNonlinearFunction): + func = jl.MOI.Nonlinear.SymbolicAD.simplify(func) + if jl.isa(func, jl.MOI.VariableIndex) or isinstance(func, float): + func = self._to_affine(func) + return func + + # -- Model building ---------------------------------------------------------- + + def add_variables(self, count): + start = len(self._variables) + self._variables.extend(self._jl.MOI.add_variables(self._optimizer, count)) + return start + + def _set(self, sense, rhs): + jl = self._jl + if sense == "<=": + return jl.MOI.LessThan(float(rhs)) + if sense == ">=": + return jl.MOI.GreaterThan(float(rhs)) + if sense == "==": + return jl.MOI.EqualTo(float(rhs)) + if sense == "binary": + return jl.MOI.ZeroOne() + if sense == "integer": + return jl.MOI.Integer() + raise ValueError(f"Unknown constraint sense: {sense}") + + def add_constraint(self, func, sense, rhs): + self._jl.MOI.Utilities.normalize_and_add_constraint( + self._optimizer, self._simplify(func), self._set(sense, rhs), ) - obj_func = _expr_to_moi(jl, all_jl_vars, model._objective.expr) - - jl.MOI.set(optimizer, jl.MOI.ObjectiveSense(), sense) - jl.MOI.set(optimizer, jl.MOI.ObjectiveFunction[jl.typeof(obj_func)](), obj_func) - - # Optimize and extract solution - jl.MOI.optimize_b(optimizer) - - # Flatten all variable blocks into one solution vector - all_vars_flat = jl.seval("vcat")(*(v for v in all_jl_vars)) - jl_solution = jl._jumpy_get_solution(optimizer, all_vars_flat) - return [float(jl_solution[i]) for i in range(len(jl_solution))] - - -def _is_linear_template(expr) -> bool: - """Check if a template expression is linear (no nonlinear functions).""" - match expr: - case Constant() | Variable() | Iterator() | IndexedVariable() | IndexedParameter(): - return True - case BinaryOp(op=op, left=left, right=right): - if op in ("+", "-", "*"): - return _is_linear_template(left) and _is_linear_template(right) - return False - case UnaryOp(op="-", arg=arg): - return _is_linear_template(arg) - case Func(): - return False - case _: - return False - - -def _add_constraint_group(jl, optimizer, all_jl_vars, model, group): - """ - Add a constraint group as a single GenOpt.FunctionGenerator. - - Python builds the template expression and iterator list, then hands - them to GenOpt. No Python-side iteration over constraint instances. - """ - # Build GenOpt iterators - genopt_iterators = jl.seval("GenOpt.Iterator[]") - iter_id_map = {} - for idx, it in enumerate(group.iterators): - jl_values = jl.seval("collect")(it.values) - jl_it = jl.GenOpt.Iterator(jl_values) - jl.push_b(genopt_iterators, jl_it) - iter_id_map[it.id] = idx + 1 # 1-based - - # Normalize: lhs - rhs in {Nonpositives, Nonnegatives, Zeros} - normalized = group.template.lhs - group.template.rhs - - # Build MOI.ScalarNonlinearFunction template with GenOpt placeholders - template_func = _expr_to_moi_template( - jl, all_jl_vars, model, normalized, genopt_iterators, iter_id_map, - ) - - # Determine target function type: affine if template is linear, else nonlinear - if _is_linear_template(group.template.lhs) and _is_linear_template(group.template.rhs): - target_type = jl.seval("MOI.ScalarAffineFunction{Float64}") - else: - target_type = jl.seval("MOI.ScalarNonlinearFunction") - - # Wrap in FunctionGenerator and add constraint — all in Julia - func_gen = jl._jumpy_make_generator(template_func, genopt_iterators, target_type) - jl._jumpy_add_constraint_group_b(optimizer, func_gen, group.template.sense) - - -def _get_jl_var(jl, all_jl_vars, var_index, model): - """Get the Julia MOI.VariableIndex for a Python Variable by its index.""" - offset = 0 - for block_idx, block in enumerate(model._var_blocks): - if var_index < offset + block.count: - local_idx = var_index - offset - return all_jl_vars[block_idx][local_idx] # PythonCall uses 0-based indexing - offset += block.count - raise IndexError(f"Variable index {var_index} out of range") - - -def _get_contiguous(jl, all_jl_vars, variable_vector, model): - """Get a GenOpt.ContiguousArrayOfVariables for a VariableVector.""" - start = variable_vector._variables[0].index - count = len(variable_vector) - return jl.seval( - f"GenOpt.ContiguousArrayOfVariables({start}, ({count},))" - ) - - -def _expr_to_moi_template(jl, all_jl_vars, model, expr, genopt_iterators, iter_id_map): - """ - Convert a Python Expr into an MOI.ScalarNonlinearFunction template - with GenOpt.IteratorIndex and ContiguousArrayOfVariables placeholders. - """ - match expr: - case Constant(value=v): - return v - case Variable(index=idx): - return _get_jl_var(jl, all_jl_vars, idx, model) - case BinaryOp(op=op, left=left, right=right): - l = _expr_to_moi_template(jl, all_jl_vars, model, left, genopt_iterators, iter_id_map) - r = _expr_to_moi_template(jl, all_jl_vars, model, right, genopt_iterators, iter_id_map) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol(op), _any_vec(l, r)) - case UnaryOp(op="-", arg=arg): - a = _expr_to_moi_template(jl, all_jl_vars, model, arg, genopt_iterators, iter_id_map) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol("-"), _any_vec(a)) - case Func(name=name, arg=arg): - a = _expr_to_moi_template(jl, all_jl_vars, model, arg, genopt_iterators, iter_id_map) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol(name), _any_vec(a)) - case Iterator() as it: - jl_idx = iter_id_map[it.id] - return jl.GenOpt.IteratorIndex(jl_idx) - case IndexedVariable() as iv: - contiguous = _get_contiguous(jl, all_jl_vars, iv.variable_vector, model) - index_expr = _expr_to_moi_template( - jl, all_jl_vars, model, iv.index_expr, genopt_iterators, iter_id_map, - ) - # 0-based Python → 1-based Julia - index_1based = jl.MOI.ScalarNonlinearFunction( - jl.Symbol("+"), _any_vec(index_expr, 1), - ) - return jl.MOI.ScalarNonlinearFunction( - jl.Symbol("getindex"), _any_vec(contiguous, index_1based), - ) - case IndexedParameter() as ip: - jl_values = jl.seval("collect")(ip.parameter.values) - index_expr = _expr_to_moi_template( - jl, all_jl_vars, model, ip.index_expr, genopt_iterators, iter_id_map, - ) - index_1based = jl.MOI.ScalarNonlinearFunction( - jl.Symbol("+"), _any_vec(index_expr, 1), - ) - return jl.MOI.ScalarNonlinearFunction( - jl.Symbol("getindex"), _any_vec(jl_values, index_1based), - ) - case _: - raise TypeError(f"Cannot convert {type(expr).__name__} to MOI template") - - -def _get_jl_var_by_index(jl, all_jl_vars, idx): - """Get Julia MOI.VariableIndex for a Python variable by global index.""" - offset = 0 - for block_idx, block_vars in enumerate(all_jl_vars): - block_size = int(jl.length(block_vars)) - if idx < offset + block_size: - return block_vars[idx - offset] # PythonCall uses 0-based indexing - offset += block_size - raise IndexError(f"Variable index {idx} out of range") - - -def _collect_linear_terms(expr, terms, sign=1.0): - """ - Try to decompose expr into linear terms: list of (coef, var_index) + constant. - Returns (success, constant). - """ - match expr: - case Constant(value=v): - return True, v * sign - case Variable(index=idx): - terms.append((sign, idx)) - return True, 0.0 - case BinaryOp(op="+", left=left, right=right): - terms_before = len(terms) - ok_l, const_l = _collect_linear_terms(left, terms, sign) - if not ok_l: - del terms[terms_before:] - return False, 0.0 - ok_r, const_r = _collect_linear_terms(right, terms, sign) - if not ok_r: - del terms[terms_before:] - return False, 0.0 - return True, const_l + const_r - case BinaryOp(op="-", left=left, right=right): - terms_before = len(terms) - ok_l, const_l = _collect_linear_terms(left, terms, sign) - if not ok_l: - del terms[terms_before:] - return False, 0.0 - ok_r, const_r = _collect_linear_terms(right, terms, -sign) - if not ok_r: - del terms[terms_before:] - return False, 0.0 - return True, const_l + const_r - case BinaryOp(op="*", left=Constant(value=v), right=right): - return _collect_linear_terms(right, terms, sign * v) - case BinaryOp(op="*", left=left, right=Constant(value=v)): - return _collect_linear_terms(left, terms, sign * v) - case UnaryOp(op="-", arg=arg): - return _collect_linear_terms(arg, terms, -sign) - case _: - return False, 0.0 - - -def _expr_to_moi_linear(jl, all_jl_vars, expr): - """ - Try to convert expr to ScalarAffineFunction. Returns None if nonlinear. - """ - terms = [] - ok, constant = _collect_linear_terms(expr, terms) - if not ok: - return None - - jl_terms = jl.seval("MOI.ScalarAffineTerm{Float64}[]") - for coef, var_idx in terms: - jl_var = _get_jl_var_by_index(jl, all_jl_vars, var_idx) - jl.push_b(jl_terms, jl.MOI.ScalarAffineTerm(float(coef), jl_var)) - - return jl.MOI.ScalarAffineFunction(jl_terms, float(constant)) - - -def _expr_to_moi(jl, all_jl_vars, expr): - """Convert a Python Expr to a concrete Julia MOI function (no iterators).""" - # Try linear first - linear = _expr_to_moi_linear(jl, all_jl_vars, expr) - if linear is not None: - return linear - - match expr: - case Constant(value=v): - return v - case Variable(index=idx): - return _get_jl_var_by_index(jl, all_jl_vars, idx) - case BinaryOp(op=op, left=left, right=right): - l = _expr_to_moi(jl, all_jl_vars, left) - r = _expr_to_moi(jl, all_jl_vars, right) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol(op), _any_vec(l, r)) - case UnaryOp(op="-", arg=arg): - a = _expr_to_moi(jl, all_jl_vars, arg) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol("-"), _any_vec(a)) - case Func(name=name, arg=arg): - a = _expr_to_moi(jl, all_jl_vars, arg) - return jl.MOI.ScalarNonlinearFunction(jl.Symbol(name), _any_vec(a)) - case _: - raise TypeError(f"Cannot convert {type(expr).__name__} to MOI") - - -def _add_individual_constraint(jl, optimizer, all_jl_vars, con): - """Add a single non-grouped constraint.""" - # Normalize: lhs - rhs in {set} - from jumpy.expressions import BinaryOp as _BinaryOp, Constant as _Constant - normalized = con.lhs - con.rhs - - func = _expr_to_moi(jl, all_jl_vars, normalized) - - if con.sense == "<=": - set_ = jl.MOI.LessThan(0.0) - elif con.sense == ">=": - set_ = jl.MOI.GreaterThan(0.0) - elif con.sense == "==": - set_ = jl.MOI.EqualTo(0.0) - else: - raise ValueError(f"Unknown constraint sense: {con.sense}") - jl.MOI.Utilities.normalize_and_add_constraint(optimizer, func, set_) + def add_constraint_group(self, func, sense, linear): + jl = self._jl + if linear: + target = "MOI.ScalarAffineFunction{Float64}" + else: + target = "MOI.ScalarNonlinearFunction" + template, iterators = jl.GenOpt.collect_iterator_refs(func) + generator = jl.seval(f"GenOpt.FunctionGenerator{{{target}}}")(template, iterators) + n = jl.MOI.output_dimension(generator) + if sense == "<=": + set_ = jl.MOI.Nonpositives(n) + elif sense == ">=": + set_ = jl.MOI.Nonnegatives(n) + elif sense == "==": + set_ = jl.MOI.Zeros(n) + else: + raise ValueError(f"Unknown constraint sense: {sense}") + jl.MOI.add_constraint(self._optimizer, generator, set_) + + def set_objective(self, sense, func): + jl = self._jl + moi_sense = jl.MOI.MIN_SENSE if sense == "min" else jl.MOI.MAX_SENSE + jl.MOI.set(self._optimizer, jl.MOI.ObjectiveSense(), moi_sense) + func = self._simplify(func) + jl.MOI.set(self._optimizer, self._objective_attr(func), func) + + def optimize(self): + jl = self._jl + jl.MOI.optimize_b(self._optimizer) + return int(jl.Integer(jl.MOI.get(self._optimizer, jl.MOI.TerminationStatus()))) + + def get_values(self, count): + jl = self._jl + return [ + float(jl.MOI.get(self._optimizer, jl.MOI.VariablePrimal(), v)) + for v in self._variables[:count] + ] diff --git a/src/jumpy/expressions.py b/src/jumpy/expressions.py index 3c59674..dd66b28 100644 --- a/src/jumpy/expressions.py +++ b/src/jumpy/expressions.py @@ -1,287 +1,220 @@ """ -Expression graph for JuMPy. - -Builds a tree of nodes that maps directly to MOI.ScalarNonlinearFunction: - ScalarNonlinearFunction(head::Symbol, args::Vector{Any}) - -Each node is either: - - Variable(index) -> MOI.VariableIndex(index) - - Constant(value) -> Float64 - - BinaryOp(op, l, r) -> MOI.ScalarNonlinearFunction(op, [l, r]) - - UnaryOp(op, arg) -> MOI.ScalarNonlinearFunction(op, [arg]) - - Func(name, arg) -> MOI.ScalarNonlinearFunction(name, [arg]) - - IteratorRef(iterator) -> GeneratorOptInterface.IteratorIndex - - IndexedVariable(vec, index) -> variable lookup resolved during expansion - - IndexedParameter(param, idx) -> data lookup resolved during expansion -""" - -from __future__ import annotations - -from typing import Union - -Numeric = Union[int, float] - - -def _wrap(other: Expr | Numeric) -> Expr: - """Wrap a numeric literal into a Constant node.""" - if isinstance(other, Expr): - return other - if isinstance(other, (int, float)): - return Constant(float(other)) - raise TypeError(f"Cannot convert {type(other).__name__} to Expr") - - -class Expr: - """Base class for all expression-graph nodes.""" - - # -- arithmetic operators -------------------------------------------------- - - def __add__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("+", self, _wrap(other)) - - def __radd__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("+", _wrap(other), self) - - def __sub__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("-", self, _wrap(other)) +Expression nodes, built eagerly as MOI functions. - def __rsub__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("-", _wrap(other), self) +Every arithmetic operation immediately performs one MOI call through the +model's `ops` object — via juliacall or the compiled library. A Node is a +thin Python handle around the resulting MOI object; there is no Python-side +expression tree and no conversion step. - def __mul__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("*", self, _wrap(other)) +Iterators (from Model.iterator) are ordinary nodes wrapping a +GenOpt.IteratorRef, so templates like `x[i] + x[i + 1] <= 10` are also +built eagerly; GenOpt discovers the iterators by identity when the group +constraint is added. +""" - def __rmul__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("*", _wrap(other), self) +from __future__ import annotations - def __truediv__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("/", self, _wrap(other)) +Numeric = (int, float) - def __rtruediv__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("/", _wrap(other), self) - def __pow__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("^", self, _wrap(other)) +def _moi(ops, value): + """The MOI object of a Node or a numeric literal.""" + if isinstance(value, Node): + return value.moi + if isinstance(value, Numeric): + return ops.constant(float(value)) + raise TypeError(f"Cannot use {type(value).__name__} in an expression") - def __rpow__(self, other: Expr | Numeric) -> BinaryOp: - return BinaryOp("^", _wrap(other), self) - def __neg__(self) -> UnaryOp: - return UnaryOp("-", self) +def _is_linear(value) -> bool: + return not isinstance(value, Node) or value.linear - def __pos__(self) -> Expr: - return self - # -- comparison operators (return Constraint objects) ---------------------- +class Node: + """A handle to an MOI expression owned by the model's backend.""" - def __le__(self, other: Expr | Numeric) -> Constraint: - return Constraint(self, "<=", _wrap(other)) + def __init__(self, ops, moi, *, linear: bool = True): + self._ops = ops + self.moi = moi + self.linear = linear - def __ge__(self, other: Expr | Numeric) -> Constraint: - return Constraint(self, ">=", _wrap(other)) + def _snf(self, head: str, args, *, linear: bool = True) -> Node: + return Node( + self._ops, + self._ops.scalar_nonlinear(head, [_moi(self._ops, a) for a in args]), + linear=linear and all(_is_linear(a) for a in args), + ) - def __eq__(self, other: Expr | Numeric) -> Constraint: - return Constraint(self, "==", _wrap(other)) + # -- arithmetic (each call is one MOI ScalarNonlinearFunction) -------------- + def __add__(self, other): + return self._snf("+", [self, other]) -class Variable(Expr): - """ - A decision variable. Maps to MOI.VariableIndex(index). + def __radd__(self, other): + return self._snf("+", [other, self]) - Users don't create these directly; they are returned by Model.variables(). - """ + def __sub__(self, other): + return self._snf("-", [self, other]) - def __init__(self, index: int, name: str | None = None): - self.index = index - self.name = name - - def __repr__(self) -> str: - if self.name: - return self.name - return f"x[{self.index}]" + def __rsub__(self, other): + return self._snf("-", [other, self]) + def __mul__(self, other): + return self._snf("*", [self, other]) -class Constant(Expr): - """A numeric constant in the expression tree.""" + def __rmul__(self, other): + return self._snf("*", [other, self]) - def __init__(self, value: float): - self.value = value + def __truediv__(self, other): + return self._snf("/", [self, other], linear=False) - def __repr__(self) -> str: - return str(self.value) + def __rtruediv__(self, other): + return self._snf("/", [other, self], linear=False) + def __pow__(self, other): + return self._snf("^", [self, other], linear=False) -class BinaryOp(Expr): - """Binary operation node: +, -, *, /, ^.""" + def __rpow__(self, other): + return self._snf("^", [other, self], linear=False) - def __init__(self, op: str, left: Expr, right: Expr): - self.op = op - self.left = left - self.right = right + def __neg__(self): + return self._snf("-", [self]) - def __repr__(self) -> str: - return f"({self.left} {self.op} {self.right})" + def __pos__(self): + return self + # -- comparisons (normalized to `self - other sense 0`) ------------------- -class UnaryOp(Expr): - """Unary operation node (currently just negation).""" + def __le__(self, other) -> Constraint: + return Constraint(self - other, "<=") - def __init__(self, op: str, arg: Expr): - self.op = op - self.arg = arg + def __ge__(self, other) -> Constraint: + return Constraint(self - other, ">=") - def __repr__(self) -> str: - return f"({self.op}{self.arg})" + def __eq__(self, other) -> Constraint: + return Constraint(self - other, "==") -class Func(Expr): - """Named function call: sin, cos, exp, log, sqrt, abs.""" +class Variable(Node): + """A single decision variable; keeps its column for solution lookup.""" - def __init__(self, name: str, arg: Expr): + def __init__(self, ops, index: int, name: str | None = None): + super().__init__(ops, ops.variable(index)) + self.index = index self.name = name - self.arg = _wrap(arg) def __repr__(self) -> str: - return f"{self.name}({self.arg})" - - -class IndexedVariable(Expr): - """ - Symbolic variable lookup: x[expr] where expr contains IteratorRefs. - - During expansion on the Julia side, the index expression is evaluated - for each iterator value to resolve to a concrete MOI.VariableIndex. - """ - - def __init__(self, variable_vector: VariableVector, index_expr: Expr): - self.variable_vector = variable_vector - self.index_expr = index_expr - - def __repr__(self) -> str: - name = self.variable_vector.name or "x" - return f"{name}[{self.index_expr}]" - - -class IndexedParameter(Expr): - """ - Symbolic data lookup: param[expr] where expr contains IteratorRefs. - - During expansion, the index expression is evaluated for each iterator - value to look up a concrete float from the parameter data. - """ - - def __init__(self, parameter: Parameter, index_expr: Expr): - self.parameter = parameter - self.index_expr = index_expr - - def __repr__(self) -> str: - name = self.parameter.name or "p" - return f"{name}[{self.index_expr}]" + return self.name or f"x[{self.index}]" class VariableVector: """ A block of decision variables returned by Model.variables(). - Supports both concrete indexing (x[0] -> Variable) and symbolic - indexing (x[i] -> IndexedVariable, where i is an Iterator/Expr). + Concrete indexing (x[0]) returns a Variable; symbolic indexing (x[i] + with an expression) builds a getindex template node over the contiguous + block. """ - def __init__(self, variables: list[Variable], name: str | None = None): - self._variables = variables + def __init__(self, ops, start: int, count: int, name: str | None = None): + self._ops = ops + self.start = start + self.count = count self.name = name - - def __getitem__(self, index: int | Expr) -> Variable | IndexedVariable: - if isinstance(index, (int,)): - return self._variables[index] - if isinstance(index, Expr): - return IndexedVariable(self, index) - raise TypeError(f"Index must be int or Expr, got {type(index).__name__}") + self._block = None # GenOpt.ContiguousArrayOfVariables, built lazily + + def __getitem__(self, index): + if isinstance(index, int): + var_name = f"{self.name}[{index}]" if self.name else None + return Variable(self._ops, self.start + index, var_name) + if isinstance(index, Node): + if self._block is None: + self._block = self._ops.contiguous_variables(self.start, self.count) + block = Node(self._ops, self._block) + # 0-based Python index -> 1-based Julia index + return block._snf("getindex", [block, index + 1]) + raise TypeError(f"Index must be int or Node, got {type(index).__name__}") def __len__(self) -> int: - return len(self._variables) + return self.count def __iter__(self): - return iter(self._variables) + return (self[k] for k in range(self.count)) def __repr__(self) -> str: - name = self.name or "x" - return f"{name}[0:{len(self._variables)}]" + return f"{self.name or 'x'}[0:{self.count}]" class Parameter: """ - A vector of constant data for use in constraint group templates. - - Supports symbolic indexing: costs[i] where i is an Iterator. + A vector of constant data returned by Model.parameter(). - Example: - costs = jp.Parameter([3.0, 1.5, 2.0]) - m.constraint_group([i], costs[i] * x[i] <= 50) + Concrete indexing (costs[0]) returns a float; symbolic indexing + (costs[i]) builds a getindex template node over the data vector. """ - def __init__(self, values: list[float], name: str | None = None): + def __init__(self, ops, values, name: str | None = None): + self._ops = ops self.values = [float(v) for v in values] self.name = name - - def __getitem__(self, index: int | Expr) -> Constant | IndexedParameter: - if isinstance(index, (int,)): - return Constant(self.values[index]) - if isinstance(index, Expr): - return IndexedParameter(self, index) - raise TypeError(f"Index must be int or Expr, got {type(index).__name__}") + self._array = None # data vector node, built lazily + + def __getitem__(self, index): + if isinstance(index, int): + return self.values[index] + if isinstance(index, Node): + if self._array is None: + self._array = self._ops.float_array(self.values) + array = Node(self._ops, self._array) + return array._snf("getindex", [array, index + 1]) + raise TypeError(f"Index must be int or Node, got {type(index).__name__}") def __len__(self) -> int: return len(self.values) def __repr__(self) -> str: - name = self.name or "param" - return f"{name}[0:{len(self.values)}]" + return f"{self.name or 'param'}[0:{len(self.values)}]" class Constraint: - """ - Represents lhs {<=, >=, ==} rhs. - - Normalized before passing to MOI as: (lhs - rhs) in {Nonpositives, Nonnegatives, Zeros}. - """ + """A normalized constraint: `func sense 0`.""" - def __init__(self, lhs: Expr, sense: str, rhs: Expr): - self.lhs = lhs + def __init__(self, func: Node, sense: str): + self.func = func self.sense = sense - self.rhs = rhs def __repr__(self) -> str: - return f"{self.lhs} {self.sense} {self.rhs}" + return f"" class Objective: """An optimization objective (minimize or maximize).""" - def __init__(self, sense: str, expr: Expr): + def __init__(self, sense: str, func: Node): self.sense = sense - self.expr = expr + self.func = func - def __repr__(self) -> str: - return f"{self.sense}({self.expr})" +# -- nonlinear functions -------------------------------------------------------- + +def _func(name: str, x: Node) -> Node: + return x._snf(name, [x], linear=False) -# -- convenience functions that return Func nodes ----------------------------- -def sin(x: Expr | Numeric) -> Func: - return Func("sin", _wrap(x)) +def sin(x): + return _func("sin", x) -def cos(x: Expr | Numeric) -> Func: - return Func("cos", _wrap(x)) +def cos(x): + return _func("cos", x) -def exp(x: Expr | Numeric) -> Func: - return Func("exp", _wrap(x)) +def exp(x): + return _func("exp", x) -def log(x: Expr | Numeric) -> Func: - return Func("log", _wrap(x)) +def log(x): + return _func("log", x) -def sqrt(x: Expr | Numeric) -> Func: - return Func("sqrt", _wrap(x)) +def sqrt(x): + return _func("sqrt", x) -def abs(x: Expr | Numeric) -> Func: - return Func("abs", _wrap(x)) +def abs(x): + return _func("abs", x) diff --git a/src/jumpy/iterators.py b/src/jumpy/iterators.py deleted file mode 100644 index eda71ec..0000000 --- a/src/jumpy/iterators.py +++ /dev/null @@ -1,41 +0,0 @@ -""" -Iterators for constraint groups. - -An Iterator IS an Expr node. When used in arithmetic (10*i + j), it builds -an expression graph. When used as an index (x[i], costs[i]), it creates -symbolic IndexedVariable / IndexedParameter nodes. - -Maps to GeneratorOptInterface.Iterator(length, values) on the Julia side. -""" - -from __future__ import annotations - -from jumpy.expressions import Expr - - -class Iterator(Expr): - """ - An index set for constraint groups. - - When used in expressions, it acts as a symbolic placeholder that - GeneratorOptInterface expands over its values during constraint generation. - - Example: - i = jp.Iterator(range(99)) - j = jp.Iterator(range(10)) - - x[i] # symbolic variable lookup - x[10*i + j] # symbolic index arithmetic - costs[i] * x[i] # symbolic data + variable lookup - """ - - _next_id: int = 0 - - def __init__(self, values): - self.values = list(values) - self.length = len(self.values) - self.id = Iterator._next_id - Iterator._next_id += 1 - - def __repr__(self) -> str: - return f"i{self.id}" diff --git a/src/jumpy/model.py b/src/jumpy/model.py index 742ecda..d1a8e0e 100644 --- a/src/jumpy/model.py +++ b/src/jumpy/model.py @@ -1,91 +1,33 @@ """ The Model class: top-level API for building optimization models in JuMPy. -A Model collects variables, constraint groups, and an objective, then hands -everything off to the compiled Julia library for expansion and solving. +The model is built eagerly: every call performs the corresponding MOI call +through the backend's ops object (juliacall or the compiled library). +optimize() is just MOI.optimize! plus solution retrieval. """ from __future__ import annotations -from dataclasses import dataclass -from typing import Callable - +from jumpy.backend import get_ops from jumpy.expressions import ( Constraint, - Expr, + Node, Objective, Parameter, Variable, VariableVector, ) -from jumpy.iterators import Iterator -from jumpy.serialize import serialize_constraint, serialize_expr -from jumpy.backend import Backend, get_backend - - -def minimize(expr: Expr) -> Objective: - return Objective("min", expr) - - -def maximize(expr: Expr) -> Objective: - return Objective("max", expr) - - -def sum_over(iterator: Iterator, expr: Expr) -> Expr: - """ - Symbolic sum over an iterator. - - Passed to Julia as a SumNode that GenOpt expands, NOT evaluated in Python. - - Example: - jp.sum_over(i, costs[i] * x[i]) - """ - return SumExpr(iterator, expr) - -@dataclass -class ConstraintGroup: - """ - A group of constraints defined by a template expression + iterators. - - Maps to GeneratorOptInterface.IteratedFunction: - IteratedFunction(func::MOI.ScalarNonlinearFunction, iterators::Vector{Iterator}) - - The template is a Constraint containing Iterator nodes as placeholders. - Expansion happens entirely on the Julia side. - """ - - template: Constraint - iterators: list[Iterator] - - def __repr__(self) -> str: - n = 1 - for it in self.iterators: - n *= it.length - return f"ConstraintGroup({n} constraints from {len(self.iterators)} iterator(s))" - - -@dataclass -class VariableBlock: - """A contiguous block of variables with shared bounds.""" - start: int - count: int - lower: float | None - upper: float | None - vector: VariableVector - binary: bool = False - integer: bool = False +# MOI.OPTIMAL in MOI.TerminationStatusCode. +OPTIMAL = 1 -class SumExpr(Expr): - """Symbolic sum over an iterator. Expanded on the Julia side.""" +def minimize(func: Node) -> Objective: + return Objective("min", func) - def __init__(self, iterator: Iterator, body: Expr): - self.iterator = iterator - self.body = body - def __repr__(self) -> str: - return f"sum({self.iterator}, {self.body})" +def maximize(func: Node) -> Objective: + return Objective("max", func) class Model: @@ -96,10 +38,10 @@ class Model: m = jp.Model() x = m.variables(100, lower=0) - i = jp.Iterator(range(99)) - m.constraint_group([i], x[i] + x[i + 1] <= 10) + i = m.iterator(range(99)) + m.constraint_group(x[i] + x[i + 1] <= 10) - m.objective = jp.minimize(sum(x)) + m.objective = jp.minimize(x[0] + x[1]) m.optimize() """ @@ -111,14 +53,22 @@ def __init__(self, backend: str = "juliac"): backend: "juliac" (default, no Julia needed) or "juliacall" (uses juliacall, installs Julia lazily if needed). """ - self._backend: Backend = get_backend(backend) - self._var_blocks: list[VariableBlock] = [] - self._constraint_groups: list[ConstraintGroup] = [] - self._individual_constraints: list[Constraint] = [] + self._ops = get_ops(backend) + self._num_vars = 0 self._objective: Objective | None = None - self._num_vars: int = 0 self._solution: list[float] | None = None + def close(self) -> None: + """Release the backend model. The model must not be used afterwards.""" + # getattr: __del__ may run when __init__ failed before setting _ops + ops = getattr(self, "_ops", None) + if ops is not None: + ops.free() + self._ops = None + + def __del__(self): + self.close() + # -- Variables ------------------------------------------------------------- def variables( @@ -131,23 +81,24 @@ def variables( binary: bool = False, integer: bool = False, ) -> VariableVector: - """ - Add a block of decision variables. - - Returns a VariableVector that supports both concrete (x[0]) - and symbolic (x[i]) indexing. - """ - start = self._num_vars - vars = [] + """Add a block of decision variables (MOI.add_variables + bounds).""" + start = self._ops.add_variables(count) + # Bounds and integrality are VariableIndex-in-set constraints, as in MOI. for k in range(count): - var_name = f"{name}[{k}]" if name else None - vars.append(Variable(start + k, var_name)) + if lower is not None: + self._ops.add_constraint( + self._ops.variable(start + k), ">=", float(lower), + ) + if upper is not None: + self._ops.add_constraint( + self._ops.variable(start + k), "<=", float(upper), + ) + if binary: + self._ops.add_constraint(self._ops.variable(start + k), "binary", 0.0) + elif integer: + self._ops.add_constraint(self._ops.variable(start + k), "integer", 0.0) self._num_vars += count - vec = VariableVector(vars, name) - self._var_blocks.append( - VariableBlock(start, count, lower, upper, vec, binary, integer) - ) - return vec + return VariableVector(self._ops, start, count, name) def variable( self, @@ -159,41 +110,41 @@ def variable( integer: bool = False, ) -> Variable: """Add a single decision variable.""" - vec = self.variables( + return self.variables( 1, lower=lower, upper=upper, name=name, binary=binary, integer=integer, - ) - return vec[0] + )[0] - # -- Constraints ----------------------------------------------------------- + # -- Template data ----------------------------------------------------------- - def constraint_group( - self, - iterators: list[Iterator], - template: Constraint, - ) -> ConstraintGroup: + def iterator(self, values) -> Node: """ - Add a constraint group. + An index set for constraint groups (a GenOpt iterator). - The template is an expression containing Iterator nodes as symbolic - placeholders. GeneratorOptInterface expands the template over all - iterator values entirely in compiled Julia. + Used in expressions, it is a symbolic placeholder that GenOpt + expands over its values when the group constraint is added. + """ + return Node(self._ops, self._ops.iterator([float(v) for v in values])) - Example: - i = jp.Iterator(range(99)) - m.constraint_group([i], x[i] + x[i + 1] <= 10) + def parameter(self, values, name: str | None = None) -> Parameter: + """A vector of constant data, symbolically indexable in templates.""" + return Parameter(self._ops, values, name) + + # -- Constraints ----------------------------------------------------------- - i = jp.Iterator(range(10)) - j = jp.Iterator(range(10)) - m.constraint_group([i, j], x[10*i + j] >= 0) + def constraint(self, con: Constraint) -> None: + """Add a single constraint (MOI.add_constraint).""" + self._ops.add_constraint(con.func.moi, con.sense, 0.0) + + def constraint_group(self, con: Constraint) -> None: """ - group = ConstraintGroup(template, iterators) - self._constraint_groups.append(group) - return group + Add a constraint group: one constraint per combination of the + values of the iterators appearing in the template. - def constraint(self, con: Constraint) -> Constraint: - """Add a single (non-grouped) constraint.""" - self._individual_constraints.append(con) - return con + Example: + i = m.iterator(range(99)) + m.constraint_group(x[i] + x[i + 1] <= 10) + """ + self._ops.add_constraint_group(con.func.moi, con.sense, con.func.linear) # -- Objective ------------------------------------------------------------- @@ -204,73 +155,21 @@ def objective(self) -> Objective | None: @objective.setter def objective(self, obj: Objective) -> None: self._objective = obj + self._ops.set_objective(obj.sense, obj.func.moi) # -- Solve ----------------------------------------------------------------- def optimize(self) -> None: - """ - Solve the model using the selected backend. - - - juliac backend: serializes to flat arrays, calls compiled shared library - - juliacall backend: builds MOI model directly in Julia via juliacall - """ - self._solution = self._backend.optimize(self) - - def _serialize(self) -> dict: - """Serialize the entire model for the Julia C ABI.""" - param_registry: dict = {} - data: dict = { - "num_vars": self._num_vars, - "var_blocks": [], - "constraint_groups": [], - "individual_constraints": [], - "objective": None, - "parameters": [], # filled at the end from param_registry - } - - for block in self._var_blocks: - data["var_blocks"].append({ - "start": block.start, - "count": block.count, - "lower": block.lower, - "upper": block.upper, - }) - - for group in self._constraint_groups: - serialized_con = serialize_constraint(group.template, param_registry) - iterators = [] - for it in group.iterators: - iterators.append({ - "id": it.id, - "length": it.length, - "values": [float(v) for v in it.values], - }) - data["constraint_groups"].append({ - **serialized_con, - "iterators": iterators, - }) - - for con in self._individual_constraints: - data["individual_constraints"].append( - serialize_constraint(con, param_registry) + """MOI.optimize!, then retrieve the solution.""" + status = self._ops.optimize() + if status != OPTIMAL: + raise RuntimeError( + f"Solve did not reach OPTIMAL (termination status {status})" ) - - if self._objective: - data["objective"] = { - "sense": self._objective.sense, - "expr": serialize_expr(self._objective.expr, param_registry), - } - - # Collect all referenced parameters - params_by_id = sorted(param_registry.values(), key=lambda p: p["id"]) - data["parameters"] = [p["values"] for p in params_by_id] - - return data - - # -- Solution retrieval ---------------------------------------------------- + self._solution = self._ops.get_values(self._num_vars) def value(self, var: Variable) -> float: """Get the solved value of a variable.""" if self._solution is None: raise RuntimeError("Model has not been solved yet. Call optimize() first.") - return self._solution[var.index] \ No newline at end of file + return self._solution[var.index] diff --git a/src/jumpy/serialize.py b/src/jumpy/serialize.py deleted file mode 100644 index 4833cf4..0000000 --- a/src/jumpy/serialize.py +++ /dev/null @@ -1,116 +0,0 @@ -""" -Serialization of Python expression graphs to a flat format for the Julia C ABI. - -The compiled Julia library (MOI + GeneratorOptInterface + Bridges + HiGHS) -receives models as flat arrays through ctypes. This module converts the Python -expression tree into that flat representation. - -Tags: - 0 = constant -> [0, float64_value] - 1 = variable -> [1, var_index] - 2 = binary op -> [2, op_code, left..., right...] - 3 = unary op -> [3, op_code, arg...] - 4 = function -> [4, func_code, arg...] - 5 = iterator -> [5, iterator_id] - 6 = indexed_var -> [6, var_block_start, var_block_count, index_expr...] - 7 = indexed_param -> [7, param_id, index_expr...] - -Op codes for binary: + = 0, - = 1, * = 2, / = 3, ^ = 4 -Op codes for unary: - = 0 -Func codes: sin = 0, cos = 1, exp = 2, log = 3, sqrt = 4, abs = 5 -""" - -from __future__ import annotations - -from jumpy.expressions import ( - BinaryOp, - Constant, - Constraint, - Expr, - Func, - IndexedParameter, - IndexedVariable, - UnaryOp, - Variable, -) -from jumpy.iterators import Iterator - -# Tag constants -TAG_CONST = 0 -TAG_VAR = 1 -TAG_BINARY = 2 -TAG_UNARY = 3 -TAG_FUNC = 4 -TAG_ITERATOR = 5 -TAG_INDEXED_VAR = 6 -TAG_INDEXED_PARAM = 7 - -BINARY_OPS = {"+": 0, "-": 1, "*": 2, "/": 3, "^": 4} -UNARY_OPS = {"-": 0} -FUNC_CODES = {"sin": 0, "cos": 1, "exp": 2, "log": 3, "sqrt": 4, "abs": 5} - - -def serialize_expr(expr: Expr, param_registry: dict | None = None) -> list[float]: - """ - Serialize an expression tree to a flat list of floats. - - This is the format passed across the ctypes boundary to the compiled - Julia library, which reconstructs it into MOI.ScalarNonlinearFunction. - """ - if param_registry is None: - param_registry = {} - buf: list[float] = [] - _serialize(expr, buf, param_registry) - return buf - - -def _serialize(node: Expr, buf: list[float], param_registry: dict) -> None: - match node: - case Constant(value=v): - buf.extend([TAG_CONST, v]) - case Variable(index=idx): - buf.extend([TAG_VAR, float(idx)]) - case BinaryOp(op=op, left=left, right=right): - buf.extend([TAG_BINARY, float(BINARY_OPS[op])]) - _serialize(left, buf, param_registry) - _serialize(right, buf, param_registry) - case UnaryOp(op=op, arg=arg): - buf.extend([TAG_UNARY, float(UNARY_OPS[op])]) - _serialize(arg, buf, param_registry) - case Func(name=name, arg=arg): - buf.extend([TAG_FUNC, float(FUNC_CODES[name])]) - _serialize(arg, buf, param_registry) - case Iterator() as it: - buf.extend([TAG_ITERATOR, float(it.id)]) - case IndexedVariable() as iv: - start = iv.variable_vector._variables[0].index - count = len(iv.variable_vector) - buf.extend([TAG_INDEXED_VAR, float(start), float(count)]) - _serialize(iv.index_expr, buf, param_registry) - case IndexedParameter() as ip: - param_id = id(ip.parameter) - if param_id not in param_registry: - param_registry[param_id] = { - "id": len(param_registry), - "values": ip.parameter.values, - } - buf.extend([TAG_INDEXED_PARAM, float(param_registry[param_id]["id"])]) - _serialize(ip.index_expr, buf, param_registry) - case _: - raise TypeError(f"Cannot serialize {type(node).__name__}") - - -def serialize_constraint(con: Constraint, param_registry: dict | None = None) -> dict: - """ - Serialize a constraint: normalized expression + sense. - - Returns: - {"expr": [...flat...], "sense": "<="|">="|"=="} - """ - if param_registry is None: - param_registry = {} - normalized = con.lhs - con.rhs # f(x) {<=,>=,==} 0 - return { - "expr": serialize_expr(normalized, param_registry), - "sense": con.sense, - } diff --git a/tests/test_expressions.py b/tests/test_expressions.py index 65e255d..1d5c87d 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -1,271 +1,141 @@ -"""Tests for the JuMPy expression graph, symbolic indexing, and serialization.""" +""" +Tests for eager expression building, using a mock ops object. + +The mock records every MOI call as a tuple, so these tests check exactly +what a backend receives — pure Python, no Julia needed. +""" import sys sys.path.insert(0, "src") -from jumpy import ( - Model, Iterator, Parameter, Variable, Constant, - sin, exp, minimize, maximize, sum_over, -) -from jumpy.expressions import BinaryOp, Constraint, IndexedVariable, IndexedParameter -from jumpy.serialize import ( - serialize_expr, - serialize_constraint, - TAG_CONST, - TAG_VAR, - TAG_BINARY, - TAG_FUNC, - TAG_ITERATOR, - TAG_INDEXED_VAR, - TAG_INDEXED_PARAM, -) - - -# -- Expression building ------------------------------------------------------- - -def test_basic_arithmetic(): - x = Variable(0, "x") - y = Variable(1, "y") - expr = x + 2 * y - assert isinstance(expr, BinaryOp) - assert repr(expr) == "(x + (2.0 * y))" - - -def test_constraint(): - x = Variable(0, "x") - y = Variable(1, "y") - con = x + y <= 10 - assert isinstance(con, Constraint) - assert con.sense == "<=" - - -def test_nonlinear(): - x = Variable(0, "x") - expr = sin(x) + exp(x) - assert repr(expr) == "(sin(x) + exp(x))" +from jumpy.expressions import Node, Parameter, Variable, VariableVector, sin -def test_negation(): - x = Variable(0, "x") - expr = -x - assert repr(expr) == "(-x)" +class MockOps: + def __init__(self): + self.constraints = [] + self.groups = [] + self.num_vars = 0 + def constant(self, v): + return v -# -- Iterator as Expr (the core idea) ------------------------------------------ + def variable(self, index): + return ("var", index) -def test_iterator_in_arithmetic(): - """Iterator objects participate in arithmetic to build index expressions.""" - i = Iterator(range(10)) - j = Iterator(range(10)) + def scalar_nonlinear(self, head, args): + return (head, *args) - expr = 10 * i + j - assert isinstance(expr, BinaryOp) - assert expr.op == "+" - assert isinstance(expr.left, BinaryOp) - assert expr.left.op == "*" + def iterator(self, values): + return ("iterator", tuple(values)) + def contiguous_variables(self, start, count): + return ("block", start, count) -def test_symbolic_variable_indexing(): - """x[i] with an Iterator produces an IndexedVariable.""" - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(99)) + def float_array(self, values): + return ("data", tuple(values)) - indexed = x[i] - assert isinstance(indexed, IndexedVariable) + def add_variables(self, count): + start = self.num_vars + self.num_vars += count + return start + def add_constraint(self, func, sense, rhs): + self.constraints.append((func, sense, rhs)) -def test_symbolic_variable_index_arithmetic(): - """x[i + 1] builds an index expression graph.""" - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(99)) + def add_constraint_group(self, func, sense, linear): + self.groups.append((func, sense, linear)) - indexed = x[i + 1] - assert isinstance(indexed, IndexedVariable) - assert isinstance(indexed.index_expr, BinaryOp) - -def test_multidim_index(): - """x[10*i + j] builds a compound index expression.""" - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(10)) - j = Iterator(range(10)) - - indexed = x[10 * i + j] - assert isinstance(indexed, IndexedVariable) - - -def test_parameter_symbolic_indexing(): - """costs[i] produces an IndexedParameter.""" - costs = Parameter([3.0, 1.5, 2.0], name="costs") - i = Iterator(range(3)) - - indexed = costs[i] - assert isinstance(indexed, IndexedParameter) +def test_arithmetic_builds_scalar_nonlinear(): + ops = MockOps() + x = Variable(ops, 0) + y = Variable(ops, 1) + expr = x + 2 * y + assert expr.moi == ("+", ("var", 0), ("*", 2.0, ("var", 1))) + assert expr.linear -def test_parameter_concrete_indexing(): - """costs[0] returns a Constant.""" - costs = Parameter([3.0, 1.5, 2.0]) - c = costs[0] - assert isinstance(c, Constant) - assert c.value == 3.0 +def test_reflected_and_unary_operators(): + ops = MockOps() + x = Variable(ops, 0) + assert (1 - x).moi == ("-", 1.0, ("var", 0)) + assert (-x).moi == ("-", ("var", 0)) + assert (2.0 / x).linear is False + assert (x**2).linear is False -# -- Full constraint group expressions ----------------------------------------- +def test_nonlinear_functions(): + ops = MockOps() + x = Variable(ops, 0) + expr = sin(x) + 1 + assert expr.moi == ("+", ("sin", ("var", 0)), 1.0) + assert not expr.linear -def test_constraint_group_expression(): - """The full expression x[i] + x[i+1] <= 10 builds a valid tree.""" - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(99)) - con = x[i] + x[i + 1] <= 10 - assert isinstance(con, Constraint) +def test_comparison_normalizes(): + ops = MockOps() + x = Variable(ops, 0) + con = x + 1 <= 10 assert con.sense == "<=" - assert isinstance(con.lhs, BinaryOp) # x[i] + x[i+1] - assert isinstance(con.lhs.left, IndexedVariable) - assert isinstance(con.lhs.right, IndexedVariable) - - -def test_nonlinear_constraint_group(): - """sin(x[i]) + exp(x[i]) <= 1.0 with symbolic indexing.""" - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(100)) - - con = sin(x[i]) + exp(x[i]) <= 1.0 - assert isinstance(con, Constraint) - - -def test_parameter_in_constraint_group(): - """costs[i] * x[i] <= 50 mixes data and variables.""" - m = Model() - x = m.variables(100, lower=0, name="x") - costs = Parameter([float(k) for k in range(100)], name="costs") - i = Iterator(range(100)) - - con = costs[i] * x[i] <= 50 - assert isinstance(con, Constraint) - assert isinstance(con.lhs, BinaryOp) - assert con.lhs.op == "*" - assert isinstance(con.lhs.left, IndexedParameter) - assert isinstance(con.lhs.right, IndexedVariable) - - -# -- Model API ----------------------------------------------------------------- - -def test_model_constraint_group(): - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(99)) - - group = m.constraint_group([i], x[i] + x[i + 1] <= 10) - assert len(m._constraint_groups) == 1 - assert "99 constraints" in repr(group) - - -def test_model_multidim_constraint_group(): - m = Model() - x = m.variables(100, lower=0, name="x") - i = Iterator(range(10)) - j = Iterator(range(10)) - - group = m.constraint_group([i, j], x[10 * i + j] >= 0) - assert len(m._constraint_groups) == 1 - assert "100 constraints" in repr(group) - - -def test_model_objective(): - m = Model() - x = m.variables(3, name="x") - m.objective = minimize(x[0] + x[1] + x[2]) - assert m.objective.sense == "min" - - -# -- Serialization ------------------------------------------------------------- - -def test_serialize_constant(): - assert serialize_expr(Constant(3.14)) == [TAG_CONST, 3.14] - - -def test_serialize_variable(): - assert serialize_expr(Variable(42)) == [TAG_VAR, 42.0] - - -def test_serialize_binary(): - x = Variable(0) - y = Variable(1) - buf = serialize_expr(x + y) - assert buf == [TAG_BINARY, 0.0, TAG_VAR, 0.0, TAG_VAR, 1.0] - - -def test_serialize_iterator(): - i = Iterator(range(10)) - buf = serialize_expr(i) - assert buf == [TAG_ITERATOR, float(i.id)] - - -def test_serialize_indexed_variable(): - m = Model() - x = m.variables(100, lower=0) - i = Iterator(range(99)) - - buf = serialize_expr(x[i]) - assert buf[0] == TAG_INDEXED_VAR - assert buf[1] == 0.0 # start index - assert buf[2] == 100.0 # count - assert TAG_ITERATOR in buf - - -def test_serialize_indexed_parameter(): - costs = Parameter([1.0, 2.0, 3.0]) - i = Iterator(range(3)) - - reg = {} - buf = serialize_expr(costs[i], reg) - assert buf[0] == TAG_INDEXED_PARAM - assert len(reg) == 1 - - -def test_serialize_full_model(): - """Serialize a complete model and verify the structure.""" - m = Model() - x = m.variables(100, lower=0, name="x") - - i = Iterator(range(99)) - m.constraint_group([i], x[i] + x[i + 1] <= 10) - - costs = Parameter([float(k) for k in range(100)], name="costs") - j = Iterator(range(100)) - m.constraint_group([j], costs[j] * x[j] <= 50) - - m.objective = minimize(x[0] + x[1]) - - data = m._serialize() - assert data["num_vars"] == 100 - assert len(data["var_blocks"]) == 1 - assert len(data["constraint_groups"]) == 2 - assert data["constraint_groups"][0]["iterators"][0]["length"] == 99 - assert data["constraint_groups"][1]["iterators"][0]["length"] == 100 - assert len(data["parameters"]) == 1 # costs - assert data["objective"]["sense"] == "min" + assert con.func.moi == ("-", ("+", ("var", 0), 1.0), 10.0) + + +def test_variable_vector_concrete_indexing(): + ops = MockOps() + start = ops.add_variables(3) + x = VariableVector(ops, start, 3, "x") + assert x[2].index == 2 + assert x[2].moi == ("var", 2) + assert len(x) == 3 + assert [v.index for v in x] == [0, 1, 2] + + +def test_variable_vector_symbolic_indexing(): + ops = MockOps() + x = VariableVector(ops, 0, 10, "x") + i = Node(ops, ops.iterator([0.0, 1.0])) + # 0-based Python index -> 1-based Julia index + assert x[i].moi == ( + "getindex", + ("block", 0, 10), + ("+", ("iterator", (0.0, 1.0)), 1.0), + ) + assert x[i].linear + + +def test_parameter_indexing(): + ops = MockOps() + p = Parameter(ops, [1.0, 2.0, 3.0], "costs") + assert p[1] == 2.0 + i = Node(ops, ops.iterator([0.0, 1.0])) + assert p[i].moi == ( + "getindex", + ("data", (1.0, 2.0, 3.0)), + ("+", ("iterator", (0.0, 1.0)), 1.0), + ) + + +def test_template_linearity_flag(): + ops = MockOps() + x = VariableVector(ops, 0, 10, "x") + i = Node(ops, ops.iterator([0.0, 1.0])) + assert (x[i] + x[i + 1] <= 10).func.linear + assert not (sin(x[i]) <= 1).func.linear if __name__ == "__main__": import traceback - tests = [v for k, v in globals().items() if k.startswith("test_")] + tests = [v for k, v in sorted(globals().items()) if k.startswith("test_")] passed = failed = 0 for test in tests: try: test() passed += 1 + print(f" PASS {test.__name__}") except Exception as e: failed += 1 - print(f"FAIL {test.__name__}: {e}") + print(f" FAIL {test.__name__}: {e}") traceback.print_exc() print(f"\n{passed} passed, {failed} failed") - import sys sys.exit(1 if failed else 0) diff --git a/tests/test_solve.py b/tests/test_solve.py index f05b11b..6eac902 100644 --- a/tests/test_solve.py +++ b/tests/test_solve.py @@ -1,17 +1,24 @@ """ -End-to-end tests that solve models with HiGHS via the juliacall backend. +End-to-end tests that solve models with HiGHS. -These require Julia + juliacall to be installed: - pip install juliacall +Backend is selected with the JUMPY_BACKEND environment variable: + JUMPY_BACKEND=juliac python tests/test_solve.py (default) + requires the compiled library (julia/README.md), no Julia needed + JUMPY_BACKEND=juliacall python tests/test_solve.py + requires Julia + juliacall: pip install juliacall """ +import os import sys sys.path.insert(0, "src") -from jumpy import Model, Iterator, Parameter, minimize, maximize +from jumpy import Model, minimize, maximize + +BACKEND = os.environ.get("JUMPY_BACKEND", "juliac") + def _model(): - return Model(backend="juliacall") + return Model(backend=BACKEND) def test_simple_lp(): @@ -35,8 +42,8 @@ def test_constraint_group_lp(): m = _model() x = m.variables(10, lower=0, name="x") - i = Iterator(range(10)) - m.constraint_group([i], x[i] >= 1) + i = m.iterator(range(10)) + m.constraint_group(x[i] >= 1) m.objective = minimize(sum(x)) m.optimize() @@ -47,13 +54,13 @@ def test_constraint_group_lp(): def test_constraint_group_consecutive(): """ - min x[0] s.t. x[i] + x[i+1] >= 2 for i in 0..8, x[i] >= 0 + min x[0]+x[1]+x[2] s.t. x[i] + x[i+1] >= 2 for i in 0..8, x[i] >= 0 """ m = _model() x = m.variables(10, lower=0, name="x") - i = Iterator(range(9)) - m.constraint_group([i], x[i] + x[i + 1] >= 2) + i = m.iterator(range(9)) + m.constraint_group(x[i] + x[i + 1] >= 2) m.objective = minimize(x[0] + x[1] + x[2]) m.optimize() @@ -70,10 +77,10 @@ def test_parameter_in_constraint_group(): """ m = _model() x = m.variables(5, lower=0, name="x") - demand = Parameter([1.0, 2.0, 3.0, 4.0, 5.0], name="demand") + demand = m.parameter([1.0, 2.0, 3.0, 4.0, 5.0], name="demand") - i = Iterator(range(5)) - m.constraint_group([i], x[i] >= demand[i]) + i = m.iterator(range(5)) + m.constraint_group(x[i] >= demand[i]) m.objective = minimize(sum(x)) m.optimize() @@ -90,9 +97,9 @@ def test_multidim_constraint_group(): m = _model() x = m.variables(9, lower=0, name="x") - i = Iterator(range(3)) - j = Iterator(range(3)) - m.constraint_group([i, j], x[3 * i + j] >= 1) + i = m.iterator(range(3)) + j = m.iterator(range(3)) + m.constraint_group(x[3 * i + j] >= 1) m.objective = minimize(sum(x)) m.optimize() @@ -101,6 +108,23 @@ def test_multidim_constraint_group(): assert abs(total - 9.0) < 1e-6 +def test_constraint_group_over_bounded_variables(): + """ + x[i] >= 0 as a group on variables that already have lower=0: must be + rows (like JuMP's @constraint), not clashing variable bounds. + """ + m = _model() + x = m.variables(4, lower=0, name="x") + + i = m.iterator(range(4)) + m.constraint_group(x[i] >= 0) + + m.objective = minimize(sum(x)) + m.optimize() + + assert abs(sum(m.value(v) for v in x)) < 1e-6 + + def test_maximize(): """max x s.t. x <= 42, x >= 0""" m = _model()