mipgen.cc

Go to the documentation of this file.

/*
 *  Copyright 2007 Martin von Gagern
 *
 *
 *  This file is part of mqn2mps.
 *
 *  mqn2mps is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  mqn2mps is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */


#include <string>
#include <vector>
#include <map>
#include <set>
#include <sstream>
#include <iostream>
#include <fstream>
#include <limits>

#include "driver.hh"
#include "bounds.hh"
#include "mipgen.hh"
#include "euclid.hh"
#include "triplet.hh"
#include "vdiff.hh"

/**
 * @file
 * Implementation of class mipgen.
 */

/**
 * Constructor.
 * @param d the driver which will be accessed read-only for settings
 *          and input data.
 */
mipgen::mipgen(const driver& d) : drv(d) {
}

/**
 * Destructor.
 * It is only present so that derived classes can be destructed
 * cleanly in the presence of polymorphism.
 */
mipgen::~mipgen() {
}

/**
 * Modulo arithmetic as the mathematician knows it.
 * The C++ idea of modulo arithmetic is a bit different from most
 * mathematicians' ideas when it comes to negative numbers. This
 * function implements the mathematician's idea of a modulo
 * operation.
 *
 * @param a the divisor, may be negative.
 * @param b the dividend (modulus), must be positive.
 * @return a number @e c such that @f$ 0 \leq c < b @f$ and
 *         @f$ a = b \cdot d + c @f$ for some integer @e d.
 */
int mipgen::realmod(int a, int b) {
  if (a >= 0) return a % b; // = a - (a / b) * b;
  else return a + ((-a - 1)/b + 1) * b;
}

/**
 * Control problem formulation by calling different methods.
 * Look at the implementation to see in which methods are called in
 * which order.
 */
void mipgen::formulate() {
  constraints = drv.constraints;
  modfields();
  if (drv.multiply_k) multiply_k();
  if (drv.gamma_constraint) append_gamma();
  comdenom();
  get_indices();
  set_bounds();
  copy_row_names();
  collect_fields();
  copy_col_names();
  // unique_names(row_names);
  // unique_names(col_names);
  make_matrix(row_names.size(), col_names.size());
  fill_matrix();
}

/**
 * Add columns to express modulo arithmetic.
 * For each constraint that has a nonzero constraint::mod, a field is
 * added to mipgen::fields to express any integral multiple of the
 * modulus.
 */
void mipgen::modfields() {
  fields.clear();
  for (std::vector<constraint>::const_iterator
         i = drv.constraints.begin(), e = drv.constraints.end(); i != e; ++i) {
    if (i->mod <= 0) continue;
    std::ostringstream ss;
    ss << "_" << i->name << "mod" << i->mod;
    fields.push_back(field());
    fields.back().name = ss.str();
    fields.back().values.resize(drv.names.size());
    fields.back().values[drv.indices.find(i->name)->second] = i->mod;
  }
  modcount = fields.size();
  fields.insert(fields.end(), drv.fields.begin(), drv.fields.end());
}

/**
 * Multiply all @c n_* by @c k.
 * Every quantum number starting in @c n_ will be multiplied by the
 * quantum number starting with @c k. driver::set_names ensures that
 * such a field @c k exists whenever driver::multiply_k is set.
 */
void mipgen::multiply_k() {
  int k = drv.indices.find("k")->second;
  for (std::vector<field>::iterator
         i = fields.begin() + modcount, e = fields.end(); i != e; ++i) {
    for (size_t j = 0; j < drv.names.size(); ++j) {
      if (drv.names[j].find("n_") == 0)
        i->values[j] *= i->values[k];
    }
  }
}

/**
 * Formulate gamma constraint.
 *
 * If driver::gamma_constraint is set, a special rule holds for the
 * field component labeled @c gamma: no solution may have a single
 * fractional value as the only nonzero value in this column. All
 * other combinations like all zeros, multiple frations or at least
 * one integer are acceptable.
 *
 * This is expressed in the MIP by assigning a new number
 * @f$ n_\gamma @f$ to each field:
 * @f[
 *   n_\gamma = \begin{cases}
 *     0 & \quad q_\gamma = 0 \\
 *     1 & \quad q_\gamma \in \mathbb Q\setminus\mathbb Z \\
 *     2 & \quad q_\gamma \in \mathbb Z\setminus\{0\}
 *   \end{cases}
 * @f]
 * Now the constraint can be reformulated thus: the sum of all the
 * @f$ n_\gamma @f$ in a solution may never be exactly 1. To express
 * this, a new binary decision variable @f$ x_\gamma \in \{0,1\} @f$
 * is introduced, where 0 stands for the all zero case, and 1 for a
 * sum of at least 2. With this variable, the following two
 * inequalities constrain the sum to be not exactly 1.
 * @f{align*}
 *  -2x_\gamma+\sum_an_\gamma^{(a)} &\geq 0 \\
 *  -\infty x_\gamma+\sum_an_\gamma^{(a)} &\leq 0
 * @f}
 *
 * This method adds two new components to each existing field, and
 * adds a new field to represent @f$ x_\gamma @f$ (called @c _gamma in
 * the output). Instead of @f$ \infty @f$, twice the maximum order is
 * used, which is infinite enough here.
 */
