update ToDos

.idea/.name

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

.idea/artifacts/BiGpairSEQ_Sim_jar.xml

@@ -1,16 +1,27 @@
 <component name="ArtifactManager">
-  <artifact type="jar" build-on-make="true" name="BiGpairSEQ_Sim:jar">
+  <artifact type="jar" name="BiGpairSEQ_Sim:jar">
     <output-path>$PROJECT_DIR$/out/artifacts/BiGpairSEQ_Sim_jar</output-path>
     <root id="archive" name="BiGpairSEQ_Sim.jar">
       <element id="directory" name="META-INF">
-        <element id="file-copy" path="$PROJECT_DIR$/src/main/java/META-INF/MANIFEST.MF" />
+        <element id="file-copy" path="$PROJECT_DIR$/META-INF/MANIFEST.MF" />
       </element>
-      <element id="module-output" name="BigPairSEQ" />
-      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.1/jgrapht-core-1.5.1.jar" path-in-jar="/" />
-      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.13/jheaps-0.13.jar" path-in-jar="/" />
-      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/commons-cli/commons-cli/1.5.0/commons-cli-1.5.0.jar" path-in-jar="/" />
-      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-csv/1.9.0/commons-csv-1.9.0.jar" path-in-jar="/" />
-      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jetbrains/annotations/23.0.0/annotations-23.0.0.jar" path-in-jar="/" />
+      <element id="module-output" name="BiGpairSEQ_Sim" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-rng-sampling/1.6/commons-rng-sampling-1.6.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-csv/1.14.0/commons-csv-1.14.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jetbrains/annotations/26.0.2/annotations-26.0.2.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-io/1.5.2/jgrapht-io-1.5.2.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-rng-simple/1.6/commons-rng-simple-1.6.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/commons-io/commons-io/2.18.0/commons-io-2.18.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-rng-core/1.6/commons-rng-core-1.6.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/commons-codec/commons-codec/1.18.0/commons-codec-1.18.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-rng-client-api/1.6/commons-rng-client-api-1.6.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/commons-cli/commons-cli/1.9.0/commons-cli-1.9.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-lang3/3.12.0/commons-lang3-3.12.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/antlr/antlr4-runtime/4.12.0/antlr4-runtime-4.12.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apfloat/apfloat/1.10.1/apfloat-1.10.1.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/apache/commons/commons-text/1.10.0/commons-text-1.10.0.jar" path-in-jar="/" />
+      <element id="extracted-dir" path="$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.14/jheaps-0.14.jar" path-in-jar="/" />
     </root>
   </artifact>
 </component>

.idea/compiler.xml

@@ -7,6 +7,7 @@
         <sourceTestOutputDir name="target/generated-test-sources/test-annotations" />
         <outputRelativeToContentRoot value="true" />
         <module name="BigPairSEQ" />
+        <module name="BiGpairSEQ_Sim" />
       </profile>
     </annotationProcessing>
   </component>

.idea/jarRepositories.xml

@@ -1,20 +1,35 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
   <component name="RemoteRepositoriesConfiguration">
+    <remote-repository>
+      <option name="id" value="my-internal-site" />
+      <option name="name" value="my-internal-site" />
+      <option name="url" value="https://myserver/repo" />
+    </remote-repository>
+    <remote-repository>
+      <option name="id" value="central" />
+      <option name="name" value="Central Repository" />
+      <option name="url" value="https://repo1.maven.org/maven2" />
+    </remote-repository>
+    <remote-repository>
+      <option name="id" value="central repo" />
+      <option name="name" value="central repo" />
+      <option name="url" value="https://repo1.maven.org/maven2/" />
+    </remote-repository>
     <remote-repository>
       <option name="id" value="central" />
       <option name="name" value="Central Repository" />
       <option name="url" value="https://repo.maven.apache.org/maven2" />
     </remote-repository>
-    <remote-repository>
-      <option name="id" value="central" />
-      <option name="name" value="Maven Central repository" />
-      <option name="url" value="https://repo1.maven.org/maven2" />
-    </remote-repository>
     <remote-repository>
       <option name="id" value="jboss.community" />
       <option name="name" value="JBoss Community repository" />
       <option name="url" value="https://repository.jboss.org/nexus/content/repositories/public/" />
     </remote-repository>
