diff --git a/doc/api.rst b/doc/api.rst index 20958857..1554ce60 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -24,6 +24,7 @@ Creating a model piecewise.segments model.Model.linexpr model.Model.remove_constraints + model.Model.copy Classes under the hook diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 35b21c67..26007b5c 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -4,6 +4,7 @@ Release Notes Upcoming Version ---------------- +* Add ``Model.copy()`` (default deep copy) with ``deep`` and ``include_solution`` options; support Python ``copy.copy`` and ``copy.deepcopy`` protocols via ``__copy__`` and ``__deepcopy__``. * Harmonize coordinate alignment for operations with subset/superset objects: - Multiplication and division fill missing coords with 0 (variable doesn't participate) - Addition and subtraction of constants fill missing coords with 0 (identity element) and pin result to LHS coords diff --git a/linopy/dual.py b/linopy/dual.py new file mode 100644 index 00000000..0a787e68 --- /dev/null +++ b/linopy/dual.py @@ -0,0 +1,581 @@ +""" +Linopy dual module. + +This module contains implementations for constructing the dual of a linear optimization problem. +""" + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING, Any + +import numpy as np +import pandas as pd +import xarray as xr + +from linopy.expressions import LinearExpression + +if TYPE_CHECKING: + from linopy.model import Model + +logger = logging.getLogger(__name__) + + +def _var_lookup(m: Model) -> dict: + """ + Build a flat label -> (var_name, coord_dict) lookup for all variables in m. + + Used to map entries in m.matrices.vlabels back to their variable name + and xarray coordinates for use in dual feasibility constraint construction. + + Skips masked entries (label == -1) and empty variables. + + Parameters + ---------- + m : Model + Primal linopy model. + + Returns + ------- + dict + Mapping from flat integer label to (var_name, coord_dict) tuple. + """ + var_lookup = {} + logger.debug("Building variable label lookup.") + for var_name, var in m.variables.items(): + labels = var.labels + flat_labels = labels.values.flatten() + + if len(flat_labels) == 0: + logger.debug(f"Skipping empty variable '{var_name}'.") + continue + if not (flat_labels != -1).any(): + logger.debug(f"Variable '{var_name}' is fully masked, skipping.") + continue + + logger.debug( + f"Creating label lookup for variable '{var_name}' with shape {labels.shape} and dims {labels.dims}." + ) + + coord_arrays = ( + np.meshgrid( + *[labels.coords[dim].values for dim in labels.dims], indexing="ij" + ) + if len(labels.dims) > 0 + else [] + ) + flat_coords = [arr.flatten() for arr in coord_arrays] + + for k, flat in enumerate(flat_labels): + if flat != -1: + var_lookup[int(flat)] = ( + var_name, + {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, + ) + + return var_lookup + + +def _con_lookup(m: Model) -> dict: + """ + Build a flat label -> (con_name, coord_dict) lookup for all constraints in m. + + Used to map entries in m.matrices.clabels back to their constraint name + and xarray coordinates for use in dual feasibility constraint construction. + + Skips masked entries (label == -1) and empty or fully-masked constraints. + + Parameters + ---------- + m : Model + Primal linopy model. + + Returns + ------- + dict + Mapping from flat integer label to (con_name, coord_dict) tuple. + """ + con_lookup = {} + logger.debug("Building constraint label lookup.") + for con_name, con in m.constraints.items(): + labels = con.labels + flat_labels = labels.values.flatten() + + if len(flat_labels) == 0: + logger.debug(f"Skipping empty constraint '{con_name}'.") + continue + if not (flat_labels != -1).any(): + logger.debug(f"Constraint '{con_name}' is fully masked, skipping.") + continue + + logger.debug( + f"Creating label lookup for constraint '{con_name}' with shape {labels.shape} and dims {labels.dims}." + ) + + coord_arrays = ( + np.meshgrid( + *[labels.coords[dim].values for dim in labels.dims], indexing="ij" + ) + if len(labels.dims) > 0 + else [] + ) + flat_coords = [arr.flatten() for arr in coord_arrays] + + for k, flat in enumerate(flat_labels): + if flat != -1: + con_lookup[int(flat)] = ( + con_name, + {dim: flat_coords[i][k] for i, dim in enumerate(labels.dims)}, + ) + + return con_lookup + + +def bounds_to_constraints(m: Model) -> None: + """ + Add explicit bound constraints for variables with bounds set directly + in the variable rather than via explicit constraints. + + Adds constraints named '{var_name}-bound-lower' and '{var_name}-bound-upper' + to distinguish from PyPSA's automatic '-fix-*' constraints. + + Also resets variable bounds to [-inf, inf] after adding constraints, + to avoid double-counting in the dual. + + Parameters + ---------- + m : Model + Linopy model to convert variable bounds to constraints. Mutates the model in-place. + """ + logger.debug("Converting variable bounds to explicit constraints.") + logger.debug("Relaxing variable bounds to [-inf, inf].") + for var_name, var in m.variables.items(): + mask = var.labels != -1 + lb = var.lower + ub = var.upper + + # lower bound + if f"{var_name}-bound-lower" not in m.constraints: + has_finite_lb = np.isfinite(lb.values[mask.values]).any() + if has_finite_lb: + m.add_constraints( + var >= lb, + name=f"{var_name}-bound-lower", + mask=mask, + ) + logger.debug(f"Added lower bound constraint for '{var_name}'.") + var.lower.values[mask.values] = -np.inf + # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. + m.variables[var_name].lower.values[mask.values] = -np.inf + else: + logger.debug( + f"Variable '{var_name}' has no finite lower bound, skipping." + ) + + # upper bound + if f"{var_name}-bound-upper" not in m.constraints: + has_finite_ub = np.isfinite(ub.values[mask.values]).any() + if has_finite_ub: + m.add_constraints( + var <= ub, + name=f"{var_name}-bound-upper", + mask=mask, + ) + logger.debug(f"Added upper bound constraint for '{var_name}'.") + var.upper.values[mask.values] = np.inf + # Remove bounds to avoid double-counting in the dual. Rely on the new constraints instead. + m.variables[var_name].upper.values[mask.values] = np.inf + else: + logger.debug( + f"Variable '{var_name}' has no finite upper bound, skipping." + ) + + +def _add_dual_variables(m: Model, m2: Model) -> dict: + """ + Add dual variables to m2 corresponding to constraints in m.. + + For each active constraint in m, adds a dual variable to m2 following + standard LP duality sign conventions. The sign of the dual variable bounds + depends on both the constraint type and the primal objective sense: + + For a minimization primal: + - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) + - <= constraints -> non-positive dual variable (lower=-inf, upper=0) + - >= constraints -> non-negative dual variable (lower=0, upper=inf) + + For a maximization primal: + - Equality constraints (=) -> free dual variable (lower=-inf, upper=inf) + - <= constraints -> non-negative dual variable (lower=0, upper=inf) + - >= constraints -> non-positive dual variable (lower=-inf, upper=0) + + Skips constraints with no active rows (empty or fully masked). + + Parameters + ---------- + m : Model + Primal linopy model containing the constraints to dualize. + m2 : Model + Dual linopy model to which dual variables are added. + + Returns + ------- + dict + Mapping from constraint name (str) to the corresponding dual + variable (linopy.Variable) in m2. + """ + primal_is_min = m.objective.sense == "min" + + dual_vars = {} + for name, con in m.constraints.items(): + sign_vals = con.sign.values.flatten() + + if len(sign_vals) == 0: + logger.warning(f"Constraint '{name}' has no sign values, skipping.") + continue + + mask = con.labels != -1 + if not mask.any(): + logger.debug(f"Constraint '{name}' is fully masked, skipping.") + continue + + if sign_vals[0] == "=": + lower, upper = -np.inf, np.inf + var_type = "free" + elif sign_vals[0] == "<=": + lower, upper = (-np.inf, 0) if primal_is_min else (0, np.inf) + var_type = "non-positive" if primal_is_min else "non-negative" + elif sign_vals[0] == ">=": + lower, upper = (0, np.inf) if primal_is_min else (-np.inf, 0) + var_type = "non-negative" if primal_is_min else "non-positive" + else: + logger.warning( + f"Constraint '{name}' has unrecognized sign '{sign_vals[0]}', skipping." + ) + continue + + logger.debug( + f"Adding {var_type} dual variable for constraint '{name}' with shape {con.shape} and dims {con.labels.dims}." + ) + coords = ( + [con.labels.coords[dim] for dim in con.labels.dims] + if con.labels.dims + else None + ) + dual_vars[name] = m2.add_variables( + lower=lower, + upper=upper, + coords=coords, + name=name, + mask=mask, + ) + + return dual_vars + + +def _build_dual_feas_terms( + m: Model, + dual_vars: dict, + var_lookup: dict, + con_lookup: dict, +) -> dict: + """ + Build dual feasibility terms for each primal variable in m. + + For each active primal variable x_j, collects the constraint matrix + entries A_ji and their corresponding constraint names and coordinates, + forming the terms of the stationarity condition: + sum_i (A_ji * lambda_i) = c_j + + Raw constraint matrix coefficients are used directly without sign + factors, as the sign convention is encoded in the dual variable bounds: + - <= constraints: lambda_i <= 0 + - >= constraints: lambda_i >= 0 + - = constraints: lambda_i free + + Parameters + ---------- + m : Model + Primal linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). Used to skip constraints + that were not dualized (e.g. empty or fully masked). + var_lookup : dict + Mapping from flat variable label to (var_name, coord_dict), + as returned by _var_lookup(). + con_lookup : dict + Mapping from flat constraint label to (con_name, coord_dict), + as returned by _con_lookup(). + + Returns + ------- + dict + Nested dict: {var_name: {flat_label: (var_coords, terms, obj_coeff)}} + where terms is a list of (con_name, con_coords, coeff) tuples. + """ + A = m.matrices.A + if A is None: + raise ValueError("Constraint matrix is None, model has no constraints.") + A_csc = A.tocsc() + c = m.matrices.c + indptr = A_csc.indptr + indices = A_csc.indices + data = A_csc.data + vlabels = m.matrices.vlabels + clabels = m.matrices.clabels + + dual_feas_terms: dict[str, dict[int, tuple]] = { + var_name: {} for var_name in m.variables + } + + logger.debug("Building dual feasibility terms for each primal variable.") + + for i in range(A_csc.shape[1]): + flat_var = vlabels[i] + if flat_var == -1: + continue + if flat_var not in var_lookup: + continue + var_name, var_coords = var_lookup[flat_var] + terms = [] + for k in range(indptr[i], indptr[i + 1]): + j = indices[k] + flat_con = clabels[j] + if flat_con == -1: + continue + if flat_con not in con_lookup: + continue + con_name, con_coords = con_lookup[flat_con] + if con_name not in dual_vars: + continue + coeff = data[k] + terms.append((con_name, con_coords, coeff)) + dual_feas_terms[var_name][flat_var] = (var_coords, terms, c[i]) + + return dual_feas_terms + + +def _add_dual_feasibility_constraints( + m: Model, + m2: Model, + dual_vars: dict, + var_lookup: dict, + con_lookup: dict, +) -> None: + """ + Add dual feasibility constraints to m2. + + For each primal variable x_j in m, adds the stationarity constraint: + sum_i (A_ji * lambda_i) = c_j + where: + - A is the primal constraint matrix + - lambda_i are the dual variables in m2 + - c_j is the objective coefficient of x_j + + Raw constraint matrix coefficients are used directly without sign factors, + because the sign convention is encoded in the dual variable bounds: + - <= constraints: lambda_i <= 0 + - >= constraints: lambda_i >= 0 + - = constraints: lambda_i free + + Skips masked variable entries (label == -1) and variables not present + in var_lookup (e.g. from empty constraints). + + Parameters + ---------- + m : Model + Primal linopy model. + m2 : Model + Dual linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). + var_lookup : dict + Mapping from flat variable label to (var_name, coord_dict), + as returned by _var_lookup(). + con_lookup : dict + Mapping from flat constraint label to (con_name, coord_dict), + as returned by _con_lookup(). + """ + dual_feas_terms = _build_dual_feas_terms(m, dual_vars, var_lookup, con_lookup) + + c = m.matrices.c + vlabels = m.matrices.vlabels + + # build objective coefficient lookup by flat variable label + c_by_label = {vlabels[i]: c[i] for i in range(len(vlabels))} + + # add dual feasibility constraints to m2 + logger.debug("Adding dual feasibility constraints to model.") + for var_name, var in m.variables.items(): + coords = [ + pd.Index(var.labels.coords[dim].values, name=dim) for dim in var.labels.dims + ] + mask = var.labels != -1 + + c_vals = xr.DataArray( + np.vectorize(lambda flat: c_by_label.get(flat, 0.0))(var.labels.values), + coords=var.labels.coords, + ) + + def rule( + m: Model, + *coord_vals: Any, + vname: str = var_name, + vdims: tuple = var.labels.dims, + ) -> LinearExpression | None: + coord_dict = dict(zip(vdims, coord_vals)) + flat = var.labels.sel(**coord_dict).item() + if flat == -1: + return None + if flat not in dual_feas_terms[vname]: + return None + _, terms, _ = dual_feas_terms[vname][flat] + if not terms: + return None + return sum( + coeff * dual_vars[con_name].at[tuple(con_coords.values())] + for con_name, con_coords, coeff in terms + ) + + lhs = LinearExpression.from_rule(m2, rule, coords) + m2.add_constraints(lhs == c_vals, name=var_name, mask=mask) + + +def _add_dual_objective( + m: Model, + m2: Model, + dual_vars: dict, + add_objective_constant: float = 0.0, +) -> None: + """ + Construct and add the dual objective to m2. + + The dual objective is sum(rhs * dual) over all constraints, added uniformly + with a + sign. The sign convention is encoded in the dual variable bounds: + - <= constraints: dual <= 0, so rhs * dual contributes negatively + - >= constraints: dual >= 0, so rhs * dual contributes positively + - = constraints: dual free + + This matches linopy's and Gurobi's native dual sign convention, allowing + direct comparison between m2.variables[con_name].solution and + m.constraints[con_name].dual without sign adjustments. + + The dual objective sense is flipped relative to the primal: + - min primal -> max dual + - max primal -> min dual + + Parameters + ---------- + m : Model + Primal linopy model. + m2 : Model + Dual linopy model. + dual_vars : dict + Mapping from constraint name to dual variable in m2, + as returned by _add_dual_variables(). + add_objective_constant : float, optional + Constant term to add to the dual objective. Use this to pass through + a primal objective constant excluded via include_objective_constant=False + during model creation. Default is 0.0. + """ + dual_obj: LinearExpression = LinearExpression(None, m2) + sense = "max" if m.objective.sense == "min" else "min" + + for name, con in m.constraints.items(): + if name not in dual_vars: + continue + + mask = con.labels != -1 + rhs_masked = con.rhs.where(mask, 0) + dual_obj += (rhs_masked * dual_vars[name]).sum() + + if add_objective_constant != 0.0: + dual_obj += add_objective_constant + logger.debug(f"Added constant {add_objective_constant} to dual objective.") + + logger.debug(f"Constructed dual objective with {len(dual_obj.coeffs)} terms.") + logger.debug("Adding dual objective to model.") + m2.add_objective(dual_obj, sense=sense, overwrite=True) + + +def dualize( + m: Model, + add_objective_constant: float = 0.0, +) -> Model: + """ + Construct the dual of a linopy LP model. + + Transforms the primal model into its dual equivalent m2 following + standard LP duality theory. The dual sense is flipped relative to the + primal (min -> max, max -> min), and dual variable bounds depend on + both constraint type and primal objective sense. + + For a minimization primal: + + Primal (min): Dual (max): + min c^T x max b_eq^T λ + b_leq^T μ + b_geq^T ν + s.t. A_eq x = b_eq : λ free s.t. A_eq^T λ + A_leq^T μ + A_geq^T ν = c + A_leq x <= b_leq : μ <= 0 λ free, μ <= 0, ν >= 0 + A_geq x >= b_geq : ν >= 0 + + For a maximization primal the dual variable bounds are flipped: + μ >= 0 for <= constraints, ν <= 0 for >= constraints. + + Variable bounds are converted to explicit constraints before dualization + via bounds_to_constraints(), so that they appear in the constraint matrix + A and are correctly reflected in the dual. + + The dual variables in m2 are named identically to their corresponding + primal constraints and are accessible via m2.variables[con_name]. + + Strong duality guarantees that at optimality: + primal objective = dual objective + + Note: The standalone dual m2 may be unbounded if the primal is degenerate. + Only linear programs (LP) are supported. + + Parameters + ---------- + m : Model + Primal linopy model to dualize. Must have a linear objective and linear constraints. + + add_objective_constant : float, optional + Constant term to add to the dual objective. Use this to pass through + a primal objective constant. Default is 0.0. + + Returns + ------- + Model + The dual linopy model. Dual variables are named after their + corresponding primal constraints. + + Examples + -------- + .. code-block:: python + + m2 = m.dualize() + m2.solve(solver_name="gurobi", Method=2, Crossover=1) + gap = abs(m.objective.value - m2.objective.value) + """ + from linopy.model import Model + + m1 = m.copy() + m2 = Model() + + if not m.variables or not m.constraints: + logger.warning( + "Primal model has no variables or constraints. Returning empty dual model." + ) + return m2 + + bounds_to_constraints(m1) + var_lup = _var_lookup(m1) + con_lup = _con_lookup(m1) + dual_vars = _add_dual_variables(m1, m2) + _add_dual_feasibility_constraints(m1, m2, dual_vars, var_lup, con_lup) + _add_dual_objective( + m1, m2, dual_vars, add_objective_constant=add_objective_constant + ) + return m2 diff --git a/linopy/io.py b/linopy/io.py index 2213cbb5..e7353b60 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -1239,3 +1239,124 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset: setattr(m, k, ds.attrs.get(k)) return m + + +def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: + """ + Return a copy of this model. + + With ``deep=True`` (default), variables, constraints, objective, + parameters, blocks, and scalar attributes are copied to a fully + independent model. With ``deep=False``, returns a shallow copy. + + :meth:`Model.copy` defaults to deep copy for workflow safety. + In contrast, ``copy.copy(model)`` is shallow via ``__copy__``, and + ``copy.deepcopy(model)`` is deep via ``__deepcopy__``. + + Solver runtime metadata (for example, ``solver_name`` and + ``solver_model``) is intentionally not copied. Solver backend state + is recreated on ``solve()``. + + Parameters + ---------- + m : Model + The model to copy. + include_solution : bool, optional + Whether to include solution and dual values in the copy. + If False (default), solve artifacts are excluded: solution/dual data, + objective value, and solve status are reset to initialized state. + If True, these values are copied when present. For unsolved models, + this has no additional effect. + deep : bool, optional + Whether to return a deep copy (default) or shallow copy. If False, + the returned model uses independent wrapper objects that share + underlying data buffers with the source model. + + Returns + ------- + Model + A deep or shallow copy of the model. + """ + from linopy.model import ( + Constraint, + Constraints, + LinearExpression, + Model, + Objective, + Variable, + Variables, + ) + + SOLVE_STATE_ATTRS = {"status", "termination_condition"} + + new_model = Model( + chunk=m._chunk, + force_dim_names=m._force_dim_names, + auto_mask=m._auto_mask, + solver_dir=str(m._solver_dir), + ) + + new_model._variables = Variables( + { + name: Variable( + var.data.copy(deep=deep) + if include_solution + else var.data[m.variables.dataset_attrs].copy(deep=deep), + new_model, + name, + ) + for name, var in m.variables.items() + }, + new_model, + ) + + new_model._constraints = Constraints( + { + name: Constraint( + con.data.copy(deep=deep) + if include_solution + else con.data[m.constraints.dataset_attrs].copy(deep=deep), + new_model, + name, + ) + for name, con in m.constraints.items() + }, + new_model, + ) + + obj_expr = LinearExpression(m.objective.expression.data.copy(deep=deep), new_model) + new_model._objective = Objective(obj_expr, new_model, m.objective.sense) + new_model._objective._value = m.objective.value if include_solution else None + + new_model._parameters = m._parameters.copy(deep=deep) + new_model._blocks = m._blocks.copy(deep=deep) if m._blocks is not None else None + + for attr in m.scalar_attrs: + if include_solution or attr not in SOLVE_STATE_ATTRS: + setattr(new_model, attr, getattr(m, attr)) + + return new_model + + +def shallowcopy(m: Model) -> Model: + """ + Support Python's ``copy.copy`` protocol for ``Model``. + + Returns a shallow copy with independent wrapper objects that share + underlying array buffers with ``m``. Solve artifacts are excluded, + matching :meth:`Model.copy` defaults. + """ + return copy(m, include_solution=False, deep=False) + + +def deepcopy(m: Model, memo: dict[int, Any]) -> Model: + """ + Support Python's ``copy.deepcopy`` protocol for ``Model``. + + Returns a deep, structurally independent copy and records it in ``memo`` + as required by Python's copy protocol. Solve artifacts are excluded, + matching :meth:`Model.copy` defaults. + """ + new_model = copy(m, include_solution=False, deep=True) + memo[id(m)] = new_model + return new_model diff --git a/linopy/model.py b/linopy/model.py index 06e814c6..404df054 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -47,12 +47,16 @@ TerminationCondition, ) from linopy.constraints import AnonymousScalarConstraint, Constraint, Constraints +from linopy.dual import bounds_to_constraints, dualize from linopy.expressions import ( LinearExpression, QuadraticExpression, ScalarLinearExpression, ) from linopy.io import ( + copy, + deepcopy, + shallowcopy, to_block_files, to_cupdlpx, to_file, @@ -1877,6 +1881,12 @@ def reset_solution(self) -> None: self.variables.reset_solution() self.constraints.reset_dual() + copy = copy + + __copy__ = shallowcopy + + __deepcopy__ = deepcopy + to_netcdf = to_netcdf to_file = to_file @@ -1890,3 +1900,7 @@ def reset_solution(self) -> None: to_cupdlpx = to_cupdlpx to_block_files = to_block_files + + bounds_to_constraints = bounds_to_constraints + + dualize = dualize diff --git a/test/test_model.py b/test/test_model.py index c363fe4c..c0988c26 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -5,6 +5,7 @@ from __future__ import annotations +import copy as pycopy from pathlib import Path from tempfile import gettempdir @@ -12,8 +13,13 @@ import pytest import xarray as xr -from linopy import EQUAL, Model -from linopy.testing import assert_model_equal +from linopy import EQUAL, Model, available_solvers +from linopy.testing import ( + assert_conequal, + assert_equal, + assert_linequal, + assert_model_equal, +) target_shape: tuple[int, int] = (10, 10) @@ -163,3 +169,167 @@ def test_assert_model_equal() -> None: m.add_objective(obj) assert_model_equal(m, m) + + +@pytest.fixture(scope="module") +def copy_test_model() -> Model: + """Small representative model used across copy tests.""" + m: Model = Model() + + lower: xr.DataArray = xr.DataArray( + np.zeros((10, 10)), coords=[range(10), range(10)] + ) + upper: xr.DataArray = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)]) + x = m.add_variables(lower, upper, name="x") + y = m.add_variables(name="y") + + m.add_constraints(1 * x + 10 * y, EQUAL, 0) + m.add_objective((10 * x + 5 * y).sum()) + + return m + + +@pytest.fixture(scope="module") +def solved_copy_test_model(copy_test_model: Model) -> Model: + """Solved representative model used across solved-copy tests.""" + m = copy_test_model.copy(deep=True) + m.solve() + return m + + +def test_model_copy_unsolved(copy_test_model: Model) -> None: + """Copy of unsolved model is structurally equal and independent.""" + m = copy_test_model.copy(deep=True) + c = m.copy(include_solution=False) + + assert_model_equal(m, c) + + # independence: mutating copy does not affect source + c.add_variables(name="z") + assert "z" not in m.variables + + +def test_model_copy_unsolved_with_solution_flag(copy_test_model: Model) -> None: + """Unsolved model with include_solution=True has no extra solve artifacts.""" + m = copy_test_model.copy(deep=True) + + c_include_solution = m.copy(include_solution=True) + c_exclude_solution = m.copy(include_solution=False) + + assert_model_equal(c_include_solution, c_exclude_solution) + assert c_include_solution.status == "initialized" + assert c_include_solution.termination_condition == "" + assert c_include_solution.objective.value is None + + +def test_model_copy_shallow(copy_test_model: Model) -> None: + """Shallow copy has independent wrappers sharing underlying data buffers.""" + m = copy_test_model.copy(deep=True) + c = m.copy(deep=False) + + assert c is not m + assert c.variables is not m.variables + assert c.constraints is not m.constraints + assert c.objective is not m.objective + + # wrappers are distinct, but shallow copy shares payload buffers + c.variables["x"].lower.values[0, 0] = 123.0 + assert m.variables["x"].lower.values[0, 0] == 123.0 + + +def test_model_deepcopy_protocol(copy_test_model: Model) -> None: + """copy.deepcopy(model) dispatches to Model.__deepcopy__ and stays independent.""" + m = copy_test_model.copy(deep=True) + c = pycopy.deepcopy(m) + + assert_model_equal(m, c) + + # Test independence: mutations to copy do not affect source + # 1. Variable mutation: add new variable + c.add_variables(name="z") + assert "z" not in m.variables + + # 2. Variable data mutation (bounds): verify buffers are independent + original_lower = m.variables["x"].lower.values[0, 0].item() + new_lower = 999 + c.variables["x"].lower.values[0, 0] = new_lower + assert c.variables["x"].lower.values[0, 0] == new_lower + assert m.variables["x"].lower.values[0, 0] == original_lower + + # 3. Constraint coefficient mutation: deep copy must not leak back + original_con_coeff = m.constraints["con0"].coeffs.values.flat[0].item() + new_con_coeff = original_con_coeff + 42 + c.constraints["con0"].coeffs.values.flat[0] = new_con_coeff + assert c.constraints["con0"].coeffs.values.flat[0] == new_con_coeff + assert m.constraints["con0"].coeffs.values.flat[0] == original_con_coeff + + # 4. Objective expression coefficient mutation: deep copy must not leak back + original_obj_coeff = m.objective.expression.coeffs.values.flat[0].item() + new_obj_coeff = original_obj_coeff + 20 + c.objective.expression.coeffs.values.flat[0] = new_obj_coeff + assert c.objective.expression.coeffs.values.flat[0] == new_obj_coeff + assert m.objective.expression.coeffs.values.flat[0] == original_obj_coeff + + # 5. Objective sense mutation + original_sense = m.objective.sense + c.objective.sense = "max" + assert c.objective.sense == "max" + assert m.objective.sense == original_sense + + +@pytest.mark.skipif(not available_solvers, reason="No solver installed") +class TestModelCopySolved: + def test_model_deepcopy_protocol_excludes_solution( + self, solved_copy_test_model: Model + ) -> None: + """copy.deepcopy on solved model drops solve state by default.""" + m = solved_copy_test_model + + c = pycopy.deepcopy(m) + + assert c.status == "initialized" + assert c.termination_condition == "" + assert c.objective.value is None + + for v in m.variables: + assert_equal( + c.variables[v].data[c.variables.dataset_attrs], + m.variables[v].data[m.variables.dataset_attrs], + ) + for con in m.constraints: + assert_conequal(c.constraints[con], m.constraints[con], strict=False) + assert_linequal(c.objective.expression, m.objective.expression) + assert c.objective.sense == m.objective.sense + + def test_model_copy_solved_with_solution( + self, solved_copy_test_model: Model + ) -> None: + """Copy with include_solution=True preserves solve state.""" + m = solved_copy_test_model + + c = m.copy(include_solution=True) + assert_model_equal(m, c) + + def test_model_copy_solved_without_solution( + self, solved_copy_test_model: Model + ) -> None: + """Copy with include_solution=False (default) drops solve state but preserves problem structure.""" + m = solved_copy_test_model + + c = m.copy(include_solution=False) + + # solve state is dropped + assert c.status == "initialized" + assert c.termination_condition == "" + assert c.objective.value is None + + # problem structure is preserved — compare only dataset_attrs to exclude solution/dual + for v in m.variables: + assert_equal( + c.variables[v].data[c.variables.dataset_attrs], + m.variables[v].data[m.variables.dataset_attrs], + ) + for con in m.constraints: + assert_conequal(c.constraints[con], m.constraints[con], strict=False) + assert_linequal(c.objective.expression, m.objective.expression) + assert c.objective.sense == m.objective.sense