Changeset 2882


Ignore:
Timestamp:
Jul 29, 2010 1:35:45 PM (9 years ago)
Author:
claird
Message:

Added nonlinear support to pyomo. This includes new intrinsic functions and new NL writer as well as corresponding tests

Location:
coopr.pyomo/trunk/coopr/pyomo
Files:
390 added
5 edited

Legend:

Unmodified
Added
Removed
  • coopr.pyomo/trunk/coopr/pyomo/base/__init__.py

    r2865 r2882  
    1111from numvalue import *
    1212from expr import *
     13from intrinsic_functions import *
    1314from numtypes import *
    1415from plugin import *
  • coopr.pyomo/trunk/coopr/pyomo/base/intrinsic_functions.py

    r2400 r2882  
    1515#        intrinsic function with the value of that expression?
    1616# TODO:  Do we need to register these expression objects?
    17 # TODO:  Verify why we need to allow Components as arguments.
    1817#
    1918
    20 __all__ = ['log']
     19__all__ = ['log', 'log10', 'sin', 'cos', 'tan', 'cosh', 'sinh', 'tanh', 'asin', 'acos', 'atan', 'exp', 'sqrt', 'asinh', 'acosh', 'atanh']
    2120
    2221import component
     
    3332    def __init__(self, args=(), name=None, operation=None):
    3433        expr.Expression.__init__(self, args=args, nargs=len(args), name=name, operation=operation, tuple_op=True)
    35 
     34        self.coef = 1.0
    3635
    3736def expr_args(*args):
     
    5150       raise "TypeError: log() takes exactly 1 argument (%d given)" % len(args)
    5251    if expr_args(*args):
    53         return _IntrinsicFuncExpression(args=args, name='log', operation=math.log)
     52        tmp = list()
     53        tmp.append(args[0])
     54        return _IntrinsicFuncExpression(args=tmp, name='log', operation=math.log)
    5455    return math.log(*args)
    5556
     57def log10(*args):
     58    if len(args) != 1:
     59       raise "TypeError: log10() takes exactly 1 argument (%d given)" % len(args)
     60    if expr_args(*args):
     61        tmp = list()
     62        tmp.append(args[0])
     63        return _IntrinsicFuncExpression(args=tmp, name='log10', operation=math.log10)
     64    return math.log10(*args)
     65
     66def sin(*args):
     67    if len(args) != 1:
     68       raise "TypeError: sin() takes exactly 1 argument (%d given)" % len(args)
     69    if expr_args(*args):
     70        tmp = list()
     71        tmp.append(args[0])
     72        return _IntrinsicFuncExpression(args=tmp, name='sin', operation=math.sin)
     73    return math.sin(*args)
     74
     75def cos(*args):
     76    if len(args) != 1:
     77       raise "TypeError: cos() takes exactly 1 argument (%d given)" % len(args)
     78    if expr_args(*args):
     79        tmp = list()
     80        tmp.append(args[0])
     81        return _IntrinsicFuncExpression(args=tmp, name='cos', operation=math.cos)
     82    return math.cos(*args)
     83
     84def tan(*args):
     85    if len(args) != 1:
     86       raise "TypeError: tan() takes exactly 1 argument (%d given)" % len(args)
     87    if expr_args(*args):
     88        tmp = list()
     89        tmp.append(args[0])
     90        return _IntrinsicFuncExpression(args=tmp, name='tan', operation=math.tan)
     91    return math.tan(*args)
     92
     93def sinh(*args):
     94    if len(args) != 1:
     95       raise "TypeError: sinh() takes exactly 1 argument (%d given)" % len(args)
     96    if expr_args(*args):
     97        tmp = list()
     98        tmp.append(args[0])
     99        return _IntrinsicFuncExpression(args=tmp, name='sinh', operation=math.sinh)
     100    return math.sinh(*args)
     101
     102def cosh(*args):
     103    if len(args) != 1:
     104       raise "TypeError: cosh() takes exactly 1 argument (%d given)" % len(args)
     105    if expr_args(*args):
     106        tmp = list()
     107        tmp.append(args[0])
     108        return _IntrinsicFuncExpression(args=tmp, name='cosh', operation=math.cosh)
     109    return math.cosh(*args)
     110
     111def tanh(*args):
     112    if len(args) != 1:
     113       raise "TypeError: tanh() takes exactly 1 argument (%d given)" % len(args)
     114    if expr_args(*args):
     115        tmp = list()
     116        tmp.append(args[0])
     117        return _IntrinsicFuncExpression(args=tmp, name='tanh', operation=math.tanh)
     118    return math.tanh(*args)
     119
     120def asin(*args):
     121    if len(args) != 1:
     122       raise "TypeError: asin() takes exactly 1 argument (%d given)" % len(args)
     123    if expr_args(*args):
     124        tmp = list()
     125        tmp.append(args[0])
     126        return _IntrinsicFuncExpression(args=tmp, name='asin', operation=math.asin)
     127    return math.asin(*args)
     128
     129def acos(*args):
     130    if len(args) != 1:
     131       raise "TypeError: acos() takes exactly 1 argument (%d given)" % len(args)
     132    if expr_args(*args):
     133        tmp = list()
     134        tmp.append(args[0])
     135        return _IntrinsicFuncExpression(args=tmp, name='acos', operation=math.acos)
     136    return math.acos(*args)
     137
     138def atan(*args):
     139    if len(args) != 1:
     140       raise "TypeError: atan() takes exactly 1 argument (%d given)" % len(args)
     141    if expr_args(*args):
     142        tmp = list()
     143        tmp.append(args[0])
     144        return _IntrinsicFuncExpression(args=tmp, name='atan', operation=math.atan)
     145    return math.atan(*args)
     146
     147def exp(*args):
     148    if len(args) != 1:
     149       raise "TypeError: exp() takes exactly 1 argument (%d given)" % len(args)
     150    if expr_args(*args):
     151        tmp = list()
     152        tmp.append(args[0])
     153        return _IntrinsicFuncExpression(args=tmp, name='exp', operation=math.exp)
     154    return math.exp(*args)
     155
     156def sqrt(*args):
     157    if len(args) != 1:
     158       raise "TypeError: sqrt() takes exactly 1 argument (%d given)" % len(args)
     159    if expr_args(*args):
     160        tmp = list()
     161        tmp.append(args[0])
     162        return _IntrinsicFuncExpression(args=tmp, name='sqrt', operation=math.sqrt)
     163    return math.sqrt(*args)
     164
     165def asinh(*args):
     166    if len(args) != 1:
     167       raise "TypeError: atan() takes exactly 1 argument (%d given)" % len(args)
     168    if expr_args(*args):
     169        tmp = list()
     170        tmp.append(args[0])
     171        return _IntrinsicFuncExpression(args=tmp, name='asinh', operation=math.asinh)
     172    return math.asinh(*args)
     173
     174def acosh(*args):
     175    if len(args) != 1:
     176       raise "TypeError: atan() takes exactly 1 argument (%d given)" % len(args)
     177    if expr_args(*args):
     178        tmp = list()
     179        tmp.append(args[0])
     180        return _IntrinsicFuncExpression(args=tmp, name='acosh', operation=math.acosh)
     181    return math.acosh(*args)
     182
     183def atanh(*args):
     184    if len(args) != 1:
     185       raise "TypeError: atan() takes exactly 1 argument (%d given)" % len(args)
     186    if expr_args(*args):
     187        tmp = list()
     188        tmp.append(args[0])
     189        return _IntrinsicFuncExpression(args=tmp, name='atanh', operation=math.atanh)
     190    return math.atanh(*args)
  • coopr.pyomo/trunk/coopr/pyomo/io/ampl.py

    r2763 r2882  
    1616from coopr.pyomo.base.numtypes import minimize, maximize
    1717from coopr.pyomo.base import *
     18from coopr.pyomo.base import expr
     19from coopr.pyomo.base import intrinsic_functions
     20from coopr.pyomo.base.var import _VarValue, Var, _VarElement, _VarBase
     21from coopr.pyomo.base import numvalue
    1822from coopr.opt.base import *
    1923from coopr.pyomo.base.param import _ParamValue
    20 from coopr.pyomo.base.var import _VarBase
     24from coopr.pyomo.base import var
     25from coopr.pyomo.base import param
    2126from pyutilib.component.core import alias
     27from coopr.pyomo.expr.canonical_repn import is_nonlinear, is_linear
    2228
    2329class ProblemWriter_nl(AbstractProblemWriter):
     
    4652           raise ValueError, "ERROR: nonconstant bound: " + str(exp)
    4753           return None
    48 
     54   
     55    def _print_nonlinear_terms_NL(self, OUTPUT, exp):
     56        if isinstance(exp, intrinsic_functions._IntrinsicFuncExpression):
     57            if exp.name == 'log':
     58                print >>OUTPUT, "o43", " #log"
     59            elif exp.name == 'log10':
     60                print >>OUTPUT, "o42", " #log10"
     61            elif exp.name == 'sin':
     62                print >>OUTPUT, "o41", " #sin"
     63            elif exp.name == 'cos':
     64                print >>OUTPUT, "o46", " #cos"
     65            elif exp.name == 'tan':
     66                print >>OUTPUT, "o38", " #tan"
     67            elif exp.name == 'sinh':
     68                print >>OUTPUT, "o40", " #sinh"
     69            elif exp.name == 'cosh':
     70                print >>OUTPUT, "o45", " #cosh"
     71            elif exp.name == 'tanh':
     72                print >>OUTPUT, "o37", " #tanh"
     73            elif exp.name == 'asin':
     74                print >>OUTPUT, "o51", " #asin"
     75            elif exp.name == 'acos':
     76                print >>OUTPUT, "o53", " #acos"
     77            elif exp.name == 'atan':
     78                print >>OUTPUT, "o49", " #atan"
     79            elif exp.name == 'exp':
     80                print >>OUTPUT, "o44", " #exp"
     81            elif exp.name == 'sqrt':
     82                print >>OUTPUT, "o39", " #sqrt"
     83            elif exp.name == 'asinh':
     84                print >>OUTPUT, "o50", " #asinh"
     85            elif exp.name == 'acosh':
     86                print >>OUTPUT, "o52", " #acosh"
     87            elif exp.name == 'atanh':
     88                print >>OUTPUT, "o47", " #atanh"
     89               
     90            for child_exp in exp._args:
     91                self._print_nonlinear_terms_NL(OUTPUT, child_exp)
     92       
     93        elif isinstance(exp, expr._SumExpression):
     94            n = len(exp._args)
     95            if exp._const != 0.0:
     96                print >>OUTPUT, "o0", " #+"
     97                print >>OUTPUT, "n" + str(exp._const)
     98               
     99            for i in range(0,n-1):
     100                print >>OUTPUT, "o0", " #+"
     101                if exp._coef[i] != 1:
     102                    print >>OUTPUT, "o2", " #*"
     103                    print >>OUTPUT, "n" + str(exp._coef[i])
     104                self._print_nonlinear_terms_NL(OUTPUT, exp._args[i])
     105
     106            if exp._coef[n-1] != 1:
     107                print >>OUTPUT, "o2", " #*"
     108                print >>OUTPUT, "n" + str(exp._coef[n-1])
     109            self._print_nonlinear_terms_NL(OUTPUT, exp._args[n-1])
     110
     111        elif isinstance(exp, list):
     112            # this is an implied summation of expressions (we did not create a new sum expression for efficiency)
     113            # this should be a list of tuples where [0] is the coeff and [1] is the expr to write
     114            n = len(exp)
     115            for i in range(0,n):
     116                assert(isinstance(exp[i],tuple))
     117                coef = exp[i][0]
     118                child_exp = exp[i][1]
     119                if i != n-1:
     120                    # need the + op if it is not the last entry in the list
     121                    print >>OUTPUT, "o0", " #+"
     122                   
     123                if coef != 1.0:
     124                    print >>OUTPUT, "o2", " #*"
     125                    print >>OUTPUT, "n" + str(coef)
     126                self._print_nonlinear_terms_NL(OUTPUT, child_exp)
     127
     128        elif isinstance(exp, expr._ProductExpression):
     129            denom_exists = False
     130            if len(exp._denominator) == 0:
     131                pass
     132            else:
     133                print >>OUTPUT, "o3 # /"
     134                denom_exists = True
     135                if len(exp._numerator) == 0:
     136                    print >>OUTPUT, "n1"
     137           
     138            if exp.coef != 1:
     139                print >>OUTPUT, "o2 # *"
     140                print >>OUTPUT, "n" + str(exp.coef)
     141
     142            # print out the numerator
     143            child_counter = 0
     144            max_count = len(exp._numerator)-1
     145            for child_exp in exp._numerator:
     146#                if child_exp != exp._numerator[-1]:
     147                if child_counter < max_count:
     148                    print >>OUTPUT, "o2", " #*"
     149                self._print_nonlinear_terms_NL(OUTPUT, child_exp)
     150                child_counter += 1
     151
     152            if denom_exists:
     153            # print out the denominator           
     154                child_counter = 0
     155                max_count = len(exp._denominator)-1
     156                for child_exp in exp._denominator:
     157                    if child_counter < max_count:
     158                        print >>OUTPUT, "o2", " #*"
     159                    self._print_nonlinear_terms_NL(OUTPUT, child_exp)
     160                    child_counter += 1
     161
     162        elif isinstance(exp, expr._PowExpression):
     163            print >>OUTPUT, "o5", " #^"
     164            for child_exp in exp._args:
     165                self._print_nonlinear_terms_NL(OUTPUT, child_exp)
     166
     167        elif isinstance(exp,var._VarValue) or isinstance(exp,var.Var):
     168            print >>OUTPUT, "v" + str(exp.ampl_var_id) + " #" + str(exp.label)
     169           
     170        elif isinstance(exp,param._ParamValue):
     171            print >>OUTPUT, "n" + str(exp.value) + " #" + str(exp.name)
     172
     173        elif isinstance(exp,numvalue.NumericConstant):
     174            print >>OUTPUT, "n" + str(exp.value) + " # numeric constant"
     175
     176        else:
     177            raise ValueError, "Unsupported expression type in _print_nonlinear_terms_NL"
     178           
    49179    def _print_model_NL(self, model, OUTPUT, verbose=False):
    50180
     
    58188        # Collect statistics
    59189        #
    60         nc = no = neqn = nrange = 0
    61190        Obj = model.active_components(Objective)
    62191        Con = model.active_components(Constraint)
    63 
    64         # Dictionary of all variables used in all expressions in this model
    65         Var = {}
    66 
    67         #
    68         # Count number of objectives
    69         #
     192        Vars = dict() # will be the entire list of all vars used in the problem
     193        #
     194        # Count number of objectives and build the ampl_repns
     195        #
     196        n_objs = 0
     197        n_nonlinear_objs = 0
     198        ObjVars = set()
     199        ObjNonlinearVars = set()
    70200        for obj in Obj:
    71201            if Obj[obj].trivial:
    72202                continue
    73203            for i in Obj[obj]:
    74                 if not is_constant(Obj[obj][i].repn):
    75                     # Merge variables from this expression
    76                     Var.update( Obj[obj][i].repn[-1] )
    77                     no += 1
    78         #
    79         # Count number of constraints
    80         #
     204                ampl_repn = generate_ampl_repn(Obj[obj][i].expr)
     205#                ampl_repn._nonlinear_expr.pprint()
     206                Obj[obj][i].ampl_repn = ampl_repn
     207               
     208                Vars.update(ampl_repn._linear_terms_var)
     209                Vars.update(ampl_repn._nonlinear_vars)
     210               
     211                ObjVars.update(ampl_repn._linear_terms_var.keys())
     212                ObjVars.update(ampl_repn._nonlinear_vars.keys())
     213                ObjNonlinearVars.update(ampl_repn._nonlinear_vars.keys())
     214                n_objs += 1
     215                if ampl_repn.is_nonlinear():
     216                    n_nonlinear_objs += 1
     217        assert(n_objs == 1)
     218
     219        #
     220        # Count number of constraints and build the ampl_repns
     221        #
     222        n_ranges = 0
     223        n_single_sided_ineq = 0
     224        n_equals = 0
     225        n_nonlinear_constraints = 0
     226        ConVars = set()
     227        ConNonlinearVars = set()
     228        nnz_grad_constraints = 0
    81229        for con in Con:
    82230            C = Con[con]
     
    84232                continue
    85233            for i in C:
    86                 if is_constant(C[i].repn):
     234                ampl_repn = generate_ampl_repn(C[i].body)
     235                C[i].ampl_repn = ampl_repn
     236                if is_constant(ampl_repn):
    87237                    continue
     238
    88239                # Merge variables from this expression
    89                 Var.update( C[i].repn[-1] )
    90                 nc += 1
    91                 if C[i].lower is not None:
    92                     L = -inf
    93                     U = inf
    94                     if C[i].lower is not None: L = self._get_bound(C[i].lower)
    95                     if C[i].upper is not None: U = self._get_bound(C[i].upper)
    96                     if (L == U):
    97                         neqn += 1
    98                 elif L > U:
     240                all_vars = ampl_repn._linear_terms_var.copy()
     241                all_vars.update(ampl_repn._nonlinear_vars)
     242                Vars.update(all_vars)
     243                nnz_grad_constraints += len(all_vars)
     244                ConVars.update(all_vars.keys())
     245                ConNonlinearVars.update(ampl_repn._nonlinear_vars.keys())
     246                if ampl_repn.is_nonlinear():
     247                    n_nonlinear_constraints += 1
     248
     249                L = None
     250                U = None
     251                if C[i].lower is not None: L = self._get_bound(C[i].lower)
     252                if C[i].upper is not None: U = self._get_bound(C[i].upper)
     253                if (L is None and U is None):
     254                    # constraint with NO bounds ???
     255                    raise ValueError, "Constraint " + str(con) +\
     256                                         "[" + str(i) + "]: lower bound = None and upper bound = None"
     257                elif (L is None or U is None):
     258                    # one sided inequality
     259                    n_single_sided_ineq += 1
     260                elif (L > U):
    99261                    msg = 'Constraint %s[%s]: lower bound greater than upper' \
    100262                          ' bound (%s > %s)'
    101263                    raise ValueError, msg % (str(con), str(i), str(L), str(U) )
    102                 elif C[i].lower is not None and C[i].upper is not None:
    103                     nrange += 1
    104         #
    105         niv = nbv = 0
     264                elif (L == U):
     265                    # equality
     266                    n_equals += 1
     267                else:
     268                    # double sided inequality
     269                    # both are not none and they are valid
     270                    n_ranges += 1
     271       
    106272        Ilist = []
    107273        Blist = []
    108274        Vlist = []
    109         keys = Var.keys()
     275        keys = Vars.keys()
    110276        keys.sort()
    111277        for ndx in keys:
    112             var = Var[ndx].var
     278            var = Vars[ndx].var
    113279            if isinstance(var.domain, IntegerSet):
    114280                Ilist.append(ndx)
    115281            elif isinstance(var.domain, BooleanSet):
    116                 L = Var[ndx].lb
    117                 U = Var[ndx].ub
    118                 if L is not None and value(L) == 0 \
    119                     and U is not None and value(U) == 1:
     282                L = Vars[ndx].lb
     283                U = Vars[ndx].ub
     284                if L is None or U is None:
     285                    raise ValueError, "Variable " + str(var.label) +\
     286                          "is binary, but does not have lb and ub set"
     287                elif value(L) == 0 and value(U) == 1:
    120288                    Blist.append(ndx)
    121289                else:
     
    123291            else:
    124292                Vlist.append(ndx)
    125         for ndx in Blist:
    126             nbv += 1
    127             Vlist.append(ndx)
    128         for ndxx in Ilist:
    129             niv += 1
    130             Vlist.append(ndx)
    131         #
    132         vctr=0
    133         var_id = {}
    134         for ndx in Vlist:
    135             var_id[ndx] = vctr
    136             vctr += 1
     293
     294        Vlist.extend(Blist)
     295        Vlist.extend(Ilist)
     296
     297        # put the ampl variable id into the variable
     298        var_id_ctr=0
     299        for var_name in Vlist:
     300            Vars[var_name].ampl_var_id = var_id_ctr
     301            var_id_ctr += 1
     302
    137303        #
    138304        # Print Header
     
    144310        # LINE 2
    145311        #
    146         print >>OUTPUT, " ",len(Var), nc, no, nrange, str(neqn) +\
    147                 "\t# vars, constraints, objectives, ranges, eqns"
     312        print >>OUTPUT, "", len(Vars), str(n_single_sided_ineq+n_ranges+n_equals), str(n_objs), str(n_ranges), str(n_equals) +\
     313                "\t# vars, constraints, objectives, general inequalities, equalities"
     314       
    148315        #
    149316        # LINE 3
    150317        #
    151         print >>OUTPUT, " 0 0\t# nonlinear constraints, objectives"
     318        print >>OUTPUT, " " + str(n_nonlinear_constraints),  str(n_nonlinear_objs) + "\t# nonlinear constraints, objectives"
    152319        #
    153320        # LINE 4
     
    157324        # LINE 5
    158325        #
    159         print >>OUTPUT, ' 0 0 0\t# nonlinear vars in constraints, '           \
    160                         'objectives, both'
     326        Nonlinear_Vars_in_Objs_and_Constraints = ObjNonlinearVars.intersection(ConNonlinearVars)
     327        print >>OUTPUT, " " + str(len(ConNonlinearVars)), \
     328              str(len(ObjNonlinearVars)), str(len(Nonlinear_Vars_in_Objs_and_Constraints))\
     329              +"\t# nonlinear vars in constraints, objectives, both"
    161330        #
    162331        # LINE 6
    163332        #
    164         print >>OUTPUT, ' 0 0 0 1\t# linear network variables; functions; '   \
    165                         'arith, flags'
     333        print >>OUTPUT, " 0 0 0 1\t# linear network variables; functions; arith, flags"
    166334        #
    167335        # LINE 7
    168336        #
    169         print >>OUTPUT, " " + str(nbv), niv, 0, 0,\
    170                 "0\t# discrete variables: binary, integer, nonlinear (b,c,o)"
     337        BIset = set(Blist).union(set(Ilist)) # maybe Blist and Ilist should be sets
     338        Discrete_Nonlinear_Vars_Both = Nonlinear_Vars_in_Objs_and_Constraints.intersection(BIset)
     339        Discrete_Nonlinear_Vars_Con_Only = ConNonlinearVars.intersection(BIset).difference(Discrete_Nonlinear_Vars_Both)
     340        Discrete_Nonlinear_Vars_Obj_Only = ObjNonlinearVars.intersection(BIset).difference(Discrete_Nonlinear_Vars_Both)
     341        n_int_nonlinear_b = len(Discrete_Nonlinear_Vars_Both)
     342        n_int_nonlinear_c = len(Discrete_Nonlinear_Vars_Con_Only)
     343        n_int_nonlinear_o = len(Discrete_Nonlinear_Vars_Obj_Only)
     344        print >>OUTPUT, " " + str(len(Blist)), str(len(Ilist)), n_int_nonlinear_b, n_int_nonlinear_c,\
     345                n_int_nonlinear_o, "\t# discrete variables: binary, integer, nonlinear (b,c,o)"
    171346        #
    172347        # LINE 8
    173348        #
    174         nsno = len(Var)
    175         ngc = ngo = 0
    176         # Compute # of nonzeros in gradient
    177         for key in Obj:
    178             if Obj[key].trivial:
    179                 continue
    180             for ondx in Obj[key]:
    181                 ngo += len(Obj[key][ondx].repn.get(1,{}))
    182         # Compute # of nonzeros in Jacobian
    183         cu = {}
    184         for i in xrange(nsno):
    185             cu[i] = 0
    186         for key in Con:
    187             C = Con[key]
    188             if C.trivial:
    189                 continue
    190             for cndx in C:
    191                 for d in C[cndx].repn.get(1,{}):
    192                     ngc += 1
    193                     cu[var_id[d.keys()[0]]] += 1
    194         print >>OUTPUT, " %s %s\t# nonzeros in Jacobian, gradients"           \
    195                         % ( str(ngc), str(ngo) )
     349        # objective info computed above
     350        print >>OUTPUT, " " + str(nnz_grad_constraints), str(len(ObjVars)) + \
     351            "\t# nonzeros in Jacobian, obj. gradient"
    196352        #
    197353        # LINE 9
     
    210366                continue
    211367            for ndx in Con[key]:
    212                 if is_constant(Con[key][ndx].repn):
     368                ampl_repn = Con[key][ndx].ampl_repn
     369                if is_constant(ampl_repn):
    213370                    continue
    214 
    215                 print >>OUTPUT, 'C%s\t#%s[%s]' % ( str(nc), str(key), str(ndx))
     371                print >>OUTPUT, "C" + str(nc) + "\t#" + str(key) + "[" + str(ndx) + "]"
    216372                nc += 1
    217                 print >>OUTPUT, "n0"
     373               
     374                if ampl_repn.is_linear():
     375                    print >>OUTPUT, "n0"
     376                else:
     377                    assert(not ampl_repn._nonlinear_expr is None)
     378                    self._print_nonlinear_terms_NL(OUTPUT, ampl_repn._nonlinear_expr)
     379                                   
    218380        #
    219381        # "O" lines
     
    225387                k = 1
    226388            for ndx in Obj[key]:
    227                 if is_constant(Obj[key][ndx].repn):
    228                     continue
    229 
    230                 print >>OUTPUT, '0%s %s\t#%s[%s]' \
    231                     % ( str(no), str(k), str(key), str(ndx) )
    232 
     389                ampl_repn = Obj[key][ndx].ampl_repn
     390                if ampl_repn.is_constant():
     391                    pass
     392#                    continue
     393                print >>OUTPUT, "O" + str(no) + " "+str(k)+"\t#" + str(key) + "[" + str(ndx) + "]"
    233394                no += 1
    234                 if 0 in Obj[key][ndx].repn:
    235                     print >>OUTPUT, "n" + str(Obj[key][ndx].repn[0][None])
     395
     396                if ampl_repn.is_linear():
     397                    print >>OUTPUT, "n" + str(ampl_repn._constant)
    236398                else:
    237                     print >>OUTPUT, "n0"
     399                    assert(not ampl_repn._nonlinear_expr is None)
     400                    if ampl_repn._constant != 0.0:
     401                        print >>OUTPUT, "o0", " #+"
     402                        print >>OUTPUT, "n" + str(ampl_repn._constant)
     403                   
     404                    self._print_nonlinear_terms_NL(OUTPUT, ampl_repn._nonlinear_expr)
     405
     406        #
     407        # "x" lines
     408        #
     409        # variable initialization
     410        n_var_inits = 0
     411        for (var_name, var) in Vars.items():
     412            if var.initial != None:
     413                n_var_inits += 1
     414
     415        print >>OUTPUT, "x" + str(n_var_inits)
     416        for var_name in Vlist:
     417            var = Vars[var_name]
     418            if var.initial != None:
     419                print >>OUTPUT, str(var.ampl_var_id) + " " + str(var.initial) + " # " + var_name + " initial"
     420       
    238421        #
    239422        # "r" lines
     
    242425        nc=0
    243426        for key in Con:
    244             if Con[key].trivial:
     427            C = Con[key]
     428
     429            if C.trivial:
    245430                continue
    246             C = Con[key]
     431
    247432            for ndx in C:
    248                 if is_constant(C[ndx].repn):
     433                ampl_repn = C[ndx].ampl_repn
     434                if ampl_repn.is_constant():
    249435                    continue
    250                 if 0 in C[ndx].repn:
    251                     offset=C[ndx].repn[0][None]
    252                 else:
    253                     offset=0
     436
     437                offset = ampl_repn._constant
     438
    254439                if C[ndx]._equality:
    255440                    print >>OUTPUT, "4", self._get_bound(C[ndx].lower)-offset,
     
    259444                            print >>OUTPUT,"3",
    260445                        else:
    261                             print >>OUTPUT,"1", \
    262                             self._get_bound(C[ndx].upper)-offset,
     446                            print >>OUTPUT,"1", self._get_bound(C[ndx].upper)-offset,
    263447                    else:
    264448                        if C[ndx].upper is None:
    265                             print >>OUTPUT,"2", \
    266                             self._get_bound(C[ndx].lower)-offset,
     449                            print >>OUTPUT,"2", self._get_bound(C[ndx].lower)-offset,
    267450                        else:
    268451                            print >>OUTPUT,"0",\
    269452                                    self._get_bound(C[ndx].lower)-offset,\
    270453                                    self._get_bound(C[ndx].upper)-offset,
    271                 print >>OUTPUT, " # c"+ str(nc)+"  "+C[ndx].label
    272                 symbol_map["c" + str(nc)] = C[ndx].label
     454                print >>OUTPUT, " # c"+ str(nc)+"  "+C[ndx].label #str(key) + "[" + str(ndx) + "]"
     455                symbol_map["c" + str(nc)] = C[ndx].label #str(key) + "[" + str(ndx) + "]"
    273456                nc += 1
    274457        #
     
    277460        print >>OUTPUT, "b"
    278461        for ndx in Vlist:
    279             vv = Var[ndx]
    280             #vi = Vlist[i]
     462            vv = Vars[ndx]
    281463            var = vv.var
    282             #ndx = vi[1]
    283464            if isinstance(var.domain, BooleanSet):
    284465                print >>OUTPUT, "0 0 1",
     
    302483            else:
    303484                print >>OUTPUT, "3",
    304             print >>OUTPUT, " # v"+str(var_id[ndx])+"  "+vv.label
    305             symbol_map["v"+str(var_id[ndx])] = vv.label
     485            print >>OUTPUT, " # v"+str(vv.ampl_var_id)+"  "+vv.label
     486            symbol_map["v"+str(vv.ampl_var_id)] = vv.label
    306487        #
    307488        # "k" lines
    308489        #
    309         n1 = len(Var) - 1
     490
     491        ktot = 0
     492        cu = {}
     493        for i in xrange(len(Vars)):
     494            cu[i] = 0
     495
     496        for key in Con:
     497            C = Con[key]
     498            if C.trivial:
     499                continue
     500            for cndx in C:
     501                con_vars = set(C[cndx].ampl_repn._linear_terms_var.keys())
     502                con_vars.update(set(C[cndx].ampl_repn._nonlinear_vars.keys()))
     503                for d in con_vars:
     504                    cu[Vars[d].ampl_var_id] += 1
     505       
     506        n1 = len(Vars) - 1
    310507        print >>OUTPUT, "k" + str(n1)
    311508        ktot = 0
     
    313510                ktot += cu[i]
    314511                print >>OUTPUT, ktot
     512
    315513        #
    316514        # "J" lines
     
    322520            con = Con[key]
    323521            for ndx in con:
    324                 if is_constant(con[ndx].repn):
    325                     continue
    326                 linear = con[ndx].repn.get(1,{})
    327                 print >>OUTPUT, "J" + str(nc), len(linear)
     522                ampl_repn = con[ndx].ampl_repn
     523                # we should probably have this set already created in ampl_repn
     524                con_vars = set(ampl_repn._linear_terms_var.keys())
     525                con_vars.update(set(ampl_repn._nonlinear_vars.keys()))
     526                print >>OUTPUT, "J" + str(nc), len(con_vars), " # ", con[ndx].label
    328527                nc += 1
    329 
    330                 for d in linear:
    331                     for id in d:
    332                         print >>OUTPUT, var_id[id], linear[d]
     528                grad_entries = {}
     529                for con_var in con_vars:
     530                    if con_var in ampl_repn._linear_terms_var:
     531                        grad_entries[Vars[con_var].ampl_var_id] = ampl_repn._linear_terms_coef[con_var]
     532                    else:
     533                        assert(con_var in ampl_repn._nonlinear_vars)
     534                        grad_entries[Vars[con_var].ampl_var_id] = 0
     535
     536                sorted_keys = grad_entries.keys()
     537                sorted_keys.sort()
     538                for var_id in sorted_keys:
     539                    print >>OUTPUT, var_id, " ", grad_entries[var_id]
    333540        #
    334541        # "G" lines
    335542        #
    336543        no = 0
     544        grad_entries = {}
    337545        for key in Obj:
    338546            if Obj[key].trivial:
     
    340548            obj = Obj[key]
    341549            for ndx in obj:
    342                 linear = obj[ndx].repn.get(1,{})
    343                 print >>OUTPUT, "G" + str(no), len(linear)
     550                ampl_repn = obj[ndx].ampl_repn
     551                obj_vars = set(ampl_repn._linear_terms_var.keys())
     552                obj_vars.update(set(ampl_repn._nonlinear_vars.keys()))
     553                if len(obj_vars) == 0:
     554                    pass
     555                else:
     556                    print >>OUTPUT, "G" + str(no), len(obj_vars)
    344557                no += 1
    345                 for d in linear:
    346                     for id in d:
    347                         print >>OUTPUT, var_id[id], linear[d]
     558
     559                for obj_var in obj_vars:
     560                    if obj_var in ampl_repn._linear_terms_var:
     561                        grad_entries[Vars[obj_var].ampl_var_id] = ampl_repn._linear_terms_coef[obj_var]
     562                    else:
     563                        assert(obj_var in ampl_repn._nonlinear_vars)
     564                        grad_entries[Vars[obj_var].ampl_var_id] = 0
     565
     566        sorted_keys = grad_entries.keys()
     567        sorted_keys.sort()
     568        for var_id in sorted_keys:
     569            print >>OUTPUT, var_id, " ", grad_entries[var_id]
    348570
    349571        return symbol_map
     572   
     573class ampl_representation:
     574    def __init__(self):
     575        self._constant = 0
     576        self._linear_terms_var = {}
     577        self._linear_terms_coef = {}
     578        self._nonlinear_expr = None
     579        self._nonlinear_vars = {}
     580
     581    def clone(self):
     582        clone = ampl_representation()
     583        clone._constant = self.constant
     584        clone._linear_terms_var = self._linear_terms_var
     585        clone._linear_terms_coef = self._linear_terms_coef
     586        clone._nonlinear_expr = self._nonlinear_expr
     587        clone._nonlinear_vars = self._nonlinear_vars
     588        return clone
     589       
     590    def is_constant(self):
     591        if len(self._linear_terms_var) == 0 and self._nonlinear_expr is None:
     592#            assert(len(self._linear_terms_coef) == 0)
     593#            assert(len(self._nonlinear_vars) == 0)
     594            return True
     595        return False
     596   
     597    def is_linear(self):
     598        if self._nonlinear_expr is None:
     599            assert(len(self._nonlinear_vars) == 0)
     600            return True
     601
     602    def is_nonlinear(self):
     603        if not self._nonlinear_expr is None:
     604            return True
     605        return False
     606
     607    def needs_sum(self):
     608        num_components = 0
     609        if self._constant != 0:
     610            num_components += 1
     611       
     612        num_components += len(self._linear_terms_var)
     613       
     614        if not self._nonlinear_expr is None:
     615            num_components += 1
     616           
     617        if num_components > 1:
     618            return True
     619        return False
     620
     621    def create_expr(self):
     622        if self.needs_sum() == True:
     623            ret_expr = expr._SumExpression()
     624            if self.is_constant():
     625                ret_expr._const = self._constant
     626           
     627            for (var_name, var) in self._linear_terms_var.items():
     628                ret_expr.add(self._linear_terms_coef[var_name], var)
     629
     630            if self.is_nonlinear():
     631                ret_expr.add(1.0, self._nonlinear_expr)
     632
     633            return ret_expr
     634        else:           
     635            # we know there is only one component present in ampl_repn
     636            if self.is_constant():
     637                return as_numeric(self._constant)
     638            elif self.is_nonlinear():
     639                return self._nonlinear_expr
     640            else:
     641                assert( len(self._linear_terms_var) == 1 )
     642                return self._linear_terms_var[self._linear_terms_var.keys()[0]]
     643           
     644           
     645       
     646def generate_ampl_repn(exp):
     647    ampl_repn = ampl_representation()
     648#    exp.pprint()
     649    #
     650    # recurse to the bottom of the identity expression
     651    #
     652    while isinstance(exp, expr._IdentityExpression):
     653        exp = exp._args[0]
     654
     655    #
     656    # Expression
     657    #
     658    if isinstance(exp,expr.Expression):
     659        #
     660        # Sum
     661        #
     662        if isinstance(exp,expr._SumExpression):
     663            ampl_repn._constant = float(exp._const)
     664            ampl_repn._nonlinear_expr = None
     665            for i in range(len(exp._args)):
     666                child_repn = generate_ampl_repn(exp._args[i])
     667                # adjust the constant
     668                ampl_repn._constant += exp._coef[i]*child_repn._constant
     669
     670                # adjust the linear terms
     671                for (var_name,var) in child_repn._linear_terms_var.items():
     672                    if var_name in ampl_repn._linear_terms_var.keys():
     673                        ampl_repn._linear_terms_coef[var_name] += exp._coef[i]*child_repn._linear_terms_coef[var_name]
     674                    else:
     675                        ampl_repn._linear_terms_var[var_name] = var
     676                        ampl_repn._linear_terms_coef[var_name] = exp._coef[i]*child_repn._linear_terms_coef[var_name]
     677           
     678                # adjust the nonlinear terms
     679                if not child_repn._nonlinear_expr is None:
     680                    if ampl_repn._nonlinear_expr is None:
     681#                        ampl_repn._nonlinear_expr = expr._SumExpression([child_repn._nonlinear_expr],[exp._coef[i]])
     682                        ampl_repn._nonlinear_expr = [(exp._coef[i], child_repn._nonlinear_expr)]
     683                    else:
     684#                        assert(isinstance(ampl_repn._nonlinear_expr, expr._SumExpression))
     685                        assert(len(ampl_repn._nonlinear_expr) > 0)
     686                        ampl_repn._nonlinear_expr.append((exp._coef[i], child_repn._nonlinear_expr))
     687
     688                # adjust the nonlinear vars
     689                ampl_repn._nonlinear_vars.update(child_repn._nonlinear_vars)
     690               
     691            return ampl_repn
     692   
     693        #
     694        # Product
     695        #
     696        elif isinstance(exp,expr._ProductExpression):
     697            #
     698            # Iterate through the denominator.  If they aren't all constants, then
     699            # simply return this expresion.
     700            #
     701            denom=1.0
     702            for e in exp._denominator:
     703                if e.fixed_value():
     704                    #print 'Z',type(e),e.fixed
     705                    denom *= e.value
     706                elif e.is_constant():
     707                    denom *= e()
     708                else:
     709                    ampl_repn._nonlinear_expr = exp
     710                    break
     711                if denom == 0.0:
     712                    print "Divide-by-zero error - offending sub-expression:"
     713                    e.pprint()
     714                    raise ZeroDivisionError
     715
     716            if not ampl_repn._nonlinear_expr is None:
     717                # we have a nonlinear expression ... build up all the vars
     718                for e in exp._denominator:
     719                    arg_repn = generate_ampl_repn(e)
     720                    ampl_repn._nonlinear_vars.update(arg_repn._linear_terms_var)
     721                    ampl_repn._nonlinear_vars.update(arg_repn._nonlinear_vars)
     722
     723                for e in exp._numerator:
     724                    arg_repn = generate_ampl_repn(e)
     725                    ampl_repn._nonlinear_vars.update(arg_repn._linear_terms_var)
     726                    ampl_repn._nonlinear_vars.update(arg_repn._nonlinear_vars)
     727                return ampl_repn
     728
     729            #
     730            # OK, the denominator is a constant.
     731            #
     732            # build up the ampl_repns for the numerator
     733            n_linear_args = 0
     734            n_nonlinear_args = 0
     735            arg_repns = list()
     736            for i in range(len(exp._numerator)):
     737                e = exp._numerator[i]
     738                e_repn = generate_ampl_repn(e)
     739                arg_repns.append(e_repn)
     740                if e_repn.is_nonlinear():
     741                    n_nonlinear_args += 1
     742                elif not e_repn.is_constant():
     743                    n_linear_args += 1
     744
     745            is_nonlinear = False
     746            if n_linear_args > 1 or n_nonlinear_args > 0:
     747                is_nonlinear = True
     748               
     749               
     750            if is_nonlinear == True:
     751                # do like AMPL and simply return the expression
     752                # without extracting the potentially linear part
     753                ampl_repn = ampl_representation()
     754                ampl_repn._nonlinear_expr = exp
     755                for repn in arg_repns:
     756                    ampl_repn._nonlinear_vars.update(repn._linear_terms_var)
     757                    ampl_repn._nonlinear_vars.update(repn._nonlinear_vars)
     758                return ampl_repn
     759
     760            else: # is linear or constant
     761                ampl_repn = current_repn = arg_repns[0]
     762                for i in range(1,len(arg_repns)):
     763                    e_repn = arg_repns[i]
     764                    ampl_repn = ampl_representation()
     765
     766                    # const_c * const_e
     767                    ampl_repn._constant = current_repn._constant * e_repn._constant
     768   
     769                    # const_e * L_c
     770                    if e_repn._constant != 0.0:
     771                        for (var_name, var) in current_repn._linear_terms_var.items():
     772                            ampl_repn._linear_terms_coef[var_name] = current_repn._linear_terms_coef[var_name] * e_repn._constant
     773                            ampl_repn._linear_terms_var[var_name] = var
     774   
     775                    # const_c * L_e
     776                    if current_repn._constant != 0.0:
     777                        for (e_var_name,e_var) in e_repn._linear_terms_var.items():
     778                            if e_var_name in ampl_repn._linear_terms_var:
     779                                ampl_repn._linear_terms_coef[e_var_name] += current_repn._constant * e_repn._linear_terms_coef[e_var_name]
     780                            else:
     781                                ampl_repn._linear_terms_coef[e_var_name] = current_repn._constant * e_repn._linear_terms_coef[e_var_name]
     782                                ampl_repn._linear_terms_var[e_var_name] = e_var
     783               
     784                current_repn = ampl_repn
     785
     786            # now deal with the product expression's coefficient that needs
     787            # to be applied to all parts of the ampl_repn
     788            ampl_repn._constant *= exp.coef/denom
     789            for var_name in ampl_repn._linear_terms_coef.keys():
     790                ampl_repn._linear_terms_coef[var_name] *= exp.coef/denom
     791            assert(ampl_repn._nonlinear_expr is None)
     792            assert(len(ampl_repn._nonlinear_vars) == 0)
     793
     794            return ampl_repn
     795
     796        elif isinstance(exp,expr._PowExpression):
     797            assert(len(exp._args) == 2)
     798            base_repn = generate_ampl_repn(exp._args[0])
     799            exponent_repn = generate_ampl_repn(exp._args[1])
     800
     801            if base_repn.is_constant() and exponent_repn.is_constant():
     802                ampl_repn._constant = base_repn._constant**exponent_repn._constant
     803            elif exponent_repn.is_constant() and exponent_repn._constant == 1.0:
     804                # not sure if this is necessary (may be able to just return base_repn
     805                ampl_repn = base_repn.clone()
     806            elif exponent_repn.is_constant() and exponent_repn._constant == 0.0:
     807                ampl_repn._constant = 1.0
     808            else:
     809                # this might be overkill, creating another sum expression
     810#                base_expr = base_repn.create_expr()
     811#                exponent_expr = exponent_repn.create_expr()
     812#                ampl_repn._nonlinear_expr = expr._PowExpression(base_expr, exponent_expr)
     813
     814                # instead, let's just return the expression we are given and only
     815                # use the ampl_repn for the vars
     816                ampl_repn._nonlinear_expr = exp
     817                ampl_repn._nonlinear_vars = base_repn._nonlinear_vars
     818                ampl_repn._nonlinear_vars.update(exponent_repn._nonlinear_vars)
     819                ampl_repn._nonlinear_vars.update(base_repn._linear_terms_var)
     820                ampl_repn._nonlinear_vars.update(exponent_repn._linear_terms_var)
     821
     822            return ampl_repn
     823             
     824        elif isinstance(exp,intrinsic_functions._IntrinsicFuncExpression):
     825            assert(exp._nargs == 1)
     826
     827            # this is inefficient since it is using much more than what we need
     828            child_repn = generate_ampl_repn(exp._args[0])
     829#            assert(child_repn.is_constant() == False)
     830           
     831            ampl_repn._nonlinear_expr = exp
     832            ampl_repn._nonlinear_vars = child_repn._nonlinear_vars
     833            ampl_repn._nonlinear_vars.update(child_repn._linear_terms_var)
     834            return ampl_repn
     835        #
     836        # ERROR
     837        #
     838        else:
     839            raise ValueError, "Unsupported expression type: "+str(exp)
     840    #
     841    # Constant
     842    #
     843    elif exp.fixed_value():
     844        ampl_repn._constant = float(exp.value)
     845        return ampl_repn
     846    #
     847    # Variable
     848    #
     849    elif type(exp) is _VarValue or exp.type() is Var:
     850        ampl_repn._linear_terms_coef = {exp.label: 1.0}
     851        ampl_repn._linear_terms_var = {exp.label: exp}
     852        return ampl_repn
     853    #
     854    # ERROR
     855    #
     856    else:
     857       raise ValueError, "Unexpected expression type: "+str(exp)
     858
  • coopr.pyomo/trunk/coopr/pyomo/scripting/pyomo.py

    r2761 r2882  
    229229        type='string',
    230230        default=None)
     231    parser.add_option("--skip-canonical-repn",
     232            help="Do not create the canonical representation. This is not necessary for solvers that do not require it.",
     233            action="store_true",
     234            dest="skip_canonical_repn",
     235            default=False)
    231236    #
    232237    # These options are depricated until we have a specific use-case for them
  • coopr.pyomo/trunk/coopr/pyomo/scripting/util.py

    r2773 r2882  
    9797        TempfileManager.tempdir = options.tempdir
    9898    #
     99    # Disable canonical repn for ASL solvers
     100    #
     101    if hasattr(options, 'solve') and options.solver == "asl":
     102        options.skip_canonical_repn = True
     103
     104    #
    99105    # Configure exception management
    100106    #
     
    246252    for ep in ExtensionPoint(IPyomoScriptPrintModel):
    247253        ep.apply( options=options, model=model )
     254       
     255    # ToDo: CDL This may not be the right place for this
     256    # Likely we need to change the framework so that canonical repn
     257    # is not assumed to be required by all solvers?
     258    # For now...
     259    if hasattr(options, 'skip_canonical_rep') and options.skip_canonical_repn == True:
     260        if 'simple_preprocessor' in model._preprocessors:
     261            model.preprocessor_ep.service('simple_preprocessor').deactivate_action('compute_canonical_repn')
     262   
    248263    #
    249264    # Create Problem Instance
Note: See TracChangeset for help on using the changeset viewer.