source: coopr.pyomo/trunk/coopr/pyomo/expr/canonical_repn.py @ 2201

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

Update to Coopr to account for changes in PyUtilib? package names.

File size: 8.1 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
13#import pyutilib.component.core
14from coopr.pyomo.base import IPyomoPresolver, IPyomoPresolveAction, Model
15from coopr.pyomo.base import expr
16from coopr.pyomo.base.var import _VarValue, Var
17import copy
18
19
20#
21# A frozen dictionary that can be hashed.  This dictionary isn't _really_
22# frozen, but it acts hashable.
23#
24class 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
34def 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
60def 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
75def 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#
128temp_const = { 0: {None:0.0} }   
129temp_var = { 1: {frozendict({None:1}):1.0} }   
130temp_nonl = { None: None }   
131
132#
133# Internal function for collecting canonical representation, which is
134# called recursively.
135#
136def 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#
234def generate_canonical_repn(expr, model):
235    return collect_canonical_repn(expr, model)
236
237def is_constant(repn):
238    """Return True if the canonical representation is a constant expression"""
239    return (0 in repn) and (len(repn) == 1)
240
241def 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
248def 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
252def 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
Note: See TracBrowser for help on using the repository browser.