Initial commit.

.gitignore

+/dist/
+*.xml
+*.properties
+*.class
+*.mf

src/models/Graph.java

+package models;
+
+import java.util.ArrayList;
+import java.util.Random;
+
+/**
+ * Graph contains an instance of a layered graph.
+ *
+ * Vertices are numbered from 0 to n-1, where n is the number of vertices in
+ * the graph. Vertices in layer i are *labeled* from 1 to n_i, where n_i is
+ * the number of vertices in layer i.
+ *
+ * @author Paul A. Rubin (rubin@msu.edu)
+ */
+public final class Graph {
+  // Dimension labels for the layers matrix.
+  private static final int SIZE = 0;   // number of vertices in the layer
+  private static final int FIRST = 1;  // index of the first vertex in the layer
+  private static final int LAST = 2;   // index of the last vertex in the layer
+
+  // Dimension labels for the edges matrix.
+  private static final int TAIL = 0;   // index of the tail node
+  private static final int HEAD = 1;   // index of the head node
+  private static final int FROM = 2;   // index of the layer containing the tail
+
+  private final int nVertices;         // number of vertices
+  private final int nLayers;           // number of layers
+  private final int nEdges;            // number of edges
+  private final int minLayerSize;      // minimum layer size
+  private final int maxLayerSize;      // maximum layer size
+  private final int[][] layers;        // info about each layer
+  private final int[][] edges;         // info about each edge
+
+  /**
+   * Constructs a random graph.
+   * @param nL the number of layers to use
+   * @param minV the minimum number of vertices in a layer
+   * @param maxV the maximum number of vertices in a layer
+   * @param p the proportion of possible arcs that actually exist
+   * @param seed the seed for the random number generator
+   */
+  @SuppressWarnings({ "checkstyle:magicnumber" })
+  public Graph(final int nL, final int minV, final int maxV, final double p,
+               final long seed) {
+    nLayers = nL;
+    // Set up a random number generator.
+    Random rng = new Random(seed);
+    // Set up the layers matrix.
+    layers = new int[nLayers][3];
+    // Assign nodes to each layer.
+    int id = 0;
+    int mv = maxV - minV + 1;
+    int s1 = Integer.MAX_VALUE;     // tracks min layer size
+    int s2 = 0;                     // tracks max layer size
+    for (int j = 0; j < nL; j++) {
+      // Randomly choose the number of nodes to put in the layer.
+      int n = minV + rng.nextInt(mv);
+      s1 = Math.min(s1, n);
+      s2 = Math.max(s2, n);
+      // Assign the nodes.
+      layers[j][SIZE] = n;
+      layers[j][FIRST] = id;
+      layers[j][LAST] = id + n - 1;
+      id += n;
+    }
+    nVertices = id;
+    minLayerSize = s1;
+    maxLayerSize = s2;
+    // Since the number of edges is initially undetermined, use two array lists
+    // to track their endpoints.
+    ArrayList<Integer> heads = new ArrayList<>();
+    ArrayList<Integer> tails = new ArrayList<>();
+    ArrayList<Integer> from = new ArrayList<>();
+    // Consider each possible edge and randomly decide whether to include it.
+    for (int i = 0; i < nLayers - 1; i++) {
+      for (int j = layers[i][FIRST]; j <= layers[i][LAST]; j++) {
+        for (int k = layers[i + 1][FIRST]; k <= layers[i + 1][LAST]; k++) {
+          if (rng.nextDouble() <= p) {
+            // Take the edge from j to k.
+            tails.add(j);
+            heads.add(k);
+            from.add(i);
+          }
+        }
+      }
+    }
+    // Convert the edge info into a matrix.
+    edges = new int[tails.size()][3];
+    for (int i = 0; i < tails.size(); i++) {
+      edges[i][TAIL] = tails.get(i);
+      edges[i][HEAD] = heads.get(i);
+      edges[i][FROM] = from.get(i);
+    }
+    nEdges = edges.length;
+  }
+
+  /**
+   * Gets the number of vertices in the graph.
+   * @return the vertex count
+   */
+  public int getNVertices() {
+    return nVertices;
+  }
+
+  /**
+   * Gets the number of layers in the graph.
+   * @return the layer count
+   */
+  public int getNLayers() {
+    return nLayers;
+  }
+
+  /**
+   * Gets the number of edges in the graph.
+   * @return the edge count
+   */
+  public int getNEdges() {
+    return nEdges;
+  }
+
+  /**
+   * Gets the maximum number of vertices in any layer.
+   * This is also the largest possible label value.
+   * @return the maximum layer size
+   */
+  public int getMaxLayerSize() {
+    return maxLayerSize;
+  }
+
+  /**
+   * Gets the size (vertex count) of a layer.
+   * @param layerIndex the index of the layer
+   * @return the number of vertices in the layer
+   */
+  public int getLayerSize(final int layerIndex) {
+    return layers[layerIndex][SIZE];
+  }
+
+  /**
+   * Gets the sizes of all layers.
+   * @return the layer sizes
+   */
+  public int[] getLayerSizes() {
+    return getColumn(layers, SIZE);
+  }
+
+  /**
+   * Gets the lowest vertex index in a layer.
+   * @param layerIndex the index of the target layer
+   * @return the index of the first vertex in the layer
+   */
+  public int getLayerStart(final int layerIndex) {
+    return layers[layerIndex][FIRST];
+  }
+
+  /**
+   * Gets the first vertex in every layer.
+   * @return the initial vertices of layers
+   */
+  public int[] getLayerStarts() {
+    return getColumn(layers, FIRST);
+  }
+
+  /**
+   * Gets the highest vertex index in a layer.
+   * @param layerIndex the index of the target layer
+   * @return the index of the last vertex in the layer
+   */
+  public int getLayerEnd(final int layerIndex) {
+    return layers[layerIndex][LAST];
+  }
+
+  /**
+   * Gets the last vertex in every layer.
+   * @return the final vertices of layers
+   */
+  public int[] getLayerEnds() {
+    return getColumn(layers, LAST);
+  }
+
+  /**
+   * Gets the head vertex of an edge.
+   * @param edgeIndex the index of the edge
+   * @return the index of the head vertex
+   */
+  public int getHead(final int edgeIndex) {
+    return edges[edgeIndex][HEAD];
+  }
+
+  /**
+   * Gets the tail vertex of an edge.
+   * @param edgeIndex the index of the edge
+   * @return the index of the tail vertex
+   */
+  public int getTail(final int edgeIndex) {
+    return edges[edgeIndex][TAIL];
+  }
+
+  /**
+   * Generates a summary of the graph.
+   * @return a string summarizing the graph
+   */
+  public String toString() {
+    StringBuilder sb = new StringBuilder("Graph dimensions:\n");
+    sb.append("\tOverall vertex count = ").append(nVertices).append(".\n")
+      .append("\tNumber of layers = ").append(nLayers).append(".\n")
+      .append("\tSmallest layer has ").append(minLayerSize)
+      .append(" vertices.\n\tLargest layer has ").append(maxLayerSize)
+      .append(" vertices.\n\tNumber of edges = ").append(nEdges).append(".");
+    return sb.toString();
+  }
+
+  /**
+   * Evaluates a candidate solution.
+   * @param labels the labels for every node
+   * @return a string report of the solution
+   */
+  public String evaluate(final int[] labels) {
+    StringBuilder sb = new StringBuilder();
+    // Verify that no label exceeds the layer maximum and that every label
+    // in a layer is used at most once.
+    for (int k = 0; k < nLayers; k++) {
+      int s = layers[k][SIZE];
+      boolean[] used = new boolean[s + 1];
+      for (int i = layers[k][FIRST]; i <= layers[k][LAST]; i++) {
+        int n = labels[i];
+        if (n < 1) {
+          sb.append("Invalid label ").append(n).append(" encountered in layer ")
+            .append(k).append("\n");
+        } else if (n > s) {
+          sb.append("Label ").append(n).append(" in layer ").append(k)
+            .append(" exceeds layer size ").append(s).append(".\n");
+        } else if (used[n]) {
+          sb.append("Label ").append(n).append(" already used in layer ")
+            .append(k).append(".\n");
+        } else {
+          used[n] = true;
+        }
+      }
+    }
+    // Evaluate the objective value.
+    sb.append("Objective value = ").append(getObjValue(labels)).append(".");
+    return sb.toString();
+  }
+
+  /**
+   * Computes the objective value of a vector of vertex labels.
+   * @param labels the vertex labels
+   * @return the objective value
+   */
+  public long getObjValue(final int[] labels) {
+    long obj = 0;
+    for (int e = 0; e < edges.length; e++) {
+      int a = labels[edges[e][TAIL]];
+      a -= labels[edges[e][HEAD]];
+      obj += a * a;
+    }
+    return obj;
+  }
+
+  /**
+   * Extracts a specified column from an integer matrix (either the edge
+   * matrix or the layer matrix).
+   * @param m the source matrix
+   * @param col the column to extract
+   * @return the extracted column
+   */
+  private int[] getColumn(final int[][] m, final int col) {
+    int[] c = new int[m.length];
+    for (int i = 0; i < m.length; i++) {
+      c[i] = m[i][col];
+    }
+    return c;
+  }
+}