+    <remote-repository>
+      <option name="id" value="34d16bdc-85f0-48ee-8e8b-144091765be1" />
+      <option name="name" value="34d16bdc-85f0-48ee-8e8b-144091765be1" />
+      <option name="url" value="https://repository.mulesoft.org/nexus/content/repositories/public/" />
+    </remote-repository>
   </component>
 </project>

.idea/libraries/apache_commons_csv.xml

@@ -1,8 +1,10 @@
 <component name="libraryTable">
   <library name="apache.commons.csv" type="repository">
-    <properties maven-id="org.apache.commons:commons-csv:1.9.0" />
+    <properties maven-id="org.apache.commons:commons-csv:1.14.0" />
     <CLASSES>
-      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-csv/1.9.0/commons-csv-1.9.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-csv/1.14.0/commons-csv-1.14.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/commons-io/commons-io/2.18.0/commons-io-2.18.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/commons-codec/commons-codec/1.18.0/commons-codec-1.18.0.jar!/" />
     </CLASSES>
     <JAVADOC />
     <SOURCES />

.idea/libraries/commons_cli.xml

@@ -1,8 +1,8 @@
 <component name="libraryTable">
   <library name="commons.cli" type="repository">
-    <properties maven-id="commons-cli:commons-cli:1.5.0" />
+    <properties maven-id="commons-cli:commons-cli:1.9.0" />
     <CLASSES>
-      <root url="jar://$MAVEN_REPOSITORY$/commons-cli/commons-cli/1.5.0/commons-cli-1.5.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/commons-cli/commons-cli/1.9.0/commons-cli-1.9.0.jar!/" />
     </CLASSES>
     <JAVADOC />
     <SOURCES />

.idea/libraries/jgrapht_core.xml

@@ -1,9 +1,10 @@
 <component name="libraryTable">
   <library name="jgrapht.core" type="repository">
-    <properties maven-id="org.jgrapht:jgrapht-core:1.5.1" />
+    <properties maven-id="org.jgrapht:jgrapht-core:1.5.2" />
     <CLASSES>
-      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.1/jgrapht-core-1.5.1.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.13/jheaps-0.13.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.14/jheaps-0.14.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/apfloat/apfloat/1.10.1/apfloat-1.10.1.jar!/" />
     </CLASSES>
     <JAVADOC />
     <SOURCES />

.idea/libraries/jgrapht_io.xml

@@ -1,13 +1,14 @@
 <component name="libraryTable">
   <library name="jgrapht.io" type="repository">
-    <properties maven-id="org.jgrapht:jgrapht-io:1.5.1" />
+    <properties maven-id="org.jgrapht:jgrapht-io:1.5.2" />
     <CLASSES>
-      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-io/1.5.1/jgrapht-io-1.5.1.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.1/jgrapht-core-1.5.1.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.13/jheaps-0.13.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/antlr/antlr4-runtime/4.8-1/antlr4-runtime-4.8-1.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-text/1.8/commons-text-1.8.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-lang3/3.9/commons-lang3-3.9.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-io/1.5.2/jgrapht-io-1.5.2.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.14/jheaps-0.14.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/apfloat/apfloat/1.10.1/apfloat-1.10.1.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/antlr/antlr4-runtime/4.12.0/antlr4-runtime-4.12.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-text/1.10.0/commons-text-1.10.0.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-lang3/3.12.0/commons-lang3-3.12.0.jar!/" />
     </CLASSES>
     <JAVADOC />
     <SOURCES />

pom.xml

@@ -5,7 +5,7 @@
     <modelVersion>4.0.0</modelVersion>
 
     <groupId>org.example</groupId>
-    <artifactId>TCellSim</artifactId>
+    <artifactId>BiGpairSEQ_Sim</artifactId>
     <version>1.0-SNAPSHOT</version>
     <build>
         <plugins>
@@ -26,8 +26,48 @@
             <version>RELEASE</version>
             <scope>compile</scope>
         </dependency>