void mipgen::append_gamma() {
  int g = drv.indices.find("gamma")->second;
  for (std::vector<field>::iterator
         i = fields.begin(), e = fields.end(); i != e; ++i) {
    rational gamma = i->values[g];
    int code;
    if (gamma.nom == 0) code = 0; // zero
    else if (gamma.nom % gamma.denom == 0) code = 2; // integral
    else code = 1; // fractional
    i->values.push_back(code);
    i->values.push_back(code);
  }
  fields.insert(fields.begin(), field());
  fields[0].name = "_gamma";
  fields[0].values.resize(drv.names.size());
  fields[0].values.push_back(-2);
  fields[0].values.push_back(-2*drv.max_order);
  ++modcount;
}

/**
 * Use common denominator.
 * The input is given in rational numbers, whereas the output will be
 * in integers. This method calculates the common denominator of each
 * row and expands all fractions to it. From now on, all calculations
 * can take place on the nominators solely.
 */
void mipgen::comdenom() {
  for (std::vector<constraint>::iterator
         i = constraints.begin(), e = constraints.end(); i != e; ++i) {
    int index = drv.indices.find(i->name)->second;
    int comdenom = i->goal.denom;
    for (std::vector<field>::iterator
           j = fields.begin(), ej = fields.end(); j != ej; ++j)
      comdenom = lcm(comdenom, j->values[index].denom);
    i->goal.set_denom(comdenom);
    for (std::vector<field>::iterator
           j = fields.begin(), ej = fields.end(); j != ej; ++j)
      j->values[index].set_denom(comdenom);
  }
}

/**
 * Build list of component indices.
 * The order of values in the field descriptions does not follow the
 * order of rows in the output. This method determines the
 * corresponding index into field::values for every output row. For
 * the first row, which constrains the maximum order, there is no
 * corresponding component. -1 is stored instead.
 */
void mipgen::get_indices() {
  indices.clear();
  indices.push_back(-1);
  if (drv.gamma_constraint) {
    indices.push_back(drv.names.size());
    indices.push_back(drv.names.size()+1);
  }
  for (std::vector<constraint>::const_iterator
         i = drv.constraints.begin(), e = drv.constraints.end(); i != e; ++i) {
    indices.push_back(drv.indices.find(i->name)->second);
  }
}

/**
 * Store bounds for all rows and columns.
 *
 * The row bounds are set to these values:
 *  -# the order is constrained between driver::min_order and
 *     driver::max_order
 *  -# the gamma rows are constrained as described in append_gamma()
 *  -# the input constraints are copied as strict equalities
 *
 * The column bounds are set to these values:
 *  -# the gamma column is constrained between 0 and 1
 *  -# the modulo columns are unconstrained
 *  -# the columns for input fields are constrained to be
 *     nonnegative.
 *
 * For most unconstrained ends of a bounds pair, an implied loose
 * bound is given as a hint, to allow derived classes to use ranges if
 * they can be handled better.
 *
 * Difference vectors are not yet encoded at the time when this method
 * is invoked.
 */
void mipgen::set_bounds() {
  row_bounds.clear();
  row_bounds.push_back(bounds(drv.min_order, drv.max_order)); // order
  if (drv.gamma_constraint) {
    row_bounds.push_back(bounds(0, 2*drv.max_order));
    row_bounds.back().upper_type = bounds::hint;
    row_bounds.push_back(bounds(-2*drv.max_order, 0));
    row_bounds.back().lower_type = bounds::hint;
  }
  for (std::vector<constraint>::iterator
         i = constraints.begin(), e = constraints.end(); i != e; ++i)
    row_bounds.push_back(bounds(i->goal.nom));

  col_bounds.clear();
  for (size_t i = 0; i < fields.size(); ++i) {
    if (drv.gamma_constraint && i == 0)
      col_bounds.push_back(bounds(0, 1));
    else if (i < modcount)
      col_bounds.push_back(bounds());
    else
      col_bounds.push_back(bounds(0, bounds::strict,
                                  drv.max_order, bounds::hint));
  }
}