src/models/MIP.java

+package models;
+
+import ilog.concert.IloException;
+import ilog.concert.IloIntVar;
+import ilog.concert.IloLinearNumExpr;
+import ilog.concert.IloNumVar;
+import ilog.concert.IloQuadNumExpr;
+import ilog.cplex.IloCplex;
+import java.util.ArrayList;
+
+/**
+ * MIP implements a mixed-integer programming model for the vertex
+ * labeling problem.
+ *
+ * The model is based on Rob Pratt's solution to the OR SE question.
+ *
+ * @author Paul A. Rubin (rubin@msu.edu)
+ */
+public final class MIP implements AutoCloseable {
+  private static final double HALF = 0.5;  // rounding cutoff
+  private final IloCplex mip;              // MIP model
+  private final IloIntVar[][] label;       // indicators for label assignments
+    // label[i][j] = 1 if vertex i is assigned label j + 1
+
+  /**
+   * Constructs the MIP model.
+   * @param g the graph to be processed
+   * @throws IloException if CPLEX balks at any step
+   */
+  public MIP(final Graph g) throws IloException {
+    // Get key dimensions.
+    int nVertices = g.getNVertices();
+    int nEdges = g.getNEdges();
+    int maxLabel = g.getMaxLayerSize();
+    int nLayers = g.getNLayers();
+    // Initialize the MIP model object.
+    mip = new IloCplex();
+    // Create indicator variables for label assignments.
+    label = new IloIntVar[nVertices][maxLabel];
+    for (int i = 0; i < nVertices; i++) {
+      for (int j = 0; j < maxLabel; j++) {
+        label[i][j] = mip.boolVar("x_" + i + "_" + (j + 1));
+      }
+    }
+    // Disallow assignments to labels greater than the layer size.
+    for (int k = 0; k < nLayers; k++) {
+      int n = g.getLayerSize(k);
+      if (n < maxLabel) {
+        for (int i = g.getLayerStart(k); i <= g.getLayerEnd(k); i++) {
+          for (int j = n; j < maxLabel; j++) {
+            label[i][j].setUB(0);
+          }
+        }
+      }
+    }
+    // The objective is to minimize the sum of the edge costs.
+    IloQuadNumExpr obj = mip.quadNumExpr();
+    for (int e = 0; e < nEdges; e++) {
+      // Get the endpoints of edge e.
+      int head = g.getHead(e);
+      int tail = g.getTail(e);
+      for (int i = 0; i < maxLabel; i++) {
+        for (int j = 0; j < maxLabel; j++) {
+          // Compute the cost of this pair of labels.
+          double c = (i - j) * (i - j);
+          // Add the cost times the product of the two indicators.
+          obj.addTerm(c, label[head][i], label[tail][j]);
+        }
+      }
+    }
+    mip.addMinimize(obj);
+    // Every vertex needs to be assigned exactly one label.
+    for (int i = 0; i < nVertices; i++) {
+      mip.addEq(mip.sum(label[i]), 1.0);
+    }
+    // Every valid label (1, ..., layer size) needs to be assigned to exactly
+    // one vertex in the layer.
+    for (int k = 0; k < nLayers; k++) {
+      // Get the size of the layer. This is the maximum valid label.
+      int n = g.getLayerSize(k);
+      for (int j = 0; j < n; j++) {
+        // Assign label j to exactly one vertex in the layer.
+        IloLinearNumExpr expr = mip.linearNumExpr();
+        for (int i = g.getLayerStart(k); i <= g.getLayerEnd(k); i++) {
+          expr.addTerm(1.0, label[i][j]);
+        }
+        mip.addEq(expr, 1.0);
+      }
+    }
+  }
+
+  /**
+   * Attempts to solve the MIP model.
+   * @param timeLimit the time limit in seconds for the solver
+   * @return the final solver status
+   * @throws IloException if CPLEX blows up
+   */
+  public IloCplex.Status solve(final double timeLimit) throws IloException {
+    mip.setParam(IloCplex.DoubleParam.TimeLimit, timeLimit);
+    mip.solve();
+    return mip.getStatus();
+  }
+
+  /**
+   * Closes the MIP model.
+   */
+  @Override
+  public void close() {
+    mip.close();
+  }
+
+  /**
+   * Gets the labels assigned by the incumbent solution to the MIP model.
+   * @return the incumbent labels
+   * @throws IloException if there is no incumbent or it cannot be accessed
+   */
+  public int[] getLabels() throws IloException {
+    IloCplex.Status status = mip.getStatus();
+    if (status == IloCplex.Status.Optimal
+        || status == IloCplex.Status.Feasible) {
+      int[] sol = new int[label.length];
+      for (int i = 0; i < label.length; i++) {
+        double[] z = mip.getValues(label[i]);
+        for (int j = 0; j < z.length; j++) {
+          if (z[j] > HALF) {
+            sol[i] = j + 1;
+            break;
+          }
+        }
+      }
+      return sol;
+    } else {
+      throw new IloException("Cannot recover solution when status = " + status);
+    }
+  }
+
+  /**
+   * Sets a starting solution for the MIP.
+   * @param labels the labels to use for the starting solution
+   * @throws IloException if the MIP start is rejected
+   */
+  public void setStart(final int[] labels) throws IloException {
+    int mx = label[0].length;  // max # of possible labels
+    ArrayList<IloNumVar> varList = new ArrayList<>();
+    ArrayList<Double> valList = new ArrayList<>();
+    for (int i = 0; i < label.length; i++) {
+      for (int j = 0; j < mx; j++) {
+        varList.add(label[i][j]);
+        double x = (labels[i] == (j + 1)) ?  1.0 : 0.0;
+        valList.add(x);
+      }
+    }
+    IloNumVar[] vars = varList.toArray(new IloNumVar[varList.size()]);
+    double[] vals = valList.stream().mapToDouble(x -> x).toArray();
+    mip.addMIPStart(vars, vals, IloCplex.MIPStartEffort.Auto);
+  }
+}

