source: coopr.pyomo/stable/2.2/coopr/pyomo/expr/canonical_repn.py @ 2230

Last change on this file since 2230 was 2230, checked in by wehart, 11 years ago

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

........

r2209 | jwatson | 2010-01-30 17:01:12 -0700 (Sat, 30 Jan 2010) | 3 lines


Fix to new LP writer to output variable label, instead of name. The name has embedded quotes, which hoses all kinds of things. The label doesn't.

........

r2210 | jwatson | 2010-01-30 20:11:21 -0700 (Sat, 30 Jan 2010) | 3 lines


Additional minor fixings to make new LP writer work with PH - which it now does.

........

r2211 | jwatson | 2010-01-30 20:17:19 -0700 (Sat, 30 Jan 2010) | 3 lines


Suppressing warning for "constant constraint" when writing LP files - this happens all the time in PH and other contexts, where variable fixing is active.

........

r2213 | wehart | 2010-01-30 23:03:26 -0700 (Sat, 30 Jan 2010) | 2 lines


Adding the ability to process multiple *.dat files on the Pyomo command line

........

r2214 | wehart | 2010-01-31 14:44:12 -0700 (Sun, 31 Jan 2010) | 3 lines


Reworking the canonical expression definition to include a
dictionary from variable ID -> _VarValue object.

........

r2217 | wehart | 2010-01-31 21:48:34 -0700 (Sun, 31 Jan 2010) | 2 lines


Update to NL file writer, after recent changes in the canonical representation.

........

r2218 | wehart | 2010-01-31 22:11:59 -0700 (Sun, 31 Jan 2010) | 7 lines


Reworking the symbol_map generated for the NL file.


Fixing the use of this within PyomoModel? when loading data.


Adding logic to specify the 'solver' option when an optimization interface
is specified (e.g. asl:cplexamp).

........

r2226 | wehart | 2010-01-31 22:27:19 -0700 (Sun, 31 Jan 2010) | 2 lines


Removing some debugging IO.

........

r2227 | wehart | 2010-02-01 14:40:08 -0700 (Mon, 01 Feb 2010) | 7 lines


This is a major rework of the ModelData? logic. I've extracted the
actual operations for data loading into different plugins. These
plugins share a common code base for processing the data into a form that
Pyomo can use, but that logic is now separate from the ModelData? class itself.
This will provide a path-forward for cleaning up this code, as well as
modularizing the parsing of AMPL *.dat files.

........

r2228 | jwatson | 2010-02-01 14:49:40 -0700 (Mon, 01 Feb 2010) | 3 lines


Final fix to LP writer to (optionally) account for the models in which a variable originates - the PySP extensive form writer now passes all tests.

........

File size: 8.7 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, 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
64def 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
89def 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#
148temp_const = { 0: {None:0.0} }   
149temp_var = { 1: {frozendict({None:1}):1.0} }   
150temp_nonl = { None: None }   
151
152#
153# Internal function for collecting canonical representation, which is
154# called recursively.
155#
156def 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#
258def generate_canonical_repn(expr):
259    return collect_canonical_repn(expr)
260
261def 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
268def 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
275def 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
279def 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
Note: See TracBrowser for help on using the repository browser.