test/loops/TSPExchangerTest.java

/*
 * Written by Doug Lea and Bill Scherer with assistance from members
 * of JCP JSR-166 Expert Group and released to the public domain, as
 * explained at http://creativecommons.org/publicdomain/zero/1.0/
 */

import java.util.*;
import java.util.concurrent.*;
import java.util.concurrent.atomic.*;
import java.util.concurrent.locks.*;

/**
 * A parallel Traveling Salesperson Problem (TSP) program based on a
 * genetic algorithm using an Exchanger.  A population of chromosomes is
 * distributed among "subpops".  Each chromosomes represents a tour,
 * and its fitness is the total tour length.
 *
 * A set of worker threads perform updates on subpops. The basic
 * update step is:
 * <ol>
 *   <li>Select a breeder b from the subpop
 *   <li>Create a strand of its tour with a random starting point and length
 *   <li>Offer the strand to the exchanger, receiving a strand from
 *       another subpop
 *   <li>Combine b and the received strand using crossing function to
 *       create new chromosome c.
 *   <li>Replace a chromosome in the subpop with c.
 * </ol>
 *
 * This continues for a given number of generations per subpop.
 * Because there are normally more subpops than threads, each worker
 * thread performs small (randomly sized) run of updates for one
 * subpop and then selects another. A run continues until there is at
 * most one remaining thread performing updates.
 *
 * See below for more details.
 */
public class TSPExchangerTest {
    static final int NCPUS = Runtime.getRuntime().availableProcessors();

    /** Runs start with two threads, increasing by two through max */
    static final int DEFAULT_MAX_THREADS = Math.max(4, NCPUS + NCPUS/2);

    /** The number of replication runs per thread value */
    static final int DEFAULT_REPLICATIONS = 3;

    /** If true, print statistics in SNAPSHOT_RATE intervals */
    static boolean verbose = true;
    static final long SNAPSHOT_RATE = 10000; // in milliseconds

    /**
     * The problem size. Each city is a random point. The goal is to
     * find a tour among them with smallest total Euclidean distance.
     */
    static final int DEFAULT_CITIES = 144;

    // Tuning parameters.

    /**
     * The number of chromosomes per subpop. Must be a power of two.
     *
     * Smaller values lead to faster iterations but poorer quality
     * results
     */
    static final int DEFAULT_SUBPOP_SIZE = 32;

    /**
     * The number of iterations per subpop. Convergence appears
     * to be roughly proportional to #cities-squared
     */
    static final int DEFAULT_GENERATIONS = DEFAULT_CITIES * DEFAULT_CITIES;

    /**
     * The number of subpops. The total population is #subpops * subpopSize,
     * which should be roughly on the order of #cities-squared
     *
     * Smaller values lead to faster total runs but poorer quality
     * results
     */
    static final int DEFAULT_NSUBPOPS = DEFAULT_GENERATIONS / DEFAULT_SUBPOP_SIZE;

    /**
     * The minimum length for a random chromosome strand.
     * Must be at least 1.
     */
    static final int MIN_STRAND_LENGTH = 3;

    /**
     * The probability mask value for creating random strands,
     * that have lengths at least MIN_STRAND_LENGTH, and grow
     * with exponential decay 2^(-(1/(RANDOM_STRAND_MASK + 1)
     * Must be 1 less than a power of two.
     */
    static final int RANDOM_STRAND_MASK = 7;

    /**
     * Probability control for selecting breeders.
     * Breeders are selected starting at the best-fitness chromosome,
     * with exponentially decaying probability
     * 1 / (subpopSize >>> BREEDER_DECAY).
     *
     * Larger values usually cause faster convergence but poorer
     * quality results
     */
    static final int BREEDER_DECAY = 1;

    /**
     * Probability control for selecting dyers.
     * Dyers are selected starting at the worst-fitness chromosome,
     * with exponentially decaying probability
     * 1 / (subpopSize >>> DYER_DECAY)
     *
     * Larger values usually cause faster convergence but poorer
     * quality results
     */
    static final int DYER_DECAY = 1;

    /**
     * The set of cities. Created once per program run, to
     * make it easier to compare solutions across different runs.
     */
    static CitySet cities;

    public static void main(String[] args) throws Exception {
        int maxThreads = DEFAULT_MAX_THREADS;
        int nCities = DEFAULT_CITIES;
        int subpopSize = DEFAULT_SUBPOP_SIZE;
        int nGen = nCities * nCities;
        int nSubpops = nCities * nCities / subpopSize;
        int nReps = DEFAULT_REPLICATIONS;

        try {
            int argc = 0;
            while (argc < args.length) {
                String option = args[argc++];
                if (option.equals("-c")) {
                    nCities = Integer.parseInt(args[argc]);
                    nGen = nCities * nCities;
                    nSubpops = nCities * nCities / subpopSize;
                }
                else if (option.equals("-p"))
                    subpopSize = Integer.parseInt(args[argc]);
                else if (option.equals("-g"))
                    nGen = Integer.parseInt(args[argc]);
                else if (option.equals("-n"))
                    nSubpops = Integer.parseInt(args[argc]);
                else if (option.equals("-q")) {
                    verbose = false;
                    argc--;
                }
                else if (option.equals("-r"))
                    nReps = Integer.parseInt(args[argc]);
                else
                    maxThreads = Integer.parseInt(option);
                argc++;
            }
        }
        catch (Exception e) {
            reportUsageErrorAndDie();
        }

        System.out.print("TSPExchangerTest");
        System.out.print(" -c " + nCities);
        System.out.print(" -g " + nGen);
        System.out.print(" -p " + subpopSize);
        System.out.print(" -n " + nSubpops);
        System.out.print(" -r " + nReps);
        System.out.print(" max threads " + maxThreads);
        System.out.println();

        cities = new CitySet(nCities);

        if (false && NCPUS > 4) {
            int h = NCPUS/2;
            System.out.printf("Threads: %4d Warmup\n", h);
            oneRun(h, nSubpops, subpopSize, nGen);
            Thread.sleep(500);
        }

        int maxt = (maxThreads < nSubpops) ? maxThreads : nSubpops;
        for (int j = 0; j < nReps; ++j) {
            for (int i = 2; i <= maxt; i += 2) {
                System.out.printf("Threads: %4d Replication: %2d\n", i, j);
                oneRun(i, nSubpops, subpopSize, nGen);
                Thread.sleep(500);
            }
        }
    }

