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 pyutilib.component.core |
---|
14 | from coopr.pyomo.base import IPyomoPresolver, IPyomoPresolveAction, Model |
---|
15 | from coopr.pyomo.base import expr |
---|
16 | from coopr.pyomo.base.var import _VarValue, Var |
---|
17 | import copy |
---|
18 | |
---|
19 | |
---|
20 | # |
---|
21 | # A frozen dictionary that can be hashed. This dictionary isn't _really_ |
---|
22 | # frozen, but it acts hashable. |
---|
23 | # |
---|
24 | class frozendict(dict): |
---|
25 | __slots__ = ('_hash',) |
---|
26 | def __hash__(self): |
---|
27 | rval = getattr(self, '_hash', None) |
---|
28 | if rval is None: |
---|
29 | rval = self._hash = hash(frozenset(self.iteritems())) |
---|
30 | return rval |
---|
31 | |
---|
32 | |
---|
33 | |
---|
34 | def as_expr(rep, vars=None, model=None, ignore_other=False): |
---|
35 | """ Convert a canonical representation into an expression. """ |
---|
36 | if isinstance(model, Model): |
---|
37 | vars = model._var |
---|
38 | id_offset=0 |
---|
39 | elif not vars is None: |
---|
40 | id_offset=1 |
---|
41 | expr = 0.0 |
---|
42 | for d in rep: |
---|
43 | if d is None: |
---|
44 | if not ignore_other: |
---|
45 | expr += rep[d] |
---|
46 | continue |
---|
47 | if d is -1: |
---|
48 | continue |
---|
49 | for v in rep[d]: |
---|
50 | if v is None: |
---|
51 | expr += rep[d][v] |
---|
52 | continue |
---|
53 | e = rep[d][v] |
---|
54 | for id in v: |
---|
55 | for i in xrange(v[id]): |
---|
56 | if vars is None: |
---|
57 | e *= rep[-1][id] |
---|
58 | else: |
---|
59 | e *= vars[id[1]+id_offset] |
---|
60 | expr += e |
---|
61 | return expr |
---|
62 | |
---|
63 | |
---|
64 | def repn_add(rep, r2, coef=1.0): |
---|
65 | for d in r2: |
---|
66 | if d is None: |
---|
67 | if d in rep: |
---|
68 | rep[d] += r2[d] |
---|
69 | else: |
---|
70 | rep[d] = r2[d] |
---|
71 | continue |
---|
72 | if d < 0: |
---|
73 | continue |
---|
74 | if not d in rep: |
---|
75 | rep[d] = {} |
---|
76 | for var in r2[d]: |
---|
77 | rep[d][var] = coef*r2[d][var] + rep[d].get(var,0.0) |
---|
78 | # |
---|
79 | # Merge the VarValue maps |
---|
80 | # |
---|
81 | if -1 in r2: |
---|
82 | if -1 in rep: |
---|
83 | rep[-1].update(r2[-1]) |
---|
84 | else: |
---|
85 | rep[-1] = r2[-1] |
---|
86 | return rep |
---|
87 | |
---|
88 | |
---|
89 | def repn_mult(r1, r2, coef=1.0): |
---|
90 | rep = {} |
---|
91 | for d1 in r1: |
---|
92 | for d2 in r2: |
---|
93 | if d1 == None or d2 == None or d1 < 0 or d2 < 0: |
---|
94 | pass |
---|
95 | else: |
---|
96 | d=d1+d2 |
---|
97 | if not d in rep: |
---|
98 | rep[d] = {} |
---|
99 | if d == 0: |
---|
100 | rep[d][None] = coef * r1[0][None] * r2[0][None] |
---|
101 | elif d1 == 0: |
---|
102 | for v2 in r2[d2]: |
---|
103 | rep[d][v2] = coef * r1[0][None] * r2[d2][v2] + rep[d].get(v2,0.0) |
---|
104 | elif d2 == 0: |
---|
105 | for v1 in r1[d1]: |
---|
106 | rep[d][v1] = coef * r1[d1][v1] * r2[0][None] + rep[d].get(v1,0.0) |
---|
107 | else: |
---|
108 | for v1 in r1[d1]: |
---|
109 | for v2 in r2[d2]: |
---|
110 | v = frozendict(v1) |
---|
111 | for id in v2: |
---|
112 | if id in v: |
---|
113 | v[id] += v2[id] |
---|
114 | else: |
---|
115 | v[id] = v2[id] |
---|
116 | rep[d][v] = coef * r1[d1][v1] * r2[d2][v2] + rep[d].get(v,0.0) |
---|
117 | # |
---|
118 | # Handle other nonlinear terms |
---|
119 | # |
---|
120 | if None in r1: |
---|
121 | rep[None] = as_expr(r2, ignore_other=True) * copy.deepcopy(r1[None]) |
---|
122 | if None in r2: |
---|
123 | rep[None] += copy.deepcopy(r1[None])*copy.deepcopy(r2[None]) |
---|
124 | if None in r2: |
---|
125 | if None in rep: |
---|
126 | rep[None] += as_expr(r1, ignore_other=True) * copy.deepcopy(r2[None]) |
---|
127 | else: |
---|
128 | rep[None] = as_expr(r1, ignore_other=True) * copy.deepcopy(r2[None]) |
---|
129 | # |
---|
130 | # Merge the VarValue maps |
---|
131 | # |
---|
132 | if -1 in r1: |
---|
133 | rep[-1] = r1[-1] |
---|
134 | if -1 in r2: |
---|
135 | if -1 in rep: |
---|
136 | rep[-1].update(r2[-1]) |
---|
137 | else: |
---|
138 | rep[-1] = r2[-1] |
---|
139 | # |
---|
140 | # Return the canonical repn |
---|
141 | # |
---|
142 | return rep |
---|
143 | |
---|
144 | |
---|
145 | # |
---|
146 | # Temporary canonical expressions |
---|
147 | # |
---|
148 | temp_const = { 0: {None:0.0} } |
---|
149 | temp_var = { 1: {frozendict({None:1}):1.0} } |
---|
150 | temp_nonl = { None: None } |
---|
151 | |
---|
152 | # |
---|
153 | # Internal function for collecting canonical representation, which is |
---|
154 | # called recursively. |
---|
155 | # |
---|
156 | def collect_canonical_repn(exp): |
---|
157 | global temp_const |
---|
158 | global temp_var |
---|
159 | global temp_nonl |
---|
160 | # |
---|
161 | # Ignore identity expressions |
---|
162 | # |
---|
163 | while isinstance(exp, expr._IdentityExpression): |
---|
164 | exp = exp._args[0] |
---|
165 | # |
---|
166 | # Expression |
---|
167 | # |
---|
168 | if isinstance(exp,expr.Expression): |
---|
169 | # |
---|
170 | # Sum |
---|
171 | # |
---|
172 | if isinstance(exp,expr._SumExpression): |
---|
173 | if exp._const != 0.0: |
---|
174 | repn = { 0: {None:exp._const} } |
---|
175 | else: |
---|
176 | repn = {} |
---|
177 | for i in range(len(exp._args)): |
---|
178 | repn = repn_add(repn, collect_canonical_repn(exp._args[i]), coef=exp._coef[i] ) |
---|
179 | return repn |
---|
180 | # |
---|
181 | # Product |
---|
182 | # |
---|
183 | elif isinstance(exp,expr._ProductExpression): |
---|
184 | # |
---|
185 | # Iterate through the denominator. If they aren't all constants, then |
---|
186 | # simply return this expresion. |
---|
187 | # |
---|
188 | denom=1.0 |
---|
189 | for e in exp._denominator: |
---|
190 | if e.fixed_value(): |
---|
191 | #print 'Z',type(e),e.fixed |
---|
192 | denom *= e.value |
---|
193 | elif e.is_constant(): |
---|
194 | denom *= e() |
---|
195 | else: |
---|
196 | temp_nonl[None] = exp |
---|
197 | return temp_nonl |
---|
198 | if denom == 0.0: |
---|
199 | print "Divide-by-zero error - offending sub-expression:" |
---|
200 | e.pprint() |
---|
201 | raise ZeroDivisionError |
---|
202 | # |
---|
203 | # OK, the denominator is a constant. |
---|
204 | # |
---|
205 | repn = { 0: {None:exp.coef / denom} } |
---|
206 | #print "Y",exp.coef/denom |
---|
207 | for e in exp._numerator: |
---|
208 | repn = repn_mult(repn, collect_canonical_repn(e)) |
---|
209 | return repn |
---|
210 | # |
---|
211 | # ERROR |
---|
212 | # |
---|
213 | else: |
---|
214 | raise ValueError, "Unsupported expression type: "+str(exp) |
---|
215 | # |
---|
216 | # Constant |
---|
217 | # |
---|
218 | elif exp.fixed_value(): |
---|
219 | temp_const[0][None] = exp.value |
---|
220 | return temp_const |
---|
221 | # |
---|
222 | # Variable |
---|
223 | # |
---|
224 | elif type(exp) is _VarValue or exp.type() is Var: |
---|
225 | if type(exp) is _VarValue: |
---|
226 | mid = id(exp.var.model) |
---|
227 | else: |
---|
228 | mid = id(exp.model) |
---|
229 | temp_var = { -1: {(mid,exp.id):exp}, 1: {frozendict({(mid,exp.id):1}):1.0} } |
---|
230 | return temp_var |
---|
231 | # |
---|
232 | # ERROR |
---|
233 | # |
---|
234 | else: |
---|
235 | raise ValueError, "Unexpected expression type: "+str(exp) |
---|
236 | |
---|
237 | |
---|
238 | # |
---|
239 | # Generate a canonical representation of an expression. |
---|
240 | # |
---|
241 | # The canonical representation is a dictionary. Each element is a mapping |
---|
242 | # from a term degree to terms in the expression. If the term degree is |
---|
243 | # None, then the map value is simply an expression of terms without a |
---|
244 | # specific degree. Otherwise, the map value is a dictionary of terms to |
---|
245 | # their coefficients. A term is represented as a frozen dictionary that |
---|
246 | # maps variable id to variable power. A constant term is represented |
---|
247 | # with None. |
---|
248 | # |
---|
249 | # Examples: |
---|
250 | # Let x[1] ... x[4] be the first 4 variables, and |
---|
251 | # y[1] ... y[4] be the next 4 variables |
---|
252 | # |
---|
253 | # 1.3 {0:{ None :1.3}} |
---|
254 | # 3.2*x[1] {1:{ {0:1} :3.2}} |
---|
255 | # 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}} |
---|
256 | # log(y[1]) + x[1]*x[1] {2:{{0:2}:1.0}, None:log(y[1])} |
---|
257 | # |
---|
258 | def generate_canonical_repn(expr): |
---|
259 | return collect_canonical_repn(expr) |
---|
260 | |
---|
261 | def is_constant(repn): |
---|
262 | """Return True if the canonical representation is a constant expression""" |
---|
263 | for key in repn: |
---|
264 | if key is None or key > 0: |
---|
265 | return False |
---|
266 | return True |
---|
267 | |
---|
268 | def is_nonlinear(repn): |
---|
269 | """Return True if the canonical representation is a nonlinear expression""" |
---|
270 | for key in repn: |
---|
271 | if not key in [-1,0,1]: |
---|
272 | return True |
---|
273 | return False |
---|
274 | |
---|
275 | def is_linear(repn): |
---|
276 | """Return True if the canonical representation is a linear expression.""" |
---|
277 | return (1 in repn) and not is_nonlinear(repn) |
---|
278 | |
---|
279 | def is_quadratic(repn): |
---|
280 | """Return True if the canonical representation is a quadratic expression.""" |
---|
281 | if not 2 in repn: |
---|
282 | return False |
---|
283 | for key in repn: |
---|
284 | if key is None or key > 2: |
---|
285 | return False |
---|
286 | return True |
---|
287 | |
---|