Update readme with CLI arguments

readme.md

@@ -12,7 +12,7 @@ Unlike pairSEQ, which calculates p-values for every TCR alpha/beta overlap and c
 against a null distribution, BiGpairSEQ does not do any statistical calculations
 directly.
 
-BiGpairSEQ creates a [simple bipartite weighted graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
+BiGpairSEQ creates a [weighted bipartite graph](https://en.wikipedia.org/wiki/Bipartite_graph) representing the sample plate.
 The distinct TCRA and TCRB sequences form the two sets of vertices. Every TCRA/TCRB pair that share a well
 are connected by an edge, with the edge weight set to the number of wells in which both sequences appear.
 (Sequences present in *all* wells are filtered out prior to creating the graph, as there is no signal in their occupancy pattern.)
@@ -29,17 +29,13 @@ Unfortunately, it's a fairly new algorithm, and not yet implemented by the graph
 So this program instead uses the Fibonacci heap-based algorithm of Fredman and Tarjan (1987), which has a worst-case
 runtime of **O(n (n log(n) + m))**. The algorithm is implemented as described in Melhorn and Näher (1999).
 
-The current version of the program uses a pairing heap instead of a Fibonacci heap for its priority queue,
-which has lower theoretical efficiency but also lower complexity overhead, and is often equivalently performant
-in practice.
-
 ## USAGE
 
 ### RUNNING THE PROGRAM
 
 [Download the current version of BiGpairSEQ_Sim.](https://gitea.ejsf.synology.me/efischer/BiGpairSEQ/releases)
 
-BiGpairSEQ_Sim is an executable .jar file. Requires Java 11 or higher. [OpenJDK 17](https://jdk.java.net/17/)
+BiGpairSEQ_Sim is an executable .jar file. Requires Java 14 or higher. [OpenJDK 17](https://jdk.java.net/17/)
 recommended.
 
 Run with the command:
@@ -52,8 +48,13 @@ For example, to run the program with 32 gigabytes of memory, use the command:
 
 `java -Xmx32G -jar BiGpairSEQ_Sim.jar`
 
-Once running, BiGpairSEQ_Sim has an interactive, menu-driven CLI for generating files and simulating TCR pairing. The
+There are a number of command line options, to allow the program to be used in shell scripts. For a full list,
-main menu looks like this:
+use the -help flag:
+
+`java -jar BiGpairSEQ_Sim.jar -help`
+
+If no command line arguments are given, BiGpairSEQ_Sim will launch with an interactive, menu-driven CLI for 
+generating files and simulating TCR pairing. The main menu looks like this:
 
 ```
 --------BiGPairSEQ SIMULATOR--------
@@ -65,21 +66,52 @@ Please select an option:
 2) Generate a sample plate of T cells
 3) Generate CDR3 alpha/beta occupancy data and overlap graph
 4) Simulate bipartite graph CDR3 alpha/beta matching (BiGpairSEQ)
+8) Options
 9) About/Acknowledgments
 0) Exit
 ```
 
-### OUTPUT
+By default, the Options menu looks like this:
+```
+--------------OPTIONS---------------
+1) Turn on cell sample file caching
+2) Turn on plate file caching
+3) Turn on graph/data file caching
+4) Turn off serialized binary graph output
+5) Turn on GraphML graph output
+6) Maximum weight matching algorithm options
+0) Return to main menu
+```
+
+
+### INPUT/OUTPUT
 
 To run the simulation, the program reads and writes 4 kinds of files:
 * Cell Sample files in CSV format
 * Sample Plate files in CSV format
-* Graph and Data files in binary object serialization format
+* Graph/Data files in binary object serialization format
 * Matching Results files in CSV format
 
-When entering filenames, it is not necessary to include the file extension (.csv or .ser). When reading or
+These files are often generated in sequence. When entering filenames, it is not necessary to include the file extension
-writing files, the program will automatically add the correct extension to any filename without one.
+(.csv or .ser). When reading or writing files, the program will automatically add the correct extension to any filename 
+without one.
 
+To save file I/O time, the most recent instance of each of these four
+files either generated or read from disk can be cached in program memory. When caching is active, subsequent uses of the 
+same data file won't need to be read in again until another file of that type is used or generated,
+or caching is turned off for that file type. The program checks whether it needs to update its cached data by comparing
+filenames as entered by the user. On encountering a new filename, the program flushes its cache and reads in the new file.
+
+(Note that cached Graph/Data files must be transformed back into their original state after a matching experiment, which
+may take some time. Whether file I/O or graph transformation takes longer for graph/data files is likely to be
+device-specific.)
+
+The program's caching behavior can be controlled in the Options menu. By default, all caching is OFF.
+
+The program can optionally output Graph/Data files in .GraphML format (.graphml) for data portability. This can be 
+turned on in the Options menu. By default, GraphML output is OFF.
+
+---
 #### Cell Sample Files
 Cell Sample files consist of any number of distinct "T cells." Every cell contains 
 four sequences: Alpha CDR3, Beta CDR3, Alpha CDR1, Beta CDR1. The sequences are represented by
@@ -97,7 +129,6 @@ Comments are preceded by `#`
 
 Structure:
 
----
     # Sample contains 1 unique CDR1 for every 4 unique CDR3s.
 | Alpha CDR3 | Beta CDR3 | Alpha CDR1 | Beta CDR1 |
 |---|---|---|---|
@@ -121,15 +152,18 @@ Options when making a Sample Plate file:
     * Standard deviation size 
   * Exponential
     * Lambda value
-      * (Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was exponential with a lambda of approximately 0.6. (Howie, et al. 2015))
+      * *(Based on the slope of the graph in Figure 4C of the pairSEQ paper, the distribution of the original experiment was approximately exponential with a lambda ~0.6. (Howie, et al. 2015))*
 * Total number of wells on the plate
-* Number of sections on plate
+* Well populations random or fixed
-* Number of T cells per well
+  * If random, minimum and maximum population sizes
-  * per section, if more than one section
+  * If fixed
+    * Number of sections on plate
+    * Number of T cells per well
+      * per section, if more than one section
 * Dropout rate
 
 Files are in CSV format. There are no header labels. Every row represents a well. 
-Every column represents an individual cell, containing four sequences, depicted as an array string:
+Every value represents an individual cell, containing four sequences, depicted as an array string:
 `[CDR3A, CDR3B, CDR1A, CDR1B]`. So a representative cell might look like this: 
 
 `[525902, 791533, -1, 866282]`
@@ -139,7 +173,6 @@ Dropout sequences are replaced with the value `-1`. Comments are preceded by `#`
 
 Structure:
 
----
 ```
 # Cell source file name:
 # Each row represents one well on the plate
@@ -155,25 +188,32 @@ Structure:
 
 ---
 
-#### Graph and Data Files
+#### Graph/Data Files
-Graph and Data files are serialized binaries of a Java object containing the weigthed bipartite graph representation of a
+Graph/Data files are serialized binaries of a Java object containing the weigthed bipartite graph representation of a
 Sample Plate, along with the necessary metadata for matching and results output. Making them requires a Cell Sample file 
 (to construct a list of correct sequence pairs for checking the accuracy of BiGpairSEQ simulations) and a 
-Sample Plate file (to construct the associated occupancy graph). These files can be several gigabytes in size.
+Sample Plate file (to construct the associated occupancy graph).
-Writing them to a file lets us generate a graph and its metadata once, then use it for multiple different BiGpairSEQ simulations.
 
-Options for creating a Graph and Data file:
+These files can be several gigabytes in size. Writing them to a file lets us generate a graph and its metadata once,
+then use it for multiple different BiGpairSEQ simulations.
+
+Options for creating a Graph/Data file:
 * The Cell Sample file to use
 * The Sample Plate file to use. (This must have been generated from the selected Cell Sample file.)
 
-These files do not have a human-readable structure, and are not portable to other programs. (Export of graphs in a 
+These files do not have a human-readable structure, and are not portable to other programs.
-portable data format may be implemented in the future. The tricky part is encoding the necessary metadata.)
+
+(For portability to other software, turn on GraphML output in the Options menu. This will produce a .graphml file
+for the weighted graph, with vertex attributes sequence, type, and occupancy data.)
 
 ---
 
 #### Matching Results Files 
-Matching results files consist of the results of a BiGpairSEQ matching simulation.
+Matching results files consist of the results of a BiGpairSEQ matching simulation. Making them requires a serialized
-Files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details.
+binary Graph/Data file (.ser). (Because .graphML files are larger than .ser files, BiGpairSEQ_Sim supports .graphML
+output only. Graph/data input must use a serialized binary.)
+
+Matching results files are in CSV format. Rows are sequence pairings with extra relevant data. Columns are pairing-specific details.
 Metadata about the matching simulation is included as comments. Comments are preceded by `#`.
 
 Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
@@ -188,7 +228,6 @@ Options when running a BiGpairSEQ simulation of CDR3 alpha/beta matching:
 
 Example output:
 
----
 ```
 # Source Sample Plate file: 4MilCellsPlate.csv
 # Source Graph and Data file: 4MilCellsPlateGraph.ser
@@ -239,26 +278,30 @@ slightly less time than the simulation itself. Real elapsed time from start to f
 ## 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
+* ~~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. If so, awesome.
+  * 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.
+* ~~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.
+    * 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._
+    * _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.
-* Re-implement command line arguments, to enable scripting and statistical simulation studies
+* ~~Enable GraphML output in addition to serialized object binaries, for data portability~~ DONE
-* Implement sample plates with random numbers of T cells per well.
+  * ~~Custom vertex type with attribute for sequence occupancy?~~ ABANDONED
-  * Possible BiGpairSEQ advantage over pairSEQ: BiGpairSEQ is resilient to variations in well population sizes on a sample plate; pairSEQ is not.
+    * Have a branch where this is implemented, but there's a bug that broke matching. Don't currently have time to fix.
-    * preliminary data suggests that BiGpairSEQ behaves roughly as though the whole plate had whatever the *average* well concentration is, but that's still speculative.
+* ~~Re-implement command line arguments, to enable scripting and statistical simulation studies~~ DONE
-* Enable GraphML output in addition to serialized object binaries, for data portability
-  * Custom vertex type with attribute for sequence occupancy?
 * Re-implement CDR1 matching method
 * Implement Duan and Su's maximum weight matching algorithm
-  * Add controllable algorithm-type parameter?
+    * Add controllable algorithm-type parameter?
-* Test whether pairing heap (currently used) or Fibonacci heap is more efficient for priority queue in current matching algorithm
+    * This would be fun and valuable, but probably take more time than I have for a hobby project.
-  * in theory Fibonacci heap should be more efficient, but complexity overhead may eliminate theoretical advantage
+* Implement Vose's alias method for arbitrary statistical distributions of cells
-  * Add controllable heap-type parameter?
-
   
   
 ## CITATIONS

src/main/java/BiGpairSEQ.java

@@ -1,8 +1,21 @@
-//main class. Only job is to choose which interface to use, and hold graph data in memory
+import java.util.Random;
+
+//main class. For choosing interface type and holding settings
 public class BiGpairSEQ {
 
+    private static final Random rand = new Random();
+    private static CellSample cellSampleInMemory = null;
+    private static String cellFilename = null;
+    private static Plate plateInMemory = null;
+    private static String plateFilename = null;
     private static GraphWithMapData graphInMemory = null;
     private static String graphFilename = null;
+    private static boolean cacheCells = false;
+    private static boolean cachePlate = false;
+    private static boolean cacheGraph = false;
+    private static String priorityQueueHeapType = "FIBONACCI";
+    private static boolean outputBinary = true;
+    private static boolean outputGraphML = false;
 
     public static void main(String[] args) {
         if (args.length == 0) {
@@ -10,33 +23,154 @@ public class BiGpairSEQ {
         }
         else {
             //This will be uncommented when command line arguments are re-implemented.
-            //CommandLineInterface.startCLI(args);
+            CommandLineInterface.startCLI(args);
-            System.out.println("Command line arguments are still being re-implemented.");
+            //System.out.println("Command line arguments are still being re-implemented.");
         }
     }
 
-    public static GraphWithMapData getGraph() {
+    public static Random getRand() {
-        return graphInMemory;
+        return rand;
     }
 
-    public static void setGraph(GraphWithMapData g) {
+    public static CellSample getCellSampleInMemory() {
+        return cellSampleInMemory;
+    }
+
+    public static void setCellSampleInMemory(CellSample cellSample, String filename) {
+        if(cellSampleInMemory != null) {
+            clearCellSampleInMemory();
+        }
+        cellSampleInMemory = cellSample;
+        cellFilename = filename;
+        System.out.println("Cell sample file " + filename + " cached.");
+    }
+
+    public static void clearCellSampleInMemory() {
+        cellSampleInMemory = null;
+        cellFilename = null;
+        System.gc();
+        System.out.println("Cell sample file cache cleared.");
+
+    }
+
+    public static String getCellFilename() {
+        return cellFilename;
+    }
+
+    public static Plate getPlateInMemory() {
+        return plateInMemory;
+    }
+
+    public static void setPlateInMemory(Plate plate, String filename) {
+        if(plateInMemory != null) {
+            clearPlateInMemory();
+        }
+        plateInMemory = plate;
+        plateFilename = filename;
+        System.out.println("Sample plate file " + filename + " cached.");
+    }
+
+    public static void clearPlateInMemory() {
+        plateInMemory = null;
+        plateFilename = null;
+        System.gc();
+        System.out.println("Sample plate file cache cleared.");
+
+    }
+
+    public static String getPlateFilename() {
+        return plateFilename;
+    }
+
+
+    public static GraphWithMapData getGraphInMemory() {return graphInMemory;
+    }
+
+    public static void setGraphInMemory(GraphWithMapData g, String filename) {
         if (graphInMemory != null) {
-            clearGraph();
+            clearGraphInMemory();
         }
         graphInMemory = g;
+        graphFilename = filename;
+        System.out.println("Graph and data file " + filename + " cached.");
     }
 
-    public static void clearGraph() {
+    public static void clearGraphInMemory() {
         graphInMemory = null;
+        graphFilename = null;
         System.gc();
+        System.out.println("Graph and data file cache cleared.");
     }
 
     public static String getGraphFilename() {
         return graphFilename;
     }
 
-    public static void setGraphFilename(String filename) {
+
-        graphFilename = filename;
+    public static boolean cacheCells() {
+        return cacheCells;
     }
 
+    public static void setCacheCells(boolean cacheCells) {
+        //if not caching, clear the memory
+        if(!cacheCells){
+            BiGpairSEQ.clearCellSampleInMemory();
+            System.out.println("Cell sample file caching: OFF.");
+        }
+        else {
+            System.out.println("Cell sample file caching: ON.");
+        }
+        BiGpairSEQ.cacheCells = cacheCells;
+    }
+
+    public static boolean cachePlate() {
+        return cachePlate;
+    }
+
+    public static void setCachePlate(boolean cachePlate) {
+        //if not caching, clear the memory
+        if(!cachePlate) {
+            BiGpairSEQ.clearPlateInMemory();
+            System.out.println("Sample plate file caching: OFF.");
+        }
+        else {
+            System.out.println("Sample plate file caching: ON.");
+        }
+        BiGpairSEQ.cachePlate = cachePlate;
+    }
+
+    public static boolean cacheGraph() {
+        return cacheGraph;
+    }
+
+    public static void setCacheGraph(boolean cacheGraph) {
+        //if not caching, clear the memory
+        if(!cacheGraph) {
+            BiGpairSEQ.clearGraphInMemory();
+            System.out.println("Graph/data file caching: OFF.");
+        }
+        else {
+            System.out.println("Graph/data file caching: ON.");
+        }
+        BiGpairSEQ.cacheGraph = cacheGraph;
+    }
+
+    public static String getPriorityQueueHeapType() {
+        return priorityQueueHeapType;
+    }
+
+    public static void setPairingHeap() {
+        priorityQueueHeapType = "PAIRING";
+    }
+
+    public static void setFibonacciHeap() {
+        priorityQueueHeapType = "FIBONACCI";
+    }
+
+    public static boolean outputBinary() {return outputBinary;}
+    public static void setOutputBinary(boolean b) {outputBinary = b;}
+
+    public static boolean outputGraphML() {return outputGraphML;}
+    public static void setOutputGraphML(boolean b) {outputGraphML = b;}
+
 }

src/main/java/CellFileReader.java

@@ -13,6 +13,7 @@ public class CellFileReader {
 
     private String filename;
     private List<Integer[]> distinctCells = new ArrayList<>();
+    private Integer cdr1Freq;
 
     public CellFileReader(String filename) {
         if(!filename.matches(".*\\.csv")){
@@ -38,19 +39,37 @@ public class CellFileReader {
                 cell[3] = Integer.valueOf(record.get("Beta CDR1"));
                 distinctCells.add(cell);
             }
+
+
         } catch(IOException ex){
             System.out.println("cell file " + filename + " not found.");
             System.err.println(ex);
         }
+
+        //get CDR1 frequency
+        ArrayList<Integer> cdr1Alphas = new ArrayList<>();
+        for (Integer[] cell : distinctCells) {
+            cdr1Alphas.add(cell[3]);
+        }
+        double count = cdr1Alphas.stream().distinct().count();
+        count = Math.ceil(distinctCells.size() / count);
+        cdr1Freq = (int) count;
+
+    }
+
+    public CellSample getCellSample() {
+        return new CellSample(distinctCells, cdr1Freq);
     }
 
     public String getFilename() { return filename;}
 
-    public List<Integer[]> getCells(){
+    //Refactor everything that uses this to have access to a Cell Sample and get the cells there instead.
+    public List<Integer[]> getListOfDistinctCellsDEPRECATED(){
         return distinctCells;
     }
 
-    public Integer getCellCount() {
+    public Integer getCellCountDEPRECATED() {
+        //Refactor everything that uses this to have access to a Cell Sample and get the count there instead.
         return distinctCells.size();
     }
 }

src/main/java/CellSample.java

@@ -1,10 +1,37 @@
+import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
+import java.util.stream.IntStream;
 
 public class CellSample {
 
     private List<Integer[]> cells;
     private Integer cdr1Freq;
 
+    public CellSample(Integer numDistinctCells, Integer cdr1Freq){
+        this.cdr1Freq = cdr1Freq;
+        List<Integer> numbersCDR3 = new ArrayList<>();
+        List<Integer> numbersCDR1 = new ArrayList<>();
+        Integer numDistCDR3s = 2 * numDistinctCells + 1;
+        IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
+        IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
+        Collections.shuffle(numbersCDR3);
+        Collections.shuffle(numbersCDR1);
+
+        //Each cell represented by 4 values
+        //two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
+        List<Integer[]> distinctCells = new ArrayList<>();
+        for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
+            Integer tmpCDR3a = numbersCDR3.get(i);
+            Integer tmpCDR3b = numbersCDR3.get(i+1);
+            Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
+            Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
+            Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
+            distinctCells.add(tmp);
+        }
+        this.cells = distinctCells;
+    }
+
     public CellSample(List<Integer[]> cells, Integer cdr1Freq){
         this.cells = cells;
         this.cdr1Freq = cdr1Freq;
@@ -18,7 +45,7 @@ public class CellSample {
         return cdr1Freq;
     }
 
-    public Integer population(){
+    public Integer getCellCount(){
         return cells.size();
     }
 

src/main/java/CommandLineInterface.java

@@ -1,5 +1,9 @@
 import org.apache.commons.cli.*;
 
+import java.io.IOException;
+import java.util.Arrays;
+import java.util.stream.Stream;
+
 /*
  * Class for parsing options passed to program from command line
  *
@@ -29,6 +33,8 @@ import org.apache.commons.cli.*;
  * cellfile : name of the cell sample file to use as input
  * platefile : name of the sample plate file to use as input
  * output : name of the output file
+ * graphml : output a graphml file
+ * binary : output a serialized binary object file
  *
  * Match flags:
  * graphFile : name of graph and data file to use as input
@@ -43,286 +49,379 @@ import org.apache.commons.cli.*;
 public class CommandLineInterface {
 
     public static void startCLI(String[] args) {
-        //These command line options are a big mess
+        //Options sets for the different modes
-        //Really, I don't think command line tools are expected to work in this many different modes
+        Options mainOptions = buildMainOptions();
-        //making cells, making plates, and matching are the sort of thing that UNIX philosophy would say
+        Options cellOptions = buildCellOptions();
-        //should be three separate programs.
+        Options plateOptions = buildPlateOptions();
-        //There might be a way to do it with option parameters?
+        Options graphOptions = buildGraphOptions();
-
+        Options matchOptions = buildMatchCDR3options();
-        //main options set
-        Options mainOptions = new Options();
-        Option makeCells = Option.builder("cells")
-                .longOpt("make-cells")
-                .desc("Makes a file of distinct cells")
-                .build();
-        Option makePlate = Option.builder("plates")
-                .longOpt("make-plates")
-                .desc("Makes a sample plate file")
-                .build();
-        Option makeGraph = Option.builder("graph")
-                .longOpt("make-graph")
-                .desc("Makes a graph and data file")
-                .build();
-        Option matchCDR3 = Option.builder("match")
-                .longOpt("match-cdr3")
-                .desc("Match CDR3s. Requires a cell sample file and any number of plate files.")
-                .build();
-        OptionGroup mainGroup = new OptionGroup();
-        mainGroup.addOption(makeCells);
-        mainGroup.addOption(makePlate);
-        mainGroup.addOption(makeGraph);
-        mainGroup.addOption(matchCDR3);
-        mainGroup.setRequired(true);
-        mainOptions.addOptionGroup(mainGroup);
-
-        //Reuse clones of this for other options groups, rather than making it lots of times
-        Option outputFile = Option.builder("o")
-                .longOpt("output-file")
-                .hasArg()
-                .argName("filename")
-                .desc("Name of output file")
-                .build();
-        mainOptions.addOption(outputFile);
-
-        //Options cellOptions = new Options();
-        Option numCells = Option.builder("nc")
-                .longOpt("num-cells")
-                .desc("The number of distinct cells to generate")
-                .hasArg()
-                .argName("number")
-                .build();
-        mainOptions.addOption(numCells);
-        Option cdr1Freq = Option.builder("d")
-                .longOpt("peptide-diversity-factor")
-                .hasArg()
-                .argName("number")
-                .desc("Number of distinct CDR3s for every CDR1")
-                .build();
-        mainOptions.addOption(cdr1Freq);
-        //Option cellOutput = (Option) outputFile.clone();
-        //cellOutput.setRequired(true);
-        //mainOptions.addOption(cellOutput);
-
-        //Options plateOptions = new Options();
-        Option inputCells = Option.builder("c")
-                .longOpt("cell-file")
-                .hasArg()
-                .argName("file")
-                .desc("The cell sample file used for filling wells")
-                .build();
-        mainOptions.addOption(inputCells);
-        Option numWells = Option.builder("w")
-                .longOpt("num-wells")
-                .hasArg()
-                .argName("number")
-                .desc("The number of wells on each plate")
-                .build();
-        mainOptions.addOption(numWells);
-        Option numPlates = Option.builder("np")
-                .longOpt("num-plates")
-                .hasArg()
-                .argName("number")
-                .desc("The number of plate files to output")
-                .build();
-        mainOptions.addOption(numPlates);
-        //Option plateOutput = (Option) outputFile.clone();
-        //plateOutput.setRequired(true);
-        //plateOutput.setDescription("Prefix for plate output filenames");
-        //mainOptions.addOption(plateOutput);
-        Option plateErr = Option.builder("err")
-                .longOpt("drop-out-rate")
-                .hasArg()
-                .argName("number")
-                .desc("Well drop-out rate. (Probability between 0 and 1)")
-                .build();
-        mainOptions.addOption(plateErr);
-        Option plateConcentrations = Option.builder("t")
-                .longOpt("t-cells-per-well")
-                .hasArgs()
-                .argName("number 1, number 2, ...")
-                .desc("Number of T cells per well for each plate section")
-                .build();
-        mainOptions.addOption(plateConcentrations);
-
-//different distributions, mutually exclusive
-        OptionGroup plateDistributions = new OptionGroup();
-        Option plateExp = Option.builder("exponential")
-                .desc("Sample from distinct cells with exponential frequency distribution")
-                .build();
-        plateDistributions.addOption(plateExp);
-        Option plateGaussian = Option.builder("gaussian")
-                .desc("Sample from distinct cells with gaussain frequency distribution")
-                .build();
-        plateDistributions.addOption(plateGaussian);
-        Option platePoisson = Option.builder("poisson")
-                .desc("Sample from distinct cells with poisson frequency distribution")
-                .build();
-        plateDistributions.addOption(platePoisson);
-        mainOptions.addOptionGroup(plateDistributions);
-
-        Option plateStdDev = Option.builder("stddev")
-                .desc("Standard deviation for gaussian distribution")
-                .hasArg()
-                .argName("number")
-                .build();
-        mainOptions.addOption(plateStdDev);
-
-        Option plateLambda = Option.builder("lambda")
-                .desc("Lambda for exponential distribution")
-                .hasArg()
-                .argName("number")
-                .build();
-        mainOptions.addOption(plateLambda);
-
-
-
-//
-//            String cellFile, String filename, Double stdDev,
-//                    Integer numWells, Integer numSections,
-//                    Integer[] concentrations, Double dropOutRate
-//
-
-        //Options matchOptions = new Options();
-        inputCells.setDescription("The cell sample file to be used for matching.");
-        mainOptions.addOption(inputCells);
-        Option lowThresh = Option.builder("low")
-                .longOpt("low-threshold")
-                .hasArg()
-                .argName("number")
-                .desc("Sets the minimum occupancy overlap to attempt matching")
-                .build();
-        mainOptions.addOption(lowThresh);
-        Option highThresh = Option.builder("high")
-                .longOpt("high-threshold")
-                .hasArg()
-                .argName("number")
-                .desc("Sets the maximum occupancy overlap to attempt matching")
-                .build();
-        mainOptions.addOption(highThresh);
-        Option occDiff = Option.builder("occdiff")
-                .longOpt("occupancy-difference")
-                .hasArg()
-                .argName("Number")
-                .desc("Maximum difference in alpha/beta occupancy to attempt matching")
-                .build();
-        mainOptions.addOption(occDiff);
-        Option overlapPer = Option.builder("ovper")
-                .longOpt("overlap-percent")
-                .hasArg()
-                .argName("Percent")
-                .desc("Minimum overlap percent to attempt matching (0 -100)")
-                .build();
-        mainOptions.addOption(overlapPer);
-        Option inputPlates = Option.builder("p")
-                .longOpt("plate-files")
-                .hasArgs()
-                .desc("Plate files to match")
-                .build();
-        mainOptions.addOption(inputPlates);
-
-
 
         CommandLineParser parser = new DefaultParser();
-        try {
+        try{
-            CommandLine line = parser.parse(mainOptions, args);
+            CommandLine line = parser.parse(mainOptions, Arrays.copyOfRange(args, 0, 1));
-            if(line.hasOption("match")){
+
-                //line = parser.parse(mainOptions, args);
+            if (line.hasOption("help")) {
-                //String cellFile = line.getOptionValue("c");
+                HelpFormatter formatter = new HelpFormatter();
-                String graphFile = line.getOptionValue("g");
+                formatter.printHelp("BiGpairSEQ_Sim", mainOptions);
-                Integer lowThreshold = Integer.valueOf(line.getOptionValue(lowThresh));
+                System.out.println();
-                Integer highThreshold = Integer.valueOf(line.getOptionValue(highThresh));
+                formatter.printHelp("BiGpairSEQ_SIM -cells", cellOptions);
-                Integer occupancyDifference = Integer.valueOf(line.getOptionValue(occDiff));
+                System.out.println();
-                Integer overlapPercent = Integer.valueOf(line.getOptionValue(overlapPer));
+                formatter.printHelp("BiGpairSEQ_Sim -plate", plateOptions);
-                for(String plate: line.getOptionValues("p")) {
+                System.out.println();
-                    matchCDR3s(graphFile, lowThreshold, highThreshold, occupancyDifference, overlapPercent);
+                formatter.printHelp("BiGpairSEQ_Sim -graph", graphOptions);
-                }
+                System.out.println();
+                formatter.printHelp("BiGpairSEQ_Sim -match", matchOptions);
             }
-            else if(line.hasOption("cells")){
+            else if (line.hasOption("cells")) {
-                //line = parser.parse(mainOptions, args);
+                line = parser.parse(cellOptions, Arrays.copyOfRange(args, 1, args.length));
+                Integer number = Integer.valueOf(line.getOptionValue("n"));
+                Integer diversity = Integer.valueOf(line.getOptionValue("d"));
                 String filename = line.getOptionValue("o");
-                Integer numDistCells = Integer.valueOf(line.getOptionValue("nc"));
+                makeCells(filename, number, diversity);
-                Integer freq = Integer.valueOf(line.getOptionValue("d"));
-                makeCells(filename, numDistCells, freq);
             }
-            else if(line.hasOption("plates")){
-                //line = parser.parse(mainOptions, args);
-                String cellFile = line.getOptionValue("c");
-                String filenamePrefix = line.getOptionValue("o");
-                Integer numWellsOnPlate = Integer.valueOf(line.getOptionValue("w"));
-                Integer numPlatesToMake = Integer.valueOf(line.getOptionValue("np"));
-                String[] concentrationsToUseString = line.getOptionValues("t");
-                Integer numSections = concentrationsToUseString.length;
 
-                Integer[] concentrationsToUse = new Integer[numSections];
+            else if (line.hasOption("plate")) {
-                for(int i = 0; i <numSections; i++){
+                line = parser.parse(plateOptions, Arrays.copyOfRange(args, 1, args.length));
-                    concentrationsToUse[i] = Integer.valueOf(concentrationsToUseString[i]);
+                //get the cells
+                String cellFilename = line.getOptionValue("c");
+                CellSample cells = getCells(cellFilename);
+                //get the rest of the parameters
+                Integer[] populations;
+                String outputFilename = line.getOptionValue("o");
+                Integer numWells = Integer.parseInt(line.getOptionValue("w"));
+                Double dropoutRate = Double.parseDouble(line.getOptionValue("err"));
+                if (line.hasOption("random")) {
+                    //Array holding values of minimum and maximum populations
+                   Integer[] min_max = Stream.of(line.getOptionValues("random"))
+                           .mapToInt(Integer::parseInt)
+                           .boxed()
+                           .toArray(Integer[]::new);
+                   populations = BiGpairSEQ.getRand().ints(min_max[0], min_max[1] + 1)
+                            .limit(numWells)
+                            .boxed()
+                            .toArray(Integer[]::new);
                 }
-                Double dropOutRate = Double.valueOf(line.getOptionValue("err"));
+                else if (line.hasOption("pop")) {
-                if(line.hasOption("exponential")){
+                    populations = Stream.of(line.getOptionValues("pop"))
-                    Double lambda = Double.valueOf(line.getOptionValue("lambda"));
+                            .mapToInt(Integer::parseInt)
-                    for(int i = 1; i <= numPlatesToMake; i++){
+                            .boxed()
-                        makePlateExp(cellFile, filenamePrefix + i, lambda, numWellsOnPlate,
+                            .toArray(Integer[]::new);
-                                concentrationsToUse,dropOutRate);
-                    }
                 }
-                else if(line.hasOption("gaussian")){
+                else{
-                    Double stdDev = Double.valueOf(line.getOptionValue("std-dev"));
+                    populations = new Integer[1];
-                    for(int i = 1; i <= numPlatesToMake; i++){
+                    populations[0] = 1;
-                        makePlate(cellFile, filenamePrefix + i, stdDev, numWellsOnPlate,
+                }
-                                concentrationsToUse,dropOutRate);
+                //make the plate
-                    }
+                Plate plate;
+                if (line.hasOption("poisson")) {
+                    Double stdDev = Math.sqrt(numWells);
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
+                }
+                else if (line.hasOption("gaussian")) {
+                    Double stdDev = Double.parseDouble(line.getOptionValue("stddev"));
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, stdDev, false);
+                }
+                else {
+                    assert line.hasOption("exponential");
+                    Double lambda = Double.parseDouble(line.getOptionValue("lambda"));
+                    plate = new Plate(cells, cellFilename, numWells, populations, dropoutRate, lambda, true);
+                }
+                PlateFileWriter writer = new PlateFileWriter(outputFilename, plate);
+                writer.writePlateFile();
+            }
 
+            else if (line.hasOption("graph")) { //Making a graph
+                line = parser.parse(graphOptions, Arrays.copyOfRange(args, 1, args.length));
+                String cellFilename = line.getOptionValue("c");
+                String plateFilename = line.getOptionValue("p");
+                String outputFilename = line.getOptionValue("o");
+                //get cells
+                CellSample cells = getCells(cellFilename);
+                //get plate
+                Plate plate = getPlate(plateFilename);
+                GraphWithMapData graph = Simulator.makeGraph(cells, plate, false);
+                if (!line.hasOption("no-binary")) { //output binary file unless told not to
+                    GraphDataObjectWriter writer = new GraphDataObjectWriter(outputFilename, graph, false);
+                    writer.writeDataToFile();
                 }
-                else if(line.hasOption("poisson")){
+                if (line.hasOption("graphml")) { //if told to, output graphml file
-                    for(int i = 1; i <= numPlatesToMake; i++){
+                    GraphMLFileWriter gmlwriter = new GraphMLFileWriter(outputFilename, graph);
-                        makePlatePoisson(cellFile, filenamePrefix + i, numWellsOnPlate,
+                    gmlwriter.writeGraphToFile();
-                                concentrationsToUse,dropOutRate);
-                    }
                 }
             }
+
+            else if (line.hasOption("match")) { //can add a flag for which match type in future, spit this in two
+                line = parser.parse(matchOptions, Arrays.copyOfRange(args, 1, args.length));
+                String graphFilename = line.getOptionValue("g");
+                String outputFilename = line.getOptionValue("o");
+                Integer minThreshold = Integer.parseInt(line.getOptionValue("min"));
+                Integer maxThreshold = Integer.parseInt(line.getOptionValue("max"));
+                Integer minOverlapPct;
+                if (line.hasOption("minpct")) { //see if this filter is being used
+                    minOverlapPct = Integer.parseInt(line.getOptionValue("minpct"));
+                }
+                else {
+                    minOverlapPct = 0;
+                }
+                Integer maxOccupancyDiff;
+                if (line.hasOption("maxdiff")) { //see if this filter is being used
+                    maxOccupancyDiff = Integer.parseInt(line.getOptionValue("maxdiff"));
+                }
+                else {
+                    maxOccupancyDiff = Integer.MAX_VALUE;
+                }
+                GraphWithMapData graph = getGraph(graphFilename);
+                MatchingResult result = Simulator.matchCDR3s(graph, graphFilename, minThreshold, maxThreshold,
+                        maxOccupancyDiff, minOverlapPct, false);
+                MatchingFileWriter writer = new MatchingFileWriter(outputFilename, result);
+                writer.writeResultsToFile();
+                //can put a bunch of ifs for outputting various things from the MatchingResult to System.out here
+                //after I put those flags in the matchOptions
+            }
         }
         catch (ParseException exp) {
             System.err.println("Parsing failed.  Reason: " + exp.getMessage());
         }
     }
 
+    private static Option outputFileOption() {
+        Option outputFile = Option.builder("o")
+                .longOpt("output-file")
+                .hasArg()
+                .argName("filename")
+                .desc("Name of output file")
+                .required()
+                .build();
+        return outputFile;
+    }
+
+    private static Options buildMainOptions() {
+        Options mainOptions = new Options();
+        Option help = Option.builder("help")
+                .desc("Displays this help menu")
+                .build();
+        Option makeCells = Option.builder("cells")
+                .longOpt("make-cells")
+                .desc("Makes a cell sample file of distinct T cells")
+                .build();
+        Option makePlate = Option.builder("plate")
+                .longOpt("make-plate")
+                .desc("Makes a sample plate file. Requires a cell sample file.")
+                .build();
+        Option makeGraph = Option.builder("graph")
+                .longOpt("make-graph")
+                .desc("Makes a graph/data file. Requires a cell sample file and a sample plate file")
+                .build();
+        Option matchCDR3 = Option.builder("match")
+                .longOpt("match-cdr3")
+                .desc("Matches CDR3s. Requires a graph/data file.")
+                .build();
+        OptionGroup mainGroup = new OptionGroup();
+        mainGroup.addOption(help);
+        mainGroup.addOption(makeCells);
+        mainGroup.addOption(makePlate);
+        mainGroup.addOption(makeGraph);
+        mainGroup.addOption(matchCDR3);
+        mainGroup.setRequired(true);
+        mainOptions.addOptionGroup(mainGroup);
+        return mainOptions;
+    }
+
+    private static Options buildCellOptions() {
+        Options cellOptions = new Options();
+        Option numCells = Option.builder("n")
+                .longOpt("num-cells")
+                .desc("The number of distinct cells to generate")
+                .hasArg()
+                .argName("number")
+                .required().build();
+        Option cdr3Diversity = Option.builder("d")
+                .longOpt("diversity-factor")
+                .desc("The factor by which unique CDR3s outnumber unique CDR1s")
+                .hasArg()
+                .argName("factor")
+                .required().build();
+        cellOptions.addOption(numCells);
+        cellOptions.addOption(cdr3Diversity);
+        cellOptions.addOption(outputFileOption());
+        return cellOptions;
+    }
+
+    private static Options buildPlateOptions() {
+        Options plateOptions = new Options();
+        Option cellFile = Option.builder("c") // add this to plate options
+                .longOpt("cell-file")
+                .desc("The cell sample file to use")
+                .hasArg()
+                .argName("filename")
+                .required().build();
+        Option numWells = Option.builder("w")// add this to plate options
+                .longOpt("wells")
+                .desc("The number of wells on the sample plate")
+                .hasArg()
+                .argName("number")
+                .required().build();
+        //options group for choosing with distribution to use
+        OptionGroup distributions = new OptionGroup();// add this to plate options
+        distributions.setRequired(true);
+        Option poisson = Option.builder("poisson")
+                .desc("Use a Poisson distribution for cell sample")
+                .build();
+        Option gaussian = Option.builder("gaussian")
+                .desc("Use a Gaussian distribution for cell sample")
+                .build();
+        Option exponential = Option.builder("exponential")
+                .desc("Use an exponential distribution for cell sample")
+                .build();
+        distributions.addOption(poisson);
+        distributions.addOption(gaussian);
+        distributions.addOption(exponential);
+        //options group for statistical distribution parameters
+        OptionGroup statParams = new OptionGroup();// add this to plate options
+        Option stdDev = Option.builder("stddev")
+                .desc("If using -gaussian flag, standard deviation for distrbution")
+                .hasArg()
+                .argName("value")
+                .build();
+        Option lambda = Option.builder("lambda")
+                .desc("If using -exponential flag, lambda value for distribution")
+                .hasArg()
+                .argName("value")
+                .build();
+        statParams.addOption(stdDev);
+        statParams.addOption(lambda);
+        //Option group for random plate or set populations
+        OptionGroup wellPopOptions = new OptionGroup(); // add this to plate options
+        wellPopOptions.setRequired(true);
+        Option randomWellPopulations = Option.builder("random")
+                .desc("Randomize well populations on sample plate. Takes two arguments: the minimum possible population and the maximum possible population.")
+                .hasArgs()
+                .numberOfArgs(2)
+                .argName("minimum maximum")
+                .build();
+        Option specificWellPopulations = Option.builder("pop")
+                .desc("The well populations for each section of the sample plate. There will be as many sections as there are populations given.")
+                .hasArgs()
+                .argName("number [number]...")
+                .build();
+        Option dropoutRate = Option.builder("err") //add this to plate options
+                .hasArg()
+                .desc("The sequence dropout rate due to amplification error. (0.0 - 1.0)")
+                .argName("rate")
+                .required()
+                .build();
+        wellPopOptions.addOption(randomWellPopulations);
+        wellPopOptions.addOption(specificWellPopulations);
+        plateOptions.addOption(cellFile);
+        plateOptions.addOption(numWells);
+        plateOptions.addOptionGroup(distributions);
+        plateOptions.addOptionGroup(statParams);
+        plateOptions.addOptionGroup(wellPopOptions);
+        plateOptions.addOption(dropoutRate);
+        plateOptions.addOption(outputFileOption());
+        return plateOptions;
+    }
+
+    private static Options buildGraphOptions() {
+        Options graphOptions = new Options();
+        Option cellFilename = Option.builder("c")
+                .longOpt("cell-file")
+                .desc("Cell sample file to use for checking accuracy")
+                .hasArg()
+                .argName("filename")
+                .required().build();
+        Option plateFilename = Option.builder("p")
+                .longOpt("plate-filename")
+                .desc("Sample plate file (made from given cell sample file) to construct graph from")
+                .hasArg()
+                .argName("filename")
+                .required().build();
+        Option outputGraphML = Option.builder("graphml")
+                .desc("Output GraphML file")
+                .build();
+        Option outputSerializedBinary = Option.builder("nb")
+                .longOpt("no-binary")
+                .desc("Don't output serialized binary file")
+                .build();
+        graphOptions.addOption(cellFilename);
+        graphOptions.addOption(plateFilename);
+        graphOptions.addOption(outputFileOption());
+        graphOptions.addOption(outputGraphML);
+        graphOptions.addOption(outputSerializedBinary);
+        return graphOptions;
+    }
+
+    private static Options buildMatchCDR3options() {
+        Options matchCDR3options = new Options();
+        Option graphFilename = Option.builder("g")
+                .longOpt("graph-file")
+                .desc("The graph/data file to use")
+                .hasArg()
+                .argName("filename")
+                .required().build();
+        Option minOccupancyOverlap = Option.builder("min")
+                .desc("The minimum number of shared wells to attempt to match a sequence pair")
+                .hasArg()
+                .argName("number")
+                .required().build();
+        Option maxOccupancyOverlap = Option.builder("max")
+                .desc("The maximum number of shared wells to attempt to match a sequence pair")
+                .hasArg()
+                .argName("number")
+                .required().build();
+        Option minOverlapPercent = Option.builder("minpct")
+                .desc("(Optional) The minimum percentage of a sequence's total occupancy shared by another sequence to attempt matching. (0 - 100) ")
+                .hasArg()
+                .argName("percent")
+                .build();
+        Option maxOccupancyDifference = Option.builder("maxdiff")
+                .desc("(Optional) The maximum difference in total occupancy between two sequences to attempt matching.")
+                .hasArg()
+                .argName("number")
+                .build();
+        matchCDR3options.addOption(graphFilename);
+        matchCDR3options.addOption(minOccupancyOverlap);
+        matchCDR3options.addOption(maxOccupancyOverlap);
+        matchCDR3options.addOption(minOverlapPercent);
+        matchCDR3options.addOption(maxOccupancyDifference);
+        matchCDR3options.addOption(outputFileOption());
+        //options for output to System.out
+       //Option printPairingErrorRate = Option.builder()
+
+        return matchCDR3options;
+    }
+
+
+
+    private static CellSample getCells(String cellFilename) {
+        assert cellFilename != null;
+        CellFileReader reader = new CellFileReader(cellFilename);
+        return reader.getCellSample();
+    }
+
+    private static Plate getPlate(String plateFilename) {
+        assert plateFilename != null;
+        PlateFileReader reader = new PlateFileReader(plateFilename);
+        return reader.getSamplePlate();
+    }
+
+    private static GraphWithMapData getGraph(String graphFilename) {
+        assert graphFilename != null;
+        try{
+            GraphDataObjectReader reader = new GraphDataObjectReader(graphFilename, false);
+            return reader.getData();
+
+        }
+        catch (IOException ex) {
+            ex.printStackTrace();
+            return null;
+        }
+    }
+
     //for calling from command line
-    public static void makeCells(String filename, Integer numCells, Integer cdr1Freq){
+    public static void makeCells(String filename, Integer numCells, Integer cdr1Freq) {
-        CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
+        CellSample sample = new CellSample(numCells, cdr1Freq);
         CellFileWriter writer = new CellFileWriter(filename, sample);
         writer.writeCellsToFile();
     }
-
-    public static void makePlateExp(String cellFile, String filename, Double lambda,
-                                    Integer numWells, Integer[] concentrations, Double dropOutRate){
-        CellFileReader cellReader = new CellFileReader(cellFile);
-        Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
-        samplePlate.fillWellsExponential(cellReader.getFilename(), cellReader.getCells(), lambda);
-        PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
-        writer.writePlateFile();
-    }
-
-    private static void makePlatePoisson(String cellFile, String filename, Integer numWells,
-                                        Integer[] concentrations, Double dropOutRate){
-        CellFileReader cellReader = new CellFileReader(cellFile);
-        Double stdDev = Math.sqrt(cellReader.getCellCount());
-        Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
-        samplePlate.fillWells(cellReader.getFilename(), cellReader.getCells(), stdDev);
-        PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
-        writer.writePlateFile();
-    }
-
-    private static void makePlate(String cellFile, String filename, Double stdDev,
-                                 Integer numWells, Integer[] concentrations, Double dropOutRate){
-        CellFileReader cellReader = new CellFileReader(cellFile);
-        Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
-        samplePlate.fillWells(cellReader.getFilename(), cellReader.getCells(), stdDev);
-        PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
-        writer.writePlateFile();
-    }
-
-    private static void matchCDR3s(String graphFile, Integer lowThreshold, Integer highThreshold,
-                                   Integer occupancyDifference, Integer overlapPercent) {
-
-    }
 }

src/main/java/Equations.java

@@ -4,10 +4,6 @@ import java.math.MathContext;
 
 public abstract class Equations {
 
-    public static int getRandomNumber(int min, int max) {
-        return (int) ((Math.random() * (max - min)) + min);
-    }
-
     //pValue calculation as described in original pairSEQ paper.
     //Included for comparison with original results.
     //Not used by BiGpairSEQ for matching.

src/main/java/GraphDataObjectReader.java

@@ -1,10 +1,12 @@
 import java.io.*;
 
 public class GraphDataObjectReader {
+
     private GraphWithMapData data;
     private String filename;
+    private boolean verbose = true;
 
-    public GraphDataObjectReader(String filename) throws IOException {
+    public GraphDataObjectReader(String filename, boolean verbose) throws IOException {
         if(!filename.matches(".*\\.ser")){
             filename = filename + ".ser";
         }

src/main/java/GraphDataObjectWriter.java

@@ -1,3 +1,5 @@
+import org.jgrapht.Graph;
+
 import java.io.BufferedOutputStream;
 import java.io.FileOutputStream;
 import java.io.IOException;
@@ -7,6 +9,7 @@ public class GraphDataObjectWriter {
 
     private GraphWithMapData data;
     private String filename;
+    private boolean verbose = true;
 
     public GraphDataObjectWriter(String filename, GraphWithMapData data) {
         if(!filename.matches(".*\\.ser")){
@@ -16,13 +19,24 @@ public class GraphDataObjectWriter {
         this.data = data;
     }
 
+    public GraphDataObjectWriter(String filename, GraphWithMapData data, boolean verbose) {
+        this.verbose = verbose;
+        if(!filename.matches(".*\\.ser")){
+            filename = filename + ".ser";
+        }
+        this.filename = filename;
+        this.data = data;
+    }
+
     public void writeDataToFile() {
         try (BufferedOutputStream bufferedOut = new BufferedOutputStream(new FileOutputStream(filename));
 
              ObjectOutputStream out = new ObjectOutputStream(bufferedOut);
         ){
-            System.out.println("Writing graph and occupancy data to file. This may take some time.");
+            if(verbose) {
-            System.out.println("File I/O time is not included in results.");
+                System.out.println("Writing graph and occupancy data to file. This may take some time.");
+                System.out.println("File I/O time is not included in results.");
+            }
             out.writeObject(data);
         } catch (IOException ex) {
             ex.printStackTrace();

src/main/java/GraphMLFileReader.java

@@ -1,35 +0,0 @@
-import org.jgrapht.graph.SimpleWeightedGraph;
-import org.jgrapht.nio.graphml.GraphMLImporter;
-
-import java.io.BufferedReader;
-import java.io.IOException;
-import java.nio.file.Files;
-import java.nio.file.Path;
-
-public class GraphMLFileReader {
-
-    private String filename;
-    private SimpleWeightedGraph graph;
-
-    public GraphMLFileReader(String filename, SimpleWeightedGraph graph) {
-        if(!filename.matches(".*\\.graphml")){
-            filename = filename + ".graphml";
-        }
-        this.filename = filename;
-        this.graph = graph;
-
-        try(//don't need to close reader bc of try-with-resources auto-closing
-            BufferedReader reader = Files.newBufferedReader(Path.of(filename));
-        ){
-            GraphMLImporter<SimpleWeightedGraph, BufferedReader> importer = new GraphMLImporter<>();
-            importer.importGraph(graph, reader);
-        }
-        catch (IOException ex) {
-            System.out.println("Graph file " + filename + " not found.");
-            System.err.println(ex);
-        }
-    }
-
-    public SimpleWeightedGraph getGraph() { return graph; }
-
-}

src/main/java/GraphMLFileWriter.java

@@ -1,4 +1,8 @@
+import org.jgrapht.graph.DefaultWeightedEdge;
 import org.jgrapht.graph.SimpleWeightedGraph;
+import org.jgrapht.nio.Attribute;
+import org.jgrapht.nio.AttributeType;
+import org.jgrapht.nio.DefaultAttribute;
 import org.jgrapht.nio.dot.DOTExporter;
 import org.jgrapht.nio.graphml.GraphMLExporter;
 
@@ -7,25 +11,69 @@ import java.io.IOException;
 import java.nio.file.Files;
 import java.nio.file.Path;
 import java.nio.file.StandardOpenOption;
+import java.util.HashMap;
+import java.util.LinkedHashMap;
+import java.util.Map;
 
 public class GraphMLFileWriter {
 
     String filename;
-    SimpleWeightedGraph graph;
+    GraphWithMapData data;
 
 
-    public GraphMLFileWriter(String filename, SimpleWeightedGraph graph) {
+    public GraphMLFileWriter(String filename, GraphWithMapData data) {
         if(!filename.matches(".*\\.graphml")){
             filename = filename + ".graphml";
         }
         this.filename = filename;
-        this.graph = graph;
+        this.data = data;
     }
 
+//    public void writeGraphToFile() {
+//        try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
+//        ){
+//            GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
+//            exporter.exportGraph(graph, writer);
+//        } catch(IOException ex){
+//            System.out.println("Could not make new file named "+filename);
+//            System.err.println(ex);
+//        }
+//    }
+
     public void writeGraphToFile() {
+        SimpleWeightedGraph graph = data.getGraph();
+        Map<Integer, Integer> vertexToAlphaMap = data.getPlateVtoAMap();
+        Map<Integer, Integer> vertexToBetaMap = data.getPlateVtoBMap();
+        Map<Integer, Integer> alphaOccs = data.getAlphaWellCounts();
+        Map<Integer, Integer> betaOccs = data.getBetaWellCounts();
         try(BufferedWriter writer = Files.newBufferedWriter(Path.of(filename), StandardOpenOption.CREATE_NEW);
         ){
-            GraphMLExporter<SimpleWeightedGraph, BufferedWriter> exporter = new GraphMLExporter<>();
+            //create exporter. Let the vertex labels be the unique ids for the vertices
+            GraphMLExporter<Integer, SimpleWeightedGraph<Vertex, DefaultWeightedEdge>> exporter = new GraphMLExporter<>(v -> v.toString());
+            //set to export weights
+            exporter.setExportEdgeWeights(true);
+            //set type, sequence, and occupancy attributes for each vertex
+            exporter.setVertexAttributeProvider( v -> {
+                Map<String, Attribute> attributes = new HashMap<>();
+                if(vertexToAlphaMap.containsKey(v)) {
+                    attributes.put("type", DefaultAttribute.createAttribute("CDR3 Alpha"));
+                    attributes.put("sequence", DefaultAttribute.createAttribute(vertexToAlphaMap.get(v)));
+                    attributes.put("occupancy", DefaultAttribute.createAttribute(
+                            alphaOccs.get(vertexToAlphaMap.get(v))));
+                }
+                else if(vertexToBetaMap.containsKey(v)) {
+                    attributes.put("type", DefaultAttribute.createAttribute("CDR3 Beta"));
+                    attributes.put("sequence", DefaultAttribute.createAttribute(vertexToBetaMap.get(v)));
+                    attributes.put("occupancy", DefaultAttribute.createAttribute(
+                            betaOccs.get(vertexToBetaMap.get(v))));
+                }
+                return attributes;
+            });
+            //register the attributes
+            exporter.registerAttribute("type", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
+            exporter.registerAttribute("sequence", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
+            exporter.registerAttribute("occupancy", GraphMLExporter.AttributeCategory.NODE, AttributeType.STRING);
+            //export the graph
             exporter.exportGraph(graph, writer);
         } catch(IOException ex){
             System.out.println("Could not make new file named "+filename);
@@ -33,3 +81,4 @@ public class GraphMLFileWriter {
         }
     }
 }
+

src/main/java/GraphModificationFunctions.java

@@ -4,61 +4,75 @@ import org.jgrapht.graph.SimpleWeightedGraph;
 import java.util.ArrayList;
 import java.util.List;
 import java.util.Map;
-import java.util.Set;
 
-public abstract class GraphModificationFunctions {
+public interface GraphModificationFunctions {
 
     //remove over- and under-weight edges
-    public static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByOverlapThresholds(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
-                                                            int low, int high) {
+                                                            int low, int high, boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
-        for(DefaultWeightedEdge e: graph.edgeSet()){
+        for (DefaultWeightedEdge e : graph.edgeSet()) {
-            if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)){
+            if ((graph.getEdgeWeight(e) > high) || (graph.getEdgeWeight(e) < low)) {
-                Integer source = graph.getEdgeSource(e);
+                if(saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer weight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, weight};
+                    Integer weight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, weight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
     //Remove edges for pairs with large occupancy discrepancy
-    public static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByRelativeOccupancy(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                                   Map<Integer, Integer> alphaWellCounts,
                                                   Map<Integer, Integer> betaWellCounts,
                                                   Map<Integer, Integer> plateVtoAMap,
                                                   Map<Integer, Integer> plateVtoBMap,
-                                                  Integer maxOccupancyDifference) {
+                                                  Integer maxOccupancyDifference, boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
             Integer betaOcc = betaWellCounts.get(plateVtoBMap.get(graph.getEdgeTarget(e)));
             if (Math.abs(alphaOcc - betaOcc) >= maxOccupancyDifference) {
-                Integer source = graph.getEdgeSource(e);
+                if (saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer weight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, weight};
+                    Integer weight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, weight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
     //Remove edges for pairs where overlap size is significantly lower than the well occupancy
-    public static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static List<Integer[]> filterByOverlapPercent(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                                          Map<Integer, Integer> alphaWellCounts,
                                                          Map<Integer, Integer> betaWellCounts,
                                                          Map<Integer, Integer> plateVtoAMap,
                                                          Map<Integer, Integer> plateVtoBMap,
-                                                         Integer minOverlapPercent) {
+                                                         Integer minOverlapPercent,
+                                                         boolean saveEdges) {
         List<Integer[]> removedEdges = new ArrayList<>();
         for (DefaultWeightedEdge e : graph.edgeSet()) {
             Integer alphaOcc = alphaWellCounts.get(plateVtoAMap.get(graph.getEdgeSource(e)));
@@ -66,20 +80,27 @@ public abstract class GraphModificationFunctions {
             double weight = graph.getEdgeWeight(e);
             double min = minOverlapPercent / 100.0;
             if ((weight / alphaOcc < min) || (weight / betaOcc < min)) {
-                Integer source = graph.getEdgeSource(e);
+                if(saveEdges) {
-                Integer target = graph.getEdgeTarget(e);
+                    Integer source = graph.getEdgeSource(e);
-                Integer intWeight = (int) graph.getEdgeWeight(e);
+                    Integer target = graph.getEdgeTarget(e);
-                Integer[] edge = {source, target, intWeight};
+                    Integer intWeight = (int) graph.getEdgeWeight(e);
-                removedEdges.add(edge);
+                    Integer[] edge = {source, target, intWeight};
+                    removedEdges.add(edge);
+                }
+                else {
+                    graph.setEdgeWeight(e, 0.0);
+                }
             }
         }
-        for (Integer[] edge : removedEdges) {
+        if(saveEdges) {
-            graph.removeEdge(edge[0], edge[1]);
+            for (Integer[] edge : removedEdges) {
+                graph.removeEdge(edge[0], edge[1]);
+            }
         }
         return removedEdges;
     }
 
-    public static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
+    static void addRemovedEdges(SimpleWeightedGraph<Integer, DefaultWeightedEdge> graph,
                                        List<Integer[]> removedEdges) {
         for (Integer[] edge : removedEdges) {
             DefaultWeightedEdge e = graph.addEdge(edge[0], edge[1]);

src/main/java/GraphWithMapData.java

@@ -11,7 +11,7 @@ public class GraphWithMapData implements java.io.Serializable {
     private String sourceFilename;
     private final SimpleWeightedGraph graph;
     private Integer numWells;
-    private Integer[] wellConcentrations;
+    private Integer[] wellPopulations;
     private Integer alphaCount;
     private Integer betaCount;
     private final Map<Integer, Integer> distCellsMapAlphaKey;
@@ -31,7 +31,7 @@ public class GraphWithMapData implements java.io.Serializable {
                             Map<Integer, Integer> betaWellCounts, Duration time) {
         this.graph = graph;
         this.numWells = numWells;
-        this.wellConcentrations = wellConcentrations;
+        this.wellPopulations = wellConcentrations;
         this.alphaCount = alphaCount;
         this.betaCount = betaCount;
         this.distCellsMapAlphaKey = distCellsMapAlphaKey;
@@ -52,8 +52,8 @@ public class GraphWithMapData implements java.io.Serializable {
         return numWells;
     }
 
-    public Integer[] getWellConcentrations() {
+    public Integer[] getWellPopulations() {
-        return wellConcentrations;
+        return wellPopulations;
     }
 
     public Integer getAlphaCount() {

src/main/java/InteractiveInterface.java

@@ -1,14 +1,15 @@
 import java.io.IOException;
-import java.util.List;
+import java.util.*;
-import java.util.Scanner;
+import java.util.regex.Matcher;
-import java.util.InputMismatchException;
+import java.util.regex.Pattern;
 
 //
 public class InteractiveInterface {
 
-    final static Scanner sc = new Scanner(System.in);
+    private static final Random rand = BiGpairSEQ.getRand();
-    static int input;
+    private static final Scanner sc = new Scanner(System.in);
-    static boolean quit = false;
+    private static int input;
+    private static boolean quit = false;
 
     public static void startInteractive() {
 
@@ -26,6 +27,7 @@ public class InteractiveInterface {
                 //Need to re-do the CDR3/CDR1 matching to correspond to new pattern
                 //System.out.println("5) Generate CDR3/CDR1 occupancy graph");
                 //System.out.println("6) Simulate CDR3/CDR1 T cell matching");
+                System.out.println("8) Options");
                 System.out.println("9) About/Acknowledgments");
                 System.out.println("0) Exit");
                 try {
@@ -36,9 +38,10 @@ public class InteractiveInterface {
                         case 3 -> makeCDR3Graph();
                         case 4 -> matchCDR3s();
                         //case 6 -> matchCellsCDR1();
+                        case 8 -> mainOptions();
                         case 9 -> acknowledge();
                         case 0 -> quit = true;
-                        default -> throw new InputMismatchException("Invalid input.");
+                        default -> System.out.println("Invalid input.");
                     }
                 } catch (InputMismatchException | IOException ex) {
                     System.out.println(ex);
@@ -71,11 +74,15 @@ public class InteractiveInterface {
             System.out.println(ex);
             sc.next();
         }
-        CellSample sample = Simulator.generateCellSample(numCells, cdr1Freq);
+        CellSample sample = new CellSample(numCells, cdr1Freq);
         assert filename != null;
+        System.out.println("Writing cells to file");
         CellFileWriter writer = new CellFileWriter(filename, sample);
         writer.writeCellsToFile();
-        System.gc();
+        System.out.println("Cell sample written to: " + filename);
+        if(BiGpairSEQ.cacheCells()) {
+            BiGpairSEQ.setCellSampleInMemory(sample, filename);
+        }
     }
 
     //Output a CSV of sample plate
@@ -85,7 +92,7 @@ public class InteractiveInterface {
         Double stdDev = 0.0;
         Integer numWells = 0;
         Integer numSections;
-        Integer[] concentrations = {1};
+        Integer[] populations = {1};
         Double dropOutRate = 0.0;
         boolean poisson = false;
         boolean exponential = false;
@@ -124,10 +131,11 @@ public class InteractiveInterface {
                 }
                 case 3 -> {
                     exponential = true;
-                    System.out.println("Please enter lambda value for exponential distribution.");
+                    System.out.print("Please enter lambda value for exponential distribution: ");
                     lambda = sc.nextDouble();
                     if (lambda <= 0.0) {
-                        throw new InputMismatchException("Value must be positive.");
+                        lambda = 0.6;
+                        System.out.println("Value must be positive. Defaulting to 0.6.");
                     }
                 }
                 default -> {
@@ -140,22 +148,57 @@ public class InteractiveInterface {
             if(numWells < 1){
                 throw new InputMismatchException("No wells on plate");
             }
-            System.out.println("\nThe plate can be evenly sectioned to allow multiple concentrations of T-cells/well");
+            //choose whether to make T cell population/well random
-            System.out.println("How many sections would you like to make (minimum 1)?");
+            boolean randomWellPopulations;
-            numSections = sc.nextInt();
+            System.out.println("Randomize number of T cells in each well? (y/n)");
-            if(numSections < 1) {
+            String ans = sc.next();
-                throw new InputMismatchException("Too few sections.");
+            Pattern pattern = Pattern.compile("(?:yes|y)", Pattern.CASE_INSENSITIVE);
+            Matcher matcher = pattern.matcher(ans);
+            if(matcher.matches()){
+                randomWellPopulations = true;
             }
-            else if (numSections > numWells) {
+            else{
-                throw new InputMismatchException("Cannot have more sections than wells.");
+                randomWellPopulations = false;
             }
-            int i = 1;
+            if(randomWellPopulations) { //if T cell population/well is random
-            concentrations = new Integer[numSections];
+                numSections = numWells;
-            while(numSections > 0) {
+                Integer minPop;
-                System.out.print("Enter number of T-cells per well in section " + i +": ");
+                Integer maxPop;
-                concentrations[i - 1] = sc.nextInt();
+                System.out.print("Please enter minimum number of T cells in a well: ");
-                i++;
+                minPop = sc.nextInt();
-                numSections--;
+                if(minPop < 1) {
+                    throw new InputMismatchException("Minimum well population must be positive");
+                }
+                System.out.println("Please enter maximum number of T cells in a well: ");
+                maxPop = sc.nextInt();
+                if(maxPop < minPop) {
+                    throw new InputMismatchException("Max well population must be greater than min well population");
+                }
+                //maximum should be inclusive, so need to add one to max of randomly generated values
+                populations = rand.ints(minPop, maxPop + 1)
+                        .limit(numSections)
+                        .boxed()
+                        .toArray(Integer[]::new);
+                System.out.print("Populations: ");
+                System.out.println(Arrays.toString(populations));
+            }
+            else{ //if T cell population/well is not random
+                System.out.println("\nThe plate can be evenly sectioned to allow different numbers of T cells per well.");
+                System.out.println("How many sections would you like to make (minimum 1)?");
+                numSections = sc.nextInt();
+                if (numSections < 1) {
+                    throw new InputMismatchException("Too few sections.");
+                } else if (numSections > numWells) {
+                    throw new InputMismatchException("Cannot have more sections than wells.");
+                }
+                int i = 1;
+                populations = new Integer[numSections];
+                while (numSections > 0) {
+                    System.out.print("Enter number of T cells per well in section " + i + ": ");
+                    populations[i - 1] = sc.nextInt();
+                    i++;
+                    numSections--;
+                }
             }
             System.out.println("\nErrors in amplification can induce a well dropout rate for sequences");
             System.out.print("Enter well dropout rate (0.0 to 1.0): ");
@@ -167,27 +210,38 @@ public class InteractiveInterface {
             System.out.println(ex);
             sc.next();
         }
-        System.out.println("Reading Cell Sample file: " + cellFile);
         assert cellFile != null;
-        CellFileReader cellReader = new CellFileReader(cellFile);
+        CellSample cells;
+        if (cellFile.equals(BiGpairSEQ.getCellFilename())){
+            cells = BiGpairSEQ.getCellSampleInMemory();
+        }
+        else {
+            System.out.println("Reading Cell Sample file: " + cellFile);
+            CellFileReader cellReader = new CellFileReader(cellFile);
+            cells = cellReader.getCellSample();
+            if(BiGpairSEQ.cacheCells()) {
+                BiGpairSEQ.setCellSampleInMemory(cells, cellFile);
+            }
+        }
+        assert filename != null;
+        Plate samplePlate;
+        PlateFileWriter writer;
         if(exponential){
-            Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
+            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, lambda, true);
-            samplePlate.fillWellsExponential(cellReader.getFilename(), cellReader.getCells(), lambda);
+            writer = new PlateFileWriter(filename, samplePlate);
-            PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
-            writer.writePlateFile();
         }
         else {
             if (poisson) {
-                stdDev = Math.sqrt(cellReader.getCellCount()); //gaussian with square root of elements approximates poisson
+                stdDev = Math.sqrt(cells.getCellCount()); //gaussian with square root of elements approximates poisson
             }
-            Plate samplePlate = new Plate(numWells, dropOutRate, concentrations);
+            samplePlate = new Plate(cells, cellFile, numWells, populations, dropOutRate, stdDev, false);
-            samplePlate.fillWells(cellReader.getFilename(), cellReader.getCells(), stdDev);
+            writer = new PlateFileWriter(filename, samplePlate);
-            assert filename != null;
+        }
-            PlateFileWriter writer = new PlateFileWriter(filename, samplePlate);
+        System.out.println("Writing Sample Plate to file");
-            System.out.println("Writing Sample Plate to file");
+        writer.writePlateFile();
-            writer.writePlateFile();
+        System.out.println("Sample Plate written to file: " + filename);
-            System.out.println("Sample Plate written to file: " + filename);
+        if(BiGpairSEQ.cachePlate()) {
-            System.gc();
+            BiGpairSEQ.setPlateInMemory(samplePlate, filename);
         }
     }
 
@@ -196,7 +250,6 @@ public class InteractiveInterface {
         String filename = null;
         String cellFile = null;
         String plateFile = null;
-
         try {
             String str = "\nGenerating bipartite weighted graph encoding occupancy overlap data ";
             str = str.concat("\nrequires a cell sample file and a sample plate file.");
@@ -212,14 +265,37 @@ public class InteractiveInterface {
             System.out.println(ex);
             sc.next();
         }
-        System.out.println("Reading Cell Sample file: " + cellFile);
+
         assert cellFile != null;
-        CellFileReader cellReader = new CellFileReader(cellFile);
+        CellSample cellSample;
-        System.out.println("Reading Sample Plate file: " + plateFile);
+        //check if cells are already in memory
+        if(cellFile.equals(BiGpairSEQ.getCellFilename()) && BiGpairSEQ.getCellSampleInMemory() != null) {
+            cellSample = BiGpairSEQ.getCellSampleInMemory();
+        }
+        else {
+            System.out.println("Reading Cell Sample file: " + cellFile);
+            CellFileReader cellReader = new CellFileReader(cellFile);
+            cellSample = cellReader.getCellSample();
+            if(BiGpairSEQ.cacheCells()) {
+                BiGpairSEQ.setCellSampleInMemory(cellSample, cellFile);
+            }
+        }
+
         assert plateFile != null;
-        PlateFileReader plateReader = new PlateFileReader(plateFile);
+        Plate plate;
-        Plate plate = new Plate(plateReader.getFilename(), plateReader.getWells());
+        //check if plate is already in memory
-        if (cellReader.getCells().size() == 0){
+        if(plateFile.equals(BiGpairSEQ.getPlateFilename())){
+            plate = BiGpairSEQ.getPlateInMemory();
+        }
+        else {
+            System.out.println("Reading Sample Plate file: " + plateFile);
+            PlateFileReader plateReader = new PlateFileReader(plateFile);
+            plate = plateReader.getSamplePlate();
+            if(BiGpairSEQ.cachePlate()) {
+                BiGpairSEQ.setPlateInMemory(plate, plateFile);
+            }
+        }
+        if (cellSample.getCells().size() == 0){
             System.out.println("No cell sample found.");
             System.out.println("Returning to main menu.");
         }
@@ -228,13 +304,22 @@ public class InteractiveInterface {
             System.out.println("Returning to main menu.");
         }
         else{
-            List<Integer[]> cells = cellReader.getCells();
+            GraphWithMapData data = Simulator.makeGraph(cellSample, plate, true);
-            GraphWithMapData data = Simulator.makeGraph(cells, plate, true);
             assert filename != null;
-            GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
+            if(BiGpairSEQ.outputBinary()) {
-            dataWriter.writeDataToFile();
+                GraphDataObjectWriter dataWriter = new GraphDataObjectWriter(filename, data);
-            System.out.println("Graph and Data file written to: " + filename);
+                dataWriter.writeDataToFile();
-            System.gc();
+                System.out.println("Serialized binary graph/data file written to: " + filename);
+            }
+            if(BiGpairSEQ.outputGraphML()) {
+                GraphMLFileWriter graphMLWriter = new GraphMLFileWriter(filename, data);
+                graphMLWriter.writeGraphToFile();
+                System.out.println("GraphML file written to: " + filename);
+            }
+            if(BiGpairSEQ.cacheGraph()) {
+                BiGpairSEQ.setGraphInMemory(data, filename);
+
+            }
         }
     }
 
@@ -256,17 +341,28 @@ public class InteractiveInterface {
             System.out.println("\nWhat is the minimum number of CDR3 alpha/beta overlap wells to attempt matching?");
             lowThreshold = sc.nextInt();
             if(lowThreshold < 1){
-                throw new InputMismatchException("Minimum value for low threshold set to 1");
+                lowThreshold = 1;
+                System.out.println("Value for low occupancy overlap threshold must be positive");
+                System.out.println("Value for low occupancy overlap threshold set to 1");
             }
             System.out.println("\nWhat is the maximum number of CDR3 alpha/beta overlap wells to attempt matching?");
             highThreshold = sc.nextInt();
-            System.out.println("\nWhat is the maximum difference in alpha/beta occupancy to attempt matching?");
+            if(highThreshold < lowThreshold) {
-            maxOccupancyDiff = sc.nextInt();
+                highThreshold = lowThreshold;
-            System.out.println("\nWell overlap percentage = pair overlap / sequence occupancy");
+                System.out.println("Value for high occupancy overlap threshold must be >= low overlap threshold");
-            System.out.println("What is the minimum well overlap percentage to attempt matching? (0 to 100)");
+                System.out.println("Value for high occupancy overlap threshold set to " + lowThreshold);
+            }
+            System.out.println("What is the minimum percentage of a sequence's wells in alpha/beta overlap to attempt matching? (0 - 100)");
             minOverlapPercent = sc.nextInt();
             if (minOverlapPercent < 0 || minOverlapPercent > 100) {
-                throw new InputMismatchException("Value outside range. Minimum percent set to 0");
+                System.out.println("Value outside range. Minimum occupancy overlap percentage set to 0");
+            }
+            System.out.println("\nWhat is the maximum difference in alpha/beta occupancy to attempt matching?");
+            maxOccupancyDiff = sc.nextInt();
+            if (maxOccupancyDiff < 0) {
+                maxOccupancyDiff = 0;
+                System.out.println("Maximum allowable difference in alpha/beta occupancy must be nonnegative");
+                System.out.println("Maximum allowable difference in alpha/beta occupancy set to 0");
             }
         } catch (InputMismatchException ex) {
             System.out.println(ex);
@@ -275,17 +371,15 @@ public class InteractiveInterface {
         assert graphFilename != null;
         //check if this is the same graph we already have in memory.
         GraphWithMapData data;
-        if(!(graphFilename.equals(BiGpairSEQ.getGraphFilename())) || BiGpairSEQ.getGraph() == null) {
+        if(graphFilename.equals(BiGpairSEQ.getGraphFilename())) {
-            BiGpairSEQ.clearGraph();
+            data = BiGpairSEQ.getGraphInMemory();
-            //read object data from file
-            GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename);
-            data = dataReader.getData();
-            //set new graph in memory and new filename
-            BiGpairSEQ.setGraph(data);
-            BiGpairSEQ.setGraphFilename(graphFilename);
         }
         else {
-            data = BiGpairSEQ.getGraph();
+            GraphDataObjectReader dataReader = new GraphDataObjectReader(graphFilename, true);
+            data = dataReader.getData();
+            if(BiGpairSEQ.cacheGraph()) {
+                BiGpairSEQ.setGraphInMemory(data, graphFilename);
+            }
         }
         //simulate matching
         MatchingResult results = Simulator.matchCDR3s(data, graphFilename, lowThreshold, highThreshold, maxOccupancyDiff,
@@ -296,7 +390,6 @@ public class InteractiveInterface {
         System.out.println("Writing results to file");
         writer.writeResultsToFile();
         System.out.println("Results written to file: " + filename);
-        System.gc();
     }
 
     ///////
@@ -403,6 +496,79 @@ public class InteractiveInterface {
 //        }
 //    }
 
+    private static void mainOptions(){
+        boolean backToMain = false;
+        while(!backToMain) {
+            System.out.println("\n--------------OPTIONS---------------");
+            System.out.println("1) Turn " + getOnOff(!BiGpairSEQ.cacheCells()) + " cell sample file caching");
+            System.out.println("2) Turn " + getOnOff(!BiGpairSEQ.cachePlate()) + " plate file caching");
+            System.out.println("3) Turn " + getOnOff(!BiGpairSEQ.cacheGraph()) + " graph/data file caching");
+            System.out.println("4) Turn " + getOnOff(!BiGpairSEQ.outputBinary()) + " serialized binary graph output");
+            System.out.println("5) Turn " + getOnOff(!BiGpairSEQ.outputGraphML()) + " GraphML graph output");
+            System.out.println("6) Maximum weight matching algorithm options");
+            System.out.println("0) Return to main menu");
+            try {
+                input = sc.nextInt();
+                switch (input) {
+                    case 1 -> BiGpairSEQ.setCacheCells(!BiGpairSEQ.cacheCells());
+                    case 2 -> BiGpairSEQ.setCachePlate(!BiGpairSEQ.cachePlate());
+                    case 3 -> BiGpairSEQ.setCacheGraph(!BiGpairSEQ.cacheGraph());
+                    case 4 -> BiGpairSEQ.setOutputBinary(!BiGpairSEQ.outputBinary());
+                    case 5 -> BiGpairSEQ.setOutputGraphML(!BiGpairSEQ.outputGraphML());
+                    case 6 -> algorithmOptions();
+                    case 0 -> backToMain = true;
+                    default -> System.out.println("Invalid input");
+                }
+            } catch (InputMismatchException ex) {
+                System.out.println(ex);
+                sc.next();
+            }
+        }
+    }
+
+    /**
+     * Helper function for printing menu items in mainOptions(). Returns a string based on the value of parameter.
+     *
+     * @param b - a boolean value
+     * @return String "on" if b is true, "off" if b is false
+     */
+    private static String getOnOff(boolean b) {
+        if (b) { return "on";}
+        else { return "off"; }
+    }
+
+    private static void algorithmOptions(){
+        boolean backToOptions = false;
+        while(!backToOptions) {
+            System.out.println("\n---------ALGORITHM OPTIONS----------");
+            System.out.println("1) Use scaling algorithm by Duan and Su.");
+            System.out.println("2) Use LEDA book algorithm with Fibonacci heap priority queue");
+            System.out.println("3) Use LEDA book algorithm with pairing heap priority queue");
+            System.out.println("0) Return to Options menu");
+            try {
+                input = sc.nextInt();
+                switch (input) {
+                    case 1 -> System.out.println("This option is not yet implemented. Choose another.");
+                    case 2 -> {
+                        BiGpairSEQ.setFibonacciHeap();
+                        System.out.println("MWM algorithm set to LEDA with Fibonacci heap");
+                        backToOptions = true;
+                    }
+                    case 3 -> {
+                        BiGpairSEQ.setPairingHeap();
+                        System.out.println("MWM algorithm set to LEDA with pairing heap");
+                        backToOptions = true;
+                    }
+                    case 0 -> backToOptions = true;
+                    default -> System.out.println("Invalid input");
+                }
+            } catch (InputMismatchException ex) {
+                System.out.println(ex);
+                sc.next();
+            }
+        }
+    }
+
     private static void acknowledge(){
         System.out.println("This program simulates BiGpairSEQ, a graph theory based adaptation");
         System.out.println("of the pairSEQ algorithm for pairing T cell receptor sequences.");

src/main/java/MatchingResult.java

@@ -21,15 +21,15 @@ public class MatchingResult {
          * well populations *
          * total alphas found *
          * total betas found *
-         * high overlap threshold
+         * high overlap threshold *
-         * low overlap threshold
+         * low overlap threshold *
-         * maximum occupancy difference
+         * maximum occupancy difference *
-         * minimum overlap percent
+         * minimum overlap percent *
-         * pairing attempt rate
+         * pairing attempt rate *
-         * correct pairing count
+         * correct pairing count *
-         * incorrect pairing count
+         * incorrect pairing count *
-         * pairing error rate
+         * pairing error rate *
-         * simulation time
+         * simulation time (seconds)
          */
         this.metadata = metadata;
         this.comments = new ArrayList<>();
@@ -91,6 +91,22 @@ public class MatchingResult {
         return Integer.parseInt(metadata.get("total beta count"));
     }
 
-    //put in the rest of these methods following the same pattern
+    public Integer getHighOverlapThreshold() { return Integer.parseInt(metadata.get("high overlap threshold"));}
+
+    public Integer getLowOverlapThreshold() { return Integer.parseInt(metadata.get("low overlap threshold"));}
+
+    public Integer getMaxOccupancyDifference() { return Integer.parseInt(metadata.get("maximum occupancy difference"));}
+
+    public Integer getMinOverlapPercent() { return Integer.parseInt(metadata.get("minimum overlap percent"));}
+
+    public Double getPairingAttemptRate() { return Double.parseDouble(metadata.get("pairing attempt rate"));}
+
+    public Integer getCorrectPairingCount() { return Integer.parseInt(metadata.get("correct pairing count"));}
+
+    public Integer getIncorrectPairingCount() { return Integer.parseInt(metadata.get("incorrect pairing count"));}
+
+    public Double getPairingErrorRate() { return Double.parseDouble(metadata.get("pairing error rate"));}
+
+    public String getSimulationTime() { return metadata.get("simulation time (seconds)"); }
 
 }

src/main/java/Plate.java

@@ -8,9 +8,11 @@ TODO: Implement discrete frequency distributions using Vose's Alias Method
 import java.util.*;
 
 public class Plate {
+    private CellSample cells;
     private String sourceFile;
+    private String filename;
     private List<List<Integer[]>> wells;
-    private Random rand = new Random();
+    private final Random rand = BiGpairSEQ.getRand();
     private int size;
     private double error;
     private Integer[] populations;
@@ -18,6 +20,25 @@ public class Plate {
     private double lambda;
     boolean exponential = false;
 
+    public Plate(CellSample cells, String cellFilename, int numWells, Integer[] populations,
+                 double dropoutRate, double stdDev_or_lambda, boolean exponential){
+        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);
+        }
+    }
+
 
     public Plate(int size, double error, Integer[] populations) {
         this.size = size;
@@ -26,8 +47,9 @@ public class Plate {
         wells = new ArrayList<>();
     }
 
-    public Plate(String sourceFileName, List<List<Integer[]>> wells) {
+    //constructor for returning a Plate from a PlateFileReader
-        this.sourceFile = sourceFileName;
+    public Plate(String filename, List<List<Integer[]>> wells) {
+        this.filename = filename;
         this.wells = wells;
         this.size = wells.size();
 
@@ -43,15 +65,13 @@ public class Plate {
         }
     }
 
-    public void fillWellsExponential(String sourceFileName, List<Integer[]> cells, double lambda){
+    private void fillWellsExponential(List<Integer[]> cells, double lambda){
         this.lambda = lambda;
         exponential = true;
-        sourceFile = sourceFileName;
         int numSections = populations.length;
         int section = 0;
         double m;
         int n;
-        int test=0;
         while (section < numSections){
             for (int i = 0; i < (size / numSections); i++) {
                 List<Integer[]> well = new ArrayList<>();
@@ -61,13 +81,6 @@ public class Plate {
                         m = (Math.log10((1 - rand.nextDouble()))/(-lambda)) * Math.sqrt(cells.size());
                     } while (m >= cells.size() || m < 0);
                     n = (int) Math.floor(m);
-                    //n = Equations.getRandomNumber(0, cells.size());
-                    // was testing generating the cell sample file with exponential dist, then sampling flat here
-                    //that would be more realistic
-                    //But would mess up other things in the simulation with how I've coded it.
-                    if(n > test){
-                        test = n;
-                    }
                     Integer[] cellToAdd = cells.get(n).clone();
                     for(int k = 0; k < cellToAdd.length; k++){
                         if(Math.abs(rand.nextDouble()) < error){//error applied to each seqeunce
@@ -80,12 +93,10 @@ public class Plate {
             }
             section++;
         }
-        System.out.println("Highest index: " +test);
     }
 
-    public void fillWells(String sourceFileName, List<Integer[]> cells, double stdDev) {
+    private void fillWells( List<Integer[]> cells, double stdDev) {
         this.stdDev = stdDev;
-        sourceFile = sourceFileName;
         int numSections = populations.length;
         int section = 0;
         double m;
@@ -168,4 +179,6 @@ public class Plate {
     public String getSourceFileName() {
         return sourceFile;
     }
+
+    public String getFilename() { return filename; }
 }

src/main/java/PlateFileReader.java

@@ -56,11 +56,8 @@ public class PlateFileReader {
 
     }
 
-    public List<List<Integer[]>> getWells() {
+    public Plate getSamplePlate() {
-        return wells;
+        return new Plate(filename, wells);
     }
 
-    public String getFilename() {
-        return filename;
-    }
 }

src/main/java/PlateFileWriter.java

@@ -16,7 +16,7 @@ public class PlateFileWriter {
     private Double error;
     private String filename;
     private String sourceFileName;
-    private Integer[] concentrations;
+    private Integer[] populations;
     private boolean isExponential = false;
 
     public PlateFileWriter(String filename, Plate plate) {
@@ -35,8 +35,8 @@ public class PlateFileWriter {
         }
         this.error = plate.getError();
         this.wells = plate.getWells();
-        this.concentrations = plate.getPopulations();
+        this.populations = plate.getPopulations();
-        Arrays.sort(concentrations);
+        Arrays.sort(populations);
     }
 
     public void writePlateFile(){
@@ -73,14 +73,12 @@ public class PlateFileWriter {
 //            rows.add(tmp);
 //        }
 
-        //get list of well populations
+        //make string out of populations array
-        List<Integer> wellPopulations = Arrays.asList(concentrations);
-        //make string out of populations list
         StringBuilder populationsStringBuilder = new StringBuilder();
-        populationsStringBuilder.append(wellPopulations.remove(0).toString());
+        populationsStringBuilder.append(populations[0].toString());
-        for(Integer i: wellPopulations){
+        for(int i = 1; i < populations.length; i++){
             populationsStringBuilder.append(", ");
-            populationsStringBuilder.append(i.toString());
+            populationsStringBuilder.append(populations[i].toString());
         }
         String wellPopulationsString = populationsStringBuilder.toString();
 

src/main/java/Simulator.java

@@ -1,9 +1,9 @@
-import org.jgrapht.Graph;
 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;
@@ -14,40 +14,19 @@ import java.time.Duration;
 import java.util.*;
 import java.util.stream.IntStream;
 
+import static java.lang.Float.*;
+
 //NOTE: "sequence" in method and variable names refers to a peptide sequence from a simulated T cell
-public class Simulator {
+public class Simulator implements GraphModificationFunctions {
     private static final int cdr3AlphaIndex = 0;
     private static final int cdr3BetaIndex = 1;
     private static final int cdr1AlphaIndex = 2;
     private static final int cdr1BetaIndex = 3;
 
-    public static CellSample generateCellSample(Integer numDistinctCells, Integer cdr1Freq) {
-        //In real T cells, CDR1s have about one third the diversity of CDR3s
-        List<Integer> numbersCDR3 = new ArrayList<>();
-        List<Integer> numbersCDR1 = new ArrayList<>();
-        Integer numDistCDR3s = 2 * numDistinctCells + 1;
-        IntStream.range(1, numDistCDR3s + 1).forEach(i -> numbersCDR3.add(i));
-        IntStream.range(numDistCDR3s + 1, numDistCDR3s + 1 + (numDistCDR3s / cdr1Freq) + 1).forEach(i -> numbersCDR1.add(i));
-        Collections.shuffle(numbersCDR3);
-        Collections.shuffle(numbersCDR1);
-
-        //Each cell represented by 4 values
-        //two CDR3s, and two CDR1s. First two values are CDR3s (alpha, beta), second two are CDR1s (alpha, beta)
-        List<Integer[]> distinctCells = new ArrayList<>();
-        for(int i = 0; i < numbersCDR3.size() - 1; i = i + 2){
-            Integer tmpCDR3a = numbersCDR3.get(i);
-            Integer tmpCDR3b = numbersCDR3.get(i+1);
-            Integer tmpCDR1a = numbersCDR1.get(i % numbersCDR1.size());
-            Integer tmpCDR1b = numbersCDR1.get((i+1) % numbersCDR1.size());
-            Integer[] tmp = {tmpCDR3a, tmpCDR3b, tmpCDR1a, tmpCDR1b};
-            distinctCells.add(tmp);
-        }
-        return new CellSample(distinctCells, cdr1Freq);
-    }
-
     //Make the graph needed for matching CDR3s
-    public static GraphWithMapData makeGraph(List<Integer[]> distinctCells, Plate samplePlate, boolean verbose) {
+    public static GraphWithMapData makeGraph(CellSample cellSample, Plate samplePlate, boolean verbose) {
         Instant start = Instant.now();
+        List<Integer[]> distinctCells = cellSample.getCells();
         int[] alphaIndex = {cdr3AlphaIndex};
         int[] betaIndex = {cdr3BetaIndex};
 
@@ -135,7 +114,7 @@ public class Simulator {
                 distCellsMapAlphaKey, plateVtoAMap, plateVtoBMap, plateAtoVMap,
                 plateBtoVMap, alphaWellCounts, betaWellCounts, time);
         //Set source file name in graph to name of sample plate
-        output.setSourceFilename(samplePlate.getSourceFileName());
+        output.setSourceFilename(samplePlate.getFilename());
         //return GraphWithMapData object
         return output;
     }
@@ -145,8 +124,8 @@ public class Simulator {
                                             Integer highThreshold, Integer maxOccupancyDifference,
                                             Integer minOverlapPercent, boolean verbose) {
         Instant start = Instant.now();
-        //Integer arrays will contain TO VERTEX, FROM VERTEX, and WEIGHT (which I'll need to cast to double)
         List<Integer[]> removedEdges = new ArrayList<>();
+        boolean saveEdges = BiGpairSEQ.cacheGraph();
         int numWells = data.getNumWells();
         Integer alphaCount = data.getAlphaCount();
         Integer betaCount = data.getBetaCount();
@@ -159,33 +138,50 @@ public class Simulator {
 
         //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");}
-        removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold));
+        removedEdges.addAll(GraphModificationFunctions.filterByOverlapThresholds(graph, lowThreshold, highThreshold, saveEdges));
         if(verbose){System.out.println("Over- and under-weight edges removed");}
 
         //remove edges between vertices with too small an overlap size, add those to removed edge list
         if(verbose){System.out.println("Eliminating edges with weights less than " + minOverlapPercent.toString() +
                 " percent of vertex occupancy value.");}
         removedEdges.addAll(GraphModificationFunctions.filterByOverlapPercent(graph, alphaWellCounts, betaWellCounts,
-                plateVtoAMap, plateVtoBMap, minOverlapPercent));
+                plateVtoAMap, plateVtoBMap, minOverlapPercent, saveEdges));
         if(verbose){System.out.println("Edges with weights too far below a vertex occupancy value removed");}
 
         //Filter by relative occupancy
         if(verbose){System.out.println("Eliminating edges between vertices with occupancy difference > "
                 + maxOccupancyDifference);}
         removedEdges.addAll(GraphModificationFunctions.filterByRelativeOccupancy(graph, alphaWellCounts, betaWellCounts,
-                plateVtoAMap, plateVtoBMap, maxOccupancyDifference));
+                plateVtoAMap, plateVtoBMap, maxOccupancyDifference, saveEdges));
         if(verbose){System.out.println("Edges between vertices of with excessively different occupancy values " +
                 "removed");}
 
         //Find Maximum Weighted Matching
         //using jheaps library class PairingHeap for improved efficiency
         if(verbose){System.out.println("Finding maximum weighted matching");}
-        //Attempting to use addressable heap to improve performance
+        MaximumWeightBipartiteMatching maxWeightMatching;
-        MaximumWeightBipartiteMatching maxWeightMatching =
+        //Use correct heap type for priority queue
-                new MaximumWeightBipartiteMatching(graph,
+        String heapType = BiGpairSEQ.getPriorityQueueHeapType();
+        switch (heapType) {
+            case "PAIRING" -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
                         plateVtoAMap.keySet(),
                         plateVtoBMap.keySet(),
                         i -> new PairingHeap(Comparator.naturalOrder()));
+            }
+            case "FIBONACCI" -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
+                        plateVtoAMap.keySet(),
+                        plateVtoBMap.keySet(),
+                        i -> new FibonacciHeap(Comparator.naturalOrder()));
+            }
+            default -> {
+                maxWeightMatching = new MaximumWeightBipartiteMatching(graph,
+                        plateVtoAMap.keySet(),
+                        plateVtoBMap.keySet());
+            }
+        }
+        //get the matching
         MatchingAlgorithm.Matching<String, DefaultWeightedEdge> graphMatching = maxWeightMatching.getMatching();
         if(verbose){System.out.println("Matching completed");}
         Instant stop = Instant.now();
@@ -241,16 +237,23 @@ public class Simulator {
         }
 
         //Metadata comments for CSV file
+        String algoType = "LEDA book with heap: " + heapType;
         int min = Math.min(alphaCount, betaCount);
         //rate of attempted matching
         double attemptRate = (double) (trueCount + falseCount) / min;
         BigDecimal attemptRateTrunc = new BigDecimal(attemptRate, mc);
         //rate of pairing error
         double pairingErrorRate = (double) falseCount / (trueCount + falseCount);
-        BigDecimal pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
+        BigDecimal pairingErrorRateTrunc;
-        //get list of well concentrations
+        if(pairingErrorRate == NaN || pairingErrorRate == POSITIVE_INFINITY || pairingErrorRate == NEGATIVE_INFINITY) {
-        Integer[] wellPopulations = data.getWellConcentrations();
+            pairingErrorRateTrunc = new BigDecimal(-1, mc);
-        //make string out of concentrations list
+        }
+        else{
+            pairingErrorRateTrunc = new BigDecimal(pairingErrorRate, mc);
+        }
+        //get list of well populations
+        Integer[] wellPopulations = data.getWellPopulations();
+        //make string out of populations list
         StringBuilder populationsStringBuilder = new StringBuilder();
         populationsStringBuilder.append(wellPopulations[0].toString());
         for(int i = 1; i < wellPopulations.length; i++){
@@ -265,18 +268,19 @@ public class Simulator {
         Map<String, String> metadata = new LinkedHashMap<>();
         metadata.put("sample plate filename", data.getSourceFilename());
         metadata.put("graph filename", dataFilename);
+        metadata.put("algorithm type", algoType);
         metadata.put("well populations", wellPopulationsString);
         metadata.put("total alphas found", alphaCount.toString());
         metadata.put("total betas found", betaCount.toString());
         metadata.put("high overlap threshold", highThreshold.toString());
         metadata.put("low overlap threshold", lowThreshold.toString());
-        metadata.put("maximum occupancy difference", maxOccupancyDifference.toString());
         metadata.put("minimum overlap percent", minOverlapPercent.toString());
+        metadata.put("maximum occupancy difference", maxOccupancyDifference.toString());
         metadata.put("pairing attempt rate", attemptRateTrunc.toString());
         metadata.put("correct pairing count", Integer.toString(trueCount));
         metadata.put("incorrect pairing count", Integer.toString(falseCount));
         metadata.put("pairing error rate", pairingErrorRateTrunc.toString());
-        metadata.put("simulation time", nf.format(time.toSeconds()));
+        metadata.put("simulation time (seconds)", nf.format(time.toSeconds()));
         //create MatchingResult object
         MatchingResult output = new MatchingResult(metadata, header, allResults, matchMap, time);
         if(verbose){
@@ -285,10 +289,11 @@ public class Simulator {
             }
         }
 
-        //put the removed edges back on the graph
+        if(saveEdges) {
-        System.out.println("Restoring removed edges to graph.");
+            //put the removed edges back on the graph
-        GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
+            System.out.println("Restoring removed edges to graph.");
-
+            GraphModificationFunctions.addRemovedEdges(graph, removedEdges);
+        }
         //return MatchingResult object
         return output;
     }
@@ -664,7 +669,7 @@ public class Simulator {
 
     private static Map<Integer, Integer> makeVertexToSequenceMap(Map<Integer, Integer> sequences, Integer startValue) {
         Map<Integer, Integer> map = new LinkedHashMap<>(); //LinkedHashMap to preserve order of entry
-        Integer index = startValue;
+        Integer index = startValue; //is this necessary? I don't think I use this.
         for (Integer k: sequences.keySet()) {
             map.put(index, k);
             index++;

src/main/java/Vertex.java

@@ -1,14 +1,20 @@
+
+
 public class Vertex {
-    private final Integer peptide;
+    private final Integer vertexLabel;
+    private final Integer sequence;
     private final Integer occupancy;
 
-    public Vertex(Integer peptide, Integer occupancy) {
+    public Vertex(Integer vertexLabel, Integer sequence, Integer occupancy) {
-        this.peptide = peptide;
+        this.vertexLabel = vertexLabel;
+        this.sequence = sequence;
         this.occupancy = occupancy;
     }
 
-    public Integer getPeptide() {
+    public Integer getVertexLabel() { return vertexLabel; }
-        return peptide;
+
+    public Integer getSequence() {
+        return sequence;
     }
 
     public Integer getOccupancy() {