source: coopr.pyomo/trunk/coopr/pyomo/io/ampl.py @ 2217

Last change on this file since 2217 was 2217, checked in by wehart, 10 years ago

Update to NL file writer, after recent changes in the canonical representation.

File size: 10.9 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the Coopr README.txt file.
9#  _________________________________________________________________________
10
11#
12# AMPL Problem Writer Plugin
13#
14
15from coopr.opt import ProblemFormat, WriterFactory
16from coopr.pyomo.base.numtypes import minimize, maximize
17from coopr.pyomo.base import *
18from coopr.opt.base import *
19from coopr.pyomo.base.param import _ParamValue
20from coopr.pyomo.base.var import _VarBase
21
22class ProblemWriter_nl(AbstractProblemWriter):
23
24    def __init__(self):
25        AbstractProblemWriter.__init__(self,ProblemFormat.nl)
26
27    def __call__(self, model, filename):
28        if filename is None:
29           filename = model.name + ".nl"
30        OUTPUT=open(filename,"w")
31        symbol_map = self._print_model_NL(model,OUTPUT)
32        OUTPUT.close()
33        return filename, symbol_map
34
35    def _get_bound(self, exp):
36        if isinstance(exp,expr._IdentityExpression):
37           return self._get_bound(exp._args[0])
38        elif isinstance(exp,NumericConstant):
39           return exp.value
40        elif isinstance(exp,_ParamValue):
41           return exp.value
42        else:
43           raise ValueError, "ERROR: nonconstant bound: " + str(exp)
44           return None
45
46    def _print_model_NL(self, model, OUTPUT, verbose=False):
47
48        # maps NL variables to the "real" variable names in the problem.
49        # it's really NL variable ordering, as there are no variable names
50        # in the NL format. however, we by convention make them go from
51        # x0 upward.
52        symbol_map = {}
53
54        #
55        # Collect statistics
56        #
57        nc = no = neqn = nrange = 0
58        Obj = model.active_components(Objective)
59        Con = model.active_components(Constraint)
60        Var = {}   # Dictionary of all variables used in all expressions in this model
61        #
62        # Count number of objectives
63        #
64        for obj in Obj:
65            if Obj[obj].trivial:
66                continue
67            for i in Obj[obj]:
68                if not is_constant(Obj[obj][i].repn):
69                    # Merge variables from this expression
70                    Var.update( Obj[obj][i].repn[-1] )
71                    no += 1
72        #
73        # Count number of constraints
74        #
75        for con in Con:
76            C = Con[con]
77            if C.trivial:
78                continue
79            for i in C:
80                if is_constant(C[i].repn):
81                    continue
82                # Merge variables from this expression
83                Var.update( C[i].repn[-1] )
84                nc += 1
85                if C[i].lower is not None:
86                    L = self._get_bound(C[i].lower) if C[i].lower is not None else -inf
87                    U = self._get_bound(C[i].upper) if C[i].upper is not None else inf
88                    if (L == U):
89                        neqn += 1
90                elif L > U:
91                    raise ValueError, "Constraint " + str(con) +\
92                                             "[" + str(i) + "]: lower bound = " + str(L) +\
93                                             " > upper bound = " + str(U)
94                elif C[i].lower is not None and C[i].upper is not None:
95                    nrange += 1
96        #
97        niv = nbv = 0
98        Ilist = []
99        Blist = []
100        Vlist = []
101        keys = Var.keys()
102        keys.sort()
103        for ndx in keys:
104            var = Var[ndx].var
105            if isinstance(var.domain, IntegerSet):
106                Ilist.append(ndx)
107            elif isinstance(var.domain, BooleanSet):
108                L = Var[ndx].lb
109                U = Var[ndx].ub
110                if L is not None and value(L) == 0 \
111                    and U is not None and value(U) == 1:
112                    Blist.append(ndx)
113                else:
114                    Ilist.append(ndx)
115            else:
116                Vlist.append(ndx)
117        for ndx in Blist:
118            nbv += 1
119            Vlist.append(ndx)
120        for ndxx in Ilist:
121            niv += 1
122            Vlist.append(ndx)
123        #
124        vctr=0
125        var_id = {}
126        for ndx in Vlist:
127            var_id[ndx] = vctr
128            vctr += 1
129        #
130        # Print Header
131        #
132        # LINE 1
133        #
134        print >>OUTPUT,"g3 1 1 0\t# problem",model.name
135        #
136        # LINE 2
137        #
138        print >>OUTPUT, " ",len(Var), nc, no, nrange, str(neqn) +\
139                "\t# vars, constraints, objectives, ranges, eqns"
140        #
141        # LINE 3
142        #
143        print >>OUTPUT, " 0 0\t# nonlinear constraints, objectives"
144        #
145        # LINE 4
146        #
147        print >>OUTPUT, " 0 0\t# network constraints: nonlinear, linear"
148        #
149        # LINE 5
150        #
151        print >>OUTPUT, " 0 0 0\t# nonlinear vars in constraints, objectives, both"
152        #
153        # LINE 6
154        #
155        print >>OUTPUT, " 0 0 0 1\t# linear network variables; functions; arith, flags"
156        #
157        # LINE 7
158        #
159        print >>OUTPUT, " " + str(nbv), niv, 0, 0,\
160                "0\t# discrete variables: binary, integer, nonlinear (b,c,o)"
161        #
162        # LINE 8
163        #
164        nsno = len(Var)
165        ngc = ngo = 0
166        # Compute # of nonzeros in gradient
167        for key in Obj:
168            if Obj[key].trivial:
169                continue
170            for ondx in Obj[key]:
171                ngo += len(Obj[key][ondx].repn.get(1,{}))
172        # Compute # of nonzeros in Jacobian
173        cu = {}
174        for i in xrange(nsno):
175            cu[i] = 0
176        for key in Con:
177            C = Con[key]
178            if C.trivial:
179                continue
180            for cndx in C:
181                for d in C[cndx].repn.get(1,{}):
182                    ngc += 1
183                    cu[var_id[d.keys()[0]]] += 1
184        print >>OUTPUT, " " + str(ngc), str(ngo) + "\t# nonzeros in Jacobian, gradients"
185        #
186        # LINE 9
187        #
188        print >>OUTPUT, " 0 0\t# max name lengths: constraints, variables"
189        #
190        # LINE 10
191        #
192        print >>OUTPUT, " 0 0 0 0 0\t# common exprs: b,c,o,c1,o1"
193        #
194        # "C" lines
195        #
196        nc = 0
197        for key in Con:
198            if Con[key].trivial:
199                continue
200            for ndx in Con[key]:
201                if is_constant(Con[key][ndx].repn):
202                    continue
203                print >>OUTPUT, "C" + str(nc) + "\t#" + str(key) + "[" + str(ndx) + "]"
204                nc += 1
205                print >>OUTPUT, "n0"
206        #
207        # "O" lines
208        #
209        no = 0
210        for key in Obj:
211            k = 0
212            if Obj[key].sense == maximize:
213                k = 1
214            for ndx in Obj[key]:
215                if is_constant(Obj[key][ndx].repn):
216                    continue
217                print >>OUTPUT, "O" + str(no) + " "+str(k)+"\t#" + str(key) + "[" + str(ndx) + "]"
218                no += 1
219                if 0 in Obj[key][ndx].repn:
220                    print >>OUTPUT, "n" + str(Obj[key][ndx].repn[0][None])
221                else:
222                    print >>OUTPUT, "n0"
223        #
224        # "r" lines
225        #
226        print >>OUTPUT, "r"
227        nc=0
228        for key in Con:
229            if Con[key].trivial:
230                continue
231            C = Con[key]
232            for ndx in C:
233                if is_constant(C[ndx].repn):
234                    continue
235                if 0 in C[ndx].repn:
236                    offset=C[ndx].repn[0][None]
237                else:
238                    offset=0
239                if C[ndx]._equality:
240                    print >>OUTPUT, "4", self._get_bound(C[ndx].lower)-offset,
241                else:
242                    if C[ndx].lower is None:
243                        if C[ndx].upper is None:
244                            print >>OUTPUT,"3",
245                        else:
246                            print >>OUTPUT,"1", self._get_bound(C[ndx].upper)-offset,
247                    else:
248                        if C[ndx].upper is None:
249                            print >>OUTPUT,"2", self._get_bound(C[ndx].lower)-offset,
250                        else:
251                            print >>OUTPUT,"0",\
252                                    self._get_bound(C[ndx].lower)-offset,\
253                                    self._get_bound(C[ndx].upper)-offset,
254                print >>OUTPUT, " # c"+ str(nc)+"  "+str(key) + "[" + str(ndx) + "]"
255                symbol_map["c" + str(nc)] = str(key) + "[" + str(ndx) + "]"
256                nc += 1
257        #
258        # "b" lines
259        #
260        print >>OUTPUT, "b"
261        nv=0
262        for ndx in Vlist:
263            vv = Var[ndx]
264            #vi = Vlist[i]
265            var = vv.var
266            #ndx = vi[1]
267            if isinstance(var.domain, BooleanSet):
268                print >>OUTPUT, "0 0 1",
269                print >>OUTPUT, " # v"+str(nv)+"  "+vv.label
270                nv += 1
271                continue
272            L = vv.lb
273            U = vv.ub
274            if L is not None:
275                Lv = str(value(L))
276                if U is not None:
277                    Uv = str(value(U))
278                    if Lv == Uv:
279                        print >>OUTPUT, "4" , Lv,
280                    else:
281                        print >>OUTPUT, "0", Lv, Uv,
282                else:
283                    print >>OUTPUT, "2", Lv,
284            elif U is not None:
285                print >>OUTPUT, "1", str(value(U)),
286            else:
287                print >>OUTPUT, "3",
288            print >>OUTPUT, " # v"+str(nv)+"  "+vv.label
289            symbol_map["v"+str(nv)] = vv.label
290
291            nv += 1
292        #
293        # "k" lines
294        #
295        n1 = len(Var) - 1
296        print >>OUTPUT, "k" + str(n1)
297        ktot = 0
298        for i in xrange(n1):
299                ktot += cu[i]
300                print >>OUTPUT, ktot
301        #
302        # "J" lines
303        #
304        nc = 0
305        for key in Con:
306            if Con[key].trivial:
307                continue
308            con = Con[key]
309            for ndx in con:
310                if is_constant(con[ndx].repn):
311                    continue
312                linear = con[ndx].repn.get(1,{})
313                print >>OUTPUT, "J" + str(nc), len(linear)
314                nc += 1
315
316                for d in linear:
317                    for id in d:
318                        print >>OUTPUT, var_id[id], linear[d]
319        #
320        # "G" lines
321        #
322        no = 0
323        for key in Obj:
324            if Obj[key].trivial:
325                continue
326            obj = Obj[key]
327            for ndx in obj:
328                linear = obj[ndx].repn.get(1,{})
329                print >>OUTPUT, "G" + str(no), len(linear)
330                no += 1
331                for d in linear:
332                    for id in d:
333                        print >>OUTPUT, var_id[id], linear[d]
334
335        return symbol_map
336
337problem.WriterRegistration(str(ProblemFormat.nl), ProblemWriter_nl)
338
Note: See TracBrowser for help on using the repository browser.