source: coopr.pyomo/trunk/coopr/pyomo/base/piecewise.py @ 5599

Revision 5599, 67.2 KB checked in by gabeh, 2 years ago (diff)

A few more updates to clean up code in the piecewise component.

Line 
1
2__all__ = ['Piecewise']
3
4from component import Component
5from constraint import Constraint, SOSConstraint, ConstraintBase
6from var import Var, _VarArray, _VarElement
7from sets import Set, _BaseSet
8from set_types import *
9from pyutilib.component.core import alias
10from pyutilib.misc import flatten_tuple
11import math
12import copy
13from itertools import product
14
15"""
16This file contains a library of functions needed to construct linear/piecewise-linear
17constraints for a Pyomo model. All piecewise types except for SOS2, BIGM_SOS1,
18BIGM_BIN were taken from the paper: Mixed-Integer Models for Non-separable Piecewise Linear
19Optimization: Unifying framework and Extensions (Vielma, Nemhauser 2008).
20
21USAGE NOTES:
22*) The BIGM methods currently determine tightest constant M values. This method is
23   implemented in such a way that binary/SOS variables are not created when this M is
24   zero. In a way, this determines convexity/concavity and can simplify constraints up to the point
25   where no binary/SOS variables may exist. The other piecewise representations have the ability to do
26   this simplification only with the check for convexity/concavity. *** not quite ready, I found a
27   test case that fails *** Note, however,for strictly affine functions the BIGM
28   method can still eliminate a few variables when determining the tightest M value.
29*) I have NOT tested what happens when Pyomo params are in f_rule/d_f_rule. I think the only safe thing to do
30   for now would be to place value() around the params if they are present. The function rules need to return a
31   numeric value, I'm not sure what happens when they don't.
32
33
34TODOS (in rough order of priority):
35_DONE_ *) insure points are sorted ???? THIS IS NEEDED FOR SOME PW_REPS
36_DONE_ *) user not providing floats can be an major issue for BIGM's and MC
37_DONE_ *) raise error when (2^n)+1 break points not givenfor LOG, DLOG piecewise reps
38_DONE_ *) check for convexity/concavity, do we really need the convex/concave kwds? ---- replaced convex/concave kwds with force_pw
39       *) check for continuity?? This is required, but should we check?
40       *) Consider another piecewise rep "SOS2_MANUAL" where we manually implement extra constraints to
41          define an SOS2 set, this would be compatible with GLPK, http://winglpk.sourceforge.net/media/glpk-sos2_02.pdf
42_DONE_ *) arbitrary number index sets of arbitrary dimension
43_DONE_  *) params in f_rules? (in particular when the they don't have value() around them)
44_DONE_ *) non-indexed variables
45       *) Clean up BIGM code if convexity check is implemented
46       *) handle when break points are given as Sets??? Right now must be dictionary( or list for non-indexed)
47       *) write tests
48_DONE_ *) check when adding elements to dicts, those elements don't already exist in model, hasattr()
49_DONE_ *) finish function documentation string
50_DONE_ *) add more description of capabilities
51_DONE_ *) clean up naming scheme for added model attributes
52_DONE_ *) remove addition of optimal BIGM values as model attribute
53_DONE_ *) tests that inputs make sense
54_DONE_ *) collect list of variables added to the model like constraints
55_DONE_ *) collect list of sets added to the model just like constraints
56       *) affine functions - BIGM_SOS1, BIGM_SOS2 *****Not quite done
57*) speed up? Search for BOTTLENECK to see a few spots worth considering
58*) double check that LOG and DLOG reps really do require 2^n points
59*) multivar linearizations
60*) multivar piecewise reps
61"""
62
63def isPowerOfTwo(x):
64    if (x <= 0):
65        return false
66    else:
67        return ( (x & (x - 1)) == 0 )
68
69def isIncreasing(l):
70    """
71    checks that list of points is
72    strictly increasing
73    """
74    for i, el in enumerate(l[1:]):
75        if el <= l[i]:
76            return False
77    return True
78
79def isDecreasing(l):
80    """
81    checks that list of points is
82    strictly decreasing
83    """
84    for i, el in enumerate(l[1:]):
85        if el >= l[i]:
86            return False
87    return True
88
89def isNonDecreasing(l):
90    """
91    checks that list of points is
92    not decreasing
93    """
94    for i, el in enumerate(l[1:]):
95        if el < l[i]:
96            return False
97    return True
98
99def isNonIncreasing(l):
100    """
101    checks that list of points is
102    not increasing
103    """
104    for i, el in enumerate(l[1:]):
105        if el > l[i]:
106            return False
107    return True
108
109def isIndexer(x):
110    """
111    Determines if object is eligable to
112    index the method
113    """
114    return isinstance(x, set) or \
115           isinstance(x,_BaseSet) or \
116           isinstance(x,list) or \
117           isinstance(x,tuple)
118
119class TupleIterator:
120    def __init__(self, underlyingIter):
121        self.underlying = underlyingIter
122
123    def __iter__(self):
124        return self
125
126    def next(self):
127        try:
128            return (self.underlying.next(),)
129        except StopIteration:
130            raise
131
132class TupleSet(set):
133    """
134    Behaves as a set but has
135    tuplize() method which
136    returns elements as a singleton
137    tuple
138    """
139   
140    def __init__(self, dimen, *args):
141        self.dimen = dimen
142        set.__init__(self,*args)
143
144    def tuplize(self):
145        if self.dimen > 1:
146            return self.__iter__()
147        elif self.dimen == 1:
148            return TupleIterator(self.__iter__())
149        else:
150            return [()]
151
152def _GrayCode(nbits):
153
154    bitset = [0 for i in xrange(nbits)]
155    graycode = [''.join(map(str,bitset))]
156
157    for i in xrange(2,(1<<nbits)+1):
158        if i%2:
159            for j in xrange(-1,-nbits,-1):
160                if bitset[j]:
161                    bitset[j-1]=bitset[j-1]^1
162                    break
163        else:
164            bitset[-1]=bitset[-1]^1
165
166        graycode.append(''.join(map(str,bitset)))
167
168    return graycode
169
170
171def _DLog_Branching_Scheme(L):
172    """
173    Branching scheme for DLOG
174    """
175
176    MAX = 2**L
177    mylists1 = {}
178    for i in xrange(1,L+1):
179        mylists1[i] = []
180        start = 1
181        step = MAX/(2**i)
182        while(start < MAX):
183            mylists1[i].extend([j for j in xrange(start,start+step)])
184            start += 2*step
185       
186    biglist = [i for i in xrange(1,MAX+1)]
187    mylists2 = {}
188    for i in sorted(mylists1.keys()):
189        mylists2[i] = []
190        for j in biglist:
191            if j not in mylists1[i]:
192                mylists2[i].append(j)
193        mylists2[i] = sorted(mylists2[i])
194   
195    return mylists1, mylists2
196
197def _Log_Branching_Scheme(n):
198    """
199    Branching scheme for LOG, requires G be a gray code
200    """
201
202    BIGL = 2**n
203    S = [j for j in xrange(1,n+1)]   
204    code = _GrayCode(n)
205    G = dict((i, [int(j) for j in code[i-1]]) for i in xrange(1,len(code)+1))
206
207    L = {}
208    R = {}
209    for s in S:
210        L[s] = []
211        R[s] = []
212        for k in range(BIGL+1):
213            if ((k == 0) or (G[k][s-1] == 1)) and ((k == BIGL) or (G[k+1][s-1] == 1)):
214                L[s].append(k+1)
215            if ((k == 0) or (G[k][s-1] == 0)) and ((k == BIGL) or (G[k+1][s-1] == 0)):
216                R[s].append(k+1)
217               
218    return S,L,R
219
220def characterize_function(f,tpl):
221   
222    # determines convexity or concavity
223    # return 1 (convex), -1 (concave), 0 (neither)
224   
225    points = tpl[-1]
226    new_tpl = flatten_tuple(tpl[:-1])
227    SLOPES = [(f( *(new_tpl+(points[i],)) )-f( *(new_tpl+(points[i-1],)) ))/(points[i]-points[i-1]) for i in xrange(1,len(points))]
228   
229    if isNonDecreasing(SLOPES):
230        # convex
231        return 1
232    elif isNonIncreasing(SLOPES):
233        # concave
234        return -1
235   
236    return 0
237
238
239class Piecewise(Component):
240    """
241    Adds linear/piecewise-linear constraints to a Pyomo model for functions of the form, y = f(x).
242   
243    Usage:
244            model.const = Piecewise(model.index,model.yvar,model.xvar,**Keywords) - for indexed variables
245            model.const = Piecewise(model.yvar,model.xvar,**Keywords) - for non-indexed variables
246           
247            where,
248                    model.index - indexing set
249                    model.yvar - Pyomo Var which represents the range values of f(x)
250                    model.xvar - Pyomo Var which represents the domain values of f(x)
251
252    Keywords:
253   
254-pw_pts={},[],() 
255          A dictionary of lists (keys are index set) or a single list (for non-indexed variables) defining
256          the set of domain vertices for a continuous piecewise linear function. Default is None
257          *** REQUIRES 'f_rule', and 'pw_constr_type' kwds, 'pw_repn' is optional ***
258
259-tan_pts={},[],()
260          A dictionary of lists (keys are index set) or a single list (for non-indexed variables) defining
261          the set of domain vertices for the linear under/over estimators.
262          *** REQUIRES 'f_rule', 'd_f_rule', and 'tan_constr_type' kwds ***
263
264-pw_repn=''
265          Indicates the type of piecewise representation to use. Each representation uses its own sets of different
266          constraints and variables. This can have a major impact on performance of the MIP solve.
267          Choices:
268                 
269                 ~ + 'SOS2' - standard representation using SOS2 variables. (Default)
270                 ~ + 'BIGM_BIN' - uses BigM constraints with binary variables. Theoretically tightest M
271                                  values are automatically determined.
272                 ~ + 'BIGM_SOS1' - uses BigM constraints with sos1 variables. Theoretically tightest M
273                                  values are automatically determined.
274                 ~*+ 'DCC' - Disaggregated convex combination model
275                 ~*+ 'DLOG' - Logarithmic Disaggregated convex combination model
276                 ~*+ 'CC' - Convex combination model
277                 ~*+ 'LOG' - Logarithmic branching convex combination
278                 ~*+ 'MC' - Multiple choice model
279                 ~*+ 'INC' - Incremental (delta) method
280
281                   * Source: "Mixed-Integer Models for Non-separable Piecewise Linear Optimization:
282                              Unifying framework and Extensions" (Vielma, Nemhauser 2008)
283                   ~ Refer to the optional 'force_pw' kwd.
284
285-pw_constr_type=''
286          Indicates whether piecewise function represents an upper bound, lower bound, or equality constraint.
287          *** REQUIRED WHEN 'pw_pts' ARE SPECIFIED ***
288          Choices:
289                 
290                   + 'UB' - y variable is bounded above by piecewise function
291                   + 'LB' - y variable is bounded below by piecewise function
292                   + 'EQ' - y variable is equal to the piecewise function
293
294-tan_constr_type
295          Indicates whether the function linearization represents an upper bound or lower bound.
296          *** REQUIRED WHEN 'tan_pts' ARE SPECIFIED ***
297          Choices:
298                 
299                   + 'UB' - y variable is bounded above by the function linearizations
300                   + 'LB' - y variable is bounded below by the function linearizations
301
302-f_rule=f(model,i,j,...,x)
303          A callable object defining the function that linearizations/piecewise constraints will be applied to.
304          First argument must be a Pyomo model. The last argument is the domain value at which the function evaluates.
305          Intermediate arguments are the corresponding indices of the Piecewise variables (model.yvar, model.xvar).
306          If these variables are non-indexed, no intermediate arguments are needed.
307          *** ALWAYS REQUIRED ***
308          Examples:
309                   + def f(model,j,x):
310                          if (j == 2):
311                             return x**2 + 1.0
312                          else:
313                             return x**2 + 5.0
314                   + f = lambda model,x: return exp(x) + value(model.p) (p is a Pyomo Param)
315                   + def f(model,x):
316                         return mydictionary[x]
317
318-d_f_rule=f(model,i,j,...,x)
319          A callable object defining the derivative of the function that linearizations will be applied to.
320          First argument must be a Pyomo model. The last argument is the domain value at which the function evaluates.
321          Intermediate arguments are the corresponding indices of the Piecewise variables (model.yvar, model.xvar).
322          If these variables are non-indexed, no intermediate arguments are needed.
323          *** REQUIRED WHEN 'tan_pts' ARE SPECIFIED ***
324          Examples:
325                   + def df(model,j,x):
326                          return 2.0*x
327                   + df = lambda model,x: return exp(x)
328                   + def df(model,x):
329                         return mydictionary[x]                     
330
331-force_pw=True/False
332          Using the given function rule and pw_pts, a check for convexity/concavity is implemented. If (1) the function is convex and
333          the piecewise constraints are lower bounds or if (2) the function is concave and the piecewise constraints are upper bounds
334          the piecewise constraints will be substituted for strictly linear constraints. Setting 'force_pw=True' will force the use
335          of the original piecewise constraints when one of these two cases applies.
336    """
337   
338    alias("Piecewise", "A Generator of linear and piecewise-linear constraints.")
339   
340    def __init__(self, *args, **kwds):
341
342        self.piecewise_choices = ['BIGM_SOS1','BIGM_BIN','SOS2','CC','DCC','DLOG','LOG','MC','INC']
343        self.PW_PTS = kwds.pop('pw_pts',None)
344        self.GRAD_PTS = kwds.pop('tan_pts',None)
345        self.PW_REP = kwds.pop('pw_repn', 'SOS2')
346        self.PW_TYPE = kwds.pop('pw_constr_type',None)
347        self.GRAD_TYPE = kwds.pop('tan_constr_type',None)
348        self.candidate_f = kwds.pop('f_rule',None)
349        self.candidate_d_f = kwds.pop('d_f_rule',None)
350        self._force_pw = kwds.pop('force_pw',False)
351        if len(kwds) != 0:
352            raise AssertionError, "WARNING: Possible invalid keywords given to LinearConstraint: "+str(kwds.keys())
353       
354        self._element_mode = False
355        self._index_args = []
356        self._index_dimen = None
357       
358        if len(args) > 2:
359            if all([isIndexer(args[i]) for i in range(len(args)-2)]) and \
360            isinstance(args[-2],_VarArray) and \
361            isinstance(args[-1],_VarArray):
362                self._index_args = [args[i] for i in range(len(args)-2)]
363                self._index_dimen = 0
364                for arg in self._index_args:
365                    if isinstance(arg,_BaseSet):
366                        self._index_dimen += arg.dimen
367                    elif isinstance(arg, set):
368                        element = iter(arg).next()
369                        if isinstance(element,tuple):
370                            self._index_dimen += len(element)
371                        else:
372                            self._index_dimen += 1
373                    elif isinstance(arg,list) or isinstance(arg,tuple):
374                        element = arg[0]
375                        if isinstance(element,tuple):
376                            self._index_dimen += len(element)
377                        else:
378                            self._index_dimen += 1
379                self.YVAR = args[-2]
380                self.XVAR = args[-1]
381                args = ()
382        elif isinstance(args[0],_VarElement) and isinstance(args[1],_VarElement):
383            self.YVAR = args[0]
384            self.XVAR = args[1]
385            self._index_dimen = 0
386            args = ()
387            self._element_mode = True
388        else:
389            raise TypeError, "Piecewise: Invalid function arguments. Syntax: \n\
390            (1) F(Set,Var,Var,**kwds) \n\
391            (2) F(set,Var,Var,**kwds) \n\
392            (3) F(list,Var,Var,**kwds) \n\
393            (4) F(tuple,Var,Var,**kwds) \n\
394            (5) F(Var,Var,**kwds)"
395       
396        self._constraints_dict = {}
397        self._vars_dict = {}
398        self._sets_dict = {}
399       
400        # Test that inputs make sense
401        if self._force_pw not in [True,False]:
402            raise AssertionError, "Piecewise: invalid kwd value for 'force_pw', must be True or False"
403        if (self.GRAD_PTS is not None) and (self.GRAD_TYPE is None):
404            raise AssertionError, "Piecewise kwd 'tan_constr_type' required when 'tan_pts' is specified"
405        if (self.PW_PTS is not None):
406            if (self.PW_TYPE is None):
407                raise AssertionError, "Piecewise kwd 'pw_constr_type' required when 'pw_pts' is specified"
408            if (self.PW_REP not in self.piecewise_choices):
409                raise AssertionError, repr(self.PW_REP)+" is not a valid type for Piecewise kwd 'pw_repn' \n\
410                                      Choices are: "+str(self.piecewise_choices)
411                                                 
412        if (self._element_mode is True) and (not (self.PW_PTS is None)) and (not isinstance(self.PW_PTS, list)) and (not isinstance(self.PW_PTS,tuple)):
413            raise AssertionError, "Piecewise kwd 'pw_pts' must be of type list or tuple when no indexing set is given"
414        if (self._element_mode is False) and (not (self.PW_PTS is None)) and (not isinstance(self.PW_PTS, dict)):
415            raise AssertionError, "Piecewise kwd 'pw_pts' must be of type dict when indexing set is given"
416        if (self._element_mode is True) and (not (self.GRAD_PTS is None)) and (not isinstance(self.GRAD_PTS, list)) and (not isinstance(self.GRAD_PTS,tuple)):
417            raise AssertionError, "Piecewise kwd 'tan_pts' must be of type list or tuple when no indexing set is given"
418        if (self._element_mode is False) and (not (self.GRAD_PTS is None)) and (not isinstance(self.GRAD_PTS, dict)):
419            raise AssertionError, "Piecewise kwd 'tan_pts' must be of type dict when indexing set is given"
420
421        self.name = "unknown"
422        # Construct parent
423        Component.__init__(self,ctype=Piecewise, **kwds)
424        self._index_set = None
425        self._index = None
426       
427       
428    def _grad_pts(self,*targs):
429        args = flatten_tuple(targs)
430        if len(args) > 1:
431            return self.GRAD_PTS[args]
432        elif len(args) == 1:
433            return self.GRAD_PTS[args[0]]
434        else:
435            return self.GRAD_PTS
436       
437    def _pw_pts(self,*targs):
438        args = flatten_tuple(targs)
439        if len(args) > 1:
440            return self.PW_PTS[args]
441        elif len(args) == 1:
442            return self.PW_PTS[args[0]]
443        else:
444            return self.PW_PTS
445       
446    def _getvar(self,model,name,*targs):
447        args = flatten_tuple(targs)
448        if len(args) > 1:
449            return getattr(model,name)[args]
450        elif len(args) == 1:
451            return getattr(model,name)[args[0]]
452        return getattr(model,name)
453       
454    def _getset(self,model,name,*targs):
455        args = flatten_tuple(targs)
456        if len(args) > 1:
457            return getattr(model,name)[args]
458        elif len(args) == 1:
459            return getattr(model,name)[args[0]]
460        return getattr(model,name)
461       
462    def pprint(self, *args, **kwds):
463        pass
464   
465    def reset(self):            #pragma:nocover
466        pass                    #pragma:nocover
467   
468    def deactivate(self):
469        """
470        Deactivates all model constraints created by this class.
471        """
472        [con.deactivate() for con in self._constraints_dict.itervalues()]
473        self.active = False
474   
475    def activate(self):
476        """
477        Activates all model constraints created by this class.
478        """
479        [con.activate() for con in self._constraints_dict.itervalues()]
480        self.active = True
481   
482    def construct(self, *args, **kwds):
483       
484        #A hack to call add after data has been loaded.
485        if self._model()._defer_construction:
486            self._model().concrete_mode()
487            self.add()
488            self._model().symbolic_mode()
489        else:
490            self.add()
491        self._constructed=True
492       
493
494    def add(self, *args):
495
496        #Mimics Constraint.add, but is only used to alert preprocessors
497        #to its existence. Ignores all arguments passed to it.
498        if len(self._index_args) > 1:
499            self.index = set([flatten_tuple(idx) for idx in product(*self._index_args)])
500        elif len(self._index_args) == 1:
501            self.index = self._index_args[0]
502        else:
503            # Non-indexed case
504            self.index = []
505       
506       
507        if self.candidate_f != None:
508            self._define_F(self.candidate_f)
509        if self.candidate_d_f != None:
510            self._define_d_F(self.candidate_d_f)
511       
512        #construct the linear over/under estimators
513        if self.GRAD_PTS != None:
514            self._construct_linear_estimators(TupleSet(self._index_dimen,self.index))
515       
516        #If a piecewise list is given then construct the piecewise constraints.
517        #These functions look at the size of each list of x points and create
518        #variables/constraints as needed. For instance if only two points are supplied
519        #as a list for piecewise, no binary variable is actually needed, therefore
520        #the constraints are created using the _construct_simple_single_bound function.
521        if self.PW_PTS != None:
522           
523            #check that lists are sorted and for special cases
524            #(pw_repn='LOG' or 'DLOG') that lists have 2^n + 1 points.
525            if self._element_mode:
526                if (self.PW_REP in ['DLOG','LOG']) and (not isPowerOfTwo(len(self._pw_pts())-1)):
527                    raise AssertionError, self.name+": 'pw_pts' list must have 2^n + 1 points in list for " \
528                                                       +self.PW_REP+" type piecewise representation"
529                if not isIncreasing(self._pw_pts()):
530                    raise AssertionError, self.name+": 'pw_pts' must be sorted in increasing order"
531            else:
532                if (self.PW_REP in ['DLOG','LOG']) and (not all((isPowerOfTwo(len(self._pw_pts(t))-1) for t in self.index))):
533                    raise AssertionError, self.name+": 'pw_pts' lists must have 2^n + 1 points in each list for " \
534                                                       +self.PW_REP+" type piecewise representation"
535                if not all((isIncreasing(self._pw_pts(t)) for t in self.index)):
536                    raise AssertionError, self.name+": 'pw_pts' must be sorted in increasing order"
537               
538            #determine convexity or concavity in terms of piecewise points
539            simplified_index = []
540            pw_index = []
541            force_simple = False
542            force_pw = False
543            if self._element_mode:
544                tpl = (self._model(),)+(self._pw_pts(),)
545                result = characterize_function(self._F, tpl)
546                if (result == -1):
547                    self.function_character = 'concave'
548                    if (self.PW_TYPE == 'UB'):
549                        force_simple = True
550                    else:
551                        force_pw = True
552                elif (result == 1):
553                    self.function_character = 'convex'
554                    if (self.PW_TYPE == 'LB'):
555                        force_simple = True
556                    else:
557                        force_pw = True
558                else:
559                    self.function_character = 'affine'
560                    force_pw = True
561            else:
562                self.function_character = {}
563                for t in self.index:
564                    tpl = (self._model(),)+(t,)+(self._pw_pts(t),)
565                    result = characterize_function(self._F, tpl)
566                    if (result == -1):
567                        self.function_character[t] = 'concave'
568                        if (self.PW_TYPE == 'UB'):
569                            simplified_index.append(t)
570                            force_simple = True
571                        else:
572                            force_pw = True
573                            pw_index.append(t)
574                    elif (result == 1):
575                        self.function_character[t] = 'convex'
576                        if (self.PW_TYPE == 'LB'):
577                            simplified_index.append(t)
578                            force_simple = True
579                        else:
580                            force_pw = True
581                            pw_index.append(t)
582                    else:
583                        self.function_character[t] = 'affine'
584                        pw_index.append(t)
585                        force_pw = True
586                       
587            # make sure user doesn't want to force the use of pw constraints
588            if self._force_pw is True:
589                force_simple = False
590                force_pw = True
591                pw_index = self.index
592                simplified_index = []
593           
594            if force_simple is True:
595                self._construct_simple_bounds(TupleSet(self._index_dimen,simplified_index))
596            if force_pw is True:
597                INDEX_SET_PW = TupleSet(self._index_dimen,[])
598                INDEX_SET_SIMPLE = TupleSet(self._index_dimen,[])
599                force_simple = False
600                force_pw = False
601                if self._element_mode:
602                    if len(self._pw_pts()) == 2:
603                        force_simple = True
604                    else:
605                        force_pw = True
606                else:
607                    for t in pw_index:
608                        if len(self._pw_pts(t)) == 2:
609                            INDEX_SET_SIMPLE.add(t)
610                            force_simple = True
611                        else:
612                            INDEX_SET_PW.add(t)
613                            force_pw = True
614                       
615                ####### SIMPLE
616                if force_simple is True:
617                    self._construct_simple_single_bound(INDEX_SET_SIMPLE)
618                   
619                ##### SOS2
620                if force_pw is True:
621                    if self.PW_REP == 'SOS2':
622                        self._construct_SOS2(INDEX_SET_PW)
623                    elif self.PW_REP in ['BIGM_SOS1','BIGM_BIN']:
624                        self._construct_BIGM(INDEX_SET_PW)
625                    elif self.PW_REP == 'DCC':
626                        self._construct_DCC(INDEX_SET_PW)
627                    elif self.PW_REP == 'DLOG':
628                        self._construct_DLOG(INDEX_SET_PW)
629                    elif self.PW_REP == 'CC':
630                        self._construct_CC(INDEX_SET_PW)
631                    elif self.PW_REP == 'LOG':
632                        self._construct_LOG(INDEX_SET_PW)
633                    elif self.PW_REP == 'MC':
634                        self._construct_MC(INDEX_SET_PW)
635                    elif self.PW_REP == 'INC':
636                        self._construct_INC(INDEX_SET_PW)
637
638   
639    def _define_F(self,candidate):
640        # check that f_rule defines a function
641        if candidate == None:
642            raise TypeError, "Piecewise kwd 'f_rule' is required."
643        else:
644            test_point = 0.0
645            if self.PW_PTS != None:
646                if self.PW_PTS.__class__ is dict:
647                    test_point = self.PW_PTS.items()[0]
648                    test_point = flatten_tuple((test_point[0],test_point[1][0]))
649                elif self.PW_PTS.__class__ is tuple or self.PW_PTS.__class__ is list:
650                    test_point = self.PW_PTS[0]
651                else:
652                    raise TypeError, "Piecewise kwd 'pw_pts' must be of type dict, list, or tuple"
653            elif self.GRAD_PTS != None:
654                if self.GRAD_PTS.__class__ is dict:
655                    test_point = self.GRAD_PTS.items()[0]
656                    test_point = flatten_tuple((test_point[0],test_point[1][0]))
657                elif self.GRAD_PTS.__class__ is tuple or self.GRAD_PTS.__class__ is list:
658                    test_point = self.GRAD_PTS[0]
659                else:
660                    raise TypeError, "Piecewise kwd 'tan_pts' must be of type dict, list, or tuple"
661            else:
662                raise AssertionError, "Must specify at least one of 'tan_pts' or 'pw_points' as kwd to Piecewise"
663            try:
664                if test_point.__class__ is tuple:
665                    candidate(self._model(),*(test_point))
666                else:
667                    candidate(self._model(),*(test_point,))
668            except Exception:
669                raise TypeError, "Piecewise kwd 'f_rule' must be callable with the form f(model,x),\n" \
670                                +"where the first argument 'model' is a Pyomo model and the second argument \n" \
671                                +"'x' is a number"
672            else:
673                self._F_imp = candidate
674
675    def _F(self,m,*args):
676        # we cast to float to avoid integer division in case the
677        # user supplies integer type points and return values
678        return float(self._F_imp(m,*flatten_tuple(args)))
679       
680    def _define_d_F(self,candidate):
681        # check that d_f_rule defines a function
682        if (candidate == None) and (self.GRAD_PTS != None):
683            raise TypeError, "Piecewise kwd 'd_f_rule' is required when kwd 'tan_pts' is specified."
684        elif (self.GRAD_PTS == None) and (candidate != None):
685            raise TypeError, "Piecewise kwd 'tan_pts' is required when kwd 'd_f_rule' is specified."
686        else:
687            test_point = 0.0
688            if self.GRAD_PTS != None:
689                if self.GRAD_PTS.__class__ is dict:
690                    test_point = self.GRAD_PTS.items()[0]
691                    test_point = flatten_tuple((test_point[0],test_point[1][0]))
692                elif self.GRAD_PTS.__class__ is tuple or self.GRAD_PTS.__class__ is list:
693                    test_point = self.GRAD_PTS[0]
694                else:
695                    raise TypeError, "Piecewise kwd 'tan_pts' must be of type dict, list, or tuple"
696            try:
697                if test_point.__class__ is tuple:
698                    candidate(self._model(),*(test_point))
699                else:
700                    candidate(self._model(),*(test_point,))
701            except Exception:
702                raise TypeError, "Piecewise kwd 'd_f_rule' must be callable with the form f(model,x),\n" \
703                                +"where the first argument 'model' is a Pyomo model and the second argument \n" \
704                                +"'x' is a number"
705            self._d_F_imp = candidate   
706
707    def _d_F(self,m,*args):
708        # we cast to float to avoid integer division in case the
709        # user supplies integer type points and return values
710        return float(self._d_F_imp(m,*flatten_tuple(args)))
711   
712    def _update_dicts(self,name,_object_):
713        if isinstance(_object_,ConstraintBase):
714            constraint_name = self.name+'_'+name
715            if hasattr(self._model(),constraint_name):
716                raise AssertionError, "Error: "+self.name+" trying to overwrite model attribute: "+constraint_name
717            setattr(self._model(),constraint_name,_object_)
718            self._constraints_dict[constraint_name] = getattr(self._model(), constraint_name)
719            return None
720        elif isinstance(_object_,_VarArray):
721            var_name = self.name+'_'+name
722            if hasattr(self._model(),var_name):
723                raise AssertionError, "Error: "+self.name+" trying to overwrite model attribute: "+var_name
724            setattr(self._model(),var_name,_object_)
725            self._vars_dict.update([(v.name,v) for (n,v) in getattr(self._model(), var_name).iteritems()])
726            return var_name
727        elif isinstance(_object_,_VarElement):
728            var_name = self.name+'_'+name
729            if hasattr(self._model(),var_name):
730                raise AssertionError, "Error: "+self.name+" trying to overwrite model attribute: "+var_name
731            setattr(self._model(),var_name,_object_)
732            self._vars_dict[var_name] = getattr(self._model(), var_name)
733            return var_name
734        elif isinstance(_object_,_BaseSet):
735            set_name = self.name+'_'+name
736            if hasattr(self._model(),set_name):
737                raise AssertionError, "Error: "+self.name+" trying to overwrite model attribute: "+set_name
738            setattr(self._model(),set_name,_object_)
739            self._sets_dict[set_name] = getattr(self._model(), set_name)
740            return set_name
741   
742   
743    def variables(self):
744        """
745        Returns a dictionary of the model variables created by this class.
746        Key, Value pairs are (var.name, var)
747        """
748        return self._vars_dict
749   
750    def constraints(self):
751        """
752        Returns a dictionary of the model constraints created by this class.
753        Key, Value pairs are (constraint.name, constraint)
754        """
755        return self._constraints_dict
756   
757    def sets(self):
758        """
759        Returns a dictionary of the model sets created by this class.
760        Key, Value pairs are (set.name, set)
761        """
762        return self._sets_dict
763   
764    def components(self):
765        """
766        Returns a dictionary of the model components created by this class.
767        Key, Value pairs are (component.name, component)
768        """
769        tmp_dict = {}
770        tmp_dict.update(self._vars_dict)
771        tmp_dict.update(self._sets_dict)
772        tmp_dict.update(self._constraints_dict)
773        return tmp_dict
774
775    def _construct_linear_estimators(self, INDEX_SET_TANGENT):
776       
777        def linear_estimators_LINEAR_indices_init(model):
778            if self._element_mode:
779                return (i for i in xrange(len(self._grad_pts())))
780            else:
781                return (t+(i,) for t in INDEX_SET_TANGENT.tuplize() \
782                              for i in xrange(len(self._grad_pts(t))))
783        tangent_indices = self._update_dicts('tangent_indices',Set(dimen=self._index_dimen+1, ordered=True, rule=linear_estimators_LINEAR_indices_init))
784
785
786        def linear_estimators_constraint_rule(model,*args):
787            t = args[:-1]
788            i = args[-1]
789            LHS = self._getvar(model,self.YVAR.name,t)
790            F_AT_XO = self._F(model,t,self._grad_pts(t)[i])
791            dF_AT_XO = self._d_F(model,t,self._grad_pts(t)[i])
792            X_MINUS_XO = (self._getvar(model,self.XVAR.name,t)-self._grad_pts(t)[i])
793            if self.GRAD_TYPE == 'LB':
794                return LHS >= F_AT_XO + dF_AT_XO*X_MINUS_XO
795            elif self.GRAD_TYPE == 'UB':
796                return LHS <= F_AT_XO + dF_AT_XO*X_MINUS_XO
797
798        self._update_dicts('tangent_constraint' , Constraint(getattr(self._model(),tangent_indices),rule=linear_estimators_constraint_rule))
799   
800       
801    def _construct_simple_single_bound(self, INDEX_SET_SIMPLE):
802       
803        def SIMPLE_LINEAR_constraint1_rule(model,*t):
804            LHS = self._getvar(model,self.YVAR.name,t)
805            F_AT_XO = self._F(model,t,min(self._pw_pts(t)))
806            dF_AT_XO = (float(self._F(model,t,max(self._pw_pts(t)))-self._F(model,t,min(self._pw_pts(t))))/float(max(self._pw_pts(t))-min(self._pw_pts(t))))
807            X_MINUS_XO = (self._getvar(model,self.XVAR.name,t)-min(self._pw_pts(t)))
808            if self.PW_TYPE == 'UB':
809                return LHS <= F_AT_XO + dF_AT_XO*X_MINUS_XO     
810            elif self.PW_TYPE == 'LB':
811                return LHS >= F_AT_XO + dF_AT_XO*X_MINUS_XO     
812            elif self.PW_TYPE == 'EQ':
813                return LHS == F_AT_XO + dF_AT_XO*X_MINUS_XO
814
815        self._update_dicts('simplified_single_line_constraint',Constraint(INDEX_SET_SIMPLE,rule=SIMPLE_LINEAR_constraint1_rule))
816       
817    def _construct_simple_bounds(self, INDEX_SET_SIMPLE_BOUNDS):
818       
819        def simple_bounds_indices_init(model):
820            if self._element_mode:
821                return (i for i in xrange(len(self._pw_pts())-1))
822            else:
823                return (t+(i,) for t in INDEX_SET_SIMPLE_BOUNDS.tuplize() \
824                              for i in xrange(len(self._pw_pts(t))-1))
825        simple_bounds_indices = self._update_dicts('simplified_bounds_indices',Set(dimen=self._index_dimen+1,ordered=True, rule=simple_bounds_indices_init))
826       
827       
828        def SIMPLE_LINEAR_constraint1_rule(model,*args):
829            t = args[:-1]
830            i = args[-1]
831            LHS = self._getvar(model,self.YVAR.name,t)
832            F_AT_XO = self._F(model,t,self._pw_pts(t)[i])
833            dF_AT_XO = (float(self._F(model,t,self._pw_pts(t)[i+1])-self._F(model,t,self._pw_pts(t)[i]))/float(self._pw_pts(t)[i+1]-self._pw_pts(t)[i]))
834            X_MINUS_XO = (self._getvar(model,self.XVAR.name,t)-self._pw_pts(t)[i])
835            if self.PW_TYPE == 'UB':
836                return LHS <= F_AT_XO + dF_AT_XO*X_MINUS_XO
837            elif self.PW_TYPE == 'LB':
838                return LHS >= F_AT_XO + dF_AT_XO*X_MINUS_XO
839
840        self._update_dicts('simplified_bounds_constraint',Constraint(getattr(self._model(),simple_bounds_indices),rule=SIMPLE_LINEAR_constraint1_rule))
841   
842    def _construct_SOS2(self, INDEX_SET_SOS2):
843       
844        def SOS_indices_init(*t):
845            if self._element_mode:
846                return (i for i in xrange(len(self._pw_pts(t))))
847            else:
848                if type(t[0]) is tuple:
849                    return (t[0]+(i,) for i in xrange(len(self._pw_pts(t))))
850                else:
851                    return ((t[0],i) for i in xrange(len(self._pw_pts(t))))
852        def y_indices_init(model):
853            if self._element_mode:
854                return (i for i in xrange(len(self._pw_pts())))
855            else:
856                return (t+(i,) for t in INDEX_SET_SOS2.tuplize() for i in xrange(len(self._pw_pts(t))))
857
858        SOS_indices = ''
859        if self._element_mode:
860            SOS_indices += self._update_dicts('SOS2_indices',Set(INDEX_SET_SOS2,dimen=self._index_dimen+1, ordered=True, \
861                                         initialize=SOS_indices_init()))
862        else:
863            SOS_indices += self._update_dicts('SOS2_indices',Set(INDEX_SET_SOS2,dimen=self._index_dimen+1, ordered=True, \
864                                         initialize=dict([(t,SOS_indices_init(t)) for t in INDEX_SET_SOS2 ])))
865        y_indices = self._update_dicts('y_SOS2_indices',Set(ordered=True, dimen=self._index_dimen+1,initialize=y_indices_init))
866       
867        #Add SOS2 weighting variable to self._model()
868        y = self._update_dicts('y_SOS2',Var(getattr(self._model(),y_indices),within=NonNegativeReals))
869       
870        #Add SOS2 constraints to self._model()
871        def SOS2_constraint1_rule(model,*t):
872            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,y,t,i)*self._pw_pts(t)[i] for i in xrange(len(self._pw_pts(t))))
873        def SOS2_constraint2_rule(model,*t):
874            LHS = self._getvar(model,self.YVAR.name,t)
875            RHS = sum(self._getvar(model,y,t,i)*self._F(model,t,self._pw_pts(t)[i]) for i in xrange(len(self._pw_pts(t))))
876            if self.PW_TYPE == 'UB':
877                return LHS <= RHS
878            elif self.PW_TYPE == 'LB':
879                return LHS >= RHS
880            elif self.PW_TYPE == 'EQ':
881                return LHS == RHS
882        def SOS2_constraint3_rule(model,*t):
883            return sum(self._getvar(model,y,t,j) for j in xrange(len(self._pw_pts(t)))) == 1
884       
885        if self._element_mode is True:
886            self._update_dicts('SOS2_constraint1',Constraint(rule=SOS2_constraint1_rule))
887            self._update_dicts('SOS2_constraint2',Constraint(rule=SOS2_constraint2_rule))
888            self._update_dicts('SOS2_constraint3',Constraint(rule=SOS2_constraint3_rule))
889            self._update_dicts('SOS2_constraint4',SOSConstraint(var=getattr(self._model(),y), sos=2))
890        else:
891            self._update_dicts('SOS2_constraint1',Constraint(INDEX_SET_SOS2,rule=SOS2_constraint1_rule))
892            self._update_dicts('SOS2_constraint2',Constraint(INDEX_SET_SOS2,rule=SOS2_constraint2_rule))
893            self._update_dicts('SOS2_constraint3',Constraint(INDEX_SET_SOS2,rule=SOS2_constraint3_rule))
894            self._update_dicts('SOS2_constraint4',SOSConstraint(INDEX_SET_SOS2, var=getattr(self._model(),y), set=getattr(self._model(),SOS_indices), sos=2))
895
896    def _construct_BIGM(self,BIGM_SET):
897       
898        #Add dictionary of optimal big M values as attribute to self._model()
899        #same index as SOS variable
900        OPT_M = {}
901        OPT_M['UB'] = {}
902        OPT_M['LB'] = {}
903
904        if self.PW_TYPE in ['UB','EQ']:
905            OPT_M['UB'] = self._find_M(self.PW_PTS, 'UB')
906        if self.PW_TYPE in ['LB','EQ']:
907            OPT_M['LB'] = self._find_M(self.PW_PTS, 'LB')
908       
909        all_keys = set(OPT_M['UB'].keys()).union(OPT_M['LB'].keys())
910        full_indices = []
911        if self._element_mode:
912            full_indices.extend([i for i in xrange(1,len(self._pw_pts()))])
913        else:
914            full_indices.extend([t+(i,) for t in BIGM_SET.tuplize() for i in xrange(1,len(self._pw_pts(t)))])
915        y_indices = ''
916        y = ''
917        if len(all_keys) > 0:
918            y_indices += self._update_dicts('y_'+self.PW_REP+'_indices',Set(ordered=True, dimen=self._index_dimen+1,initialize=all_keys))
919           
920            #Add indicator variable to self._model()
921            def y_domain():
922                if self.PW_REP == 'BIGM_BIN':
923                    return Binary
924                elif self.PW_REP == 'BIGM_SOS1':
925                    return NonNegativeReals
926            y += self._update_dicts('y_'+self.PW_REP,Var(getattr(self._model(),y_indices),within=y_domain()))
927       
928        def BIGM_LINEAR_constraint1_rule(model,*args):
929            t = args[:-1]
930            i = args[-1]
931            if self._element_mode:
932                targs = i
933            else:
934                targs = flatten_tuple(args)
935            if (self.PW_TYPE == 'UB') or (self.PW_TYPE == 'EQ'):
936                rhs = 1.0
937                if targs not in OPT_M['UB'].keys():
938                    rhs *= 0.0
939                else:
940                    rhs *= OPT_M['UB'][targs]*(1-self._getvar(model,y,t,i))
941                return self._getvar(model,self.YVAR.name,t) - (self._F(model,t,self._pw_pts(t)[i-1]) + ((self._F(model,t,self._pw_pts(t)[i])-self._F(model,t,self._pw_pts(t)[i-1]))/(self._pw_pts(t)[i]-self._pw_pts(t)[i-1]))*(self._getvar(model,self.XVAR.name,t)-self._pw_pts(t)[i-1])) <= rhs
942            elif self.PW_TYPE == 'LB':
943                rhs = 1.0
944                if targs not in OPT_M['LB'].keys():
945                    rhs *= 0.0
946                else:
947                    rhs *= OPT_M['LB'][targs]*(1-self._getvar(model,y,t,i))
948                return self._getvar(model,self.YVAR.name,t) - (self._F(model,t,self._pw_pts(t)[i-1]) + ((self._F(model,t,self._pw_pts(t)[i])-self._F(model,t,self._pw_pts(t)[i-1]))/(self._pw_pts(t)[i]-self._pw_pts(t)[i-1]))*(self._getvar(model,self.XVAR.name,t)-self._pw_pts(t)[i-1])) >= rhs
949        def BIGM_LINEAR_constraint2_rule(model,*t):
950            if self._element_mode:
951                expr = [self._getvar(model,y,t,i) for i in xrange(1,len(self._pw_pts())) if i in all_keys]
952            else:
953                expr = [self._getvar(model,y,t,i) for i in xrange(1,len(self._pw_pts(t))) if flatten_tuple((t,i)) in all_keys]
954            if len(expr) > 0:
955                return sum(expr) == 1
956            else:
957                return Constraint.Skip
958        def BIGM_LINEAR_constraintAFF_rule(model,*args):
959            t = args[:-1]
960            i = args[-1]
961            targs = 0
962            if self._element_mode:
963                targs = i
964            else:
965                targs = flatten_tuple(args)
966            rhs = 1.0
967            if targs not in OPT_M['LB'].keys():
968                rhs *= 0.0
969            else:
970                rhs *= OPT_M['LB'][targs]*(1-self._getvar(model,y,t,i))
971            return self._getvar(model,self.YVAR.name,t) - (self._F(model,t,self._pw_pts(t)[i-1]) + ((self._F(model,t,self._pw_pts(t)[i])-self._F(model,t,self._pw_pts(t)[i-1]))/(self._pw_pts(t)[i]-self._pw_pts(t)[i-1]))*(self._getvar(model,self.XVAR.name,t)-self._pw_pts(t)[i-1])) >= rhs
972
973        self._update_dicts(self.PW_REP+'_constraint1',Constraint(full_indices,rule=BIGM_LINEAR_constraint1_rule))
974        if len(all_keys) > 0:
975            if self._element_mode:
976                self._update_dicts(self.PW_REP+'_constraint2',Constraint(rule=BIGM_LINEAR_constraint2_rule))
977            else:
978                self._update_dicts(self.PW_REP+'_constraint2',Constraint(BIGM_SET,rule=BIGM_LINEAR_constraint2_rule))
979        if self.PW_TYPE == 'EQ':
980                self._update_dicts(self.PW_REP+'_constraint3',Constraint(full_indices,rule=BIGM_LINEAR_constraintAFF_rule))
981       
982        if len(all_keys) > 0:
983            if (self.PW_REP == 'BIGM_SOS1'):
984                def BIGM_SOS1_indices_init(*t):
985                    if self._element_mode is True:
986                        return [i for i in xrange(len(self._pw_pts())) if i in all_keys]
987                    else:
988                        if type(t[0]) is tuple:
989                            return [t[0]+(i,) for i in xrange(1,len(self._pw_pts(t))) if t[0]+(i,) in all_keys]
990                        else:
991                            return [(t[0],i) for i in xrange(1,len(self._pw_pts(t))) if (t[0],i) in all_keys]
992           
993                if self._element_mode:
994                    self._update_dicts('BIGM_SOS1_constraintSOS',SOSConstraint(var=getattr(self._model(),y), sos=1))
995                else:
996                    SOS_indices = self._update_dicts('BIGM_SOS1_indices',Set(list(BIGM_SET),dimen=self._index_dimen+1, ordered=True, \
997                                                                            initialize=dict([(t,BIGM_SOS1_indices_init(t)) for t in BIGM_SET])))
998                    self._update_dicts('BIGM_SOS1_constraintSOS',SOSConstraint(BIGM_SET, var=getattr(self._model(),y), set=getattr(self._model(),SOS_indices), sos=1))
999
1000    def _construct_DCC(self,DCC_SET):
1001       
1002        def DCC_POLYTOPES_init(model):
1003            if self._element_mode:
1004                return (i for i in xrange(1,len(self._pw_pts())))
1005
1006            else:
1007                return (t+(i,) for t in DCC_SET.tuplize() \
1008                              for i in xrange(1,len(self._pw_pts(t))))
1009        def DCC_poly_init(*t):
1010            return (i for i in xrange(1,len(self._pw_pts(t))))
1011        def DCC_VERTICES_init(model):
1012            if self._element_mode:
1013                return (i for i in xrange(1,len(self._pw_pts())+1))
1014            else:
1015                return (t+(i,) for t in DCC_SET.tuplize() \
1016                              for i in xrange(1,len(self._pw_pts(t))+1))
1017        def DCC_vert_init(args):
1018            try:
1019                t = args[:-1]
1020                p = args[-1]
1021            except:
1022                p = args
1023            return (i for i in xrange(p,p+2))
1024        def DCC_lamba_set_init(model):
1025            if self._element_mode is True:
1026                return ((p,v) for p in xrange(1,len(self._pw_pts())) \
1027                                for v in xrange(1,len(self._pw_pts())+1))
1028            else:
1029                # BOTTLENECK: speed/memory, this set can be extremely large and time consuming to
1030                #             build when number of indices and breakpoints is large.
1031                return (t+(p,v) for t in DCC_SET.tuplize() \
1032                                for p in xrange(1,len(self._pw_pts(t))) \
1033                                for v in xrange(1,len(self._pw_pts(t))+1))
1034       
1035        DCC_POLYTOPES = self._update_dicts('DCC_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=DCC_POLYTOPES_init))
1036        if self._element_mode:
1037            DCC_poly = self._update_dicts('DCC_poly',Set(ordered=True, initialize=DCC_poly_init()))
1038        else:
1039            DCC_poly = self._update_dicts('DCC_poly',Set(DCC_SET,ordered=True, \
1040                                                         initialize=dict([(t,DCC_poly_init(t)) for t in DCC_SET])))
1041        DCC_VERTICES = self._update_dicts('DCC_VERTICES',Set(ordered=True,dimen=self._index_dimen+1,initialize=DCC_VERTICES_init))
1042        DCC_vert = self._update_dicts('DCC_vert',Set(getattr(self._model(),DCC_POLYTOPES),ordered=True, \
1043                                                     initialize=dict([(args,DCC_vert_init(args)) for args in getattr(self._model(),DCC_POLYTOPES)])))   
1044        DCC_LAMBDA_SET = self._update_dicts('DCC_LAMBDA_SET',Set(ordered=True,dimen=self._index_dimen+2,initialize=DCC_lamba_set_init))
1045       
1046        lmda = self._update_dicts('DCC_LAMBDA',Var(getattr(self._model(),DCC_LAMBDA_SET),within=PositiveReals))
1047       
1048       
1049        y = self._update_dicts('DCC_BIN_y',Var(getattr(self._model(),DCC_POLYTOPES),within=Binary))
1050       
1051        def DCC_LINEAR_constraint1_rule(model,*t):
1052            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,lmda,t,p,v)*self._pw_pts(t)[v-1] for p in self._getset(model,DCC_poly,t) for v in self._getset(model,DCC_vert,t,p))
1053        def DCC_LINEAR_constraint2_rule(model,*t):
1054            LHS = self._getvar(model,self.YVAR.name,t)
1055            RHS = sum(self._getvar(model,lmda,t,p,v)*self._F(model,t,self._pw_pts(t)[v-1]) for p in self._getset(model,DCC_poly,t) for v in self._getset(model,DCC_vert,t,p))
1056            if self.PW_TYPE == 'UB':
1057                return LHS <= RHS
1058            elif self.PW_TYPE == 'LB':
1059                return LHS >= RHS
1060            elif self.PW_TYPE == 'EQ':
1061                return LHS == RHS
1062        def DCC_LINEAR_constraint3_rule(model,*args):
1063            t = args[:-1]
1064            p = args[-1]
1065            return self._getvar(model,y,t,p) == sum(self._getvar(model,lmda,t,p,v) for v in self._getset(model,DCC_vert,t,p))
1066        def DCC_LINEAR_constraint4_rule(model,*t):
1067            return sum(self._getvar(model,y,t,p) for p in self._getset(model,DCC_poly,t)) == 1
1068
1069        if self._element_mode:
1070            self._update_dicts('DCC_constraint1',Constraint(rule=DCC_LINEAR_constraint1_rule))
1071            self._update_dicts('DCC_constraint2',Constraint(rule=DCC_LINEAR_constraint2_rule))
1072            self._update_dicts('DCC_constraint3',Constraint(getattr(self._model(),DCC_POLYTOPES),rule=DCC_LINEAR_constraint3_rule))
1073            self._update_dicts('DCC_constraint4',Constraint(rule=DCC_LINEAR_constraint4_rule))
1074        else:
1075            self._update_dicts('DCC_constraint1',Constraint(DCC_SET,rule=DCC_LINEAR_constraint1_rule))
1076            self._update_dicts('DCC_constraint2',Constraint(DCC_SET,rule=DCC_LINEAR_constraint2_rule))
1077            self._update_dicts('DCC_constraint3',Constraint(getattr(self._model(),DCC_POLYTOPES),rule=DCC_LINEAR_constraint3_rule))
1078            self._update_dicts('DCC_constraint4',Constraint(DCC_SET,rule=DCC_LINEAR_constraint4_rule))
1079
1080
1081    def _construct_DLOG(self,DLOG_SET):
1082       
1083        L = {}
1084        L_i = 0
1085        B_ZERO = {}
1086        B_ZERO_l = []
1087        B_ONE = {}
1088        B_ONE_l = []
1089        if self._element_mode:
1090            L_i += int(math.log(len(self._pw_pts())-1,2))
1091            B_ZERO_l,B_ONE_l = _DLog_Branching_Scheme(L_i)
1092        else:
1093            for t in DLOG_SET:
1094                L[t] = int(math.log(len(self._pw_pts(t))-1,2))
1095                L[(t,)] = L[t]
1096            for t in sorted(L.keys()):
1097                B_ZERO[t],B_ONE[t] = _DLog_Branching_Scheme(L[t])
1098                B_ZERO[(t,)] = B_ZERO[t]
1099                B_ONE[(t,)] = B_ONE[t]
1100           
1101        def DLOG_POLYTOPES_init(model):
1102            if self._element_mode:
1103                return (i for i in xrange(1,len(self._pw_pts())))
1104
1105            else:
1106                return (t+(i,) for t in DLOG_SET.tuplize() \
1107                              for i in xrange(1,len(self._pw_pts(t))))
1108        def DLOG_LENGTH_POLY_init(model):
1109            if self._element_mode:
1110                return (i for i in xrange(1,L_i+1))
1111            else:
1112                return (t+(i,) for t in DLOG_SET.tuplize() \
1113                              for i in xrange(1,L[t]+1))
1114        def DLOG_poly_init(*t):
1115            return (i for i in xrange(1,len(self._pw_pts(t))))
1116        def DLOG_poly_one_init(args):
1117            if self._element_mode:
1118                l = args
1119                return B_ZERO_l[l]
1120            else:
1121                t = args[:-1]
1122                l = args[-1]
1123                return B_ZERO[t][l]
1124        def DLOG_poly_zero_init(args):
1125            if self._element_mode:
1126                l = args
1127                return B_ONE_l[l]
1128            else:
1129                t = args[:-1]
1130                l = args[-1]
1131                return B_ONE[t][l]
1132        def DLOG_VERTICES_init(model):
1133            if self._element_mode:
1134                return (i for i in xrange(1,len(self._pw_pts())+1))
1135            else:
1136                return (t+(i,) for t in DLOG_SET.tuplize() \
1137                              for i in xrange(1,len(self._pw_pts(t))+1))
1138        def DLOG_vert_init(args):
1139            try:
1140                t = args[:-1]
1141                p = args[-1]
1142            except:
1143                p = args
1144            return (i for i in xrange(p,p+2))
1145        def DLOG_lamba_set_init(model):
1146            if self._element_mode:
1147                return ((p,v) for p in xrange(1,len(self._pw_pts())) \
1148                              for v in xrange(1,len(self._pw_pts())+1))
1149            else:
1150                # BOTTLENECK: speed/memory, this set can be extremely large and time consuming to
1151                #             build when number of indices and breakpoints is large.
1152                return (t+(p,v) for t in DLOG_SET.tuplize() \
1153                                for p in xrange(1,len(self._pw_pts(t))) \
1154                                for v in xrange(1,len(self._pw_pts(t))+1))
1155       
1156        DLOG_POLYTOPES = self._update_dicts('DLOG_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=DLOG_POLYTOPES_init))
1157        DLOG_LENGTH_POLY = self._update_dicts('DLOG_LENGTH_POLY',Set(ordered=True,dimen=self._index_dimen+1,initialize=DLOG_LENGTH_POLY_init))
1158        if self._element_mode:
1159            DLOG_poly = self._update_dicts('DLOG_poly',Set(ordered=True, \
1160                                                           initialize=DLOG_poly_init()))
1161        else:
1162            DLOG_poly = self._update_dicts('DLOG_poly',Set(DLOG_SET,ordered=True, \
1163                                                           initialize=dict([(t,DLOG_poly_init(t)) for t in DLOG_SET])))
1164        DLOG_poly_one = self._update_dicts('DLOG_poly_one',Set(getattr(self._model(),DLOG_LENGTH_POLY),ordered=True, \
1165                                                               initialize=dict([(args,DLOG_poly_one_init(args)) for args in getattr(self._model(),DLOG_LENGTH_POLY)])))
1166        DLOG_poly_zero = self._update_dicts('DLOG_poly_zero',Set(getattr(self._model(),DLOG_LENGTH_POLY),ordered=True, \
1167                                                                 initialize=dict([(args,DLOG_poly_zero_init(args)) for args in getattr(self._model(),DLOG_LENGTH_POLY)])))
1168        DLOG_VERTICES = self._update_dicts('DLOG_VERTICES',Set(ordered=True,dimen=self._index_dimen+1,initialize=DLOG_VERTICES_init))
1169        DLOG_vert = self._update_dicts('DLOG_vert',Set(getattr(self._model(),DLOG_POLYTOPES),ordered=True, \
1170                                                       initialize=dict([(args,DLOG_vert_init(args)) for args in getattr(self._model(),DLOG_POLYTOPES)])))       
1171        DLOG_LAMBDA_SET = self._update_dicts('DLOG_LAMBDA_SET',Set(ordered=True,dimen=self._index_dimen+2,initialize=DLOG_lamba_set_init))
1172       
1173        lmda = self._update_dicts('DLOG_LAMBDA',Var(getattr(self._model(),DLOG_LAMBDA_SET),within=PositiveReals))
1174        y = self._update_dicts('DLOG_BIN_y',Var(getattr(self._model(),DLOG_LENGTH_POLY),within=Binary))
1175       
1176        def DLOG_LINEAR_constraint1_rule(model,*t):
1177            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,lmda,t,p,v)*self._pw_pts(t)[v-1] for p in self._getset(model,DLOG_poly,t) for v in self._getset(model,DLOG_vert,t,p))
1178        def DLOG_LINEAR_constraint2_rule(model,*t):
1179            LHS = self._getvar(model,self.YVAR.name,t)
1180            RHS = sum(self._getvar(model,lmda,t,p,v)*self._F(model,t,self._pw_pts(t)[v-1]) for p in self._getset(model,DLOG_poly,t) for v in self._getset(model,DLOG_vert,t,p))
1181            if self.PW_TYPE == 'UB':
1182                return LHS <= RHS
1183            elif self.PW_TYPE == 'LB':
1184                return LHS >= RHS
1185            elif self.PW_TYPE == 'EQ':
1186                return LHS == RHS
1187        def DLOG_LINEAR_constraint3_rule(model,*t):
1188            return 1 == sum(self._getvar(model,lmda,t,p,v) for p in self._getset(model,DLOG_poly,t) for v in self._getset(model,DLOG_vert,t,p))
1189        def DLOG_LINEAR_constraint4_rule(model,*args):
1190            t = args[:-1]
1191            l = args[-1]
1192            return sum(self._getvar(model,lmda,t,p,v) for p in self._getvar(model,DLOG_poly_one,t,l) for v in self._getset(model,DLOG_vert,t,p)) <= self._getvar(model,y,t,l)
1193        def DLOG_LINEAR_constraint5_rule(model,*args):
1194            t = args[:-1]
1195            l = args[-1]
1196            return sum(self._getvar(model,lmda,t,p,v) for p in self._getvar(model,DLOG_poly_zero,t,l) for v in self._getset(model,DLOG_vert,t,p)) <= (1-self._getvar(model,y,t,l))
1197
1198        if self._element_mode:
1199            self._update_dicts('DLOG_constraint1',Constraint(rule=DLOG_LINEAR_constraint1_rule))
1200            self._update_dicts('DLOG_constraint2',Constraint(rule=DLOG_LINEAR_constraint2_rule))
1201            self._update_dicts('DLOG_constraint3',Constraint(rule=DLOG_LINEAR_constraint3_rule))
1202            self._update_dicts('DLOG_constraint4',Constraint(getattr(self._model(),DLOG_LENGTH_POLY),rule=DLOG_LINEAR_constraint4_rule))
1203            self._update_dicts('DLOG_constraint5',Constraint(getattr(self._model(),DLOG_LENGTH_POLY),rule=DLOG_LINEAR_constraint5_rule))
1204        else:
1205            self._update_dicts('DLOG_constraint1',Constraint(DLOG_SET,rule=DLOG_LINEAR_constraint1_rule))
1206            self._update_dicts('DLOG_constraint2',Constraint(DLOG_SET,rule=DLOG_LINEAR_constraint2_rule))
1207            self._update_dicts('DLOG_constraint3',Constraint(DLOG_SET,rule=DLOG_LINEAR_constraint3_rule))
1208            self._update_dicts('DLOG_constraint4',Constraint(getattr(self._model(),DLOG_LENGTH_POLY),rule=DLOG_LINEAR_constraint4_rule))
1209            self._update_dicts('DLOG_constraint5',Constraint(getattr(self._model(),DLOG_LENGTH_POLY),rule=DLOG_LINEAR_constraint5_rule))
1210
1211    def _construct_CC(self,CC_SET):
1212       
1213        def CC_POLYTOPES_init(model):
1214            if self._element_mode:
1215                return (i for i in xrange(1,len(self._pw_pts())))
1216            else:
1217                return (t+(i,) for t in CC_SET.tuplize() \
1218                              for i in xrange(1,len(self._pw_pts(t))))
1219        def CC_poly_init(args):
1220            if self._element_mode:
1221                v = args
1222                if v == 1:
1223                    return [v]
1224                elif v == len(self._pw_pts()):
1225                    return [v-1]
1226                else:
1227                    return [v-1,v]
1228            else:
1229                t = args[:-1]
1230                v = args[-1]
1231                if v == 1:
1232                    return [v]
1233                elif v == len(self._pw_pts(t)):
1234                    return [v-1]
1235                else:
1236                    return [v-1,v]
1237        def CC_VERTICES_init(model):
1238            if self._element_mode:
1239                return (i for i in xrange(1,len(self._pw_pts())+1))
1240            else:
1241                return (t+(i,) for t in CC_SET.tuplize() \
1242                              for i in xrange(1,len(self._pw_pts(t))+1))
1243        def CC_vert_init(*t):
1244            return (i for i in xrange(1,len(self._pw_pts(t))+1))
1245        def CC_lamba_set_init(model):
1246            if self._element_mode:
1247                return (v for v in xrange(1,len(self._pw_pts())+1))
1248            else:
1249                return (t+(v,) for t in CC_SET.tuplize() \
1250                              for v in xrange(1,len(self._pw_pts(t))+1))
1251       
1252        CC_VERTICES = self._update_dicts('CC_VERTICES',Set(ordered=True,dimen=self._index_dimen+1,initialize=CC_VERTICES_init))
1253        if self._element_mode:
1254            CC_vert = self._update_dicts('CC_vert',Set(ordered=True, initialize=CC_vert_init()))
1255        else:
1256            CC_vert = self._update_dicts('CC_vert',Set(CC_SET,ordered=True,\
1257                                                       initialize=dict([(t,CC_vert_init(t)) for t in CC_SET])))
1258        CC_POLYTOPES = self._update_dicts('CC_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=CC_POLYTOPES_init))
1259        CC_poly = self._update_dicts('CC_poly',Set(getattr(self._model(),CC_VERTICES),ordered=True, \
1260                                                   initialize=dict([(args,CC_poly_init(args)) for args in getattr(self._model(),CC_VERTICES)])))
1261        CC_LAMBDA_SET = self._update_dicts('CC_LAMBDA_SET',Set(ordered=True,dimen=self._index_dimen+1,initialize=CC_lamba_set_init))
1262       
1263        lmda = self._update_dicts('CC_LAMBDA',Var(getattr(self._model(),CC_LAMBDA_SET),within=NonNegativeReals))
1264        y = self._update_dicts('CC_BIN_y',Var(getattr(self._model(),CC_POLYTOPES),within=Binary))
1265
1266       
1267        def CC_LINEAR_constraint1_rule(model,*t):
1268            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,lmda,t,v)*self._pw_pts(t)[v-1] for v in self._getset(model,CC_vert,t))
1269        def CC_LINEAR_constraint2_rule(model,*t):
1270            LHS = self._getvar(model,self.YVAR.name,t)
1271            RHS = sum(self._getvar(model,lmda,t,v)*self._F(model,t,self._pw_pts(t)[v-1]) for v in self._getset(model,CC_vert,t))
1272            if self.PW_TYPE == 'UB':
1273                return LHS <= RHS
1274            elif self.PW_TYPE == 'LB':
1275                return LHS >= RHS
1276            elif self.PW_TYPE == 'EQ':
1277                return LHS == RHS
1278        def CC_LINEAR_constraint3_rule(model,*t):
1279            return 1 == sum(self._getvar(model,lmda,t,v) for v in self._getset(model,CC_vert,t))
1280        def CC_LINEAR_constraint4_rule(model,*args):
1281            t = args[:-1]
1282            v = args[-1]
1283            return self._getvar(model,lmda,t,v) <= sum(self._getvar(model,y,t,p) for p in self._getset(model,CC_poly,t,v))
1284        def CC_LINEAR_constraint5_rule(model,*t):
1285            return sum(self._getvar(model,y,t,p) for p in xrange(1,len(self._pw_pts(t)))) == 1
1286
1287        if self._element_mode:
1288            self._update_dicts('CC_constraint1',Constraint(rule=CC_LINEAR_constraint1_rule))
1289            self._update_dicts('CC_constraint2',Constraint(rule=CC_LINEAR_constraint2_rule))
1290            self._update_dicts('CC_constraint3',Constraint(rule=CC_LINEAR_constraint3_rule))
1291            self._update_dicts('CC_constraint4',Constraint(getattr(self._model(),CC_VERTICES),rule=CC_LINEAR_constraint4_rule))
1292            self._update_dicts('CC_constraint5',Constraint(rule=CC_LINEAR_constraint5_rule))
1293        else:
1294            self._update_dicts('CC_constraint1',Constraint(CC_SET,rule=CC_LINEAR_constraint1_rule))
1295            self._update_dicts('CC_constraint2',Constraint(CC_SET,rule=CC_LINEAR_constraint2_rule))
1296            self._update_dicts('CC_constraint3',Constraint(CC_SET,rule=CC_LINEAR_constraint3_rule))
1297            self._update_dicts('CC_constraint4',Constraint(getattr(self._model(),CC_VERTICES),rule=CC_LINEAR_constraint4_rule))
1298            self._update_dicts('CC_constraint5',Constraint(CC_SET,rule=CC_LINEAR_constraint5_rule))
1299
1300    def _construct_LOG(self,LOG_SET):
1301       
1302        L = {}
1303        Li = 0
1304        S = {}
1305        Si = 0
1306        LEFT = {}
1307        LEFT_l = []
1308        RIGHT = {}
1309        RIGHT_l = []
1310        if self._element_mode:
1311            Li += int(math.log(len(self._pw_pts())-1,2))
1312            Si,LEFT_l,RIGHT_l = _Log_Branching_Scheme(Li)
1313        else:
1314            for t in LOG_SET:
1315                L[t] = int(math.log(len(self._pw_pts(t))-1,2))
1316                S[t],LEFT[t],RIGHT[t] = _Log_Branching_Scheme(L[t])
1317                S[(t,)] = S[t]
1318                LEFT[(t,)] = LEFT[t]
1319                RIGHT[(t,)] = RIGHT[t]
1320       
1321       
1322        def LOG_POLYTOPES_init(model):
1323            if self._element_mode:
1324                return (i for i in xrange(1,len(self._pw_pts())))
1325            else:
1326                return (t+(i,) for t in LOG_SET.tuplize() \
1327                              for i in xrange(1,len(self._pw_pts(t))))
1328        def LOG_poly_init(args):
1329            if self._element_mode:
1330                v = args
1331                if v == 1:
1332                    return [v]
1333                elif v == len(self._pw_pts()):
1334                    return [v-1]
1335                else:
1336                    return [v-1,v]
1337            else:
1338                t = args[:-1]
1339                v = args[-1]
1340                if v == 1:
1341                    return [v]
1342                elif v == len(self._pw_pts(t)):
1343                    return [v-1]
1344                else:
1345                    return [v-1,v]
1346
1347        def LOG_VERTICES_init(model):
1348            if self._element_mode:
1349                return (i for i in xrange(1,len(self._pw_pts())+1))
1350            else:
1351                return (t+(i,) for t in LOG_SET.tuplize() \
1352                              for i in xrange(1,len(self._pw_pts(t))+1))
1353        def LOG_vert_init(*t):
1354            return (i for i in xrange(1,len(self._pw_pts(t))+1))
1355        def LOG_lamba_set_init(model):
1356            if self._element_mode:
1357                return (v for v in xrange(1,len(self._pw_pts())+1))
1358            else:
1359                return (t+(v,) for t in LOG_SET.tuplize() \
1360                              for v in xrange(1,len(self._pw_pts(t))+1))
1361        def LOG_BRANCHING_SCHEME_init(model):
1362            if self._element_mode:
1363                return (s for s in Si)
1364            else:
1365                return (t+(s,) for t in LOG_SET.tuplize() \
1366                              for s in S[t])
1367        def LOG_BRANCHING_LEFT_init(args):
1368            if self._element_mode:
1369                s = args
1370                return LEFT_l[s]
1371            else:
1372                t = args[:-1]
1373                s = args[-1]
1374                return LEFT[t][s]
1375        def LOG_BRANCHING_RIGHT_init(args):
1376            if self._element_mode:
1377                s = args
1378                return RIGHT_l[s]
1379            else:
1380                t = args[:-1]
1381                s = args[-1]
1382                return RIGHT[t][s]
1383       
1384        LOG_VERTICES = self._update_dicts('LOG_VERTICES',Set(ordered=True,dimen=self._index_dimen+1,initialize=LOG_VERTICES_init))
1385        if self._element_mode:
1386            LOG_vert = self._update_dicts('LOG_vert',Set(ordered=True, initialize=LOG_vert_init()))
1387        else:
1388            LOG_vert = self._update_dicts('LOG_vert',Set(LOG_SET,ordered=True, \
1389                                                         initialize=dict([(t,LOG_vert_init(t)) for t in LOG_SET])))
1390        LOG_POLYTOPES = self._update_dicts('LOG_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=LOG_POLYTOPES_init))
1391        LOG_poly = self._update_dicts('LOG_poly',Set(getattr(self._model(),LOG_VERTICES),ordered=True, \
1392                                                             initialize=dict([(args,LOG_poly_init(args)) for args in getattr(self._model(),LOG_VERTICES)])))
1393        LOG_LAMBDA_SET = self._update_dicts('LOG_LAMBDA_SET',Set(ordered=True,dimen=self._index_dimen+1,initialize=LOG_lamba_set_init))
1394        LOG_BRANCHING_SCHEME = self._update_dicts('LOG_BRANCHING_SCHEME',Set(ordered=True,dimen=self._index_dimen+1,initialize=LOG_BRANCHING_SCHEME_init))
1395        LOG_BRANCHING_LEFT = self._update_dicts('LOG_BRANCHING_LEFT',Set(getattr(self._model(),LOG_BRANCHING_SCHEME),ordered=True, \
1396                                                                         initialize=dict([(args,LOG_BRANCHING_LEFT_init(args)) for args in getattr(self._model(),LOG_BRANCHING_SCHEME)])))
1397        LOG_BRANCHING_RIGHT = self._update_dicts('LOG_BRANCHING_RIGHT',Set(getattr(self._model(),LOG_BRANCHING_SCHEME),ordered=True,\
1398                                                                           initialize=dict([(args,LOG_BRANCHING_RIGHT_init(args)) for args in getattr(self._model(),LOG_BRANCHING_SCHEME)])))
1399       
1400        lmda = self._update_dicts('LOG_LAMBDA',Var(getattr(self._model(),LOG_LAMBDA_SET),within=NonNegativeReals))
1401        y = self._update_dicts('LOG_BIN_y',Var(getattr(self._model(),LOG_BRANCHING_SCHEME),within=Binary))
1402
1403        def LOG_LINEAR_constraint1_rule(model,*t):
1404            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,lmda,t,v)*self._pw_pts(t)[v-1] for v in self._getset(model,LOG_vert,t))
1405        def LOG_LINEAR_constraint2_rule(model,*t):
1406            LHS = self._getvar(model,self.YVAR.name,t)
1407            RHS = sum(self._getvar(model,lmda,t,v)*self._F(model,t,self._pw_pts(t)[v-1]) for v in self._getset(model,LOG_vert,t))
1408            if self.PW_TYPE == 'UB':
1409                return LHS <= RHS
1410            elif self.PW_TYPE == 'LB':
1411                return LHS >= RHS
1412            elif self.PW_TYPE == 'EQ':
1413                return LHS == RHS
1414        def LOG_LINEAR_constraint3_rule(model,*t):
1415            return 1 == sum(self._getvar(model,lmda,t,v) for v in self._getset(model,LOG_vert,t))
1416        def LOG_LINEAR_constraint4_rule(model,*args):
1417            t = args[:-1]
1418            s = args[-1]
1419            return sum(self._getvar(model,lmda,t,v) for v in self._getset(model,LOG_BRANCHING_LEFT,t,s)) <= self._getvar(model,y,t,s)
1420        def LOG_LINEAR_constraint5_rule(model,*args):
1421            t = args[:-1]
1422            s = args[-1]
1423            return sum(self._getvar(model,lmda,t,v) for v in self._getset(model,LOG_BRANCHING_RIGHT,t,s)) <= (1-self._getvar(model,y,t,s))
1424
1425        if self._element_mode:
1426            self._update_dicts('LOG_constraint1',Constraint(rule=LOG_LINEAR_constraint1_rule))
1427            self._update_dicts('LOG_constraint2',Constraint(rule=LOG_LINEAR_constraint2_rule))
1428            self._update_dicts('LOG_constraint3',Constraint(rule=LOG_LINEAR_constraint3_rule))
1429            self._update_dicts('LOG_constraint4',Constraint(getattr(self._model(),LOG_BRANCHING_SCHEME),rule=LOG_LINEAR_constraint4_rule))
1430            self._update_dicts('LOG_constraint5',Constraint(getattr(self._model(),LOG_BRANCHING_SCHEME),rule=LOG_LINEAR_constraint5_rule))
1431        else:
1432            self._update_dicts('LOG_constraint1',Constraint(LOG_SET,rule=LOG_LINEAR_constraint1_rule))
1433            self._update_dicts('LOG_constraint2',Constraint(LOG_SET,rule=LOG_LINEAR_constraint2_rule))
1434            self._update_dicts('LOG_constraint3',Constraint(LOG_SET,rule=LOG_LINEAR_constraint3_rule))
1435            self._update_dicts('LOG_constraint4',Constraint(getattr(self._model(),LOG_BRANCHING_SCHEME),rule=LOG_LINEAR_constraint4_rule))
1436            self._update_dicts('LOG_constraint5',Constraint(getattr(self._model(),LOG_BRANCHING_SCHEME),rule=LOG_LINEAR_constraint5_rule))
1437
1438    def _construct_MC(self,MC_SET):
1439       
1440        SLOPE = {}
1441        INTERSEPT = {}
1442        if self._element_mode:
1443            for i in xrange(1,len(self._pw_pts())):
1444                SLOPE[i] = (self._F(self._model(),self._pw_pts()[i])-self._F(self._model(),self._pw_pts()[i-1]))/(self._pw_pts()[i]-self._pw_pts()[i-1])
1445                INTERSEPT[i] = self._F(self._model(),self._pw_pts()[i-1]) - (SLOPE[i]*self._pw_pts()[i-1])
1446        else:
1447            for t in MC_SET:
1448                SLOPE[t] = {}
1449                INTERSEPT[t] = {}
1450                for i in xrange(1,len(self._pw_pts(t))):
1451                    SLOPE[t][i] = (self._F(self._model(),t,self._pw_pts(t)[i])-self._F(self._model(),t,self._pw_pts(t)[i-1]))/(self._pw_pts(t)[i]-self._pw_pts(t)[i-1])
1452                    INTERSEPT[t][i] = self._F(self._model(),t,self._pw_pts(t)[i-1]) - (SLOPE[t][i]*self._pw_pts(t)[i-1])
1453                SLOPE[(t,)] = SLOPE[t]
1454                INTERSEPT[(t,)] = INTERSEPT[t]
1455       
1456        def MC_POLYTOPES_init(model):
1457            if self._element_mode:
1458                return (i for i in xrange(1,len(self._pw_pts())))
1459            else:
1460                return (t+(i,) for t in MC_SET.tuplize() \
1461                              for i in xrange(1,len(self._pw_pts(t))))
1462        def MC_poly_init(*t):
1463            return (i for i in xrange(1,len(self._pw_pts(t))))
1464       
1465        MC_POLYTOPES = self._update_dicts('MC_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=MC_POLYTOPES_init))
1466        if self._element_mode:
1467            MC_poly = self._update_dicts('MC_poly',Set(ordered=True, initialize=MC_poly_init()))
1468        else:
1469            MC_poly = self._update_dicts('MC_poly',Set(MC_SET,ordered=True, \
1470                                                       initialize=dict([(t,MC_poly_init(t)) for t in MC_SET])))
1471       
1472        xp = self._update_dicts('MC_X_POLY',Var(getattr(self._model(),MC_POLYTOPES)))
1473        y = self._update_dicts('MC_BIN_y',Var(getattr(self._model(),MC_POLYTOPES),within=Binary))
1474       
1475        def MC_LINEAR_constraint1_rule(model,*t):
1476            return self._getvar(model,self.XVAR.name,t) == sum(self._getvar(model,xp,t,p) for p in self._getset(model,MC_poly,t))
1477        def MC_LINEAR_constraint2_rule(model,*t):
1478            LHS = self._getvar(model,self.YVAR.name,t)
1479            RHS = 0
1480            if self._element_mode:
1481                RHS += sum((self._getvar(model,xp,t,p)*SLOPE[p])+(self._getvar(model,y,t,p)*INTERSEPT[p]) for p in self._getset(model,MC_poly,t))
1482            else:
1483                RHS += sum((self._getvar(model,xp,t,p)*SLOPE[t][p])+(self._getvar(model,y,t,p)*INTERSEPT[t][p]) for p in self._getset(model,MC_poly,t))
1484            if self.PW_TYPE == 'UB':
1485                return LHS <= RHS
1486            elif self.PW_TYPE == 'LB':
1487                return LHS >= RHS
1488            elif self.PW_TYPE == 'EQ':
1489                return LHS == RHS
1490        def MC_LINEAR_constraint3_rule(model,*args):
1491            t = args[:-1]
1492            p = args[-1]
1493            return self._getvar(model,y,t,p)*self._pw_pts(t)[p-1] <= self._getvar(model,xp,t,p)
1494        def MC_LINEAR_constraint4_rule(model,*args):
1495            t = args[:-1]
1496            p = args[-1]
1497            return self._getvar(model,xp,t,p)  <= self._getvar(model,y,t,p)*self._pw_pts(t)[p]
1498        def MC_LINEAR_constraint5_rule(model,*t):
1499            return sum(self._getvar(model,y,t,p) for p in self._getset(model,MC_poly,t)) == 1
1500           
1501        if self._element_mode:
1502            self._update_dicts('MC_constraint1',Constraint(rule=MC_LINEAR_constraint1_rule))
1503            self._update_dicts('MC_constraint2',Constraint(rule=MC_LINEAR_constraint2_rule))
1504            self._update_dicts('MC_constraint3',Constraint(getattr(self._model(),MC_POLYTOPES),rule=MC_LINEAR_constraint3_rule))
1505            self._update_dicts('MC_constraint4',Constraint(getattr(self._model(),MC_POLYTOPES),rule=MC_LINEAR_constraint4_rule))           
1506            self._update_dicts('MC_constraint5',Constraint(rule=MC_LINEAR_constraint5_rule))
1507        else:
1508            self._update_dicts('MC_constraint1',Constraint(MC_SET,rule=MC_LINEAR_constraint1_rule))
1509            self._update_dicts('MC_constraint2',Constraint(MC_SET,rule=MC_LINEAR_constraint2_rule))
1510            self._update_dicts('MC_constraint3',Constraint(getattr(self._model(),MC_POLYTOPES),rule=MC_LINEAR_constraint3_rule))
1511            self._update_dicts('MC_constraint4',Constraint(getattr(self._model(),MC_POLYTOPES),rule=MC_LINEAR_constraint4_rule))           
1512            self._update_dicts('MC_constraint5',Constraint(MC_SET,rule=MC_LINEAR_constraint5_rule))
1513
1514    def _construct_INC(self,INC_SET):
1515       
1516        def INC_POLYTOPES_init(model):
1517            if self._element_mode:
1518                return (i for i in xrange(1,len(self._pw_pts())))
1519            else:
1520                return (t+(i,) for t in INC_SET.tuplize() \
1521                              for i in xrange(1,len(self._pw_pts(t))))
1522        def INC_poly_init(*t):
1523            return (i for i in xrange(1,len(self._pw_pts(t))))
1524        def INC_VERTICES_init(model):
1525            if self._element_mode:
1526                return (i for i in xrange(1,len(self._pw_pts())+1))
1527            else:
1528                return (t+(i,) for t in INC_SET.tuplize() \
1529                              for i in xrange(1,len(self._pw_pts(t))+1))
1530        def INC_Y_init(model):
1531            if self._element_mode:
1532                return (i for i in xrange(1,len(self._pw_pts())-1))
1533            else:
1534                return (t+(i,) for t in INC_SET.tuplize() \
1535                              for i in xrange(1,len(self._pw_pts(t))-1))
1536       
1537        INC_POLYTOPES = self._update_dicts('INC_POLYTOPES',Set(ordered=True,dimen=self._index_dimen+1,initialize=INC_POLYTOPES_init))
1538        if self._element_mode:
1539            INC_poly = self._update_dicts('INC_poly',Set(ordered=True, initialize=INC_poly_init()))
1540        else:
1541            INC_poly = self._update_dicts('INC_poly',Set(INC_SET,ordered=True, \
1542                                                         initialize=dict([(t,INC_poly_init(t)) for t in INC_SET])))
1543        INC_VERTICES = self._update_dicts('INC_VERTICES',Set(ordered=True,dimen=self._index_dimen+1,initialize=INC_VERTICES_init))
1544        INC_YSET = self._update_dicts('INC_YSET',Set(ordered=True,dimen=self._index_dimen+1,initialize=INC_Y_init))
1545       
1546        delta = self._update_dicts('INC_DELTA',Var(getattr(self._model(),INC_POLYTOPES)))
1547        if self._element_mode:
1548            self._getvar(self._model(),delta,1).setub(1)
1549            self._getvar(self._model(),delta,len(self._pw_pts())-1).setlb(0)
1550        else:
1551            for t in INC_SET:
1552                self._getvar(self._model(),delta,t,1).setub(1)
1553                self._getvar(self._model(),delta,t,len(self._pw_pts(t))-1).setlb(0)
1554        y = self._update_dicts('INC_BIN_y',Var(getattr(self._model(),INC_YSET),within=Binary))
1555       
1556        def INC_LINEAR_constraint1_rule(model,*t):
1557            return self._getvar(model,self.XVAR.name,t) == self._pw_pts(t)[0] + sum(self._getvar(model,delta,t,p)*(self._pw_pts(t)[p]-self._pw_pts(t)[p-1]) for p in self._getset(model,INC_poly,t))
1558        def INC_LINEAR_constraint2_rule(model,*t):
1559            LHS = self._getvar(model,self.YVAR.name,t)
1560            RHS = self._F(model,t,self._pw_pts(t)[0]) + sum(self._getvar(model,delta,t,p)*(self._F(model,t,self._pw_pts(t)[p])-self._F(model,t,self._pw_pts(t)[p-1])) for p in self._getset(model,INC_poly,t))
1561            if self.PW_TYPE == 'UB':
1562                return LHS <= RHS
1563            elif self.PW_TYPE == 'LB':
1564                return LHS >= RHS
1565            elif self.PW_TYPE == 'EQ':
1566                return LHS == RHS
1567        def INC_LINEAR_constraint3_rule(model,*args):
1568            t = args[:-1]
1569            p = args[-1]
1570            if p != self._getset(model,INC_poly,t).last():
1571                return self._getvar(model,delta,t,p+1) <= self._getvar(model,y,t,p)
1572            else:
1573                return Constraint.Skip
1574        def INC_LINEAR_constraint4_rule(model,*args):
1575            t = args[:-1]
1576            p = args[-1]
1577            if p != self._getset(model,INC_poly,t).last():
1578                return self._getvar(model,y,t,p) <= self._getvar(model,delta,t,p)
1579            else:
1580                return Constraint.Skip
1581       
1582        if self._element_mode:
1583            self._update_dicts('INC_constraint1',Constraint(rule=INC_LINEAR_constraint1_rule))
1584            self._update_dicts('INC_constraint2',Constraint(rule=INC_LINEAR_constraint2_rule))
1585            self._update_dicts('INC_constraint3',Constraint(getattr(self._model(),INC_POLYTOPES),rule=INC_LINEAR_constraint3_rule))
1586            self._update_dicts('INC_constraint4',Constraint(getattr(self._model(),INC_POLYTOPES),rule=INC_LINEAR_constraint4_rule))
1587        else:
1588            self._update_dicts('INC_constraint1',Constraint(INC_SET,rule=INC_LINEAR_constraint1_rule))
1589            self._update_dicts('INC_constraint2',Constraint(INC_SET,rule=INC_LINEAR_constraint2_rule))
1590            self._update_dicts('INC_constraint3',Constraint(getattr(self._model(),INC_POLYTOPES),rule=INC_LINEAR_constraint3_rule))
1591            self._update_dicts('INC_constraint4',Constraint(getattr(self._model(),INC_POLYTOPES),rule=INC_LINEAR_constraint4_rule))
1592
1593   
1594   
1595    ###### Find tightest possible BIGM values for convex/concave BIG constraints
1596    def _M_func(self,a,b,c,*args):
1597        tpl = flatten_tuple((self._model(),)+args)
1598        return self._F(*(tpl+(a,))) - self._F(*(tpl+(b,))) - ((a-b) * (float(self._F(*(tpl+(c,)))-self._F(*(tpl+(b,)))) / float(c-b)))
1599
1600    def _find_M(self,PTS,BOUND):
1601       
1602        M_final = {}
1603        if self._element_mode:
1604            for j in xrange(1,len(PTS)):
1605                index = j
1606                if (BOUND == 'LB'):
1607                    M_final[index] = min( [0.0, min([self._M_func(PTS[k],PTS[j-1],PTS[j]) for k in xrange(len(PTS))])] )
1608                elif (BOUND == 'UB'):
1609                    M_final[index] = max( [0.0, max([self._M_func(PTS[k],PTS[j-1],PTS[j]) for k in xrange(len(PTS))])] )
1610                if M_final[index] == 0.0:
1611                    del M_final[index]
1612        else:
1613            for t in PTS.keys():
1614                for j in xrange(1,len(PTS[t])):
1615                    index = flatten_tuple((t,j))
1616                    if (BOUND == 'LB'):
1617                        M_final[index] = min( [0.0, min([self._M_func(PTS[t][k],PTS[t][j-1],PTS[t][j],t) for k in xrange(len(PTS[t]))])] )
1618                    elif (BOUND == 'UB'):
1619                        M_final[index] = max( [0.0, max([self._M_func(PTS[t][k],PTS[t][j-1],PTS[t][j],t) for k in xrange(len(PTS[t]))])] )
1620                    if M_final[index] == 0.0:
1621                        del M_final[index]
1622        return M_final
1623    ########
Note: See TracBrowser for help on using the repository browser.