    static void reportUsageErrorAndDie() {
        System.out.print("usage: TSPExchangerTest");
        System.out.print(" [-c #cities]");
        System.out.print(" [-p #subpopSize]");
        System.out.print(" [-g #generations]");
        System.out.print(" [-n #subpops]");
        System.out.print(" [-r #replications]");
        System.out.print(" [-q <quiet>]");
        System.out.print(" #threads]");
        System.out.println();
        System.exit(0);
    }

    /**
     * Performs one run with the given parameters.  Each run completes
     * when there are fewer than 2 active threads.  When there is
     * only one remaining thread, it will have no one to exchange
     * with, so it is terminated (via interrupt).
     */
    static void oneRun(int nThreads, int nSubpops, int subpopSize, int nGen)
        throws InterruptedException {
        Population p = new Population(nThreads, nSubpops, subpopSize, nGen);
        ProgressMonitor mon = null;
        if (verbose) {
            p.printSnapshot(0);
            mon = new ProgressMonitor(p);
            mon.start();
        }
        long startTime = System.nanoTime();
        p.start();
        p.awaitDone();
        long stopTime = System.nanoTime();
        if (mon != null)
            mon.interrupt();
        p.shutdown();
        //        Thread.sleep(100);

        long elapsed = stopTime - startTime;
        double secs = (double) elapsed / 1000000000.0;
        p.printSnapshot(secs);
    }

    /**
     * A Population creates the subpops, subpops, and threads for a run
     * and has control methods to start, stop, and report progress.
     */
    static final class Population {
        final Worker[] threads;
        final Subpop[] subpops;
        final Exchanger<Strand> exchanger;
        final CountDownLatch done;
        final int nGen;
        final int subpopSize;
        final int nThreads;

        Population(int nThreads, int nSubpops, int subpopSize, int nGen) {
            this.nThreads = nThreads;
            this.nGen = nGen;
            this.subpopSize = subpopSize;
            this.exchanger = new Exchanger<Strand>();
            this.done = new CountDownLatch(nThreads - 1);

            this.subpops = new Subpop[nSubpops];
            for (int i = 0; i < nSubpops; i++)
                subpops[i] = new Subpop(this);

            this.threads = new Worker[nThreads];
            int maxExchanges = nGen * nSubpops / nThreads;
            for (int i = 0; i < nThreads; ++i) {
                threads[i] = new Worker(this, maxExchanges);
            }

        }

        void start() {
            for (int i = 0; i < nThreads; ++i) {
                threads[i].start();
            }
        }

        /** Stop the tasks */
        void shutdown() {
            for (int i = 0; i < threads.length; ++ i)
                threads[i].interrupt();
        }

        void threadDone() {
            done.countDown();
        }

        /** Wait for tasks to complete */
        void awaitDone() throws InterruptedException {
            done.await();
        }

        int totalExchanges() {
            int xs = 0;
            for (int i = 0; i < threads.length; ++i)
                xs += threads[i].exchanges;
            return xs;
        }

        /**
         * Prints statistics, including best and worst tour lengths
         * for points scaled in [0,1), scaled by the square root of
         * number of points. This simplifies checking results.  The
         * expected optimal TSP for random points is believed to be
         * around 0.76 * sqrt(N). For papers discussing this, see
         * http://www.densis.fee.unicamp.br/~moscato/TSPBIB_home.html
         */
        void printSnapshot(double secs) {
            int xs = totalExchanges();
            long rate = (xs == 0) ? 0L : (long) ((secs * 1000000000.0) / xs);
            Chromosome bestc = subpops[0].chromosomes[0];
            Chromosome worstc = bestc;
            for (int k = 0; k < subpops.length; ++k) {
                Chromosome[] cs = subpops[k].chromosomes;
                if (cs[0].fitness < bestc.fitness)
                    bestc = cs[0];
                int w = cs[cs.length-1].fitness;
                if (cs[cs.length-1].fitness > worstc.fitness)
                    worstc = cs[cs.length-1];
            }
            double sqrtn = Math.sqrt(cities.length);
            double best = bestc.unitTourLength() / sqrtn;
            double worst = worstc.unitTourLength() / sqrtn;
            System.out.printf("N:%4d T:%8.3f B:%6.3f W:%6.3f X:%9d R:%7d\n",
                              nThreads, secs, best, worst, xs, rate);
            //            exchanger.printStats();
            //            System.out.print(" s: " + exchanger.aveSpins());
            //            System.out.print(" p: " + exchanger.aveParks());
        }
    }

    /**
     * Worker threads perform updates on subpops.
     */
    static final class Worker extends Thread {
        final Population pop;
        final int maxExchanges;
        int exchanges;
        final RNG rng = new RNG();

        Worker(Population pop, int maxExchanges) {
            this.pop = pop;
            this.maxExchanges = maxExchanges;
        }

        /**
         * Repeatedly, find a subpop that is not being updated by
         * another thread, and run a random number of updates on it.
         */
        public void run() {
            try {
                int len = pop.subpops.length;
                int pos = (rng.next() & 0x7FFFFFFF) % len;
                while (exchanges < maxExchanges) {
                    Subpop s = pop.subpops[pos];
                    AtomicBoolean busy = s.busy;
                    if (!busy.get() && busy.compareAndSet(false, true)) {
                        exchanges += s.runUpdates();
                        busy.set(false);
                        pos = (rng.next() & 0x7FFFFFFF) % len;
                    }
                    else if (++pos >= len)
                        pos = 0;
                }
                pop.threadDone();
            } catch (InterruptedException fallthrough) {
            }
        }
    }

