# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC # Under the terms of Contract DE-NA0003525 with National Technology and # Engineering Solutions of Sandia, LLC, the U.S. Government retains certain # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ # # Problem Writer for (Free) MPS Format Files # import logging from six import iteritems, StringIO from six.moves import xrange from pyomo.common.gc_manager import PauseGC from pyomo.opt import ProblemFormat from pyomo.opt.base import AbstractProblemWriter, WriterFactory from pyomo.core.base import \ (SymbolMap, TextLabeler, NumericLabeler, Constraint, SortComponents, Var, value, SOSConstraint, Objective, ComponentMap, is_fixed) from pyomo.repn import generate_standard_repn logger = logging.getLogger('pyomo.core') def _no_negative_zero(val): """Make sure -0 is never output. Makes diff tests easier.""" if val == 0: return 0 return val def _get_bound(exp): if exp is None: return None if is_fixed(exp): return value(exp) raise ValueError("non-fixed bound or weight: " + str(exp)) @WriterFactory.register('mps', 'Generate the corresponding MPS file') class ProblemWriter_mps(AbstractProblemWriter): def __init__(self): AbstractProblemWriter.__init__(self, ProblemFormat.mps) # the MPS writer is responsible for tracking which variables are # referenced in constraints, so that one doesn't end up with a # zillion "unreferenced variables" warning messages. stored at # the object level to avoid additional method arguments. # dictionary of id(_VarData)->_VarData. self._referenced_variable_ids = {} # Keven Hunter made a nice point about using %.16g in his attachment # to ticket #4319. I am adjusting this to %.17g as this mocks the # behavior of using %r (i.e., float('%r'%) == ) with # the added benefit of outputting (+/-). The only case where this # fails to mock the behavior of %r is for large (long) integers (L), # which is a rare case to run into and is probably indicative of # other issues with the model. # *** NOTE ***: If you use 'r' or 's' here, it will break code that # relies on using '%+' before the formatting character # and you will need to go add extra logic to output # the number's sign. self._precision_string = '.17g' def __call__(self, model, output_filename, solver_capability, io_options): # Make sure not to modify the user's dictionary, # they may be reusing it outside of this call io_options = dict(io_options) # Skip writing constraints whose body section is # fixed (i.e., no variables) skip_trivial_constraints = \ io_options.pop("skip_trivial_constraints", False) # Use full Pyomo component names in the MPS file rather # than shortened symbols (slower, but useful for debugging). symbolic_solver_labels = \ io_options.pop("symbolic_solver_labels", False) output_fixed_variable_bounds = \ io_options.pop("output_fixed_variable_bounds", False) # If False, unused variables will not be included in # the MPS file. Otherwise, include all variables in # the bounds sections. include_all_variable_bounds = \ io_options.pop("include_all_variable_bounds", False) labeler = io_options.pop("labeler", None) # How much effort do we want to put into ensuring the # MPS file is written deterministically for a Pyomo model: # 0 : None # 1 : sort keys of indexed components (default) # 2 : sort keys AND sort names (over declaration order) file_determinism = io_options.pop("file_determinism", 1) # user defined orderings for variable and constraint # output row_order = io_options.pop("row_order", None) column_order = io_options.pop("column_order", None) # make sure the ONE_VAR_CONSTANT variable appears in # the objective even if the constant part of the # objective is zero force_objective_constant = \ io_options.pop("force_objective_constant", False) # Whether or not to include the OBJSENSE section in # the MPS file. Some solvers, like GLPK and CBC, # either throw an error or flat out ignore this # section (I assume the default is to minimize) skip_objective_sense = \ io_options.pop("skip_objective_sense", False) if len(io_options): raise ValueError( "ProblemWriter_mps passed unrecognized io_options:\n\t" + "\n\t".join("%s = %s" % (k,v) for k,v in iteritems(io_options))) if symbolic_solver_labels and (labeler is not None): raise ValueError("ProblemWriter_mps: Using both the " "'symbolic_solver_labels' and 'labeler' " "I/O options is forbidden") if symbolic_solver_labels: labeler = TextLabeler() elif labeler is None: labeler = NumericLabeler('x') # clear the collection of referenced variables. self._referenced_variable_ids.clear() if output_filename is None: output_filename = model.name + ".mps" # when sorting, there are a non-trivial number of # temporary objects created. these all yield # non-circular references, so disable GC - the # overhead is non-trivial, and because references # are non-circular, everything will be collected # immediately anyway. with PauseGC() as pgc: with open(output_filename, "w") as output_file: symbol_map = self._print_model_MPS( model, output_file, solver_capability, labeler, output_fixed_variable_bounds=output_fixed_variable_bounds, file_determinism=file_determinism, row_order=row_order, column_order=column_order, skip_trivial_constraints=skip_trivial_constraints, force_objective_constant=force_objective_constant, include_all_variable_bounds=include_all_variable_bounds, skip_objective_sense=skip_objective_sense) self._referenced_variable_ids.clear() return output_filename, symbol_map def _extract_variable_coefficients( self, row_label, repn, column_data, quadratic_data, variable_to_column): # # Linear # if len(repn.linear_coefs) > 0: for vardata, coef in zip(repn.linear_vars, repn.linear_coefs): self._referenced_variable_ids[id(vardata)] = vardata column_data[variable_to_column[vardata]].append((row_label, coef)) # # Quadratic # if len(repn.quadratic_coefs) > 0: quad_terms = [] for vardata, coef in zip(repn.quadratic_vars, repn.quadratic_coefs): self._referenced_variable_ids[id(vardata[0])] = vardata[0] self._referenced_variable_ids[id(vardata[1])] = vardata[1] quad_terms.append( (vardata, coef) ) quadratic_data.append((row_label, quad_terms)) # # Return the constant # return repn.constant def _printSOS(self, symbol_map, labeler, variable_symbol_map, soscondata, output_file): if hasattr(soscondata, 'get_items'): sos_items = list(soscondata.get_items()) else: sos_items = list(soscondata.items()) if len(sos_items) == 0: return output_file.write("SOS\n") level = soscondata.level # I think there are many flavors to the SOS # section in the Free MPS format. I'm going with # what Cplex and Gurobi seem to recognize output_file.write(" S%d %s\n" % (level, symbol_map.getSymbol(soscondata,labeler))) sos_template_string = " %s %"+self._precision_string+"\n" for vardata, weight in sos_items: weight = _get_bound(weight) if weight < 0: raise ValueError( "Cannot use negative weight %f " "for variable %s is special ordered " "set %s " % (weight, vardata.name, soscondata.name)) if vardata.fixed: raise RuntimeError( "SOSConstraint '%s' includes a fixed variable '%s'. This is " "currently not supported. Deactive this constraint in order to " "proceed." % (soscondata.name, vardata.name)) self._referenced_variable_ids[id(vardata)] = vardata output_file.write(sos_template_string % (variable_symbol_map.getSymbol(vardata), weight)) def _print_model_MPS(self, model, output_file, solver_capability, labeler, output_fixed_variable_bounds=False, file_determinism=1, row_order=None, column_order=None, skip_trivial_constraints=False, force_objective_constant=False, include_all_variable_bounds=False, skip_objective_sense=False): symbol_map = SymbolMap() variable_symbol_map = SymbolMap() # NOTE: we use createSymbol instead of getSymbol because we # know whether or not the symbol exists, and don't want # to the overhead of error/duplicate checking. # cache frequently called functions extract_variable_coefficients = self._extract_variable_coefficients create_symbol_func = SymbolMap.createSymbol create_symbols_func = SymbolMap.createSymbols alias_symbol_func = SymbolMap.alias variable_label_pairs = [] sortOrder = SortComponents.unsorted if file_determinism >= 1: sortOrder = sortOrder | SortComponents.indices if file_determinism >= 2: sortOrder = sortOrder | SortComponents.alphabetical # # Create variable symbols (and cache the block list) # all_blocks = [] variable_list = [] for block in model.block_data_objects(active=True, sort=sortOrder): all_blocks.append(block) for vardata in block.component_data_objects( Var, active=True, sort=sortOrder, descend_into=False): variable_list.append(vardata) variable_label_pairs.append( (vardata,create_symbol_func(symbol_map, vardata, labeler))) variable_symbol_map.addSymbols(variable_label_pairs) # and extract the information we'll need for rapid labeling. object_symbol_dictionary = symbol_map.byObject variable_symbol_dictionary = variable_symbol_map.byObject # sort the variable ordering by the user # column_order ComponentMap if column_order is not None: variable_list.sort(key=lambda _x: column_order[_x]) # prepare to hold the sparse columns variable_to_column = ComponentMap( (vardata, i) for i, vardata in enumerate(variable_list)) # add one position for ONE_VAR_CONSTANT column_data = [[] for i in xrange(len(variable_list)+1)] quadobj_data = [] quadmatrix_data = [] # constraint rhs rhs_data = [] # print the model name and the source, so we know # roughly where output_file.write("* Source: Pyomo MPS Writer\n") output_file.write("* Format: Free MPS\n") output_file.write("*\n") output_file.write("NAME %s\n" % (model.name,)) # # ROWS section # objective_label = None numObj = 0 onames = [] for block in all_blocks: gen_obj_repn = \ getattr(block, "_gen_obj_repn", True) # Get/Create the ComponentMap for the repn if not hasattr(block,'_repn'): block._repn = ComponentMap() block_repn = block._repn for objective_data in block.component_data_objects( Objective, active=True, sort=sortOrder, descend_into=False): numObj += 1 onames.append(objective_data.name) if numObj > 1: raise ValueError( "More than one active objective defined for input " "model '%s'; Cannot write legal MPS file\n" "Objectives: %s" % (model.name, ' '.join(onames))) objective_label = create_symbol_func(symbol_map, objective_data, labeler) symbol_map.alias(objective_data, '__default_objective__') if not skip_objective_sense: output_file.write("OBJSENSE\n") if objective_data.is_minimizing(): output_file.write(" MIN\n") else: output_file.write(" MAX\n") # This section is not recognized by the COIN-OR # MPS reader #output_file.write("OBJNAME\n") #output_file.write(" %s\n" % (objective_label)) output_file.write("ROWS\n") output_file.write(" N %s\n" % (objective_label)) if gen_obj_repn: repn = \ generate_standard_repn(objective_data.expr) block_repn[objective_data] = repn else: repn = block_repn[objective_data] degree = repn.polynomial_degree() if degree == 0: logger.warning("Constant objective detected, replacing " "with a placeholder to prevent solver failure.") force_objective_constant = True elif degree is None: raise RuntimeError( "Cannot write legal MPS file. Objective '%s' " "has nonlinear terms that are not quadratic." % objective_data.name) constant = extract_variable_coefficients( objective_label, repn, column_data, quadobj_data, variable_to_column) if force_objective_constant or (constant != 0.0): # ONE_VAR_CONSTANT column_data[-1].append((objective_label, constant)) if numObj == 0: raise ValueError( "Cannot write legal MPS file: No objective defined " "for input model '%s'." % str(model)) assert objective_label is not None # Constraints def constraint_generator(): for block in all_blocks: gen_con_repn = \ getattr(block, "_gen_con_repn", True) # Get/Create the ComponentMap for the repn if not hasattr(block,'_repn'): block._repn = ComponentMap() block_repn = block._repn for constraint_data in block.component_data_objects( Constraint, active=True, sort=sortOrder, descend_into=False): if (not constraint_data.has_lb()) and \ (not constraint_data.has_ub()): assert not constraint_data.equality continue # non-binding, so skip if constraint_data._linear_canonical_form: repn = constraint_data.canonical_form() elif gen_con_repn: repn = generate_standard_repn(constraint_data.body) block_repn[constraint_data] = repn else: repn = block_repn[constraint_data] yield constraint_data, repn if row_order is not None: sorted_constraint_list = list(constraint_generator()) sorted_constraint_list.sort(key=lambda x: row_order[x[0]]) def yield_all_constraints(): for constraint_data, repn in sorted_constraint_list: yield constraint_data, repn else: yield_all_constraints = constraint_generator for constraint_data, repn in yield_all_constraints(): degree = repn.polynomial_degree() # Write constraint if degree == 0: if skip_trivial_constraints: continue elif degree is None: raise RuntimeError( "Cannot write legal MPS file. Constraint '%s' " "has nonlinear terms that are not quadratic." % constraint_data.name) # Create symbol con_symbol = create_symbol_func(symbol_map, constraint_data, labeler) if constraint_data.equality: assert value(constraint_data.lower) == \ value(constraint_data.upper) label = 'c_e_' + con_symbol + '_' alias_symbol_func(symbol_map, constraint_data, label) output_file.write(" E %s\n" % (label)) offset = extract_variable_coefficients( label, repn, column_data, quadmatrix_data, variable_to_column) bound = constraint_data.lower bound = _get_bound(bound) - offset rhs_data.append((label, _no_negative_zero(bound))) else: if constraint_data.has_lb(): if constraint_data.has_ub(): label = 'r_l_' + con_symbol + '_' else: label = 'c_l_' + con_symbol + '_' alias_symbol_func(symbol_map, constraint_data, label) output_file.write(" G %s\n" % (label)) offset = extract_variable_coefficients( label, repn, column_data, quadmatrix_data, variable_to_column) bound = constraint_data.lower bound = _get_bound(bound) - offset rhs_data.append((label, _no_negative_zero(bound))) else: assert constraint_data.has_ub() if constraint_data.has_ub(): if constraint_data.has_lb(): label = 'r_u_' + con_symbol + '_' else: label = 'c_u_' + con_symbol + '_' alias_symbol_func(symbol_map, constraint_data, label) output_file.write(" L %s\n" % (label)) offset = extract_variable_coefficients( label, repn, column_data, quadmatrix_data, variable_to_column) bound = constraint_data.upper bound = _get_bound(bound) - offset rhs_data.append((label, _no_negative_zero(bound))) else: assert constraint_data.has_lb() if len(column_data[-1]) > 0: # ONE_VAR_CONSTANT = 1 output_file.write(" E c_e_ONE_VAR_CONSTANT\n") column_data[-1].append(("c_e_ONE_VAR_CONSTANT",1)) rhs_data.append(("c_e_ONE_VAR_CONSTANT",1)) # # COLUMNS section # column_template = " %s %s %"+self._precision_string+"\n" output_file.write("COLUMNS\n") cnt = 0 for vardata in variable_list: col_entries = column_data[variable_to_column[vardata]] cnt += 1 if len(col_entries) > 0: var_label = variable_symbol_dictionary[id(vardata)] for i, (row_label, coef) in enumerate(col_entries): output_file.write(column_template % (var_label, row_label, _no_negative_zero(coef))) elif include_all_variable_bounds: # the column is empty, so add a (0 * var) # term to the objective # * Note that some solvers (e.g., Gurobi) # will accept an empty column as a line # with just the column name. This doesn't # seem to work for CPLEX 12.6, so I am # doing it this way so that it will work for both var_label = variable_symbol_dictionary[id(vardata)] output_file.write(column_template % (var_label, objective_label, 0)) assert cnt == len(column_data)-1 if len(column_data[-1]) > 0: col_entries = column_data[-1] var_label = "ONE_VAR_CONSTANT" for i, (row_label, coef) in enumerate(col_entries): output_file.write(column_template % (var_label, row_label, _no_negative_zero(coef))) # # RHS section # rhs_template = " RHS %s %"+self._precision_string+"\n" output_file.write("RHS\n") for i, (row_label, rhs) in enumerate(rhs_data): # note: we have already converted any -0 to 0 by this point output_file.write(rhs_template % (row_label, rhs)) # SOS constraints SOSlines = StringIO() sos1 = solver_capability("sos1") sos2 = solver_capability("sos2") for block in all_blocks: for soscondata in block.component_data_objects( SOSConstraint, active=True, sort=sortOrder, descend_into=False): create_symbol_func(symbol_map, soscondata, labeler) level = soscondata.level if (level == 1 and not sos1) or \ (level == 2 and not sos2) or \ (level > 2): raise ValueError( "Solver does not support SOS level %s constraints" % (level)) # This updates the referenced_variable_ids, just in case # there is a variable that only appears in an # SOSConstraint, in which case this needs to be known # before we write the "bounds" section (Cplex does not # handle this correctly, Gurobi does) self._printSOS(symbol_map, labeler, variable_symbol_map, soscondata, SOSlines) # # BOUNDS section # entry_template = "%s %"+self._precision_string+"\n" output_file.write("BOUNDS\n") for vardata in variable_list: if include_all_variable_bounds or \ (id(vardata) in self._referenced_variable_ids): var_label = variable_symbol_dictionary[id(vardata)] if vardata.fixed: if not output_fixed_variable_bounds: raise ValueError( "Encountered a fixed variable (%s) inside an active " "objective or constraint expression on model %s, which is " "usually indicative of a preprocessing error. Use the " "IO-option 'output_fixed_variable_bounds=True' to suppress " "this error and fix the variable by overwriting its bounds " "in the MPS file." % (vardata.name, model.name)) if vardata.value is None: raise ValueError("Variable cannot be fixed to a value of None.") output_file.write((" FX BOUND "+entry_template) % (var_label, _no_negative_zero(value(vardata.value)))) continue # convert any -0 to 0 to make baseline diffing easier vardata_lb = _no_negative_zero(_get_bound(vardata.lb)) vardata_ub = _no_negative_zero(_get_bound(vardata.ub)) unbounded_lb = not vardata.has_lb() unbounded_ub = not vardata.has_ub() treat_as_integer = False if vardata.is_binary(): if (vardata_lb == 0) and (vardata_ub == 1): output_file.write(" BV BOUND %s\n" % (var_label)) continue else: # so we can add bounds treat_as_integer = True if treat_as_integer or vardata.is_integer(): # Indicating unbounded integers is tricky because # the only way to indicate a variable is integer # is using the bounds section. Thus, we signify # infinity with a large number (10E20) # * Note: Gurobi allows values like inf and -inf # but CPLEX 12.6 does not, so I am just # using a large value if not unbounded_lb: output_file.write((" LI BOUND "+entry_template) % (var_label, vardata_lb)) else: output_file.write(" LI BOUND %s -10E20\n" % (var_label)) if not unbounded_ub: output_file.write((" UI BOUND "+entry_template) % (var_label, vardata_ub)) else: output_file.write(" UI BOUND %s 10E20\n" % (var_label)) else: assert vardata.is_continuous() if unbounded_lb and unbounded_ub: output_file.write(" FR BOUND %s\n" % (var_label)) else: if not unbounded_lb: output_file.write((" LO BOUND "+entry_template) % (var_label, vardata_lb)) else: output_file.write(" MI BOUND %s\n" % (var_label)) if not unbounded_ub: output_file.write((" UP BOUND "+entry_template) % (var_label, vardata_ub)) # # SOS section # output_file.write(SOSlines.getvalue()) # Formatting of the next two sections comes from looking # at Gurobi and Cplex output # # QUADOBJ section # if len(quadobj_data) > 0: assert len(quadobj_data) == 1 # it looks like the COIN-OR MPS Reader only # recognizes QUADOBJ (Gurobi and Cplex seem to # be okay with this) output_file.write("QUADOBJ\n") #output_file.write("QMATRIX\n") label, quad_terms = quadobj_data[0] assert label == objective_label # sort by the sorted tuple of symbols (or column assignments) # for the variables appearing in the term quad_terms = sorted(quad_terms, key=lambda _x: \ sorted((variable_to_column[_x[0][0]], variable_to_column[_x[0][1]]))) for term, coef in quad_terms: # sort the term for consistent output var1, var2 = sorted(term, key=lambda _x: variable_to_column[_x]) var1_label = variable_symbol_dictionary[id(var1)] var2_label = variable_symbol_dictionary[id(var2)] # Don't forget that a quadratic objective is always # assumed to be divided by 2 if var1_label == var2_label: output_file.write(column_template % (var1_label, var2_label, _no_negative_zero(coef * 2))) else: # the matrix needs to be symmetric so split # the coefficient (but remember it is divided by 2) output_file.write(column_template % (var1_label, var2_label, _no_negative_zero(coef))) output_file.write(column_template % (var2_label, var1_label, _no_negative_zero(coef))) # # QCMATRIX section # if len(quadmatrix_data) > 0: for row_label, quad_terms in quadmatrix_data: output_file.write("QCMATRIX %s\n" % (row_label)) # sort by the sorted tuple of symbols (or # column assignments) for the variables # appearing in the term quad_terms = sorted(quad_terms, key=lambda _x: \ sorted((variable_to_column[_x[0][0]], variable_to_column[_x[0][1]]))) for term, coef in quad_terms: # sort the term for consistent output var1, var2 = sorted(term, key=lambda _x: variable_to_column[_x]) var1_label = variable_symbol_dictionary[id(var1)] var2_label = variable_symbol_dictionary[id(var2)] if var1_label == var2_label: output_file.write(column_template % (var1_label, var2_label, _no_negative_zero(coef))) else: # the matrix needs to be symmetric so split # the coefficient output_file.write(column_template % (var1_label, var2_label, _no_negative_zero(coef * 0.5))) output_file.write(column_template % (var2_label, var1_label, coef * 0.5)) output_file.write("ENDATA\n") # Clean up the symbol map to only contain variables referenced # in the active constraints **Note**: warm start method may # rely on this for choosing the set of potential warm start # variables vars_to_delete = set(variable_symbol_map.byObject.keys()) - \ set(self._referenced_variable_ids.keys()) sm_byObject = symbol_map.byObject sm_bySymbol = symbol_map.bySymbol var_sm_byObject = variable_symbol_map.byObject for varid in vars_to_delete: symbol = var_sm_byObject[varid] del sm_byObject[varid] del sm_bySymbol[symbol] del variable_symbol_map return symbol_map