Even more efficient graph creation (my initial scheme, but this time without accidentally changing what's in the sequence records)

.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="module-output" name="BiGpairSEQ_Sim" />
-      <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/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.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$/org/apache/commons/commons-rng-sampling/1.6/commons-rng-sampling-1.6.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.14.0/commons-csv-1.14.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/26.0.2/annotations-26.0.2.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="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/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.13/jheaps-0.13.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-io/1.5.2/jgrapht-io-1.5.2.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/jgrapht/jgrapht-core/1.5.2/jgrapht-core-1.5.2.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.13/jheaps-0.13.jar!/" />
+      <root url="jar://$MAVEN_REPOSITORY$/org/jheaps/jheaps/0.14/jheaps-0.14.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/apfloat/apfloat/1.10.1/apfloat-1.10.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/antlr/antlr4-runtime/4.12.0/antlr4-runtime-4.12.0.jar!/" />
-      <root url="jar://$MAVEN_REPOSITORY$/org/apache/commons/commons-lang3/3.9/commons-lang3-3.9.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)
+   4. [INPUT/OUTPUT](#input-output)
-      1. Cell Sample Files
+      1. [Cell Sample Files](#cell-sample-files)
-      2. Sample Plate Files
+      2. [Sample Plate Files](#sample-plate-files)
-      3. Graph/Data Files
+      3. [Graph/Data Files](#graph-data-files)
-      4. Matching Results Files
+      4. [Matching Results Files](#matching-results-files)
 5. [RESULTS](#results)
-   1. SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL
+   1. [SAMPLE PLATES WITH VARYING NUMBERS OF CELLS PER WELL](#sample-plates-with-varying-numbers-of-cells-per-well)
-   2. SIMULATING EXPERIMENTS FROM pairSEQ PAPER
+   2. [SIMULATING EXPERIMENTS FROM THE 2015 pairSEQ PAPER](#simulating-experiments-from-the-2015-pairseq-paper)
-6. [TODO](#todo)
+      1. [EXPERIMENT 1](#experiment-1)
-7. [CITATIONS](#citations)
+      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,67 +601,13 @@ underlying frequency distribution drastically affect the results. The real distr
 than the simulated exponential distribution. Implementing a way to exert finer control over the sampling distribution from 
 the file of distinct cells may enable better simulated replication of this experiment.
   
-## TODO
-
-* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
-* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
-  * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
-  * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
-  * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific.
-* ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
-  * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
-  * ~~Add controllable heap-type parameter?~~
-    * Parameter implemented. Fibonacci heap the current default.
-* ~~Implement sample plates with random numbers of T cells per well.~~ DONE
-  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
-    * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
-* ~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows.~~
-  * ~~Problem is variable number of cells in a well~~
-  * ~~Apache Commons CSV library writes entries a row at a time~~ 
-    * Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
-* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
-  * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
-* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
-* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
-  * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
-  * This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
-* ~~Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.~~ DONE
-  * Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
-* ~~Add read depth simulation options to CLI~~ DONE
-* ~~Update graphml output to reflect current Vertex class attributes~~ DONE
-  * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
-* ~~Implement simulation of sequences being misread as other real sequence~~ DONE
-* Implement redistributive heap for LEDA matching algorithm to achieve theoretical worst case of O(n(m + n log C)) where C is highest edge weight.
-* Update matching metadata output options in CLI
-* Add frequency distribution details to metadata output
-  * need to make an enum for the different distribution types and refactor the Plate class and user interfaces, also add the necessary fields to GraphWithMapData and then call if from Simulator
-* Update performance data in this readme
-* Add section to ReadMe describing data filtering methods.
-* 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?
-    * This would be fun and valuable, but probably take more time than I have for a hobby project.
-* ~~Implement an auction algorithm for maximum weight matching~~ DONE
-* Implement a forward/reverse auction algorithm for maximum weight matching
-* Implement an algorithm for approximating a maximum weight matching
-  * Some of these run in linear or near-linear time
-  * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
-* Implement Vose's alias method for arbitrary statistical distributions of cells
-  * Should probably refactor to use apache commons rng for this
-* Use commons JCS for caching
-* Parameterize pre-filtering options
-  
-  
 ## 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)
+* 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
+* 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
@@ -663,14 +620,68 @@ BiGpairSEQ was conceived in collaboration with the author's spouse, Dr. Alice Ma
 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.
+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 has been 
+of the published pairSEQ protocol while she was interviewing for a position at Adaptive Biotechnologies. She was 
-employed at Adaptive Biotechnologies since 2022.
+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 BiGpairSEQ
+paper and other academic publications. He continues to work on the BiGpairSEQ simulator recreationally, as it has been  
-recreationally, as it involves some very beautiful math.
+a means of exploring some very beautiful math.
+
+## TODO
+
+* Update CLI option text in this readme to include Zipf distribution options
+* ~~Try invoking GC at end of workloads to reduce paging to disk~~ DONE
+* ~~Hold graph data in memory until another graph is read-in? ABANDONED UNABANDONED~~ DONE
+    * ~~*No, this won't work, because BiGpairSEQ simulations alter the underlying graph based on filtering constraints. Changes would cascade with multiple experiments.*~~
+    * Might have figured out a way to do it, by taking edges out and then putting them back into the graph. This may actually be possible.
+    * It is possible, though the modifications to the graph incur their own performance penalties. Need testing to see which option is best. It may be computer-specific.
+* ~~Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm~~ DONE
+    * ~~in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage~~
+    * ~~Add controllable heap-type parameter?~~
+        * Parameter implemented. Pairing heap the current default.
+* ~~Implement sample plates with random numbers of T cells per well.~~ DONE
+    * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not due to nature of probability calculations.
+        * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
+* ~~See if there's a reasonable way to reformat Sample Plate files so that wells are columns instead of rows.~~
+    * ~~Problem is variable number of cells in a well~~
+    * ~~Apache Commons CSV library writes entries a row at a time~~
+        * Got this working, but at the cost of a profoundly strange bug in graph occupancy filtering. Have reverted the repo until I can figure out what caused that. Given how easily Thingiverse transposes CSV matrices in R, might not even be worth fixing.
+* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
+    * ~~Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.~~
+* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
+* ~~Implement custom Vertex class to simplify code and make it easier to implement different MWM algorithms~~ DONE
+    * Advantage: would eliminate the need to use maps to associate vertices with sequences, which would make the code easier to understand.
+    * This also seems to be faster when using the same algorithm than the version with lots of maps, which is a nice bonus!
+* ~~Implement simulation of read depth, and of read errors. Pre-filter graph for difference in read count to eliminate spurious sequences.~~ DONE
+    * Pre-filtering based on comparing (read depth) * (occupancy) to (read count) for each sequence works extremely well
+* ~~Add read depth simulation options to CLI~~ DONE
+* ~~Update graphml output to reflect current Vertex class attributes~~ DONE
+    * Individual well data from the SequenceRecords could be included, if there's ever a reason for it
+* ~~Implement simulation of sequences being misread as other real sequence~~ DONE
+* Implement redistributive heap for LEDA matching algorithm to achieve theoretical worst case of O(n(m + n log C)) where C is highest edge weight.
+* Update matching metadata output options in CLI
+* Add frequency distribution details to metadata output
+    * need to make an enum for the different distribution types and refactor the Plate class and user interfaces, also add the necessary fields to GraphWithMapData and then call if from Simulator
+* Update performance data in this readme
+* ~~Add section to ReadMe describing data filtering methods.~~ DONE, now part of algorithm description
+* Re-implement CDR1 matching method
+* ~~Refactor simulator code to collect all needed data in a single scan of the plate~~ DONE
+    * ~~Currently it scans once for the vertices and then again for the edge weights. This made simulating read depth awkward, and incompatible with caching of plate files.~~
+    * ~~This would be a fairly major rewrite of the simulator code, but could make things faster, and would definitely make them cleaner.~~
+* Implement Duan and Su's maximum weight matching algorithm
+    * ~~Add controllable algorithm-type parameter?~~ DONE
+    * This would be fun and valuable, but probably take more time than I have for a hobby project.
+* ~~Implement an auction algorithm for maximum weight matching~~ DONE
+* Implement a forward/reverse auction algorithm for maximum weight matching
+* Implement an algorithm for approximating a maximum weight matching
+    * Some of these run in linear or near-linear time
+    * given that the underlying biological samples have many, many sources of error, this would probably be the most useful option in practice. It seems less mathematically elegant, though, and so less fun for me.
+* Implement Vose's alias method for arbitrary statistical distributions of cells
+    * Should probably refactor to use apache commons rng for this
+* Use commons JCS for caching
+* Parameterize pre-filtering options

src/main/java/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.*;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
 
 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) {
+                                                              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[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
-                    Vertex target = graph.getEdgeTarget(e);
+                    removedEdges.put(e, vertices);
-                    Integer weight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, weight);
-                }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
                 }
+                edgesToRemove.add(e);
             }
         }
-        if(saveEdges) {
+        edgesToRemove.forEach(graph::removeEdge);
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
-            }
-        }
         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[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
-                    Vertex target = graph.getEdgeTarget(e);
+                    removedEdges.put(e, vertices);
-                    Integer weight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, weight);
-                }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
                 }
+                edgesToRemove.add(e);
             }
         }
-        if(saveEdges) {
+        edgesToRemove.forEach(graph::removeEdge);
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
-            }
-        }
         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[] vertices = {graph.getEdgeSource(e), graph.getEdgeTarget(e)};
-                    Vertex target = graph.getEdgeTarget(e);
+                    removedEdges.put(e, vertices);
-                    Integer intWeight = (int) graph.getEdgeWeight(e);
-                    Vertex[] edge = {source, target};
-                    removedEdges.put(edge, intWeight);
-                }
-                else {
-                    graph.setEdgeWeight(e, 0.0);
                 }
+                edgesToRemove.add(e);
             }
         }