    /**
     * A Subpop maintains a set of chromosomes.
     */
    static final class Subpop {
        /** The chromosomes, kept in sorted order */
        final Chromosome[] chromosomes;
        /** The parent population */
        final Population pop;
        /** Reservation bit for worker threads */
        final AtomicBoolean busy;
        /** The common exchanger, same for all subpops */
        final Exchanger<Strand> exchanger;
        /** The current strand being exchanged */
        Strand strand;
        /** Bitset used in cross */
        final int[] inTour;
        final RNG rng;
        final int subpopSize;

        Subpop(Population pop) {
            this.pop = pop;
            this.subpopSize = pop.subpopSize;
            this.exchanger = pop.exchanger;
            this.busy = new AtomicBoolean(false);
            this.rng = new RNG();
            int length = cities.length;
            this.strand = new Strand(length);
            this.inTour = new int[(length >>> 5) + 1];
            this.chromosomes = new Chromosome[subpopSize];
            for (int j = 0; j < subpopSize; ++j)
                chromosomes[j] = new Chromosome(length, rng);
            Arrays.sort(chromosomes);
        }

        /**
         * Run a random number of updates.  The number of updates is
         * at least 1 and no more than subpopSize.  This
         * controls the granularity of multiplexing subpop updates on
         * to threads. It is small enough to balance out updates
         * across tasks, but large enough to avoid having runs
         * dominated by subpop selection. It is randomized to avoid
         * long runs where pairs of subpops exchange only with each
         * other.  It is hardwired because small variations of it
         * don't matter much.
         *
         * @param g the first generation to run
         */
        int runUpdates() throws InterruptedException {
            int n = 1 + (rng.next() & ((subpopSize << 1) - 1));
            for (int i = 0; i < n; ++i)
                update();
            return n;
        }

        /**
         * Chooses a breeder, exchanges strand with another subpop, and
         * crosses them to create new chromosome to replace a chosen
         * dyer.
         */
        void update() throws InterruptedException {
            int b = chooseBreeder();
            int d = chooseDyer(b);
            Chromosome breeder = chromosomes[b];
            Chromosome child = chromosomes[d];
            chooseStrand(breeder);
            strand = exchanger.exchange(strand);
            cross(breeder, child);
            fixOrder(child, d);
        }

        /**
         * Chooses a breeder, with exponentially decreasing probability
         * starting at best.
         * @return index of selected breeder
         */
        int chooseBreeder() {
            int mask = (subpopSize >>> BREEDER_DECAY) - 1;
            int b = 0;
            while ((rng.next() & mask) != mask) {
                if (++b >= subpopSize)
                    b = 0;
            }
            return b;
        }

        /**
         * Chooses a chromosome that will be replaced, with
         * exponentially decreasing probability starting at
         * worst, ignoring the excluded index.
         * @param exclude index to ignore; use -1 to not exclude any
         * @return index of selected dyer
         */
        int chooseDyer(int exclude) {
            int mask = (subpopSize >>> DYER_DECAY)  - 1;
            int d = subpopSize - 1;
            while (d == exclude || (rng.next() & mask) != mask) {
                if (--d < 0)
                    d = subpopSize - 1;
            }
            return d;
        }

        /**
         * Select a random strand of b's.
         * @param breeder the breeder
         */
        void chooseStrand(Chromosome breeder) {
            int[] bs = breeder.alleles;
            int length = bs.length;
            int strandLength = MIN_STRAND_LENGTH;
            while (strandLength < length &&
                   (rng.next() & RANDOM_STRAND_MASK) != RANDOM_STRAND_MASK)
                strandLength++;
            strand.strandLength = strandLength;
            int[] ss = strand.alleles;
            int k = (rng.next() & 0x7FFFFFFF) % length;
            for (int i = 0; i < strandLength; ++i) {
                ss[i] = bs[k];
                if (++k >= length) k = 0;
            }
        }

        /**
         * Copies current strand to start of c's, and then appends all
         * remaining b's that aren't in the strand.
         * @param breeder the breeder
         * @param child the child
         */
        void cross(Chromosome breeder, Chromosome child) {
            for (int k = 0; k < inTour.length; ++k) // clear bitset
                inTour[k] = 0;

            // Copy current strand to c
            int[] cs = child.alleles;
            int ssize = strand.strandLength;
            int[] ss = strand.alleles;
            int i;
            for (i = 0; i < ssize; ++i) {
                int x = ss[i];
                cs[i] = x;
                inTour[x >>> 5] |= 1 << (x & 31); // record in bit set
            }

            // Find index of matching origin in b
            int first = cs[0];
            int j = 0;
            int[] bs = breeder.alleles;
            while (bs[j] != first)
                ++j;

            // Append remaining b's that aren't already in tour
            while (i < cs.length) {
                if (++j >= bs.length) j = 0;
                int x = bs[j];
                if ((inTour[x >>> 5] & (1 << (x & 31))) == 0)
                    cs[i++] = x;
            }

        }

        /**
         * Fixes the sort order of a changed Chromosome c at position k.
         * @param c the chromosome
         * @param k the index
         */
        void fixOrder(Chromosome c, int k) {
            Chromosome[] cs = chromosomes;
            int oldFitness = c.fitness;
            c.recalcFitness();
            int newFitness = c.fitness;
            if (newFitness < oldFitness) {
                int j = k;
                int p = j - 1;
                while (p >= 0 && cs[p].fitness > newFitness) {
                    cs[j] = cs[p];
                    j = p--;
                }
                cs[j] = c;
            } else if (newFitness > oldFitness) {
                int j = k;
                int n = j + 1;
                while (n < cs.length && cs[n].fitness < newFitness) {
                    cs[j] = cs[n];
                    j = n++;
                }
                cs[j] = c;
            }
        }
    }

