source: coopr.pyomo/trunk/coopr/pyomo/presolve/collect_linear_terms.py @ 2163

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

Adding diagnostic output.

File size: 9.5 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
12import pyutilib.plugin.core
13from coopr.pyomo.base import expr, Var, Constraint, Objective
14from coopr.pyomo.base.numvalue import NumericConstant
15from coopr.pyomo.base.var import _VarValue
16from coopr.pyomo.base import IPyomoPresolver, IPyomoPresolveAction
17
18import string
19import StringIO
20
21
22class CollectLinearTerms(pyutilib.plugin.core.SingletonPlugin):
23    """
24    This plugin processes a Pyomo Model instance to collect
25    linear terms.
26    """
27
28    pyutilib.plugin.core.implements(IPyomoPresolveAction)
29
30    def __init__(self, **kwds):
31        kwds['name'] = "collect_linear_terms"
32        pyutilib.plugin.core.Plugin.__init__(self, **kwds)
33
34    def rank(self):
35        return 10000
36
37    #
38    # Identify variables and confirm that the expression is linear
39    #
40    def _Collect2(self, exp, x, scale=1.0):
41        #
42        # Ignore identity expressions
43        #
44        while isinstance(exp,expr._IdentityExpression):
45            exp = exp._args[0]
46        #
47        # Expression
48        #
49        if isinstance(exp,expr.Expression):
50           #
51           # Sum
52           #
53           if isinstance(exp,expr._SumExpression):
54              for i in xrange(len(exp._args)):
55                x = self._Collect2(exp._args[i], x, scale*exp._coef[i])
56              c = exp._const * scale
57              if "0" in x:
58                   xv = x["0"]
59                   x["0"] = (xv[0]+c,xv[1])
60              else:
61                   x["0"] = (c,"0")
62           #
63           # Product
64           #
65           elif isinstance(exp,expr._ProductExpression):
66              c = scale * exp.coef
67              ve = v = "0"
68              for e in exp._numerator:
69                 while isinstance(e,expr._IdentityExpression):
70                        e = e._args[0]
71                 if e.fixed_value():
72                     # at this point, we have checked that the expression has a fixed value.
73                     # however, that check isn't particularly strong, as the "fixedness" of
74                     # something is typically a function of the class, e.g., numeric values or
75                     # parameter values. it is currently possible to specify expressions for
76                     # parameter values, which is something to warn on and completely chokes
77                     # the presolve.
78                     if isinstance(e.value, expr.Expression):
79                        value_as_string = StringIO.StringIO()
80                        e.value.pprint(ostream=value_as_string)
81                        raise AttributeError, "Expression encountered where fixed value (probably a parameter) was expected - offending subexpression="+string.strip(value_as_string.getvalue())
82                     c *= e.value
83                 else:
84                     if not v is "0":
85                         raise ValueError, "Nonlinear expression: found two non-fixed variables in ProductExpression"
86                     if hasattr(e, "label") is False:
87                        expression_string = StringIO.StringIO()
88                        e.pprint(ostream=expression_string)
89                        msg = "No label for expression object="+string.strip(expression_string.getvalue())+" of type="+str(type(e))
90                        raise AttributeError, msg
91                     v = e.label
92                     ve = e
93              for e in exp._denominator:
94                   if e.fixed_value():
95                       c /= e.value
96                   elif e.is_constant():
97                      # This is experimental, based on issues with logically valid
98                      # expressions being flagged as "nonlinear, with variables".
99                      try:
100                         c /= e()
101                      except ZeroDivisionError:
102                         print "Divide-by-zero error - offending sub-expression:"
103                         e.pprint()
104                         raise ZeroDivisionError                       
105                   else:
106                       print "***Error identified with expression=",e
107                       raise ValueError, "Nonlinear expression: found variable in denominator"
108              if v in x:
109                 xv = x[v]
110                 x[v] = (xv[0]+c,xv[1])
111              else:
112                 x[v] = (c,ve)
113           #
114           # ERROR
115           #
116           else:
117              raise ValueError, "Unsupported expression type: "+str(exp)
118        #
119        # Constant
120        #
121        elif exp.fixed_value():
122           c = exp.value * scale
123           if "0" in x:
124              xv = x["0"]
125              x["0"] = (xv[0]+c,xv[1])
126           else:
127              x["0"] = (c,"0")
128        #
129        # Variable
130        #
131        elif isinstance(exp,_VarValue):
132          v = exp.label
133          if v in x:
134               xv = x[v]
135               x[v] = (xv[0] + scale, xv[1])
136          else:
137               x[v] = (scale, exp)
138        #
139        # ERROR
140        #
141        else:
142           raise ValueError, "Unexpected expression type in _Collect2:"+str(exp)
143        return x
144
145    def _Collect3(self, exp):
146        """
147        The highest-level linear term collection process for an expression.
148        Pretty much just a simple filter over _Collect2.
149        """
150        # the 'x' is just a dictonary mapping variable names to (coefficient, VarValue) pairs.
151        x = self._Collect2(exp,{})
152        # the 'y' is just a reduction of 'x' that eliminates entries with coefficients = 0.0.
153        y = {}
154        for i in x:
155           if x[i][0] != 0.:
156              y[i] = x[i]
157        return y
158
159    def presolve(self,model):
160        """
161        The main routine to perform the presolve
162        """
163        #Vars = model.active_components(Var)
164        Con = model.active_components(Constraint)
165        Obj = model.active_components(Objective)
166        #
167        # Collect the linear terms
168        #
169        # Objectives
170        #
171        Obj1 = []
172        for key in Obj.keys():
173            obj = Obj[key]
174            obj._linterm = {}
175            Onz = []
176            nt = 0
177            for ondx in obj._data:
178                if not obj._data[ondx].active:
179                    continue
180                if obj._data[ondx].expr is None:
181                    raise ValueError, "No expression has been defined for objective %s" % str(key)
182                try:
183                    t = obj._linterm[ondx] = self._Collect3(obj._data[ondx].expr)
184                except Exception, e:
185                    print "Error collecting linear terms for constraint %s (index %s)" % (str(key), str(ondx))
186                    raise e
187                lt = len(t)
188                if lt > 0:
189                    Onz.append(ondx)
190                    nt += 1
191            if nt > 0:
192                Obj1.append(key)
193            obj._Onz = Onz
194        model.Onontriv = Obj1
195        #
196        # Constraints
197        #
198
199        # there are currently expressions within pyomo that the collect routines
200        # have issues with - most notably, _SumExpression objects not having labels
201        # (which of course they don't). we are working to fix this issue, but in
202        # the meantime, we are providing a bit more context regarding which specific
203        # constraint is causing the issue.
204
205        # a model maintains a list of constraints with at least one non-empty expression index.
206        # attribute is "Cnontriv" - see set below.
207        Con1 = [] 
208        for key in Con.keys():
209           C = Con[key]
210           # for each constraint, a map between each index and the associated linear term construct.             
211           C._linterm = {}
212           # each constraint tracks a list of the indices with non-empty expressions.
213           # attribute is _Cnz - see set below.             
214           Cnz = [] 
215           nt = 0 # track number of constraint indicies with non-trivial bodies
216           for cndx in C.keys():
217              if not C._data[cndx].active:
218                    continue
219                 
220#              print "Collecting linear terms for constraint="+key+"["+str(cndx)+"]"
221#              junk = StringIO.StringIO()
222#              C._data[cndx].body.pprint(ostream=junk)                 
223#              print "CONSTRAINT BODY:"+string.strip(junk.getvalue())
224
225              try:                 
226                 t = C._linterm[cndx] = self._Collect3(C._data[cndx].body)
227              except AttributeError, msg:
228                 print "***Error encountered in class CollectLinearTerms, method=presolve, encountered when invoking _Collect3 on constraint="+key+"["+str(cndx)+"]"
229                 constraint_as_string = StringIO.StringIO()
230                 C._data[cndx].body.pprint(ostream=constraint_as_string)                                 
231                 print "***Constraint body:"+string.strip(constraint_as_string.getvalue())
232                 print "***Problem: "+str(msg)
233                 raise RuntimeError
234                 
235              if len(t) > 0:
236                 Cnz.append(cndx)
237                 nt += 1
238                   
239           if nt > 0:
240              Con1.append(key)
241           C._Cnz = Cnz
242        model.Cnontriv = Con1
243
244        #
245        # Return modified instance
246        #
247        return model
248
Note: See TracBrowser for help on using the repository browser.