-        if(saveEdges) {
+        edgesToRemove.forEach(graph::removeEdge);
-            for (Vertex[] edge : removedEdges.keySet()) {
-                graph.removeEdge(edge[0], edge[1]);
-            }
-        }
         return removedEdges;
     }
 
@@ -126,10 +99,10 @@ public interface GraphModificationFunctions {
     }
 
     static void addRemovedEdges(SimpleWeightedGraph<Vertex, DefaultWeightedEdge> graph,
-                                Map<Vertex[], Integer> removedEdges) {
+                                Map<DefaultWeightedEdge, Vertex[]> removedEdges) {
-        for (Vertex[] edge : removedEdges.keySet()) {
+        for (DefaultWeightedEdge edge : removedEdges.keySet()) {
-            DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);
+            Vertex[] vertices = removedEdges.get(edge);
-            graph.setEdgeWeight(e, 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("4) Zipf");
-//            System.out.println("(lambda value approximated from slope of log-log graph in figure 4c)");
+
             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();
+                    parameter = sc.nextDouble();
-                    if (stdDev <= 0.0) {
+                    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();
+                    parameter = sc.nextDouble();
-                    if (lambda <= 0.0) {
+                    if (parameter <= 0.0) {
-                        lambda = 0.6;
+                        parameter = 1.4;
-                        System.out.println("Value must be positive. Defaulting to 0.6.");
+                        System.out.println("Value must be positive. Defaulting to 1.4.");
+                    }
+                }
+                case 4 -> {
+                    BiGpairSEQ.setDistributionType(DistributionType.ZIPF);
+                    System.out.print("Please enter exponent value for Zipf distribution: ");
+                    parameter = sc.nextDouble();
+                    if (parameter <= 0.0) {
+                        parameter = 1.4;
+                        System.out.println("Value must be positive. Defaulting to 1.4.");
                     }
                 }
                 default -> {
                     System.out.println("Invalid input. Defaulting to exponential.");
-                    exponential = true;
+                    parameter = 1.4;
+                    BiGpairSEQ.setDistributionType(DistributionType.EXPONENTIAL);
                 }
             }
             System.out.print("\nNumber of wells on plate: ");
@@ -226,16 +237,17 @@ public class InteractiveInterface {
         assert filename != null;
         Plate samplePlate;
         PlateFileWriter writer;
-        if(exponential){
+        DistributionType type = BiGpairSEQ.getDistributionType();
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
+        switch(type) {
-            writer = new PlateFileWriter(filename, samplePlate);
+            case POISSON -> {
-        }
+                parameter = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
-        else {
+                samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, parameter);
-            if (poisson) {
+                writer = new PlateFileWriter(filename, samplePlate);
-                stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
+            }
+            default -> {
+                samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, parameter);
+                writer = new PlateFileWriter(filename, samplePlate);
             }
-            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
-            writer = new PlateFileWriter(filename, samplePlate);
         }
         System.out.println("Writing Sample Plate to file");
         writer.writePlateFile();
@@ -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;
+        this.stdDev = parameter;
-        if (this.exponential) {
+        this.lambda = parameter;
-            this.lambda = stdDev_or_lambda;
+        this.zipfExponent = parameter;
-            fillWellsExponential(cells.getCells(), this.lambda);
+        this.distributionType = BiGpairSEQ.getDistributionType();
-        }
+        fillWells(cells.getCells());
-        else {
-            this.stdDev = stdDev_or_lambda;
-            fillWells(cells.getCells(), this.stdDev);
-        }
     }
 
 
@@ -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();
+        this.distributionType = plate.getDistributionType();
-        if(isExponential) {
+        switch(distributionType) {
-            this.lambda = plate.getLambda();
+            case POISSON, GAUSSIAN -> {
-        }
+                this.stdDev = plate.getStdDev();
-        else{
+            }
-            this.stdDev = plate.getStdDev();
+            case EXPONENTIAL -> {
+                this.lambda = plate.getLambda();
+            }
+            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){
+            switch (distributionType) {
-                printer.printComment("Lambda: " + lambda);
+                case POISSON -> {
-            }
+                    printer.printComment("Cell frequency distribution: POISSON");
-            else {
+                }
-                printer.printComment("Std. dev.: " + stdDev);
+                case GAUSSIAN -> {
+                    printer.printComment("Cell frequency distribution: GAUSSIAN");
+                    printer.printComment("--Standard deviation: " + stdDev);
+                }
+                case EXPONENTIAL -> {
+                    printer.printComment("Cell frequency distribution: EXPONENTIAL");
+                    printer.printComment("--Lambda: " + lambda);
+                }
+                case ZIPF -> {
+                    printer.printComment("Cell frequency distribution: ZIPF");
+                    printer.printComment("--Exponent: " + zipfExponent);
+                }
             }
             printer.printRecords(wellsAsStrings);
         } catch(IOException ex){

src/main/java/Simulator.java

@@ -1,9 +1,7 @@
 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 +68,102 @@ 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");}
+         * The commented out code below works beautifully for small enough graphs. However, after implementing a
-        //For the SimpleWeightedBipartiteGraphMatrixGenerator, all vertices must have
+         * Zipf distribution and attempting to simulate Experiment 3 from the paper again, I discovered that
-        //distinct numbers associated with them. Since I'm using a 2D array, that means
+         * this method uses too much memory. Even a 120GB heap is not enough to build this adjacency matrix.
-        //distinct indices between the rows and columns. vertexStartValue lets me track where I switch
+         * So I'm going to attempt to build this graph directly and see if that is less memory intensive
-        //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
+//        //construct the graph. For simplicity, going to make
-        int vertexStartValue = 0;
+//        if(verbose){System.out.println("Making vertex maps");}
-        //keys are sequential integer vertices, values are alphas
+//        //For the SimpleWeightedBipartiteGraphMatrixGenerator, all vertices must have
-        Map<String, Integer> plateAtoVMap = makeSequenceToVertexMap(alphaSequences, vertexStartValue);
+//        //distinct numbers associated with them. Since I'm using a 2D array, that means
-        //new start value for vertex to beta map should be one more than final vertex value in alpha map
+//        //distinct indices between the rows and columns. vertexStartValue lets me track where I switch
-        vertexStartValue += plateAtoVMap.size();
+//        //from numbering rows to columns, so I can assign unique numbers to every vertex, and then
-        //keys are betas, values are sequential integers
+//        //subtract the vertexStartValue from betas to use their vertex labels as array indices
-        Map<String, Integer> plateBtoVMap = makeSequenceToVertexMap(betaSequences, vertexStartValue);
+//        int vertexStartValue = 0;
-        if(verbose){System.out.println("Vertex maps made");}
+//        //keys are sequential integer vertices, values are alphas
-        //make adjacency matrix for bipartite graph generator
+//        Map<String, Integer> plateAtoVMap = makeSequenceToVertexMap(alphaSequences, vertexStartValue);
-        //(technically this is only 1/4 of an adjacency matrix, but that's all you need
+//        //new start value for vertex to beta map should be one more than final vertex value in alpha map
-        //for a bipartite graph, and all the SimpleWeightedBipartiteGraphMatrixGenerator class expects.)
+//        vertexStartValue += plateAtoVMap.size();
-        if(verbose){System.out.println("Making adjacency matrix");}
+//        //keys are betas, values are sequential integers
-        double[][] weights = new double[plateAtoVMap.size()][plateBtoVMap.size()];
+//        Map<String, Integer> plateBtoVMap = makeSequenceToVertexMap(betaSequences, vertexStartValue);
-        fillAdjacencyMatrix(weights, vertexStartValue, alphaSequences, betaSequences, plateAtoVMap, plateBtoVMap);
+//        if(verbose){System.out.println("Vertex maps made");}
-        if(verbose){System.out.println("Adjacency matrix 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
+        int vertexLabelValue = 0;
-        SimpleWeightedBipartiteGraphMatrixGenerator graphGenerator = new SimpleWeightedBipartiteGraphMatrixGenerator();
+        //create and add alpha sequence vertices
-        //the list of alpha vertices
         List<Vertex> alphaVertices = new ArrayList<>();
-        for (String seq : plateAtoVMap.keySet()) {
+        for (Map.Entry<String, SequenceRecord> entry: alphaSequences.entrySet()) {
-            Vertex alphaVertex = new Vertex(alphaSequences.get(seq), plateAtoVMap.get(seq));
+            alphaVertices.add(new Vertex(entry.getValue(), vertexLabelValue));
-            alphaVertices.add(alphaVertex);
+            vertexLabelValue++;
         }
-        //Sort to make sure the order of vertices in list matches the order of the adjacency matrix
+        alphaVertices.forEach(graph::addVertex);
-        Collections.sort(alphaVertices);
+        //add beta sequence vertices
-        //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()) {
+        for (Map.Entry<String, SequenceRecord> entry: betaSequences.entrySet()) {
-            Vertex betaVertex = new Vertex(betaSequences.get(seq), plateBtoVMap.get(seq));
+            betaVertices.add(new Vertex(entry.getValue(), vertexLabelValue));
-            betaVertices.add(betaVertex);
+            vertexLabelValue++;
+        }
+        betaVertices.forEach(graph::addVertex);
+        //add edges
+        for(Vertex a: alphaVertices) {
+            for(Vertex b: betaVertices) {
+                Set<Integer> sharedWells = new HashSet<>(a.getRecord().getWells());
+                sharedWells.retainAll(b.getRecord().getWells());
+                double weight = (double) sharedWells.size();
+                if (weight != 0.0) {
+                    System.out.println("Edge weight: " + weight);
+                    DefaultWeightedEdge edge = graph.addEdge(a, b);
+                    graph.setEdgeWeight(edge, weight);
+                }
+                else {
+                    System.out.println("No overlap");
+
+                }
+            }
         }
-        //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 +187,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 +205,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 +225,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 +378,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());