From c999721253f6f59c93212ccaab403001645bf5ff Mon Sep 17 00:00:00 2001 From: Felipe Serrano Date: Wed, 28 Feb 2018 17:25:22 +0100 Subject: [PATCH 1/3] add support general nonlinear cons inspired in SCIP.jl, pyomo and biogeme's expressions --- src/pyscipopt/__init__.py | 3 + src/pyscipopt/expr.pxi | 435 +++++++++++++++++++++++++++++++++++++- src/pyscipopt/scip.pyx | 98 ++++++++- tests/test_expr.py | 189 +++++++++++++++++ tests/test_linexpr.py | 10 - 5 files changed, 715 insertions(+), 20 deletions(-) create mode 100644 tests/test_expr.py diff --git a/src/pyscipopt/__init__.py b/src/pyscipopt/__init__.py index c5d7f4c26..a612d09be 100644 --- a/src/pyscipopt/__init__.py +++ b/src/pyscipopt/__init__.py @@ -13,6 +13,9 @@ from pyscipopt.scip import Sepa from pyscipopt.scip import LP from pyscipopt.scip import quicksum +from pyscipopt.scip import exp +from pyscipopt.scip import log +from pyscipopt.scip import sqrt from pyscipopt.scip import PY_SCIP_RESULT as SCIP_RESULT from pyscipopt.scip import PY_SCIP_PARAMSETTING as SCIP_PARAMSETTING from pyscipopt.scip import PY_SCIP_PARAMEMPHASIS as SCIP_PARAMEMPHASIS diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi index 2b0dd4fab..66e7a5b50 100644 --- a/src/pyscipopt/expr.pxi +++ b/src/pyscipopt/expr.pxi @@ -1,3 +1,47 @@ +# In this file we implemenet the handling of expressions +# We have two types of expressions: Expr and GenExpr. +# The Expr can only handle polynomial expressions. +# In addition, one can recover easily information from them. +# A polynomial is a dictionary between `terms` and coefficients. +# A `term` is a tuple of variables +# For examples, 2*x*x*y*z - 1.3 x*y*y + 1 is stored as a +# {Term(x,x,y,z) : 2, Term(x,y,y) : -1.3, Term() : 1} +# Addition of common terms and expansion of exponents occur automatically. +# Given the way `Expr`s are stored, it is easy to access the terms: e.g. +# expr = 2*x*x*y*z - 1.3 x*y*y + 1 +# expr[Term(x,x,y,z)] returns 1.3 +# expr[Term(x)] returns 0.0 +# +# On the other hand, when dealing with expressions more general than polynomials, +# that is, absolute values, exp, log, sqrt or any general exponent, we use GenExpr. +# GenExpr stores expression trees in a rudimentary way. +# Basically, it stores the operator and the list of children. +# We have different types of general expressions that in addition +# to the operation and list of children stores +# SumExpr: coefficients and constant +# ProdExpr: constant +# Constant: constant +# VarExpr: variable +# PowExpr: exponent +# UnaryExpr: nothing +# We do not provide any way of accessing the internal information of the expression tree, +# nor we simplify common terms or do any other type of simplification. +# The `GenExpr` is pass as is to SCIP and SCIP will do what it see fits during presolving. +# +# TODO: All this is very complicated, so we might wanna unify Expr and GenExpr. +# Maybe when consexpr is released it makes sense to revisit this. +# TODO: We have to think about the operations that we define: __isub__, __add__, etc +# and when to copy expressions and when to not copy them. +# For example: when creating a ExprCons from an Expr expr, we store the expression expr +# and then we normalize. When doing the normalization, we do +# ``` +# c = self.expr[CONST] +# self.expr -= c +# ``` +# 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. + def _is_number(e): try: @@ -11,21 +55,21 @@ def _is_number(e): def _expr_richcmp(self, other, op): if op == 1: # <= - if isinstance(other, Expr): + if isinstance(other, Expr) or isinstance(other, GenExpr): return (self - other) <= 0.0 elif _is_number(other): return ExprCons(self, rhs=float(other)) else: raise NotImplementedError elif op == 5: # >= - if isinstance(other, Expr): + if isinstance(other, Expr) or isinstance(other, GenExpr): return (self - other) >= 0.0 elif _is_number(other): return ExprCons(self, lhs=float(other)) else: raise NotImplementedError elif op == 2: # == - if isinstance(other, Expr): + if isinstance(other, Expr) or isinstance(other, GenExpr): return (self - other) == 0.0 elif _is_number(other): return ExprCons(self, lhs=float(other), rhs=float(other)) @@ -66,6 +110,32 @@ class Term: CONST = Term() +# helper function +def buildGenExprObj(expr): + if isinstance(expr, int) or isinstance(expr, float): + return Constant(expr) + elif isinstance(expr, Expr): + # loop over terms and create a sumexpr with the sum of each term + # each term is either a variable (which gets transformed into varexpr) + # or a product of variables (which gets tranformed into a prod) + sumexpr = SumExpr() + for vars, coef in expr.terms.items(): + if len(vars) == 0: + sumexpr += coef + elif len(vars) == 1: + varexpr = VarExpr(vars[0]) + sumexpr += coef * varexpr + else: + prodexpr = ProdExpr() + for v in vars: + varexpr = VarExpr(v) + prodexpr *= varexpr + sumexpr += coef * prodexpr + return sumexpr + else: + assert isinstance(expr, GenExpr) + return expr + cdef class Expr: '''Polynomial expressions of variables with operator overloading.''' cdef public terms @@ -84,6 +154,9 @@ cdef class Expr: key = Term(key) return self.terms.get(key, 0.0) + def __abs__(self): + return abs(buildGenExprObj(self)) + def __add__(self, other): left = self right = other @@ -100,6 +173,8 @@ cdef class Expr: elif _is_number(right): c = float(right) terms[CONST] = terms.get(CONST, 0.0) + c + elif isinstance(right, GenExpr): + return buildGenExprObj(left) + right else: raise NotImplementedError return Expr(terms) @@ -111,6 +186,11 @@ cdef class Expr: elif _is_number(other): c = float(other) self.terms[CONST] = self.terms.get(CONST, 0.0) + c + elif isinstance(other, GenExpr): + # is no longer in place, might affect performance? + # can't do `self = buildGenExprObj(self) + other` since I get + # TypeError: Cannot convert pyscipopt.scip.SumExpr to pyscipopt.scip.Expr + return buildGenExprObj(self) + other else: raise NotImplementedError return self @@ -129,14 +209,25 @@ cdef class Expr: v = v1 + v2 terms[v] = terms.get(v, 0.0) + c1 * c2 return Expr(terms) + elif isinstance(other, GenExpr): + return buildGenExprObj(self) * other else: raise NotImplementedError + def __div__(self, other): + ''' transforms Expr into GenExpr''' + divisor = buildGenExprObj(other) + # we can't divide by 0 + if divisor.getOp() == Operator.const: + assert divisor.number != 0.0 + + return buildGenExprObj(self) * divisor**(-1) + def __pow__(self, other, modulo): if float(other).is_integer() and other >= 0: exp = int(other) - else: - raise NotImplementedError + else: # need to transform to GenExpr + return buildGenExprObj(self)**other res = 1 for _ in range(exp): @@ -189,15 +280,20 @@ cdef class ExprCons: def normalize(self): '''move constant terms in expression to bounds''' - c = self.expr[CONST] - self.expr -= c + if isinstance(self.expr, Expr): + c = self.expr[CONST] + self.expr -= c + assert self.expr[CONST] == 0.0 + self.expr.normalize() + else: + assert isinstance(self.expr, GenExpr) + return + if not self.lhs is None: self.lhs -= c if not self.rhs is None: self.rhs -= c - assert self.expr[CONST] == 0.0 - self.expr.normalize() def __richcmp__(self, other, op): '''turn it into a constraint''' @@ -247,3 +343,324 @@ def quicksum(termlist): for term in termlist: result += term return result + + +class Op: + const = 'const' + varidx = 'var' + exp, log, sqrt = 'exp','log', 'sqrt' + plus, minus, mul, div, power = '+', '-', '*', '/', '**' + add = 'sum' + prod = 'prod' + fabs = 'abs' + operatorIndexDic={ + varidx:SCIP_EXPR_VARIDX, + const:SCIP_EXPR_CONST, + plus:SCIP_EXPR_PLUS, + minus:SCIP_EXPR_MINUS, + mul:SCIP_EXPR_MUL, + div:SCIP_EXPR_DIV, + sqrt:SCIP_EXPR_SQRT, + power:SCIP_EXPR_REALPOWER, + exp:SCIP_EXPR_EXP, + log:SCIP_EXPR_LOG, + fabs:SCIP_EXPR_ABS, + add:SCIP_EXPR_SUM, + prod:SCIP_EXPR_PRODUCT + } + def getOpIndex(self, op): + return Op.operatorIndexDic[op]; + +Operator = Op() + +cdef class GenExpr: + '''General expressions of variables with operator overloading. + + Notes: + - this expressions are not smart enough to identify equal terms + - in constrast to polynomial expressions, __getitem__ is not implemented + so expr[x] will generate an error instead of returning the coefficient of x + ''' + cdef public operatorIndex + cdef public op + cdef public children + + + def __init__(self): # do we need it + ''' ''' + + def __abs__(self): + ans = UnaryExpr() + ans.op = Operator.fabs + ans.operatorIndex = Operator.operatorIndexDic[self.op] + ans.children.append(self) + return ans + + def __add__(self, other): + left = buildGenExprObj(self) + right = buildGenExprObj(other) + ans = SumExpr() + + # add left term + if left.getOp() == Operator.add: + ans.coefs.extend(left.coefs) + ans.children.extend(left.children) + ans.constant += left.constant + elif left.getOp() == Operator.const: + ans.constant += left.number + else: + ans.coefs.append(1.0) + ans.children.append(left) + + # add right term + if right.getOp() == Operator.add: + ans.coefs.extend(right.coefs) + ans.children.extend(right.children) + ans.constant += right.constant + elif right.getOp() == Operator.const: + ans.constant += right.number + else: + ans.coefs.append(1.0) + ans.children.append(right) + + return ans + + #def __iadd__(self, other): + #''' in-place addition, i.e., expr += other ''' + # assert isinstance(self, Expr) + # right = buildGenExprObj(other) + # + # # transform self into sum + # if self.getOp() != Operator.add: + # newsum = SumExpr() + # if self.getOp() == Operator.const: + # newsum.constant += self.number + # else: + # newsum.coefs.append(1.0) + # newsum.children.append(self.copy()) # TODO: what is copy? + # self = newsum + # # add right term + # if right.getOp() == Operator.add: + # self.coefs.extend(right.coefs) + # self.children.extend(right.children) + # self.constant += right.constant + # elif right.getOp() == Operator.const: + # self.constant += right.number + # else: + # self.coefs.append(1.0) + # self.children.append(right) + # return self + + def __mul__(self, other): + left = buildGenExprObj(self) + right = buildGenExprObj(other) + ans = ProdExpr() + + # multiply left factor + if left.getOp() == Operator.prod: + ans.children.extend(left.children) + ans.constant *= left.constant + elif left.getOp() == Operator.const: + ans.constant *= left.number + else: + ans.children.append(left) + + # multiply right factor + if right.getOp() == Operator.prod: + ans.children.extend(right.children) + ans.constant *= right.constant + elif right.getOp() == Operator.const: + ans.constant *= right.number + else: + ans.children.append(right) + + return ans + + #def __imul__(self, other): + #''' in-place multiplication, i.e., expr *= other ''' + # assert isinstance(self, Expr) + # right = buildGenExprObj(other) + # # transform self into prod + # if self.getOp() != Operator.prod: + # newprod = ProdExpr() + # if self.getOp() == Operator.const: + # newprod.constant *= self.number + # else: + # newprod.children.append(self.copy()) # TODO: what is copy? + # self = newprod + # # multiply right factor + # if right.getOp() == Operator.prod: + # self.children.extend(right.children) + # self.constant *= right.constant + # elif right.getOp() == Operator.const: + # self.constant *= right.number + # else: + # self.children.append(right) + # return self + + def __pow__(self, other, modulo): + expo = buildGenExprObj(other) + if expo.getOp() != Operator.const: + raise NotImplementedError("exponents must be numbers") + if self.getOp() == Operator.const: + return Constant(self.number**expo.number) + ans = PowExpr() + ans.children.append(self) + ans.expo = expo.number + + return ans + + #TODO: ipow, idiv, true_div, etc + def __div__(self, other): + divisor = buildGenExprObj(other) + # we can't divide by 0 + if divisor.getOp() == Operator.const: + assert divisor.number != 0.0 + + return self * divisor**(-1) + + def __neg__(self): + return -1.0 * self + + def __sub__(self, other): + return self + (-other) + + def __radd__(self, other): + return self.__add__(other) + + def __rmul__(self, other): + return self.__mul__(other) + + def __rsub__(self, other): + return -1.0 * self + other + + def __richcmp__(self, other, op): + '''turn it into a constraint''' + return _expr_richcmp(self, other, op) + + def degree(self): + return float('inf') # none of these expressions should be polynomial + + def getOp(self): + return self.op + + +# Sum Expressions +cdef class SumExpr(GenExpr): + + cdef public constant + cdef public coefs + + def __init__(self): + self.constant = 0.0 + self.coefs = [] + self.children = [] + self.op = Operator.add + self.operatorIndex = Operator.operatorIndexDic[self.op] + def __repr__(self): + return self.op + "(" + str(self.constant) + "," + ",".join(map(lambda child : child.__repr__(), self.children)) + ")" + +# Prod Expressions +cdef class ProdExpr(GenExpr): + cdef public constant + def __init__(self): + self.constant = 1.0 + self.children = [] + self.op = Operator.prod + self.operatorIndex = Operator.operatorIndexDic[self.op] + def __repr__(self): + return self.op + "(" + str(self.constant) + "," + ",".join(map(lambda child : child.__repr__(), self.children)) + ")" + +# Var Expressions +cdef class VarExpr(GenExpr): + cdef public var + def __init__(self, var): + self.children = [var] + self.op = Operator.varidx + self.operatorIndex = Operator.operatorIndexDic[self.op] + def __repr__(self): + return self.children[0].__repr__() + +# Pow Expressions +cdef class PowExpr(GenExpr): + cdef public expo + def __init__(self): + self.expo = 1.0 + self.children = [] + self.op = Operator.power + self.operatorIndex = Operator.operatorIndexDic[self.op] + def __repr__(self): + return self.op + "(" + self.children[0].__repr__() + "," + str(self.expo) + ")" + +# Exp, Log, Sqrt Expressions +cdef class UnaryExpr(GenExpr): + def __init__(self): + self.children = [] + def __repr__(self): + return self.op + "(" + self.children[0].__repr__() + ")" + +# class for constant expressions +cdef class Constant(GenExpr): + cdef public number + def __init__(self,number): + self.number = number + self.op = Operator.const + self.operatorIndex = Operator.operatorIndexDic[self.op] + + def __repr__(self): + return str(self.number) + +def exp(expr): + ans = UnaryExpr() + ans.op = Operator.exp + ans.operatorIndex = Operator.operatorIndexDic[ans.op] + ans.children.append(buildGenExprObj(expr)) + return ans +def log(expr): + ans = UnaryExpr() + ans.op = Operator.log + ans.operatorIndex = Operator.operatorIndexDic[ans.op] + ans.children.append(buildGenExprObj(expr)) + return ans +def sqrt(expr): + ans = UnaryExpr() + ans.op = Operator.sqrt + ans.operatorIndex = Operator.operatorIndexDic[ans.op] + ans.children.append(buildGenExprObj(expr)) + return ans + +# transforms tree to an array of nodes. each node is an operator and the position of the +# children of that operator (i.e. the other nodes) in the array +def expr_to_nodes(expr): + assert isinstance(expr, GenExpr) + nodes = [] + expr_to_array(expr, nodes) + return nodes + +def value_to_array(val, nodes): + nodes.append(tuple(['const', [val]])) + return len(nodes) - 1 + +# there many hacky things here: value_to_array is trying to mimick +# the multiple dispatch of julia. Also that we have to ask which expression is which +# in order to get the constants correctly +# also, for sums, we are not considering coefficients, because basically all coefficients are 1 +# haven't even consider substractions, but I guess we would interpret them as a - b = a + (-1) * b +def expr_to_array(expr, nodes): + op = expr.op + if op != Operator.varidx: + indices = [] + nchildren = len(expr.children) + for child in expr.children: + pos = expr_to_array(child, nodes) # position of child in the final array of nodes, 'nodes' + indices.append(pos) + if op == Operator.power: + pos = value_to_array(expr.expo, nodes) + indices.append(pos) + elif (op == Operator.add and expr.constant != 0.0) or (op == Operator.prod and expr.constant != 1.0): + pos = value_to_array(expr.constant, nodes) + indices.append(pos) + nodes.append( tuple( [op, indices] ) ) + else: # var + nodes.append( tuple( [op, expr.children] ) ) + return len(nodes) - 1 diff --git a/src/pyscipopt/scip.pyx b/src/pyscipopt/scip.pyx index a65fada2e..7de6d8d63 100644 --- a/src/pyscipopt/scip.pyx +++ b/src/pyscipopt/scip.pyx @@ -891,6 +891,8 @@ cdef class Model: return self._addLinCons(cons, **kwargs) elif deg <= 2: return self._addQuadCons(cons, **kwargs) + elif deg == float('inf'): # general nonlinear + return self._addGenNonlinearCons(cons, **kwargs) else: return self._addNonlinearCons(cons, **kwargs) @@ -1008,6 +1010,100 @@ cdef class Model: free(varexprs) return PyCons + def _addGenNonlinearCons(self, ExprCons cons, **kwargs): + cdef SCIP_EXPR** childrenexpr + cdef SCIP_EXPR** scipexprs + cdef SCIP_EXPRTREE* exprtree + cdef SCIP_CONS* scip_cons + cdef int nchildren + + # get arrays from python's expression tree + expr = cons.expr + nodes = expr_to_nodes(expr) + op2idx = Operator.operatorIndexDic + + # in nodes we have a list of tuples: each tuple is of the form + # (operator, [indices]) where indices are the indices of the tuples + # that are the children of this operator. This is sorted, + # so we are going to do is: + # loop over the nodes and create the expression of each + # Note1: when the operator is SCIP_EXPR_CONST, [indices] stores the value + # Note2: we need to compute the number of variable operators to find out + # how many variables are there. + nvars = 0 + for node in nodes: + if op2idx[node[0]] == SCIP_EXPR_VARIDX: + nvars += 1 + vars = malloc(nvars * sizeof(SCIP_VAR*)) + + varpos = 0 + scipexprs = malloc(len(nodes) * sizeof(SCIP_EXPR*)) + for i,node in enumerate(nodes): + op = node[0] + opidx = op2idx[op] + if opidx == SCIP_EXPR_VARIDX: + assert len(node[1]) == 1 + pyvar = node[1][0] # for vars we store the actual var! + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], opidx, varpos) ) + vars[varpos] = (pyvar).var + varpos += 1 + continue + if opidx == SCIP_EXPR_CONST: + assert len(node[1]) == 1 + value = node[1][0] + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], opidx, value) ) + continue + if opidx == SCIP_EXPR_SUM or opidx == SCIP_EXPR_PRODUCT: + nchildren = len(node[1]) + childrenexpr = malloc(nchildren * sizeof(SCIP_EXPR*)) + for c, pos in enumerate(node[1]): + childrenexpr[c] = scipexprs[pos] + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], opidx, nchildren, childrenexpr) ) + + free(childrenexpr); + continue + if opidx == SCIP_EXPR_REALPOWER: + # the second child is the exponent which is a const + valuenode = nodes[node[1][1]] + assert op2idx[valuenode[0]] == SCIP_EXPR_CONST + exponent = valuenode[1][0] + if float(exponent).is_integer(): + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], SCIP_EXPR_INTPOWER, scipexprs[node[1][0]], exponent) ) + else: + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], opidx, scipexprs[node[1][0]], exponent) ) + continue + if opidx == SCIP_EXPR_EXP or opidx == SCIP_EXPR_LOG or opidx == SCIP_EXPR_SQRT or opidx == SCIP_EXPR_ABS: + assert len(node[1]) == 1 + PY_SCIP_CALL( SCIPexprCreate(SCIPblkmem(self._scip), &scipexprs[i], opidx, scipexprs[node[1][0]]) ) + continue + # default: + raise NotImplementedError + assert varpos == nvars + + # create expression tree + PY_SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(self._scip), &exprtree, scipexprs[len(nodes) - 1], nvars, 0, NULL) ); + PY_SCIP_CALL( SCIPexprtreeSetVars(exprtree, nvars, vars) ); + + # create nonlinear constraint for exprtree + PY_SCIP_CALL( SCIPcreateConsNonlinear( + self._scip, &scip_cons, str_conversion(kwargs['name']), + 0, NULL, NULL, # linear + 1, &exprtree, NULL, # nonlinear + kwargs['lhs'], kwargs['rhs'], + kwargs['initial'], kwargs['separate'], kwargs['enforce'], + kwargs['check'], kwargs['propagate'], kwargs['local'], + kwargs['modifiable'], kwargs['dynamic'], kwargs['removable'], + kwargs['stickingatnode']) ) + PY_SCIP_CALL(SCIPaddCons(self._scip, scip_cons)) + PyCons = Constraint.create(scip_cons) + PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons)) + PY_SCIP_CALL( SCIPexprtreeFree(&exprtree) ) + + # free more memory + free(vars) + + return PyCons + def addConsCoeff(self, Constraint cons, Variable var, coeff): """Add coefficient to the linear constraint (if non-zero). @@ -1177,7 +1273,7 @@ cdef class Model: raise ValueError("expected inequality that has either only a left or right hand side") if cons.expr.degree() > 1: - raise ValueError("expected linear inequality, expression has degree %d" % cons.expr.degree) + raise ValueError("expected linear inequality, expression has degree %d" % cons.expr.degree()) assert cons.expr.degree() <= 1 diff --git a/tests/test_expr.py b/tests/test_expr.py new file mode 100644 index 000000000..20d5e9912 --- /dev/null +++ b/tests/test_expr.py @@ -0,0 +1,189 @@ +import pytest + +from pyscipopt import Model, sqrt, log, exp +from pyscipopt.scip import Expr, GenExpr, ExprCons, Term, quicksum + +@pytest.fixture(scope="module") +def model(): + m = Model() + x = m.addVar("x") + y = m.addVar("y") + z = m.addVar("z") + return m, x, y, z + +CONST = Term() + +def test_upgrade(model): + m, x, y, z = model + expr = x + y + assert isinstance(expr, Expr) + expr += exp(z) + assert isinstance(expr, GenExpr) + + expr = x + y + assert isinstance(expr, Expr) + expr -= exp(z) + assert isinstance(expr, GenExpr) + + expr = x + y + assert isinstance(expr, Expr) + expr /= x + assert isinstance(expr, GenExpr) + + expr = x + y + assert isinstance(expr, Expr) + expr *= sqrt(x) + assert isinstance(expr, GenExpr) + + expr = x + y + assert isinstance(expr, Expr) + expr **= 1.5 + assert isinstance(expr, GenExpr) + + expr = x + y + assert isinstance(expr, Expr) + assert isinstance(expr + exp(x), GenExpr) + assert isinstance(expr - exp(x), GenExpr) + assert isinstance(expr/x, GenExpr) + assert isinstance(expr * x**1.2, GenExpr) + assert isinstance(sqrt(expr), GenExpr) + assert isinstance(abs(expr), GenExpr) + assert isinstance(log(expr), GenExpr) + assert isinstance(exp(expr), GenExpr) + +def test_genexpr_op_expr(model): + m, x, y, z = model + genexpr = x**1.5 + y + assert isinstance(genexpr, GenExpr) + genexpr += x**2 + assert isinstance(genexpr, GenExpr) + genexpr += 1 + assert isinstance(genexpr, GenExpr) + genexpr += x + assert isinstance(genexpr, GenExpr) + genexpr += 2 * y + assert isinstance(genexpr, GenExpr) + genexpr -= x**2 + assert isinstance(genexpr, GenExpr) + genexpr -= 1 + assert isinstance(genexpr, GenExpr) + genexpr -= x + assert isinstance(genexpr, GenExpr) + genexpr -= 2 * y + assert isinstance(genexpr, GenExpr) + genexpr *= x + y + assert isinstance(genexpr, GenExpr) + genexpr *= 2 + assert isinstance(genexpr, GenExpr) + genexpr /= 2 + assert isinstance(genexpr, GenExpr) + genexpr /= x + y + assert isinstance(genexpr, GenExpr) + assert isinstance(x**1.2 + x + y, GenExpr) + assert isinstance(x**1.2 - x, GenExpr) + assert isinstance(x**1.2 *(x+y), GenExpr) + +def test_genexpr_op_genexpr(model): + m, x, y, z = model + genexpr = x**1.5 + y + assert isinstance(genexpr, GenExpr) + genexpr **= 2.2 + assert isinstance(genexpr, GenExpr) + genexpr += exp(x) + assert isinstance(genexpr, GenExpr) + genexpr -= exp(x) + assert isinstance(genexpr, GenExpr) + genexpr /= log(x + 1) + assert isinstance(genexpr, GenExpr) + genexpr *= (x + y)**1.2 + assert isinstance(genexpr, GenExpr) + genexpr /= exp(2) + assert isinstance(genexpr, GenExpr) + genexpr /= x + y + assert isinstance(genexpr, GenExpr) + genexpr = x**1.5 + y + assert isinstance(genexpr, GenExpr) + assert isinstance(sqrt(x) + genexpr, GenExpr) + assert isinstance(exp(x) + genexpr, GenExpr) + assert isinstance(1/x + genexpr, GenExpr) + assert isinstance(1/x**1.5 - genexpr, GenExpr) + assert isinstance(y/x - exp(genexpr), GenExpr) + # sqrt(2) is not a constant expression and + # we can only power to constant expressions! + with pytest.raises(NotImplementedError): + genexpr **= sqrt(2) + +def test_degree(model): + m, x, y, z = model + expr = GenExpr() + assert expr.degree() == float('inf') + +# In contrast to Expr inequalities, we can't expect much of the sides +def test_inequality(model): + m, x, y, z = model + + expr = x + 2*y + assert isinstance(expr, Expr) + cons = expr <= x**1.2 + assert isinstance(cons, ExprCons) + assert isinstance(cons.expr, GenExpr) + assert cons.lhs is None + assert cons.rhs == 0.0 + + assert isinstance(expr, Expr) + cons = expr >= x**1.2 + assert isinstance(cons, ExprCons) + assert isinstance(cons.expr, GenExpr) + assert cons.lhs == 0.0 + assert cons.rhs is None + + assert isinstance(expr, Expr) + cons = expr >= 1 + x**1.2 + assert isinstance(cons, ExprCons) + assert isinstance(cons.expr, GenExpr) + assert cons.lhs == 0.0 # NOTE: the 1 is pass the the other side because of the way GenExprs work + assert cons.rhs is None + + assert isinstance(expr, Expr) + cons = exp(expr) <= 1 + x**1.2 + assert isinstance(cons, ExprCons) + assert isinstance(cons.expr, GenExpr) + assert cons.rhs == 0.0 + assert cons.lhs is None + + +def test_equation(model): + m, x, y, z = model + equat = 2*x**1.2 - 3*sqrt(y) == 1 + assert isinstance(equat, ExprCons) + assert equat.lhs == equat.rhs + assert equat.lhs == 1.0 + + equat = exp(x+2*y) == 1 + x**1.2 + assert isinstance(equat, ExprCons) + assert isinstance(equat.expr, GenExpr) + assert equat.lhs == equat.rhs + assert equat.lhs == 0.0 + + equat = x == 1 + x**1.2 + assert isinstance(equat, ExprCons) + assert isinstance(equat.expr, GenExpr) + assert equat.lhs == equat.rhs + assert equat.lhs == 0.0 + +def test_objective(model): + m, x, y, z = model + + # setting linear objective + m.setObjective(x + y) + + # using quicksum + m.setObjective(quicksum(2*v for v in [x, y, z])) + + # setting nonlinear objective + with pytest.raises(ValueError): + m.setObjective(x**2 - y*z) + + # setting affine objective + with pytest.raises(ValueError): + m.setObjective(x + y + 1) diff --git a/tests/test_linexpr.py b/tests/test_linexpr.py index 18416d64d..ab54dcc09 100644 --- a/tests/test_linexpr.py +++ b/tests/test_linexpr.py @@ -109,16 +109,6 @@ def test_operations_poly(model): assert expr[Term(y,y)] == 2.0 assert expr.terms == (x**3 + 2*y**2).terms -def test_invalid_power(model): - m, x, y, z = model - assert (x + (y + 1)**0).terms == (x + 1).terms - - with pytest.raises(NotImplementedError): - expr = (x + 1)**0.5 - - with pytest.raises(NotImplementedError): - expr = (x + 1)**(-1) - def test_degree(model): m, x, y, z = model expr = Expr() From 4a8c024af5ed5342f5ec79cff02bc128e06ee3ab Mon Sep 17 00:00:00 2001 From: Felipe Serrano Date: Fri, 2 Mar 2018 12:51:40 +0100 Subject: [PATCH 2/3] make it python3 compatible; remove old test --- src/pyscipopt/expr.pxi | 30 +++++++++++++++++++++++++++++- tests/test_expr.py | 17 ----------------- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi index 66e7a5b50..4f8e2913f 100644 --- a/src/pyscipopt/expr.pxi +++ b/src/pyscipopt/expr.pxi @@ -223,6 +223,20 @@ cdef class Expr: return buildGenExprObj(self) * divisor**(-1) + def __rdiv__(self, other): + ''' other / self ''' + otherexpr = buildGenExprObj(other) + return otherexpr.__div__(self) + + def __truediv__(self,other): + selfexpr = buildGenExprObj(self) + return selfexpr.__div__(other) + + def __rtruediv__(self, other): + ''' other / self ''' + otherexpr = buildGenExprObj(other) + return otherexpr.__div__(self) + def __pow__(self, other, modulo): if float(other).is_integer() and other >= 0: exp = int(other) @@ -510,7 +524,7 @@ cdef class GenExpr: return ans - #TODO: ipow, idiv, true_div, etc + #TODO: ipow, idiv, etc def __div__(self, other): divisor = buildGenExprObj(other) # we can't divide by 0 @@ -519,6 +533,20 @@ cdef class GenExpr: return self * divisor**(-1) + def __rdiv__(self, other): + ''' other / self ''' + otherexpr = buildGenExprObj(other) + return otherexpr.__div__(self) + + def __truediv__(self,other): + selfexpr = buildGenExprObj(self) + return selfexpr.__div__(other) + + def __rtruediv__(self, other): + ''' other / self ''' + otherexpr = buildGenExprObj(other) + return otherexpr.__div__(self) + def __neg__(self): return -1.0 * self diff --git a/tests/test_expr.py b/tests/test_expr.py index 20d5e9912..5539b24f0 100644 --- a/tests/test_expr.py +++ b/tests/test_expr.py @@ -170,20 +170,3 @@ def test_equation(model): assert isinstance(equat.expr, GenExpr) assert equat.lhs == equat.rhs assert equat.lhs == 0.0 - -def test_objective(model): - m, x, y, z = model - - # setting linear objective - m.setObjective(x + y) - - # using quicksum - m.setObjective(quicksum(2*v for v in [x, y, z])) - - # setting nonlinear objective - with pytest.raises(ValueError): - m.setObjective(x**2 - y*z) - - # setting affine objective - with pytest.raises(ValueError): - m.setObjective(x + y + 1) From e6227a4f80c759246a70d85ddb7d082b81fad470 Mon Sep 17 00:00:00 2001 From: Felipe Serrano Date: Fri, 2 Mar 2018 14:48:42 +0100 Subject: [PATCH 3/3] make python3 work second try --- src/pyscipopt/expr.pxi | 56 ++++++++++++++++-------------------------- tests/test_expr.py | 3 +++ 2 files changed, 24 insertions(+), 35 deletions(-) diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi index 4f8e2913f..726cb620d 100644 --- a/src/pyscipopt/expr.pxi +++ b/src/pyscipopt/expr.pxi @@ -216,12 +216,8 @@ cdef class Expr: def __div__(self, other): ''' transforms Expr into GenExpr''' - divisor = buildGenExprObj(other) - # we can't divide by 0 - if divisor.getOp() == Operator.const: - assert divisor.number != 0.0 - - return buildGenExprObj(self) * divisor**(-1) + selfexpr = buildGenExprObj(self) + return selfexpr.__div__(other) def __rdiv__(self, other): ''' other / self ''' @@ -230,12 +226,12 @@ cdef class Expr: def __truediv__(self,other): selfexpr = buildGenExprObj(self) - return selfexpr.__div__(other) + return selfexpr.__truediv__(other) def __rtruediv__(self, other): ''' other / self ''' otherexpr = buildGenExprObj(other) - return otherexpr.__div__(self) + return otherexpr.__truediv__(self) def __pow__(self, other, modulo): if float(other).is_integer() and other >= 0: @@ -404,11 +400,7 @@ cdef class GenExpr: ''' ''' def __abs__(self): - ans = UnaryExpr() - ans.op = Operator.fabs - ans.operatorIndex = Operator.operatorIndexDic[self.op] - ans.children.append(self) - return ans + return UnaryExpr(Operator.fabs, self) def __add__(self, other): left = buildGenExprObj(self) @@ -528,9 +520,8 @@ cdef class GenExpr: def __div__(self, other): divisor = buildGenExprObj(other) # we can't divide by 0 - if divisor.getOp() == Operator.const: - assert divisor.number != 0.0 - + if divisor.getOp() == Operator.const and divisor.number == 0.0: + raise ValueError("cannot divide by 0") return self * divisor**(-1) def __rdiv__(self, other): @@ -539,13 +530,17 @@ cdef class GenExpr: return otherexpr.__div__(self) def __truediv__(self,other): - selfexpr = buildGenExprObj(self) - return selfexpr.__div__(other) + print("truediv of GenExpr", self, other) + divisor = buildGenExprObj(other) + # we can't divide by 0 + if divisor.getOp() == Operator.const and divisor.number == 0.0: + raise ValueError("cannot divide by 0") + return self * divisor**(-1) def __rtruediv__(self, other): ''' other / self ''' otherexpr = buildGenExprObj(other) - return otherexpr.__div__(self) + return otherexpr.__truediv__(self) def __neg__(self): return -1.0 * self @@ -622,8 +617,11 @@ cdef class PowExpr(GenExpr): # Exp, Log, Sqrt Expressions cdef class UnaryExpr(GenExpr): - def __init__(self): + def __init__(self, op, expr): self.children = [] + self.children.append(expr) + self.op = op + self.operatorIndex = Operator.operatorIndexDic[op] def __repr__(self): return self.op + "(" + self.children[0].__repr__() + ")" @@ -639,23 +637,11 @@ cdef class Constant(GenExpr): return str(self.number) def exp(expr): - ans = UnaryExpr() - ans.op = Operator.exp - ans.operatorIndex = Operator.operatorIndexDic[ans.op] - ans.children.append(buildGenExprObj(expr)) - return ans + return UnaryExpr(Operator.exp, buildGenExprObj(expr)) def log(expr): - ans = UnaryExpr() - ans.op = Operator.log - ans.operatorIndex = Operator.operatorIndexDic[ans.op] - ans.children.append(buildGenExprObj(expr)) - return ans + return UnaryExpr(Operator.log, buildGenExprObj(expr)) def sqrt(expr): - ans = UnaryExpr() - ans.op = Operator.sqrt - ans.operatorIndex = Operator.operatorIndexDic[ans.op] - ans.children.append(buildGenExprObj(expr)) - return ans + return UnaryExpr(Operator.sqrt, buildGenExprObj(expr)) # transforms tree to an array of nodes. each node is an operator and the position of the # children of that operator (i.e. the other nodes) in the array diff --git a/tests/test_expr.py b/tests/test_expr.py index 5539b24f0..a9dca9698 100644 --- a/tests/test_expr.py +++ b/tests/test_expr.py @@ -51,6 +51,9 @@ def test_upgrade(model): assert isinstance(log(expr), GenExpr) assert isinstance(exp(expr), GenExpr) + with pytest.raises(ValueError): + expr /= 0.0 + def test_genexpr_op_expr(model): m, x, y, z = model genexpr = x**1.5 + y