    /**
     * A Chromosome is a candidate TSP tour.
     */
    static final class Chromosome implements Comparable {
        /** Index of cities in tour order */
        final int[] alleles;
        /** Total tour length */
        int fitness;

        /**
         * Initializes to random tour.
         */
        Chromosome(int length, RNG random) {
            alleles = new int[length];
            for (int i = 0; i < length; i++)
                alleles[i] = i;
            for (int i = length - 1; i > 0; i--) {
                int idx = (random.next() & 0x7FFFFFFF) % alleles.length;
                int tmp = alleles[i];
                alleles[i] = alleles[idx];
                alleles[idx] = tmp;
            }
            recalcFitness();
        }

        public int compareTo(Object x) { // to enable sorting
            int xf = ((Chromosome) x).fitness;
            int f = fitness;
            return ((f == xf) ? 0 :((f < xf) ? -1 : 1));
        }

        void recalcFitness() {
            int[] a = alleles;
            int len = a.length;
            int p = a[0];
            long f = cities.distanceBetween(a[len-1], p);
            for (int i = 1; i < len; i++) {
                int n = a[i];
                f += cities.distanceBetween(p, n);
                p = n;
            }
            fitness = (int) (f / len);
        }

        /**
         * Returns tour length for points scaled in [0, 1).
         */
        double unitTourLength() {
            int[] a = alleles;
            int len = a.length;
            int p = a[0];
            double f = cities.unitDistanceBetween(a[len-1], p);
            for (int i = 1; i < len; i++) {
                int n = a[i];
                f += cities.unitDistanceBetween(p, n);
                p = n;
            }
            return f;
        }

        /**
         * Checks that this tour visits each city.
         */
        void validate() {
            int len = alleles.length;
            boolean[] used = new boolean[len];
            for (int i = 0; i < len; ++i)
                used[alleles[i]] = true;
            for (int i = 0; i < len; ++i)
                if (!used[i])
                    throw new Error("Bad tour");
        }

    }

    /**
     * A Strand is a random sub-sequence of a Chromosome.  Each subpop
     * creates only one strand, and then trades it with others,
     * refilling it on each iteration.
     */
    static final class Strand {
        final int[] alleles;
        int strandLength;
        Strand(int length) { alleles = new int[length]; }
    }

    /**
     * A collection of (x,y) points that represent cities.
     */
    static final class CitySet {

        final int length;
        final int[] xPts;
        final int[] yPts;
        final int[][] distances;

