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