source: coopr.pyomo/stable/coopr/pyomo/expr/canonical_repn.py @ 3285

Last change on this file since 3285 was 3285, checked in by wehart, 10 years ago

Merged revisions 3184-3284 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pyomo/trunk

........

r3188 | jwatson | 2010-10-29 08:09:35 -0600 (Fri, 29 Oct 2010) | 3 lines


Eliminating initial domain check when constructing numeric constants and the default domain (Reals) is specified.

........

r3203 | jwatson | 2010-10-29 14:46:18 -0600 (Fri, 29 Oct 2010) | 3 lines


Fixing bugs in has_discrete_variables() method of PyomoModel?.

........

r3211 | jdsiiro | 2010-11-01 14:49:11 -0600 (Mon, 01 Nov 2010) | 4 lines


bugfixes for Blocks:

  • avoid infinite loop when adding a block to a model
  • support pretty printing of user-defined components

........

r3212 | jdsiiro | 2010-11-02 15:17:13 -0600 (Tue, 02 Nov 2010) | 5 lines


  • cleaning up the management of Block._parent_block and Component.model attributes. Adding & removing blocks now updates the model attribute on all children
  • renaming Block._setattr_exec -> Block._add_component

........

r3213 | jdsiiro | 2010-11-03 10:47:06 -0600 (Wed, 03 Nov 2010) | 2 lines


Bugfix to the PyomoLogHandler? for python 2.4 compatibility

........

r3214 | jdsiiro | 2010-11-03 15:20:47 -0600 (Wed, 03 Nov 2010) | 3 lines


There is no point logging a warning when the problem is encountered
generating a logging.info message: log the warning at the info level.

........

r3215 | wehart | 2010-11-04 23:05:17 -0600 (Thu, 04 Nov 2010) | 3 lines


Adding a simple knapsack example to illustrate the difference between
a concrete and abstract model.

........

r3216 | jwatson | 2010-11-05 09:39:19 -0600 (Fri, 05 Nov 2010) | 3 lines


Fixing error diagnostic when indexing a variable with a bad index.

........

r3219 | jwatson | 2010-11-05 16:01:23 -0600 (Fri, 05 Nov 2010) | 3 lines


Supressing a validation test with NumericConstant?. If the user specifies a value, we are (now) assuming it is actually a numeric value - otherwise, the domain check significantly inflates the run-time associated with expression tree creation. This needs to be revisited in the Coopr 2.5 re-write.

........

r3226 | wehart | 2010-11-06 21:32:59 -0600 (Sat, 06 Nov 2010) | 2 lines


Setting up example, which was never converted.

........

r3233 | wehart | 2010-11-12 15:56:28 -0700 (Fri, 12 Nov 2010) | 4 lines


Migrating OS-specific functionality into coopr.os


Adding coopr.os to the dev.ini config file.

........

r3242 | wehart | 2010-11-13 01:28:57 -0700 (Sat, 13 Nov 2010) | 4 lines


Type fix.


Updating error message.

........

r3244 | wehart | 2010-11-13 10:44:26 -0700 (Sat, 13 Nov 2010) | 2 lines


Skipping OSiL writer when not defined.

........

r3246 | wehart | 2010-11-13 10:54:15 -0700 (Sat, 13 Nov 2010) | 2 lines


bug fix.

........

r3247 | wehart | 2010-11-13 11:14:01 -0700 (Sat, 13 Nov 2010) | 2 lines


Bug fix.

........

r3248 | jwatson | 2010-11-17 13:51:47 -0700 (Wed, 17 Nov 2010) | 3 lines


Interim fixes to output of quadratic terms in LP writer - more to do, but at least the basic examples now work.

........

r3254 | jwatson | 2010-11-19 13:19:19 -0700 (Fri, 19 Nov 2010) | 3 lines


Fixed bug in LP writer involving quadratic terms involving two distinct variables. Added two new quadratic examples.

........

r3257 | jwatson | 2010-11-19 13:59:35 -0700 (Fri, 19 Nov 2010) | 3 lines