        CitySet(int n) {
            this.length = n;
            this.xPts = new int[n];
            this.yPts = new int[n];
            this.distances = new int[n][n];

            RNG random = new RNG();
            for (int i = 0; i < n; i++) {
                xPts[i] = (random.next() & 0x7FFFFFFF);
                yPts[i] = (random.next() & 0x7FFFFFFF);
            }

            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    double dx = (double) xPts[i] - (double) xPts[j];
                    double dy = (double) yPts[i] - (double) yPts[j];
                    double dd = Math.hypot(dx, dy) / 2.0;
                    long ld = Math.round(dd);
                    distances[i][j] = (ld >= Integer.MAX_VALUE) ?
                        Integer.MAX_VALUE : (int) ld;
                }
            }
        }

        /**
         * Returns the cached distance between a pair of cities.
         */
        int distanceBetween(int i, int j) {
            return distances[i][j];
        }

        // Scale ints to doubles in [0,1)
        static final double PSCALE = (double) 0x80000000L;

        /**
         * Returns distance for points scaled in [0,1). This simplifies
         * checking results.  The expected optimal TSP for random
         * points is believed to be around 0.76 * sqrt(N). For papers
         * discussing this, see
         * http://www.densis.fee.unicamp.br/~moscato/TSPBIB_home.html
         */
        double unitDistanceBetween(int i, int j) {
            double dx = ((double) xPts[i] - (double) xPts[j]) / PSCALE;
            double dy = ((double) yPts[i] - (double) yPts[j]) / PSCALE;
            return Math.hypot(dx, dy);
        }

    }

    /**
     * Cheap XorShift random number generator
     */
    static final class RNG {
        /** Seed generator for XorShift RNGs */
        static final Random seedGenerator = new Random();

        int seed;
        RNG(int seed) { this.seed = seed; }
        RNG()         { this.seed = seedGenerator.nextInt() | 1; }

        int next() {
            int x = seed;
            x ^= x << 6;
            x ^= x >>> 21;
            x ^= x << 7;
            seed = x;
            return x;
        }
    }

    static final class ProgressMonitor extends Thread {
        final Population pop;
        ProgressMonitor(Population p) { pop = p; }
        public void run() {
            double time = 0;
            try {
                while (!Thread.interrupted()) {
                    sleep(SNAPSHOT_RATE);
                    time += SNAPSHOT_RATE;
                    pop.printSnapshot(time / 1000.0);
                }
            } catch (InterruptedException ie) {}
        }
    }
}
Revision:	1.19
Committed:	Sun Sep 13 16:28:14 2015 UTC (8 years, 7 months ago) by jsr166
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.18:	+7 -7 lines
Log Message:	consistent style for <li> tags, removing </li> end tags
#	Content
1	/*
2	* Written by Doug Lea and Bill Scherer with assistance from members
3	* of JCP JSR-166 Expert Group and released to the public domain, as
4	* explained at http://creativecommons.org/publicdomain/zero/1.0/
5	*/
6
7	import java.util.*;
8	import java.util.concurrent.*;
9	import java.util.concurrent.atomic.*;
10	import java.util.concurrent.locks.*;
11
12	/**
13	* A parallel Traveling Salesperson Problem (TSP) program based on a
14	* genetic algorithm using an Exchanger. A population of chromosomes is
15	* distributed among "subpops". Each chromosomes represents a tour,
16	* and its fitness is the total tour length.
17	*
18	* A set of worker threads perform updates on subpops. The basic
19	* update step is:
20	* <ol>
21	* <li>Select a breeder b from the subpop
22	* <li>Create a strand of its tour with a random starting point and length
23	* <li>Offer the strand to the exchanger, receiving a strand from
24	* another subpop
25	* <li>Combine b and the received strand using crossing function to
26	* create new chromosome c.
27	* <li>Replace a chromosome in the subpop with c.
28	* </ol>
29	*
30	* This continues for a given number of generations per subpop.
31	* Because there are normally more subpops than threads, each worker
32	* thread performs small (randomly sized) run of updates for one
33	* subpop and then selects another. A run continues until there is at
34	* most one remaining thread performing updates.
35	*
36	* See below for more details.
37	*/
38	public class TSPExchangerTest {
39	static final int NCPUS = Runtime.getRuntime().availableProcessors();
40
41	/** Runs start with two threads, increasing by two through max */
42	static final int DEFAULT_MAX_THREADS = Math.max(4, NCPUS + NCPUS/2);
43
44	/** The number of replication runs per thread value */
45	static final int DEFAULT_REPLICATIONS = 3;
46
47	/** If true, print statistics in SNAPSHOT_RATE intervals */
48	static boolean verbose = true;
49	static final long SNAPSHOT_RATE = 10000; // in milliseconds
50
51	/**
52	* The problem size. Each city is a random point. The goal is to
53	* find a tour among them with smallest total Euclidean distance.
54	*/
55	static final int DEFAULT_CITIES = 144;
56
57	// Tuning parameters.
58
59	/**
60	* The number of chromosomes per subpop. Must be a power of two.
61	*
62	* Smaller values lead to faster iterations but poorer quality
63	* results
64	*/
65	static final int DEFAULT_SUBPOP_SIZE = 32;
66
67	/**
68	* The number of iterations per subpop. Convergence appears
69	* to be roughly proportional to #cities-squared
70	*/
71	static final int DEFAULT_GENERATIONS = DEFAULT_CITIES * DEFAULT_CITIES;
72
73	/**
74	* The number of subpops. The total population is #subpops * subpopSize,
75	* which should be roughly on the order of #cities-squared
76	*
77	* Smaller values lead to faster total runs but poorer quality
78	* results
79	*/
80	static final int DEFAULT_NSUBPOPS = DEFAULT_GENERATIONS / DEFAULT_SUBPOP_SIZE;
81
82	/**
83	* The minimum length for a random chromosome strand.
84	* Must be at least 1.
85	*/
86	static final int MIN_STRAND_LENGTH = 3;
87
88	/**
89	* The probability mask value for creating random strands,
90	* that have lengths at least MIN_STRAND_LENGTH, and grow
91	* with exponential decay 2^(-(1/(RANDOM_STRAND_MASK + 1)
92	* Must be 1 less than a power of two.
93	*/
94	static final int RANDOM_STRAND_MASK = 7;
95
96	/**
97	* Probability control for selecting breeders.
98	* Breeders are selected starting at the best-fitness chromosome,
99	* with exponentially decaying probability
100	* 1 / (subpopSize >>> BREEDER_DECAY).
101	*
102	* Larger values usually cause faster convergence but poorer
103	* quality results
104	*/
105	static final int BREEDER_DECAY = 1;
106
107	/**
108	* Probability control for selecting dyers.
109	* Dyers are selected starting at the worst-fitness chromosome,
110	* with exponentially decaying probability
111	* 1 / (subpopSize >>> DYER_DECAY)
112	*
113	* Larger values usually cause faster convergence but poorer
114	* quality results
115	*/
116	static final int DYER_DECAY = 1;
117
118	/**
119	* The set of cities. Created once per program run, to
120	* make it easier to compare solutions across different runs.
121	*/
122	static CitySet cities;
123
124	public static void main(String[] args) throws Exception {
125	int maxThreads = DEFAULT_MAX_THREADS;
126	int nCities = DEFAULT_CITIES;
127	int subpopSize = DEFAULT_SUBPOP_SIZE;
128	int nGen = nCities * nCities;
129	int nSubpops = nCities * nCities / subpopSize;
130	int nReps = DEFAULT_REPLICATIONS;
131
132	try {
133	int argc = 0;
134	while (argc < args.length) {
135	String option = args[argc++];
136	if (option.equals("-c")) {
137	nCities = Integer.parseInt(args[argc]);
138	nGen = nCities * nCities;
139	nSubpops = nCities * nCities / subpopSize;
140	}
141	else if (option.equals("-p"))
142	subpopSize = Integer.parseInt(args[argc]);
143	else if (option.equals("-g"))
144	nGen = Integer.parseInt(args[argc]);
145	else if (option.equals("-n"))
146	nSubpops = Integer.parseInt(args[argc]);
147	else if (option.equals("-q")) {
148	verbose = false;
149	argc--;
150	}
151	else if (option.equals("-r"))
152	nReps = Integer.parseInt(args[argc]);
153	else
154	maxThreads = Integer.parseInt(option);
155	argc++;
156	}
157	}
158	catch (Exception e) {
159	reportUsageErrorAndDie();
160	}
161
162	System.out.print("TSPExchangerTest");
163	System.out.print(" -c " + nCities);
164	System.out.print(" -g " + nGen);
165	System.out.print(" -p " + subpopSize);
166	System.out.print(" -n " + nSubpops);
167	System.out.print(" -r " + nReps);
168	System.out.print(" max threads " + maxThreads);
169	System.out.println();
170
171	cities = new CitySet(nCities);
172
173	if (false && NCPUS > 4) {
174	int h = NCPUS/2;
175	System.out.printf("Threads: %4d Warmup\n", h);
176	oneRun(h, nSubpops, subpopSize, nGen);
177	Thread.sleep(500);
178	}
179
180	int maxt = (maxThreads < nSubpops) ? maxThreads : nSubpops;
181	for (int j = 0; j < nReps; ++j) {
182	for (int i = 2; i <= maxt; i += 2) {
183	System.out.printf("Threads: %4d Replication: %2d\n", i, j);
184	oneRun(i, nSubpops, subpopSize, nGen);
185	Thread.sleep(500);
186	}
187	}
188	}
189
190	static void reportUsageErrorAndDie() {
191	System.out.print("usage: TSPExchangerTest");
192	System.out.print(" [-c #cities]");
193	System.out.print(" [-p #subpopSize]");
194	System.out.print(" [-g #generations]");
195	System.out.print(" [-n #subpops]");
196	System.out.print(" [-r #replications]");
197	System.out.print(" [-q <quiet>]");
198	System.out.print(" #threads]");
199	System.out.println();
200	System.exit(0);
201	}
202
203	/**
204	* Performs one run with the given parameters. Each run completes
205	* when there are fewer than 2 active threads. When there is
206	* only one remaining thread, it will have no one to exchange
207	* with, so it is terminated (via interrupt).
208	*/
209	static void oneRun(int nThreads, int nSubpops, int subpopSize, int nGen)
210	throws InterruptedException {
211	Population p = new Population(nThreads, nSubpops, subpopSize, nGen);
212	ProgressMonitor mon = null;
213	if (verbose) {
214	p.printSnapshot(0);
215	mon = new ProgressMonitor(p);
216	mon.start();
217	}
218	long startTime = System.nanoTime();
219	p.start();
220	p.awaitDone();
221	long stopTime = System.nanoTime();
222	if (mon != null)
223	mon.interrupt();
224	p.shutdown();
225	// Thread.sleep(100);
226
227	long elapsed = stopTime - startTime;
228	double secs = (double) elapsed / 1000000000.0;
229	p.printSnapshot(secs);
230	}
231
232	/**
233	* A Population creates the subpops, subpops, and threads for a run
234	* and has control methods to start, stop, and report progress.
235	*/
236	static final class Population {
237	final Worker[] threads;
238	final Subpop[] subpops;
239	final Exchanger<Strand> exchanger;
240	final CountDownLatch done;
241	final int nGen;
242	final int subpopSize;
243	final int nThreads;
244
245	Population(int nThreads, int nSubpops, int subpopSize, int nGen) {
246	this.nThreads = nThreads;
247	this.nGen = nGen;
248	this.subpopSize = subpopSize;
249	this.exchanger = new Exchanger<Strand>();
250	this.done = new CountDownLatch(nThreads - 1);
251
252	this.subpops = new Subpop[nSubpops];
253	for (int i = 0; i < nSubpops; i++)
254	subpops[i] = new Subpop(this);
255
256	this.threads = new Worker[nThreads];
257	int maxExchanges = nGen * nSubpops / nThreads;
258	for (int i = 0; i < nThreads; ++i) {
259	threads[i] = new Worker(this, maxExchanges);
260	}
261
262	}
263
264	void start() {
265	for (int i = 0; i < nThreads; ++i) {
266	threads[i].start();
267	}
268	}
269
270	/** Stop the tasks */
271	void shutdown() {
272	for (int i = 0; i < threads.length; ++ i)
273	threads[i].interrupt();
274	}
275
276	void threadDone() {
277	done.countDown();
278	}
279
280	/** Wait for tasks to complete */
281	void awaitDone() throws InterruptedException {
282	done.await();
283	}
284
285	int totalExchanges() {
286	int xs = 0;
287	for (int i = 0; i < threads.length; ++i)
288	xs += threads[i].exchanges;
289	return xs;
290	}
291
292	/**
293	* Prints statistics, including best and worst tour lengths
294	* for points scaled in [0,1), scaled by the square root of
295	* number of points. This simplifies checking results. The
296	* expected optimal TSP for random points is believed to be
297	* around 0.76 * sqrt(N). For papers discussing this, see
298	* http://www.densis.fee.unicamp.br/~moscato/TSPBIB_home.html
299	*/
300	void printSnapshot(double secs) {
301	int xs = totalExchanges();
302	long rate = (xs == 0) ? 0L : (long) ((secs * 1000000000.0) / xs);
303	Chromosome bestc = subpops[0].chromosomes[0];
304	Chromosome worstc = bestc;
305	for (int k = 0; k < subpops.length; ++k) {
306	Chromosome[] cs = subpops[k].chromosomes;
307	if (cs[0].fitness < bestc.fitness)
308	bestc = cs[0];
309	int w = cs[cs.length-1].fitness;
310	if (cs[cs.length-1].fitness > worstc.fitness)
311	worstc = cs[cs.length-1];
312	}
313	double sqrtn = Math.sqrt(cities.length);
314	double best = bestc.unitTourLength() / sqrtn;
315	double worst = worstc.unitTourLength() / sqrtn;
316	System.out.printf("N:%4d T:%8.3f B:%6.3f W:%6.3f X:%9d R:%7d\n",
317	nThreads, secs, best, worst, xs, rate);
318	// exchanger.printStats();
319	// System.out.print(" s: " + exchanger.aveSpins());
320	// System.out.print(" p: " + exchanger.aveParks());
321	}
322	}
323
324	/**
325	* Worker threads perform updates on subpops.
326	*/
327	static final class Worker extends Thread {
328	final Population pop;
329	final int maxExchanges;
330	int exchanges;
331	final RNG rng = new RNG();
332
333	Worker(Population pop, int maxExchanges) {
334	this.pop = pop;
335	this.maxExchanges = maxExchanges;
336	}
337
338	/**
339	* Repeatedly, find a subpop that is not being updated by
340	* another thread, and run a random number of updates on it.
341	*/
342	public void run() {
343	try {
344	int len = pop.subpops.length;
345	int pos = (rng.next() & 0x7FFFFFFF) % len;
346	while (exchanges < maxExchanges) {
347	Subpop s = pop.subpops[pos];
348	AtomicBoolean busy = s.busy;
349	if (!busy.get() && busy.compareAndSet(false, true)) {
350	exchanges += s.runUpdates();
351	busy.set(false);
352	pos = (rng.next() & 0x7FFFFFFF) % len;
353	}
354	else if (++pos >= len)
355	pos = 0;
356	}
357	pop.threadDone();
358	} catch (InterruptedException fallthrough) {
359	}
360	}
361	}
362
363	/**
364	* A Subpop maintains a set of chromosomes.
365	*/
366	static final class Subpop {
367	/** The chromosomes, kept in sorted order */
368	final Chromosome[] chromosomes;
369	/** The parent population */
370	final Population pop;
371	/** Reservation bit for worker threads */
372	final AtomicBoolean busy;
373	/** The common exchanger, same for all subpops */
374	final Exchanger<Strand> exchanger;
375	/** The current strand being exchanged */
376	Strand strand;
377	/** Bitset used in cross */
378	final int[] inTour;
379	final RNG rng;
380	final int subpopSize;
381
382	Subpop(Population pop) {
383	this.pop = pop;
384	this.subpopSize = pop.subpopSize;
385	this.exchanger = pop.exchanger;
386	this.busy = new AtomicBoolean(false);
387	this.rng = new RNG();
388	int length = cities.length;
389	this.strand = new Strand(length);
390	this.inTour = new int[(length >>> 5) + 1];
391	this.chromosomes = new Chromosome[subpopSize];
392	for (int j = 0; j < subpopSize; ++j)
393	chromosomes[j] = new Chromosome(length, rng);
394	Arrays.sort(chromosomes);
395	}
396
397	/**
398	* Run a random number of updates. The number of updates is
399	* at least 1 and no more than subpopSize. This
400	* controls the granularity of multiplexing subpop updates on
401	* to threads. It is small enough to balance out updates
402	* across tasks, but large enough to avoid having runs
403	* dominated by subpop selection. It is randomized to avoid
404	* long runs where pairs of subpops exchange only with each
405	* other. It is hardwired because small variations of it
406	* don't matter much.
407	*
408	* @param g the first generation to run
409	*/
410	int runUpdates() throws InterruptedException {
411	int n = 1 + (rng.next() & ((subpopSize << 1) - 1));
412	for (int i = 0; i < n; ++i)
413	update();
414	return n;
415	}
416
417	/**
418	* Chooses a breeder, exchanges strand with another subpop, and
419	* crosses them to create new chromosome to replace a chosen
420	* dyer.
421	*/
422	void update() throws InterruptedException {
423	int b = chooseBreeder();
424	int d = chooseDyer(b);
425	Chromosome breeder = chromosomes[b];
426	Chromosome child = chromosomes[d];
427	chooseStrand(breeder);
428	strand = exchanger.exchange(strand);
429	cross(breeder, child);
430	fixOrder(child, d);
431	}
432
433	/**
434	* Chooses a breeder, with exponentially decreasing probability
435	* starting at best.
436	* @return index of selected breeder
437	*/
438	int chooseBreeder() {
439	int mask = (subpopSize >>> BREEDER_DECAY) - 1;
440	int b = 0;
441	while ((rng.next() & mask) != mask) {
442	if (++b >= subpopSize)
443	b = 0;
444	}
445	return b;
446	}
447
448	/**
449	* Chooses a chromosome that will be replaced, with
450	* exponentially decreasing probability starting at
451	* worst, ignoring the excluded index.
452	* @param exclude index to ignore; use -1 to not exclude any
453	* @return index of selected dyer
454	*/
455	int chooseDyer(int exclude) {
456	int mask = (subpopSize >>> DYER_DECAY) - 1;
457	int d = subpopSize - 1;
458	while (d == exclude \|\| (rng.next() & mask) != mask) {
459	if (--d < 0)
460	d = subpopSize - 1;
461	}
462	return d;
463	}
464
465	/**
466	* Select a random strand of b's.
467	* @param breeder the breeder
468	*/
469	void chooseStrand(Chromosome breeder) {
470	int[] bs = breeder.alleles;
471	int length = bs.length;
472	int strandLength = MIN_STRAND_LENGTH;
473	while (strandLength < length &&
474	(rng.next() & RANDOM_STRAND_MASK) != RANDOM_STRAND_MASK)
475	strandLength++;
476	strand.strandLength = strandLength;
477	int[] ss = strand.alleles;
478	int k = (rng.next() & 0x7FFFFFFF) % length;
479	for (int i = 0; i < strandLength; ++i) {
480	ss[i] = bs[k];
481	if (++k >= length) k = 0;
482	}
483	}
484
485	/**
486	* Copies current strand to start of c's, and then appends all
487	* remaining b's that aren't in the strand.
488	* @param breeder the breeder
489	* @param child the child
490	*/
491	void cross(Chromosome breeder, Chromosome child) {
492	for (int k = 0; k < inTour.length; ++k) // clear bitset
493	inTour[k] = 0;
494
495	// Copy current strand to c
496	int[] cs = child.alleles;
497	int ssize = strand.strandLength;
498	int[] ss = strand.alleles;
499	int i;
500	for (i = 0; i < ssize; ++i) {
501	int x = ss[i];
502	cs[i] = x;
503	inTour[x >>> 5] \|= 1 << (x & 31); // record in bit set
504	}
505
506	// Find index of matching origin in b
507	int first = cs[0];
508	int j = 0;
509	int[] bs = breeder.alleles;
510	while (bs[j] != first)
511	++j;
512
513	// Append remaining b's that aren't already in tour
514	while (i < cs.length) {
515	if (++j >= bs.length) j = 0;
516	int x = bs[j];
517	if ((inTour[x >>> 5] & (1 << (x & 31))) == 0)
518	cs[i++] = x;
519	}
520
521	}
522
523	/**
524	* Fixes the sort order of a changed Chromosome c at position k.
525	* @param c the chromosome
526	* @param k the index
527	*/
528	void fixOrder(Chromosome c, int k) {
529	Chromosome[] cs = chromosomes;
530	int oldFitness = c.fitness;
531	c.recalcFitness();
532	int newFitness = c.fitness;
533	if (newFitness < oldFitness) {
534	int j = k;
535	int p = j - 1;
536	while (p >= 0 && cs[p].fitness > newFitness) {
537	cs[j] = cs[p];
538	j = p--;
539	}
540	cs[j] = c;
541	} else if (newFitness > oldFitness) {
542	int j = k;
543	int n = j + 1;
544	while (n < cs.length && cs[n].fitness < newFitness) {
545	cs[j] = cs[n];
546	j = n++;
547	}
548	cs[j] = c;
549	}
550	}
551	}
552
553	/**
554	* A Chromosome is a candidate TSP tour.
555	*/
556	static final class Chromosome implements Comparable {
557	/** Index of cities in tour order */
558	final int[] alleles;
559	/** Total tour length */
560	int fitness;
561
562	/**
563	* Initializes to random tour.
564	*/
565	Chromosome(int length, RNG random) {
566	alleles = new int[length];
567	for (int i = 0; i < length; i++)
568	alleles[i] = i;
569	for (int i = length - 1; i > 0; i--) {
570	int idx = (random.next() & 0x7FFFFFFF) % alleles.length;
571	int tmp = alleles[i];
572	alleles[i] = alleles[idx];
573	alleles[idx] = tmp;
574	}
575	recalcFitness();
576	}
577
578	public int compareTo(Object x) { // to enable sorting
579	int xf = ((Chromosome) x).fitness;
580	int f = fitness;
581	return ((f == xf) ? 0 :((f < xf) ? -1 : 1));
582	}
583
584	void recalcFitness() {
585	int[] a = alleles;
586	int len = a.length;
587	int p = a[0];
588	long f = cities.distanceBetween(a[len-1], p);
589	for (int i = 1; i < len; i++) {
590	int n = a[i];
591	f += cities.distanceBetween(p, n);
592	p = n;
593	}
594	fitness = (int) (f / len);
595	}
596
597	/**
598	* Returns tour length for points scaled in [0, 1).
599	*/
600	double unitTourLength() {
601	int[] a = alleles;
602	int len = a.length;
603	int p = a[0];
604	double f = cities.unitDistanceBetween(a[len-1], p);
605	for (int i = 1; i < len; i++) {
606	int n = a[i];
607	f += cities.unitDistanceBetween(p, n);
608	p = n;
609	}
610	return f;
611	}
612
613	/**
614	* Checks that this tour visits each city.
615	*/
616	void validate() {
617	int len = alleles.length;
618	boolean[] used = new boolean[len];
619	for (int i = 0; i < len; ++i)
620	used[alleles[i]] = true;
621	for (int i = 0; i < len; ++i)
622	if (!used[i])
623	throw new Error("Bad tour");
624	}
625
626	}
627
628	/**
629	* A Strand is a random sub-sequence of a Chromosome. Each subpop
630	* creates only one strand, and then trades it with others,
631	* refilling it on each iteration.
632	*/
633	static final class Strand {
634	final int[] alleles;
635	int strandLength;
636	Strand(int length) { alleles = new int[length]; }
637	}
638
639	/**
640	* A collection of (x,y) points that represent cities.
641	*/
642	static final class CitySet {
643
644	final int length;
645	final int[] xPts;
646	final int[] yPts;
647	final int[][] distances;
648
649	CitySet(int n) {
650	this.length = n;
651	this.xPts = new int[n];
652	this.yPts = new int[n];
653	this.distances = new int[n][n];
654
655	RNG random = new RNG();
656	for (int i = 0; i < n; i++) {
657	xPts[i] = (random.next() & 0x7FFFFFFF);
658	yPts[i] = (random.next() & 0x7FFFFFFF);
659	}
660
661	for (int i = 0; i < n; i++) {
662	for (int j = 0; j < n; j++) {
663	double dx = (double) xPts[i] - (double) xPts[j];
664	double dy = (double) yPts[i] - (double) yPts[j];
665	double dd = Math.hypot(dx, dy) / 2.0;
666	long ld = Math.round(dd);
667	distances[i][j] = (ld >= Integer.MAX_VALUE) ?
668	Integer.MAX_VALUE : (int) ld;
669	}
670	}
671	}
672
673	/**
674	* Returns the cached distance between a pair of cities.
675	*/
676	int distanceBetween(int i, int j) {
677	return distances[i][j];
678	}
679
680	// Scale ints to doubles in [0,1)
681	static final double PSCALE = (double) 0x80000000L;
682
683	/**
684	* Returns distance for points scaled in [0,1). This simplifies
685	* checking results. The expected optimal TSP for random
686	* points is believed to be around 0.76 * sqrt(N). For papers
687	* discussing this, see
688	* http://www.densis.fee.unicamp.br/~moscato/TSPBIB_home.html
689	*/
690	double unitDistanceBetween(int i, int j) {
691	double dx = ((double) xPts[i] - (double) xPts[j]) / PSCALE;
692	double dy = ((double) yPts[i] - (double) yPts[j]) / PSCALE;
693	return Math.hypot(dx, dy);
694	}
695
696	}
697
698	/**
699	* Cheap XorShift random number generator
700	*/
701	static final class RNG {
702	/** Seed generator for XorShift RNGs */
703	static final Random seedGenerator = new Random();
704
705	int seed;
706	RNG(int seed) { this.seed = seed; }
707	RNG() { this.seed = seedGenerator.nextInt() \| 1; }
708
709	int next() {
710	int x = seed;
711	x ^= x << 6;
712	x ^= x >>> 21;
713	x ^= x << 7;
714	seed = x;
715	return x;
716	}
717	}
718
719	static final class ProgressMonitor extends Thread {
720	final Population pop;
721	ProgressMonitor(Population p) { pop = p; }
722	public void run() {
723	double time = 0;
724	try {
725	while (!Thread.interrupted()) {
726	sleep(SNAPSHOT_RATE);
727	time += SNAPSHOT_RATE;
728	pop.printSnapshot(time / 1000.0);
729	}
730	} catch (InterruptedException ie) {}
731	}
732	}
733	}