From ced2c580e7d7e55e1c76d1c72f85fca91c1b2acf Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:04:19 -0500 Subject: [PATCH 01/17] Various fixes --- src/pyoframe/_arithmetic.py | 3 ++- src/pyoframe/core.py | 41 +++++++++++++++++++++++++++++++++---- src/pyoframe/util.py | 14 +++++++++++-- tests/test_arithmetic.py | 24 +++++++++++----------- tests/test_solver.py | 16 +++++++-------- 5 files changed, 71 insertions(+), 27 deletions(-) diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index 00b10c9..3491582 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -116,7 +116,7 @@ def get_indices(expr): not Config.disable_unmatched_checks ), "This code should not be reached when unmatched checks are disabled." outer_join = get_indices(left).join( - get_indices(right), how="outer", on=dims + get_indices(right), how="full", on=dims ) if outer_join.get_column(dims[0]).null_count() > 0: raise PyoframeError( @@ -188,6 +188,7 @@ def _add_dimension(self: "Expression", target: "Expression") -> "Expression": return self if not set(missing_dims) <= set(self.allowed_new_dims): + # TODO actually suggest using e.g. .add_dim("a", "b") instead of just "use .add_dim()" raise PyoframeError( f"Dataframe has missing dimensions {missing_dims}. If this is intentional, use .add_dim()\n{self.data}" ) diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 294eb19..fcdcca9 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -126,6 +126,32 @@ def __rmul__(self, other): def __radd__(self, other): return self.to_expr() + other + def __truediv__(self, other): + """ + Support division. + >>> from pyoframe import Variable + >>> var = Variable({"dim1": [1,2,3]}) + >>> var / 2 + + [1]: 0.5 x1 + [2]: 0.5 x2 + [3]: 0.5 x3 + """ + return self.to_expr() * (1 / other) + + def __rsub__(self, other): + """ + Support right subtraction. + >>> from pyoframe import Variable + >>> var = Variable({"dim1": [1,2,3]}) + >>> 1 - var + + [1]: 1 - x1 + [2]: 1 - x2 + [3]: 1 - x3 + """ + return other + (-self.to_expr()) + def __le__(self, other): """Equality constraint. Examples @@ -248,7 +274,7 @@ def __add__(self, other): return self._new( pl.concat([self.data, other.data]).unique(maintain_order=True) ) - except pl.ShapeError as e: + except pl.exceptions.ShapeError as e: if "unable to vstack, column names don't match" in str(e): raise PyoframeError( f"Failed to add sets '{self.friendly_name}' and '{other.friendly_name}' because dimensions do not match ({self.dimensions} != {other.dimensions}) " @@ -295,6 +321,10 @@ def _set_to_polars(set: "SetTypes") -> pl.DataFrame: df = set.to_frame() elif isinstance(set, Set): df = set.data + elif isinstance(set, range): + raise ValueError( + "Cannot convert a range to a set without a dimension name. Try Set(dim_name=range(...))" + ) else: raise ValueError(f"Cannot convert type {type(set)} to a polars DataFrame") @@ -507,7 +537,9 @@ def rolling_sum(self, over: str, window_size: int): [ df.with_columns(pl.col(over).max()) for _, df in self.data.rolling( - index_column=over, period=f"{window_size}i", by=remaining_dims + index_column=over, + period=f"{window_size}i", + group_by=remaining_dims, ) ] ) @@ -659,7 +691,7 @@ def _add_const(self, const: int | float) -> Expression: .unique(maintain_order=True) .with_columns(pl.lit(CONST_TERM).alias(VAR_KEY).cast(VAR_TYPE)) ) - data = data.join(keys, on=dim + [VAR_KEY], how="outer_coalesce") + data = data.join(keys, on=dim + [VAR_KEY], how="full", coalesce=True) data = data.with_columns(pl.col(COEF_KEY).fill_null(0.0)) data = data.with_columns( @@ -678,7 +710,8 @@ def constant_terms(self): return constant_terms.join( self.data.select(dims).unique(maintain_order=True), on=dims, - how="outer_coalesce", + how="full", + coalesce=True, ).with_columns(pl.col(COEF_KEY).fill_null(0.0)) else: if len(constant_terms) == 0: diff --git a/src/pyoframe/util.py b/src/pyoframe/util.py index 37846dd..7293ad6 100644 --- a/src/pyoframe/util.py +++ b/src/pyoframe/util.py @@ -47,10 +47,20 @@ def parse_inputs_as_iterable( def _is_iterable(input: Union[Any, Iterable[Any]]) -> bool: - # Inspired from the polars library + # Inspired from the polars library, TODO: Consider using opposite check, i.e. equals list or tuple return isinstance(input, Iterable) and not isinstance( input, - (str, bytes, pl.DataFrame, pl.Series, pd.DataFrame, pd.Series, pd.Index, dict), + ( + str, + bytes, + pl.DataFrame, + pl.Series, + pd.DataFrame, + pd.Series, + pd.Index, + dict, + range, + ), ) diff --git a/tests/test_arithmetic.py b/tests/test_arithmetic.py index 82b9cce..e525ecf 100644 --- a/tests/test_arithmetic.py +++ b/tests/test_arithmetic.py @@ -51,7 +51,7 @@ def test_multiplication_no_common_dimensions(): VAR_KEY: [CONST_TERM] * 6, } ), - check_dtype=False, + check_dtypes=False, ) @@ -71,7 +71,7 @@ def test_within_set(): VAR_KEY: [1, 4], } ), - check_dtype=False, + check_dtypes=False, ) @@ -82,7 +82,7 @@ def test_filter_expression(): assert_frame_equal( result.data, pl.DataFrame({"dim1": [2], COEF_KEY: [2], VAR_KEY: [CONST_TERM]}), - check_dtype=False, + check_dtypes=False, ) @@ -92,7 +92,7 @@ def test_filter_constraint(): assert_frame_equal( result, pl.DataFrame({"dim1": [2], COEF_KEY: [2], VAR_KEY: [CONST_TERM]}), - check_dtype=False, + check_dtypes=False, ) @@ -106,7 +106,7 @@ def test_filter_variable(): def test_filter_set(): s = Set(x=[1, 2, 3]) result = s.filter(x=2) - assert_frame_equal(result.data, pl.DataFrame({"x": [2]}), check_dtype=False) + assert_frame_equal(result.data, pl.DataFrame({"x": [2]}), check_dtypes=False) def test_drops_na(): @@ -139,7 +139,7 @@ def test_add_expressions(): assert_frame_equal( result.data, pl.DataFrame({VAR_KEY: [CONST_TERM], COEF_KEY: [2]}), - check_dtype=False, + check_dtypes=False, check_column_order=False, ) @@ -150,7 +150,7 @@ def test_add_expressions_with_vars(): assert_frame_equal( result.data, pl.DataFrame({VAR_KEY: [1, 2], COEF_KEY: [2, 4]}), - check_dtype=False, + check_dtypes=False, check_column_order=False, ) @@ -167,7 +167,7 @@ def test_add_expressions_with_vars_and_dims(): pl.DataFrame( {"dim1": [1, 1, 2, 2], VAR_KEY: [1, 2, 1, 2], COEF_KEY: [2, 4, 6, 8]} ), - check_dtype=False, + check_dtypes=False, check_column_order=False, ) @@ -226,7 +226,7 @@ def test_add_expression_with_vars_and_add_dim(): assert_frame_equal( result.data, expected_result, - check_dtype=False, + check_dtypes=False, check_column_order=False, check_row_order=False, ) @@ -236,7 +236,7 @@ def test_add_expression_with_vars_and_add_dim(): assert_frame_equal( result.data, expected_result, - check_dtype=False, + check_dtypes=False, check_column_order=False, check_row_order=False, ) @@ -391,7 +391,7 @@ def test_three_way_add(): pl.DataFrame( {"dim1": [1, 2], VAR_KEY: [CONST_TERM, CONST_TERM], COEF_KEY: [9, 4]} ), - check_dtype=False, + check_dtypes=False, check_column_order=False, ) @@ -402,7 +402,7 @@ def test_three_way_add(): assert_frame_equal( result.data, pl.DataFrame({"dim1": [1], VAR_KEY: [CONST_TERM], COEF_KEY: [9]}), - check_dtype=False, + check_dtypes=False, check_column_order=False, ) diff --git a/tests/test_solver.py b/tests/test_solver.py index 67fecfb..ba7c8e8 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -44,7 +44,7 @@ def test_retrieving_duals_vectorized(): assert_frame_equal( m.X.solution, pl.DataFrame({"t": [1, 2], "solution": [45, 10]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( @@ -52,19 +52,19 @@ def test_retrieving_duals_vectorized(): pl.DataFrame( {"t": [1, 2], "RC": [0, 0]} ), # Somehow the reduced cost is 0 since we are no longer using a bound. - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( m.max_AB.dual, pl.DataFrame({"c": [1, 2], "dual": [0.1, 0]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( m.max_AB.slack, pl.DataFrame({"c": [1, 2], "slack": [0, 50]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) @@ -87,25 +87,25 @@ def test_support_variable_attributes(): assert_frame_equal( m.X.solution, pl.DataFrame({"t": [1, 2], "solution": [45, 10]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( m.X.RC, pl.DataFrame({"t": [1, 2], "RC": [0.0, 1.9]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( m.max_AB.dual, pl.DataFrame({"c": [1, 2], "dual": [0.1, 0]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) assert_frame_equal( m.max_AB.slack, pl.DataFrame({"c": [1, 2], "slack": [0, 50]}), - check_dtype=False, + check_dtypes=False, check_row_order=False, ) From 5e7852854db8af32de9120748c14a3db09761dae Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:08:53 -0500 Subject: [PATCH 02/17] Add test for sets with range --- src/pyoframe/core.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index fcdcca9..564814b 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -195,6 +195,15 @@ def __eq__(self, value: object): class Set(ModelElement, SupportsMath, SupportPolarsMethodMixin): + """ + A set which can then be used to index variables. + + Examples: + >>> import pyoframe as pf + >>> pf.Set(x=range(2), y=range(3)) + + [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)] + """ def __init__(self, *data: SetTypes | Iterable[SetTypes], **named_data): data_list = list(data) for name, set in named_data.items(): From e0188c264edece0e957c07235f94e61f5cd4739c Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:15:22 -0500 Subject: [PATCH 03/17] Branch .join statement for different versions --- src/pyoframe/_arithmetic.py | 5 ++++- src/pyoframe/core.py | 20 ++++++++++++-------- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index 3491582..e1a754f 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -8,6 +8,7 @@ UnmatchedStrategy, Config, PyoframeError, + POLARS_VERSION, ) if TYPE_CHECKING: # pragma: no cover @@ -116,7 +117,9 @@ def get_indices(expr): not Config.disable_unmatched_checks ), "This code should not be reached when unmatched checks are disabled." outer_join = get_indices(left).join( - get_indices(right), how="full", on=dims + get_indices(right), + how="full" if POLARS_VERSION.major >= 1 else "outer", + on=dims, ) if outer_join.get_column(dims[0]).null_count() > 0: raise PyoframeError( diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 564814b..6868440 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -197,13 +197,14 @@ def __eq__(self, value: object): class Set(ModelElement, SupportsMath, SupportPolarsMethodMixin): """ A set which can then be used to index variables. - + Examples: >>> import pyoframe as pf >>> pf.Set(x=range(2), y=range(3)) [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)] """ + def __init__(self, *data: SetTypes | Iterable[SetTypes], **named_data): data_list = list(data) for name, set in named_data.items(): @@ -700,7 +701,10 @@ def _add_const(self, const: int | float) -> Expression: .unique(maintain_order=True) .with_columns(pl.lit(CONST_TERM).alias(VAR_KEY).cast(VAR_TYPE)) ) - data = data.join(keys, on=dim + [VAR_KEY], how="full", coalesce=True) + if POLARS_VERSION.major >= 1: + data = data.join(keys, on=dim + [VAR_KEY], how="full", coalesce=True) + else: + data = data.join(keys, on=dim + [VAR_KEY], how="outer_coalesce") data = data.with_columns(pl.col(COEF_KEY).fill_null(0.0)) data = data.with_columns( @@ -716,12 +720,12 @@ def constant_terms(self): dims = self.dimensions constant_terms = self.data.filter(pl.col(VAR_KEY) == CONST_TERM).drop(VAR_KEY) if dims is not None: - return constant_terms.join( - self.data.select(dims).unique(maintain_order=True), - on=dims, - how="full", - coalesce=True, - ).with_columns(pl.col(COEF_KEY).fill_null(0.0)) + dims_df = self.data.select(dims).unique(maintain_order=True) + if POLARS_VERSION.major >= 1: + df = constant_terms.join(dims_df, on=dims, how="full", coalesce=True) + else: + df = constant_terms.join(dims_df, on=dims, how="outer_coalesce") + return df.with_columns(pl.col(COEF_KEY).fill_null(0.0)) else: if len(constant_terms) == 0: return pl.DataFrame( From 33ee7ae0161ce5ccbca04dee6e6e3d3b46ff3874 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:19:46 -0500 Subject: [PATCH 04/17] Support previous polars versions --- tests/test_arithmetic.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/tests/test_arithmetic.py b/tests/test_arithmetic.py index e525ecf..b0aefe8 100644 --- a/tests/test_arithmetic.py +++ b/tests/test_arithmetic.py @@ -6,10 +6,11 @@ from polars.testing import assert_frame_equal import polars as pl -from pyoframe.constants import COEF_KEY, CONST_TERM, VAR_KEY +from pyoframe.constants import COEF_KEY, CONST_TERM, VAR_KEY, POLARS_VERSION from pyoframe import Variable, Model, sum, Set, Config, Expression, VType from .util import csvs_to_expr +check_dtypes_false = {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": True} def test_set_multiplication(): dim1 = [1, 2, 3] @@ -51,7 +52,7 @@ def test_multiplication_no_common_dimensions(): VAR_KEY: [CONST_TERM] * 6, } ), - check_dtypes=False, + **check_dtypes_false, ) @@ -71,7 +72,7 @@ def test_within_set(): VAR_KEY: [1, 4], } ), - check_dtypes=False, + **check_dtypes_false, ) @@ -82,7 +83,7 @@ def test_filter_expression(): assert_frame_equal( result.data, pl.DataFrame({"dim1": [2], COEF_KEY: [2], VAR_KEY: [CONST_TERM]}), - check_dtypes=False, + **check_dtypes_false, ) @@ -92,7 +93,7 @@ def test_filter_constraint(): assert_frame_equal( result, pl.DataFrame({"dim1": [2], COEF_KEY: [2], VAR_KEY: [CONST_TERM]}), - check_dtypes=False, + **check_dtypes_false, ) @@ -139,7 +140,7 @@ def test_add_expressions(): assert_frame_equal( result.data, pl.DataFrame({VAR_KEY: [CONST_TERM], COEF_KEY: [2]}), - check_dtypes=False, + **check_dtypes_false, check_column_order=False, ) @@ -150,7 +151,7 @@ def test_add_expressions_with_vars(): assert_frame_equal( result.data, pl.DataFrame({VAR_KEY: [1, 2], COEF_KEY: [2, 4]}), - check_dtypes=False, + **check_dtypes_false, check_column_order=False, ) @@ -167,7 +168,7 @@ def test_add_expressions_with_vars_and_dims(): pl.DataFrame( {"dim1": [1, 1, 2, 2], VAR_KEY: [1, 2, 1, 2], COEF_KEY: [2, 4, 6, 8]} ), - check_dtypes=False, + **check_dtypes_false, check_column_order=False, ) @@ -226,7 +227,7 @@ def test_add_expression_with_vars_and_add_dim(): assert_frame_equal( result.data, expected_result, - check_dtypes=False, + **check_dtypes_false, check_column_order=False, check_row_order=False, ) @@ -236,7 +237,7 @@ def test_add_expression_with_vars_and_add_dim(): assert_frame_equal( result.data, expected_result, - check_dtypes=False, + **check_dtypes_false, check_column_order=False, check_row_order=False, ) @@ -391,7 +392,7 @@ def test_three_way_add(): pl.DataFrame( {"dim1": [1, 2], VAR_KEY: [CONST_TERM, CONST_TERM], COEF_KEY: [9, 4]} ), - check_dtypes=False, + **check_dtypes_false, check_column_order=False, ) @@ -402,7 +403,7 @@ def test_three_way_add(): assert_frame_equal( result.data, pl.DataFrame({"dim1": [1], VAR_KEY: [CONST_TERM], COEF_KEY: [9]}), - check_dtypes=False, + **check_dtypes_false, check_column_order=False, ) From 8769b227f03342c145532fd8b7ad3cafec98cf0f Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:20:24 -0500 Subject: [PATCH 05/17] Fix lint --- tests/test_arithmetic.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_arithmetic.py b/tests/test_arithmetic.py index b0aefe8..de355cf 100644 --- a/tests/test_arithmetic.py +++ b/tests/test_arithmetic.py @@ -10,7 +10,10 @@ from pyoframe import Variable, Model, sum, Set, Config, Expression, VType from .util import csvs_to_expr -check_dtypes_false = {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": True} +check_dtypes_false = ( + {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": True} +) + def test_set_multiplication(): dim1 = [1, 2, 3] From 7412a634fcaba2edec8f8f3464731b761fdaecaa Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:23:38 -0500 Subject: [PATCH 06/17] Fix tests --- tests/test_arithmetic.py | 4 ++-- tests/test_solver.py | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/tests/test_arithmetic.py b/tests/test_arithmetic.py index de355cf..8ad8689 100644 --- a/tests/test_arithmetic.py +++ b/tests/test_arithmetic.py @@ -11,7 +11,7 @@ from .util import csvs_to_expr check_dtypes_false = ( - {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": True} + {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": False} ) @@ -110,7 +110,7 @@ def test_filter_variable(): def test_filter_set(): s = Set(x=[1, 2, 3]) result = s.filter(x=2) - assert_frame_equal(result.data, pl.DataFrame({"x": [2]}), check_dtypes=False) + assert_frame_equal(result.data, pl.DataFrame({"x": [2]}), **check_dtypes_false) def test_drops_na(): diff --git a/tests/test_solver.py b/tests/test_solver.py index ba7c8e8..b572158 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -1,8 +1,13 @@ import pyoframe as pf +from pyoframe.constants import POLARS_VERSION import polars as pl from polars.testing import assert_frame_equal import pytest +check_dtypes_false = ( + {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": False} +) + def test_retrieving_duals(): m = pf.Model("max") @@ -44,28 +49,28 @@ def test_retrieving_duals_vectorized(): assert_frame_equal( m.X.solution, pl.DataFrame({"t": [1, 2], "solution": [45, 10]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) assert_frame_equal( m.X.RC, pl.DataFrame( {"t": [1, 2], "RC": [0, 0]} ), # Somehow the reduced cost is 0 since we are no longer using a bound. - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) assert_frame_equal( m.max_AB.dual, pl.DataFrame({"c": [1, 2], "dual": [0.1, 0]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) assert_frame_equal( m.max_AB.slack, pl.DataFrame({"c": [1, 2], "slack": [0, 50]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) @@ -87,26 +92,26 @@ def test_support_variable_attributes(): assert_frame_equal( m.X.solution, pl.DataFrame({"t": [1, 2], "solution": [45, 10]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) assert_frame_equal( m.X.RC, pl.DataFrame({"t": [1, 2], "RC": [0.0, 1.9]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) assert_frame_equal( m.max_AB.dual, pl.DataFrame({"c": [1, 2], "dual": [0.1, 0]}), - check_dtypes=False, + **check_dtypes_false, check_row_order=False, ) assert_frame_equal( m.max_AB.slack, pl.DataFrame({"c": [1, 2], "slack": [0, 50]}), - check_dtypes=False, check_row_order=False, + **check_dtypes_false, ) From d6b2b4bf1ab6316d1f9562b4a6aaf5bb72d69281 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Fri, 8 Nov 2024 13:31:12 -0500 Subject: [PATCH 07/17] Fix tests (try 2) --- tests/test_arithmetic.py | 2 +- tests/test_solver.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_arithmetic.py b/tests/test_arithmetic.py index 8ad8689..9f941f6 100644 --- a/tests/test_arithmetic.py +++ b/tests/test_arithmetic.py @@ -11,7 +11,7 @@ from .util import csvs_to_expr check_dtypes_false = ( - {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": False} + {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtype": False} ) diff --git a/tests/test_solver.py b/tests/test_solver.py index b572158..3be1ee3 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -5,7 +5,7 @@ import pytest check_dtypes_false = ( - {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtypes": False} + {"check_dtypes": False} if POLARS_VERSION.major >= 1 else {"check_dtype": False} ) From e9ae0f7d3bc1ccb51076a8a68a33532218f298fe Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Sun, 10 Nov 2024 18:18:58 -0500 Subject: [PATCH 08/17] Add support for quadratics --- docs/quadratic_expressions.md | 16 +++ src/pyoframe/_arithmetic.py | 146 ++++++++++++++++++++++- src/pyoframe/constants.py | 16 ++- src/pyoframe/core.py | 219 +++++++++++++++++++++++----------- src/pyoframe/io_mappers.py | 17 ++- src/pyoframe/model_element.py | 9 +- 6 files changed, 342 insertions(+), 81 deletions(-) create mode 100644 docs/quadratic_expressions.md diff --git a/docs/quadratic_expressions.md b/docs/quadratic_expressions.md new file mode 100644 index 0000000..27c316c --- /dev/null +++ b/docs/quadratic_expressions.md @@ -0,0 +1,16 @@ +# Quadratic Expressions + +Quadratic expressions work as you'd expect. Simply multiply two linear expression together (or square an expression with `**2`) and you'll get a quadratic. The quadratic can then be used in constraints or the objective. + +## Example + +## TODO + +```python3 +``` + +## Note for developers: Internal Representation of Quadratics + +Internally, Pyoframe's `Expression` object is used for both linear and quadratic expressions. When the dataframe within an `Expression` object (i.e. `Expression.data`) contains an additional column (named `__quadratic_variable_id`) we know that the expression is a quadratic. + +This extra column stores the ID of the second variable in quadratic terms. For terms with only one variable, this column contains ID `0` (a reserved variable ID which can thought of as equalling `1`). The variable ID in the `__variable_id` column is always greater or equal than the variable ID in the `__quadratic_variable_id`. This means that linear terms always have the variable id in the first column and `0` in the second column. Also, a `0` in the first column implies that the second column must also be `0` and therefore the term is a constant. \ No newline at end of file diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index e1a754f..729b30e 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -3,6 +3,9 @@ from pyoframe.constants import ( COEF_KEY, + CONST_TERM, + QUAD_VAR_KEY, + VAR_TYPE, RESERVED_COL_KEYS, VAR_KEY, UnmatchedStrategy, @@ -15,6 +18,20 @@ from pyoframe.core import Expression +def _multiply_expressions(*expressions: "Expression") -> "Expression": + try: + return _multiply_expressions_core(*expressions) + except PyoframeError as error: + raise PyoframeError( + "Failed to multiply expressions:\n" + + " * ".join( + e.to_str(include_header=True, include_data=False) for e in expressions + ) + + "\nDue to error:\n" + + str(error) + ) from error + + def _add_expressions(*expressions: "Expression") -> "Expression": try: return _add_expressions_core(*expressions) @@ -29,6 +46,106 @@ def _add_expressions(*expressions: "Expression") -> "Expression": ) from error +def _multiply_expressions_core(*expressions: "Expression") -> "Expression": + if len(expressions) > 2: + result = expressions[0] + for expr in expressions[1:]: + result = _multiply_expressions(result, expr) + return result + + self, other = expressions + + self_degree, other_degree = self.degree(), other.degree() + if self_degree + other_degree > 2: + # We know one of the two must be a quadratic since 1 + 1 is not greater than 2. + raise PyoframeError("Cannot multiply a quadratic expression by a non-constant.") + if self_degree < other_degree: + self, other = other, self + self_degree, other_degree = other_degree, self_degree + if other_degree == 1: + assert ( + self_degree == 1 + ), "This should always be true since the sum of degrees must be <=2." + return _quadratic_multiplication(self, other) + + assert ( + other_degree == 0 + ), "This should always be true since other cases have already been handled." + multiplier = other.data.drop( + VAR_KEY + ) # QUAD_VAR_KEY doesn't need to be dropped since we know it doesn't exist + + dims = self.dimensions_unsafe + other_dims = other.dimensions_unsafe + dims_in_common = [dim for dim in dims if dim in other_dims] + + data = ( + self.data.join( + multiplier, + on=dims_in_common if len(dims_in_common) > 0 else None, + how="inner" if dims_in_common else "cross", + ) + .with_columns(pl.col(COEF_KEY) * pl.col(COEF_KEY + "_right")) + .drop(COEF_KEY + "_right") + ) + + return self._new(data) + + +def _quadratic_multiplication(self: "Expression", other: "Expression") -> "Expression": + """ + Multiply two expressions of degree 1. + + Examples: + >>> import polars as pl + >>> from pyoframe import Variable + >>> df = pl.DataFrame({"dim": [1, 2, 3], "value": [1, 2, 3]}) + >>> x1 = Variable() + >>> x2 = Variable() + >>> expr1 = df * x1 + >>> expr2 = df * x2 * 2 + 4 + >>> expr1 * expr2 + + [1]: 4 x1 +2 x2*x1 + [2]: 8 x1 +8 x2*x1 + [3]: 12 x1 +18 x2*x1 + >>> (expr1 * expr2) - df * x1 * df * x2 * 2 + + [1]: 4 x1 + [2]: 8 x1 + [3]: 12 x1 + """ + dims = self.dimensions_unsafe + other_dims = other.dimensions_unsafe + dims_in_common = [dim for dim in dims if dim in other_dims] + + data = ( + self.data.join( + other.data, + on=dims_in_common if len(dims_in_common) > 0 else None, + how="inner" if dims_in_common else "cross", + ) + .with_columns(pl.col(COEF_KEY) * pl.col(COEF_KEY + "_right")) + .drop(COEF_KEY + "_right") + .rename({VAR_KEY + "_right": QUAD_VAR_KEY}) + # Swap VAR_KEY and QUAD_VAR_KEY so that VAR_KEy is always the larger one + .with_columns( + pl.when(pl.col(VAR_KEY) < pl.col(QUAD_VAR_KEY)) + .then(pl.col(QUAD_VAR_KEY)) + .otherwise(pl.col(VAR_KEY)) + .alias(VAR_KEY), + pl.when(pl.col(VAR_KEY) < pl.col(QUAD_VAR_KEY)) + .then(pl.col(VAR_KEY)) + .otherwise(pl.col(QUAD_VAR_KEY)) + .alias(QUAD_VAR_KEY), + ) + ) + + data = _sum_like_terms(data) + + return self._new(data) + + def _add_expressions_core(*expressions: "Expression") -> "Expression": # Mapping of how a sum of two expressions should propogate the unmatched strategy propogatation_strategies = { @@ -162,11 +279,22 @@ def get_indices(expr): propogate_strat = expressions[0].unmatched_strategy expr_data = [expr.data for expr in expressions] + # Add quadratic column if it is needed and doesn't already exist + if any(QUAD_VAR_KEY in df.columns for df in expr_data): + expr_data = [ + ( + df.with_columns(pl.lit(CONST_TERM).alias(QUAD_VAR_KEY).cast(VAR_TYPE)) + if QUAD_VAR_KEY not in df.columns + else df + ) + for df in expr_data + ] + # Sort columns to allow for concat - expr_data = [e.select(sorted(e.columns)) for e in expr_data] + expr_data = [e.select(dims + [c for c in e.columns if c not in dims]) for e in expr_data] data = pl.concat(expr_data, how="vertical_relaxed") - data = data.group_by(dims + [VAR_KEY], maintain_order=True).sum() + data = _sum_like_terms(data) new_expr = expressions[0]._new(data) new_expr.unmatched_strategy = propogate_strat @@ -214,6 +342,20 @@ def _add_dimension(self: "Expression", target: "Expression") -> "Expression": return self._new(result) +def _sum_like_terms(df: pl.DataFrame) -> pl.DataFrame: + """Combines terms with the same variables. Removes quadratic column if they all happen to cancel.""" + dims = [c for c in df.columns if c not in RESERVED_COL_KEYS] + var_cols = [VAR_KEY] + ([QUAD_VAR_KEY] if QUAD_VAR_KEY in df.columns else []) + df = ( + df.group_by(dims + var_cols, maintain_order=True) + .sum() + .filter(pl.col(COEF_KEY) != 0) + ) + if QUAD_VAR_KEY in df.columns and (df.get_column(QUAD_VAR_KEY) == CONST_TERM).all(): + df = df.drop(QUAD_VAR_KEY) + return df + + def _get_dimensions(df: pl.DataFrame) -> Optional[List[str]]: """ Returns the dimensions of the DataFrame. Reserved columns do not count as dimensions. diff --git a/src/pyoframe/constants.py b/src/pyoframe/constants.py index 949ef68..6b0e3b7 100644 --- a/src/pyoframe/constants.py +++ b/src/pyoframe/constants.py @@ -20,17 +20,31 @@ COEF_KEY = "__coeff" VAR_KEY = "__variable_id" +QUAD_VAR_KEY = "__quadratic_variable_id" CONSTRAINT_KEY = "__constraint_id" SOLUTION_KEY = "solution" DUAL_KEY = "dual" RC_COL = "RC" SLACK_COL = "slack" -CONST_TERM = 0 +COL_DTYPES = { + COEF_KEY: pl.Float64, + VAR_KEY: pl.UInt32, + QUAD_VAR_KEY: pl.UInt32, + CONSTRAINT_KEY: pl.UInt32, + SOLUTION_KEY: pl.Float64, + DUAL_KEY: pl.Float64, + RC_COL: pl.Float64, + SLACK_COL: pl.Float64, +} +VAR_TYPE = COL_DTYPES[VAR_KEY] + +CONST_TERM = 0 # 0 is a reserved value which makes it easy to detect. It must be zero for e.g. ensuring consistency during quadratic multiplication. RESERVED_COL_KEYS = ( COEF_KEY, VAR_KEY, + QUAD_VAR_KEY, CONSTRAINT_KEY, SOLUTION_KEY, DUAL_KEY, diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 6868440..9516624 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -15,9 +15,12 @@ import pandas as pd import polars as pl -from packaging import version -from pyoframe._arithmetic import _add_expressions, _get_dimensions +from pyoframe._arithmetic import ( + _add_expressions, + _get_dimensions, + _multiply_expressions, +) from pyoframe.constants import ( COEF_KEY, CONST_TERM, @@ -29,6 +32,8 @@ SLACK_COL, SOLUTION_KEY, VAR_KEY, + QUAD_VAR_KEY, + VAR_TYPE, Config, ConstraintSense, ObjSense, @@ -55,8 +60,6 @@ if TYPE_CHECKING: # pragma: no cover from pyoframe.model import Model -VAR_TYPE = pl.UInt32 - def _forward_to_expression(func_name: str): def wrapper(self: "SupportsMath", *args, **kwargs) -> "Expression": @@ -98,6 +101,18 @@ def to_expr(self) -> "Expression": ... sum = _forward_to_expression("sum") map = _forward_to_expression("map") + def __pow__(self, power: int): + """ + Support squaring expressions: + >>> from pyoframe import Variable + >>> Variable() ** 2 + + x1*x1 + """ + if power == 2: + return self * self + return NotImplemented + def __neg__(self): res = self.to_expr() * -1 # Negating a constant term should keep the unmatched strategy @@ -353,7 +368,7 @@ def _set_to_polars(set: "SetTypes") -> pl.DataFrame: class Expression(ModelElement, SupportsMath, SupportPolarsMethodMixin): - """A linear expression.""" + """A linear or quadratic expression.""" def __init__(self, data: pl.DataFrame): """ @@ -433,10 +448,17 @@ def sum(self, over: Union[str, Iterable[str]]): return self._new( self.data.drop(over) - .group_by(remaining_dims + [VAR_KEY], maintain_order=True) + .group_by(remaining_dims + self._variable_columns, maintain_order=True) .sum() ) + @property + def _variable_columns(self) -> List[str]: + if self.is_quadratic: + return [VAR_KEY, QUAD_VAR_KEY] + else: + return [VAR_KEY] + def map(self, mapping_set: SetTypes, drop_shared_dims: bool = True): """ Replaces the dimensions that are shared with mapping_set with the other dimensions found in mapping_set. @@ -585,6 +607,48 @@ def within(self, set: "SetTypes") -> Expression: by_dims = df.select(dims_in_common).unique(maintain_order=True) return self._new(self.data.join(by_dims, on=dims_in_common)) + @property + def is_quadratic(self) -> bool: + """ + Returns True if the expression is quadratic, False otherwise. + + Computes in O(1) since expressions are quadratic if and + only if self.data contain the QUAD_VAR_KEY column. + + Examples: + >>> import pandas as pd + >>> from pyoframe import Variable + >>> expr = pd.DataFrame({"dim1": [1, 2, 3], "value": [1, 2, 3]}) * Variable() + >>> expr *= Variable() + >>> expr.is_quadratic + True + """ + return QUAD_VAR_KEY in self.data.columns + + def degree(self) -> int: + """ + Returns the degree of the expression (0=constant, 1=linear, 2=quadratic). + + Examples: + >>> import pandas as pd + >>> from pyoframe import Variable + >>> expr = pd.DataFrame({"dim1": [1, 2, 3], "value": [1, 2, 3]}).to_expr() + >>> expr.degree() + 0 + >>> expr *= Variable() + >>> expr.degree() + 1 + >>> expr += (Variable() ** 2).add_dim("dim1") + >>> expr.degree() + 2 + """ + if self.is_quadratic: + return 2 + elif (self.data.get_column(VAR_KEY) != CONST_TERM).any(): + return 1 + else: + return 0 + def __add__(self, other): """ Examples: @@ -639,31 +703,7 @@ def __mul__( other = other.to_expr() self._learn_from_other(other) - - if (other.data.get_column(VAR_KEY) != CONST_TERM).any(): - self, other = other, self - - if (other.data.get_column(VAR_KEY) != CONST_TERM).any(): - raise ValueError( - "Multiplication of two expressions with variables is non-linear and not supported." - ) - multiplier = other.data.drop(VAR_KEY) - - dims = self.dimensions_unsafe - other_dims = other.dimensions_unsafe - dims_in_common = [dim for dim in dims if dim in other_dims] - - data = ( - self.data.join( - multiplier, - on=dims_in_common if len(dims_in_common) > 0 else None, - how="inner" if dims_in_common else "cross", - ) - .with_columns(pl.col(COEF_KEY) * pl.col(COEF_KEY + "_right")) - .drop(COEF_KEY + "_right") - ) - - return self._new(data) + return _multiply_expressions(self, other) def to_expr(self) -> Expression: return self @@ -680,19 +720,32 @@ def _new(self, data: pl.DataFrame) -> Expression: return e def _add_const(self, const: int | float) -> Expression: + """ + >>> Variable() + 5 + + x1 +5 + >>> Variable() ** 2 + 5 + + x2*x2 +5 + >>> Variable() ** 2 + Variable() + 5 + + x3*x3 + x4 +5 + """ dim = self.dimensions data = self.data # Fill in missing constant terms if not dim: if CONST_TERM not in data[VAR_KEY]: + const_df = pl.DataFrame( + {COEF_KEY: [0.0], VAR_KEY: [CONST_TERM]}, + schema={COEF_KEY: pl.Float64, VAR_KEY: VAR_TYPE}, + ) + if self.is_quadratic: + const_df = const_df.with_columns( + pl.lit(CONST_TERM).alias(QUAD_VAR_KEY).cast(VAR_TYPE) + ) data = pl.concat( - [ - data, - pl.DataFrame( - {COEF_KEY: [0.0], VAR_KEY: [CONST_TERM]}, - schema={COEF_KEY: pl.Float64, VAR_KEY: VAR_TYPE}, - ), - ], + [data, const_df], how="vertical_relaxed", ) else: @@ -701,10 +754,18 @@ def _add_const(self, const: int | float) -> Expression: .unique(maintain_order=True) .with_columns(pl.lit(CONST_TERM).alias(VAR_KEY).cast(VAR_TYPE)) ) + if self.is_quadratic: + keys = keys.with_columns( + pl.lit(CONST_TERM).alias(QUAD_VAR_KEY).cast(VAR_TYPE) + ) if POLARS_VERSION.major >= 1: - data = data.join(keys, on=dim + [VAR_KEY], how="full", coalesce=True) + data = data.join( + keys, on=dim + self._variable_columns, how="full", coalesce=True + ) else: - data = data.join(keys, on=dim + [VAR_KEY], how="outer_coalesce") + data = data.join( + keys, on=dim + self._variable_columns, how="outer_coalesce" + ) data = data.with_columns(pl.col(COEF_KEY).fill_null(0.0)) data = data.with_columns( @@ -719,6 +780,8 @@ def _add_const(self, const: int | float) -> Expression: def constant_terms(self): dims = self.dimensions constant_terms = self.data.filter(pl.col(VAR_KEY) == CONST_TERM).drop(VAR_KEY) + if self.is_quadratic: + constant_terms = constant_terms.drop(QUAD_VAR_KEY) if dims is not None: dims_df = self.data.select(dims).unique(maintain_order=True) if POLARS_VERSION.major >= 1: @@ -776,19 +839,23 @@ def value(self) -> pl.DataFrame: "Can't obtain value of expression since the model has not been solved." ) - df = ( - self.data.join(self._model.result.solution.primal, on=VAR_KEY, how="left") - .with_columns( - ( - pl.when(pl.col(VAR_KEY) == CONST_TERM) - .then(1) - .otherwise(pl.col(SOLUTION_KEY)) - * pl.col(COEF_KEY) - ).alias(SOLUTION_KEY) + df = self.data + for var_col in self._variable_columns: + df = ( + df.join(self._model.result.solution.primal, on=var_col, how="left") + .with_columns( + ( + pl.when(pl.col(var_col) == CONST_TERM) + .then(1) + .otherwise(pl.col(SOLUTION_KEY)) + * pl.col(COEF_KEY) + ).alias(COEF_KEY) + ) + .drop(var_col) + .drop(SOLUTION_KEY) ) - .drop(VAR_KEY) - .drop(COEF_KEY) - ) + + df = df.rename({COEF_KEY: SOLUTION_KEY}) dims = self.dimensions if dims is not None: @@ -807,24 +874,34 @@ def to_str_table( data = self.data if include_const_term else self.variable_terms data = cast_coef_to_string(data, float_precision=float_precision) - if var_map is not None: - data = var_map.apply(data, to_col="str_var") - elif self._model is not None and self._model.var_map is not None: - var_map = self._model.var_map - data = var_map.apply(data, to_col="str_var") - else: - data = data.with_columns( - pl.concat_str(pl.lit("x"), VAR_KEY).alias("str_var") - ) - if include_const_variable: - data = data.drop(VAR_KEY).rename({"str_var": VAR_KEY}) - else: + for var_column in self._variable_columns: + temp_var_column = f"{var_column}_temp" + if var_map is not None: + data = var_map.apply(data, to_col=temp_var_column, id_col=var_column) + elif self._model is not None and self._model.var_map is not None: + var_map = self._model.var_map + data = var_map.apply(data, to_col=temp_var_column, id_col=var_column) + else: + data = data.with_columns( + pl.concat_str(pl.lit("x"), var_column).alias(temp_var_column) + ) + if include_const_variable and var_column == VAR_KEY: + data = data.drop(var_column).rename({temp_var_column: var_column}) + else: + data = data.with_columns( + pl.when(pl.col(var_column) == CONST_TERM) + .then(pl.lit("")) + .otherwise(temp_var_column) + .alias(var_column) + ).drop(temp_var_column) + + if self.is_quadratic: data = data.with_columns( - pl.when(pl.col(VAR_KEY) == CONST_TERM) - .then(pl.lit("")) - .otherwise("str_var") + pl.when(pl.col(QUAD_VAR_KEY) == "") + .then(pl.col(VAR_KEY)) + .otherwise(pl.concat_str(VAR_KEY, pl.lit("*"), pl.col(QUAD_VAR_KEY))) .alias(VAR_KEY) - ).drop("str_var") + ).drop(QUAD_VAR_KEY) dimensions = self.dimensions @@ -888,7 +965,11 @@ def to_str( result = "" if include_header: result += get_obj_repr( - self, size=len(self), dimensions=self.shape, terms=len(self.data) + self, + size=len(self), + dimensions=self.shape, + terms=len(self.data), + degree=2 if self.degree() == 2 else None, ) if include_header and include_data: result += "\n" diff --git a/src/pyoframe/io_mappers.py b/src/pyoframe/io_mappers.py index 0c3ad83..17e793f 100644 --- a/src/pyoframe/io_mappers.py +++ b/src/pyoframe/io_mappers.py @@ -49,15 +49,21 @@ def apply( self, df: pl.DataFrame, to_col: Optional[str] = None, + id_col: Optional[str] = None, ) -> pl.DataFrame: if df.height == 0: return df + if id_col is None: + id_col = self._ID_COL result = df.join( - self.mapping_registry, on=self._ID_COL, how="left", validate="m:1" + self.mapping_registry, how="left", validate="m:1", left_on=id_col, right_on=self._ID_COL ) - if to_col is None: + if id_col != self._ID_COL: result = result.drop(self._ID_COL) - to_col = self._ID_COL + if to_col is None: + # Drop self._ID_COL so we can replace it with the result + result = result.drop(id_col) + to_col = id_col return result.rename({Mapper.NAME_COL: to_col}) def undo(self, df: pl.DataFrame) -> pl.DataFrame: @@ -125,13 +131,16 @@ def apply( self, df: pl.DataFrame, to_col: Optional[str] = None, + id_col: Optional[str] = None, ) -> pl.DataFrame: if df.height == 0: return df + if id_col is None: + id_col = self._ID_COL query = pl.concat_str( pl.lit(self._prefix), - pl.col(self._ID_COL).map_batches( + pl.col(id_col).map_batches( Base36Mapper._to_base36, return_dtype=pl.String, is_elementwise=True, diff --git a/src/pyoframe/model_element.py b/src/pyoframe/model_element.py index fb4f0e4..73190e1 100644 --- a/src/pyoframe/model_element.py +++ b/src/pyoframe/model_element.py @@ -5,7 +5,7 @@ import polars as pl from typing import TYPE_CHECKING -from pyoframe.constants import COEF_KEY, RESERVED_COL_KEYS, VAR_KEY +from pyoframe.constants import COEF_KEY, RESERVED_COL_KEYS, VAR_KEY, QUAD_VAR_KEY, COL_DTYPES from pyoframe._arithmetic import _get_dimensions from pyoframe.user_defined import AttrContainerMixin @@ -29,10 +29,9 @@ def __init__(self, data: pl.DataFrame, **kwargs) -> None: data = data.select(cols) # Cast to proper dtype - if COEF_KEY in data.columns: - data = data.cast({COEF_KEY: pl.Float64}) - if VAR_KEY in data.columns: - data = data.cast({VAR_KEY: pl.UInt32}) + for c in [COEF_KEY, VAR_KEY, QUAD_VAR_KEY]: + if c in data.columns: + data = data.cast({c: COL_DTYPES[c]}) self._data = data self._model: Optional[Model] = None From 9eef395b8de99d10b0a9f29c725b537209df1e82 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Sun, 10 Nov 2024 18:20:09 -0500 Subject: [PATCH 09/17] Apply linter --- src/pyoframe/_arithmetic.py | 4 +++- src/pyoframe/io_mappers.py | 6 +++++- src/pyoframe/model_element.py | 8 +++++++- 3 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index 729b30e..a6eefb8 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -291,7 +291,9 @@ def get_indices(expr): ] # Sort columns to allow for concat - expr_data = [e.select(dims + [c for c in e.columns if c not in dims]) for e in expr_data] + expr_data = [ + e.select(dims + [c for c in e.columns if c not in dims]) for e in expr_data + ] data = pl.concat(expr_data, how="vertical_relaxed") data = _sum_like_terms(data) diff --git a/src/pyoframe/io_mappers.py b/src/pyoframe/io_mappers.py index 17e793f..9498c25 100644 --- a/src/pyoframe/io_mappers.py +++ b/src/pyoframe/io_mappers.py @@ -56,7 +56,11 @@ def apply( if id_col is None: id_col = self._ID_COL result = df.join( - self.mapping_registry, how="left", validate="m:1", left_on=id_col, right_on=self._ID_COL + self.mapping_registry, + how="left", + validate="m:1", + left_on=id_col, + right_on=self._ID_COL, ) if id_col != self._ID_COL: result = result.drop(self._ID_COL) diff --git a/src/pyoframe/model_element.py b/src/pyoframe/model_element.py index 73190e1..a5539fb 100644 --- a/src/pyoframe/model_element.py +++ b/src/pyoframe/model_element.py @@ -5,7 +5,13 @@ import polars as pl from typing import TYPE_CHECKING -from pyoframe.constants import COEF_KEY, RESERVED_COL_KEYS, VAR_KEY, QUAD_VAR_KEY, COL_DTYPES +from pyoframe.constants import ( + COEF_KEY, + RESERVED_COL_KEYS, + VAR_KEY, + QUAD_VAR_KEY, + COL_DTYPES, +) from pyoframe._arithmetic import _get_dimensions from pyoframe.user_defined import AttrContainerMixin From 062006ce0467b80d45117b05ba90b5ccf63310e1 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Sun, 10 Nov 2024 18:28:11 -0500 Subject: [PATCH 10/17] Update documentation --- docs/quadratic_expressions.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/quadratic_expressions.md b/docs/quadratic_expressions.md index 27c316c..6e866fb 100644 --- a/docs/quadratic_expressions.md +++ b/docs/quadratic_expressions.md @@ -9,8 +9,10 @@ Quadratic expressions work as you'd expect. Simply multiply two linear expressio ```python3 ``` -## Note for developers: Internal Representation of Quadratics +## Note for Pyoframe developers: Internal Representation of Quadratics Internally, Pyoframe's `Expression` object is used for both linear and quadratic expressions. When the dataframe within an `Expression` object (i.e. `Expression.data`) contains an additional column (named `__quadratic_variable_id`) we know that the expression is a quadratic. -This extra column stores the ID of the second variable in quadratic terms. For terms with only one variable, this column contains ID `0` (a reserved variable ID which can thought of as equalling `1`). The variable ID in the `__variable_id` column is always greater or equal than the variable ID in the `__quadratic_variable_id`. This means that linear terms always have the variable id in the first column and `0` in the second column. Also, a `0` in the first column implies that the second column must also be `0` and therefore the term is a constant. \ No newline at end of file +This extra column stores the ID of the second variable in quadratic terms. For terms with only one variable, this column contains ID `0` (a reserved variable ID which can thought of as meaning 'no variable'). The variables in a quadratic are rearranged such that the ID in the `__variable_id` column is always greater or equal than the variable ID in the `__quadratic_variable_id` (recall: a*b=b*a). This rearranging not only ensures that a*b+b*a=2a*b but also generates a useful property: If the variable ID in the first column (`__variable_id`) is `0` we know the variable ID in the second must also be `0` and therefore the term must be a constant. + +The additional quadratic variable ID column is automatically dropped if through arithmetic the quadratic terms cancel out. \ No newline at end of file From 251b7244619bd13e3c6f21add21ab59c72d8f570 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Sun, 10 Nov 2024 19:33:44 -0500 Subject: [PATCH 11/17] Start creating quadratic example --- src/pyoframe/io_mappers.py | 5 +- tests/examples/facility_location/model.py | 55 ++++++++ .../results/pyoframe-problem.lp | 124 ++++++++++++++++++ .../facility_location/results/solution.sol | 93 +++++++++++++ 4 files changed, 274 insertions(+), 3 deletions(-) create mode 100644 tests/examples/facility_location/model.py create mode 100644 tests/examples/facility_location/results/pyoframe-problem.lp create mode 100644 tests/examples/facility_location/results/solution.sol diff --git a/src/pyoframe/io_mappers.py b/src/pyoframe/io_mappers.py index 9498c25..4263f14 100644 --- a/src/pyoframe/io_mappers.py +++ b/src/pyoframe/io_mappers.py @@ -62,8 +62,6 @@ def apply( left_on=id_col, right_on=self._ID_COL, ) - if id_col != self._ID_COL: - result = result.drop(self._ID_COL) if to_col is None: # Drop self._ID_COL so we can replace it with the result result = result.drop(id_col) @@ -74,8 +72,9 @@ def undo(self, df: pl.DataFrame) -> pl.DataFrame: if df.height == 0: return df df = df.rename({self._ID_COL: Mapper.NAME_COL}) + # We use an inner join since sometimes Gurobi returns additional variables (i.e. quadratic variables) return df.join( - self.mapping_registry, on=Mapper.NAME_COL, how="left", validate="m:1" + self.mapping_registry, on=Mapper.NAME_COL, how="inner", validate="m:1" ).drop(Mapper.NAME_COL) diff --git a/tests/examples/facility_location/model.py b/tests/examples/facility_location/model.py new file mode 100644 index 0000000..9bb969c --- /dev/null +++ b/tests/examples/facility_location/model.py @@ -0,0 +1,55 @@ +import pyoframe as pf +import polars as pl +import os +from pathlib import Path + + +def main(F, G, **kwargs): + model = pf.Model("min") + + model.facilities = pf.Set(f=range(F)) + model.customers = pf.Set(c=range(G)) + + model.facility_position_x = pf.Variable(model.facilities, lb=0, ub=1) + model.facility_position_y = pf.Variable(model.facilities, lb=0, ub=1) + model.customer_position_x = pl.DataFrame( + {"c": range(G), "x": [step / G for step in range(G)]} + ).to_expr() + model.customer_position_y = pl.DataFrame( + {"c": range(G), "y": [step / G for step in range(G)]} + ).to_expr() + + model.max_distance = pf.Variable(lb=0) + model.objective = model.max_distance + + model.is_closest = pf.Variable(model.customers, model.facilities, vtype="binary") + model.con_only_one_closest = pf.sum("f", model.is_closest) == 1 + + model.dist_x = pf.Variable(model.customers, model.facilities) + model.dist_y = pf.Variable(model.customers, model.facilities) + model.con_dist_x = model.dist_x == model.customer_position_x.add_dim( + "f" + ) - model.facility_position_x.add_dim("c") + model.con_dist_y = model.dist_y == model.customer_position_y.add_dim( + "f" + ) - model.facility_position_y.add_dim("c") + model.dist = pf.Variable(model.customers, model.facilities, lb=0) + model.con_dist = model.dist**2 >= model.dist_x**2 + model.dist_y**2 + + M = 2 * 1.414 # Twice the max distance which ensures that when is_closest is 0, the constraint is not binding. + model.con_max_distance = model.max_distance.add_dim("c", "f") == model.dist - M * ( + 1 - model.is_closest + ) + + model.solve(**kwargs) + + +if __name__ == "__main__": + working_dir = Path(os.path.dirname(os.path.realpath(__file__))) + main( + 3, + 4, + directory=working_dir / "results", + use_var_names=True, + solution_file=working_dir / "results" / "solution.sol", + ) diff --git a/tests/examples/facility_location/results/pyoframe-problem.lp b/tests/examples/facility_location/results/pyoframe-problem.lp new file mode 100644 index 0000000..eaa10b5 --- /dev/null +++ b/tests/examples/facility_location/results/pyoframe-problem.lp @@ -0,0 +1,124 @@ +minimize + +obj: + +max_distance + +s.t. + +con_only_one_closest[0]: is_closest[0,0] + is_closest[0,1] + is_closest[0,2] = 1 +con_only_one_closest[1]: is_closest[1,0] + is_closest[1,1] + is_closest[1,2] = 1 +con_only_one_closest[2]: is_closest[2,0] + is_closest[2,1] + is_closest[2,2] = 1 +con_only_one_closest[3]: is_closest[3,0] + is_closest[3,1] + is_closest[3,2] = 1 +con_dist_x[0,0]: dist_x[0,0] + facility_position_x[0] = 0 +con_dist_x[0,1]: dist_x[0,1] + facility_position_x[1] = 0 +con_dist_x[0,2]: dist_x[0,2] + facility_position_x[2] = 0 +con_dist_x[1,0]: dist_x[1,0] + facility_position_x[0] = 0.25 +con_dist_x[1,1]: dist_x[1,1] + facility_position_x[1] = 0.25 +con_dist_x[1,2]: dist_x[1,2] + facility_position_x[2] = 0.25 +con_dist_x[2,0]: dist_x[2,0] + facility_position_x[0] = 0.5 +con_dist_x[2,1]: dist_x[2,1] + facility_position_x[1] = 0.5 +con_dist_x[2,2]: dist_x[2,2] + facility_position_x[2] = 0.5 +con_dist_x[3,0]: dist_x[3,0] + facility_position_x[0] = 0.75 +con_dist_x[3,1]: dist_x[3,1] + facility_position_x[1] = 0.75 +con_dist_x[3,2]: dist_x[3,2] + facility_position_x[2] = 0.75 +con_dist_y[0,0]: dist_y[0,0] + facility_position_y[0] = 0 +con_dist_y[0,1]: dist_y[0,1] + facility_position_y[1] = 0 +con_dist_y[0,2]: dist_y[0,2] + facility_position_y[2] = 0 +con_dist_y[1,0]: dist_y[1,0] + facility_position_y[0] = 0.25 +con_dist_y[1,1]: dist_y[1,1] + facility_position_y[1] = 0.25 +con_dist_y[1,2]: dist_y[1,2] + facility_position_y[2] = 0.25 +con_dist_y[2,0]: dist_y[2,0] + facility_position_y[0] = 0.5 +con_dist_y[2,1]: dist_y[2,1] + facility_position_y[1] = 0.5 +con_dist_y[2,2]: dist_y[2,2] + facility_position_y[2] = 0.5 +con_dist_y[3,0]: dist_y[3,0] + facility_position_y[0] = 0.75 +con_dist_y[3,1]: dist_y[3,1] + facility_position_y[1] = 0.75 +con_dist_y[3,2]: dist_y[3,2] + facility_position_y[2] = 0.75 +con_dist[0,0]: dist[0,0]*dist[0,0] - dist_x[0,0]*dist_x[0,0] - dist_y[0,0]*dist_y[0,0] >= 0 +con_dist[0,1]: dist[0,1]*dist[0,1] - dist_x[0,1]*dist_x[0,1] - dist_y[0,1]*dist_y[0,1] >= 0 +con_dist[0,2]: dist[0,2]*dist[0,2] - dist_x[0,2]*dist_x[0,2] - dist_y[0,2]*dist_y[0,2] >= 0 +con_dist[1,0]: dist[1,0]*dist[1,0] - dist_x[1,0]*dist_x[1,0] - dist_y[1,0]*dist_y[1,0] >= 0 +con_dist[1,1]: dist[1,1]*dist[1,1] - dist_x[1,1]*dist_x[1,1] - dist_y[1,1]*dist_y[1,1] >= 0 +con_dist[1,2]: dist[1,2]*dist[1,2] - dist_x[1,2]*dist_x[1,2] - dist_y[1,2]*dist_y[1,2] >= 0 +con_dist[2,0]: dist[2,0]*dist[2,0] - dist_x[2,0]*dist_x[2,0] - dist_y[2,0]*dist_y[2,0] >= 0 +con_dist[2,1]: dist[2,1]*dist[2,1] - dist_x[2,1]*dist_x[2,1] - dist_y[2,1]*dist_y[2,1] >= 0 +con_dist[2,2]: dist[2,2]*dist[2,2] - dist_x[2,2]*dist_x[2,2] - dist_y[2,2]*dist_y[2,2] >= 0 +con_dist[3,0]: dist[3,0]*dist[3,0] - dist_x[3,0]*dist_x[3,0] - dist_y[3,0]*dist_y[3,0] >= 0 +con_dist[3,1]: dist[3,1]*dist[3,1] - dist_x[3,1]*dist_x[3,1] - dist_y[3,1]*dist_y[3,1] >= 0 +con_dist[3,2]: dist[3,2]*dist[3,2] - dist_x[3,2]*dist_x[3,2] - dist_y[3,2]*dist_y[3,2] >= 0 +con_max_distance[0,0]: max_distance - dist[0,0] -2.828 is_closest[0,0] = -2.828 +con_max_distance[0,1]: max_distance - dist[0,1] -2.828 is_closest[0,1] = -2.828 +con_max_distance[0,2]: max_distance - dist[0,2] -2.828 is_closest[0,2] = -2.828 +con_max_distance[1,0]: max_distance - dist[1,0] -2.828 is_closest[1,0] = -2.828 +con_max_distance[1,1]: max_distance - dist[1,1] -2.828 is_closest[1,1] = -2.828 +con_max_distance[1,2]: max_distance - dist[1,2] -2.828 is_closest[1,2] = -2.828 +con_max_distance[2,0]: max_distance - dist[2,0] -2.828 is_closest[2,0] = -2.828 +con_max_distance[2,1]: max_distance - dist[2,1] -2.828 is_closest[2,1] = -2.828 +con_max_distance[2,2]: max_distance - dist[2,2] -2.828 is_closest[2,2] = -2.828 +con_max_distance[3,0]: max_distance - dist[3,0] -2.828 is_closest[3,0] = -2.828 +con_max_distance[3,1]: max_distance - dist[3,1] -2.828 is_closest[3,1] = -2.828 +con_max_distance[3,2]: max_distance - dist[3,2] -2.828 is_closest[3,2] = -2.828 + + +bounds + +facility_position_x[0] <= 1 +facility_position_x[1] <= 1 +facility_position_x[2] <= 1 +facility_position_y[0] <= 1 +facility_position_y[1] <= 1 +facility_position_y[2] <= 1 +is_closest[0,0] <= 1 +is_closest[0,1] <= 1 +is_closest[0,2] <= 1 +is_closest[1,0] <= 1 +is_closest[1,1] <= 1 +is_closest[1,2] <= 1 +is_closest[2,0] <= 1 +is_closest[2,1] <= 1 +is_closest[2,2] <= 1 +is_closest[3,0] <= 1 +is_closest[3,1] <= 1 +is_closest[3,2] <= 1 +-inf <= dist_x[0,0] +-inf <= dist_x[0,1] +-inf <= dist_x[0,2] +-inf <= dist_x[1,0] +-inf <= dist_x[1,1] +-inf <= dist_x[1,2] +-inf <= dist_x[2,0] +-inf <= dist_x[2,1] +-inf <= dist_x[2,2] +-inf <= dist_x[3,0] +-inf <= dist_x[3,1] +-inf <= dist_x[3,2] +-inf <= dist_y[0,0] +-inf <= dist_y[0,1] +-inf <= dist_y[0,2] +-inf <= dist_y[1,0] +-inf <= dist_y[1,1] +-inf <= dist_y[1,2] +-inf <= dist_y[2,0] +-inf <= dist_y[2,1] +-inf <= dist_y[2,2] +-inf <= dist_y[3,0] +-inf <= dist_y[3,1] +-inf <= dist_y[3,2] + + +binary + +is_closest[0,0] +is_closest[0,1] +is_closest[0,2] +is_closest[1,0] +is_closest[1,1] +is_closest[1,2] +is_closest[2,0] +is_closest[2,1] +is_closest[2,2] +is_closest[3,0] +is_closest[3,1] +is_closest[3,2] + +end diff --git a/tests/examples/facility_location/results/solution.sol b/tests/examples/facility_location/results/solution.sol new file mode 100644 index 0000000..0a372c6 --- /dev/null +++ b/tests/examples/facility_location/results/solution.sol @@ -0,0 +1,93 @@ +# Solution for model obj +# Objective value = 0 +max_distance 0 +is_closest[0,0] 0 +is_closest[0,1] 0 +is_closest[0,2] 1 +is_closest[1,0] 1 +is_closest[1,1] 0 +is_closest[1,2] 0 +is_closest[2,0] 0 +is_closest[2,1] 0 +is_closest[2,2] 1 +is_closest[3,0] 0 +is_closest[3,1] 0 +is_closest[3,2] 1 +dist_x[0,0] 0 +facility_position_x[0] 0 +dist_x[0,1] 0 +facility_position_x[1] 0 +dist_x[0,2] 0 +facility_position_x[2] 0 +dist_x[1,0] 0.25 +dist_x[1,1] 0.25 +dist_x[1,2] 0.25 +dist_x[2,0] 0.5 +dist_x[2,1] 0.5 +dist_x[2,2] 0.5 +dist_x[3,0] 0.75 +dist_x[3,1] 0.75 +dist_x[3,2] 0.75 +dist_y[0,0] 0 +facility_position_y[0] 0 +dist_y[0,1] 0 +facility_position_y[1] 0 +dist_y[0,2] 0 +facility_position_y[2] 0 +dist_y[1,0] 0.25 +dist_y[1,1] 0.25 +dist_y[1,2] 0.25 +dist_y[2,0] 0.5 +dist_y[2,1] 0.5 +dist_y[2,2] 0.5 +dist_y[3,0] 0.75 +dist_y[3,1] 0.75 +dist_y[3,2] 0.75 +dist[0,0]*dist[0,0] 2000000000 +dist_x[0,0]*dist_x[0,0] 0 +dist_y[0,0]*dist_y[0,0] 0 +dist[0,1]*dist[0,1] 2000000000 +dist_x[0,1]*dist_x[0,1] 0 +dist_y[0,1]*dist_y[0,1] 0 +dist[0,2]*dist[0,2] 2000000000 +dist_x[0,2]*dist_x[0,2] 0 +dist_y[0,2]*dist_y[0,2] 0 +dist[1,0]*dist[1,0] 2000000000 +dist_x[1,0]*dist_x[1,0] 0 +dist_y[1,0]*dist_y[1,0] 0 +dist[1,1]*dist[1,1] 2000000000 +dist_x[1,1]*dist_x[1,1] 0 +dist_y[1,1]*dist_y[1,1] 0 +dist[1,2]*dist[1,2] 2000000000 +dist_x[1,2]*dist_x[1,2] 0 +dist_y[1,2]*dist_y[1,2] 0 +dist[2,0]*dist[2,0] 2000000000 +dist_x[2,0]*dist_x[2,0] 0 +dist_y[2,0]*dist_y[2,0] 0 +dist[2,1]*dist[2,1] 2000000000 +dist_x[2,1]*dist_x[2,1] 0 +dist_y[2,1]*dist_y[2,1] 0 +dist[2,2]*dist[2,2] 2000000000 +dist_x[2,2]*dist_x[2,2] 0 +dist_y[2,2]*dist_y[2,2] 0 +dist[3,0]*dist[3,0] 2000000000 +dist_x[3,0]*dist_x[3,0] 0 +dist_y[3,0]*dist_y[3,0] 0 +dist[3,1]*dist[3,1] 2000000000 +dist_x[3,1]*dist_x[3,1] 0 +dist_y[3,1]*dist_y[3,1] 0 +dist[3,2]*dist[3,2] 2000000000 +dist_x[3,2]*dist_x[3,2] 0 +dist_y[3,2]*dist_y[3,2] 0 +dist[0,0] 2.828 +dist[0,1] 2.828 +dist[0,2] 0 +dist[1,0] 0 +dist[1,1] 2.828 +dist[1,2] 2.828 +dist[2,0] 2.828 +dist[2,1] 2.828 +dist[2,2] 0 +dist[3,0] 2.828 +dist[3,1] 2.828 +dist[3,2] 0 From 2518b3bfa612f188332aebcfb69724d3d38b7e35 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Mon, 11 Nov 2024 23:14:16 -0500 Subject: [PATCH 12/17] Make quadratic model work --- docs/quadratic_expressions.md | 4 +- src/pyoframe/core.py | 60 +++- src/pyoframe/io.py | 4 +- src/pyoframe/io_mappers.py | 3 +- src/pyoframe/model_element.py | 19 ++ tests/examples/facility_location/model.py | 90 +++-- .../results/pyoframe-problem.lp | 320 +++++++++++++----- .../facility_location/results/solution.sol | 216 +++++++----- tests/test_io.py | 33 ++ 9 files changed, 551 insertions(+), 198 deletions(-) diff --git a/docs/quadratic_expressions.md b/docs/quadratic_expressions.md index 6e866fb..76010f2 100644 --- a/docs/quadratic_expressions.md +++ b/docs/quadratic_expressions.md @@ -15,4 +15,6 @@ Internally, Pyoframe's `Expression` object is used for both linear and quadratic This extra column stores the ID of the second variable in quadratic terms. For terms with only one variable, this column contains ID `0` (a reserved variable ID which can thought of as meaning 'no variable'). The variables in a quadratic are rearranged such that the ID in the `__variable_id` column is always greater or equal than the variable ID in the `__quadratic_variable_id` (recall: a*b=b*a). This rearranging not only ensures that a*b+b*a=2a*b but also generates a useful property: If the variable ID in the first column (`__variable_id`) is `0` we know the variable ID in the second must also be `0` and therefore the term must be a constant. -The additional quadratic variable ID column is automatically dropped if through arithmetic the quadratic terms cancel out. \ No newline at end of file +The additional quadratic variable ID column is automatically dropped if through arithmetic the quadratic terms cancel out. + +One final detail: the `.to_str()` method on constraints and the objective has a special `quadratic_divider` parameter which is used when writing to a .lp file--a file format that requires quadratic terms to be enclosed by brackets and, in the case of the objective function, divided by 2. See [Gurobi Documentation](https://www.gurobi.com/documentation/current/refman/lp_format.html) \ No newline at end of file diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 9516624..1c57a47 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -870,8 +870,21 @@ def to_str_table( include_const_variable=False, var_map=None, float_precision=None, + quadratic_divider=None, ): + assert ( + quadratic_divider is None + or quadratic_divider == 1 + or quadratic_divider == 2 + ) data = self.data if include_const_term else self.variable_terms + if self.is_quadratic and quadratic_divider == 2: + data = data.with_columns( + pl.when(pl.col(QUAD_VAR_KEY) == CONST_TERM) + .then(pl.col(COEF_KEY)) + .otherwise(pl.col(COEF_KEY) * quadratic_divider) + ) + data = cast_coef_to_string(data, float_precision=float_precision) for var_column in self._variable_columns: @@ -896,10 +909,14 @@ def to_str_table( ).drop(temp_var_column) if self.is_quadratic: + if quadratic_divider is not None: + data = data.sort( + by=QUAD_VAR_KEY + ) # This ensures all elements in the brackets are grouped together data = data.with_columns( pl.when(pl.col(QUAD_VAR_KEY) == "") .then(pl.col(VAR_KEY)) - .otherwise(pl.concat_str(VAR_KEY, pl.lit("*"), pl.col(QUAD_VAR_KEY))) + .otherwise(pl.concat_str(VAR_KEY, pl.lit(" * "), pl.col(QUAD_VAR_KEY))) .alias(VAR_KEY) ).drop(QUAD_VAR_KEY) @@ -915,6 +932,43 @@ def to_str_table( ).drop(COEF_KEY, VAR_KEY) # Combine terms into one string + if self.is_quadratic and quadratic_divider is not None: + query_first = pl.col("expr").str.contains("*", literal=True).cum_sum() + query_last = pl.col("expr").str.contains("*", literal=True).cum_sum(reverse=True) + if dimensions is not None: + query_first = query_first.over(dimensions) + query_last = query_last.over(dimensions) + + data = data.with_columns( + pl.when( + pl.col("expr").str.contains("*", literal=True), + query_first == 1, + ) + .then( + pl.concat_str( + pl.lit("+ [ "), pl.col("expr").str.strip_chars(characters=" +") + ) + ) + .otherwise(pl.col("expr")) + .alias("expr") + ).with_columns( + pl.when( + pl.col("expr").str.contains("*", literal=True), query_last == 1 + ) + .then( + pl.concat_str( + pl.col("expr"), + pl.lit( + " ]" + if quadratic_divider == 1 + else f" ] / {quadratic_divider}" + ), + ) + ) + .otherwise(pl.col("expr")) + .alias("expr") + ) + if dimensions is not None: data = data.group_by(dimensions, maintain_order=True).agg( pl.col("expr").str.concat(delimiter=" ") @@ -961,6 +1015,7 @@ def to_str( include_header=False, include_data=True, float_precision=None, + quadratic_divider=None, ): result = "" if include_header: @@ -981,6 +1036,7 @@ def to_str( include_const_variable=include_const_variable, var_map=var_map, float_precision=float_precision, + quadratic_divider=quadratic_divider, ) if include_prefix: str_table = self.to_str_create_prefix(str_table) @@ -1224,6 +1280,7 @@ def to_str( var_map=None, float_precision=None, const_map=None, + quadratic_divider=None, ) -> str: dims = self.dimensions str_table = self.lhs.to_str_table( @@ -1231,6 +1288,7 @@ def to_str( max_rows=max_rows, include_const_term=False, var_map=var_map, + quadratic_divider=quadratic_divider, ) str_table = self.to_str_create_prefix(str_table, const_map=const_map) rhs = self.lhs.constant_terms.with_columns(pl.col(COEF_KEY) * -1) diff --git a/src/pyoframe/io.py b/src/pyoframe/io.py index c442de9..9ef9daf 100644 --- a/src/pyoframe/io.py +++ b/src/pyoframe/io.py @@ -86,7 +86,7 @@ def objective_to_file(m: "Model", f: TextIOWrapper, var_map): objective_sense = "minimize" if m.sense == ObjSense.MIN else "maximize" f.write(f"{objective_sense}\n\nobj:\n\n") result = m.objective.to_str( - var_map=var_map, include_prefix=False, include_const_variable=True + var_map=var_map, include_prefix=False, include_const_variable=True, quadratic_divider=2 ) f.write(result) @@ -99,7 +99,7 @@ def constraints_to_file(m: "Model", f: TextIOWrapper, var_map, const_map): f, "s.t.", ): - f.write(constraint.to_str(var_map=var_map, const_map=const_map) + "\n") + f.write(constraint.to_str(var_map=var_map, const_map=const_map, quadratic_divider=1) + "\n") def bounds_to_file(m: "Model", f, var_map): diff --git a/src/pyoframe/io_mappers.py b/src/pyoframe/io_mappers.py index 4263f14..82d3b3c 100644 --- a/src/pyoframe/io_mappers.py +++ b/src/pyoframe/io_mappers.py @@ -72,9 +72,8 @@ def undo(self, df: pl.DataFrame) -> pl.DataFrame: if df.height == 0: return df df = df.rename({self._ID_COL: Mapper.NAME_COL}) - # We use an inner join since sometimes Gurobi returns additional variables (i.e. quadratic variables) return df.join( - self.mapping_registry, on=Mapper.NAME_COL, how="inner", validate="m:1" + self.mapping_registry, on=Mapper.NAME_COL, how="left", validate="m:1" ).drop(Mapper.NAME_COL) diff --git a/src/pyoframe/model_element.py b/src/pyoframe/model_element.py index a5539fb..8183c91 100644 --- a/src/pyoframe/model_element.py +++ b/src/pyoframe/model_element.py @@ -140,6 +140,25 @@ def _new(self, data: pl.DataFrame): @abstractmethod def data(self): ... + + def select(self, **kwargs): + """ + Filter elements by the given criteria and then drop the filtered dimensions. + + Example: + >>> from pyoframe.core import Variable + >>> v = Variable([{"hour": ["00:00", "06:00", "12:00", "18:00"]}, {"city": ["Toronto", "Berlin", "Paris"]}]) + >>> v.select(hour="06:00") + + [Toronto]: x4 + [Berlin]: x5 + [Paris]: x6 + >>> v.select(hour="06:00", city="Toronto") + + x4 + """ + return self._new(self.data.filter(**kwargs).drop(kwargs.keys())) + class ModelElementWithId(ModelElement, AttrContainerMixin): """ diff --git a/tests/examples/facility_location/model.py b/tests/examples/facility_location/model.py index 9bb969c..0ec69a8 100644 --- a/tests/examples/facility_location/model.py +++ b/tests/examples/facility_location/model.py @@ -4,52 +4,102 @@ from pathlib import Path -def main(F, G, **kwargs): +def main(G, F, **kwargs): model = pf.Model("min") + g_range = range(G) model.facilities = pf.Set(f=range(F)) - model.customers = pf.Set(c=range(G)) + model.x_axis = pf.Set(x=g_range) + model.y_axis = pf.Set(y=g_range) + model.axis = pf.Set(d=[1, 2]) + model.customers = model.x_axis * model.y_axis - model.facility_position_x = pf.Variable(model.facilities, lb=0, ub=1) - model.facility_position_y = pf.Variable(model.facilities, lb=0, ub=1) + model.facility_position = pf.Variable(model.facilities, model.axis, lb=0, ub=1) model.customer_position_x = pl.DataFrame( - {"c": range(G), "x": [step / G for step in range(G)]} + {"x": g_range, "x_pos": [step / (G - 1) for step in g_range]} ).to_expr() model.customer_position_y = pl.DataFrame( - {"c": range(G), "y": [step / G for step in range(G)]} + {"y": g_range, "y_pos": [step / (G - 1) for step in g_range]} ).to_expr() model.max_distance = pf.Variable(lb=0) - model.objective = model.max_distance model.is_closest = pf.Variable(model.customers, model.facilities, vtype="binary") model.con_only_one_closest = pf.sum("f", model.is_closest) == 1 - model.dist_x = pf.Variable(model.customers, model.facilities) - model.dist_y = pf.Variable(model.customers, model.facilities) + model.dist_x = pf.Variable(model.x_axis, model.facilities) + model.dist_y = pf.Variable(model.y_axis, model.facilities) model.con_dist_x = model.dist_x == model.customer_position_x.add_dim( "f" - ) - model.facility_position_x.add_dim("c") + ) - model.facility_position.select(d=1).add_dim("x") model.con_dist_y = model.dist_y == model.customer_position_y.add_dim( "f" - ) - model.facility_position_y.add_dim("c") - model.dist = pf.Variable(model.customers, model.facilities, lb=0) - model.con_dist = model.dist**2 >= model.dist_x**2 + model.dist_y**2 + ) - model.facility_position.select(d=2).add_dim("y") + model.dist = pf.Variable(model.x_axis, model.y_axis, model.facilities, lb=0) + model.con_dist = model.dist**2 == (model.dist_x**2).add_dim("y") + ( + model.dist_y**2 + ).add_dim("x") - M = 2 * 1.414 # Twice the max distance which ensures that when is_closest is 0, the constraint is not binding. - model.con_max_distance = model.max_distance.add_dim("c", "f") == model.dist - M * ( - 1 - model.is_closest - ) + M = ( + 2 * 1.414 + ) # Twice the max distance which ensures that when is_closest is 0, the constraint is not binding. + model.con_max_distance = model.max_distance.add_dim( + "x", "y", "f" + ) >= model.dist - M * (1 - model.is_closest) + + model.objective = model.max_distance model.solve(**kwargs) + return model + + +def draw_results(model, G, T): + import tkinter + + root = tkinter.Tk() + scale = 500 + padding = 20 + canvas = tkinter.Canvas(root, width=scale + 2 * padding, height=scale + 2 * padding) + size = 5 + max_dist = model.max_distance.solution + for x in range(G): + for y in range(G): + canvas.create_rectangle( + x / (G - 1) * scale - size + padding, + y / (G - 1) * scale - size + padding, + x / (G - 1) * scale + size + padding, + y / (G - 1) * scale + size + padding, + fill="black", + ) + for f, x, y in model.facility_position.solution.pivot( + on="d", values="solution", index="f" + ).iter_rows(): + canvas.create_rectangle( + x * scale - size + padding, + y * scale - size + padding, + x * scale + size + padding, + y * scale + size + padding, + fill="red", + ) + canvas.create_oval( + (x - max_dist) * scale + padding, + (y - max_dist) * scale + padding, + (x + max_dist) * scale + padding, + (y + max_dist) * scale + padding, + outline="red", + ) + canvas.pack() + root.mainloop() if __name__ == "__main__": + G, F = 4, 3 working_dir = Path(os.path.dirname(os.path.realpath(__file__))) - main( - 3, - 4, + model = main( + G, + F, directory=working_dir / "results", use_var_names=True, solution_file=working_dir / "results" / "solution.sol", ) + draw_results(model, G, F) diff --git a/tests/examples/facility_location/results/pyoframe-problem.lp b/tests/examples/facility_location/results/pyoframe-problem.lp index eaa10b5..16ef277 100644 --- a/tests/examples/facility_location/results/pyoframe-problem.lp +++ b/tests/examples/facility_location/results/pyoframe-problem.lp @@ -6,80 +6,200 @@ max_distance s.t. -con_only_one_closest[0]: is_closest[0,0] + is_closest[0,1] + is_closest[0,2] = 1 -con_only_one_closest[1]: is_closest[1,0] + is_closest[1,1] + is_closest[1,2] = 1 -con_only_one_closest[2]: is_closest[2,0] + is_closest[2,1] + is_closest[2,2] = 1 -con_only_one_closest[3]: is_closest[3,0] + is_closest[3,1] + is_closest[3,2] = 1 -con_dist_x[0,0]: dist_x[0,0] + facility_position_x[0] = 0 -con_dist_x[0,1]: dist_x[0,1] + facility_position_x[1] = 0 -con_dist_x[0,2]: dist_x[0,2] + facility_position_x[2] = 0 -con_dist_x[1,0]: dist_x[1,0] + facility_position_x[0] = 0.25 -con_dist_x[1,1]: dist_x[1,1] + facility_position_x[1] = 0.25 -con_dist_x[1,2]: dist_x[1,2] + facility_position_x[2] = 0.25 -con_dist_x[2,0]: dist_x[2,0] + facility_position_x[0] = 0.5 -con_dist_x[2,1]: dist_x[2,1] + facility_position_x[1] = 0.5 -con_dist_x[2,2]: dist_x[2,2] + facility_position_x[2] = 0.5 -con_dist_x[3,0]: dist_x[3,0] + facility_position_x[0] = 0.75 -con_dist_x[3,1]: dist_x[3,1] + facility_position_x[1] = 0.75 -con_dist_x[3,2]: dist_x[3,2] + facility_position_x[2] = 0.75 -con_dist_y[0,0]: dist_y[0,0] + facility_position_y[0] = 0 -con_dist_y[0,1]: dist_y[0,1] + facility_position_y[1] = 0 -con_dist_y[0,2]: dist_y[0,2] + facility_position_y[2] = 0 -con_dist_y[1,0]: dist_y[1,0] + facility_position_y[0] = 0.25 -con_dist_y[1,1]: dist_y[1,1] + facility_position_y[1] = 0.25 -con_dist_y[1,2]: dist_y[1,2] + facility_position_y[2] = 0.25 -con_dist_y[2,0]: dist_y[2,0] + facility_position_y[0] = 0.5 -con_dist_y[2,1]: dist_y[2,1] + facility_position_y[1] = 0.5 -con_dist_y[2,2]: dist_y[2,2] + facility_position_y[2] = 0.5 -con_dist_y[3,0]: dist_y[3,0] + facility_position_y[0] = 0.75 -con_dist_y[3,1]: dist_y[3,1] + facility_position_y[1] = 0.75 -con_dist_y[3,2]: dist_y[3,2] + facility_position_y[2] = 0.75 -con_dist[0,0]: dist[0,0]*dist[0,0] - dist_x[0,0]*dist_x[0,0] - dist_y[0,0]*dist_y[0,0] >= 0 -con_dist[0,1]: dist[0,1]*dist[0,1] - dist_x[0,1]*dist_x[0,1] - dist_y[0,1]*dist_y[0,1] >= 0 -con_dist[0,2]: dist[0,2]*dist[0,2] - dist_x[0,2]*dist_x[0,2] - dist_y[0,2]*dist_y[0,2] >= 0 -con_dist[1,0]: dist[1,0]*dist[1,0] - dist_x[1,0]*dist_x[1,0] - dist_y[1,0]*dist_y[1,0] >= 0 -con_dist[1,1]: dist[1,1]*dist[1,1] - dist_x[1,1]*dist_x[1,1] - dist_y[1,1]*dist_y[1,1] >= 0 -con_dist[1,2]: dist[1,2]*dist[1,2] - dist_x[1,2]*dist_x[1,2] - dist_y[1,2]*dist_y[1,2] >= 0 -con_dist[2,0]: dist[2,0]*dist[2,0] - dist_x[2,0]*dist_x[2,0] - dist_y[2,0]*dist_y[2,0] >= 0 -con_dist[2,1]: dist[2,1]*dist[2,1] - dist_x[2,1]*dist_x[2,1] - dist_y[2,1]*dist_y[2,1] >= 0 -con_dist[2,2]: dist[2,2]*dist[2,2] - dist_x[2,2]*dist_x[2,2] - dist_y[2,2]*dist_y[2,2] >= 0 -con_dist[3,0]: dist[3,0]*dist[3,0] - dist_x[3,0]*dist_x[3,0] - dist_y[3,0]*dist_y[3,0] >= 0 -con_dist[3,1]: dist[3,1]*dist[3,1] - dist_x[3,1]*dist_x[3,1] - dist_y[3,1]*dist_y[3,1] >= 0 -con_dist[3,2]: dist[3,2]*dist[3,2] - dist_x[3,2]*dist_x[3,2] - dist_y[3,2]*dist_y[3,2] >= 0 -con_max_distance[0,0]: max_distance - dist[0,0] -2.828 is_closest[0,0] = -2.828 -con_max_distance[0,1]: max_distance - dist[0,1] -2.828 is_closest[0,1] = -2.828 -con_max_distance[0,2]: max_distance - dist[0,2] -2.828 is_closest[0,2] = -2.828 -con_max_distance[1,0]: max_distance - dist[1,0] -2.828 is_closest[1,0] = -2.828 -con_max_distance[1,1]: max_distance - dist[1,1] -2.828 is_closest[1,1] = -2.828 -con_max_distance[1,2]: max_distance - dist[1,2] -2.828 is_closest[1,2] = -2.828 -con_max_distance[2,0]: max_distance - dist[2,0] -2.828 is_closest[2,0] = -2.828 -con_max_distance[2,1]: max_distance - dist[2,1] -2.828 is_closest[2,1] = -2.828 -con_max_distance[2,2]: max_distance - dist[2,2] -2.828 is_closest[2,2] = -2.828 -con_max_distance[3,0]: max_distance - dist[3,0] -2.828 is_closest[3,0] = -2.828 -con_max_distance[3,1]: max_distance - dist[3,1] -2.828 is_closest[3,1] = -2.828 -con_max_distance[3,2]: max_distance - dist[3,2] -2.828 is_closest[3,2] = -2.828 +con_only_one_closest[0,0]: is_closest[0,0,0] + is_closest[0,0,1] + is_closest[0,0,2] = 1 +con_only_one_closest[0,1]: is_closest[0,1,0] + is_closest[0,1,1] + is_closest[0,1,2] = 1 +con_only_one_closest[0,2]: is_closest[0,2,0] + is_closest[0,2,1] + is_closest[0,2,2] = 1 +con_only_one_closest[0,3]: is_closest[0,3,0] + is_closest[0,3,1] + is_closest[0,3,2] = 1 +con_only_one_closest[1,0]: is_closest[1,0,0] + is_closest[1,0,1] + is_closest[1,0,2] = 1 +con_only_one_closest[1,1]: is_closest[1,1,0] + is_closest[1,1,1] + is_closest[1,1,2] = 1 +con_only_one_closest[1,2]: is_closest[1,2,0] + is_closest[1,2,1] + is_closest[1,2,2] = 1 +con_only_one_closest[1,3]: is_closest[1,3,0] + is_closest[1,3,1] + is_closest[1,3,2] = 1 +con_only_one_closest[2,0]: is_closest[2,0,0] + is_closest[2,0,1] + is_closest[2,0,2] = 1 +con_only_one_closest[2,1]: is_closest[2,1,0] + is_closest[2,1,1] + is_closest[2,1,2] = 1 +con_only_one_closest[2,2]: is_closest[2,2,0] + is_closest[2,2,1] + is_closest[2,2,2] = 1 +con_only_one_closest[2,3]: is_closest[2,3,0] + is_closest[2,3,1] + is_closest[2,3,2] = 1 +con_only_one_closest[3,0]: is_closest[3,0,0] + is_closest[3,0,1] + is_closest[3,0,2] = 1 +con_only_one_closest[3,1]: is_closest[3,1,0] + is_closest[3,1,1] + is_closest[3,1,2] = 1 +con_only_one_closest[3,2]: is_closest[3,2,0] + is_closest[3,2,1] + is_closest[3,2,2] = 1 +con_only_one_closest[3,3]: is_closest[3,3,0] + is_closest[3,3,1] + is_closest[3,3,2] = 1 +con_dist_x[0,0]: dist_x[0,0] + facility_position[0,1] = 0 +con_dist_x[0,1]: dist_x[0,1] + facility_position[1,1] = 0 +con_dist_x[0,2]: dist_x[0,2] + facility_position[2,1] = 0 +con_dist_x[1,0]: dist_x[1,0] + facility_position[0,1] = 0.3333333333333333 +con_dist_x[1,1]: dist_x[1,1] + facility_position[1,1] = 0.3333333333333333 +con_dist_x[1,2]: dist_x[1,2] + facility_position[2,1] = 0.3333333333333333 +con_dist_x[2,0]: dist_x[2,0] + facility_position[0,1] = 0.6666666666666666 +con_dist_x[2,1]: dist_x[2,1] + facility_position[1,1] = 0.6666666666666666 +con_dist_x[2,2]: dist_x[2,2] + facility_position[2,1] = 0.6666666666666666 +con_dist_x[3,0]: dist_x[3,0] + facility_position[0,1] = 1 +con_dist_x[3,1]: dist_x[3,1] + facility_position[1,1] = 1 +con_dist_x[3,2]: dist_x[3,2] + facility_position[2,1] = 1 +con_dist_y[0,0]: dist_y[0,0] + facility_position[0,2] = 0 +con_dist_y[0,1]: dist_y[0,1] + facility_position[1,2] = 0 +con_dist_y[0,2]: dist_y[0,2] + facility_position[2,2] = 0 +con_dist_y[1,0]: dist_y[1,0] + facility_position[0,2] = 0.3333333333333333 +con_dist_y[1,1]: dist_y[1,1] + facility_position[1,2] = 0.3333333333333333 +con_dist_y[1,2]: dist_y[1,2] + facility_position[2,2] = 0.3333333333333333 +con_dist_y[2,0]: dist_y[2,0] + facility_position[0,2] = 0.6666666666666666 +con_dist_y[2,1]: dist_y[2,1] + facility_position[1,2] = 0.6666666666666666 +con_dist_y[2,2]: dist_y[2,2] + facility_position[2,2] = 0.6666666666666666 +con_dist_y[3,0]: dist_y[3,0] + facility_position[0,2] = 1 +con_dist_y[3,1]: dist_y[3,1] + facility_position[1,2] = 1 +con_dist_y[3,2]: dist_y[3,2] + facility_position[2,2] = 1 +con_dist[0,0,0]: [ dist[0,0,0] * dist[0,0,0] - dist_x[0,0] * dist_x[0,0] - dist_y[0,0] * dist_y[0,0] ] = 0 +con_dist[0,0,1]: [ dist[0,0,1] * dist[0,0,1] - dist_x[0,1] * dist_x[0,1] - dist_y[0,1] * dist_y[0,1] ] = 0 +con_dist[0,0,2]: [ dist[0,0,2] * dist[0,0,2] - dist_x[0,2] * dist_x[0,2] - dist_y[0,2] * dist_y[0,2] ] = 0 +con_dist[0,1,0]: [ dist[0,1,0] * dist[0,1,0] - dist_x[0,0] * dist_x[0,0] - dist_y[1,0] * dist_y[1,0] ] = 0 +con_dist[0,1,1]: [ dist[0,1,1] * dist[0,1,1] - dist_x[0,1] * dist_x[0,1] - dist_y[1,1] * dist_y[1,1] ] = 0 +con_dist[0,1,2]: [ dist[0,1,2] * dist[0,1,2] - dist_x[0,2] * dist_x[0,2] - dist_y[1,2] * dist_y[1,2] ] = 0 +con_dist[0,2,0]: [ dist[0,2,0] * dist[0,2,0] - dist_x[0,0] * dist_x[0,0] - dist_y[2,0] * dist_y[2,0] ] = 0 +con_dist[0,2,1]: [ dist[0,2,1] * dist[0,2,1] - dist_x[0,1] * dist_x[0,1] - dist_y[2,1] * dist_y[2,1] ] = 0 +con_dist[0,2,2]: [ dist[0,2,2] * dist[0,2,2] - dist_x[0,2] * dist_x[0,2] - dist_y[2,2] * dist_y[2,2] ] = 0 +con_dist[0,3,0]: [ dist[0,3,0] * dist[0,3,0] - dist_x[0,0] * dist_x[0,0] - dist_y[3,0] * dist_y[3,0] ] = 0 +con_dist[0,3,1]: [ dist[0,3,1] * dist[0,3,1] - dist_x[0,1] * dist_x[0,1] - dist_y[3,1] * dist_y[3,1] ] = 0 +con_dist[0,3,2]: [ dist[0,3,2] * dist[0,3,2] - dist_x[0,2] * dist_x[0,2] - dist_y[3,2] * dist_y[3,2] ] = 0 +con_dist[1,0,0]: [ dist[1,0,0] * dist[1,0,0] - dist_x[1,0] * dist_x[1,0] - dist_y[0,0] * dist_y[0,0] ] = 0 +con_dist[1,0,1]: [ dist[1,0,1] * dist[1,0,1] - dist_x[1,1] * dist_x[1,1] - dist_y[0,1] * dist_y[0,1] ] = 0 +con_dist[1,0,2]: [ dist[1,0,2] * dist[1,0,2] - dist_x[1,2] * dist_x[1,2] - dist_y[0,2] * dist_y[0,2] ] = 0 +con_dist[1,1,0]: [ dist[1,1,0] * dist[1,1,0] - dist_x[1,0] * dist_x[1,0] - dist_y[1,0] * dist_y[1,0] ] = 0 +con_dist[1,1,1]: [ dist[1,1,1] * dist[1,1,1] - dist_x[1,1] * dist_x[1,1] - dist_y[1,1] * dist_y[1,1] ] = 0 +con_dist[1,1,2]: [ dist[1,1,2] * dist[1,1,2] - dist_x[1,2] * dist_x[1,2] - dist_y[1,2] * dist_y[1,2] ] = 0 +con_dist[1,2,0]: [ dist[1,2,0] * dist[1,2,0] - dist_x[1,0] * dist_x[1,0] - dist_y[2,0] * dist_y[2,0] ] = 0 +con_dist[1,2,1]: [ dist[1,2,1] * dist[1,2,1] - dist_x[1,1] * dist_x[1,1] - dist_y[2,1] * dist_y[2,1] ] = 0 +con_dist[1,2,2]: [ dist[1,2,2] * dist[1,2,2] - dist_x[1,2] * dist_x[1,2] - dist_y[2,2] * dist_y[2,2] ] = 0 +con_dist[1,3,0]: [ dist[1,3,0] * dist[1,3,0] - dist_x[1,0] * dist_x[1,0] - dist_y[3,0] * dist_y[3,0] ] = 0 +con_dist[1,3,1]: [ dist[1,3,1] * dist[1,3,1] - dist_x[1,1] * dist_x[1,1] - dist_y[3,1] * dist_y[3,1] ] = 0 +con_dist[1,3,2]: [ dist[1,3,2] * dist[1,3,2] - dist_x[1,2] * dist_x[1,2] - dist_y[3,2] * dist_y[3,2] ] = 0 +con_dist[2,0,0]: [ dist[2,0,0] * dist[2,0,0] - dist_x[2,0] * dist_x[2,0] - dist_y[0,0] * dist_y[0,0] ] = 0 +con_dist[2,0,1]: [ dist[2,0,1] * dist[2,0,1] - dist_x[2,1] * dist_x[2,1] - dist_y[0,1] * dist_y[0,1] ] = 0 +con_dist[2,0,2]: [ dist[2,0,2] * dist[2,0,2] - dist_x[2,2] * dist_x[2,2] - dist_y[0,2] * dist_y[0,2] ] = 0 +con_dist[2,1,0]: [ dist[2,1,0] * dist[2,1,0] - dist_x[2,0] * dist_x[2,0] - dist_y[1,0] * dist_y[1,0] ] = 0 +con_dist[2,1,1]: [ dist[2,1,1] * dist[2,1,1] - dist_x[2,1] * dist_x[2,1] - dist_y[1,1] * dist_y[1,1] ] = 0 +con_dist[2,1,2]: [ dist[2,1,2] * dist[2,1,2] - dist_x[2,2] * dist_x[2,2] - dist_y[1,2] * dist_y[1,2] ] = 0 +con_dist[2,2,0]: [ dist[2,2,0] * dist[2,2,0] - dist_x[2,0] * dist_x[2,0] - dist_y[2,0] * dist_y[2,0] ] = 0 +con_dist[2,2,1]: [ dist[2,2,1] * dist[2,2,1] - dist_x[2,1] * dist_x[2,1] - dist_y[2,1] * dist_y[2,1] ] = 0 +con_dist[2,2,2]: [ dist[2,2,2] * dist[2,2,2] - dist_x[2,2] * dist_x[2,2] - dist_y[2,2] * dist_y[2,2] ] = 0 +con_dist[2,3,0]: [ dist[2,3,0] * dist[2,3,0] - dist_x[2,0] * dist_x[2,0] - dist_y[3,0] * dist_y[3,0] ] = 0 +con_dist[2,3,1]: [ dist[2,3,1] * dist[2,3,1] - dist_x[2,1] * dist_x[2,1] - dist_y[3,1] * dist_y[3,1] ] = 0 +con_dist[2,3,2]: [ dist[2,3,2] * dist[2,3,2] - dist_x[2,2] * dist_x[2,2] - dist_y[3,2] * dist_y[3,2] ] = 0 +con_dist[3,0,0]: [ dist[3,0,0] * dist[3,0,0] - dist_x[3,0] * dist_x[3,0] - dist_y[0,0] * dist_y[0,0] ] = 0 +con_dist[3,0,1]: [ dist[3,0,1] * dist[3,0,1] - dist_x[3,1] * dist_x[3,1] - dist_y[0,1] * dist_y[0,1] ] = 0 +con_dist[3,0,2]: [ dist[3,0,2] * dist[3,0,2] - dist_x[3,2] * dist_x[3,2] - dist_y[0,2] * dist_y[0,2] ] = 0 +con_dist[3,1,0]: [ dist[3,1,0] * dist[3,1,0] - dist_x[3,0] * dist_x[3,0] - dist_y[1,0] * dist_y[1,0] ] = 0 +con_dist[3,1,1]: [ dist[3,1,1] * dist[3,1,1] - dist_x[3,1] * dist_x[3,1] - dist_y[1,1] * dist_y[1,1] ] = 0 +con_dist[3,1,2]: [ dist[3,1,2] * dist[3,1,2] - dist_x[3,2] * dist_x[3,2] - dist_y[1,2] * dist_y[1,2] ] = 0 +con_dist[3,2,0]: [ dist[3,2,0] * dist[3,2,0] - dist_x[3,0] * dist_x[3,0] - dist_y[2,0] * dist_y[2,0] ] = 0 +con_dist[3,2,1]: [ dist[3,2,1] * dist[3,2,1] - dist_x[3,1] * dist_x[3,1] - dist_y[2,1] * dist_y[2,1] ] = 0 +con_dist[3,2,2]: [ dist[3,2,2] * dist[3,2,2] - dist_x[3,2] * dist_x[3,2] - dist_y[2,2] * dist_y[2,2] ] = 0 +con_dist[3,3,0]: [ dist[3,3,0] * dist[3,3,0] - dist_x[3,0] * dist_x[3,0] - dist_y[3,0] * dist_y[3,0] ] = 0 +con_dist[3,3,1]: [ dist[3,3,1] * dist[3,3,1] - dist_x[3,1] * dist_x[3,1] - dist_y[3,1] * dist_y[3,1] ] = 0 +con_dist[3,3,2]: [ dist[3,3,2] * dist[3,3,2] - dist_x[3,2] * dist_x[3,2] - dist_y[3,2] * dist_y[3,2] ] = 0 +con_max_distance[0,0,0]: max_distance - dist[0,0,0] -2.828 is_closest[0,0,0] >= -2.828 +con_max_distance[0,0,1]: max_distance - dist[0,0,1] -2.828 is_closest[0,0,1] >= -2.828 +con_max_distance[0,0,2]: max_distance - dist[0,0,2] -2.828 is_closest[0,0,2] >= -2.828 +con_max_distance[0,1,0]: max_distance - dist[0,1,0] -2.828 is_closest[0,1,0] >= -2.828 +con_max_distance[0,1,1]: max_distance - dist[0,1,1] -2.828 is_closest[0,1,1] >= -2.828 +con_max_distance[0,1,2]: max_distance - dist[0,1,2] -2.828 is_closest[0,1,2] >= -2.828 +con_max_distance[0,2,0]: max_distance - dist[0,2,0] -2.828 is_closest[0,2,0] >= -2.828 +con_max_distance[0,2,1]: max_distance - dist[0,2,1] -2.828 is_closest[0,2,1] >= -2.828 +con_max_distance[0,2,2]: max_distance - dist[0,2,2] -2.828 is_closest[0,2,2] >= -2.828 +con_max_distance[0,3,0]: max_distance - dist[0,3,0] -2.828 is_closest[0,3,0] >= -2.828 +con_max_distance[0,3,1]: max_distance - dist[0,3,1] -2.828 is_closest[0,3,1] >= -2.828 +con_max_distance[0,3,2]: max_distance - dist[0,3,2] -2.828 is_closest[0,3,2] >= -2.828 +con_max_distance[1,0,0]: max_distance - dist[1,0,0] -2.828 is_closest[1,0,0] >= -2.828 +con_max_distance[1,0,1]: max_distance - dist[1,0,1] -2.828 is_closest[1,0,1] >= -2.828 +con_max_distance[1,0,2]: max_distance - dist[1,0,2] -2.828 is_closest[1,0,2] >= -2.828 +con_max_distance[1,1,0]: max_distance - dist[1,1,0] -2.828 is_closest[1,1,0] >= -2.828 +con_max_distance[1,1,1]: max_distance - dist[1,1,1] -2.828 is_closest[1,1,1] >= -2.828 +con_max_distance[1,1,2]: max_distance - dist[1,1,2] -2.828 is_closest[1,1,2] >= -2.828 +con_max_distance[1,2,0]: max_distance - dist[1,2,0] -2.828 is_closest[1,2,0] >= -2.828 +con_max_distance[1,2,1]: max_distance - dist[1,2,1] -2.828 is_closest[1,2,1] >= -2.828 +con_max_distance[1,2,2]: max_distance - dist[1,2,2] -2.828 is_closest[1,2,2] >= -2.828 +con_max_distance[1,3,0]: max_distance - dist[1,3,0] -2.828 is_closest[1,3,0] >= -2.828 +con_max_distance[1,3,1]: max_distance - dist[1,3,1] -2.828 is_closest[1,3,1] >= -2.828 +con_max_distance[1,3,2]: max_distance - dist[1,3,2] -2.828 is_closest[1,3,2] >= -2.828 +con_max_distance[2,0,0]: max_distance - dist[2,0,0] -2.828 is_closest[2,0,0] >= -2.828 +con_max_distance[2,0,1]: max_distance - dist[2,0,1] -2.828 is_closest[2,0,1] >= -2.828 +con_max_distance[2,0,2]: max_distance - dist[2,0,2] -2.828 is_closest[2,0,2] >= -2.828 +con_max_distance[2,1,0]: max_distance - dist[2,1,0] -2.828 is_closest[2,1,0] >= -2.828 +con_max_distance[2,1,1]: max_distance - dist[2,1,1] -2.828 is_closest[2,1,1] >= -2.828 +con_max_distance[2,1,2]: max_distance - dist[2,1,2] -2.828 is_closest[2,1,2] >= -2.828 +con_max_distance[2,2,0]: max_distance - dist[2,2,0] -2.828 is_closest[2,2,0] >= -2.828 +con_max_distance[2,2,1]: max_distance - dist[2,2,1] -2.828 is_closest[2,2,1] >= -2.828 +con_max_distance[2,2,2]: max_distance - dist[2,2,2] -2.828 is_closest[2,2,2] >= -2.828 +con_max_distance[2,3,0]: max_distance - dist[2,3,0] -2.828 is_closest[2,3,0] >= -2.828 +con_max_distance[2,3,1]: max_distance - dist[2,3,1] -2.828 is_closest[2,3,1] >= -2.828 +con_max_distance[2,3,2]: max_distance - dist[2,3,2] -2.828 is_closest[2,3,2] >= -2.828 +con_max_distance[3,0,0]: max_distance - dist[3,0,0] -2.828 is_closest[3,0,0] >= -2.828 +con_max_distance[3,0,1]: max_distance - dist[3,0,1] -2.828 is_closest[3,0,1] >= -2.828 +con_max_distance[3,0,2]: max_distance - dist[3,0,2] -2.828 is_closest[3,0,2] >= -2.828 +con_max_distance[3,1,0]: max_distance - dist[3,1,0] -2.828 is_closest[3,1,0] >= -2.828 +con_max_distance[3,1,1]: max_distance - dist[3,1,1] -2.828 is_closest[3,1,1] >= -2.828 +con_max_distance[3,1,2]: max_distance - dist[3,1,2] -2.828 is_closest[3,1,2] >= -2.828 +con_max_distance[3,2,0]: max_distance - dist[3,2,0] -2.828 is_closest[3,2,0] >= -2.828 +con_max_distance[3,2,1]: max_distance - dist[3,2,1] -2.828 is_closest[3,2,1] >= -2.828 +con_max_distance[3,2,2]: max_distance - dist[3,2,2] -2.828 is_closest[3,2,2] >= -2.828 +con_max_distance[3,3,0]: max_distance - dist[3,3,0] -2.828 is_closest[3,3,0] >= -2.828 +con_max_distance[3,3,1]: max_distance - dist[3,3,1] -2.828 is_closest[3,3,1] >= -2.828 +con_max_distance[3,3,2]: max_distance - dist[3,3,2] -2.828 is_closest[3,3,2] >= -2.828 bounds -facility_position_x[0] <= 1 -facility_position_x[1] <= 1 -facility_position_x[2] <= 1 -facility_position_y[0] <= 1 -facility_position_y[1] <= 1 -facility_position_y[2] <= 1 -is_closest[0,0] <= 1 -is_closest[0,1] <= 1 -is_closest[0,2] <= 1 -is_closest[1,0] <= 1 -is_closest[1,1] <= 1 -is_closest[1,2] <= 1 -is_closest[2,0] <= 1 -is_closest[2,1] <= 1 -is_closest[2,2] <= 1 -is_closest[3,0] <= 1 -is_closest[3,1] <= 1 -is_closest[3,2] <= 1 +facility_position[0,1] <= 1 +facility_position[0,2] <= 1 +facility_position[1,1] <= 1 +facility_position[1,2] <= 1 +facility_position[2,1] <= 1 +facility_position[2,2] <= 1 +is_closest[0,0,0] <= 1 +is_closest[0,0,1] <= 1 +is_closest[0,0,2] <= 1 +is_closest[0,1,0] <= 1 +is_closest[0,1,1] <= 1 +is_closest[0,1,2] <= 1 +is_closest[0,2,0] <= 1 +is_closest[0,2,1] <= 1 +is_closest[0,2,2] <= 1 +is_closest[0,3,0] <= 1 +is_closest[0,3,1] <= 1 +is_closest[0,3,2] <= 1 +is_closest[1,0,0] <= 1 +is_closest[1,0,1] <= 1 +is_closest[1,0,2] <= 1 +is_closest[1,1,0] <= 1 +is_closest[1,1,1] <= 1 +is_closest[1,1,2] <= 1 +is_closest[1,2,0] <= 1 +is_closest[1,2,1] <= 1 +is_closest[1,2,2] <= 1 +is_closest[1,3,0] <= 1 +is_closest[1,3,1] <= 1 +is_closest[1,3,2] <= 1 +is_closest[2,0,0] <= 1 +is_closest[2,0,1] <= 1 +is_closest[2,0,2] <= 1 +is_closest[2,1,0] <= 1 +is_closest[2,1,1] <= 1 +is_closest[2,1,2] <= 1 +is_closest[2,2,0] <= 1 +is_closest[2,2,1] <= 1 +is_closest[2,2,2] <= 1 +is_closest[2,3,0] <= 1 +is_closest[2,3,1] <= 1 +is_closest[2,3,2] <= 1 +is_closest[3,0,0] <= 1 +is_closest[3,0,1] <= 1 +is_closest[3,0,2] <= 1 +is_closest[3,1,0] <= 1 +is_closest[3,1,1] <= 1 +is_closest[3,1,2] <= 1 +is_closest[3,2,0] <= 1 +is_closest[3,2,1] <= 1 +is_closest[3,2,2] <= 1 +is_closest[3,3,0] <= 1 +is_closest[3,3,1] <= 1 +is_closest[3,3,2] <= 1 -inf <= dist_x[0,0] -inf <= dist_x[0,1] -inf <= dist_x[0,2] @@ -108,17 +228,53 @@ is_closest[3,2] <= 1 binary -is_closest[0,0] -is_closest[0,1] -is_closest[0,2] -is_closest[1,0] -is_closest[1,1] -is_closest[1,2] -is_closest[2,0] -is_closest[2,1] -is_closest[2,2] -is_closest[3,0] -is_closest[3,1] -is_closest[3,2] +is_closest[0,0,0] +is_closest[0,0,1] +is_closest[0,0,2] +is_closest[0,1,0] +is_closest[0,1,1] +is_closest[0,1,2] +is_closest[0,2,0] +is_closest[0,2,1] +is_closest[0,2,2] +is_closest[0,3,0] +is_closest[0,3,1] +is_closest[0,3,2] +is_closest[1,0,0] +is_closest[1,0,1] +is_closest[1,0,2] +is_closest[1,1,0] +is_closest[1,1,1] +is_closest[1,1,2] +is_closest[1,2,0] +is_closest[1,2,1] +is_closest[1,2,2] +is_closest[1,3,0] +is_closest[1,3,1] +is_closest[1,3,2] +is_closest[2,0,0] +is_closest[2,0,1] +is_closest[2,0,2] +is_closest[2,1,0] +is_closest[2,1,1] +is_closest[2,1,2] +is_closest[2,2,0] +is_closest[2,2,1] +is_closest[2,2,2] +is_closest[2,3,0] +is_closest[2,3,1] +is_closest[2,3,2] +is_closest[3,0,0] +is_closest[3,0,1] +is_closest[3,0,2] +is_closest[3,1,0] +is_closest[3,1,1] +is_closest[3,1,2] +is_closest[3,2,0] +is_closest[3,2,1] +is_closest[3,2,2] +is_closest[3,3,0] +is_closest[3,3,1] +is_closest[3,3,2] end diff --git a/tests/examples/facility_location/results/solution.sol b/tests/examples/facility_location/results/solution.sol index 0a372c6..def939c 100644 --- a/tests/examples/facility_location/results/solution.sol +++ b/tests/examples/facility_location/results/solution.sol @@ -1,93 +1,129 @@ # Solution for model obj -# Objective value = 0 -max_distance 0 -is_closest[0,0] 0 -is_closest[0,1] 0 -is_closest[0,2] 1 -is_closest[1,0] 1 -is_closest[1,1] 0 -is_closest[1,2] 0 -is_closest[2,0] 0 -is_closest[2,1] 0 -is_closest[2,2] 1 -is_closest[3,0] 0 -is_closest[3,1] 0 -is_closest[3,2] 1 +# Objective value = 4.9999957731881350e-01 +max_distance 4.9999957731881350e-01 +is_closest[0,0,0] 0 +is_closest[0,0,1] 1 +is_closest[0,0,2] 0 +is_closest[0,1,0] 1 +is_closest[0,1,1] 0 +is_closest[0,1,2] 0 +is_closest[0,2,0] 1 +is_closest[0,2,1] 0 +is_closest[0,2,2] 0 +is_closest[0,3,0] 1 +is_closest[0,3,1] 0 +is_closest[0,3,2] 0 +is_closest[1,0,0] 0 +is_closest[1,0,1] 1 +is_closest[1,0,2] 0 +is_closest[1,1,0] 1 +is_closest[1,1,1] 0 +is_closest[1,1,2] 0 +is_closest[1,2,0] 1 +is_closest[1,2,1] 0 +is_closest[1,2,2] 0 +is_closest[1,3,0] 1 +is_closest[1,3,1] 0 +is_closest[1,3,2] 0 +is_closest[2,0,0] 0 +is_closest[2,0,1] 1 +is_closest[2,0,2] 0 +is_closest[2,1,0] 0 +is_closest[2,1,1] 0 +is_closest[2,1,2] 1 +is_closest[2,2,0] 0 +is_closest[2,2,1] 0 +is_closest[2,2,2] 1 +is_closest[2,3,0] 0 +is_closest[2,3,1] 0 +is_closest[2,3,2] 1 +is_closest[3,0,0] 0 +is_closest[3,0,1] 1 +is_closest[3,0,2] 0 +is_closest[3,1,0] 0 +is_closest[3,1,1] 0 +is_closest[3,1,2] 1 +is_closest[3,2,0] 0 +is_closest[3,2,1] 0 +is_closest[3,2,2] 1 +is_closest[3,3,0] 0 +is_closest[3,3,1] 0 +is_closest[3,3,2] 1 dist_x[0,0] 0 -facility_position_x[0] 0 -dist_x[0,1] 0 -facility_position_x[1] 0 -dist_x[0,2] 0 -facility_position_x[2] 0 -dist_x[1,0] 0.25 -dist_x[1,1] 0.25 -dist_x[1,2] 0.25 -dist_x[2,0] 0.5 -dist_x[2,1] 0.5 -dist_x[2,2] 0.5 -dist_x[3,0] 0.75 -dist_x[3,1] 0.75 -dist_x[3,2] 0.75 -dist_y[0,0] 0 -facility_position_y[0] 0 +facility_position[0,1] 0 +dist_x[0,1] -4.9999957731881595e-01 +facility_position[1,1] 4.9999957731881595e-01 +dist_x[0,2] -9.4444444444444442e-01 +facility_position[2,1] 9.4444444444444442e-01 +dist_x[1,0] 3.3333333333333331e-01 +dist_x[1,1] -1.6666624398548263e-01 +dist_x[1,2] -6.1111111111111116e-01 +dist_x[2,0] 6.6666666666666663e-01 +dist_x[2,1] 1.6666708934785068e-01 +dist_x[2,2] -2.7777777777777779e-01 +dist_x[3,0] 1 +dist_x[3,1] 5.0000042268118461e-01 +dist_x[3,2] 5.5555555555555552e-02 +dist_y[0,0] -6.2732143666390394e-01 +facility_position[0,2] 6.2732143666390394e-01 dist_y[0,1] 0 -facility_position_y[1] 0 -dist_y[0,2] 0 -facility_position_y[2] 0 -dist_y[1,0] 0.25 -dist_y[1,1] 0.25 -dist_y[1,2] 0.25 -dist_y[2,0] 0.5 -dist_y[2,1] 0.5 -dist_y[2,2] 0.5 -dist_y[3,0] 0.75 -dist_y[3,1] 0.75 -dist_y[3,2] 0.75 -dist[0,0]*dist[0,0] 2000000000 -dist_x[0,0]*dist_x[0,0] 0 -dist_y[0,0]*dist_y[0,0] 0 -dist[0,1]*dist[0,1] 2000000000 -dist_x[0,1]*dist_x[0,1] 0 -dist_y[0,1]*dist_y[0,1] 0 -dist[0,2]*dist[0,2] 2000000000 -dist_x[0,2]*dist_x[0,2] 0 -dist_y[0,2]*dist_y[0,2] 0 -dist[1,0]*dist[1,0] 2000000000 -dist_x[1,0]*dist_x[1,0] 0 -dist_y[1,0]*dist_y[1,0] 0 -dist[1,1]*dist[1,1] 2000000000 -dist_x[1,1]*dist_x[1,1] 0 -dist_y[1,1]*dist_y[1,1] 0 -dist[1,2]*dist[1,2] 2000000000 -dist_x[1,2]*dist_x[1,2] 0 -dist_y[1,2]*dist_y[1,2] 0 -dist[2,0]*dist[2,0] 2000000000 -dist_x[2,0]*dist_x[2,0] 0 -dist_y[2,0]*dist_y[2,0] 0 -dist[2,1]*dist[2,1] 2000000000 -dist_x[2,1]*dist_x[2,1] 0 -dist_y[2,1]*dist_y[2,1] 0 -dist[2,2]*dist[2,2] 2000000000 -dist_x[2,2]*dist_x[2,2] 0 -dist_y[2,2]*dist_y[2,2] 0 -dist[3,0]*dist[3,0] 2000000000 -dist_x[3,0]*dist_x[3,0] 0 -dist_y[3,0]*dist_y[3,0] 0 -dist[3,1]*dist[3,1] 2000000000 -dist_x[3,1]*dist_x[3,1] 0 -dist_y[3,1]*dist_y[3,1] 0 -dist[3,2]*dist[3,2] 2000000000 -dist_x[3,2]*dist_x[3,2] 0 -dist_y[3,2]*dist_y[3,2] 0 -dist[0,0] 2.828 -dist[0,1] 2.828 -dist[0,2] 0 -dist[1,0] 0 -dist[1,1] 2.828 -dist[1,2] 2.828 -dist[2,0] 2.828 -dist[2,1] 2.828 -dist[2,2] 0 -dist[3,0] 2.828 -dist[3,1] 2.828 -dist[3,2] 0 +facility_position[1,2] 0 +dist_y[0,2] -6.6666666666666663e-01 +facility_position[2,2] 6.6666666666666663e-01 +dist_y[1,0] -2.9398810333057063e-01 +dist_y[1,1] 3.3333333333333331e-01 +dist_y[1,2] -3.3333333333333331e-01 +dist_y[2,0] 3.9345230002764076e-02 +dist_y[2,1] 6.6666666666666663e-01 +dist_y[2,2] 0 +dist_y[3,0] 3.7267856333609484e-01 +dist_y[3,1] 1 +dist_y[3,2] 3.3333333333333331e-01 +dist[0,0,0] 6.2732143666390272e-01 +dist[0,0,1] 4.9999957731881473e-01 +dist[0,0,2] 1.1560361710998912e+00 +dist[0,1,0] 2.9398810333056835e-01 +dist[0,1,1] 6.0092486088537456e-01 +dist[0,1,2] 1.0015419589602610e+00 +dist[0,2,0] 3.9345230002757123e-02 +dist[0,2,1] 8.3333307972469206e-01 +dist[0,2,2] 9.4444444444444442e-01 +dist[0,3,0] 3.7267856333609767e-01 +dist[0,3,1] 1.1180337997211873e+00 +dist[0,3,2] 1.0015419589602470e+00 +dist[1,0,0] 7.1038249979090418e-01 +dist[1,0,1] 1.6666624398547991e-01 +dist[1,0,2] 9.0437887757640079e-01 +dist[1,1,0] 4.4445485261274381e-01 +dist[1,1,1] 3.7267780722138583e-01 +dist[1,1,2] 6.9610905817480417e-01 +dist[1,2,0] 3.3564737185784826e-01 +dist[1,2,1] 6.8718416842115460e-01 +dist[1,2,2] 6.1111111111111116e-01 +dist[1,3,0] 5.0000042268118638e-01 +dist[1,3,1] 1.0137936855614309e+00 +dist[1,3,2] 6.9610905817478963e-01 +dist[2,0,0] 9.1541063427431679e-01 +dist[2,0,1] 1.6666708934784868e-01 +dist[2,0,2] 7.2222222222222221e-01 +dist[2,1,0] 7.2861062944782018e-01 +dist[2,1,1] 3.7267818527893048e-01 +dist[2,1,2] 4.3390271768024286e-01 +dist[2,2,0] 6.6782669276423312e-01 +dist[2,2,1] 6.8718437345164507e-01 +dist[2,2,2] 2.7777777777777779e-01 +dist[2,3,0] 7.6376289253583285e-01 +dist[2,3,1] 1.0137938245381486e+00 +dist[2,3,2] 4.3390271768022504e-01 +dist[3,0,0] 1.1804796418820878e+00 +dist[3,0,1] 5.0000042268118350e-01 +dist[3,0,2] 6.6897748205985585e-01 +dist[3,1,0] 1.0423190513944873e+00 +dist[3,1,1] 6.0092556426938215e-01 +dist[3,1,2] 3.3793126249254807e-01 +dist[3,2,0] 1.0007737242373866e+00 +dist[3,2,1] 8.3333358694211312e-01 +dist[3,2,2] 5.5555555555555552e-02 +dist[3,3,0] 1.0671875709406753e+00 +dist[3,3,1] 1.1180341777787315e+00 +dist[3,3,2] 3.3793126249252342e-01 diff --git a/tests/test_io.py b/tests/test_io.py index 337af48..2e5185b 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -49,5 +49,38 @@ def test_expression_with_const_to_str(): assert str(expr) == "2 x1 +5" +def test_quadratic_objective_to_str(): + kwargs = dict( + include_prefix=False, include_const_variable=True, quadratic_divider=2 + ) + expr = ( + 5 * Variable() - 2 * Variable() ** 2 + 3 + 2 * Variable() + 4 * Variable() ** 2 + ) + assert expr.to_str(**kwargs) == "5 x1 +3 x0 +2 x3 + [ -4 x2 * x2 +8 x4 * x4 ] / 2" + + +def test_quadratic_constraints_to_str(): + kwargs = dict(quadratic_divider=1) + constraint = ( + 5 * Variable() - 2 * Variable() ** 2 + 3 + 2 * Variable() + 4 * Variable() ** 2 + == 0 + ) + assert constraint.to_str(**kwargs) == "5 x1 +2 x3 + [ -2 x2 * x2 +4 x4 * x4 ] = -3" + df = pl.DataFrame({"x": [1, 2]}) + constraint = ( + 5 * Variable(df) + - 2 * Variable(df) ** 2 + + 3 + + 2 * Variable(df) + + 4 * Variable(df) ** 2 + == 0 + ) + assert ( + constraint.to_str(**kwargs) + == """[1]: 5 x5 +2 x9 + [ 4 x11 * x11 -2 x7 * x7 ] = -3 +[2]: 5 x6 +2 x10 + [ 4 x12 * x12 -2 x8 * x8 ] = -3""" + ) + + if __name__ == "__main__": pytest.main([__file__]) From 37b8342a82c912917695315bb54b74e3632dbff7 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Tue, 12 Nov 2024 09:43:12 -0500 Subject: [PATCH 13/17] lint --- src/pyoframe/core.py | 8 ++++---- src/pyoframe/io.py | 10 ++++++++-- src/pyoframe/model_element.py | 1 - 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 1c57a47..ddfce91 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -934,7 +934,9 @@ def to_str_table( # Combine terms into one string if self.is_quadratic and quadratic_divider is not None: query_first = pl.col("expr").str.contains("*", literal=True).cum_sum() - query_last = pl.col("expr").str.contains("*", literal=True).cum_sum(reverse=True) + query_last = ( + pl.col("expr").str.contains("*", literal=True).cum_sum(reverse=True) + ) if dimensions is not None: query_first = query_first.over(dimensions) query_last = query_last.over(dimensions) @@ -952,9 +954,7 @@ def to_str_table( .otherwise(pl.col("expr")) .alias("expr") ).with_columns( - pl.when( - pl.col("expr").str.contains("*", literal=True), query_last == 1 - ) + pl.when(pl.col("expr").str.contains("*", literal=True), query_last == 1) .then( pl.concat_str( pl.col("expr"), diff --git a/src/pyoframe/io.py b/src/pyoframe/io.py index 9ef9daf..99d4689 100644 --- a/src/pyoframe/io.py +++ b/src/pyoframe/io.py @@ -86,7 +86,10 @@ def objective_to_file(m: "Model", f: TextIOWrapper, var_map): objective_sense = "minimize" if m.sense == ObjSense.MIN else "maximize" f.write(f"{objective_sense}\n\nobj:\n\n") result = m.objective.to_str( - var_map=var_map, include_prefix=False, include_const_variable=True, quadratic_divider=2 + var_map=var_map, + include_prefix=False, + include_const_variable=True, + quadratic_divider=2, ) f.write(result) @@ -99,7 +102,10 @@ def constraints_to_file(m: "Model", f: TextIOWrapper, var_map, const_map): f, "s.t.", ): - f.write(constraint.to_str(var_map=var_map, const_map=const_map, quadratic_divider=1) + "\n") + f.write( + constraint.to_str(var_map=var_map, const_map=const_map, quadratic_divider=1) + + "\n" + ) def bounds_to_file(m: "Model", f, var_map): diff --git a/src/pyoframe/model_element.py b/src/pyoframe/model_element.py index 8183c91..4ccb90c 100644 --- a/src/pyoframe/model_element.py +++ b/src/pyoframe/model_element.py @@ -140,7 +140,6 @@ def _new(self, data: pl.DataFrame): @abstractmethod def data(self): ... - def select(self, **kwargs): """ Filter elements by the given criteria and then drop the filtered dimensions. From 49edbf0167274b2e410de1b3b518ba2028fdcb11 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Tue, 12 Nov 2024 10:05:02 -0500 Subject: [PATCH 14/17] Add facility_location as a valid test --- docs/quadratic_expressions.md | 16 +++++++++++++++- src/pyoframe/core.py | 6 +++--- tests/examples/facility_location/model.py | 6 +++--- .../{solution.sol => pyoframe-problem.sol} | 0 tests/test_examples.py | 19 ++++++++++++++----- 5 files changed, 35 insertions(+), 12 deletions(-) rename tests/examples/facility_location/results/{solution.sol => pyoframe-problem.sol} (100%) diff --git a/docs/quadratic_expressions.md b/docs/quadratic_expressions.md index 76010f2..66b04ae 100644 --- a/docs/quadratic_expressions.md +++ b/docs/quadratic_expressions.md @@ -4,10 +4,24 @@ Quadratic expressions work as you'd expect. Simply multiply two linear expressio ## Example -## TODO +### Maximize area of box +Here's a short example that shows that a square maximizes the area of any box with a fixed perimeter. ```python3 +import pyoframe as pf +model = pf.Model("max") +model.w = pf.Variable(lb=0) +model.h = pf.Variable(lb=0) +model.limit_perimter = 2 * (model.w + model.h) <= 20 +model.objective = model.w * model.h +model.solve() +print(f"It's a square: {model.w.solution==model.h.solution}") + +# Outputs: It's a square: True ``` +### Facility Location Problem + +See [examples/facility_location](../tests/examples/facility_location/). ## Note for Pyoframe developers: Internal Representation of Quadratics diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index ddfce91..f5cf637 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -107,7 +107,7 @@ def __pow__(self, power: int): >>> from pyoframe import Variable >>> Variable() ** 2 - x1*x1 + x1 * x1 """ if power == 2: return self * self @@ -726,10 +726,10 @@ def _add_const(self, const: int | float) -> Expression: x1 +5 >>> Variable() ** 2 + 5 - x2*x2 +5 + x2 * x2 +5 >>> Variable() ** 2 + Variable() + 5 - x3*x3 + x4 +5 + x3 * x3 + x4 +5 """ dim = self.dimensions data = self.data diff --git a/tests/examples/facility_location/model.py b/tests/examples/facility_location/model.py index 0ec69a8..373f247 100644 --- a/tests/examples/facility_location/model.py +++ b/tests/examples/facility_location/model.py @@ -4,7 +4,7 @@ from pathlib import Path -def main(G, F, **kwargs): +def main(G=4, F=3, **kwargs): model = pf.Model("min") g_range = range(G) @@ -53,7 +53,7 @@ def main(G, F, **kwargs): return model -def draw_results(model, G, T): +def draw_results(model, G, F): import tkinter root = tkinter.Tk() @@ -100,6 +100,6 @@ def draw_results(model, G, T): F, directory=working_dir / "results", use_var_names=True, - solution_file=working_dir / "results" / "solution.sol", + solution_file=working_dir / "results" / "pyoframe-problem.sol", ) draw_results(model, G, F) diff --git a/tests/examples/facility_location/results/solution.sol b/tests/examples/facility_location/results/pyoframe-problem.sol similarity index 100% rename from tests/examples/facility_location/results/solution.sol rename to tests/examples/facility_location/results/pyoframe-problem.sol diff --git a/tests/test_examples.py b/tests/test_examples.py index 51cb142..a15a515 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -27,6 +27,11 @@ class Example: many_valid_solutions=True, has_gurobi_version=False, ), + Example( + "facility_location", + has_gurobi_version=False, + many_valid_solutions=True, + ), ] @@ -49,15 +54,19 @@ def test_examples(example: Example): symbolic_solution_file = symbolic_output_dir / "pyoframe-problem.sol" dense_solution_file = dense_output_dir / "pyoframe-problem.sol" - symbolic_model = main_module.main( - input_dir, + symbolic_kwargs = dict( directory=symbolic_output_dir, solution_file=symbolic_solution_file, use_var_names=True, ) - dense_model = main_module.main( - input_dir, directory=dense_output_dir, solution_file=dense_solution_file - ) + dense_kwargs = dict(directory=dense_output_dir, solution_file=dense_solution_file) + + if input_dir.exists(): + symbolic_kwargs["input_dir"] = input_dir + dense_kwargs["input_dir"] = input_dir + + symbolic_model = main_module.main(**symbolic_kwargs) + dense_model = main_module.main(**dense_kwargs) if example.check_params is not None: for param, value in example.check_params.items(): From dd8ca2943e4dd2c23d9e43d212c37ab6193443dd Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Tue, 12 Nov 2024 10:10:15 -0500 Subject: [PATCH 15/17] Fix tests --- src/pyoframe/_arithmetic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index a6eefb8..e37ca96 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -106,9 +106,9 @@ def _quadratic_multiplication(self: "Expression", other: "Expression") -> "Expre >>> expr2 = df * x2 * 2 + 4 >>> expr1 * expr2 - [1]: 4 x1 +2 x2*x1 - [2]: 8 x1 +8 x2*x1 - [3]: 12 x1 +18 x2*x1 + [1]: 4 x1 +2 x2 * x1 + [2]: 8 x1 +8 x2 * x1 + [3]: 12 x1 +18 x2 * x1 >>> (expr1 * expr2) - df * x1 * df * x2 * 2 [1]: 4 x1 From 41613e1678b0c5ab6f22e2d0f747b85510ca8fbc Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Tue, 12 Nov 2024 10:28:20 -0500 Subject: [PATCH 16/17] Increase test coverage --- src/pyoframe/_arithmetic.py | 37 +++++++++++++++++++++++++------------ src/pyoframe/core.py | 6 +++++- 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/pyoframe/_arithmetic.py b/src/pyoframe/_arithmetic.py index e37ca96..a60be70 100644 --- a/src/pyoframe/_arithmetic.py +++ b/src/pyoframe/_arithmetic.py @@ -18,14 +18,35 @@ from pyoframe.core import Expression -def _multiply_expressions(*expressions: "Expression") -> "Expression": +def _multiply_expressions(self: "Expression", other: "Expression") -> "Expression": + """ + Multiply two or more expressions together. + + Examples: + >>> import pyoframe as pf + >>> m = pf.Model("min") + >>> m.x1 = pf.Variable() + >>> m.x2 = pf.Variable() + >>> m.x3 = pf.Variable() + >>> result = 5 * m.x1 * m.x2 + >>> result + + 5 x2 * x1 + >>> result * m.x3 + Traceback (most recent call last): + ... + pyoframe.constants.PyoframeError: Failed to multiply expressions: + * + Due to error: + Cannot multiply a quadratic expression by a non-constant. + """ try: - return _multiply_expressions_core(*expressions) + return _multiply_expressions_core(self, other) except PyoframeError as error: raise PyoframeError( "Failed to multiply expressions:\n" + " * ".join( - e.to_str(include_header=True, include_data=False) for e in expressions + e.to_str(include_header=True, include_data=False) for e in [self, other] ) + "\nDue to error:\n" + str(error) @@ -46,15 +67,7 @@ def _add_expressions(*expressions: "Expression") -> "Expression": ) from error -def _multiply_expressions_core(*expressions: "Expression") -> "Expression": - if len(expressions) > 2: - result = expressions[0] - for expr in expressions[1:]: - result = _multiply_expressions(result, expr) - return result - - self, other = expressions - +def _multiply_expressions_core(self: "Expression", other: "Expression") -> "Expression": self_degree, other_degree = self.degree(), other.degree() if self_degree + other_degree > 2: # We know one of the two must be a quadratic since 1 + 1 is not greater than 2. diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index f5cf637..0b93fe7 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -108,10 +108,14 @@ def __pow__(self, power: int): >>> Variable() ** 2 x1 * x1 + >>> Variable() ** 3 + Traceback (most recent call last): + ... + ValueError: Raising an expressions to **3 is not supported. Expressions can only be squared (**2). """ if power == 2: return self * self - return NotImplemented + raise ValueError(f"Raising an expressions to **{power} is not supported. Expressions can only be squared (**2).") def __neg__(self): res = self.to_expr() * -1 From e537d1ed0fac416b3e7973994ed90838f7d86805 Mon Sep 17 00:00:00 2001 From: Martin Staadecker Date: Tue, 12 Nov 2024 10:30:48 -0500 Subject: [PATCH 17/17] Lint --- src/pyoframe/core.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pyoframe/core.py b/src/pyoframe/core.py index 0b93fe7..a7e8617 100644 --- a/src/pyoframe/core.py +++ b/src/pyoframe/core.py @@ -115,7 +115,9 @@ def __pow__(self, power: int): """ if power == 2: return self * self - raise ValueError(f"Raising an expressions to **{power} is not supported. Expressions can only be squared (**2).") + raise ValueError( + f"Raising an expressions to **{power} is not supported. Expressions can only be squared (**2)." + ) def __neg__(self): res = self.to_expr() * -1