diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..ad37434 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,59 @@ +name: CI + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + test-expressions: + name: Python expression tests + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13"] + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Run expression tests + run: python tests/test_expressions.py + + test-solve: + name: End-to-end solve (juliacall + HiGHS) + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.11", "3.12"] + julia-version: ["1.11", "1.12"] + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.julia-version }} + + - 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 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..1b5ea51 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +__pycache__/ +*.pyc +*.egg-info/ +dist/ +build/ diff --git a/pyproject.toml b/pyproject.toml index ac2f484..0538ad6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,9 +8,10 @@ version = "0.1.0" description = "A Python interface to JuMP's MathOptInterface via GeneratorOptInterface" requires-python = ">=3.10" license = "MIT" -dependencies = [ - "numpy", -] +dependencies = [] + +[project.optional-dependencies] +juliacall = ["juliacall>=0.9.23"] [tool.hatch.build.targets.wheel] packages = ["src/jumpy"] diff --git a/src/jumpy/__pycache__/__init__.cpython-311.pyc b/src/jumpy/__pycache__/__init__.cpython-311.pyc deleted file mode 100644 index 6dd6b2d..0000000 Binary files a/src/jumpy/__pycache__/__init__.cpython-311.pyc and /dev/null differ diff --git a/src/jumpy/__pycache__/expressions.cpython-311.pyc b/src/jumpy/__pycache__/expressions.cpython-311.pyc deleted file mode 100644 index 56290c6..0000000 Binary files a/src/jumpy/__pycache__/expressions.cpython-311.pyc and /dev/null differ diff --git a/src/jumpy/__pycache__/iterators.cpython-311.pyc b/src/jumpy/__pycache__/iterators.cpython-311.pyc deleted file mode 100644 index ecaf839..0000000 Binary files a/src/jumpy/__pycache__/iterators.cpython-311.pyc and /dev/null differ diff --git a/src/jumpy/__pycache__/model.cpython-311.pyc b/src/jumpy/__pycache__/model.cpython-311.pyc deleted file mode 100644 index f41eba2..0000000 Binary files a/src/jumpy/__pycache__/model.cpython-311.pyc and /dev/null differ diff --git a/src/jumpy/__pycache__/serialize.cpython-311.pyc b/src/jumpy/__pycache__/serialize.cpython-311.pyc deleted file mode 100644 index 0099ae7..0000000 Binary files a/src/jumpy/__pycache__/serialize.cpython-311.pyc and /dev/null differ diff --git a/src/jumpy/backend.py b/src/jumpy/backend.py new file mode 100644 index 0000000..dfff732 --- /dev/null +++ b/src/jumpy/backend.py @@ -0,0 +1,140 @@ +""" +Solver backend abstraction. + +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). +""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from jumpy.model import Model + + +class Backend(ABC): + """Abstract solver backend.""" + + @abstractmethod + def optimize(self, model: Model) -> list[float]: + """Solve the model and return the solution vector.""" + ... + + +class JuliacBackend(Backend): + """ + Default backend: calls a precompiled Julia shared library via ctypes. + + The library is built with juliac from: + MOI + GenOpt + Bridges + HiGHS + + No Julia installation required. + """ + + 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 + 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" + " 1. Install the pre-built wheel: pip install jumpy\n" + " 2. 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." + ) + + +class JuliaCallBackend(Backend): + """ + Optional backend: calls Julia directly through juliacall. + + Requires `pip install jumpy[juliacall]`. Julia is installed lazily + by juliacall on first use if not already present. + + This backend has full flexibility — it can use any solver or MOI + feature, not just what's compiled into the juliac library. + """ + + 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())}" + ) + return cls() diff --git a/src/jumpy/bridge_juliacall.py b/src/jumpy/bridge_juliacall.py new file mode 100644 index 0000000..2126c8a --- /dev/null +++ b/src/jumpy/bridge_juliacall.py @@ -0,0 +1,391 @@ +""" +Bridge between JuMPy's Python expression graph and MathOptInterface via juliacall. + +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. +""" + +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_scalar_affine(terms, constant) + return MOI.ScalarAffineFunction(terms, constant) + end + + function _jumpy_affine_term(coef, var) + return MOI.ScalarAffineTerm(coef, var) + end + + function _jumpy_set_objective!(optimizer, sense, func) + MOI.set(optimizer, MOI.ObjectiveSense(), sense) + F = typeof(func) + MOI.set(optimizer, MOI.ObjectiveFunction{F}(), func) + end + """) + jl.seval(""" + function _jumpy_add_variables!(optimizer, count, lower, upper) + 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 + 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") + 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 + + function _jumpy_get_solution(optimizer, vars) + return [MOI.get(optimizer, MOI.VariablePrimal(), v) for v in vars] + end + """) + _HELPERS_DEFINED = True + + +def build_moi_model(jl, model: Model) -> list[float]: + """ + Build an MOI model in Julia from a JuMPy Model and solve it. + + Constraint groups are passed as GenOpt.FunctionGenerator objects + so that expansion happens entirely in Julia. + """ + _define_helpers(jl) + + optimizer = jl._jumpy_create_optimizer() + + # 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, + ) + 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 + ) + obj_func = _expr_to_moi(jl, all_jl_vars, model._objective.expr) + jl._jumpy_set_objective_b(optimizer, sense, 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._jumpy_affine_term(float(coef), jl_var)) + + return jl._jumpy_scalar_affine(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_) diff --git a/src/jumpy/model.py b/src/jumpy/model.py index 1024fea..e7dce37 100644 --- a/src/jumpy/model.py +++ b/src/jumpy/model.py @@ -20,6 +20,7 @@ ) 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: @@ -101,7 +102,15 @@ class Model: m.optimize() """ - def __init__(self): + def __init__(self, backend: str = "juliac"): + """ + Create a new model. + + Args: + 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] = [] @@ -191,17 +200,12 @@ def objective(self, obj: Objective) -> None: def optimize(self) -> None: """ - Serialize the model and send it to the compiled Julia library for solving. - - The Julia side: - 1. Reconstructs MOI.ScalarNonlinearFunction trees from flat arrays - 2. Wraps constraint groups in GeneratorOptInterface.IteratedFunction - 3. Adds bridges - 4. Calls HiGHS - 5. Returns the solution vector + 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 """ - model_data = self._serialize() - self._solution = _call_julia_solver(model_data) + self._solution = self._backend.optimize(self) def _serialize(self) -> dict: """Serialize the entire model for the Julia C ABI.""" @@ -263,13 +267,3 @@ def value(self, var: Variable) -> float: return self._solution[var.index] -def _call_julia_solver(model_data: dict) -> list[float]: - """ - Call the compiled Julia shared library via ctypes. - - TODO: Implement once the Julia library is compiled with juliac. - """ - raise NotImplementedError( - "Julia solver backend not yet available. " - "Compile the Julia library with juliac and place it on the library path." - ) diff --git a/tests/test_expressions.py b/tests/test_expressions.py index a70b87c..65e255d 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -267,3 +267,5 @@ def test_serialize_full_model(): 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 new file mode 100644 index 0000000..effc976 --- /dev/null +++ b/tests/test_solve.py @@ -0,0 +1,130 @@ +""" +End-to-end tests that solve models with HiGHS via the juliacall backend. + +These require Julia + juliacall to be installed: + pip install juliacall +""" + +import sys +sys.path.insert(0, "src") + +from jumpy import Model, Iterator, Parameter, minimize, maximize, sin, exp + + +def _model(): + return Model(backend="juliacall") + + +def test_simple_lp(): + """min x + y s.t. x + y >= 10, x >= 0, y >= 0""" + m = _model() + x = m.variable(lower=0, name="x") + y = m.variable(lower=0, name="y") + + m.constraint(x + y >= 10) + m.objective = minimize(x + y) + m.optimize() + + assert abs(m.value(x) + m.value(y) - 10.0) < 1e-6 + + +def test_constraint_group_lp(): + """ + min sum(x) s.t. x[i] >= 1 for i in 0..9, x[i] >= 0 + Optimal: all x[i] = 1, obj = 10 + """ + m = _model() + x = m.variables(10, lower=0, name="x") + + i = Iterator(range(10)) + m.constraint_group([i], x[i] >= 1) + + m.objective = minimize(sum(x)) + m.optimize() + + total = sum(m.value(v) for v in x) + assert abs(total - 10.0) < 1e-6 + + +def test_constraint_group_consecutive(): + """ + min x[0] 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) + + m.objective = minimize(x[0] + x[1] + x[2]) + m.optimize() + + for k in range(9): + assert m.value(x[k]) + m.value(x[k + 1]) >= 2.0 - 1e-6 + + +def test_parameter_in_constraint_group(): + """ + min sum(x) s.t. x[i] >= demand[i], x[i] >= 0 + demand = [1, 2, 3, 4, 5] + Optimal: x[i] = demand[i], obj = 15 + """ + m = _model() + x = m.variables(5, lower=0, name="x") + demand = 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]) + + m.objective = minimize(sum(x)) + m.optimize() + + total = sum(m.value(v) for v in x) + assert abs(total - 15.0) < 1e-6 + + +def test_multidim_constraint_group(): + """ + 2D indexing: x[3*i + j] >= 1 for i in 0..2, j in 0..2 + 9 variables, all >= 1 + """ + 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) + + m.objective = minimize(sum(x)) + m.optimize() + + total = sum(m.value(v) for v in x) + assert abs(total - 9.0) < 1e-6 + + +def test_maximize(): + """max x s.t. x <= 42, x >= 0""" + m = _model() + x = m.variable(lower=0, upper=42, name="x") + + m.objective = maximize(x) + m.optimize() + + assert abs(m.value(x) - 42.0) < 1e-6 + + +if __name__ == "__main__": + 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 + import traceback + print(f" FAIL {test.__name__}: {e}") + traceback.print_exc() + print(f"\n{passed} passed, {failed} failed") + sys.exit(1 if failed else 0)