/**
 * Fill list of row names.
 *
 * The names of the @c _MaxOrder row and the gamma rows are hardcoded
 * into this method. All other names come from the cooresponding
 * constraints.
 *
 * Difference vectors are not yet encoded at the time when this method
 * is invoked.
 */
void mipgen::copy_row_names() {  
  row_names.clear();
  row_names.push_back("_MaxOrder");
  if (drv.gamma_constraint) {
    row_names.push_back("_gamma1");
    row_names.push_back("_gamma2");
  }
  for (std::vector<constraint>::const_iterator
         i = constraints.begin(), e = constraints.end(); i != e; ++i)
    row_names.push_back(i->name);
}

/**
 * Transform using difference vectors.
 *
 * The input often contains multiple representations of the same
 * field, a so called multiplet. These multipletts are represented as
 * different field descriptions sharing a common name. For the
 * solution it is unimportant, which element of a multiplet is
 * actually used, as long as a valid solution does exist. Together
 * with the fact that multiplets tend to follow certain patterns, the
 * following encoding gan improve the problem description.
 *
 * For every multiplet, its first description is used as
 * representative. The others are encoded as differences in relation
 * to this representative. If another multiplet follows the same
 * pattern and thus results in the same set of difference vecotrs, the
 * same difference vectors can be shared between the two multiplets. A
 * new line for each set of difference vectors ensures that no
 * solution will contain more difference vectors than matching
 * multiplet representatives.
 *
 * At the time this method is invoked, bounds have been set and row
 * names filled in. This method will update the stored information.
 */
void mipgen::collect_fields() {
  // sort fields by name, collecting multiple fields with the same
  // name into a single vector.
  std::map<std::string, std::vector<const field*> > mf;
  for (std::vector<field>::iterator
         i = fields.begin() + modcount, e = fields.end(); i != e; ++i)
    mf[i->name].push_back(&*i);

  std::map<vdiff, int> md; // map difference sets to index of their row
  std::vector<field> of; // original offset fields
  std::vector<field> df; // difference fields
  int cnt = 0;
  int rows = row_names.size();
  for (std::map<std::string, std::vector<const field*> >::iterator
         i = mf.begin(), e = mf.end(); i != e; ++i) {
    if (i->second.size() == 1) { // this is a single vector, no multiplet
      of.push_back(*i->second[0]);
      continue;
    }
    // now we have a real multiplet, and let vdiff calculate and
    // compare the difference vectors.
    std::pair<std::map<vdiff, int>::iterator, bool> ires =
      md.insert(std::make_pair(vdiff(i->second), rows));
    if (ires.second) {
      // this is a new difference set, we need to create a row
      std::ostringstream ss;
      ss << "_d" << ++cnt;
      ++rows;
      indices.push_back(rows - 1);
      row_names.push_back(ss.str());
      row_bounds.push_back(bounds(0, bounds::strict, // non-negative
                                  drv.max_order, bounds::hint));
      for (size_t col = 0; col < ires.first->first.num_cols(); ++col) {
        std::ostringstream ssi;
        ssi << ss.str() << "_" << col+1;
        df.push_back(ires.first->first[col]);
        df.back().name = ssi.str();
        df.back().values.resize(rows);
        // a difference row has -1 for every associated difference vector:
        df.back().values[rows-1] = -1;
      }
    }
    else {
      // we already know this difference set, share difference vectors
      row_bounds[ires.first->second].upper_value += drv.max_order;
    }
    of.push_back(*i->second[0]);
    of.back().values.resize(rows);
    // a difference row has +1 for every associated representative:
    of.back().values[ires.first->second] = 1;
  }

  // Now we have the original fields of and the differences df.
  // Time to combine them into a single list of fields again.
  fields.resize(modcount); // remove old field descriptions
  fields.insert(fields.end(), df.begin(), df.end()); // add all of df
  fields.insert(fields.end(), of.begin(), of.end()); // add all of of
  for (std::vector<field>::iterator
         i = fields.begin(), e = fields.end(); i != e; ++i)
    i->values.resize(rows); // ensure all fields have correct size
  col_bounds.resize(modcount); // remove old column bounds
  for (std::vector<field>::iterator
         i = df.begin(), e = df.end(); i != e; ++i) {
    size_t row;
    for (row = 0; !i->values[row]; ++row) ; // find associated row
    // set bounds for difference cols to nonnegative with hint
    col_bounds.push_back(bounds(0, bounds::strict,
                                row_bounds[row].upper_value, bounds::hint));
  }
  for (std::vector<field>::iterator
         i = of.begin(), e = of.end(); i != e; ++i) {
    // set bounds for representatives to nonnegative with hint
    col_bounds.push_back(bounds(0, bounds::strict,
                                drv.max_order, bounds::hint));
  }
  modcount += df.size(); // the df are generated columns as well
}

