#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2022
#  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 io import StringIO

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'%<number>) == <number>) 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 io_options.items()))

        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 range(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
