diff --git a/CHANGELOG.md b/CHANGELOG.md index ae37c3270..746e9c7bd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,8 @@ - Made `test_markDoNotAggrVar_and_getStatus` robust to SCIP presolve changes by discovering the aggregated/multi-aggregated variables instead of hardcoding them ### Changed - Move magic methods (`__radd__`, `__sub__`, `__rsub__`, `__rmul__`, `__richcmp__`, `__neg__`, and `__rtruediv__`) to `ExprLike` base class (#1204) -- Speed up `Expr.__add__` and `Expr.__iadd__` via the C-level API +- Speed up `Expr.__add__` and `Expr.__iadd__` via the C-level API (#1205) +- Replace Python math with C-level math functions and refactor unary expressions (#1224) - Extended `structured_optimization_trace` recipe to support context-managed JSONL tracing with final `run_end` records, alongside the existing attach-style in-memory tracing. ### Removed diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi index 3b232fea2..22f1eb5dd 100644 --- a/src/pyscipopt/expr.pxi +++ b/src/pyscipopt/expr.pxi @@ -42,7 +42,6 @@ # which should, in princple, modify the expr. However, since we do not implement __isub__, __sub__ # gets called (I guess) and so a copy is returned. # Modifying the expression directly would be a bug, given that the expression might be re-used by the user. -import math from typing import TYPE_CHECKING, Literal, Union import numpy as np @@ -54,6 +53,13 @@ from cpython.number cimport PyNumber_Check from cpython.object cimport Py_LE, Py_EQ, Py_GE, Py_TYPE from cpython.ref cimport PyObject from cpython.tuple cimport PyTuple_GET_ITEM +from libc.math cimport cos as c_cos +from libc.math cimport exp as c_exp +from libc.math cimport fabs as c_fabs +from libc.math cimport INFINITY +from libc.math cimport log as c_log +from libc.math cimport sqrt as c_sqrt +from libc.math cimport sin as c_sin cimport numpy as cnp from pyscipopt.scip cimport Variable, Solution @@ -281,23 +287,28 @@ cdef class ExprLike: def __pos__(self, /) -> Union[Expr, GenExpr]: return self.copy() - def __abs__(self) -> GenExpr: - return UnaryExpr(Operator.fabs, buildGenExprObj(self)) + def __abs__(self, /) -> AbsExpr: + return AbsExpr(Operator.fabs, buildGenExprObj(self)) - def exp(self) -> GenExpr: - return UnaryExpr(Operator.exp, buildGenExprObj(self)) + def exp(self, /) -> ExpExpr: + return ExpExpr(Operator.exp, buildGenExprObj(self)) - def log(self) -> GenExpr: - return UnaryExpr(Operator.log, buildGenExprObj(self)) + def log(self, /) -> LogExpr: + return LogExpr(Operator.log, buildGenExprObj(self)) - def sqrt(self) -> GenExpr: - return UnaryExpr(Operator.sqrt, buildGenExprObj(self)) + def sqrt(self, /) -> SqrtExpr: + return SqrtExpr(Operator.sqrt, buildGenExprObj(self)) - def sin(self) -> GenExpr: - return UnaryExpr(Operator.sin, buildGenExprObj(self)) + def sin(self, /) -> SinExpr: + return SinExpr(Operator.sin, buildGenExprObj(self)) - def cos(self) -> GenExpr: - return UnaryExpr(Operator.cos, buildGenExprObj(self)) + def cos(self, /) -> CosExpr: + return CosExpr(Operator.cos, buildGenExprObj(self)) + + cpdef double _evaluate(self, Solution sol) except *: + raise NotImplementedError( + f"{self.__class__.__name__!s} need to implement _evaluate() method" + ) cdef ExprLike copy(self, bint copy=True): raise NotImplementedError( @@ -711,7 +722,7 @@ cdef class GenExpr(ExprLike): def degree(self): '''Note: none of these expressions should be polynomial''' - return float('inf') + return INFINITY def getOp(self): '''returns operator of GenExpr''' @@ -828,24 +839,66 @@ cdef class PowExpr(GenExpr): return res -# Exp, Log, Sqrt, Sin, Cos Expressions cdef class UnaryExpr(GenExpr): + def __init__(self, op, expr): self.children = [] self.children.append(expr) self._op = op - def __abs__(self) -> UnaryExpr: - if self._op == "abs": - return self.copy() - return UnaryExpr(Operator.fabs, self) - - def __repr__(self): + def __repr__(self) -> str: return self._op + "(" + self.children[0].__repr__() + ")" + +cdef class AbsExpr(UnaryExpr): + + def __abs__(self) -> AbsExpr: + return self.copy() + + cpdef double _evaluate(self, Solution sol) except *: + return c_fabs((self.children[0])._evaluate(sol)) + + +cdef class ExpExpr(UnaryExpr): + + cpdef double _evaluate(self, Solution sol) except *: + return c_exp((self.children[0])._evaluate(sol)) + + +cdef class LogExpr(UnaryExpr): + + cpdef double _evaluate(self, Solution sol) except *: + cdef double val = (self.children[0])._evaluate(sol) + if val <= 0.0: + raise ValueError("math domain error") + return c_log(val) + + +cdef class SqrtExpr(UnaryExpr): + + cpdef double _evaluate(self, Solution sol) except *: + cdef double val = (self.children[0])._evaluate(sol) + if val < 0.0: + raise ValueError("math domain error") + return c_sqrt(val) + + +cdef class SinExpr(UnaryExpr): + + cpdef double _evaluate(self, Solution sol) except *: + cdef double val = (self.children[0])._evaluate(sol) + if c_fabs(val) == INFINITY: + raise ValueError("math domain error") + return c_sin(val) + + +cdef class CosExpr(UnaryExpr): + cpdef double _evaluate(self, Solution sol) except *: - cdef double res = (self.children[0])._evaluate(sol) - return math.fabs(res) if self._op == "abs" else getattr(math, self._op)(res) + cdef double val = (self.children[0])._evaluate(sol) + if c_fabs(val) == INFINITY: + raise ValueError("math domain error") + return c_cos(val) # class for constant expressions diff --git a/src/pyscipopt/scip.pxd b/src/pyscipopt/scip.pxd index 3a126f4dc..4144f91ca 100644 --- a/src/pyscipopt/scip.pxd +++ b/src/pyscipopt/scip.pxd @@ -2156,13 +2156,12 @@ cdef extern from "tpi/tpi.h": cdef class ExprLike: + cpdef double _evaluate(self, Solution sol) cdef ExprLike copy(self, bint copy=*) cdef class Expr(ExprLike): cdef public terms - cpdef double _evaluate(self, Solution sol) - cdef class Event: cdef SCIP_EVENT* event # can be used to store problem data diff --git a/src/pyscipopt/scip.pyi b/src/pyscipopt/scip.pyi index 1e3164109..ca9da7d77 100644 --- a/src/pyscipopt/scip.pyi +++ b/src/pyscipopt/scip.pyi @@ -355,12 +355,12 @@ class ExprLike: **kwargs: Incomplete, ) -> Incomplete: ... def __pos__(self, /) -> Self: ... - def __abs__(self, /) -> UnaryExpr: ... - def exp(self) -> UnaryExpr: ... - def log(self) -> UnaryExpr: ... - def sqrt(self) -> UnaryExpr: ... - def sin(self) -> UnaryExpr: ... - def cos(self) -> UnaryExpr: ... + def __abs__(self, /) -> AbsExpr: ... + def exp(self, /) -> ExpExpr: ... + def log(self, /) -> LogExpr: ... + def sqrt(self, /) -> SqrtExpr: ... + def sin(self, /) -> SinExpr: ... + def cos(self, /) -> CosExpr: ... @overload def __le__(self, other: float | ExprLike, /) -> ExprCons: ... @overload @@ -588,7 +588,7 @@ class MatrixConstraint(np.ndarray): def isStickingAtNode(self) -> Incomplete: ... class MatrixExpr(np.ndarray): - def _evaluate(self, sol: Incomplete) -> Incomplete: ... + def _evaluate(self, sol: Solution) -> Incomplete: ... def __array_ufunc__( self, ufunc: np.ufunc, @@ -2342,6 +2342,13 @@ class Term: class UnaryExpr(GenExpr): def __init__(self, *args: Incomplete, **kwargs: Incomplete) -> None: ... +class AbsExpr(UnaryExpr): ... +class ExpExpr(UnaryExpr): ... +class LogExpr(UnaryExpr): ... +class SqrtExpr(UnaryExpr): ... +class SinExpr(UnaryExpr): ... +class CosExpr(UnaryExpr): ... + @disjoint_base class VarExpr(GenExpr): var: Incomplete diff --git a/tests/test_expr.py b/tests/test_expr.py index 1b51e4f2a..36341e518 100644 --- a/tests/test_expr.py +++ b/tests/test_expr.py @@ -212,8 +212,16 @@ def test_getVal_with_GenExpr(): assert m.getVal(y / x) == 2 # test "**(prod(1.0,**(sum(0.0,prod(1.0,x)),-1)),2)" assert m.getVal((1 / x) ** 2) == 1 - # test "sin(sum(0.0,prod(1.0,x)))" + + # test C-level math functions + assert m.getVal(abs(x)) == 1 + assert m.getVal(abs(-x)) == 1 + assert m.getVal(abs(abs(-x))) == 1 + assert round(m.getVal(exp(x)), 6) == round(math.exp(1), 6) + assert round(m.getVal(log(x)), 6) == round(math.log(1), 6) + assert round(m.getVal(sqrt(x)), 6) == round(math.sqrt(1), 6) assert round(m.getVal(sin(x)), 6) == round(math.sin(1), 6) + assert round(m.getVal(cos(x)), 6) == round(math.cos(1), 6) with pytest.raises(TypeError): m.getVal(1) @@ -221,6 +229,33 @@ def test_getVal_with_GenExpr(): with pytest.raises(ZeroDivisionError): m.getVal(1 / z) + # math domain errors match the math module + with pytest.raises(ValueError, match="math domain error"): + m.getVal(log(z)) # log(0) + + with pytest.raises(ValueError, match="math domain error"): + m.getVal(log(-y)) # log(-2) + + with pytest.raises(ValueError, match="math domain error"): + m.getVal(sqrt(-y)) # sqrt(-2) + + # sqrt(0) is inside the domain, like math.sqrt(0) + assert m.getVal(sqrt(z)) == 0 + + # +inf is inside log's domain, like math.log(inf) -> inf + assert m.getVal(log(math.inf)) == math.inf + + # sin and cos reject infinite arguments, like math.sin(inf) + with pytest.raises(ValueError, match="math domain error"): + m.getVal(sin(math.inf)) + + with pytest.raises(ValueError, match="math domain error"): + m.getVal(cos(-math.inf)) + + # nested unary expressions propagate the inner domain error + with pytest.raises(ValueError, match="math domain error"): + m.getVal(exp(log(-x))) + def test_unary_ufunc(model): m, x, y, z = model