MIPModel.java

package solvers;

import ilog.concert.IloException;
import ilog.concert.IloIntVar;
import ilog.concert.IloLinearNumExpr;
import ilog.concert.IloNumVar;
import ilog.cplex.IloCplex;
import java.util.ArrayList;
import java.util.Collection;
import java.util.HashSet;
import xmtrclusters.Problem;

/**
 * MIPModel implements a mixed-integer linear program for the clustering
 * problem.
 *
 * @author Paul A. Rubin (rubin@msu.edu)
 */
public final class MIPModel {
  private static final double HALF = 0.5;  // for rounding binary values

  private final Problem problem;  // the parent problem

  // Model components.
  private final IloCplex mip;     // the MIP model
  private final IloIntVar[][] x;  // x[i][j] = 1 if trans i & j are in the same
                                  // cluster; x[i][i] = 1 if i anchors a cluster
  private final IloIntVar[][] y;  // y[i][j] = 1 if trans i & user j are
                                  // in the same cluster
  private final IloNumVar[][] z;  // z[i][j] is the objective contribution of
                                  // the trans i / user j pairing
  private final IloNumVar[] q;    // q[j] is the solution quality for user j
  private final IloNumVar[][] v;  // v[i][j] = 1 if transmitters i and j
                                  // are in the same cluster

  /**
   * Constructor.
   * @param prob the problem to solve
   * @throws IloException if CPLEX throws an error building the model
   */
  public MIPModel(final Problem prob) throws IloException {
    problem = prob;
    int nU = problem.getNUsers();
    int nT = problem.getNTrans();
    double[] qMax = problem.getMaxQuality();
    // Initialize the model.
    mip = new IloCplex();
    // Add the variables.
    x = new IloIntVar[nT][nT];
    y = new IloIntVar[nT][nU];
    z = new IloNumVar[nT][nU];
    q = new IloNumVar[nU];
    v = new IloNumVar[nT][nT];
    for (int t = 0; t < nT; t++) {
      for (int t1 = 0; t1 < nT; t1++) {
        x[t][t1] = mip.boolVar("x_" + t + "_" + t1);
        v[t][t1] = mip.numVar(0, 1, "v_" + t + "_" + t1);
      }
      for (int u = 0; u < nU; u++) {
        y[t][u] = mip.boolVar("y_" + t + "_" + u);
        z[t][u] = mip.numVar(0, qMax[u], "z_" + t + "_" + u);
      }
    }
    for (int u = 0; u < nU; u++) {
      q[u] = mip.numVar(0, qMax[u], "q_" + u);
    }
    // The objective is to maximize the total userQuality.
    mip.addMaximize(mip.sum(q));
    // Make the x variables symmetric.
    for (int t = 0; t < nT; t++) {
      for (int t1 = t + 1; t1 < nT; t1++) {
        mip.addEq(x[t][t1], x[t1][t], "symmetry_" + t + "_" + t1);
      }
    }
    // Limit the number of anchors (hence the number of clusters).
    IloLinearNumExpr expr = mip.linearNumExpr();
    for (int t = 0; t < nT; t++) {
      expr.addTerm(1.0, x[t][t]);
    }
    mip.addLe(expr, problem.getMaxClusters(), "max_cluster_count");
    // Limit the size of every cluster.
    int maxS = problem.getMaxSize();
    for (int t = 0; t < nT; t++) {
      mip.addLe(mip.diff(mip.sum(x[t]), x[t][t]), maxS - 1,
                "cluster_size_" + t);
    }
    // Avoid multiple anchors in a cluster.
    // (These constraints may be unnecessary.)
    for (int t = 0; t < nT; t++) {
      for (int t1 = t + 1; t1 < nT; t1++) {
        mip.addLe(mip.sum(x[t][t], x[t][t1], x[t1][t1]), 2.0,
                  "one_anchor_" + t + "_" + t1);
      }
    }
    // Force each transmitter to be in a cluster with an anchor.
    for (int t = 0; t < nT; t++) {
      expr = mip.linearNumExpr();
      expr.addTerm(1.0, x[t][t]);
      for (int t1 = 0; t1 < nT; t1++) {
        if (t != t1) {
          mip.addLe(v[t][t1], x[t][t1], "v_def1_" + t + "_" + t1);
          mip.addLe(v[t][t1], x[t1][t1], "v_def2_" + t + "_" + t1);
          expr.addTerm(1.0, v[t][t1]);
        }
      }
      mip.addEq(expr, 1.0, "in_cluster_" + t);
    }
    // Cluster membership is transitive.
    for (int t = 0; t < nT; t++) {
      for (int t1 = 0; t1 < nT; t1++) {
        if (t != t1) {
          for (int t2 = 0; t2 < nT; t2++) {
            if (t2 != t && t2 != t1) {
              mip.addGe(x[t][t2], mip.diff(mip.sum(x[t][t1], x[t1][t2]), 1.0),
                        "transitive_" + t + "_" + t1 + "_" + t2);
            }
          }
        }
      }
    }
    // Antisymmetry: lowest index transmitter is the anchor.
    for (int t = 0; t < nT; t++) {
      for (int t1 = t + 1; t1 < nT; t1++) {
        mip.addLe(x[t1][t1], mip.diff(1.0, x[t][t1]),
                  "asymmetry_" + t + "_" + t1);
      }
    }
    // Put each user in with their favorite transmitter.
    int[] tau = problem.getBest();
    for (int u = 0; u < nU; u++) {
      mip.addEq(y[tau[u]][u], 1.0, "favorite_transmitter_" + u);
    }
    // Determine other transmitter-user pairings.
    for (int u = 0; u < nU; u++) {
      for (int t = 0; t < nT; t++) {
        if (t != tau[u]) {
          mip.addEq(y[t][u], x[tau[u]][t], "pairing_" + t + "_" + u);
        }
      }
    }
    // Define the z variables.
    for (int t = 0; t < nT; t++) {
      for (int u = 0; u < nU; u++) {
        mip.addLe(z[t][u], mip.prod(qMax[u], y[t][u]), "zdef1_" + t + "_" + u);
        expr = mip.linearNumExpr();
        expr.addTerm(1.0, z[t][u]);
        expr.addTerm(-1.0, q[u]);
        expr.addTerm(qMax[u], y[t][u]);
        mip.addLe(expr, qMax[u], "zdef2_" + t + "_" + u);
        expr = mip.linearNumExpr();
        expr.addTerm(1.0, z[t][u]);
        expr.addTerm(-1.0, q[u]);
        expr.addTerm(-qMax[u], y[t][u]);
        mip.addGe(expr, -qMax[u], "zdef3_" + t + "_" + u);
      }
    }
    // Linearize q.
    for (int u = 0; u < nU; u++) {
      expr = mip.linearNumExpr();
      double[] w = problem.getWeights(u);
      for (int i = 0; i < nT; i++) {
        expr.addTerm(w[i], q[u]);
        expr.addTerm(-w[i], z[i][u]);
        expr.addTerm(-w[i], y[i][u]);
      }
      mip.addEq(expr, 0.0, "qdef_" + u);
    }
  }

