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

Revision 5462, 71.4 KB checked in by gabeh, 2 years ago (diff)

Updating the required syntax for f_rule and d_f_rule keywords passed to the Piecewise component. Also updating the piecewise tests and examples directory to comply with the new syntax. Any other tests (if any) which make use of this component may need to be updated (only if using an _indexed_ Piecewise component). The new syntax allows for more flexibility in defining the Piecewise component over an index, e.g, model.C = Piecewise(model.INDEX, ...). All function rules supplied to a Piecewise component having an index set must now accept as arguments each of the indices of the set as well as a pyomo model and domain point to evaluate the function, e.g.

def f(model, i, j, x)

...

model.C = Piecewise([(1,1),(2,2)], model.y, model.x, f_rule=f)

Here i,j are the indices and x is the domain point. The domain point should always be the last argument. This makes the look of Piecewise function rules more agreeable with the current syntax of Pyomo Constraint rules, and allows the modeler to approximate different functions with a single indexed Piecewise component using the indices passed into the function rule.

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