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

Revision 5597, 71.5 KB checked in by wehart, 2 years ago (diff)

A first code review, focusing on Component and
IndexedComponent?. There were many changes because I
changed the names of some component attributes to clarify what
were public and private attributes.

A private attribute refers to an attribute that an end-user is
not expected to use.

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