source: coopr.pyomo/stable/coopr/pyomo/base/var.py @ 3182

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

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

........

r3116 | jwatson | 2010-10-18 14:08:00 -0600 (Mon, 18 Oct 2010) | 3 lines


First of a few commits involving an experimental linear compaction of the expression tree associated with a constraint, in order to save both memory and time.

........

r3117 | jwatson | 2010-10-18 14:09:13 -0600 (Mon, 18 Oct 2010) | 3 lines


Renaming linear representation module.

........

r3119 | jwatson | 2010-10-18 21:06:55 -0600 (Mon, 18 Oct 2010) | 3 lines


More exhaustive additions relating to linear expression encodings. No functional change to pyomo unless the --linearize-expressions options is specified. All tests are passing.

........

r3121 | jwatson | 2010-10-20 12:02:31 -0600 (Wed, 20 Oct 2010) | 3 lines


Fixing white-space issue when writing quadratic terms in LP files. Gurobi 4.0 beta wants more white-space than CPLEX, which is easy enough to accomodate.

........

r3123 | jwatson | 2010-10-20 13:02:06 -0600 (Wed, 20 Oct 2010) | 3 lines


Changing "integer" section label with the more technically correct "general" in the LP file writer. The former is not recognized by Gurobi, and actually does not appear in the CPLEX documentation either. I'm not sure where we got this, but the fix corrects/bypasses the issue entirely.

........

r3130 | wehart | 2010-10-20 20:44:17 -0600 (Wed, 20 Oct 2010) | 2 lines


Update of coopr.pyomo CHANGELOG for the 2.4 release.

........

r3132 | prsteel | 2010-10-21 15:30:40 -0600 (Thu, 21 Oct 2010) | 2 lines


Adding myself as a developer.

........

r3133 | jwatson | 2010-10-21 15:41:58 -0600 (Thu, 21 Oct 2010) | 5 lines


Various updates, mostly relating to linear expression representations debugging. Initial tests indicate very significant speed and memory reductions (at least for PH), with regression testing performed against a few initial benchmarks. No functional change if linear expressions are not enabled, which is the default


On the variable front, I have added output of variable Status in the pprint(), to facilitate debugging.

........

r3136 | khunter | 2010-10-22 10:19:32 -0600 (Fri, 22 Oct 2010) | 2 lines


Add name as a developer.

........

r3140 | jwatson | 2010-10-22 15:47:33 -0600 (Fri, 22 Oct 2010) | 3 lines


When constructing a linear representation, dump the canonical representation on the ConstraintData? class if it exists - these take up a large amount of space as well.

........

r3141 | jwatson | 2010-10-22 16:54:49 -0600 (Fri, 22 Oct 2010) | 3 lines


Added a "nochecking" keyword argument to Params. Intended strictly for use by algorithm developers, it by-passes the call to _valid_indexed_component, which was taking a huge amount of time in various PySP algorithms. Use sparingly and judiciously! May also apply in the near future to other aspects of setitem.

........

r3142 | jwatson | 2010-10-22 17:07:45 -0600 (Fri, 22 Oct 2010) | 3 lines


Suppression of index validation in Param::setitem when nochecking is enabled.

........

r3144 | jwatson | 2010-10-22 22:15:52 -0600 (Fri, 22 Oct 2010) | 3 lines


Update of test output baseline to reflect addition of variable status output in pretty-print.

........

r3149 | wehart | 2010-10-23 22:59:06 -0600 (Sat, 23 Oct 2010) | 3 lines


Changes to allow RangeSet? instances to be defined with
floating point start, stop and step arguments.

........

r3150 | wehart | 2010-10-23 22:59:31 -0600 (Sat, 23 Oct 2010) | 2 lines


Allowing string values to be elements of a Boolean set.

........

r3151 | wehart | 2010-10-23 23:00:47 -0600 (Sat, 23 Oct 2010) | 2 lines


Removing the name_str() method, which I didn't mean to commit.

........

r3153 | jwatson | 2010-10-24 22:14:15 -0600 (Sun, 24 Oct 2010) | 5 lines


Adding logic for LP writer when dealing with linearized expressions involving fixed variables.


Added "has_discrete_variables" method to PyomoModel?, to support query by MIP plugins - which often need this method to determine the type of solution information (e.g., duals) to extract.

........

r3158 | claird | 2010-10-25 12:01:21 -0600 (Mon, 25 Oct 2010) | 1 line


Added exception for multiple objectives - not currently handled by the NL writier

........

r3159 | claird | 2010-10-25 12:03:18 -0600 (Mon, 25 Oct 2010) | 1 line


Commented out second objective function so there are not multiple objectives in the example

........

r3166 | jwatson | 2010-10-26 19:42:59 -0600 (Tue, 26 Oct 2010) | 3 lines


Various performance-related enhancements, including a rework of NumericValue? and NumericConstant? internals.

........

r3167 | jwatson | 2010-10-26 20:05:25 -0600 (Tue, 26 Oct 2010) | 3 lines


Restoring older version of rangeset.py - the latest commit caused all kinds of test failures, for reasons I haven't had time to explore.

........

r3171 | jwatson | 2010-10-28 21:02:40 -0600 (Thu, 28 Oct 2010) | 3 lines


Fixing bug observed by John regarding passing of keyword arguments to the Constraint base class initializer.

........

r3181 | wehart | 2010-10-28 21:45:53 -0600 (Thu, 28 Oct 2010) | 2 lines


Updating changelog

........

