update readme and add Zipf exponent option to CLI

.idea/.name

@@ -0,0 +1 @@
+BiGpairSEQ

.idea/libraries/commons_rng_1.xml

@@ -0,0 +1,13 @@
+<component name="libraryTable">
+  <library name="commons-rng-1">
+    <CLASSES>
+      <root url="file://$USER_HOME$/Downloads/commons-rng-1.6" />
+    </CLASSES>
+    <JAVADOC />
+    <SOURCES>
+      <root url="file://$USER_HOME$/Downloads/commons-rng-1.6" />
+    </SOURCES>
+    <jarDirectory url="file://$USER_HOME$/Downloads/commons-rng-1.6" recursive="false" />
+    <jarDirectory url="file://$USER_HOME$/Downloads/commons-rng-1.6" recursive="false" type="SOURCES" />
+  </library>
+</component>

readme.md

@@ -1,30 +1,34 @@
 # BiGpairSEQ SIMULATOR
 
 ## CONTENTS
-1. ABOUT
-2. THEORY
-3. THE BiGpairSEQ ALGORITHM
-4. USAGE
-   1. RUNNING THE PROGRAM
-   2. COMMAND LINE OPTIONS
-   3. INTERACTIVE INTERFACE
-   4. INPUT/OUTPUT
-      1. Cell Sample Files
-      2. Sample Plate Files
-      3. Graph/Data Files
-      4. Matching Results Files
-5. RESULTS
-   1. SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL
-   2. SIMULATING EXPERIMENTS FROM pairSEQ PAPER
-6. TODO
-7. CITATIONS
-8. ACKNOWLEDGEMENTS
-9. AUTHOR
+1. [ABOUT](#about)
+2. [THEORY](#theory)
+3. [THE BiGpairSEQ ALGORITHM](#the-bigpairseq-algorithm)
+4. [USAGE](#usage)
+   1. [RUNNING THE PROGRAM](#running-the-program)
+   2. [COMMAND LINE OPTIONS](#command-line-options)
+   3. [INTERACTIVE INTERFACE](#interactive-interface)
+   4. [INPUT/OUTPUT](#input-output)
+      1. [Cell Sample Files](#cell-sample-files)
+      2. [Sample Plate Files](#sample-plate-files)
+      3. [Graph/Data Files](#graph-data-files)
+      4. [Matching Results Files](#matching-results-files)
+5. [RESULTS](#results)
+   1. [SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL](#sample-plates-with-varying-numbers-of-cells-per-well)
+   2. [SIMULATING EXPERIMENTS FROM THE 2015 pairSEQ PAPER](#simulating-experiments-from-the-2015-pairseq-paper)
+      1. [EXPERIMENT 1](#experiment-1)
+      2. [EXPERIMENT 3](#experiment-3)
+6. [CITATIONS](#citations)
+7. [EXTERNAL LIBRARIES USED](#external-libraries-used)
+8. [ACKNOWLEDGEMENTS](#acknowledgements)
+9. [AUTHOR](#author)
+10. [DISCLOSURE](#disclosure)
+11. [TODO](#todo)
 
 ## ABOUT
 
 This program simulates BiGpairSEQ (Bipartite Graph pairSEQ), a graph theory-based adaptation
-of the pairSEQ algorithm (Howie, et al. 2015) for pairing T cell receptor sequences.
+of the pairSEQ algorithm ([Howie, et al. 2015](#citations)) for pairing T cell receptor sequences.
 
 ## THEORY
 
@@ -50,21 +54,17 @@ matching (MWM) on a bipartite graph--the subset of vertex-disjoint edges whose w
 
 This is a well-studied combinatorial optimization problem, with many known algorithms that produce 
 provably-optimal solutions. The most theoretically efficient algorithm known to the author for maximum weight matching of a bipartite
-graph with strictly integral weights is from Duan and Su (2012). For a graph with m edges, n vertices per side, 
+graph with strictly integral weights is from [Duan and Su (2012)](#citations). For a graph with m edges, n vertices per side, 
 and maximum integer edge weight N, their algorithm runs in **O(m sqrt(n) log(N))** time. As the graph representation of 
-a pairSEQ experiment is bipartite with integer weights, this algorithm seems ideal for BiGpairSEQ. Unfortunately, it's a 
-fairly new algorithm, and not yet implemented by the graph theory library used in this simulator (JGraphT), nor has the author had 
-time to implement it himself.
-
-There have been some studies which show that [auction algorithms](https://en.wikipedia.org/wiki/Auction_algorithm) for the assignment problem can have superior performance in 
-real-world implementations, due to their simplicity, than more complex algorithms with better theoretical asymptotic 
-performance. But, again, there is no such algorithms implemented by JGraphT, nor has the author yet had time to implement one.
+a pairSEQ experiment is bipartite with integer weights, this algorithm seems ideal for BiGpairSEQ. Unfortunately, it is not 
+implemented by the graph theory library used in this simulator (JGraphT), and the author has not yet had time to write a
+full, optimized implementation himself for testing.
 
 So this program instead uses the [Fibonacci heap](https://en.wikipedia.org/wiki/Fibonacci_heap) based algorithm of Fredman and Tarjan (1987) (essentially
-[the Hungarian algorithm](https://en.wikipedia.org/wiki/Hungarian_algorithm) augmented with a more efficeint priority queue) which has a worst-case
-runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in Melhorn and Näher (1999). (The simulator
-allows the substitution of a [pairing heap](https://en.wikipedia.org/wiki/Pairing_heap) for a Fibonacci heap, though the relative performance difference of the two
-has not yet been thoroughly tested.)
+[the Hungarian algorithm](https://en.wikipedia.org/wiki/Hungarian_algorithm) augmented with a more efficient priority queue) which has a worst-case
+runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in [Melhorn and Näher (1999)](#citations). (The simulator can use either a
+Fibonacci heap or a [pairing heap](https://en.wikipedia.org/wiki/Pairing_heap) as desired. By default, a pairing heap is used,
+as in practice they often offer superior performance.)
 
 One possible advantage of this less efficient algorithm is that the Hungarian algorithm and its variations work with both the balanced and the unbalanced assignment problem
 (that is, cases where both sides of the bipartite graph have the same number of vertices and those in which they don't.)
@@ -74,6 +74,15 @@ be balanced assignment problems, in practice sequence dropout can cause them to
 the Hungarian algorithm, graph doubling--which could be challenging with the computational resources available to the 
 author--has not yet been necessary.
 
+There have been some studies which show that [auction algorithms](https://en.wikipedia.org/wiki/Auction_algorithm) for the assignment problem can have superior performance in
+real-world implementations, due to their simplicity, than more complex algorithms with better theoretical asymptotic
+performance. The author has implemented a basic forward auction algorithm, which produces optimal assignment for unbalanced bipartite graphs with 
+integer weights. To allow for unbalanced assignment, this algorithm eschews epsilon-scaling, 
+and as a result is prone to "bidding-wars" which increase run time, making it less efficient than the implementation of 
+the Fredman-Tarjan algorithm in JGraphT. A forward/reverse auction algorithm as developed by Bertsekas and Castañon
+should be able to handle unbalanced (or, as they call it, asymmetric) assignment much more efficiently, but has yet to be 
+implemented.
+
 The relative time/space efficiencies of BiGpairSEQ when backed by different MWM algorithms remains an open problem.
 
 ## THE BiGpairSEQ ALGORITHM
@@ -127,7 +136,7 @@ There are a number of command line options, to allow the program to be used in s
 `java -jar BiGpairSEQ_Sim.jar -help`
 
 ```
-usage: BiGpairSEQ_Sim.jar
+ usage: BiGpairSEQ_Sim.jar
  -cells,--make-cells   Makes a cell sample file of distinct T cells
  -graph,--make-graph   Makes a graph/data file. Requires a cell sample
                        file and a sample plate file
@@ -147,6 +156,8 @@ usage: BiGpairSEQ_Sim.jar -plate
  -c,--cell-file <filename>     The cell sample file to use
  -d,--dropout-rate <rate>      The sequence dropout rate due to
                                amplification error. (0.0 - 1.0)
+ -exp <value>                  If using -zipf flag, exponent value for
+                               distribution
  -exponential                  Use an exponential distribution for cell
                                sample
  -gaussian                     Use a Gaussian distribution for cell sample
@@ -164,6 +175,7 @@ usage: BiGpairSEQ_Sim.jar -plate
  -stddev <value>               If using -gaussian flag, standard deviation
                                for distrbution
  -w,--wells <number>           The number of wells on the sample plate
+ -zipf                         Use a Zipf distribution for cell sample
 
 usage: BiGpairSEQ_Sim.jar -graph
  -c,--cell-file <filename>                Cell sample file to use for
@@ -225,7 +237,6 @@ usage: BiGpairSEQ_Sim.jar -match
                                to stdout.
  -pv,--p-value                 (Optional) Calculate p-values for sequence
                                pairs.
-
 ```
 
 ### INTERACTIVE INTERFACE
@@ -331,6 +342,8 @@ Options when making a Sample Plate file:
     * Standard deviation size 
   * Exponential
     * Lambda value
+  * Zipf
+    * Exponent value
 * Total number of wells on the plate
 * Well populations random or fixed
   * If random, minimum and maximum population sizes
@@ -474,8 +487,9 @@ Several BiGpairSEQ simulations were performed on a home computer with the follow
 * Linux Mint 21 (5.15 kernel)
 
 ### SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL
-NOTE: these results were obtained with an earlier version of BiGpairSEQ_Sim, and should be re-run with the current version.
-The observed behavior is not believed to be likely to change, however.
+
+The probability calculations used by pairSEQ require that every well on the sample plate contain the same number of T cells.
+BiGpairSEQ does not share this limitation; it is robust to variations in the number of cells per well.
 
 A series of BiGpairSEQ simulations were conducted using a cell sample file of 3.5 million unique T cells. From these cells,
 10 sample plate files were created. All of these sample plates had 96 wells, used an exponential distribution with a lambda of 0.6, and
@@ -492,6 +506,9 @@ The well populations of the plates were:
 All BiGpairSEQ simulations were run with a low overlap threshold of 3 and a high overlap threshold of 94.
 No optional filters were used, so pairing was attempted for all sequences with overlaps within the threshold values.
 
+NOTE: these results were obtained with an earlier version of BiGpairSEQ_Sim, and should be re-run with the current version.
+The observed behavior is not believed to be likely to change, however.
+
 Constant well population plate results:
 
 | |1000 Cell/Well Plate|2000 Cell/Well Plate|3000 Cell/Well Plate|4000 Cell/Well Plate|5000 Cell/Well Plate
@@ -584,63 +601,13 @@ underlying frequency distribution drastically affect the results. The real distr
 than the simulated exponential distribution. Implementing a way to exert finer control over the sampling distribution from 
 the file of distinct cells may enable better simulated replication of this experiment.
   
-## TODO
-
-* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
-* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
-  * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
-  * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
-  * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific.
-* ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
-  * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
-  * ~~Add controllable heap-type parameter?~~
-    * Parameter implemented. Fibonacci heap the current default.
-* ~~Implement sample plates with random numbers of T cells per well.~~ DONE
-  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
-    * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
-* ~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows.~~
-  * ~~Problem is variable number of cells in a well~~
-  * ~~Apache Commons CSV library writes entries a row at a time~~ 
-    * Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
-* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
-  * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
-* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
-* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
-  * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
-  * This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
-* ~~Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.~~ DONE
-  * Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
-* ~~Add read depth simulation options to CLI~~ DONE
-* ~~Update graphml output to reflect current Vertex class attributes~~ DONE
-  * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
-* ~~Implement simulation of sequences being misread as other real sequence~~ DONE
-* Update matching metadata output options in CLI
-* Add frequency distribution details to metadata output
-  * need to make an enum for the different distribution types and refactor the Plate class and user interfaces, also add the necessary fields to GraphWithMapData and then call if from Simulator
-* Update performance data in this readme
-* Add section to ReadMe describing data filtering methods.
-* Re-implement CDR1 matching method
-* Refactor simulator code to collect all needed data in a single scan of the plate
-  * Currently it scans once for the vertices and then again for the edge weights. This made simulating read depth awkward, and incompatible with caching of plate files.
-  * This would be a fairly major rewrite of the simulator code, but could make things faster, and would definitely make them cleaner.
-* Implement Duan and Su's maximum weight matching algorithm
-    * Add controllable algorithm-type parameter?
-    * This would be fun and valuable, but probably take more time than I have for a hobby project.
-* Implement an auction algorithm for maximum weight matching
-* Implement an algorithm for approximating a maximum weight matching
-  * Some of these run in linear or near-linear time
-  * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
-* Implement Vose's alias method for arbitrary statistical distributions of cells
-  * Should probably refactor to use apache commons rng for this
-* Use commons JCS for caching
-* Parameterize pre-filtering options
-  
-  
 ## CITATIONS
 * Howie, B., Sherwood, A. M., et al. ["High-throughput pairing of T cell receptor alpha and beta sequences."](https://pubmed.ncbi.nlm.nih.gov/26290413/) Sci. Transl. Med. 7, 301ra131 (2015)
 * Duan, R., Su H. ["A Scaling Algorithm for Maximum Weight Matching in Bipartite Graphs."](https://web.eecs.umich.edu/~pettie/matching/Duan-Su-scaling-bipartite-matching.pdf) Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms, p. 1413-1424. (2012)
 * Melhorn, K., Näher, St. [The LEDA Platform of Combinatorial and Geometric Computing.](https://people.mpi-inf.mpg.de/~mehlhorn/LEDAbook.html) Cambridge University Press. Chapter 7, Graph Algorithms; p. 132-162 (1999)
 * Fredman, M., Tarjan, R. ["Fibonacci heaps and their uses in improved network optimization algorithms."](https://www.cl.cam.ac.uk/teaching/1011/AlgorithII/1987-FredmanTar-fibonacci.pdf) J. ACM, 34(3):596–615 (1987))
+* Bertsekas, D., Castañon, D. ["A forward/reverse auction algorithm for asymmetric assignment problems."](https://www.mit.edu/~dimitrib/For_Rev_Asym_Auction.pdf) Computational Optimization and Applications 1, 277-297 (1992)
+* Dimitrios Michail, Joris Kinable, Barak Naveh, and John V. Sichi. ["JGraphT—A Java Library for Graph Data Structures and Algorithms."](https://dl.acm.org/doi/10.1145/3381449) ACM Trans. Math. Softw. 46, 2, Article 16 (2020)
 
 ## EXTERNAL LIBRARIES USED
 * [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
@@ -649,8 +616,72 @@ the file of distinct cells may enable better simulated replication of this exper
 * [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
 
 ## ACKNOWLEDGEMENTS
-BiGpairSEQ was conceived in collaboration with Dr. Alice MacQueen, who brought the original 
+BiGpairSEQ was conceived in collaboration with the author's spouse, Dr. Alice MacQueen, who brought the original 
 pairSEQ paper to the author's attention and explained all the biology terms he didn't know.
 
 ## AUTHOR
-BiGpairSEQ algorithm and simulation by Eugene Fischer, 2021. Improvements and documentation, 2022.
+BiGpairSEQ algorithm and simulation by Eugene Fischer, 2021. Improvements and documentation, 2022–2025.
+
+## DISCLOSURE
+The earliest versions of the BiGpairSEQ simulator were written in 2021 to let Dr. MacQueen test hypothetical extensions 
+of the published pairSEQ protocol while she was interviewing for a position at Adaptive Biotechnologies. She was 
+employed at Adaptive Biotechnologies starting in 2022.
+
+The author has worked on this BiGpairSEQ simulator since 2021 without Dr. MacQueen's involvement, since she has had
+access to related, proprietary technologies. The author has had no such access, relying exclusively on the 2015 pairSEQ
+paper and other academic publications. He continues to work on the BiGpairSEQ simulator recreationally, as it has been  
+a means of exploring some very beautiful math.
+
+## TODO
+
+* Update CLI option text in this readme to include Zipf distribution options
+* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
+* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
+    * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
+    * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
+    * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific.
+* ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
+    * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
+    * ~~Add controllable heap-type parameter?~~
+        * Parameter implemented. Pairing heap the current default.
+* ~~Implement sample plates with random numbers of T cells per well.~~ DONE
+    * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
+        * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
+* ~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows.~~
+    * ~~Problem is variable number of cells in a well~~
+    * ~~Apache Commons CSV library writes entries a row at a time~~
+        * Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
+* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
+    * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
+* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
+* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
+    * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
+    * This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
+* ~~Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.~~ DONE
+    * Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
+* ~~Add read depth simulation options to CLI~~ DONE
+* ~~Update graphml output to reflect current Vertex class attributes~~ DONE
+    * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
+* ~~Implement simulation of sequences being misread as other real sequence~~ DONE
+* Implement redistributive heap for LEDA matching algorithm to achieve theoretical worst case of O(n(m + n log C)) where C is highest edge weight.
+* Update matching metadata output options in CLI
+* Add frequency distribution details to metadata output
+    * need to make an enum for the different distribution types and refactor the Plate class and user interfaces, also add the necessary fields to GraphWithMapData and then call if from Simulator
+* Update performance data in this readme
+* ~~Add section to ReadMe describing data filtering methods.~~ DONE, now part of algorithm description
+* Re-implement CDR1 matching method
+* ~~Refactor simulator code to collect all needed data in a single scan of the plate~~ DONE
+    * ~~Currently it scans once for the vertices and then again for the edge weights. This made simulating read depth awkward, and incompatible with caching of plate files.~~
+    * ~~This would be a fairly major rewrite of the simulator code, but could make things faster, and would definitely make them cleaner.~~
+* Implement Duan and Su's maximum weight matching algorithm
+    * ~~Add controllable algorithm-type parameter?~~ DONE
+    * This would be fun and valuable, but probably take more time than I have for a hobby project.
+* ~~Implement an auction algorithm for maximum weight matching~~ DONE
+* Implement a forward/reverse auction algorithm for maximum weight matching
+* Implement an algorithm for approximating a maximum weight matching
+    * Some of these run in linear or near-linear time
+    * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
+* Implement Vose's alias method for arbitrary statistical distributions of cells
+    * Should probably refactor to use apache commons rng for this
+* Use commons JCS for caching
+* Parameterize pre-filtering options

src/main/java/AlgorithmType.java

@@ -0,0 +1,5 @@
+public enum AlgorithmType {
+    HUNGARIAN, //Hungarian algorithm
+    AUCTION, //Forward auction algorithm
+    INTEGER_WEIGHT_SCALING, //integer weight scaling algorithm of Duan and Su
+}

src/main/java/BiGpairSEQ.java

@@ -13,7 +13,9 @@ public class BiGpairSEQ {
     private static boolean cacheCells = false;
     private static boolean cachePlate = false;
     private static boolean cacheGraph = false;
-    private static HeapType priorityQueueHeapType = HeapType.FIBONACCI;
+    private static AlgorithmType matchingAlgoritmType = AlgorithmType.HUNGARIAN;
+    private static HeapType priorityQueueHeapType = HeapType.PAIRING;
+    private static DistributionType distributionType = DistributionType.ZIPF;
     private static boolean outputBinary = true;
     private static boolean outputGraphML = false;
     private static boolean calculatePValue = false;
@@ -59,6 +61,10 @@ public class BiGpairSEQ {
         return cellFilename;
     }
 
+    public static DistributionType getDistributionType() {return distributionType;}
+
+    public static void setDistributionType(DistributionType type) {distributionType = type;}
+
     public static Plate getPlateInMemory() {
         return plateInMemory;
     }
@@ -108,7 +114,6 @@ public class BiGpairSEQ {
         return graphFilename;
     }
 
-
     public static boolean cacheCells() {
         return cacheCells;
     }
@@ -157,10 +162,18 @@ public class BiGpairSEQ {
         BiGpairSEQ.cacheGraph = cacheGraph;
     }
 
-    public static String getPriorityQueueHeapType() {
-        return priorityQueueHeapType.name();
+    public static HeapType getPriorityQueueHeapType() {
+        return priorityQueueHeapType;
     }
 
+    public static AlgorithmType getMatchingAlgoritmType() { return matchingAlgoritmType; }
+
+    public static void setHungarianAlgorithm() { matchingAlgoritmType = AlgorithmType.HUNGARIAN; }
+
+    public static void setIntegerWeightScalingAlgorithm() { matchingAlgoritmType = AlgorithmType.INTEGER_WEIGHT_SCALING; }
+
+    public static void setAuctionAlgorithm() { matchingAlgoritmType = AlgorithmType.AUCTION; }
+
     public static void setPairingHeap() {
         priorityQueueHeapType = HeapType.PAIRING;
     }

src/main/java/CommandLineInterface.java

@@ -123,16 +123,20 @@ public class CommandLineInterface {
                 Plate plate;
                 if (line.hasOption("poisson")) {
                     Double stdDev = Math.sqrt(numWells);
-                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev);
                 }
                 else if (line.hasOption("gaussian")) {
                     Double stdDev = Double.parseDouble(line.getOptionValue("stddev"));
-                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev);
+                }
+                else if (line.hasOption("zipf")) {
+                    Double zipfExponent = Double.parseDouble(line.getOptionValue("exp"));
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, zipfExponent);
                 }
                 else {
                     assert line.hasOption("exponential");
                     Double lambda = Double.parseDouble(line.getOptionValue("lambda"));
-                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, lambda, true);
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, lambda);
                 }
                 PlateFileWriter writer = new PlateFileWriter(outputFilename, plate);
                 writer.writePlateFile();
@@ -340,9 +344,13 @@ public class CommandLineInterface {
         Option exponential = Option.builder("exponential")
                 .desc("Use an exponential distribution for cell sample")
                 .build();
+        Option zipf = Option.builder("zipf")
+                .desc("Use a Zipf distribution for cell sample")
+                .build();
         distributions.addOption(poisson);
         distributions.addOption(gaussian);
         distributions.addOption(exponential);
+        distributions.addOption(zipf);
         //options group for statistical distribution parameters
         OptionGroup statParams = new OptionGroup();// add this to plate options
         Option stdDev = Option.builder("stddev")
@@ -355,6 +363,11 @@ public class CommandLineInterface {
                 .hasArg()
                 .argName("value")
                 .build();
+        Option zipfExponent = Option.builder("exp")
+                .desc("If using -zipf flag, exponent value for distribution")
+                .hasArg()
+                .argName("value")
+                .build();
         statParams.addOption(stdDev);
         statParams.addOption(lambda);
         //Option group for random plate or set populations
@@ -386,6 +399,7 @@ public class CommandLineInterface {
         plateOptions.addOptionGroup(statParams);
         plateOptions.addOptionGroup(wellPopOptions);
         plateOptions.addOption(dropoutRate);
+        plateOptions.addOption(zipfExponent);
         plateOptions.addOption(outputFileOption());
         return plateOptions;
     }

src/main/java/DistributionType.java

@@ -0,0 +1,6 @@
+public enum DistributionType {
+    POISSON,
+    GAUSSIAN,
+    EXPONENTIAL,
+    ZIPF
+}

src/main/java/InteractiveInterface.java

@@ -89,14 +89,12 @@ public class InteractiveInterface {
     private static void makePlate() {
         String cellFile = null;
         String filename = null;
-        Double stdDev = 0.0;
+        Double parameter = 0.0;
         Integer numWells = 0;
         Integer numSections;
         Integer[] populations = {1};
         Double dropOutRate = 0.0;
-        boolean poisson = false;
-        boolean exponential = false;
-        double lambda = 1.5;
+;
         try {
             System.out.println("\nSimulated sample plates consist of:");
             System.out.println("* a number of wells");
@@ -114,33 +112,46 @@ public class InteractiveInterface {
             System.out.println("1) Poisson");
             System.out.println("2) Gaussian");
             System.out.println("3) Exponential");
-//            System.out.println("(Note: approximate distribution in original paper is exponential, lambda = 0.6)");
-//            System.out.println("(lambda value approximated from slope of log-log graph in figure 4c)");
+            System.out.println("4) Zipf");
+
             System.out.println("(Note: wider distributions are more memory intensive to match)");
             System.out.print("Enter selection value: ");
             input = sc.nextInt();
             switch (input) {
-                case 1 -> poisson = true;
+                case 1 -> {
+                    BiGpairSEQ.setDistributionType(DistributionType.POISSON);
+                }
                 case 2 -> {
+                    BiGpairSEQ.setDistributionType(DistributionType.GAUSSIAN);
                     System.out.println("How many distinct T-cells within one standard deviation of peak frequency?");
                     System.out.println("(Note: wider distributions are more memory intensive to match)");
-                    stdDev = sc.nextDouble();
-                    if (stdDev <= 0.0) {
+                    parameter = sc.nextDouble();
+                    if (parameter <= 0.0) {
                         throw new InputMismatchException("Value must be positive.");
                     }
                 }
                 case 3 -> {
-                    exponential = true;
+                    BiGpairSEQ.setDistributionType(DistributionType.EXPONENTIAL);
                     System.out.print("Please enter lambda value for exponential distribution: ");
-                    lambda = sc.nextDouble();
-                    if (lambda <= 0.0) {
-                        lambda = 0.6;
-                        System.out.println("Value must be positive. Defaulting to 0.6.");
+                    parameter = sc.nextDouble();
+                    if (parameter <= 0.0) {
+                        parameter = 1.4;
+                        System.out.println("Value must be positive. Defaulting to 1.4.");
+                    }
+                }
+                case 4 -> {
+                    BiGpairSEQ.setDistributionType(DistributionType.ZIPF);
+                    System.out.print("Please enter exponent value for Zipf distribution: ");
+                    parameter = sc.nextDouble();
+                    if (parameter <= 0.0) {
+                        parameter = 1.4;
+                        System.out.println("Value must be positive. Defaulting to 1.4.");
                     }
                 }
                 default -> {
                     System.out.println("Invalid input. Defaulting to exponential.");
-                    exponential = true;
+                    parameter = 1.4;
+                    BiGpairSEQ.setDistributionType(DistributionType.EXPONENTIAL);
                 }
             }
             System.out.print("\nNumber of wells on plate: ");
@@ -226,16 +237,17 @@ public class InteractiveInterface {
         assert filename != null;
         Plate samplePlate;
         PlateFileWriter writer;
-        if(exponential){
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
-            writer = new PlateFileWriter(filename, samplePlate);
-        }
-        else {
-            if (poisson) {
-                stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
+        DistributionType type = BiGpairSEQ.getDistributionType();
+        switch(type) {
+            case POISSON -> {
+                parameter = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
+                samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, parameter);
+                writer = new PlateFileWriter(filename, samplePlate);
+            }
+            default -> {
+                samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, parameter);
+                writer = new PlateFileWriter(filename, samplePlate);
             }
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
-            writer = new PlateFileWriter(filename, samplePlate);
         }
         System.out.println("Writing Sample Plate to file");
         writer.writePlateFile();
@@ -582,24 +594,37 @@ public class InteractiveInterface {
         boolean backToOptions = false;
         while(!backToOptions) {
             System.out.println("\n---------ALGORITHM OPTIONS----------");
-            System.out.println("1) Use scaling algorithm by Duan and Su.");
-            System.out.println("2) Use LEDA book algorithm with Fibonacci heap priority queue");
-            System.out.println("3) Use LEDA book algorithm with pairing heap priority queue");
+            System.out.println("1) Use Hungarian algorithm with Fibonacci heap priority queue");
+            System.out.println("2) Use Hungarian algorithm with pairing heap priority queue");
+            System.out.println("3) Use auction algorithm");
+            System.out.println("4) Use integer weight scaling algorithm by Duan and Su. (buggy, not yet fully implemented!)");
             System.out.println("0) Return to Options menu");
             try {
                 input = sc.nextInt();
                 switch (input) {
-                    case 1 -> System.out.println("This option is not yet implemented. Choose another.");
-                    case 2 -> {
+                    case 1 -> {
+                        BiGpairSEQ.setHungarianAlgorithm();
                         BiGpairSEQ.setFibonacciHeap();
-                        System.out.println("MWM algorithm set to LEDA with Fibonacci heap");
+                        System.out.println("MWM algorithm set to Hungarian with Fibonacci heap");
+                        backToOptions = true;
+                    }
+                    case 2 -> {
+                        BiGpairSEQ.setHungarianAlgorithm();
+                        BiGpairSEQ.setPairingHeap();
+                        System.out.println("MWM algorithm set to Hungarian with pairing heap");
                         backToOptions = true;
                     }
                     case 3 -> {
-                        BiGpairSEQ.setPairingHeap();
-                        System.out.println("MWM algorithm set to LEDA with pairing heap");
+                        BiGpairSEQ.setAuctionAlgorithm();
+                        System.out.println("MWM algorithm set to auction");
                         backToOptions = true;
                     }
+                    case 4 -> {
+                        System.out.println("Scaling integer weight MWM algorithm not yet fully implemented. Sorry.");
+//                        BiGpairSEQ.setIntegerWeightScalingAlgorithm();
+//                        System.out.println("MWM algorithm set to integer weight scaling algorithm of Duan and Su");
+//                        backToOptions = true;
+                    }
                     case 0 -> backToOptions = true;
                     default -> System.out.println("Invalid input");
                 }

src/main/java/MaximumIntegerWeightBipartiteAuctionMatching.java

@@ -0,0 +1,177 @@
+import org.jgrapht.Graph;
+import org.jgrapht.GraphTests;
+import org.jgrapht.alg.interfaces.MatchingAlgorithm;
+
+import java.math.BigDecimal;
+import java.util.*;
+
+/**
+ * Maximum weight matching in bipartite graphs with strictly integer edge weights, using a forward auction algorithm.
+ * This implementation uses the Gauss-Seidel version of the forward auction algorithm, in which bids are submitted
+ * one at a time. For any weighted bipartite graph with n vertices in the smaller partition, this algorithm will produce
+ * a matching that is within n*epsilon of being optimal. Using an epsilon = 1/(n+1) ensures that this matching differs
+ * from an optimal matching by <1. Thus, for a bipartite graph with strictly integer weights, this algorithm returns
+ * a maximum weight matching.
+ *
+ * See:
+ * "Towards auction algorithms for large dense assignment problems"
+ * Libor Buš and Pavel Tvrdík, Comput Optim Appl (2009) 43:411-436
+ * https://link.springer.com/article/10.1007/s10589-007-9146-5
+ *
+ * See also:
+ * Many books and papers by Dimitri Bertsekas, including chapter 4 of Linear Network Optimization:
+ * https://web.mit.edu/dimitrib/www/LNets_Full_Book.pdf
+ *
+ * @param <V> the graph vertex type
+ * @param <E> the graph edge type
+ *
+ * @author Eugene Fischer
+ */
+
+public class MaximumIntegerWeightBipartiteAuctionMatching<V, E> implements MatchingAlgorithm<V, E> {
+
+    private final Graph<V, E> graph;
+    private final Set<V> partition1;
+    private final Set<V> partition2;
+    private final BigDecimal epsilon;
+    private final Set<E> matching;
+    private BigDecimal matchingWeight;
+
+    private boolean swappedPartitions = false;
+
+    public MaximumIntegerWeightBipartiteAuctionMatching(Graph<V, E> graph, Set<V> partition1, Set<V> partition2) {
+        this.graph = GraphTests.requireUndirected(graph);
+        this.partition1 = Objects.requireNonNull(partition1, "Partition 1 cannot be null");
+        this.partition2 = Objects.requireNonNull(partition2, "Partition 2 cannot be null");
+        int n = Math.max(partition1.size(), partition2.size());
+        this.epsilon = BigDecimal.valueOf(1 / ((double) n + 1)); //The minimum price increase of a bid
+        this.matching = new LinkedHashSet<>();
+        this.matchingWeight = BigDecimal.ZERO;
+    }
+
+
+    /*
+    Method coded using MaximumWeightBipartiteMatching.class from JgraphT as a model
+     */
+    @Override
+    public Matching<V, E> getMatching() {
+
+        /*
+         * Test input instance
+         */
+        if (!GraphTests.isSimple(graph)) {
+            throw new IllegalArgumentException("Only simple graphs supported");
+        }
+        if (!GraphTests.isBipartitePartition(graph, partition1, partition2)) {
+            throw new IllegalArgumentException("Graph partition is not bipartite");
+        }
+
+        /*
+        If the two partitions are different sizes, the bidders must be the smaller of the two partitions.
+         */
+        Set<V> items;
+        Set<V> bidders;
+        if (partition2.size() >= partition1.size()) {
+            bidders = partition1;
+            items = partition2;
+        }
+        else {
+            bidders = partition2;
+            items = partition1;
+            swappedPartitions = true;
+        }
+
+        /*
+        Create a map to track the owner of each item, which is initially null,
+        and a map to track the price of each item, which is initially 0. An
+        Initial price of 0 allows for asymmetric assignment (though does mean
+        that this form of the algorithm cannot take advantage of epsilon-scaling).
+         */
+        Map<V, V> owners = new HashMap<>();
+        Map<V, BigDecimal> prices = new HashMap<>();
+        for(V item: items) {
+            owners.put(item, null);
+            prices.put(item, BigDecimal.ZERO);
+        }
+
+        //Create a queue of bidders that don't currently own an item, which is initially all of them
+        Queue<V> unmatchedBidders = new ArrayDeque<>();
+        for(V bidder: bidders) {
+            unmatchedBidders.offer(bidder);
+        }
+
+        //Run the auction while there are remaining unmatched bidders
+        while (unmatchedBidders.size() > 0) {
+            V bidder = unmatchedBidders.poll();
+            V item = null;
+            BigDecimal bestValue = BigDecimal.valueOf(-1.0);
+            BigDecimal runnerUpValue = BigDecimal.valueOf(-1.0);
+            /*
+            Find the items that offer the best and second-best value for the bidder,
+            then submit a bid equal to the price of the best-valued item plus the marginal value over
+            the second-best-valued item plus epsilon.
+             */
+            for (E edge: graph.edgesOf(bidder)) {
+                double weight = graph.getEdgeWeight(edge);
+                if(weight == 0.0) {
+                    continue;
+                }
+                V tmp = getItem(edge);
+                BigDecimal value = BigDecimal.valueOf(weight).subtract(prices.get(tmp));
+                if (value.compareTo(bestValue) >= 0) {
+                    runnerUpValue = bestValue;
+                    bestValue = value;
+                    item = tmp;
+                }
+                else if (value.compareTo(runnerUpValue) >= 0) {
+                    runnerUpValue = value;
+                }
+            }
+            if(bestValue.compareTo(BigDecimal.ZERO) >= 0) {
+                V formerOwner = owners.get(item);
+                BigDecimal price = prices.get(item);
+                BigDecimal bid = price.add(bestValue).subtract(runnerUpValue).add(epsilon);
+                if (formerOwner != null) {
+                    unmatchedBidders.offer(formerOwner);
+                }
+                owners.put(item, bidder);
+                prices.put(item, bid);
+            }
+        }
+        //Add all edges between items and their owners to the matching
+        for (V item: owners.keySet()) {
+            if (owners.get(item) != null) {
+                matching.add(graph.getEdge(item, owners.get(item)));
+            }
+        }
+        //Sum the edges of the matching to obtain the matching weight
+        for(E edge: matching) {
+            this.matchingWeight = this.matchingWeight.add(BigDecimal.valueOf(graph.getEdgeWeight(edge)));
+        }
+
+        return new MatchingImpl<>(graph, matching, matchingWeight.doubleValue());
+    }
+
+    private V getItem(E edge) {
+        if (swappedPartitions) {
+            return graph.getEdgeSource(edge);
+        }
+        else {
+            return graph.getEdgeTarget(edge);
+        }
+    }
+
+//    //method for implementing a forward-reverse auction algorithm, not used here
+//    private V getBidder(E edge) {
+//        if (swappedPartitions) {
+//            return graph.getEdgeTarget(edge);
+//        }
+//        else {
+//            return graph.getEdgeSource(edge);
+//        }
+//    }
+
+    public BigDecimal getMatchingWeight() {
+        return matchingWeight;
+    }
+}

src/main/java/MaximumWeightBipartiteLookBackAuctionMatching.java

@@ -0,0 +1,212 @@
+import org.jgrapht.Graph;
+import org.jgrapht.GraphTests;
+import org.jgrapht.alg.interfaces.MatchingAlgorithm;
+import org.jgrapht.alg.util.Pair;
+
+import java.math.BigDecimal;
+import java.util.*;
+
+/*
+Maximum weight matching in bipartite graphs with strictly integer edge weights, found using the
+unscaled look-back auction algorithm
+ */
+
+public class MaximumWeightBipartiteLookBackAuctionMatching<V, E> implements MatchingAlgorithm<V, E> {
+
+    private final Graph<V, E> graph;
+    private final Set<V> partition1;
+    private final Set<V> partition2;
+    private final BigDecimal delta;
+    private final Set<E> matching;
+    private BigDecimal matchingWeight;
+    private boolean swappedPartitions = false;
+
+    public MaximumWeightBipartiteLookBackAuctionMatching(Graph<V, E> graph, Set<V> partition1, Set<V> partition2) {
+        this.graph = GraphTests.requireUndirected(graph);
+        this.partition1 = Objects.requireNonNull(partition1, "Partition 1 cannot be null");
+        this.partition2 = Objects.requireNonNull(partition2, "Partition 2 cannot be null");
+        int n = Math.max(partition1.size(), partition2.size());
+        this.delta = BigDecimal.valueOf(1 / ((double) n + 1));
+        this.matching = new LinkedHashSet<>();
+        this.matchingWeight = BigDecimal.ZERO;
+    }
+
+
+    /*
+    Method coded using MaximumWeightBipartiteMatching.class from JgraphT as a model
+     */
+    @Override
+    public Matching<V, E> getMatching() {
+
+        /*
+         * Test input instance
+         */
+        if (!GraphTests.isSimple(graph)) {
+            throw new IllegalArgumentException("Only simple graphs supported");
+        }
+        if (!GraphTests.isBipartitePartition(graph, partition1, partition2)) {
+            throw new IllegalArgumentException("Graph partition is not bipartite");
+        }
+
+        /*
+        If the two partitions are different sizes, the bidders must be the smaller of the two partitions.
+         */
+        Set<V> items;
+        Set<V> bidders;
+        if (partition2.size() >= partition1.size()) {
+            bidders = partition1;
+            items = partition2;
+        }
+        else {
+            bidders = partition2;
+            items = partition1;
+            swappedPartitions = true;
+        }
+
+        /*
+        Create a map to track the owner of each item, which is initially null,
+        and a map to track the price of each item, which is initially 0.
+         */
+        Map<V, V> owners = new HashMap<>();
+
+        /*
+        Create a map to track the prices of the objects
+         */
+        Map<V, BigDecimal> prices = new HashMap<>();
+        for(V item: items) {
+            owners.put(item, null);
+            prices.put(item, BigDecimal.ZERO);
+        }
+
+        /*
+        Create a map to track the most valuable object for a bidder
+         */
+        Map<V, V> mostValuableItems = new HashMap<>();
+
+        /*
+        Create a map to track the second most valuable object for a bidder
+         */
+        Map<V, V> runnerUpItems = new HashMap<>();
+
+        /*
+        Create a map to track the bidder value thresholds
+         */
+        Map<V, BigDecimal> valueThresholds = new HashMap<>();
+
+
+        //Initialize queue of all bidders that don't currently own an item
+        Queue<V> unmatchedBidders = new ArrayDeque<>();
+        for(V bidder: bidders) {
+            unmatchedBidders.offer(bidder);
+            valueThresholds.put(bidder, BigDecimal.ZERO);
+            mostValuableItems.put(bidder, null);
+            runnerUpItems.put(bidder, null);
+        }
+
+        while (unmatchedBidders.size() > 0) {
+            V bidder = unmatchedBidders.poll();
+//            BigDecimal valueThreshold = valueThresholds.get(bidder);
+            BigDecimal bestValue = BigDecimal.ZERO;
+            BigDecimal runnerUpValue = BigDecimal.ZERO;
+            boolean reinitialize = true;
+//            if (mostValuableItems.get(bidder) != null && runnerUpItems.get(bidder) != null) {
+//                reinitialize = false;
+//                //get the weight of the edge between the bidder and the best valued item
+//                V bestItem = mostValuableItems.get(bidder);
+//                BigDecimal bestItemWeight = BigDecimal.valueOf(graph.getEdgeWeight(graph.getEdge(bidder, bestItem)));
+//                bestValue =  bestItemWeight.subtract(prices.get(bestItem));
+//                V runnerUpItem = runnerUpItems.get(bidder);
+//                BigDecimal runnerUpWeight = BigDecimal.valueOf(graph.getEdgeWeight(graph.getEdge(bidder, runnerUpItem)));
+//                runnerUpValue = runnerUpWeight.subtract(prices.get(runnerUpItem));
+//                //if both values are still above the threshold
+//                if (bestValue.compareTo(valueThreshold) >= 0 && runnerUpValue.compareTo(valueThreshold) >= 0) {
+//                    if (bestValue.compareTo(runnerUpValue) < 0) { //if best value is lower than runner up
+//                        BigDecimal tmp = bestValue;
+//                        bestValue = runnerUpValue;
+//                        runnerUpValue = tmp;
+//                        mostValuableItems.put(bidder, runnerUpItem);
+//                        runnerUpItems.put(bidder, bestItem);
+//                    }
+//                    BigDecimal newValueThreshold = bestValue.min(runnerUpValue);
+//                    valueThresholds.put(bidder, newValueThreshold);
+//                    System.out.println("lookback successful");
+//                }
+//                else {
+//                    reinitialize = true; //lookback failed
+//                }
+//            }
+            if (reinitialize){
+                bestValue = BigDecimal.ZERO;
+                runnerUpValue = BigDecimal.ZERO;
+                for (E edge: graph.edgesOf(bidder)) {
+                    double weight = graph.getEdgeWeight(edge);
+                    if (weight == 0.0) {
+                        continue;
+                    }
+                    V tmpItem = getItem(bidder, edge);
+                    BigDecimal tmpValue = BigDecimal.valueOf(weight).subtract(prices.get(tmpItem));
+                    if (tmpValue.compareTo(bestValue) >= 0) {
+                        runnerUpValue = bestValue;
+                        bestValue = tmpValue;
+                        runnerUpItems.put(bidder, mostValuableItems.get(bidder));
+                        mostValuableItems.put(bidder, tmpItem);
+                    }
+                    else if (tmpValue.compareTo(runnerUpValue) >= 0) {
+                        runnerUpValue = tmpValue;
+                        runnerUpItems.put(bidder, tmpItem);
+                    }
+                }
+                valueThresholds.put(bidder, runnerUpValue);
+            }
+            //Should now have initialized the maps to make look back possible
+            //skip this bidder if the best value is still zero
+            if (BigDecimal.ZERO.equals(bestValue)) {
+                continue;
+            }
+            V mostValuableItem = mostValuableItems.get(bidder);
+            BigDecimal price = prices.get(mostValuableItem);
+            BigDecimal bid = price.add(bestValue).subtract(runnerUpValue).add(this.delta);
+            V formerOwner = owners.get(mostValuableItem);
+            if (formerOwner != null) {
+                    unmatchedBidders.offer(formerOwner);
+            }
+            owners.put(mostValuableItem, bidder);
+            prices.put(mostValuableItem, bid);
+        }
+
+        for (V item: owners.keySet()) {
+            if (owners.get(item) != null) {
+                matching.add(graph.getEdge(item, owners.get(item)));
+            }
+        }
+
+        for(E edge: matching) {
+            this.matchingWeight = this.matchingWeight.add(BigDecimal.valueOf(graph.getEdgeWeight(edge)));
+        }
+
+
+        return new MatchingImpl<>(graph, matching, matchingWeight.doubleValue());
+    }
+
+    private V getItem(V bidder, E edge) {
+        if (swappedPartitions) {
+            return graph.getEdgeSource(edge);
+        }
+        else {
+            return graph.getEdgeTarget(edge);
+        }
+    }
+
+    private V getBidder(V item, E edge) {
+        if (swappedPartitions) {
+            return graph.getEdgeTarget(edge);
+        }
+        else {
+            return graph.getEdgeSource(edge);
+        }
+    }
+
+    public BigDecimal getMatchingWeight() {
+        return matchingWeight;
+    }
+}

src/main/java/Plate.java

@@ -13,6 +13,11 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method
  */
 
 
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.core.BaseProvider;
+import org.apache.commons.rng.sampling.distribution.RejectionInversionZipfSampler;
+import org.apache.commons.rng.simple.JDKRandomWrapper;
+
 import java.util.*;
 
 public class Plate {
@@ -26,25 +31,22 @@ public class Plate {
     private Integer[] populations;
     private double stdDev;
     private double lambda;
-    boolean exponential = false;
+    private double zipfExponent;
+    private DistributionType distributionType;
 
     public Plate(CellSample cells, String cellFilename, int numWells, Integer[] populations,
-                 double dropoutRate, double stdDev_or_lambda, boolean exponential){
+                 double dropoutRate, double parameter){
         this.cells = cells;
         this.sourceFile = cellFilename;
         this.size = numWells;
         this.wells = new ArrayList<>();
         this.error = dropoutRate;
         this.populations = populations;
-        this.exponential = exponential;
-        if (this.exponential) {
-            this.lambda = stdDev_or_lambda;
-            fillWellsExponential(cells.getCells(), this.lambda);
-        }
-        else {
-            this.stdDev = stdDev_or_lambda;
-            fillWells(cells.getCells(), this.stdDev);
-        }
+        this.stdDev = parameter;
+        this.lambda = parameter;
+        this.zipfExponent = parameter;
+        this.distributionType = BiGpairSEQ.getDistributionType();
+        fillWells(cells.getCells());
     }
 
 
@@ -61,21 +63,57 @@ public class Plate {
         this.wells = wells;
         this.size = wells.size();
 
+        double totalCellCount = 0.0;
+        double totalDropoutCount = 0.0;
         List<Integer> concentrations = new ArrayList<>();
         for (List<String[]> w: wells) {
             if(!concentrations.contains(w.size())){
                 concentrations.add(w.size());
             }
+            for (String[] cell: w) {
+                totalCellCount += 1.0;
+                for (String sequence: cell) {
+                    if("-1".equals(sequence)) {
+                        totalDropoutCount += 1.0;
+                    }
+                }
+            }
         }
+        double totalSequenceCount = totalCellCount * 4;
+        this.error = totalDropoutCount / totalSequenceCount;
         this.populations = new Integer[concentrations.size()];
         for (int i = 0; i < this.populations.length; i++) {
             this.populations[i] = concentrations.get(i);
         }
     }
 
+    private void fillWellsZipf(List<String[]> cells, double exponent) {
+        int numSections = populations.length;
+        int section = 0;
+        int n;
+        RejectionInversionZipfSampler zipfSampler = new RejectionInversionZipfSampler(new JDKRandomWrapper(rand), cells.size(), exponent);
+        while (section < numSections){
+            for (int i = 0; i < (size / numSections); i++) {
+                List<String[]> well = new ArrayList<>();
+                for (int j = 0; j < populations[section]; j++) {
+                    do {
+                        n = zipfSampler.sample();
+                    } while (n >= cells.size() || n < 0);
+                    String[] cellToAdd = cells.get(n).clone();
+                    for(int k = 0; k < cellToAdd.length; k++){
+                        if(Math.abs(rand.nextDouble()) < error){//error applied to each sequence
+                            cellToAdd[k] = "-1";
+                        }
+                    }
+                    well.add(cellToAdd);
+                }
+                wells.add(well);
+            }
+            section++;
+        }
+    }
+
     private void fillWellsExponential(List<String[]> cells, double lambda){
-        this.lambda = lambda;
-        exponential = true;
         int numSections = populations.length;
         int section = 0;
         double m;
@@ -131,6 +169,24 @@ public class Plate {
         }
     }
 
+    private void fillWells(List<String[]> cells){
+        DistributionType type = BiGpairSEQ.getDistributionType();
+        switch (type) {
+            case POISSON, GAUSSIAN -> {
+                fillWells(cells, getStdDev());
+                break;
+            }
+            case EXPONENTIAL -> {
+                fillWellsExponential(cells, getLambda());
+                break;
+            }
+            case ZIPF -> {
+                fillWellsZipf(cells, getZipfExponent());
+                break;
+            }
+        }
+    }
+
     public Integer[] getPopulations(){
         return populations;
     }
@@ -143,10 +199,12 @@ public class Plate {
         return stdDev;
     }
 
-    public boolean isExponential(){return exponential;}
+    public DistributionType getDistributionType() { return distributionType;}
 
     public double getLambda(){return lambda;}
 
+    public double getZipfExponent(){return zipfExponent;}
+
     public double getError() {
         return error;
     }

src/main/java/PlateFileWriter.java

@@ -13,11 +13,13 @@ public class PlateFileWriter {
     private List<List<String[]>> wells;
     private double stdDev;
     private double lambda;
+    private double zipfExponent;
+    private DistributionType distributionType;
     private Double error;
     private String filename;
     private String sourceFileName;
     private Integer[] populations;
-    private boolean isExponential = false;
+
 
     public PlateFileWriter(String filename, Plate plate) {
         if(!filename.matches(".*\\.csv")){
@@ -26,12 +28,17 @@ public class PlateFileWriter {
         this.filename = filename;
         this.sourceFileName = plate.getSourceFileName();
         this.size = plate.getSize();
-        this.isExponential = plate.isExponential();
-        if(isExponential) {
-            this.lambda = plate.getLambda();
-        }
-        else{
-            this.stdDev = plate.getStdDev();
+        this.distributionType = plate.getDistributionType();
+        switch(distributionType) {
+            case POISSON, GAUSSIAN -> {
+                this.stdDev = plate.getStdDev();
+            }
+            case EXPONENTIAL -> {
+                this.lambda = plate.getLambda();
+            }
+            case ZIPF -> {
+                this.zipfExponent = plate.getZipfExponent();
+            }
         }
         this.error = plate.getError();
         this.wells = plate.getWells();
@@ -93,13 +100,24 @@ public class PlateFileWriter {
             printer.printComment("Cell source file name: " + sourceFileName);
             printer.printComment("Each row represents one well on the plate.");
             printer.printComment("Plate size: " + size);
-            printer.printComment("Error rate: " + error);
             printer.printComment("Well populations: " + wellPopulationsString);
-            if(isExponential){
-                printer.printComment("Lambda: " + lambda);
-            }
-            else {
-                printer.printComment("Std. dev.: " + stdDev);
+            printer.printComment("Error rate: " + error);
+            switch (distributionType) {
+                case POISSON -> {
+                    printer.printComment("Cell frequency distribution: POISSON");
+                }
+                case GAUSSIAN -> {
+                    printer.printComment("Cell frequency distribution: GAUSSIAN");
+                    printer.printComment("--Standard deviation: " + stdDev);
+                }
+                case EXPONENTIAL -> {
+                    printer.printComment("Cell frequency distribution: EXPONENTIAL");
+                    printer.printComment("--Lambda: " + lambda);
+                }
+                case ZIPF -> {
+                    printer.printComment("Cell frequency distribution: ZIPF");
+                    printer.printComment("--Exponent: " + zipfExponent);
+                }
             }
             printer.printRecords(wellsAsStrings);
         } catch(IOException ex){

src/main/java/Simulator.java

@@ -183,32 +183,36 @@ public class Simulator implements GraphModificationFunctions {
                 "removed");}
 
         //Find Maximum Weight Matching
-        //using jheaps library class PairingHeap for improved efficiency
         if(verbose){System.out.println("Finding maximum weight matching");}
-        MaximumWeightBipartiteMatching maxWeightMatching;
-        //Use correct heap type for priority queue
-        String heapType = BiGpairSEQ.getPriorityQueueHeapType();
-        switch (heapType) {
-            case "PAIRING" -> {
-                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        alphas,
-                        betas,
-                        i -> new PairingHeap(Comparator.naturalOrder()));
+        //The matching object
+        MatchingAlgorithm<Vertex, DefaultWeightedEdge> maxWeightMatching;
+        //Determine algorithm type
+        AlgorithmType algorithm = BiGpairSEQ.getMatchingAlgoritmType();
+        switch (algorithm) { //Only two options now, but I have room to add more algorithms in the future this way
+            case AUCTION -> {
+                //create a new MaximumIntegerWeightBipartiteAuctionMatching
+                maxWeightMatching = new MaximumIntegerWeightBipartiteAuctionMatching<>(graph, alphas, betas);
             }
-            case "FIBONACCI" -> {
-                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        alphas,
-                        betas,
-                        i -> new FibonacciHeap(Comparator.naturalOrder()));
+            case INTEGER_WEIGHT_SCALING -> {
+                maxWeightMatching = new MaximumIntegerWeightBipartiteMatching<>(graph, alphas, betas, new BigDecimal(highThreshold));
             }
-            default -> {
-                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
-                        alphas,
-                        betas);
+            default -> { //HUNGARIAN
+                //use selected heap type for priority queue
+                HeapType heap = BiGpairSEQ.getPriorityQueueHeapType();
+                if(HeapType.PAIRING.equals(heap)) {
+                    maxWeightMatching = new MaximumWeightBipartiteMatching<Vertex, DefaultWeightedEdge>(graph,
+                            alphas,
+                            betas,
+                            i -> new PairingHeap(Comparator.naturalOrder()));
+                }
+                else {//Fibonacci is the default, and what's used in the JGraphT implementation
+                    maxWeightMatching = new MaximumWeightBipartiteMatching<Vertex, DefaultWeightedEdge>(graph,
+                            alphas,
+                            betas);
+                }
             }
         }
-        //get the matching
-        MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
+        MatchingAlgorithm.Matching<Vertex, DefaultWeightedEdge> matching = maxWeightMatching.getMatching();
         if(verbose){System.out.println("Matching completed");}
         Instant stop = Instant.now();
 
@@ -226,7 +230,7 @@ public class Simulator implements GraphModificationFunctions {
         List<List<String>> allResults = new ArrayList<>();
         NumberFormat nf = NumberFormat.getInstance(Locale.US);
         MathContext mc = new MathContext(3);
-        Iterator<DefaultWeightedEdge> weightIter = graphMatching.iterator();
+        Iterator<DefaultWeightedEdge> weightIter = matching.iterator();
         DefaultWeightedEdge e;
         int trueCount = 0;
         int falseCount = 0;
@@ -267,10 +271,22 @@ public class Simulator implements GraphModificationFunctions {
         }
 
         //Metadata comments for CSV file
-        String algoType = "LEDA book with heap: " + heapType;
+        String algoType;
+        switch(algorithm) {
+            case AUCTION -> {
+                algoType = "Auction algorithm";
+            }
+            case INTEGER_WEIGHT_SCALING -> {
+                algoType = "Integer weight scaling algorithm from Duan and Su (not yet perfectly implemented)";
+            }
+            default -> { //HUNGARIAN
+                algoType = "Hungarian algorithm with heap: " + BiGpairSEQ.getPriorityQueueHeapType().name();
+            }
+        }
+
         int min = Math.min(graphAlphaCount, graphBetaCount);
         //matching weight
-        BigDecimal totalMatchingWeight = maxWeightMatching.getMatchingWeight();
+        Double matchingWeight = matching.getWeight();
         //rate of attempted matching
         double attemptRate = (double) (trueCount + falseCount) / min;
         BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
@@ -309,7 +325,7 @@ public class Simulator implements GraphModificationFunctions {
         metadata.put("sequence dropout rate", data.getDropoutRate().toString());
         metadata.put("graph filename", dataFilename);
         metadata.put("MWM algorithm type", algoType);
-        metadata.put("matching weight", totalMatchingWeight.toString());
+        metadata.put("matching weight", matchingWeight.toString());
         metadata.put("well populations", wellPopulationsString);
         metadata.put("sequence read depth", data.getReadDepth().toString());
         metadata.put("sequence read error rate", data.getReadErrorRate().toString());
@@ -347,6 +363,7 @@ public class Simulator implements GraphModificationFunctions {
         return output;
     }
 
+
     //Commented out CDR1 matching until it's time to re-implement it
 //    //Simulated matching of CDR1s to CDR3s. Requires MatchingResult from prior run of matchCDR3s.
 //    public static MatchingResult[] matchCDR1s(List<Integer[]> distinctCells,

src/main/java/Vertex.java

@@ -74,4 +74,12 @@ public class Vertex implements Serializable, Comparable<Vertex> {
     public int compareTo(Vertex other) {
         return this.vertexLabel - other.getVertexLabel();
     }
+
+    public Double getPotential() {
+        return potential;
+    }
+
+    public void setPotential(Double potential) {
+        this.potential = potential;
+    }
 }