Fixing diagnostic error message when attempting to solve quadratic programs with GLPK - code for generating message was not syntatically legal.

........

r3268 | jwatson | 2010-12-01 15:08:28 -0700 (Wed, 01 Dec 2010) | 3 lines


Fixing issues with the Piecewise construct when breakpoints and slopes are generated via rules. Works now (on a sample of size 1 - the newly added example5.py) for non-indexed rules, likely broken for indexed breakpoint/slope rules.

........

r3272 | jwatson | 2010-12-02 13:53:51 -0700 (Thu, 02 Dec 2010) | 3 lines


Adding omitted pprint() method for SOS constraints - identified while debugging a piecewise issue.

........

r3274 | jwatson | 2010-12-02 16:32:29 -0700 (Thu, 02 Dec 2010) | 3 lines


Adding example of Piecewise construct using breakpoint and slope rules, as opposed to explicit/direct lists.

........

r3276 | jwatson | 2010-12-03 14:06:40 -0700 (Fri, 03 Dec 2010) | 3 lines


Some progress toward functional indexed Piecewise components.

........

File size: 10.4 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the FAST README.txt file.
9#  _________________________________________________________________________
10
11__all__ = ['generate_canonical_repn', 'as_expr', 'is_constant', 'is_linear', 'is_quadratic', 'is_nonlinear']
12
13import logging
14
15#import pyutilib.component.core
16from coopr.pyomo.base import IPyomoPresolver, IPyomoPresolveAction, Model, \
17                             Constraint, Objective
18from coopr.pyomo.base import expr
19from coopr.pyomo.base.var import _VarValue, Var
20import copy
21
22
23Logger = logging.getLogger('Pyomo')
24#
25# A frozen dictionary that can be hashed.  This dictionary isn't _really_
26# frozen, but it acts hashable.
27#
28class frozendict(dict):
29    __slots__ = ('_hash',)
30    def __hash__(self):
31        rval = getattr(self, '_hash', None)
32        if rval is None:
33            rval = self._hash = hash(frozenset(self.iteritems()))
34        return rval
35
36#
37# Generate a canonical representation of an expression.
38#
39# The canonical representation is a dictionary.  Each element is a mapping
40# from a term degree to terms in the expression.  If the term degree is
41# None, then the map value is simply an expression of terms without a
42# specific degree.  Otherwise, the map value is a dictionary of terms to
43# their coefficients.  A term is represented as a frozen dictionary that
44# maps variable id to variable power.  A constant term is represented
45# with None.
46#
47# Examples:
48#  Let x[1] ... x[4] be the first 4 variables, and
49#      y[1] ... y[4] be the next 4 variables
50#
51# 1.3                           {0:{ None :1.3}}
52# 3.2*x[1]                      {1:{ {0:1} :3.2}}
53# 2*x[1]*y[2] + 3*x[2] + 4      {0:{None:4.0}, 1:{{1:1}:3.0}, 2:{{0:1, 5:1}:2.0}}
54# 2*x[1]**4*y[2]    + 4         {0:{None:4.0}, 1:{{1:1}:3.0}, 5:{{0:4, 5:1 }:2.0}}
55# log(y[1]) + x[1]*x[1]         {2:{{0:2}:1.0}, None:log(y[1])}
56#
57def generate_canonical_repn(expr, **kwargs):
58    context = kwargs.get('context', None)
59    key     = kwargs.get('key', None)
60
61    # Prettify output of user notification "What's it doin' now?!"
62    if ( context is not None ):
63        name = " '%s" % context.name
64        if ( isinstance( context, Constraint ) ):
65            item = 'Constraint'
66        elif ( isinstance( context, Objective ) ):
67            item = 'Objective'
68    else:
69        item = 'unknown object'
70        name = ''
71        Logger.info("Argument 'context' not passed")
72
73    if ( key is not None ):
74        index = str(key)
75        if ( type(key) is tuple ):
76            index = index[1:-1] # remove tuple
77        index = "[%s]'" % index
78    else:
79        index = "'"
80        if ( '' == name ):
81            index = ''
82
83    Logger.info("Generating expression for %(item)s%(name)s%(index)s" %
84        { 'item' : item, 'name' : name, 'index' : index } )
85
86    return collect_canonical_repn(expr)
87
88
89def as_expr(rep, vars=None, model=None, ignore_other=False):
90    """ Convert a canonical representation into an expression. """
91    if isinstance(model, Model):
92        vars = model._var
93        id_offset=0
94    elif not vars is None:
95        id_offset=1
96    expr = 0.0
97    for d in rep:
98        if d is None:
99            if not ignore_other:
100                expr += rep[d]
101            continue
102        if d is -1:
103            continue
104        for v in rep[d]:
105            if v is None:
106                expr += rep[d][v]
107                continue
108            e = rep[d][v]
109            for id in v:
110                for i in xrange(v[id]):
111                    if vars is None:
112                        e *= rep[-1][id]
113                    else:
114                        e *= vars[id[1]+id_offset]
115            expr += e
116    return expr
117   
118def repn_add(lhs, rhs, coef=1.0):
119    """
120    lhs and rhs are the expressions being added together.
121    'lhs' and 'rhs' are left-hand and right-hand side operands
122    See generate_canonical_repn for explanation of pyomo expressions
123    """
124    for order in rhs:
125        # For each order term, first-order might be 3*x,
126        # second order 4*x*y or 5*x**2
127        if order is None:
128            # i.e., (currently) order is a constant or logarithm
129            if order in lhs:
130                lhs[order] += rhs[order]
131            else:
132                lhs[order] = rhs[order]
133            continue
134        if order < 0:
135            # ignore now, handled below
136            continue
137        if not order in lhs:
138            lhs[order] = {}
139        for var in rhs[order]:
140            # Add coefficients of variables in this order (e.g., third power)
141            lhs[order][var] = coef*rhs[order][var] + lhs[order].get(var,0.0)
142    #
143    # Merge the VarValue maps
144    #
145    if -1 in rhs:
146        if -1 in lhs:
147            lhs[-1].update(rhs[-1])
148        else:
149            lhs[-1] = rhs[-1]
150    return lhs
151
152
153def repn_mult(r1, r2, coef=1.0):
154    rep = {}
155    for d1 in r1:
156        for d2 in r2:
157            if d1 == None or d2 == None or d1 < 0 or d2 < 0:
158                pass
159            else:
160                d=d1+d2
161                if not d in rep:
162                    rep[d] = {}
163                if d == 0:
164                        rep[d][None] = coef * r1[0][None] * r2[0][None]
165                elif d1 == 0:
166                        for v2 in r2[d2]:
167                            rep[d][v2]  = coef * r1[0][None] * r2[d2][v2] + rep[d].get(v2,0.0)
168                elif d2 == 0:
169                        for v1 in r1[d1]:
170                            rep[d][v1]  = coef * r1[d1][v1] * r2[0][None] + rep[d].get(v1,0.0)
171                else:
172                        for v1 in r1[d1]:
173                            for v2 in r2[d2]:
174                                v = frozendict(v1)
175                                for id in v2:
176                                    if id in v:
177                                        v[id] += v2[id]
178                                    else:
179                                        v[id]  = v2[id]
180                                rep[d][v] = coef * r1[d1][v1] * r2[d2][v2] + rep[d].get(v,0.0)
181    #
182    # Handle other nonlinear terms
183    #
184    if None in r1:
185        rep[None] = as_expr(r2, ignore_other=True) * copy.deepcopy(r1[None])
186        if None in r2:
187            rep[None] += copy.deepcopy(r1[None])*copy.deepcopy(r2[None])
188    if None in r2:
189        if None in rep:
190            rep[None] += as_expr(r1, ignore_other=True) * copy.deepcopy(r2[None])
191        else:
192            rep[None]  = as_expr(r1, ignore_other=True) * copy.deepcopy(r2[None])
193    #
194    # Merge the VarValue maps
195    #
196    if -1 in r1:
197        rep[-1] = r1[-1]
198    if -1 in r2:
199        if -1 in rep:
200            rep[-1].update(r2[-1])
201        else:
202            rep[-1] = r2[-1]
203    #
204    # Return the canonical repn
205    #
206    return rep
207
208
209#
210# Temporary canonical expressions
211#
212# temp_const = { 0: {None:0.0} }   
213# temp_var = { 1: {frozendict({None:1}):1.0} }   
214# temp_nonl = { None: None }
215
216#
217# Internal function for collecting canonical representation, which is
218# called recursively.
219#
220def collect_canonical_repn(exp):
221#     global temp_const
222#     global temp_var
223#     global temp_nonl
224    temp_const = { 0: {None:0.0} }   
225    temp_var = { 1: {frozendict({None:1}):1.0} }   
226    temp_nonl = { None: None }   
227    #
228    # Ignore identity expressions
229    #
230    while isinstance(exp, expr._IdentityExpression):
231        exp = exp._args[0]
232    #
233    # Expression
234    #
235    if isinstance(exp,expr.Expression):
236        #
237        # Sum
238        #
239        if isinstance(exp,expr._SumExpression):
240            if exp._const != 0.0:
241                repn = { 0: {None:exp._const} }   
242            else:
243                repn = {}
244            for i in range(len(exp._args)):
245                repn = repn_add(repn, collect_canonical_repn(exp._args[i]), coef=exp._coef[i] )
246            return repn
247        #
248        # Product
249        #
250        elif isinstance(exp,expr._ProductExpression):
251            #
252            # Iterate through the denominator.  If they aren't all constants, then
253            # simply return this expresion.
254            #
255            denom=1.0
256            for e in exp._denominator:
257                if e.fixed_value():
258                    #print 'Z',type(e),e.fixed
259                    denom *= e.value
260                elif e.is_constant():
261                    denom *= e()
262                else:
263                    temp_nonl[None] = exp
264                    return temp_nonl
265                if denom == 0.0:
266                    print "Divide-by-zero error - offending sub-expression:"
267                    e.pprint()
268                    raise ZeroDivisionError                       
269            #
270            # OK, the denominator is a constant.
271            #
272            repn = { 0: {None:exp.coef / denom} }   
273            #print "Y",exp.coef/denom
274            for e in exp._numerator:
275                repn = repn_mult(repn, collect_canonical_repn(e))
276            return repn
277        #
278        # ERROR
279        #
280        else:
281            raise ValueError, "Unsupported expression type: "+str(exp)
282    #
283    # Constant
284    #
285    elif exp.fixed_value():
286        temp_const[0][None] = exp.value
287        return temp_const
288    #
289    # Variable
290    #
291    elif type(exp) is _VarValue or exp.type() is Var:
292        if type(exp) is _VarValue:
293            mid = id(exp.var.model)
294        else:
295            mid = id(exp.model)
296        temp_var = { -1: {(mid,exp.id):exp}, 1: {frozendict({(mid,exp.id):1}):1.0} }   
297        return temp_var
298    #
299    # ERROR
300    #
301    else:
302       raise ValueError, "Unexpected expression type: "+str(exp)
303
304
305
306def is_constant(repn):
307    """Return True if the canonical representation is a constant expression"""
308    for key in repn:
309        if key is None or key > 0:
310            return False
311    return True
312
313def is_nonlinear(repn):
314    """Return True if the canonical representation is a nonlinear expression"""
315    for key in repn:
316        if not key in [-1,0,1]:
317            return True
318    return False
319
320def is_linear(repn):
321    """Return True if the canonical representation is a linear expression."""
322    return (1 in repn) and not is_nonlinear(repn)
323
324def is_quadratic(repn):
325    """Return True if the canonical representation is a quadratic expression."""
326    if not 2 in repn:
327        return False
328    for key in repn:
329        if key is None or key > 2:
330            return False
331    return True
332
Note: See TracBrowser for help on using the repository browser.