+        <!-- https://mvnrepository.com/artifact/org.apache.commons/commons-rng-simple -->
+        <dependency>
+            <groupId>org.apache.commons</groupId>
+            <artifactId>commons-rng-simple</artifactId>
+            <version>1.6</version>
+        </dependency>
+        <dependency>
+            <groupId>org.apache.commons</groupId>
+            <artifactId>commons-rng-sampling</artifactId>
+            <version>1.6</version>
+        </dependency>
+        <!-- https://mvnrepository.com/artifact/org.apache.commons/commons-csv -->
+        <dependency>
+            <groupId>org.apache.commons</groupId>
+            <artifactId>commons-csv</artifactId>
+            <version>1.14.0</version>
+        </dependency>
+        <!-- https://mvnrepository.com/artifact/org.jgrapht/jgrapht-core -->
+        <dependency>
+            <groupId>org.jgrapht</groupId>
+            <artifactId>jgrapht-core</artifactId>
+            <version>1.5.2</version>
+        </dependency>
+        <!-- https://mvnrepository.com/artifact/org.jgrapht/jgrapht-io -->
+        <dependency>
+            <groupId>org.jgrapht</groupId>
+            <artifactId>jgrapht-io</artifactId>
+            <version>1.5.2</version>
+        </dependency>
+        <!-- https://mvnrepository.com/artifact/org.jheaps/jheaps -->
+        <dependency>
+            <groupId>org.jheaps</groupId>
+            <artifactId>jheaps</artifactId>
+            <version>0.14</version>
+        </dependency>
+        <!-- https://mvnrepository.com/artifact/commons-cli/commons-cli -->
+        <dependency>
+            <groupId>commons-cli</groupId>
+            <artifactId>commons-cli</artifactId>
+            <version>1.9.0</version>
+        </dependency>
     </dependencies>
-
     <properties>
         <maven.compiler.source>11</maven.compiler.source>
         <maven.compiler.target>11</maven.compiler.target>

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 
@@ -133,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
@@ -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,43 @@ 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
 
+* Consider whether a graph database might be a better option than keeping things in memory.
+* Look at fastUtil for more performant maps and arrays. Note that there is an optional jGraphT library to work with fastUtil (see FastutilMapIntVertexGraph, for example).
+* Consider implementing an option to use the jGrapht sparse graph representation for a lower memory cost with very large graphs (tens or hundreds of thousands of distinct sequences).
+* ~~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 +646,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 +671,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 +688,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

@@ -13,8 +13,9 @@ public class BiGpairSEQ {
     private static boolean cacheCells = false;
     private static boolean cachePlate = false;
     private static boolean cacheGraph = false;
-    private static AlgorithmType matchingAlgoritmType = AlgorithmType.HUNGARIAN;
+    private static AlgorithmType matchingAlgorithmType = 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;
     }
@@ -161,13 +166,13 @@ public class BiGpairSEQ {
         return priorityQueueHeapType;
     }
 
-    public static AlgorithmType getMatchingAlgoritmType() { return matchingAlgoritmType; }
+    public static AlgorithmType getMatchingAlgorithmType() { return matchingAlgorithmType; }
 
-    public static void setHungarianAlgorithm() { matchingAlgoritmType = AlgorithmType.HUNGARIAN; }
+    public static void setHungarianAlgorithm() { matchingAlgorithmType = AlgorithmType.HUNGARIAN; }
 
-    public static void setIntegerWeightScalingAlgorithm() { matchingAlgoritmType = AlgorithmType.INTEGER_WEIGHT_SCALING; }
+    public static void setIntegerWeightScalingAlgorithm() { matchingAlgorithmType = AlgorithmType.INTEGER_WEIGHT_SCALING; }
 
