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

@@ -8,19 +8,22 @@
    1. [RUNNING THE PROGRAM](#running-the-program)
    2. [COMMAND LINE OPTIONS](#command-line-options)
    3. [INTERACTIVE INTERFACE](#interactive-interface)
-   4. [INPUT/OUTPUT](#inputoutput)
-      1. Cell Sample Files
-      2. Sample Plate Files
-      3. Graph/Data Files
-      4. Matching Results Files
+   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
-   2. SIMULATING EXPERIMENTS FROM pairSEQ PAPER
-6. [TODO](#todo)
-7. [CITATIONS](#citations)
+   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
 
@@ -74,7 +77,7 @@ 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 algorithim eschews epsilon-scaling, 
+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 
@@ -153,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
@@ -170,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
@@ -231,7 +237,6 @@ usage: BiGpairSEQ_Sim.jar -match
                                to stdout.
  -pv,--p-value                 (Optional) Calculate p-values for sequence
                                pairs.
-
 ```
 
 ### INTERACTIVE INTERFACE
@@ -337,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
@@ -480,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
@@ -498,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
@@ -590,8 +601,40 @@ 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.
   
+## 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
+* [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
+* [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
+* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
+
+## ACKNOWLEDGEMENTS
+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–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.*~~
@@ -600,7 +643,7 @@ the file of distinct cells may enable better simulated replication of this exper
 * ~~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.
+        * 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.
@@ -625,13 +668,13 @@ the file of distinct cells may enable better simulated replication of this exper
 * 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.
+* ~~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?
+    * ~~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
@@ -642,35 +685,3 @@ the file of distinct cells may enable better simulated replication of this exper
     * 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. 2020. JGraphT—A Java Library for Graph Data Structures and Algorithms. ACM Trans. Math. Softw. 46, 2, Article 16
-
-## EXTERNAL LIBRARIES USED
-* [JGraphT](https://jgrapht.org) -- Graph theory data structures and algorithms
-* [JHeaps](https://www.jheaps.org) -- For pairing heap priority queue used in maximum weight matching algorithm
-* [Apache Commons CSV](https://commons.apache.org/proper/commons-csv/) -- For CSV file output
-* [Apache Commons CLI](https://commons.apache.org/proper/commons-cli/) -- To enable command line arguments for scripting.
-
-## ACKNOWLEDGEMENTS
-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–2023.
-
-## 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 has been 
-employed at Adaptive Biotechnologies since 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 BiGpairSEQ
-recreationally, as it involves some very beautiful math.

src/main/java/BiGpairSEQ.java

@@ -15,6 +15,7 @@ public class BiGpairSEQ {
     private static boolean cacheGraph = false;
     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;
@@ -60,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;
     }

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,17 +237,18 @@ public class InteractiveInterface {
         assert filename != null;
         Plate samplePlate;
         PlateFileWriter writer;
-        if(exponential){
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
+        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);
             }
-        else {
-            if (poisson) {
-                stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
-            }
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
+            default -> {
+                samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, parameter);
                 writer = new PlateFileWriter(filename, samplePlate);
             }
+        }
         System.out.println("Writing Sample Plate to file");
         writer.writePlateFile();
         System.out.println("Sample Plate written to file: " + filename);
@@ -605,12 +617,13 @@ public class InteractiveInterface {
                     case 3 -> {
                         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;
+//                        backToOptions = true;
                     }
                     case 0 -> backToOptions = true;
                     default -> System.out.println("Invalid input");

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());
     }
 
 
@@ -85,9 +87,33 @@ public class Plate {
         }
     }
 
+    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;
@@ -143,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;
     }
@@ -155,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.distributionType = plate.getDistributionType();
+        switch(distributionType) {
+            case POISSON, GAUSSIAN -> {
+                this.stdDev = plate.getStdDev();
+            }
+            case EXPONENTIAL -> {
                 this.lambda = plate.getLambda();
             }
-        else{
-            this.stdDev = plate.getStdDev();
+            case ZIPF -> {
+                this.zipfExponent = plate.getZipfExponent();
+            }
         }
         this.error = plate.getError();
         this.wells = plate.getWells();
@@ -95,11 +102,22 @@ public class PlateFileWriter {
             printer.printComment("Plate size: " + size);
             printer.printComment("Well populations: " + wellPopulationsString);
             printer.printComment("Error rate: " + error);
-            if(isExponential){
-                printer.printComment("Lambda: " + lambda);
+            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);
                 }
-            else {
-                printer.printComment("Std. dev.: " + stdDev);
             }
             printer.printRecords(wellsAsStrings);
         } catch(IOException ex){