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

Revision 4744, 70.1 KB checked in by wehart, 3 years ago (diff)

Documentation update.

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