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 | |
---|
13 | import logging |
---|
14 | |
---|
15 | #import pyutilib.component.core |
---|
16 | from coopr.pyomo.base import IPyomoPresolver, IPyomoPresolveAction, Model, \ |
---|
17 | Constraint, Objective |
---|
18 | from coopr.pyomo.base import expr |
---|
19 | from coopr.pyomo.base.var import _VarValue, Var |
---|
20 | import copy |
---|
21 | |
---|
22 | |
---|
23 | Logger = logging.getLogger('Pyomo') |
---|
24 | # |
---|
25 | # A frozen dictionary that can be hashed. This dictionary isn't _really_ |
---|
26 | # frozen, but it acts hashable. |
---|
27 | # |
---|
28 | class 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 | # |
---|
57 | def 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 | |
---|
89 | def 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 | |
---|
118 | def 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 | |
---|
153 | def 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 | # |
---|
220 | def 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 | |
---|
306 | def 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 | |
---|
313 | def 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 | |
---|
320 | def 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 | |
---|
324 | def 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 | |
---|