/**
 * Fill list of column names.
 * These names are simply copied from mipgen::fields.
 */
void mipgen::copy_col_names() {  
  col_names.clear();
  for (std::vector<field>::iterator
         i = fields.begin(), e = fields.end(); i != e; ++i)
    col_names.push_back(i->name);
}

/**
 * Ensure uniqueness of names in a list.
 *
 * Before mipgen::collect_fields() was introduced, this method used to
 * ensure unique output names for the different columns of the same
 * multiplet. At the moment this method is unused.
 *
 * @param[in,out] names the list of names to be made unique.
 */
void mipgen::unique_names(std::vector<std::string>& names) {
  std::map<std::string, int> count;
  for (std::vector<std::string>::iterator
         i = names.begin(), e = names.end(); i != e; ++i) {
    std::ostringstream ss;
    ss << *i << "_" << ++count[*i];
    *i = ss.str();
  }
}

/**
 * Create a matrix data structure.
 *
 * Derived classes may wish to override this method.
 * The base implementation prepares mipgen::matrix by clearing any
 * data it might contain and adjusting its size according to the
 * request.
 *
 * @param rows the number of rows the coefficient matrix has.
 * @param cols the number of columns the coefficient matrix has.
 */
void mipgen::make_matrix(size_t rows, size_t cols) {
  matrix.clear();
  matrix.resize(rows, std::vector<int>(cols));
}

/**
 * Copy elements to matrix data structure.
 *
 * This class iterates over all matrix elements and calls
 * mipgen::set_element() for each one.
 */
void mipgen::fill_matrix() {
  for (size_t col = 0; col < col_names.size(); ++col) {
    set_element(0, col, col < modcount ? 0 : 1);
    for (size_t row = 1; row < row_names.size(); ++row) {
      set_element(row, col, fields[col].values[indices[row]].nom);
    }
  }
}

/**
 * Set a single matrix element.
 *
 * Derived classes may wish to override this method.
 * The base implementation stores the value in mipgen::matrix.
 * If a derived classes wishes to omit zero entries from its data
 * structure, it should check that case itself. This method will be
 * called for all matrix elements, zero and nonzero, in column major
 * order.
 *
 * @param row the row index of the matrix element.
 * @param col the column index of the matrix element.
 * @param val the value of the matrix element.
 */
void mipgen::set_element(size_t row, size_t col, int val) {
  matrix[row][col]=val;
}

/**
 * Write data to filename.
 *
 * The default implementation relies on write(std::ostream&) to do the
 * actual output. A derived class must override either method.
 *
 * @param filename the output file name, or 0 for standard output.
 */
void mipgen::write(const char* filename) {
  if (filename == 0) {
    write(std::cout);
  }
  else {
    std::ofstream out(filename);
    write(out);
    out.close();
  }
}

/**
 * Write data to stream.
 *
 * The default implementation prints the problem in a tabulated form
 * that can be opened using spreadsheet applications. This
 * functionality is accessibly using the @c --format=cvs command line
 * switch. For this use, most data is read from the driver and not
 * from the transformed data structures stored inside mipgen itself.
 */
void mipgen::write(std::ostream& out) {
  out << "Name:";
  for (std::vector<std::string>::const_iterator
         i = drv.names.begin(), e = drv.names.end(); i != e; ++i) {
    out << '\t' << *i;
  }
  out << std::endl << "Goal:";
  for (std::vector<std::string>::const_iterator
         i = drv.names.begin(), e = drv.names.end(); i != e; ++i) {
    out << '\t';
    for (std::vector<constraint>::const_iterator
           j = drv.constraints.begin(), je = drv.constraints.end();
         j != je; ++j) {
      if (*i == j->name) {
        out << j->goal;
        break;
      }
    }
  }
  out << std::endl << "Mod:";
  for (std::vector<std::string>::const_iterator
         i = drv.names.begin(), e = drv.names.end(); i != e; ++i) {
    out << '\t';
    for (std::vector<constraint>::const_iterator
           j = drv.constraints.begin(), je = drv.constraints.end();
         j != je; ++j) {
      if (*i == j->name) {
        if (j->mod) out << j->mod;
        break;
      }
    }
  }
  out << std::endl;
  for (std::vector<field>::const_iterator
         i = drv.fields.begin(), e = drv.fields.end(); i != e; ++i) {
    out << i->name;
    for (size_t j = 0; j < drv.names.size(); ++j)
      out << '\t' << i->values[j];
    out << std::endl;
  }
}

---