-    public static void setAuctionAlgorithm() { matchingAlgoritmType = AlgorithmType.AUCTION; }
+    public static void setAuctionAlgorithm() { matchingAlgorithmType = 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/GraphModificationFunctions.java

@@ -1,72 +1,54 @@
 import org.jgrapht.graph.DefaultWeightedEdge;
 import org.jgrapht.graph.SimpleWeightedGraph;
 
-import java.util.ArrayList;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
+import java.util.*;
 
 public interface GraphModificationFunctions {
 
     //remove over- and under-weight edges, return removed edges
-    static Map<Vertex[], Integer> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
+    static Map<DefaultWeightedEdge, Vertex[]> filterByOverlapThresholds(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                               int low, int high, boolean saveEdges) {
-        Map<Vertex[], Integer> removedEdges = new HashMap<>();
+        Map<DefaultWeightedEdge, Vertex[]> removedEdges = new HashMap<>();
+        Set<DefaultWeightedEdge> edgesToRemove = new HashSet<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
                 if(saveEdges) {
-                    Vertex source = graph.getEdgeSource(e);
-                    Vertex target = graph.getEdgeTarget(e);
-                    Integer weight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, weight);
+                    Vertex[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
+                    removedEdges.put(e, vertices);
                 }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
-                }
-            }
-        }
-        if(saveEdges) {
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
+                edgesToRemove.add(e);
             }
         }
+        edgesToRemove.forEach(graph::removeEdge);
         return removedEdges;
     }
 
     //Remove edges for pairs with large occupancy discrepancy, return removed edges
-    static Map<Vertex[], Integer> filterByRelativeOccupancy(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
+    static Map<DefaultWeightedEdge, Vertex[]> filterByRelativeOccupancy(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                   Integer maxOccupancyDifference, boolean saveEdges) {
-        Map<Vertex[], Integer> removedEdges = new HashMap<>();
+        Map<DefaultWeightedEdge, Vertex[]> removedEdges = new HashMap<>();
+        Set<DefaultWeightedEdge> edgesToRemove = new HashSet<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
             Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
             if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
                 if (saveEdges) {
-                    Vertex source = graph.getEdgeSource(e);
-                    Vertex target = graph.getEdgeTarget(e);
-                    Integer weight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, weight);
+                    Vertex[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
+                    removedEdges.put(e, vertices);
                 }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
-                }
-            }
-        }
-        if(saveEdges) {
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
+                edgesToRemove.add(e);
             }
         }
+        edgesToRemove.forEach(graph::removeEdge);
         return removedEdges;
     }
 
     //Remove edges for pairs where overlap size is significantly lower than the well occupancy, return removed edges
-    static Map<Vertex[], Integer> filterByOverlapPercent(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
+    static Map<DefaultWeightedEdge, Vertex[]> filterByOverlapPercent(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
                                                          Integer minOverlapPercent,
                                                          boolean saveEdges) {
-        Map<Vertex[], Integer> removedEdges = new HashMap<>();
+        Map<DefaultWeightedEdge, Vertex[]> removedEdges = new HashMap<>();
+        Set<DefaultWeightedEdge> edgesToRemove = new HashSet<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = graph.getEdgeSource(e).getOccupancy();
             Integer betaOcc = graph.getEdgeTarget(e).getOccupancy();
@@ -74,22 +56,13 @@ public interface GraphModificationFunctions {
             double min = minOverlapPercent / 100.0;
             if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
                 if (saveEdges) {
-                    Vertex source = graph.getEdgeSource(e);
-                    Vertex target = graph.getEdgeTarget(e);
-                    Integer intWeight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, intWeight);
+                    Vertex[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
+                    removedEdges.put(e, vertices);
                 }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
-                }
-            }
-        }
-        if(saveEdges) {
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
+                edgesToRemove.add(e);
             }
         }
+        edgesToRemove.forEach(graph::removeEdge);
         return removedEdges;
     }
 
@@ -126,10 +99,10 @@ public interface GraphModificationFunctions {
     }
 
     static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