src/models/RKGA.java

+package models;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Random;
+import org.uncommons.maths.random.MersenneTwisterRNG;
+import org.uncommons.maths.random.Probability;
+import org.uncommons.watchmaker.framework.CandidateFactory;
+import org.uncommons.watchmaker.framework.EvolutionaryOperator;
+import org.uncommons.watchmaker.framework.FitnessEvaluator;
+import org.uncommons.watchmaker.framework.PopulationData;
+import org.uncommons.watchmaker.framework.factories.AbstractCandidateFactory;
+import org.uncommons.watchmaker.framework.islands.IslandEvolution;
+import org.uncommons.watchmaker.framework.islands.IslandEvolutionObserver;
+import org.uncommons.watchmaker.framework.islands.RingMigration;
+import org.uncommons.watchmaker.framework.operators.DoubleArrayCrossover;
+import org.uncommons.watchmaker.framework.operators.EvolutionPipeline;
+import org.uncommons.watchmaker.framework.selection.RouletteWheelSelection;
+import org.uncommons.watchmaker.framework.termination.ElapsedTime;
+import org.uncommons.watchmaker.framework.termination.Stagnation;
+import org.uncommons.watchmaker.framework.termination.TargetFitness;
+import valueorderedmap.ValueOrderedMap;
+
+/**
+ * RKGA implements a random key genetic algorithm for the vertex numbering
+ * problem.
+ *
+ * @author Paul A. Rubin (rubin@msu.edu)
+ */
+public final class RKGA implements FitnessEvaluator<double[]> {
+  private static final long SECTOMS = 1000;
+  private final Graph graph;    // the graph to label
+  private final int nVertices;  // number of vertices in the graph
+  private final int nLayers;    // number of layers in the graph
+  private final int[] first;    // the first vertex of every layer
+  private final int[] last;     // the last vertex of every layer
+  private final IslandEvolution<double[]> engine;  // the evolutionary engine
+  private final MersenneTwisterRNG rng;            // the random no. generator
+
+  /**
+   * Constructs the GA instance.
+   * @param g the graph to label
+   * @param nIslands the number of islands to use
+   * @param pMut the mutation probability
+   * @param nCross number of crossover points
+   */
+  public RKGA(final Graph g, final int nIslands,
+              final double pMut, final int nCross) {
+    graph = g;
+    nVertices = g.getNVertices();
+    nLayers = g.getNLayers();
+    first = g.getLayerStarts();
+    last = g.getLayerEnds();
+    // Create the random number generator.
+    rng = new MersenneTwisterRNG();
+    // Create the chromosome factory.
+    CandidateFactory<double[]> factory = new ChromosomeFactory(nVertices);
+    // Create the operators.
+    List<EvolutionaryOperator<double[]>> operators =
+      new ArrayList<EvolutionaryOperator<double[]>>(nVertices);
+    operators.add(new DoubleArrayCrossover(nCross));
+    operators.add(new Mutator(new Probability(pMut)));
+    EvolutionaryOperator<double[]> pipeline =
+      new EvolutionPipeline<double[]>(operators);
+    // Set up the evolutionary engine.
+    engine = new IslandEvolution<>(nIslands,
+                                   new RingMigration(),
+                                   factory,
+                                   pipeline,
+                                   this,
+                                   new RouletteWheelSelection(),
+                                   rng);
+    // Add a logger.
+    engine.addEvolutionObserver(new Logger());
+  }
+
+  /**
+   * Converts a chromosome into a vector of vertex labels.
+   * @param chromosome the chromosome to decode
+   * @return the corresponding labels
+   */
+  private int[] decode(final double[] chromosome) {
+    ValueOrderedMap<Integer, Double> map = new ValueOrderedMap<>();
+    int[] labels = new int[nVertices];
+    // Proceed layer by layer.
+    for (int k = 0; k < nLayers; k++) {
+      map.clear();
+      // Map the indices in the layer to their chromosome values.
+      for (int i = first[k]; i <= last[k]; i++) {
+        map.put(i, chromosome[i]);
+      }
+      // Generate a list of vertices in the layer in chromosome order.
+      int[] sorted = map.ascendingKeyStream().mapToInt(i -> i).toArray();
+      // Label the vertices.
+      int j = 1;
+      for (int i = 0; i < sorted.length; i++) {
+        labels[sorted[i]] = j;
+        j += 1;
+      }
+    }
+    return labels;
+  }
+
+  /**
+   * Evaluates the fitness of a candidate.
+   * @param candidate the candidate solution
+   * @param population the population of solutions (ignored)
+   * @return the objective value of the candidate
+   */
+  @Override
+  public double getFitness(final double[] candidate,
+                           final List<? extends double[]> population) {
+    int[] labels = decode(candidate);
+    return graph.getObjValue(labels);
+  }
+
+  /**
+   * Checks whether the fitness function is "natural" (fitter candidates
+   * have higher scores) or not.
+   * @return false (lower fitness values are preferable)
+   */
+  @Override
+  public boolean isNatural() {
+    return false;
+  }
+
+  /**
+   * Seeds the random number generator.
+   * @param seed the new seed
+   */
+  public void setSeed(final long seed) {
+    rng.setSeed(seed);
+  }
+
+  /**
+   * Runs the evolutionary engine.
+   * @param popSize the population size per island to use
+   * @param elite the number of elite candidates to retain per island
+   * @param epoch the number of epochs between migrations
+   * @param migrations the number of migrants per island
+   * @param timeLimit time limit in seconds
+   * @param stagnation max generations without improvement
+   * @return the best candidate found
+   */
+  public int[] evolve(final int popSize, final int elite, final int epoch,
+                      final int migrations, final double timeLimit,
+                      final int stagnation) {
+    long tl = Math.round(SECTOMS * timeLimit);
+    // Run the evolutionary engine.
+    double[] sol =
+      engine.evolve(popSize, elite, epoch, migrations,
+                    new TargetFitness(0, false),       // stop on fitness = 0
+                    new ElapsedTime(tl),               // stop at time limit
+                    new Stagnation(stagnation, false)  // stop on stagnation
+      );
+    // Decode and return the best solution.
+    return decode(sol);
+  }
+
+  /**
+   * ChromosomeFactory is a factory for generating candidate solutions
+   * (chromosomes).
+   */
+  private static class ChromosomeFactory
+                 extends AbstractCandidateFactory<double[]>
+                 implements CandidateFactory<double[]> {
+
+    private final int length; // chromosome length
+
+    /**
+     * Constructs the factory.
+     * @param dimension the dimension of a chromosome vector (# of vertices)
+     */
+    ChromosomeFactory(final int dimension) {