source: coopr.plugins/trunk/coopr/plugins/mip/glpk_direct.py @ 3664

Last change on this file since 3664 was 3664, checked in by wehart, 9 years ago

Adding conditional import of pyutilib.logging
File size: 19.8 KB
Line 
1# vim: set fileencoding=utf-8
2# _________________________________________________________________________
3#
4#  Coopr: A COmmon Optimization Python Repository
5#  Copyright (c) 2010 Sandia Corporation.
6#  This software is distributed under the BSD License.
7#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8#  the U.S. Government retains certain rights in this software.
9#  For more information, see the Coopr README.txt file.
10#  _________________________________________________________________________
11
12
13try:
14    import pyutilib.logging as logging
15except ImportError:
16    import logging
17
18from coopr.opt.base import *
19from coopr.opt.results import *
20from coopr.opt.solver import *
21
22import mockmip
23from pyutilib.misc import Options
24from pyutilib.component.core import alias
25
26Logger = logging.getLogger('coopr.plugins')
27
28try:
29   # import all the glp_* functions
30   from glpk.glpkpi import *
31   glpk_python_api_exists = True
32except ImportError:
33   glpk_python_api_exists = False
34
35
36class GLPKDirect ( OptSolver ):
37   """The GLPK LP/MIP solver (direct API plugin)
38
39The glpk_direct plugin offers an API interface to the GLPK.  It requires the
40Python-GLPK API interface provided via the SWIG interface available through
41most Linux distributions repositories.  For more information, see João Pedro
42Pedroso's (the author) page at http://www.dcc.fc.up.pt/~jpp/
43
44Because of the direct connection with the GLPK API, no temporary files need be
45written or read.  That ostensibly makes this a faster plugin than the
46file-based glpk plugin.  However, you will likely not notice any speed up
47unless you are using the GLPK solver with PySP problems (due to the rapid
48re-solves).
49
50One downside to the lack of temporary files, is that there is no LP file to
51inspect for clues while debugging a model.  For that, use the write_lp solver
52option:
53
54$ pyomo model.{py,dat} \
55  --solver=glpk_direct \
56  --solver-options  write_lp=/path/to/some/file.lp
57
58One can also specify the particular GLPK algorithm to use with the 'algorithm'
59solver option.  There are 4 algorithms:
60
61  simplex  - the default algorithm for non-MIP problems (primal simplex)
62  intopt   - the default algorithm for MIP problems
63  exact    - tentative implementation of two-phase primal simplex with exact
64             arithmetic internally.
65  interior - only aplicable to LP problems
66
67$ pyomo model.{py,dat} \
68  --solver glpk_direct \
69  --solver-options  algorithm=exact
70
71For more information on available algorithms, see the GLPK documentation.
72   """
73
74   alias('glpk_direct', doc='Direct Python interface to the GLPK LP/MIP solver.')
75
76
77   def __init__(self, **kwds):
78      #
79      # Call base class constructor
80      #
81      kwds['type'] = 'glpk_direct'
82      OptSolver.__init__(self, **kwds)
83
84      # NOTE: eventually both of the following attributes should be migrated
85      # to a common base class.  Is the current solve warm-started?  A
86      # transient data member to communicate state information across the
87      # _presolve, _apply_solver, and _postsolve methods.
88      self.warm_start_solve = False
89      self._timelimit = None
90
91      # Note: Undefined capabilities default to 'None'
92      self._capabilities = Options()
93      self._capabilities.linear = True
94      self._capabilities.integer = True
95
96
97   def _populate_glpk_instance ( self, model ):
98      from coopr.pyomo.base import Var, VarStatus, Objective, Constraint, \
99                                   IntegerSet, BooleanSet
100      from coopr.pyomo.expr import is_constant
101
102      try:
103         lp = glp_create_prob()
104      except Exception, e:
105         msg = 'Unable to create GLPK problem instance.  Have you installed' \
106         '\n       the Python bindings for GLPK?\n\n\tError message: %s'
107         raise Exception, msg % e
108
109      objective = sorted( model.active_components( Objective ).values() )[0]
110      sense = GLP_MAX
111      if objective.is_minimizing(): sense = GLP_MIN
112
113      constraint_list = model.active_components( Constraint )
114      variable_list   = model.active_components( Var )
115      num_constraints = model.statistics.number_of_constraints
116      num_variables   = model.statistics.number_of_variables
117
118      glp_set_prob_name( lp, model.name )
119
120      glp_set_obj_dir( lp, sense )
121      glp_add_rows( lp, num_constraints )
122      glp_add_cols( lp, num_variables )
123
124      # 1 extra because GLPK's arrays in this context are 1-based, not 0-based
125      coef_count = num_constraints * num_variables + 1
126      Ai = intArray( coef_count )
127      Aj = intArray( coef_count )
128      Ar = doubleArray( coef_count )
129
130      row = col = coef_count = 0
131      colvar_map = dict()
132      rowvar_map = dict()
133
134      # In matrix parlance, variables are columns
135      for name in variable_list.keys():
136         var_set = variable_list[ name ]
137         for ii in var_set.keys():
138            var = var_set[ ii ]
139            if (not var.active) or (var.status is VarStatus.unused) \
140               or (var.fixed is True):
141               continue
142
143            lb = ub = 0.0
144            if var.lb is None and var.ub is None:
145               var_type = GLP_FR
146            elif var.lb is None:
147               var_type = GLP_UB
148               ub = var.ub()
149            elif var.ub is None:
150               var_type = GLP_LO
151               lb = var.lb()
152            else:
153               var_type = GLP_DB
154               lb = var.lb()
155               ub = var.ub()
156
157            col += 1
158            colvar_map[ var.label ] = col
159
160            # the name is perhaps not necessary, but for completeness ...
161            glp_set_col_name( lp, col, var.label )
162            glp_set_col_bnds( lp, col, var_type, lb, ub )
163
164            # Be sure to impart the integer and binary nature of any variables
165            if isinstance(var.domain, IntegerSet):
166               glp_set_col_kind( lp, col, GLP_IV )
167            elif isinstance(var.domain, BooleanSet):
168               glp_set_col_kind( lp, col, GLP_BV )
169            else:
170               glp_set_col_kind( lp, col, GLP_CV )   # continuous
171
172      for name in constraint_list.keys():
173         constraint_set = constraint_list[ name ]
174         if constraint_set.trivial: continue
175
176         for ii in constraint_set.keys():
177            constraint = constraint_set[ ii ]
178            if not constraint.active: continue
179            elif constraint.lower is None and constraint.upper is None:
180               continue
181
182            offset = 0.0
183            if 0 in constraint.repn:
184               offset = constraint.repn[0][None]
185
186            lbound = ubound = -offset
187
188            if constraint._equality:
189               var_type = GLP_FX    # Fixed
190               lbound = ubound = constraint.lower() - offset
191            elif constraint.lower is None:
192               var_type = GLP_UP    # Upper bounded only
193               ubound += constraint.upper()
194            elif constraint.upper is None:
195               var_type = GLP_LO    # Lower bounded only
196               lbound += constraint.lower()
197            else:
198               var_type = GLP_DB    # Double bounded
199               lbound += constraint.lower()
200               ubound += constraint.upper()
201
202            row += 1
203            rowvar_map[ constraint.label ] = row
204
205            # just as with variables, set the name just for completeness ...
206            glp_set_row_name( lp, row, constraint.label )
207            glp_set_row_bnds( lp, row, var_type, lbound, ubound )
208
209            expression = constraint.repn
210            if 1 in expression: # first-order terms
211               keys = sorted( expression[1].keys() )
212               for var_key in keys:
213                  index = var_key.keys()[0]
214                  var = expression[-1][ index ]
215                  coef  = expression[ 1][ var_key ]
216                  col = colvar_map[ var.label ]
217
218                  coef_count += 1
219                  Ai[ coef_count ] = row
220                  Aj[ coef_count ] = col
221                  Ar[ coef_count ] = coef
222
223      # with the rows and columns named and bounded, load the coefficients
224      glp_load_matrix( lp, coef_count, Ai, Aj, Ar )
225
226      for key in objective:
227         expression = objective[ key ].repn
228         if is_constant( expression ):
229            msg = "Ignoring objective '%s[%s]' which is constant"
230            Logger.warning( msg % (str(objective), str(key)) )
231            continue
232
233         if 1 in expression: # first-order terms
234            keys = sorted( expression[1].keys() )
235            for var_key in keys:
236               index = var_key.keys()[0]
237               label = expression[-1][ index ].label
238               coef  = expression[ 1][ var_key ]
239               col = colvar_map[ label ]
240               glp_set_obj_coef( lp, col, coef )
241
242         elif -1 in expression:
243            pass
244         else:
245            msg = "Nonlinear objective to GLPK.  GLPK can only handle "       \
246                  "linear problems."
247            raise RuntimeError( msg )
248
249
250      self._glpk_instance = lp
251      self._glpk_rowvar_map = rowvar_map
252      self._glpk_colvar_map = colvar_map
253
254
255   def warm_start_capable(self):
256      # Note to dev who gets back here: GLPK has the notion of presolving, but
257      # it's "built-in" to each optimization function.  To disable it, make use
258      # of the second argument of glp_smcp type.  See the GLPK documentation
259      # PDF for further information.
260      msg = "GLPK has the ability to use warmstart solutions.  However, it "  \
261            "has not yet been implemented into the Coopr glpk_direct plugin."
262      Logger.info( msg )
263      return False
264
265
266   def warm_start(self, instance):
267      pass
268
269
270   def _presolve(self, *args, **kwargs):
271      from coopr.pyomo.base.PyomoModel import Model
272
273      self.warm_start_solve = kwargs.pop( 'warmstart', False )
274
275      model = args[0]
276      if len(args) != 1:
277         msg = "The glpk_direct plugin method '_presolve' must be supplied "  \
278               "a single problem instance - %s were supplied"
279         raise ValueError, msg % len(args)
280      elif not isinstance(model, Model):
281         msg = "The problem instance supplied to the glpk_direct plugin "     \
282               "'_presolve' method must be of type 'Model'"
283         raise ValueError, msg
284
285      self._populate_glpk_instance( model )
286      lp = self._glpk_instance
287      self.is_integer = ( glp_get_num_int( lp ) > 0 and True or False )
288
289      if 'write_lp' in self.options:
290         fname = self.options.write_lp
291         glp_write_lp( lp, None, fname )
292
293      self.algo = 'Simplex (primal)'
294      algorithm = glp_simplex
295      if self.is_integer > 0:
296         self.algo = 'Mixed Integer'
297         algorithm = glp_intopt
298
299      if 'algorithm' in self.options:
300         if 'simplex' == self.options.algorithm:
301            self.algo = 'Simplex (primal)'
302            algorithm = glp_simplex
303         elif 'exact' == self.options.algorithm:
304            self.algo = 'Simplex (two-phase primal)'
305            algorithm = glp_exact
306         elif 'interior' == self.options.algorithm:
307            self.algo = 'Interior Point'
308            algorithm = glp_interior
309         elif 'intopt' == self.options.algorithm:
310            self.alpo = 'Mixed Integer'
311            algorithm = glp_intopt
312         else:
313            msg = "Unknown solver specified\n  Unknown: %s\n  Using:   %s\n"
314            Logger.warn( msg % (self.options.algorithm, self.algo) )
315      self._algorithm = algorithm
316
317      if 'Simplex (primal)' == self.algo:
318         parm = glp_smcp()
319         glp_init_smcp( parm )
320      elif 'Simplex (two-phase primal)' == self.algo:
321         parm = glp_smcp()
322         glp_init_smcp( parm )
323      elif 'Interior Point' == self.algo:
324         parm = glp_iptcp()
325         glp_init_iptcp( parm )
326      elif 'Mixed Integer' == self.algo:
327         parm = glp_iocp()
328         glp_init_iocp( parm )
329         if self.mipgap:
330            parm.mip_gap = self.mipgap
331
332      if self._timelimit and self._timelimit > 0.0:
333         parm.tm_lim = self._timelimit
334
335      parm.msg_lev = GLP_MSG_OFF
336      parm.presolve = GLP_ON
337
338      self._solver_params = parm
339
340      # Scaffolding in place: I /believe/ GLPK can do warmstarts; just supply
341      # the basis solution and turn of the presolver
342      if self.warm_start_solve is True:
343
344         if len(args) != 1:
345            msg = "The glpk_direct _presolve method can only handle a single "\
346                  "problem instance - %s were supplied"
347            raise ValueError, msg % len(args)
348
349         parm.presolve = GLP_OFF
350         self.warm_start( model )
351
352
353   def _apply_solver(self):
354      lp = self._glpk_instance
355      parm = self._solver_params
356      algorithm = self._algorithm
357
358      # Actually solve the problem.
359      try:
360         beg = glp_time()
361         self.solve_return_code = algorithm( self._glpk_instance, parm )
362         end = glp_time()
363         self._glpk_solve_time = glp_difftime( end, beg )
364      except Exception, e:
365         msg = str(e)
366         if 'algorithm' in self.options:
367            msg = "Unexpected error using '%s' algorithm.  Is it the correct "\
368                  "correct algorithm for the problem type?"
369            msg %= self.options.algorithm
370         Logger.error( msg )
371         raise
372
373
374   def _glpk_return_code_to_message ( self ):
375      code = self.solve_return_code
376      if 0 == code:
377         return "Algorithm completed successfully.  (This does not "         \
378                "necessarily mean an optimal solution was found.)"
379      elif GLP_EBADB == code:
380         return "Unable to start the search, because the initial basis "     \
381            "specified in the problem object is invalid -- the number of "   \
382            "basic (auxiliary and structural) variables is not the same as " \
383            "the number of rows in the problem object."
384      elif GLP_ESING == code:
385         return "Unable to start the search, because the basis matrix "      \
386            "corresponding to the initial basis is singular within the "     \
387            "working precision."
388      elif GLP_ECOND == code:
389        return "Unable to start the search, because the basis matrix "       \
390           "corresponding to the initial basis is ill-conditioned, i.e. its "\
391           "condition number is too large."
392      elif GLP_EBOUND == code:
393        return "Unable to start the search, because some double-bounded "    \
394           "(auxiliary or structural) variables have incorrect bounds."
395      elif GLP_EFAIL == code:
396        return "The search was prematurely terminated due to the solver "    \
397           "failure."
398      elif GLP_EOBJLL == code:
399        return "The search was prematurely terminated, because the "         \
400           "objective function being maximized has reached its lower limit " \
401           "and continues decreasing (the dual simplex only)."
402      elif GLP_EOBJUL == code:
403        return "The search was prematurely terminated, because the "         \
404           "objective function being minimized has reached its upper limit " \
405           "and continues increasing (the dual simplex only)."
406      elif GLP_EITLIM == code:
407        return "The search was prematurely terminated, because the simplex " \
408           "iteration limit has been exceeded."
409      elif GLP_ETMLIM == code:
410        return "The search was prematurely terminated, because the time "    \
411           "limit has been exceeded."
412      elif GLP_ENOPFS == code:
413        return "The LP problem instance has no primal feasible solution "    \
414           "(only if the LP presolver is used)."
415      elif GLP_ENODFS == code:
416         return "The LP problem instance has no dual feasible solution "     \
417            "(only if the LP presolver is used)."
418      else:
419         return "Unexpected error condition.  Please consider remitting "    \
420            "this problem to the Coopr developers and/or the GLPK project "  \
421            "so they can improve their softwares."
422
423
424   def _glpk_get_solution_status ( self ):
425      getstatus = glp_get_status
426      if self.is_integer: getstatus = glp_mip_status
427
428      status = getstatus( self._glpk_instance )
429      if   GLP_OPT    == status: return SolutionStatus.optimal
430      elif GLP_FEAS   == status: return SolutionStatus.feasible
431      elif GLP_INFEAS == status: return SolutionStatus.infeasible
432      elif GLP_NOFEAS == status: return SolutionStatus.other
433      elif GLP_UNBND  == status: return SolutionStatus.other
434      elif GLP_UNDEF  == status: return SolutionStatus.other
435      raise RuntimeError, "Unknown solution status returned by GLPK solver"
436
437
438   def _glpk_get_solver_status ( self ):
439      rc = self.solve_return_code
440      if 0 == rc:            return SolverStatus.ok
441      elif GLP_EBADB  == rc: return SolverStatus.error
442      elif GLP_ESING  == rc: return SolverStatus.error
443      elif GLP_ECOND  == rc: return SolverStatus.error
444      elif GLP_EBOUND == rc: return SolverStatus.error
445      elif GLP_EFAIL  == rc: return SolverStatus.aborted
446      elif GLP_EOBJLL == rc: return SolverStatus.aborted
447      elif GLP_EOBJUL == rc: return SolverStatus.aborted
448      elif GLP_EITLIM == rc: return SolverStatus.aborted
449      elif GLP_ETMLIM == rc: return SolverStatus.aborted
450      elif GLP_ENOPFS == rc: return SolverStatus.warning
451      elif GLP_ENODFS == rc: return SolverStatus.warning
452      else: return SolverStatus.unkown
453
454
455   def _postsolve(self):
456      lp = self._glpk_instance
457      num_variables = glp_get_num_cols( lp )
458      bin_variables = glp_get_num_bin( lp )
459      int_variables = glp_get_num_int( lp )
460
461      tpeak = glp_long()
462      glp_mem_usage( None, None, None, tpeak )
463      # black magic trickery, thanks to Python's lack of pointers and SWIG's
464      # automatic API conversion
465      peak_mem = tpeak.lo
466
467      results = SolverResults()
468      soln = Solution()
469      prob = results.problem
470      solv = results.solver
471
472      solv.name = "GLPK " + glp_version()
473      solv.status = self._glpk_get_solver_status()
474      solv.return_code = self.solve_return_code
475      solv.message = self._glpk_return_code_to_message()
476      solv.algorithm = self.algo
477      solv.memory_used = "%d bytes, (%d KiB)" % (peak_mem, peak_mem/1024)
478      # solv.user_time = None
479      # solv.system_time = None
480      solv.wallclock_time = self._glpk_solve_time
481      # solv.termination_condition = None
482      # solv.termination_message = None
483
484      prob.name = glp_get_prob_name( lp )
485      prob.number_of_constraints = glp_get_num_rows( lp )
486      prob.number_of_nonzeros = glp_get_num_nz( lp )
487      prob.number_of_variables = num_variables
488      prob.number_of_binary_variables = bin_variables
489      prob.number_of_integer_variables = int_variables
490      prob.number_of_continuous_variables = num_variables - int_variables
491      prob.number_of_objectives = 1
492
493      prob.sense = ProblemSense.minimize
494      if GLP_MAX == glp_get_obj_dir( lp ):
495         prob.sense = ProblemSense.maximize
496
497      soln.status = self._glpk_get_solution_status()
498
499      if soln.status in ( SolutionStatus.optimal, SolutionStatus.feasible ):
500         get_col_prim = glp_get_col_prim
501         get_row_prim = glp_get_row_prim
502         get_obj_val  = glp_get_obj_val
503         if self.is_integer:
504            get_col_prim = glp_mip_col_val
505            get_row_prim = glp_mip_row_val
506            get_obj_val  = glp_mip_obj_val
507
508         obj_val = get_obj_val( lp )
509         if prob.sense == ProblemSense.minimize:
510            prob.lower_bound = obj_val
511         else:
512            prob.upper_bound = obj_val
513
514         soln.objective['f'].value = obj_val
515
516         colvar_map = self._glpk_colvar_map
517         rowvar_map = self._glpk_rowvar_map
518
519         for var_label in colvar_map.keys():
520            col = colvar_map[ var_label ]
521            soln.variable[ var_label ] = get_col_prim( lp, col )
522
523         for row_label in rowvar_map.keys():
524            row = rowvar_map[ row_label ]
525            soln.constraint[ row_label ] = get_row_prim( lp, row )
526
527
528      results.solution.insert(soln)
529
530      self.results = results
531
532      # All done with the GLPK object, so free up some memory.
533      glp_free( lp )
534      del self._glpk_instance, lp
535
536      # let the base class deal with returning results.
537      return OptSolver._postsolve(self)
538
539
540# TODO: add MockGLPKDirect class
541
542if not glpk_python_api_exists:
543   SolverFactory().deactivate('glpk_direct')
544   # SolverFactory().deactivate('_mock_glpk_direct')
Note: See TracBrowser for help on using the repository browser.