  /**
   * Exports the model to a file.
   * @param target the target file path/name
   * @throws IloException if CPLEX encounters an error during export
   */
  public void export(final String target) throws IloException {
    mip.exportModel(target);
  }

  /**
   * Solves the model.
   * @param sec time limit in seconds
   * @return the final objective value
   * @throws IloException if CPLEX encounters an error
   */
  public double solve(final double sec) throws IloException {
    mip.setParam(IloCplex.DoubleParam.TimeLimit, sec);
    mip.setParam(IloCplex.Param.Emphasis.MIP, 3);
    mip.solve();
    return mip.getObjValue();
  }

  /**
   * Gets the transmitter clusters in the final solution.
   * @return the collection of clusters
   * @throws IloException if CPLEX cannot extract the solution.
   */
  public Collection<Collection<Integer>> getClusters() throws IloException {
    // Make sure a solution exists.
    IloCplex.Status status = mip.getStatus();
    if (status != IloCplex.Status.Optimal
        && status != IloCplex.Status.Feasible) {
      throw new IllegalArgumentException("Cannot extract a"
                                         + " nonexistent solution!");
    }
    int nT = problem.getNTrans();
    // Extract the values of x[][] in the solution.
    boolean[][] xx = new boolean[nT][nT];
    for (int t = 0; t < nT; t++) {
      for (int t1 = t; t1 < nT; t1++) {
        xx[t][t1] = mip.getValue(x[t][t1]) > HALF;
        xx[t1][t] = xx[t][t1];
      }
    }
    // Identify the clusters.
    ArrayList<Collection<Integer>> clusters = new ArrayList<>();
    for (int t = 0; t < nT; t++) {
      // Look for anchors.
      if (xx[t][t]) {
        // t is an anchor; add the cluster containing it (ncluding itself).
        HashSet<Integer> c = new HashSet<>();
        for (int t1 = 0; t1 < nT; t1++) {
          if (xx[t][t1]) {
            c.add(t1);
          }
        }
        clusters.add(c);
      }
    }
    return clusters;
  }
}