-                                Map<Vertex[], Integer> removedEdges) {
-        for (Vertex[] edge : removedEdges.keySet()) {
-            DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
-            graph.setEdgeWeight(e, removedEdges.get(edge));
+                                Map<DefaultWeightedEdge, Vertex[]> removedEdges) {
+        for (DefaultWeightedEdge edge : removedEdges.keySet()) {
+            Vertex[] vertices = removedEdges.get(edge);
+            graph.addEdge(vertices[0], vertices[1], edge);
         }
     }
 

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,10 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method
  */
 
 
+
+import org.apache.commons.rng.sampling.distribution.RejectionInversionZipfSampler;
+import org.apache.commons.rng.simple.JDKRandomWrapper;
+
 import java.util.*;
 
 public class Plate {
@@ -26,25 +30,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 +86,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 +168,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 +198,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;
     }
@@ -196,7 +241,7 @@ public class Plate {
                                     sequencesAndMisreads.put(currentSequence, new ArrayList<>());
                                 }
                                 //The specific misread hasn't happened before
-                                if (rand.nextDouble() >= errorCollisionRate || sequencesAndMisreads.get(currentSequence).size() == 0) {
+                                if (rand.nextDouble() >= errorCollisionRate || sequencesAndMisreads.get(currentSequence).isEmpty()) {
                                     //The misread doesn't collide with a real sequence already on the plate and some sequences have already been read
                                     if(rand.nextDouble() >= realSequenceCollisionRate || !sequenceMap.isEmpty()){
                                         StringBuilder spurious = new StringBuilder(currentSequence);

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){

src/main/java/Simulator.java

@@ -1,9 +1,8 @@
+import org.jgrapht.Graphs;
 import org.jgrapht.alg.interfaces.MatchingAlgorithm;
 import org.jgrapht.alg.matching.MaximumWeightBipartiteMatching;
-import org.jgrapht.generate.SimpleWeightedBipartiteGraphMatrixGenerator;
 import org.jgrapht.graph.DefaultWeightedEdge;
 import org.jgrapht.graph.SimpleWeightedGraph;
-import org.jheaps.tree.FibonacciHeap;
 import org.jheaps.tree.PairingHeap;
 
 import java.math.BigDecimal;
@@ -70,58 +69,104 @@ public class Simulator implements GraphModificationFunctions {
             if(verbose){System.out.println("Total beta sequence wells removed: " + betaWellsRemoved);}
         }
 
-        //construct the graph. For simplicity, going to make
-        if(verbose){System.out.println("Making vertex maps");}
-        //For the SimpleWeightedBipartiteGraphMatrixGenerator, all vertices must have
-        //distinct numbers associated with them. Since I'm using a 2D array, that means
-        //distinct indices between the rows and columns. vertexStartValue lets me track where I switch
-        //from numbering rows to columns, so I can assign unique numbers to every vertex, and then
-        //subtract the vertexStartValue from betas to use their vertex labels as array indices
-        int vertexStartValue = 0;
-        //keys are sequential integer vertices, values are alphas
-        Map<String, Integer> plateAtoVMap = makeSequenceToVertexMap(alphaSequences, vertexStartValue);
-        //new start value for vertex to beta map should be one more than final vertex value in alpha map
-        vertexStartValue += plateAtoVMap.size();
-        //keys are betas, values are sequential integers
-        Map<String, Integer> plateBtoVMap = makeSequenceToVertexMap(betaSequences, vertexStartValue);
-        if(verbose){System.out.println("Vertex maps made");}
-        //make adjacency matrix for bipartite graph generator
-        //(technically this is only 1/4 of an adjacency matrix, but that's all you need
-        //for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
-        if(verbose){System.out.println("Making adjacency matrix");}
-        double[][] weights = new double[plateAtoVMap.size()][plateBtoVMap.size()];
-        fillAdjacencyMatrix(weights, vertexStartValue, alphaSequences, betaSequences, plateAtoVMap, plateBtoVMap);
-        if(verbose){System.out.println("Adjacency matrix made");}
+        /*
+         * The commented out code below works beautifully for small enough graphs. However, after implementing a
+         * Zipf distribution and attempting to simulate Experiment 3 from the paper again, I discovered that
+         * this method uses too much memory. Even a 120GB heap is not enough to build this adjacency matrix.
+         * So I'm going to attempt to build this graph directly and see if that is less memory intensive
+         */
+//        //construct the graph. For simplicity, going to make
+//        if(verbose){System.out.println("Making vertex maps");}
+//        //For the SimpleWeightedBipartiteGraphMatrixGenerator, all vertices must have
+//        //distinct numbers associated with them. Since I'm using a 2D array, that means
+//        //distinct indices between the rows and columns. vertexStartValue lets me track where I switch
+//        //from numbering rows to columns, so I can assign unique numbers to every vertex, and then
+//        //subtract the vertexStartValue from betas to use their vertex labels as array indices
+//        int vertexStartValue = 0;
+//        //keys are sequential integer vertices, values are alphas
+//        Map<String, Integer> plateAtoVMap = makeSequenceToVertexMap(alphaSequences, vertexStartValue);
+//        //new start value for vertex to beta map should be one more than final vertex value in alpha map
+//        vertexStartValue += plateAtoVMap.size();
+//        //keys are betas, values are sequential integers
+//        Map<String, Integer> plateBtoVMap = makeSequenceToVertexMap(betaSequences, vertexStartValue);
+//        if(verbose){System.out.println("Vertex maps made");}
+//        //make adjacency matrix for bipartite graph generator
+//        //(technically this is only 1/4 of an adjacency matrix, but that's all you need
+//        //for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
+//        if(verbose){System.out.println("Making adjacency matrix");}
+//        double[][] weights = new double[plateAtoVMap.size()][plateBtoVMap.size()];
+//        fillAdjacencyMatrix(weights, vertexStartValue, alphaSequences, betaSequences, plateAtoVMap, plateBtoVMap);
+//        if(verbose){System.out.println("Adjacency matrix made");}
+//        //make bipartite graph
+//        if(verbose){System.out.println("Making bipartite weighted graph");}
+//        //the graph object
+//        SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph =
+//                new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
+//        //the graph generator
+//        SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
+//        //the list of alpha vertices
+//        List<Vertex> alphaVertices = new ArrayList<>();
+//        for (String seq : plateAtoVMap.keySet()) {
+//            Vertex alphaVertex = new Vertex(alphaSequences.get(seq), plateAtoVMap.get(seq));
+//            alphaVertices.add(alphaVertex);
+//        }
+//        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
+//        Collections.sort(alphaVertices);
+//        //Add ordered list of vertices to the graph
+//        graphGenerator.first(alphaVertices);
+//        //the list of beta vertices
+//        List<Vertex> betaVertices = new ArrayList<>();
+//        for (String seq : plateBtoVMap.keySet()) {
+//            Vertex betaVertex = new Vertex(betaSequences.get(seq), plateBtoVMap.get(seq));
+//            betaVertices.add(betaVertex);
+//        }
+//        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
+//        Collections.sort(betaVertices);
+//        //Add ordered list of vertices to the graph
+//        graphGenerator.second(betaVertices);
+//        //use adjacency matrix of weight created previously
+//        graphGenerator.weights(weights);
+//        graphGenerator.generateGraph(graph);
+
         //make bipartite graph
         if(verbose){System.out.println("Making bipartite weighted graph");}
         //the graph object
         SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph =
                 new SimpleWeightedGraph<>(DefaultWeightedEdge.class);
-        //the graph generator
-        SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
-        //the list of alpha vertices
+        int vertexLabelValue = 0;
+        //create and add alpha sequence vertices
         List<Vertex> alphaVertices = new ArrayList<>();
-        for (String seq : plateAtoVMap.keySet()) {
-            Vertex alphaVertex = new Vertex(alphaSequences.get(seq), plateAtoVMap.get(seq));
-            alphaVertices.add(alphaVertex);
+        for (Map.Entry<String, SequenceRecord> entry: alphaSequences.entrySet()) {
+            alphaVertices.add(new Vertex(entry.getValue(), vertexLabelValue));
+            vertexLabelValue++;
         }
-        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
-        Collections.sort(alphaVertices);
-        //Add ordered list of vertices to the graph
-        graphGenerator.first(alphaVertices);
-        //the list of beta vertices
+        alphaVertices.forEach(graph::addVertex);
+        //add beta sequence vertices
         List<Vertex> betaVertices = new ArrayList<>();
-        for (String seq : plateBtoVMap.keySet()) {
-            Vertex betaVertex = new Vertex(betaSequences.get(seq), plateBtoVMap.get(seq));
-            betaVertices.add(betaVertex);
+        for (Map.Entry<String, SequenceRecord> entry: betaSequences.entrySet()) {
+            betaVertices.add(new Vertex(entry.getValue(), vertexLabelValue));
+            vertexLabelValue++;
+        }
+        betaVertices.forEach(graph::addVertex);
+        //add edges (best so far)
+        int edgesAddedCount = 0;
+        for(Vertex a: alphaVertices) {
+            Set<Integer> a_wells = a.getRecord().getWells();
+            for(Vertex b: betaVertices) {
+                Set<Integer> sharedWells = new HashSet<>(a_wells);
+                sharedWells.retainAll(b.getRecord().getWells());
+                if (!sharedWells.isEmpty()) {
+                    Graphs.addEdge(graph, a, b, (double) sharedWells.size());
+                }
+                edgesAddedCount++;
+                if (edgesAddedCount % 10000000 == 0) { //collect garbage every 10,000,000 edges
+                    System.out.println(edgesAddedCount + " edges added");
+                    //request garbage collection
+                    System.gc();
+                    System.out.println("Garbage collection requested");
+                }
+            }
         }
-        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
-        Collections.sort(betaVertices);
-        //Add ordered list of vertices to the graph
-        graphGenerator.second(betaVertices);
-        //use adjacency matrix of weight created previously
-        graphGenerator.weights(weights);
-        graphGenerator.generateGraph(graph);
         if(verbose){System.out.println("Graph created");}
         //stop timing
         Instant stop = Instant.now();
@@ -145,7 +190,7 @@ public class Simulator implements GraphModificationFunctions {
                                             Integer minOverlapPercent, boolean verbose, boolean calculatePValue) {
         Instant start = Instant.now();
         SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph = data.getGraph();
-        Map<Vertex[], Integer> removedEdges = new HashMap<>();
+        Map<DefaultWeightedEdge, Vertex[]> removedEdges = new HashMap<>();
         boolean saveEdges = BiGpairSEQ.cacheGraph();
         int numWells = data.getNumWells();
         //Integer alphaCount = data.getAlphaCount();
@@ -163,6 +208,7 @@ public class Simulator implements GraphModificationFunctions {
         }
         Integer graphAlphaCount = alphas.size();
         Integer graphBetaCount = betas.size();
+        Integer graphEdgeCount = graph.edgeSet().size();
 
         //remove edges with weights outside given overlap thresholds, add those to removed edge list
         if(verbose){System.out.println("Eliminating edges with weights outside overlap threshold values");}
@@ -182,12 +228,14 @@ public class Simulator implements GraphModificationFunctions {
         if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
                 "removed");}
 
+        Integer filteredGraphEdgeCount = graph.edgeSet().size();
+
         //Find Maximum Weight Matching
         if(verbose){System.out.println("Finding maximum weight matching");}
         //The matching object
         MatchingAlgorithm<Vertex, DefaultWeightedEdge> maxWeightMatching;
         //Determine algorithm type
-        AlgorithmType algorithm = BiGpairSEQ.getMatchingAlgoritmType();
+        AlgorithmType algorithm = BiGpairSEQ.getMatchingAlgorithmType();
         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
@@ -333,8 +381,10 @@ public class Simulator implements GraphModificationFunctions {
         metadata.put("real sequence collision rate", data.getRealSequenceCollisionRate().toString());
         metadata.put("total alphas read from plate", data.getAlphaCount().toString());
         metadata.put("total betas read from plate", data.getBetaCount().toString());
+        metadata.put("initial edges in graph", graphEdgeCount.toString());
         metadata.put("alphas in graph (after pre-filtering)", graphAlphaCount.toString());
         metadata.put("betas in graph (after pre-filtering)", graphBetaCount.toString());
+        metadata.put("final edges in graph (after pre-filtering)", filteredGraphEdgeCount.toString());
         metadata.put("high overlap threshold for pairing", highThreshold.toString());
         metadata.put("low overlap threshold for pairing", lowThreshold.toString());
         metadata.put("minimum overlap percent for pairing", minOverlapPercent.toString());