File size: 41.6 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 Coopr README.txt file.
9#  _________________________________________________________________________
10
11__all__ = ['Var', 'VarStatus', '_VarBase', "Piecewise"]
12
13from component import Component
14from numvalue import *
15from numvalue import create_name
16from pyutilib.component.core import alias
17import pyomo
18import types
19from set_types import *
20from indexed_component import IndexedComponent
21import pyutilib.math
22import pyutilib.misc
23from pyutilib.enum import Enum
24import sys
25from param import _ParamValue
26from coopr.pyomo.base.util import isfunctor
27
28VarStatus = Enum( 'undefined', 'fixed_by_presolve', 'fixed_by_optimizer',
29                  'unused', 'used', )
30
31class _VarValue(NumericValue):
32    """Holds the numeric value of a variable, along with suffix info"""
33
34    def __init__(self,**kwds):
35        """Constructor"""
36        if kwds.pop('_deep_copying', None):
37            # Hack for Python 2.4 compatibility
38            # Deep copy will copy all items as necessary, so no need to
39            # complete parsing
40            return
41
42        # NOTE: we currently use the name (in the base class), id, and label attributes,
43        #       but don't document their differences anywhere!
44        # NOTE: the "within" keyword takes precedence over the "domain" keyword.
45        NumericValue.__init__(self, \
46                              name=kwds.get('name', None), \
47                              domain=kwds.get('within', kwds.get('domain', Reals)))
48
49        self.var = None # the "parent" variable.
50        self.index = None # the index of this variable within the "parent"
51
52        # we currently have too many names/labels floating around - the
53        # two below and the name in the base numeric value class.
54        self.id = None
55        self.label = None
56
57        self.initial = None
58
59        self.active = True
60
61        # IMPORTANT: in contrast to the initial value, the lower and upper
62        #             bounds of a variable can in principle be either numeric constants,
63        #             parameter values, or expressions - anything is OK, as long as
64        #             invoking () is legal.
65        self.lb = None
66        self.ub = None
67
68        self.fixed = False
69
70        self.status = VarStatus.undefined
71
72        # cached for efficiency purposes - isinstance is not cheap.
73        self._is_binary = isinstance(self.domain, BooleanSet)
74        self._is_integer = isinstance(self.domain, IntegerSet)
75        self._is_continuous = not (self._is_binary or self._is_integer)
76
77    def __str__(self):
78        return self.name
79
80    def activate(self):
81        self.active=True
82
83    def deactivate(self):
84        self.active=False
85
86    def setlb(self, value):
87        # python only has 3 numeric built-in type that we deal with (we ignore
88        # complex numbers).
89        if isinstance(value, (int, long, float)):
90            self.lb = NumericConstant(value=value)
91        elif value is None:
92            self.lb = None
93        elif isinstance(value, _ParamValue):
94            self.lb = value
95        else:
96            msg = "Unknown type '%s' supplied as variable lower bound - "     \
97                  'legal types are numeric constants or parameters'
98            raise ValueError, msg % str( type(value) )
99
100    def setub(self, value):
101        # python only has 3 numeric built-in type that we deal with (we ignore
102        # complex numbers).
103        if isinstance(value, (int, long, float)):
104            self.ub = NumericConstant(value=value)
105        elif value is None:
106            self.ub = None
107        elif isinstance(value, _ParamValue):
108            self.ub = value
109        else:
110            msg = "Unknown type '%s' supplied as variable lower bound - "     \
111                  'legal types are numeric constants or parameters'
112            raise ValueError, msg % str( type(value) )
113
114
115    def fixed_value(self):
116        if self.fixed:
117            return True
118        return False
119
120    def is_constant(self):
121        if self.fixed:
122            return True
123        return False
124
125    def is_binary(self):
126        return self._is_binary
127
128    def is_integer(self):
129        return self._is_integer
130
131    def is_continuous(self):
132        return self._is_continuous
133
134    def simplify(self, model):
135        return self                 #pragma:nocover
136
137    def pprint(self, ostream=None):
138        if ostream is None:
139           ostream = sys.stdout
140        print >>ostream, str(self),
141
142    def _tighten_bounds(self):
143        """
144        Attempts to tighten the lower- and upper-bounds on this variable
145        by using the domain.
146        """
147        try:
148            dbounds = self.domain.bounds()
149        except:
150            # No explicit domain bounds
151            return
152
153        if self.lb is None or \
154               (dbounds[0] is not None and dbounds[0] > self.lb.value):
155            self.setlb(dbounds[0])
156
157        if self.ub is None or \
158               (dbounds[1] is not None and dbounds[1] < self.ub.value):
159            self.setub(dbounds[1])
160
161
162#
163# Variable attributes:
164#
165# Single variable:
166#
167# x = Var()
168# x.fixed = True
169# x.domain = Reals
170# x.name = "x"
171# x.setlb(-1.0)
172# x.setub(1.0)
173# x.initial = 0.0
174# x = -0.4
175#
176# Array of variables:
177#
178# y = Var(set1,set2,...,setn)
179# y.fixed = True (fixes all variables)
180# y.domain = Reals
181# y.name = "y"
182# y[i,j,...,k] (returns value of a given index)
183#
184# y[i,j,...,k].fixed
185# y[i,j,...,k].lb()
186# y[i,j,...,k].ub()
187# y[i,j,...,k].initial
188#
189
190class _VarBase(IndexedComponent):
191    """A numeric variable, which may be defined over a index"""
192
193    """ Constructor
194        Arguments:
195           name         The name of this variable
196           index        The index set that defines the distinct variables.
197                          By default, this is None, indicating that there
198                          is a single variable.
199           domain       A set that defines the type of values that
200                          each parameter must be.
201           default      A set that defines default values for this
202                          variable.
203           bounds       A rule for defining bounds values for this
204                          variable.
205           rule         A rule for setting up this parameter with
206                          existing model data
207    """
208    # TODO: default and rule keywords are not used?  Need to talk to Bill ...?
209    def __init__(self, *args, **kwd):
210        if kwd.pop('_deep_copying', None):
211            # Hack for Python 2.4 compatibility
212            # Deep copy will copy all items as necessary, so no need to
213            # complete parsing
214            return
215
216        tkwd = {'ctype':Var}
217        IndexedComponent.__init__(self, *args, **tkwd)
218
219        # Default keyword values
220        #
221        self._attr_declarations = {}
222        self._varval = {}
223        self._initialize = kwd.pop('initialize', None )
224        self.name   = kwd.pop('name', 'unknown')
225        self.domain = kwd.pop('within', Reals )
226        self.domain = kwd.pop('domain', self.domain )
227        self.bounds = kwd.pop('bounds', None )
228        self.doc    = kwd.pop('doc', None )
229
230        self._binary_keys = []
231        self._continuous_keys = []
232        self._integer_keys = []
233
234        # Check for domain rules
235        if isfunctor(self.domain):
236            self._domain_rule = self.domain
237
238            # Individual variables will be restricted
239            self.domain = None
240        else:
241            self._domain_rule = None
242
243        if kwd:
244            msg = "Var constructor: unknown keyword '%s'"
245            raise ValueError, msg % kwd.keys()[0]
246
247    def as_numeric(self):
248        if None in self._varval:
249            return self._varval[None]
250        return self
251
252    def keys(self):
253        return self._varval.keys()
254
255    def binary_keys(self):
256        """ Returns the keys of all binary variables """
257        return self._binary_keys
258
259    def continuous_keys(self):
260        """ Returns the keys of all continuous variables """
261        return self._continuous_keys
262
263    def integer_keys(self):
264        """ Returns the keys of all integer variables """
265        return self._integer_keys
266
267    def __iter__(self):
268        return self._varval.keys().__iter__()
269
270    def __contains__(self,ndx):
271        return ndx in self._varval
272
273    def dim(self):
274        return self._ndim
275
276    def reset(self):
277        for key in self._varval:
278            if not self._varval[key].initial is None:
279                self._varval[key].set_value( self._varval[key].initial )
280
281    def __len__(self):
282        return len(self._varval)
283
284    def __setitem__(self,ndx,val):
285        #print "HERE",ndx,val, self._valid_value(val,False), self.domain
286        if (None in self._index) and (ndx is not None):
287            # allow for None indexing if this is truly a singleton
288            msg = "Cannot set an array value in singleton variable '%s'"
289            raise KeyError, msg % self.name
290
291        if ndx not in self._index:
292            msg = "Cannot set the value of array variable '%s' with invalid " \
293                  "index '%s'"
294            raise KeyError, msg % ( self.name, str(ndx) )
295
296        if not self._valid_value(val,False):
297            msg = "Cannot set variable '%s[%s]' with invalid value: %s"
298            raise ValueError, msg % ( self.name, str(ndx), str(val) )
299        self._varval[ndx].value = val
300
301    def __getitem__(self,ndx):
302        """This method returns a _VarValue object.  This object can be
303           coerced to a numeric value using the value() function, or using
304           explicity coercion with float().
305        """
306        try:
307            return self._varval[ndx]
308        except KeyError: # thrown if the supplied index is hashable, but not defined.
309            msg = "Unknown index '%s' in variable %s;" % (str(ndx), self.name)
310            if (isinstance(ndx, (tuple, list)) and len(ndx) != self.dim()) or \
311                   (ndx != self.dim()):
312                msg += "    Expecting %i-dimensional indices" % self.dim()
313            else:
314                msg += "    Make sure the correct index sets were used.\n"
315                msg += "    Is the ordering of the indices correct?"
316            raise KeyError, msg
317        except TypeError, msg: # thrown if the supplied index is not hashable
318            msg2 = "Unable to index variable %s using supplied index with " % self.name
319            msg2 += str(msg)
320            raise TypeError, msg2
321     
322
323    # TODO: Convert to 'has_XXX'
324#     def is_binary(self):
325#         return isinstance(self.domain, BooleanSet)
326
327#     def is_integer(self):
328#         return isinstance(self.domain, IntegerSet)
329
330#     def is_continuous(self):
331#         return not (self.is_binary() or self.is_integer())
332
333    def simplify(self, model):
334        return self
335
336    def _add_domain_key(self, ndx, domain):
337        """ Register an index with a specific set of keys """
338        if isinstance(domain, BooleanSet):
339            self._binary_keys.append(ndx)
340        elif isinstance(domain, IntegerSet):
341            self._integer_keys.append(ndx)
342        else:
343            self._continuous_keys.append(ndx)
344
345    def construct(self, data=None):
346        if pyomo.debug("verbose"):      #pragma:nocover
347           print "Constructing Variable, name="+self.name+", from data="+`data`
348        if self._constructed:
349            return
350        self._constructed=True
351        #
352        # Construct _VarValue() objects for all index values
353        #
354        if self._ndim > 0:
355            if type(self._initialize) is dict:
356                self._index = self._initialize.keys()
357
358            for ndx in self._index:
359                # Check for domain rules
360                if self._domain_rule is not None:
361                    domain = self._domain_rule(ndx, self.model)
362                    if isinstance(domain, BooleanSet):
363                        self._binary_keys.append(ndx)
364                    elif isinstance(domain, IntegerSet):
365                        self._integer_keys.append(ndx)
366                    else:
367                        self._continuous_keys.append(ndx)
368                else:
369                    domain = self.domain
370                    self._add_domain_key(ndx, domain)
371
372                self._varval[ndx] = _VarValue(
373                      name = create_name(self.name,ndx),
374                      domain = domain,
375                      )
376
377                self._varval[ndx].var = self
378                self._varval[ndx].index = ndx
379
380        #
381        # Define the _XXX_keys objects if domain isn't a rule;
382        # they were defined individually above in the case of
383        # a rule
384        #
385        if not self._domain_rule is not None:
386            if isinstance(self.domain, BooleanSet):
387                self._binary_keys = self._varval.keys()
388            elif isinstance(self.domain, IntegerSet):
389                self._integer_keys = self._varval.keys()
390            else:
391                self._continuous_keys = self._varval.keys()
392
393        #
394        # Initialize values with a dictionary if provided
395        #
396        if self._initialize is not None and \
397           type(self._initialize) is not types.FunctionType:
398
399           if type(self._initialize) is dict:
400              for key in self._initialize:
401                self._varval[key].value = None
402                self._varval[key].initial = self._initialize[key]
403                self._valid_value(self._initialize[key],True)
404           else:
405              for key in self._index:
406                self._varval[key].value = None
407                self._varval[key].initial = self._initialize
408              self._valid_value(self._initialize,True)
409        #
410        # Initialize values with the _rule function if provided
411        #
412        elif type(self._initialize) is types.FunctionType:
413           for key in self._varval:
414             if isinstance(key,tuple):
415                tmp = list(key)
416             elif key is None:
417                tmp = []
418             else:
419                tmp = [key]
420             tmp.append(self.model)
421             tmp = tuple(tmp)
422             self._varval[key].value = None
423             self._varval[key].initial = self._initialize(*tmp)
424             self._valid_value(self._varval[key].initial,True)
425        #
426        # Initialize bounds with the bounds function if provided
427        #
428        if self.bounds is not None:
429           # bounds are specified via a tuple
430           if type(self.bounds) is tuple:
431             for key in self._varval:
432               (lb, ub) = self.bounds
433               self._varval[key].setlb(lb)
434               self._varval[key].setub(ub)
435           else:
436             # bounds are specified via a function
437             for key in self._varval:
438               if isinstance(key,tuple):
439                  tmp = list(key)
440               elif key is None:
441                  tmp = []
442               else:
443                  tmp = [key]
444               tmp.append(self.model)
445               tmp = tuple(tmp)
446               (lb, ub) = self.bounds(*tmp)
447               self._varval[key].setlb(lb)
448               self._varval[key].setub(ub)
449           for key in self._varval:
450             if self._varval[key].lb is not None and \
451                not pyutilib.math.is_finite(self._varval[key].lb()):
452                self._varval[key].setlb(None)
453             if self._varval[key].ub is not None and \
454                not pyutilib.math.is_finite(self._varval[key].ub()):
455                self._varval[key].setub(None)
456        #
457        # Iterate through all variables, and tighten the bounds based on
458        # the domain bounds information.
459        #
460        # Only done if self.domain is not a rule. If it is, _VarArray level
461        # bounds become meaningless, since the individual _VarElement objects
462        # likely have more restricted domains.
463        #
464        if self._domain_rule is None:
465            dbounds = self.domain.bounds()
466            if not dbounds is None and dbounds != (None,None):
467                for key in self._varval:
468                    if not dbounds[0] is None:
469                        if self._varval[key].lb is None or                     \
470                            dbounds[0] > self._varval[key].lb():
471                            self._varval[key].setlb(dbounds[0])
472
473                    if not dbounds[1] is None:
474                        if self._varval[key].ub is None or                     \
475                            dbounds[1] < self._varval[key].ub():
476                            self._varval[key].setub(dbounds[1])
477        #
478        # Setup declared attributes for all variables
479        #
480        for attr in self._attr_declarations:
481            for key in self._varval:
482                setattr(self._varval[key],attr,self._attr_declarations[attr][0])
483
484    def pprint(self, ostream=None):
485        if ostream is None:
486           ostream = sys.stdout
487        print >>ostream, "  ",self.name,":",
488        print >>ostream, "\tSize="+str(len(self)),
489        if self.domain is not None:
490            print >>ostream, "\tDomain="+self.domain.name
491        else:
492            print >>ostream, "\tDomain=None"
493        if self._index_set is not None:
494            print >>ostream, "\tIndicies: ",
495            for idx in self._index_set:
496               print >>ostream, str(idx.name)+", ",
497            print ""
498        if None in self._varval:
499           print >>ostream, "\tInitial Value : Lower Bound : Upper Bound : "  \
500                            "Current Value: Fixed: Status"
501           lb_value = None
502           if self._varval[None].lb is not None:
503              lb_value = self._varval[None].lb()
504           ub_value = None
505           if self._varval[None].ub is not None:
506              ub_value = self._varval[None].ub()
507
508           print >>ostream, "\t %s : %s : %s : %s : %s : %s" % (
509             str(self._varval[None].initial),
510             str(lb_value),
511             str(ub_value),
512             str(self._varval[None].value),
513             str(self._varval[None].fixed),
514             str(self._varval[None].status)
515           )
516        else:
517           print >>ostream, "\tKey : Initial Value : Lower Bound : "          \
518                            "Upper Bound : Current Value: Fixed: Status"
519           tmp=self._varval.keys()
520           tmp.sort()
521           for key in tmp:
522             initial_val = self._varval[key].initial
523             lb_value = None
524             if self._varval[key].lb is not None:
525                lb_value = self._varval[key].lb()
526             ub_value = None
527             if self._varval[key].ub is not None:
528                ub_value = self._varval[key].ub()
529
530             print >>ostream, "\t%s : %s : %s : %s : %s : %s : %s" % (
531               str(key),
532               str(initial_val),
533               str(lb_value),
534               str(ub_value),
535               str(value(self._varval[key].value)),
536               str(value(self._varval[key].fixed)),
537               str(self._varval[key].status)
538             )
539
540
541    def display(self, prefix="", ostream=None):
542        if ostream is None:
543           ostream = sys.stdout
544        print >>ostream, prefix+"Variable "+self.name,":",
545        print >>ostream, "  Size="+str(len(self)),
546        print >>ostream, "Domain="+self.domain.name
547        if None in self._varval:
548           print >>ostream, "%s  Value=%s" % (
549             prefix,
550             pyutilib.misc.format_io(self._varval[None].value)
551           )
552        else:
553           for key in self._varval:
554             val = self._varval[key].value
555             print >>ostream, prefix+"  "+str(key)+" : "+str(val)
556
557    def declare_attribute(self, name, default=None):
558        """
559        Declare a user-defined attribute.
560        """
561        if name[0] == "_":
562            msg = "Cannot define an attribute that begins with  '_'"
563            raise AttributeError, msg
564        if name in self._attr_declarations:
565            raise AttributeError, "Attribute %s is already defined" % name
566        self._attr_declarations[name] = (default,)
567        #if not default is None:
568            #self._valid_value(default)
569        #
570        # If this variable has been constructed, then
571        # generate this attribute for all variables
572        #
573        if len(self._varval) > 0:
574            for key in self._varval:
575                setattr(self._varval[key],name,default)
576
577
578class _VarElement(_VarBase,_VarValue):
579
580    def __init__(self, *args, **kwd):
581        if kwd.pop('_deep_copying', None):
582            # Hack for Python 2.4 compatibility
583            # Deep copy will copy all items as necessary, so no need to
584            # complete parsing
585            return
586
587        _VarValue.__init__(self, **kwd)
588        _VarBase.__init__(self, *args, **kwd)
589        self._varval[None] = self
590        self._varval[None].var = self
591        self._varval[None].index = None
592
593    def __call__(self, exception=True):
594        if None in self._varval:
595            return self._varval[None].value
596        return None
597
598    def is_constant(self):
599        return _VarValue.is_constant(self)
600
601class _VarArray(_VarBase):
602
603    def __init__(self, *args, **kwd):
604        if kwd.pop('_deep_copying', None):
605            # Hack for Python 2.4 compatibility
606            # Deep copy will copy all items as necessary, so no need to
607            # complete parsing
608            return
609
610        _VarBase.__init__(self, *args, **kwd)
611        self._dummy_val = _VarValue(**kwd)
612
613    def __float__(self):
614        raise TypeError, "Cannot access the value of array variable "+self.name
615
616    def __int__(self):
617        raise TypeError, "Cannot access the value of array variable "+self.name
618
619    def _valid_value(self,value,use_exception=True):
620        #print "VALID",self._dummy_val.domain
621        return self._dummy_val._valid_value(value,use_exception)
622
623    def set_value(self, value):
624        msg = "Cannot specify the value of array variable '%s'"
625        raise ValueError, msg % self.name
626
627    def __str__(self):
628        return self.name
629
630    def construct(self, data=None):
631        _VarBase.construct(self, data)
632
633        for (ndx, var) in self._varval.iteritems():
634            var._tighten_bounds()
635
636
637class Var(Component):
638    """
639    Variable objects that are used to construct Pyomo models.
640    """
641
642    alias("Var", "Decision variables in a model.")
643
644    def __new__(cls, *args, **kwds):
645        if kwds.pop('_deep_copying', None):
646            # Hack for Python 2.4 compatibility
647            # Deep copy will copy all items as necessary, so no need to
648            # complete parsing
649            return
650
651        if args == ():
652            self = _VarElement(*args, **kwds)
653        else:
654            self = _VarArray(*args, **kwds)
655        return self
656
657
658class Piecewise(_VarElement):
659    """
660    A piecewise linear expression.
661
662    Usage:
663
664    model.X = Var()
665
666    # Let breakpoint_rule and slope_rule define a piecewise linear
667    # function f(x) when indexed over sets (...)
668    #
669    # Then model.XF is a variable constrained by model.XF = f(model.X)
670    model.XF = Piecewise(
671                         breakpoint_rule=breakPointRule,
672                         slope_rule=slopeRule,
673                         breakpoint_sets=(...),
674                         slope_sets=(...),
675                         offset=0,
676                         value=model.X,
677                         forcesos2=True
678                         )
679
680    # Evaluate the piecewise expression at 7.0
681    model.C2.eval(7.0)
682
683    Breakpoints refer to the value on the domain that separate the
684    domain into distinct intervals. All intervals are assumed to be
685    closed on the lower bound and open on the upper bound, with the
686    exception of the interval extending to negative infinity on the
687    lower bound. The first slope specified is assumed to be the slope
688    of an interval coming from negative infinity to the first
689    breakpoint; all other slopes are assumed to correspond to the
690    intervals between the current and next breakpoint. The offset
691    keyword is used to indicate the y-intercept of the expression
692    (default is zero).
693
694    Example:
695
696    The concave function
697
698           / 2x + 3, x < 0
699    f(x) = | 3,      0 < x < 1
700           \ 3 - x,  1 < x
701
702    would be modeled as:
703
704    breakpoints = [0,1]
705    slopes = [2, 0, -1]
706    offset = 3
707
708    If the function should have a bounded domain, specify more
709    breakpoints than slopes.
710
711    Example:
712
713    The concave function
714
715           / 2x + 3, -2 < x < 0
716    f(x) = | 3,       0 < x < 1
717           \ 3 - x,   1 < x < 4
718
719    would be modeled as:
720
721    breakpoints = [-2, 0, 1, 4]
722    slopes = [2, 0, -1]
723    offset = 3
724
725    If we model a bounded domain, the variable passed to model.value
726    will be bounded by this domain. Finally, for bounded domains, the
727    offset is the vertical offset of the leftmost breakpoint.
728
729    The forcesos2 keyword is simply used to always generate the sos2
730    transformation, even if the slopes specify a concave/convex function.
731
732    """
733
734    alias("Piecewise", "A variable constrained by a piecewise function.")
735
736    def __init__(self, *args, **kwargs):
737        """
738        Create a new Piecewise variable
739        """
740
741        # ---------------------------------------------------------------------
742        # Process arguments
743
744        # We can either get explicit breakpoints and slopes, or rules
745        # defining them
746        self._breakpoints = kwargs.pop('breakpoints', None)
747        self._slopes = kwargs.pop('slopes', None)
748
749        self._breakpoint_rule = kwargs.pop('breakpoint_rule', None)
750        self._slope_rule = kwargs.pop('slope_rule', None)
751
752        _breakpoint_sets = kwargs.pop('breakpoint_sets', None)
753        _slope_sets = kwargs.pop('slope_sets', None)
754
755        self._forcesos2 = kwargs.pop('forcesos2', False)
756
757        # Domain variable
758        self.variable = kwargs.pop('value', None)
759
760        # y-intercept
761        self.offset = kwargs.pop('offset', 0)
762
763        # Don't pop the name, we want to pass this to superclass constructor
764        tmpname = kwargs.get('name', 'unknown')
765
766        # ---------------------------------------------------------------------
767        # Error checking
768
769        # Make sure a domain variable was passed
770        if self.variable is None:
771            raise TypeError, "Specify a Variable object for the 'value' " \
772                  "attribute for Piecewise object '%s'" % tmpname
773
774        # Make sure repetitive arguments aren't specified
775        if self._breakpoints is not None and self._breakpoint_rule is not None:
776            raise TypeError, "Specify only one of 'breakpoints' and " \
777                  "'breakpointrule' in Piecewise object '%s'" % tmpname
778        if self._slopes is not None and self._slope_rule is not None:
779            raise TypeError, "Specify only one of 'slopes' and " \
780                  "'sloperule' in Piecewise object '%s'" % tmpname
781
782        # Make sure enough arguments were specified
783        if self._breakpoints is None and self._breakpoint_rule is None:
784            raise TypeError, "Specify exactly one of 'breakpoints' or " \
785                  "'breakpointrule' in Piecewise object '%s'" % tmpname
786        if self._slopes is None and self._slope_rule is None:
787            raise TypeError, "Specify exactly one of 'slopes' or " \
788                  "'sloperule' in Piecewise object '%s'" % tmpname
789
790        # ---------------------------------------------------------------------
791        # Initialization
792
793        # Call superclass
794        _VarElement.__init__(self, *args, **kwargs)
795
796        # Create indexed_components for the breakpoints and slopes
797        if _breakpoint_sets is not None:
798            self._breakpoint_component = IndexedComponent(
799                _breakpoint_sets,
800                ctype="IndexedComponent"
801                )
802        else:
803            self._breakpoint_component = None
804        if _slope_sets is not None:
805            self._slope_component = IndexedComponent(
806                _slope_sets,
807                ctype="IndexedComponent"
808                )
809        else:
810            self._slope_component = None
811
812        # If index sets were passed to either breakpoint or slopes,
813        # make sure the appropriate rule is defined
814        if self._breakpoint_component is not None and \
815               self._breakpoint_rule is None:
816            raise TypeError, "Cannot specify 'breakpoint_sets' without " \
817                  "'breakpoint_rule' in Piecewise object '%s'" % tmpname
818        if self._slope_component is not None and \
819               self._slope_rule is None:
820            raise TypeError, "Cannot specify 'slope_sets' without " \
821                  "'slope_rule' in Piecewise object '%s'" % tmpname
822
823
824        # Have we added implicit SOS constraints?
825        self._constraints_added = False
826
827    def construct(self, *args):
828        """
829        Construct the breakpoint and slope expressions
830        """
831        # Construct the break points
832        self._breakvals = list()
833        if self._breakpoints is not None:
834            # Unpack breakpoints from an iterable
835            for exp in self._breakpoints:
836                if exp is not None:
837                    self._breakvals.append(exp)
838        else:
839            # Get breakpoints from a rule
840            if self._breakpoint_component._ndim == 0:
841                # Rule isn't indexed
842                exp = self._breakpoint_rule((self.model,))
843                if exp is not None:
844                    self._breakvals.append(exp)
845            else:
846                # Rule is indexed
847                for index in self._breakpoint_component._index:
848                    # Make index a list, and then append self.model
849                    if isinstance(index, (list, tuple)):
850                        index = list(index)
851                    else:
852                        index = [index]
853                    index.append(self.model)
854
855                    # Get the next breakpoint
856                    exp = self._breakpoint_rule(*index)
857                    if exp is not None:
858                        self._breakvals.append(exp)
859
860        # Construct the slopes
861        self._slopevals = list()
862        if self._slopes is not None:
863            # Unpack slopes from an iterable
864            for exp in self._slopes:
865                if exp is not None:
866                    self._slopevals.append(exp)
867        else:
868            # Get slopes from a rule
869            if self._slope_component._ndim == 0:
870                # Rule isn't indexed
871                exp = self._slope_rule((self.model,))
872                if exp is not None:
873                    self._slopevals.append(exp)
874            else:
875                # Rule is indexed
876                for index in self._slope_component._index:
877                    # Make index a list, and then append self.model
878                    if isinstance(index, (list, tuple)):
879                        index = list(index)
880                    else:
881                        index = [index]
882                    index.append(self.model)
883
884                    # Get the next breakpoint
885                    exp = self._slope_rule(*index)
886                    if exp is not None:
887                        self._slopevals.append(exp)
888
889        # If N breakpoints are specified, then there can be N-1 or N+1
890        # slopes
891        if abs(len(self._slopevals) - len(self._breakvals)) != 1:
892            msg = "Error constructing Piecewise object '%(name)s'; "          \
893                  "%(breaks)i breakpoints and %(slopes)i slopes were "        \
894                  "specified.\n\tGiven %(breaks)i breakpoints, there should " \
895                  "be %(breaks)i-1 or %(breaks)i+1 slopes."
896            raise ValueError, msg % {
897              "breaks" : len(self._breakvals),
898              "slopes" : len(self._slopevals),
899              "name"   : self.name
900            }
901
902        # If they have specified an unbounded problem, add "None" to both
903        # ends of self._breakvals
904        if len(self._breakvals) < len(self._slopevals):
905            self._breakvals.insert(0, None)
906            self._breakvals.append(None)
907
908        self._generate_constraints()
909
910    def _concave_or_convex(self):
911        """
912        Determines if the expressions forms a concave or convex expression.
913        Returns -1 if concave, 1 if convex, 2 if affine, and 0 otherwise.
914        """
915
916        # Must be convex if there are less than three segments
917        if len(self._slopevals) <= 2:
918            if len(self._slopevals) == 1:
919                # Affine
920                return 2
921            elif self._slopevals[0] < self._slopevals[1]:
922                # Convex
923                return 1
924            elif self._slopevals[0] == self._slopevals[1]:
925                # Affine
926                return 2
927            else:
928                # Concave
929                return -1
930
931        # Find first segment that doesn't have the same slope,
932        # get initial direction
933        start = 1
934        while start < len(self._slopevals):
935            if value(self._slopevals[start]) == \
936                   value(self._slopevals[start-1]):
937                start += 1
938            else:
939                direction = value(self._slopevals[start]) > \
940                            value(self._slopevals[start-1])
941                break
942
943        # If we didn't find a change in slope, we're affine
944        if start == len(self._slopevals) - 1:
945            # Affine
946            return 2
947
948        # Make sure the rest of the segments follow the same direction
949        for i in xrange(start, len(self._slopevals)):
950            if not (direction and value(self._slopevals[i]) >= \
951                    value(self._slopevals[i-1])):
952                return 0
953
954        if direction:
955            # Convex
956            return 1
957        else:
958            # Concave
959            return -1
960
961    def is_convex(self):
962        res = self._concave_or_convex()
963        return res == 1 or res == 2
964
965    def is_concave(self):
966        res = self._concave_or_convex()
967        return res == -1 or res == 2
968
969    def is_affine(self):
970        return self._concave_or_convex() == 2
971
972    def is_concave_or_convex(self):
973        return self._concave_or_convex() != 0
974
975    def _generate_constraints(self):
976        """
977        Generate the constraints on the variable. If this is a convex
978        (or concave) expression, then we simply minimize (or maximize)
979        over a number of linear constraints. If not, we introduce SOS2
980        variables to constrain the expression.
981        """
982
983        # Error if imported at module level
984        from constraint import Constraint, SOSConstraint
985        from rangeset import RangeSet
986
987        if self._constraints_added:
988            # We've already made the implicit constraints
989            return
990
991        # The root of all implicit constraint, variable, and set names
992        baseName = str(self.name) + "_implicit_"
993
994        # Make a partial closure that returns a rule
995        # Arguments passed become closures
996        def makeRule(slope, hOffset, vOffset):
997            expr = slope*(self.variable - hOffset) + vOffset
998            if convex:
999                def rule(model):
1000                    return self >= expr
1001            else:
1002                def rule(model):
1003                    return self <= expr
1004            return rule
1005
1006        def makeBoundRule(value, greater=True):
1007            if greater:
1008                def rule(model):
1009                    return self.variable >= value
1010            else:
1011                def rule(model):
1012                    return self.variable <= value
1013            return rule
1014
1015
1016        # Bound model.value if domain is bounded
1017        if (self._breakvals[0] is not None):
1018            self.model.__setattr__(
1019                baseName + "value_lower_bound",
1020                Constraint(rule=makeBoundRule(self._breakvals[0], True))
1021                )
1022        if (self._breakvals[-1] is not None):
1023            self.model.__setattr__(
1024                baseName + "value_upper_bound",
1025                Constraint(rule=makeBoundRule(self._breakvals[-1], False))
1026                )
1027
1028        convex = self.is_convex()
1029        concave = self.is_concave()
1030        if (convex or concave) and (self._forcesos2 is False):
1031            # Maximize (convex) or minimize (concave) over all
1032            # constraints We first find the interval containing zero
1033
1034            # Handle the case where the domain has no lower bound
1035            if self._breakvals[0] is None:
1036                start = 1
1037
1038                # Line segment with a lower domain bound of -inf
1039                self.model.__setattr__(
1040                    baseName + "constraint[%i]" % 0,
1041                    Constraint(rule=makeRule(self._slopevals[0],
1042                                             self._breakvals[1],
1043                                             self.offset
1044                                             ))
1045                    )
1046            else:
1047                start = 0
1048
1049            # Create the rest of the constraints. Walk from the lower
1050            # to upper bounds of the domain, keeping track of the
1051            # vertical offset.
1052            offset = self.offset
1053            i = start
1054            stop = len(self._breakvals)-1
1055            if (self._breakvals[-1] is None):
1056                stop -= 1
1057
1058            while i <= stop:
1059                slope = self._slopevals[i]
1060                # Make constraint
1061                self.model.__setattr__(
1062                    baseName + "constraint[%i]" % i,
1063                    Constraint(rule=makeRule(slope,
1064                                             self._breakvals[i],
1065                                             offset))
1066                    )
1067
1068                # Update offset
1069                if i < stop:
1070                    offset += slope*(self._breakvals[i+1] - self._breakvals[i])
1071
1072                i += 1
1073        else:
1074            # Introduce SOS2 constraints
1075
1076            # IMPORTANT
1077            #
1078            # For now, we assume that the domain is bounded. It will
1079            # fail otherwise. Support for extended domains coming
1080            # soon.
1081
1082            # Hack to artificially bound the problem, if necessary
1083            if self._breakvals[0] is None:
1084                self._breakvals = self._breakvals[1:]
1085                self._slopevals = self._slopevals[1:]
1086            if self._breakvals[-1] is None:
1087                self._breakvals = self._breakvals[:-1]
1088                self._slopevals = self._slopevals[:-1]
1089
1090            # Walk the breakpoints of the function and find the associated
1091            # images
1092            images = [self.offset]
1093            i = 0
1094            stop = len(self._slopevals) - 1
1095            while i <= stop:
1096                images.append(
1097                    images[-1] + self._slopevals[i] * \
1098                    (self._breakvals[i+1] - self._breakvals[i])
1099                    )
1100                i += 1
1101
1102            # Make the new set indexing the variables
1103            setName = baseName + "SOS_index_set"
1104            self.model.__setattr__(
1105                setName,
1106                RangeSet(0, len(self._breakvals)-1)
1107                )
1108
1109            # Make the new variables
1110            varName = baseName + "SOS_variables"
1111            self.model.__setattr__(
1112                varName,
1113                Var(self.model.__getattribute__(setName))
1114                )
1115
1116            # Constrain variables
1117
1118            # It would be nice to just bound all the SOS variables in the
1119            # Var() class, but then if the variable doesn't appear anywhere
1120            # else in the LP file CPLEX fails, so we explicity constrain
1121            # them here
1122
1123            def makeBoundsRule(varName):
1124                def boundsRule(i, model):
1125                    return (0, model.__getattribute__(varName)[i], 1)
1126                return boundsRule
1127
1128            self.model.__setattr__(
1129                baseName + "bounds_constraint",
1130                Constraint(self.model.__getattribute__(setName),
1131                           rule=makeBoundsRule(varName))
1132                )
1133
1134            # The SOS2 variables sum to 1
1135            def makeSumRule(varName, setName):
1136                def sumRule(model):
1137                    sum(model.__getattribute__(varName)[i] for i in \
1138                        model.__getattribute__(setName)) == 1
1139                return sumRule
1140
1141            self.model.__setattr__(
1142                baseName + "sum_constraint",
1143                Constraint(rule=makeSumRule(varName, setName)))
1144
1145            # The variable self is a convex combination of the images
1146
1147            def makeImageConvexCombRule(varName, setName, images):
1148                def imageConvexCombRule(model):
1149                    return self == sum(
1150                        model.__getattribute__(varName)[i] * images[i] \
1151                        for i in model.__getattribute__(setName))
1152
1153                return imageConvexCombRule
1154
1155            self.model.__setattr__(
1156                baseName + "image_convex_combination_constraint",
1157                Constraint(rule=makeImageConvexCombRule(varName,
1158                                                        setName,
1159                                                        images))
1160                )
1161
1162            # The variable model.value is a convex combination of the
1163            # breakpoints
1164
1165            def makeDomainConvexCombRule(varname, setName, breakpoints):
1166                def domainConvexCombRule(model):
1167                    return self.variable == sum(
1168                        model.__getattribute__(varName)[i] * breakpoints[i] \
1169                        for i in model.__getattribute__(setName))
1170
1171                return domainConvexCombRule
1172
1173            self.model.__setattr__(
1174                baseName + "domain_convex_combination_constraint",
1175                Constraint(rule=makeDomainConvexCombRule(varName,
1176                                                         setName,
1177                                                         self._breakvals
1178                                                         )))
1179
1180            # Finally, nonzero auxialliary variables must be adjacent
1181            self.model.__setattr__(
1182                baseName + "SOS2_constraint",
1183                SOSConstraint(var=self.model.__getattribute__(varName),
1184                              set=self.model.__getattribute__(setName),
1185                              sos=2))
1186
1187        # Don't make the constraints again
1188        self._constraints_added = True
Note: See TracBrowser for help on using the repository browser.