class SimulationSetup {
  String seqTot; // String of all bases
  int lenTot; // Total number of bases
  int[] starts; // Start index of each strand (starting on 0)
  String[] sequences; // String of each strand's bases
  int[] seqIds;  // Strand ID of each base
  ArrayList<int[]> foldSteps = new ArrayList<int[]>();
  int[] strandOrder; // Optimal order of strands to minimize knotting
  int[] residueOrder; // New residue order based off of strandOrder
  Short[] pairTable; // Holds index of base paired with the baseId at each index (-1 if unpaired)
  ArrayList<Short> pairTableOptimal = new ArrayList<Short>(); // Stores the position in the new order that each newly ordered residue is paired with
  //ArrayList<String> nonCanonicals = new ArrayList<String>(); // example: "id1 id2 cHS"
  HashMap<String, String> nonCanonicals = new HashMap<String, String>(); // "id1 id2", "cww" ID1 ALWAYS LESS THAN ID2!
  String errorMessage = "";
  String warningMessage = "";
  
  SimulationSetup cloneFunction() {
    SimulationSetup result = new SimulationSetup();
    result.nonCanonicals.putAll(this.nonCanonicals);
    result.pairTableOptimal.addAll(this.pairTableOptimal);
    result.seqTot = this.seqTot;
    result.lenTot = this.lenTot;
    result.starts = intArrayClone(this.starts);
    result.sequences = stringArrayClone(this.sequences);
    result.seqIds = intArrayClone(this.seqIds);
    result.foldSteps.addAll(this.foldSteps);
    result.strandOrder = intArrayClone(this.strandOrder);
    result.residueOrder = intArrayClone(this.residueOrder);
    result.pairTable = shortArrayClone(this.pairTable);
    result.errorMessage = this.errorMessage;
    result.warningMessage = this.warningMessage;
    
    return result;
  }
  
  /** Determines file type, calls initFileType */
  SimulationState initFile(String fileName, String fileType) {
    if (!BROWSER) {
      if (debugPrints) { println("\nStarting SimulationSetup.initFile for input file " + fileName); }
      
      File loadFile = new File(fileName);
      if (!loadFile.isFile()) {
        println("Could not find file with name " + fileName); 
        errorMessage = "File  " + fileName + "  could not be loaded";
        return new SimulationState();
      }
      
      String userFileName = userFile.getName();
      surface.setTitle(userFileName); // Set the title of the window
      
      if (fileType == null) { // default
        fileType = "";
        for (int i = userFileName.length(); i > 0; --i) {
          if ( ".".equals(userFileName.substring(i-1, i)) ) {
            if (i < userFileName.length()) {
              fileType = userFileName.substring(i);
            }
            break;
          }
        }
      }
      if (debugPrints) { println("File type is " + fileType); }
      return initFileType(fileName, fileType);
    }
    else {
      //println("\nStarting SimulationSetup.initFile in browser!");
      if (javascript != null) {
        if (fileType == null) { // default
          fileType = javascript.getFileTypeBrowser();
        }
        if (debugPrints) { println("File type is " + fileType); }
        return initFileType(null, fileType);
      }
      else {
        if (debugPrints) { println("JAVASCRIPT NOT BOUND!"); }
        errorMessage = "Error in communication with web browser.";
        return new SimulationState();
      }
    }
  }
  
  
  /** Determines which file parser to call, and which setup method */
  SimulationState initFileType(String fileName, String fileType) {
    if ("ct".equals(fileType) || "bpseq".equals(fileType)) {
      initColumnFile(fileName, fileType);
    }
    else if ("dbn".equals(fileType)) {
      initBracketFile(fileName);
    }
    else if ("rs".equals(fileType)) {
      SimulationState simState = initSaveFile(fileName);
      if ("".equals(errorMessage)) {
        checkLengthForRigidLoops = false;
        return simState;
      }
      else {
        if (debugPrints) {
          println("InitFileException caught from rs file type in initFile:");
          println(errorMessage);
        }
        return new SimulationState();
      }
    }
    else if ("nts".equals(fileType)) {
      initNtsPairs(fileName);
    }
    else if ("seq".equals(fileType)) {
      initSequenceFile(fileName);
    }
    else {
      // Prompt user to select file type
      String[] fileTypes = {"ct", "bpseq", "dbn", "rs", "nts", "seq"};
      int itemHeight = 40;
      float listHeight = fileTypes.length * itemHeight;
      float yStart = 100;
      if (listHeight + yStart > (height - 135)) {
        listHeight = int((height - 135 - yStart)/itemHeight) * itemHeight;
        if (listHeight < itemHeight) {
          listHeight = (float)itemHeight;
        }
      }
      fileTypeSelector = new Listbox(30, yStart, 150, listHeight, itemHeight, fileTypes, ++displaySelectorId);
      selectingFileType = true;
      if (debugPrints) { println("Please input file type"); }
      errorMessage += "Could not identify file type automatically";
      return new SimulationState();
    }
    
    if ("".equals(errorMessage)) {
      initTables();
      
      if (debugPrints) { setupDebug(); }
      ds = new DisplaySettings();
      if (lenTot > 1300) {
        ds.labelMode = false;
        if (lenTot > 2500) {
          ds.outlineMode = false;
        }
      }
      
      SphereList spheres = this.initDefSphereList();
      return createRadialState(spheres);
    }
    else {
      if (debugPrints) {
        println("InitFileException caught in initFile:");
        println(errorMessage);
      }
      return new SimulationState();
    }
  }
  
  void initFoldSteps() {
    ArrayList<int[]> pairs = new ArrayList<int[]>();
    for (int i = 0; i < pairTable.length; ++i) {
      int pairId = pairTable[i];
      if (pairId > i) {
        int[] pair = new int[2];
        pair[0] = i;
        pair[1] = pairId;
        pairs.add(pair);
      }
    }
    foldSteps = pairs;
  }
  
  void initTables() {
    initFoldSteps();
    strandOrder = orderStrands();
    residueOrder = initOrder(this);
    generatePairTableOptimal();
  }
  
  SphereList initDefSphereList() { // Outputs the base spherelist without coordinates
    SphereList sl = new SphereList();
    for (int i = 0; i < lenTot; ++i) {
      sl.add(new Sphere2D(i, seqIds[i], seqTot.substring(i, i+1), new String[lenTot]));
      
      if ((i > 0) && (seqIds[i-1] == seqIds[i])) {  // connect them if they are in the same strand
        sl.get(i-1).threeP = sl.get(i);
        sl.get(i).fiveP = sl.get(i-1);
      }
    }
    
    // Initialize relative IDs
    int currStrandId = sl.get(0).strandId;
    int idOffset = 0;
    for (int i = 0; i < lenTot; ++i) {
      int strandId = sl.get(i).strandId;
      if (strandId != currStrandId) {
        idOffset = i;
        currStrandId = strandId;
      }
      int relId = i - idOffset;
      sl.get(i).relId = relId;
    }
    
    // Initialize bond arrays based on pair table and nonCanonicals
    for (int i = 0; i < pairTable.length; ++i) {
      int pairId = pairTable[i];
      if (pairId > i) {
        String type = "cWW";
        if (nonCanonicals.containsKey(i + " " + pairId)) {
          type = nonCanonicals.get(i + " " + pairId);
        }
        setForcePair(sl.get(i), sl.get(pairId), type, this);
      }
    }
    for (Map.Entry<String, String> nc : nonCanonicals.entrySet()) {
      String[] pair = nc.getKey().split(" ");
      int id1 = int(pair[0]);
      int id2 = int(pair[1]);
      String bondType = nc.getValue();
      sl.get(id1).bonds[id2] = bondType;
      sl.get(id2).bonds[id1] = bondType;
    }

    return sl;
  }
  
  
  /** Parses ct or bpseq file. */
  void initColumnFile(String fileName, String fileType) {
    int numColumns = 6;
    int pairedCol = 4;
    if ("bpseq".equals(fileType)) {
      numColumns = 3;
      pairedCol = 2;
    }
    String[] lines;
    
    if (!BROWSER) {
      lines = loadStrings(fileName);
    }
    else {
      if (javascript != null) {
        String inputText = javascript.getTextInput("inputtext");
        lines = inputText.split("\\r\\n|\\n|\\r");
      }
      else {
        if (debugPrints) { println("ERROR: JAVASCRIPT NOT INITIALIZED!"); }
        errorMessage = "ERROR: JAVASCRIPT NOT INITIALIZED!";
        return;
      }
    }
    
    // Take out extra whitespace
    for (int i = 0; i < lines.length; ++i) {
      lines[i] = lines[i].trim().replaceAll("\\s+", " ");
    }
    
    // First line not preceded by # is header
    String[] header = null;
    int headerIndex = -1;
    for (int i = 0; i < lines.length; ++i) {
      if (!lines[i].equals("") && !lines[i].substring(0, 1).equals("#")) {
        header = split(lines[i], " ");
        headerIndex = i;
        break;
      }
    }
    if (header == null) {
      if (debugPrints) { println("No valid lines"); }
      errorMessage = "No valid lines!";
      return;
    }
    // If the first line is data
    if (header.length == numColumns && header[0].matches("^\\d+$") && int(header[0]) == 1 && header[1].matches("[a-zA-Z]") && header[pairedCol].matches("^\\d+$")) {
      if ( "bpseq".equals(fileType) || ("ct".equals(fileType) && header[2].matches("^\\d+$") && header[3].matches("^\\d+$") && header[5].matches("^\\d+$")) ) {
        header = null;
        headerIndex--;
      }
    }
    if (debugPrints) { println("headerIndex: " + headerIndex); }
    
    int numSeqs = 1; // Number of strands
    starts = new int[1]; // Array of start index of each strand (starting on 0)
    starts[0] = 0;
    sequences = new String[1];
    int dataLength = lines.length - (headerIndex + 1);
    if (debugPrints) { println("dataLength: " + dataLength); }
    boolean multipleStrands = false;
    
    // Check to see if it is in Shapiro Lab header format
    if (header != null && header.length > 3) {
      // Check if all arguments are numbers
      boolean flag = false;
      for (String argument : header) {
        if (!(argument.matches("^\\d+$"))) {
          flag = true;
          break;
        }
      }
      // Check if second argument equals number of strand start arguments, and first strand starts on 1
      if ((!flag) && int(header[1]) == (header.length - 2) && int(header[2]) == 1) {
        numSeqs = int(header[1]);
        starts = new int[numSeqs];
        starts[0] = 0;
        
        // Check if strand start arguments are in ascending order
        flag = false;
        for (int i = 1; i < numSeqs; ++i) { // for each strand
          if (int(header[i+2]) > int(header[i+1])) {
            starts[i] = int(header[i+2]) - 1; // set start
          } else {
            flag = true;
            if (debugPrints) { println("The strand start indexes (arguments 3 and up in header) must be in ascending order!"); }
            break;
          }
        }
        
        // Check if first and last arguments are in-bounds and compatible with the other arguments
        if (!flag) {
          int arg1 = int(header[0]);
          if (arg1 >= 2 && (arg1 <= dataLength) && starts[starts.length-1] < arg1 && arg1 >= numSeqs) {
            multipleStrands = true;
            sequences = new String[numSeqs];
          }
        }
      }
    }
    if (debugPrints && multipleStrands) { println("multipleStrands = true"); }
    
    if (!multipleStrands) {
      numSeqs = 1;
      starts = new int[1];
      starts[0] = 0;
      sequences = new String[1];
    }
    
    /**
    // Sets lenInHeader to header value if it is valid. If not, sets to default if not multiple strands, or throws error if there are multiple strands.
    int dataLength;
    if (header != null && header[0].matches("^\\d+$")) {
      dataLength = int(header[0]);
      if (dataLength < 1 || dataLength >= foldSteps.length - headerIndex) {
        if (debugPrints) { println("Header out of bounds"); }
        if (multipleStrands) {
          errorMessage = "\"Total length\" argument (first # in the header) is out of bounds!";
          return;
        } else {
          dataLength = foldSteps.length - (headerIndex + 1);
        }
      }
    } else {
      dataLength = foldSteps.length - (headerIndex + 1);
    }
    
    if (multipleStrands && starts[starts.length-1] >= dataLength) {
      if (debugPrints) {
        println("Last strand start (last argument in header) is greater than the length of the sequence (first argument), which is impossible!");
      }
      multipleStrands = false;
    }

    if (multipleStrands && header != null) {
      numSeqs = int(header[1]);
      
      // Check to see if numSeqs matches second argument
      int validNums = 2;
      for (int i = 4; i < header.length; i++) {
        if (!header[i].matches("^\\d+$")) { break; }
        validNums++;
      }
      if (numSeqs != validNums) {
        if (debugPrints) {
          println("\"Number of strands\" argument (2nd number in the file header) is not equal to the amount of strand starts provided!");
        }
        errorMessage = "\"Number of strands\" (2nd argument in the file header) is not equal to the amount of strand start indexes provided!\n\n";
        return;
      }
      if (int(header[2]) != 1) {
        if (debugPrints) { println("The first strand start index (3rd argument in the file header) must be \"1\" by nature!"); }
        errorMessage = "The first strand start index (3rd argument in the file header) must be \"1\" by nature!";
        return;
      }
    } 
    
    if (multipleStrands) {
      if (debugPrints) { println("The number of strands is " + str(numSeqs)); }
      
      // Check to see if last start is in bounds
      if (starts[starts.length-1] >= dataLength) {
        if (debugPrints) {
          println("Last strand start (last argument in header) is greater than the length of the sequence (first argument), which is impossible!");
        }
        errorMessage = "Last strand start (last argument in header) is greater than the length of the sequence (first argument), which is impossible!";
        return;
      }
    }
    */
    
    String[] columns;
    seqTot = ""; // String of all bases
    ArrayList<Short> pairedIdCol = new ArrayList<Short>();
    
    // Native header format
    if ("bpseq".equals(fileType) || multipleStrands) {
      int index = headerIndex + 1; // index in actual file
      boolean endOfData = false;
      for (int i = 0; i < numSeqs; ++i) { // for each strand
        if (endOfData) { break; }
        int endIndex; // end id for strand in the strucutre 
        if (i+1 < starts.length) {
          endIndex = starts[i+1];
        } else { // Last strand
          endIndex = dataLength;
        }
        sequences[i] = "";
        
        for (int j = starts[i]; j < endIndex; index++) { // index is actual index, j is the # of found valid lines
          if (index >= lines.length) { // out of bounds, so break
            if (debugPrints) { println("End of file"); }
            break;
          }
          
          if (lines[index].equals("") || lines[index].substring(0, 1).equals("#")) { // Ignore empty lines and comments
            continue;
          }
          
          columns = split(lines[index], " ");
          if (columns.length != numColumns || !columns[0].matches("^\\d+$") || (int(columns[0]) != j+1) || !columns[1].matches("[a-zA-Z]") || columns[1].length() != 1 || !columns[pairedCol].matches("^\\d+$") || int(columns[pairedCol]) > 32765) {
            if (j < dataLength) {
              if (debugPrints) { println("CAUTION, INVALID LINE: " + str(index+1)); }
              warningMessage += "Caution: File line " + str(index + 1) + " (info for base " + str(j + 1) + ") is invalid.\n";
            }
            endOfData = true;
            break;
          }
          
          j++; // found valid line
          seqTot += columns[1];
          sequences[i] += columns[1];
          int pairedId = int(columns[pairedCol]);
          
          pairedIdCol.add( (short) (pairedId - 1) );
        }
      }
    }
    
    // Looking for 0s in the 3rd column, which some labs use to mark a new strand start.
    else {
      if (debugPrints) { println("Looking for strand starts as 0s in the 3rd column of ct"); }
      ArrayList<String> sequencesList = new ArrayList<String>();
      ArrayList<Integer> startsList = new ArrayList<Integer>();
      int baseCount = 0;
      String currentSeq = "";
      
      for (int i = headerIndex+1; i < lines.length && baseCount < dataLength; i++) {
        if (lines[i].equals("") || lines[i].substring(0, 1).equals("#")) { // Ignore empty lines and comments
          continue;
        }
        
        columns = split(lines[i], " ");
        if (columns.length != numColumns || !columns[0].matches("^\\d+$") || (int(columns[0]) != baseCount+1) || !columns[1].matches("[a-zA-Z]") || columns[1].length() != 1 || !columns[pairedCol].matches("^\\d+$") || int(columns[pairedCol]) > 32765 || !columns[2].matches("^\\d+$") || !columns[3].matches("^\\d+$")) {
          if (baseCount < dataLength) {
            if (debugPrints) { println("CAUTION: INVALID LINE: " + str(i+1)); }
            warningMessage += "Caution: File line " + str(i+1) + " (info for base " + str(seqTot.length() + 1) + ") is invalid.\n";
          }
          break;
        }
        
        pairedIdCol.add( (short) (int(columns[pairedCol]) - 1) );
        
        if ("0".equals(columns[2])) {
          startsList.add(baseCount);
          if (baseCount != 0) {
            sequencesList.add(currentSeq);
            currentSeq = "";
          }
        }
        
        currentSeq += columns[1];
        seqTot += columns[1];
        
        baseCount++;
      }
      sequencesList.add(currentSeq);
      
      // copy to starts
      starts = new int[startsList.size()];
      Iterator<Integer> iterator = startsList.iterator();
      for (int i = 0; i < starts.length; i++) {
        starts[i] = iterator.next(); // Removed .intValue()
      }
      
      sequences = sequencesList.toArray(new String[sequencesList.size()]); // copy sequences
      
      numSeqs = starts.length;
    }
    
    pairTable = pairedIdCol.toArray(new Short[pairedIdCol.size()]); // copy
    
    if (debugPrints) { println("Number of pairs: " + lines.length); }
    
    lenTot = seqTot.length();
    
    if (lenTot <= 0 || starts.length == 0) {
      errorMessage = "No valid lines found!";
      return;
    }
    
    if ( starts[starts.length-1] >= lenTot || (multipleStrands && header != null && lenTot != int(header[0])) ) {
      if (debugPrints) { println("First arg incorrect"); }
      warningMessage += "Caution: The first argument in the file header (total sequence length) does not equal the # of valid lines.\n\n";
      numSeqs = 1;
      starts = new int[1];
      starts[0] = 0;
      sequences = new String[1];
      sequences[0] = seqTot;
      multipleStrands = false;
    }
    
    checkPairTable(pairTable, lenTot);
    if (!"".equals(errorMessage)) { return; }
    
    // Creates int array seqIds that stores the strand number (starting at 0) for each residue starting at index 0
    seqIds = new int[lenTot];
    int sId = 0; // Strand ID
    int pc = 0; // residue index
    for (int i = 0; i < numSeqs; ++i) { // each strand
      for (int j = 0; j < sequences[i].length(); ++j) { // for the length of each strand
        seqIds[pc] = sId;
        pc++;
      }
      sId++;
    }
  }
  
  
  /** Parses dot bracket notation file */
  void initBracketFile(String fileName) {
    String[] lines;
  
    if (!BROWSER) {
      lines = loadStrings(fileName);
    }
    else {
      if (javascript != null) {
        String inputText = javascript.getTextInput("inputtext");
        lines = inputText.split("\\r\\n|\\n|\\r");
      }
      else {
        if (debugPrints) { println("ERROR: JAVASCRIPT NOT INITIALIZED!"); }
        errorMessage = "ERROR: JAVASCRIPT NOT INITIALIZED!";
        return;
      }
    }
    
    if (debugPrints) {
      println("Raw input:");
      for (String s : lines) {
        println(s);
      }
    }
    
    // Take out extra whitespace
    for (int i = 0; i < lines.length; ++i) {
      lines[i] = lines[i].trim().replaceAll("\\s+", " ");
    }
    
    int sequenceIndex = -1;
    int structureIndex = -1;
    
    for (int i = 0; i < lines.length; ++i) {
      String line = lines[i];
      if (line.length() > 0) {
        String seqChar = line.substring(0, 1);
        if ( (!BROWSER && line.matches("[a-zA-Z& ]+")) || (BROWSER && !(seqChar.equals("#") || seqChar.equals(" ") || seqChar.equals(">"))) ) {
          sequenceIndex = i;
          for (int j = i+1; j < lines.length; ++j) {
            if (lines[j].length() < 1) { continue; }
            String startChar = lines[j].substring(0,1);
            if ( !(startChar.equals("#") || startChar.equals(" ")) ) {
              structureIndex = j;
            }
          }
          break;
        }
      }
    }
    
    if (debugPrints) {
      println("sequenceIndex: " + sequenceIndex);
      println("structureIndex: " + structureIndex);
    }
  
    // Assert both indices in bounds
    if (sequenceIndex < 0 || sequenceIndex >= lines.length-1 || structureIndex < 1 || structureIndex >= lines.length) {
      errorMessage = "Could not read input \"dot bracket\" file";
      return;
    }
    
    String sequenceInput = lines[sequenceIndex];
    String structureInput = lines[structureIndex];
    
    ArrayList<Integer> startsList = new ArrayList<Integer>();
    startsList.add(0);
  
    // Build strings without spaces or ampersands
    seqTot = "";
    for (int i = 0, l = 0; i < sequenceInput.length(); ++i) {
      String currentChar = sequenceInput.substring(i, i+1);
      if ( !(" ".equals(currentChar) || "&".equals(currentChar)) ) {
        seqTot += currentChar;
        ++l;
        if (l > 1) {
          String prevChar = sequenceInput.substring(i-1, i);
          if (" ".equals(prevChar) || "&".equals(prevChar)) {
            startsList.add(l-1);
          }
        }
      }
    }
    
    // copy
    starts = new int[startsList.size()];
    Iterator<Integer> iterator = startsList.iterator();
    for (int i = 0; i < starts.length; i++) {
      starts[i] = iterator.next(); // Removed .intValue()
    }
    
    // Build structure string without spaces or ampersands
    String catStruct = "";
    for (int i = 0; i < structureInput.length(); ++i) {
      String currentChar = structureInput.substring(i, i+1);
      if ( !(" ".equals(currentChar) || "&".equals(currentChar)) ) {
        catStruct += currentChar;
      }
    }
  
    // Assert seqTot and catStruct are same length
    if (seqTot.length() != catStruct.length()) {
      if (debugPrints) { println("seqTot and catStruct not same length"); }
      errorMessage = "Sequence and structure input lines are not the same length";
      return;
    }
    
    lenTot = seqTot.length();
    
    if (lenTot > 32765) {
      if (debugPrints) { println("Input structure (" + lenTot + " nucleotides) is too large."); }
      errorMessage = "Input structure (" + lenTot + " nucleotides) is too large.";
      return;
    }
    
    pairTable = new Short[lenTot];
    
    String[] pairChars = {"()", "[]", "{}", "<>"};
    ArrayList<OpenPair> pairs = new ArrayList<OpenPair>();
  
    if (debugPrints) { println("\n***Starting structure building loop***"); }
  
    for (int i = 0; i < lenTot; ++i) { // For each open side of a helix
      String currentChar = catStruct.substring(i, i+1);
      boolean match = false;
      for (int j = 0; j < pairChars.length; ++j) {
        if (pairChars[j].substring(0, 1).equals(currentChar)) {
          match = true;
          break;
        }
      }
      if (match) { // Found 5' open end of a pair
        pairs.add(new OpenPair(i, currentChar));
        if (debugPrints) { println("Added new open pair " + i + ", " + currentChar); }
      }
      else {
        pairTable[i] = -1;
        String pairChar = "";
        for (int j = 0; j < pairChars.length; ++j) {
          if (pairChars[j].substring(1, 2).equals(currentChar)) {
            match = true;
            pairChar = pairChars[j].substring(0, 1);
            break;
          }
        }
        if (match) { // Found 3' closing end of a pair
          if (debugPrints) { println("Found closing " + i + ", " + currentChar); }
          for (int j = pairs.size()-1; j >= 0; --j) {
            OpenPair op = pairs.get(j);
            if ( pairChar.equals(op.openChar) ) {
              // Add bond
              pairTable[op.index] = (short) i;
              pairTable[i] = (short) op.index;
              if (debugPrints) { println("Adding bond " + i + ", " + op.index); }
              pairs.remove(j);
              break;
            }
          }
          if (pairTable[i] == -1) {
            if (debugPrints) { println("Unpaired " + currentChar + " at position " + str(i+1)); }
            errorMessage = "Unpaired \"" + currentChar + "\" at position " + str(i+1) + " in structure file";
            return;
          }
        }
      }
    }
    
    checkPairTable(pairTable, lenTot);
    if (!"".equals(errorMessage)) { return; }
    
    // Initialize sequences
    int numSeqs = starts.length;
    sequences = new String[numSeqs];
    for (int i = 0, l = numSeqs-1; i < l; ++i) {
      sequences[i] = seqTot.substring(starts[i], starts[i+1]);
    }
    sequences[numSeqs-1] = seqTot.substring(starts[numSeqs-1], lenTot);
    
    // Creates int array seqIds that stores the strand number (starting at 0) for each residue starting at index 0
    seqIds = new int[lenTot];
    int sId = 0; // Strand ID
    int pc = 0; // residue index
    for (int i = 0; i < numSeqs; ++i) { // each strand
      for (int j = 0; j < sequences[i].length(); ++j) { // for the length of each strand
        seqIds[pc] = sId;
        pc++;
      }
      sId++;
    }
  }

  /** Parses sequence file (code is derived from parser of dot bracket notation file without the bracket part */
  void initSequenceFile(String fileName) {
    String[] lines;
  
    if (!BROWSER) {
      lines = loadStrings(fileName);
    }
    else {
      if (javascript != null) {
        String inputText = javascript.getTextInput("inputtext");
        lines = inputText.split("\\r\\n|\\n|\\r");
      }
      else {
        if (debugPrints) { println("ERROR: JAVASCRIPT NOT INITIALIZED!"); }
        errorMessage = "ERROR: JAVASCRIPT NOT INITIALIZED!";
        return;
      }
    }
    
    if (debugPrints) {
      println("Raw input:");
      for (String s : lines) {
        println(s);
      }
    }
    
    // Take out extra whitespace
    for (int i = 0; i < lines.length; ++i) {
      lines[i] = lines[i].trim().replaceAll("\\s+", " ");
    }
    
    int sequenceIndex = -1;
    int structureIndex = -1;
    
    for (int i = 0; i < lines.length; ++i) {
      String line = lines[i];
      if (line.length() > 0) {
        String seqChar = line.substring(0, 1);
        if ( (!BROWSER && line.matches("[a-zA-Z& ]+")) || (BROWSER && !(seqChar.equals("#") || seqChar.equals(" ") || seqChar.equals(">"))) ) {
          sequenceIndex = i;
          for (int j = i+1; j < lines.length; ++j) {
            if (lines[j].length() < 1) { continue; }
            String startChar = lines[j].substring(0,1);
            if ( !(startChar.equals("#") || startChar.equals(" ")) ) {
              structureIndex = j;
            }
          }
          break;
        }
      }
    }
    
    if (debugPrints) {
      println("sequenceIndex: " + sequenceIndex);
      println("structureIndex: " + structureIndex);
    }
  
    // Assert both indices in bounds
    if (sequenceIndex < 0 || sequenceIndex >= lines.length-1) {
      errorMessage = "Could not read input \"dot bracket\" file";
      return;
    }
    
    String sequenceInput = lines[sequenceIndex];
    
    ArrayList<Integer> startsList = new ArrayList<Integer>();
    startsList.add(0);
  
    // Build strings without spaces or ampersands
    seqTot = "";
    for (int i = 0, l = 0; i < sequenceInput.length(); ++i) {
      String currentChar = sequenceInput.substring(i, i+1);
      if ( !(" ".equals(currentChar) || "&".equals(currentChar)) ) {
        seqTot += currentChar;
        ++l;
        if (l > 1) {
          String prevChar = sequenceInput.substring(i-1, i);
          if (" ".equals(prevChar) || "&".equals(prevChar)) {
            startsList.add(l-1);
          }
        }
      }
    }
    
    // copy
    starts = new int[startsList.size()];
    Iterator<Integer> iterator = startsList.iterator();
    for (int i = 0; i < starts.length; i++) {
      starts[i] = iterator.next(); // Removed .intValue()
    }
        
    lenTot = seqTot.length();
    
    if (lenTot > 32765) {
      if (debugPrints) { println("Input structure (" + lenTot + " nucleotides) is too large."); }
      errorMessage = "Input structure (" + lenTot + " nucleotides) is too large.";
      return;
    }
    
    if (!"".equals(errorMessage)) { return; }
    
    // Initialize sequences
    int numSeqs = starts.length;
    sequences = new String[numSeqs];
    for (int i = 0, l = numSeqs-1; i < l; ++i) {
      sequences[i] = seqTot.substring(starts[i], starts[i+1]);
    }
    sequences[numSeqs-1] = seqTot.substring(starts[numSeqs-1], lenTot);
    
    // Creates int array seqIds that stores the strand number (starting at 0) for each residue starting at index 0
    seqIds = new int[lenTot];
    int sId = 0; // Strand ID
    int pc = 0; // residue index
    for (int i = 0; i < numSeqs; ++i) { // each strand
      for (int j = 0; j < sequences[i].length(); ++j) { // for the length of each strand
        seqIds[pc] = sId;
        pc++;
      }
      sId++;
    }
  }
  
  /** Helper data structure for bracket parser */
  class OpenPair {
    int index = -1;
    String openChar = "";
  
    OpenPair(int id, String oc) {
      this.index = id;
      this.openChar = oc;
    }
  }
  
  
  /** Parse native nucleotide pairs format file */
  void initNtsPairs(String fileName) {
    String[] lines;
  
    if (!BROWSER) {
      lines = loadStrings(fileName);
    } else {
      if (javascript != null) {
        String inputText = javascript.getTextInput("inputtext");
        lines = inputText.split("\\r\\n|\\n|\\r");
      } else {
        if (debugPrints) { println("ERROR: JAVASCRIPT NOT INITIALIZED!"); }
        errorMessage = "ERROR: JAVASCRIPT NOT INITIALIZED!";
        return;
      }
    }
    
    // Take out extra whitespace
    for (int i = 0; i < lines.length; ++i) {
      lines[i] = lines[i].trim().replaceAll("\\s+", " ");
    }
    
    int ntsIndex = -1;
    int pairsIndex = -1;
    
    for (int i = 0; i < lines.length; ++i) {
      String line = lines[i];
      if ("NTS:".equals(line) || "nts:".equals(line) || "NUCLEOTIDES:".equals(line) || "nucleotides:".equals(line) || "Nucleotides:".equals(line)) {
        ntsIndex = i;
      }
      else if ("PAIRS:".equals(line) || "pairs:".equals(line) || "Pairs:".equals(line)) {
        pairsIndex = i;
      }
      
      if (ntsIndex != -1 && pairsIndex != -1) {
        break;
      }
    }
    
    if (ntsIndex < 0 || ntsIndex >= lines.length || pairsIndex < 0 || pairsIndex >= lines.length || ntsIndex == pairsIndex) {
      if (debugPrints) { println("Could not read Nucleotides / Pairs file!"); }
      errorMessage = "Could not read Nucleotides / Pairs file!";
      return;
    }
    
    int ntsEndIndex, pairsEndIndex;
    
    if (ntsIndex < pairsIndex) {
      ntsEndIndex = pairsIndex -1;
      pairsEndIndex = lines.length;
    }
    else {
      pairsEndIndex = ntsIndex -1;
      ntsEndIndex = lines.length;
    }
    
    if (debugPrints) {
      println("Starts: Nts=" + ntsIndex + ", pairs=" + pairsIndex);
      println("Ends: Nts=" + ntsEndIndex + ", pairs=" + pairsEndIndex);
    }
    
    HashMap<String, Short> ntCodes = new HashMap<String, Short>(); // Key: string code. Value: index of residue
    ArrayList<Integer> startsList = new ArrayList<Integer>();
    
    seqTot = "";
    
    for (int i = ntsIndex+2, count = 0; i < ntsEndIndex; ++i) {
      String[] line = lines[i].split(" ");
      if ( line.length == 0 || line[0].trim().length() == 0 || "#".equals(line[0].trim().substring(0, 1)) ) {
        continue;
      }
      
      // Check if it is data. else, break
      if ( line.length != 8 || !line[2].matches("^\\d+$") || line[6].length() != 3 ) {
        break;
      }
      
      if ("1".equals(line[2])) {
        startsList.add(count);
      }
      seqTot += line[6].substring(1, 2);
      ntCodes.put( line[7].substring(1, line[7].length()), (short) count );
      
      ++count;
    }
    
    // copy to starts
    starts = new int[startsList.size()];
    Iterator<Integer> iterator = startsList.iterator();
    for (int i = 0; i < starts.length; i++) {
      starts[i] = iterator.next(); // Removed .intValue()
    }
    
    if (starts.length == 0 || starts[0] != 0) {
      if (debugPrints) { println("Strand indecies in Nucleotides / Pairs file not formatted correctly!"); }
      errorMessage = "Strand indecies in Nucleotides / Pairs file not formatted correctly!";
      return;
    }
    
    lenTot = seqTot.length();
    
    if (lenTot > 32765) {
      if (debugPrints) { println("Input structure (" + lenTot + " nucleotides) is too large."); }
      errorMessage = "Input structure (" + lenTot + " nucleotides) is too large.";
      return;
    }
    
    
    int numSeqs = starts.length;
    sequences = new String[numSeqs];
    for (int i = 0, l = numSeqs-1; i < l; ++i) {
      sequences[i] = seqTot.substring(starts[i], starts[i+1]);
    }
    sequences[numSeqs-1] = seqTot.substring(starts[numSeqs-1], lenTot);
    
    pairTable = new Short[lenTot];
    for (int i = 0, l = pairTable.length; i < l; ++i) {
      pairTable[i] = -1;
    }
    
    for (int i = pairsIndex+2; i < pairsEndIndex; ++i) {
      String[] line = lines[i].split(" ");
      if ( line.length == 0 || "#".equals(lines[i].trim().substring(0, 1)) ) {
        continue;
      }
      
      if ( line.length != 9 || line[7].length() != 5 ) {
        break;
      }
      
      String ntCode1 = line[2].substring( 1, line[2].length() );
      String ntCode2 = line[3].substring( 1, line[3].length() );
      
      int id1 = ntCodes.get(ntCode1);
      int id2 = ntCodes.get(ntCode2);
      
      if (id1 > id2) {
        int tempId;
        tempId = id1;
        id1 = id2;
        id2 = tempId;
      }
      
      boolean firstBond = false;
      if (pairTable[id1] == -1 && pairTable[id2] == -1) {
        pairTable[id1] = (short)id2;
        pairTable[id2] = (short)id1;
        firstBond = true;
      }
      if (!"\"cWW\"".equals(line[7]) || !firstBond) {
        nonCanonicals.put(id1 + " " + id2, line[7].substring(1, 4));
      }
    }
    
    checkPairTable(pairTable, lenTot);
    if (!"".equals(errorMessage)) { return; }
    
    // Creates int array seqIds that stores the strand number (starting at 0) for each residue starting at index 0
    seqIds = new int[lenTot];
    int sId = 0; // Strand ID
    int pc = 0; // residue index
    for (int i = 0; i < numSeqs; ++i) { // each strand
      for (int j = 0; j < sequences[i].length(); ++j) { // for the length of each strand
        seqIds[pc] = sId;
        pc++;
      }
      sId++;
    }
    
    // Initialize foldSteps
    /*
    ArrayList<String> foldStepsList = new ArrayList<String>();
    for (int i = 0; i < pairTable.length; ++i) {
      if (i < pairTable[i]) {
        foldStepsList.add("F " + str(i+1) + " " + str(pairTable[i]+1));
      }
    }
    foldSteps = foldStepsList.toArray(new String[foldStepsList.size()]); // copy
    */
  }
  
  
  /** Parses native save file format, and sets up state. */
  SimulationState initSaveFile(String loadFileName) {
    SimulationState simState = new SimulationState();
    SphereList spheres = new SphereList();
    simState.sim = new SimulationSetup();
    SimulationSetup s = simState.sim;
    ForceField tempff = new ForceField();
    DisplaySettings tempds = new DisplaySettings();
    
    //boolean tempRigidHelices;
    //boolean tempRigidHairpins;
    
    boolean noRigidLoopsEntry = true;
    
    String[] coordinates = null;
    
    ArrayList<Double> colorData = new ArrayList<Double>();
    
    String[] input;
    if (!BROWSER) {
      input = loadStrings(loadFileName);
    }
    else {
      if (javascript != null) {
        String inputText = javascript.getTextInput("inputtext");
        input = inputText.split("\\r\\n|\\n|\\r");
      }
      else {
        if (debugPrints) { println("ERROR: JAVASCRIPT NOT INITIALIZED!"); }
        errorMessage = "ERROR: JAVASCRIPT NOT INITIALIZED!";
        return new SimulationState();
      }
    }
    
    for (int i = 0; i < input.length; ++i) {
      if (input[i].equals("") || input[i].substring(0, 1).equals("#")) { // Ignore empty lines and comments
        continue;
      }
      String[] words = input[i].split("=");
      if (words.length >= 2) {
        if (words[0].equals("sequences")) {
          s.sequences = words[1].split("; ");
        }
        else if (words[0].equals("starts")) {
          String[] words2 = words[1].split(",");
          s.starts = new int[words2.length];
          for (int j = 0; j < words2.length; j++) {
            if (words2[j].matches("^\\d+$")) {
              s.starts[j] = int(words2[j]);
            }
            else {
              errorMessage = "Start " + (j+1) + " is not a valid number!";
              return new SimulationState();
            }
          }
        }
        else if (words[0].equals("positions")) {
          coordinates = words[1].split("; ");
        }
        else if (words[0].equals("radius")) {
          if (words[1].matches("[-+]?\\d*\\.?\\d*")) {
            tempff.radius = stringToDouble(words[1]);
            
            if (tempff.radius < tempff.minRadius) {
              tempff.radius = tempff.minRadius;
            } else if (tempff.radius > tempff.minRadius + tempff.radiusScale) {
              tempff.radius = tempff.minRadius + tempff.radiusScale;
            }
          }
          else {
            errorMessage = "Radius is not a valid number!";
            return new SimulationState();
          }
        }
        else if (words[0].equals("backbone_distance")) {
          if (words[1].matches("[-+]?\\d*\\.?\\d*")) {
            tempff.backboneDist = stringToDouble(words[1]);
            
            if (tempff.backboneDist < tempff.minBackboneDist) {
              tempff.backboneDist = tempff.minBackboneDist;
            } else if (tempff.backboneDist > (tempff.minBackboneDist + tempff.backboneDistScale)) {
              tempff.backboneDist = tempff.minBackboneDist + tempff.backboneDistScale;
            }
          }
          else {
            errorMessage = "Backbone distance is not a valid number!";
            return new SimulationState();
          }
        }
        else if (words[0].equals("basepair_distance")) {
          if (words[1].matches("[-+]?\\d*\\.?\\d*")) {
            tempff.bpDist = stringToDouble(words[1]);
            
            if (tempff.bpDist < tempff.minBpDist) {
              tempff.bpDist = tempff.minBpDist;
            } else if (tempff.bpDist > (tempff.minBpDist + tempff.bpDistScale)) {
              tempff.bpDist = tempff.minBpDist + tempff.bpDistScale;
            }
          }
          else {
            errorMessage = "Basepair distance is not a valid number!";
            return new SimulationState();
          }
        }
        /*
        else if (words[0].equals("padding")) {
          if (words[1].matches("[-+]?\\d*\\.?\\d*")) {
            tempds.padding = stringToDouble(words[1]);
            
            if (tempds.padding < 0) {
              tempds.padding = 0;
            } else if (tempds.padding > tempds.paddingScale) {
              tempds.padding = tempds.paddingScale;
            }
          }
          else {
            errorMessage = "Padding is not a valid number!";
            return new SimulationState();
          }
        }
        */
        else if (words[0].equals("color_scheme")) {
          if ( words[1].matches("^\\d+$") && int(words[1]) < DisplaySettings.COLOR_MODE_MAX ) {
            tempds.colorMode = int(words[1]);
          }
        }
        else if (words[0].equals("outline")) {
          tempds.outlineMode = stringToBoolean(words[1]);
        }
        else if (words[0].equals("movement")) {
          simulationMode = stringToBoolean(words[1]);
        }
        else if (words[0].equals("simulation_type")) {
          simulateAllMode = stringToBoolean(words[1]);
        }
        else if (words[0].equals("labels")) {
          tempds.labelMode = stringToBoolean(words[1]);
        }
        else if (words[0].equals("rigid_helices")) {
          rigidHelices = stringToBoolean(words[1]);
        }
        else if (words[0].equals("rigid_hairpins")) {
          rigidHairpins = stringToBoolean(words[1]);
          wasRigidHairpins = rigidHairpins;
        }
        else if (words[0].equals("rigid_loops")) {
          noRigidLoopsEntry = false;
          rigidLoops = stringToBoolean(words[1]);
          wasRigidLoops = rigidLoops;
        }
        else if (words[0].equals("relaxed_ids")) {
          String[] ids = words[1].split(" ");
          for (String id : ids) {
            if (id.matches("^\\d+$")) {
              spheres.relaxedIds.put(int(id), 1);
            }
          }
        }
        else if (words[0].equals("display_order")) {
          String[] ids = words[1].split(",");
          for (String id : ids) {
            if (id.matches("^\\d+$")) {
              spheres.displayOrder.add(int(id));
            }
          }
        }
        else if (words[0].equals("pair_table")) {
          String[] words2 = words[1].split(",");
          s.pairTable = new Short[words2.length];
          for (int j = 0; j < words2.length; j++) {
            String pairedW = words2[j];
            if (pairedW.matches("^-?\\d+$")) {
              if (int(pairedW) < 32765) {
                s.pairTable[j] = (short) int(pairedW);
              }
              else {
                errorMessage = "pair_table value " + (j+1) + " is greater than the max value (32764)!";
                return new SimulationState();
              }
            }
            else {
              errorMessage = "pair_table value " + (j+1) + " is not a valid number";
              return new SimulationState();
            }
          }
        }
        else if (words[0].equals("non_canonicals")) {
          String[] words2 = words[1].split(", ");
          s.nonCanonicals = new HashMap<String, String>(words2.length);
          for (int j = 0; j < words2.length; ++j) {
            String word = words2[j];
            String[] entries = word.split(" ");
            if (entries.length == 3 && entries[0].matches("^\\d+$") && entries[1].matches("^\\d+$") && entries[2].length() == 3) {
              s.nonCanonicals.put(entries[0] + " " + entries[1], entries[2]);
            } else if (debugPrints) {
              println("Non Canonical entry number " + (j+1) + " in save file incorrectly formatted");
            }
          }
        }
        else if (words[0].equals("base_colors")) {
          String[] nums = words[1].split(",");
          for (String num : nums) {
            if (num.matches("-?\\d+(\\.\\d+)?")) {
              colorData.add(stringToDouble(num));
            }
            else {
              break;
            }
          }
        }
      }
    }
    
    // File corruption checks
    if (s.starts == null) {
      if (debugPrints) { println("Starts is null!"); }
      errorMessage = "Missing starts argument in load file";
      return new SimulationState();
    }
    else if (s.starts[0] != 0) {
      if (debugPrints) { println("First start is not 0!"); }
      errorMessage = "First sequence start is not 0 (cooresponding to first residue)!";
      return new SimulationState();
    }
    for (int i = 1; i < s.starts.length; i++) {
      if (s.starts[i] <= s.starts[i-1]) {
        if (debugPrints) { println("Strand starts not in ascending order!"); }
        errorMessage = "Strand starts not in ascending order!";
        return new SimulationState();
      }
    }
    
    if (s.sequences == null) {
      if (debugPrints) { println("Sequences is null!"); }
      errorMessage = "Missing sequences argument in load file";
      return new SimulationState();
    }
    else if (s.starts.length != s.sequences.length) {
      if (debugPrints) { println("Sequences and starts arrays not of same length!"); }
      errorMessage = "Sequences and starts arrays not of same length!";
      return new SimulationState();
    }
    for (int i = 1; i < s.starts.length; i++) {
      int size = s.starts[i] - s.starts[i-1];
      if (size != s.sequences[i-1].length()) {
        if (debugPrints) { println("Start " + (i-1) + " to " + i + " does not match sequence " + i + " length!"); }
        errorMessage = "Start " + (i-1) + " to " + i + " does not match sequence " + i + " length!";
        return new SimulationState();
      }
    }
    
    // Create seqTot, lenTot, and seqIds
    s.seqTot = "";
    for (String seq : s.sequences) {
      s.seqTot += seq;
    }
    s.lenTot = s.seqTot.length();
    
    s.seqIds = new int[s.lenTot];
    int sId = 0; // Strand ID
    int pc = 0; // residue index
    for (int i = 0; i < s.sequences.length; ++i) { // each strand
      for (int j = 0; j < s.sequences[i].length(); ++j) { // for the length of each strand
        s.seqIds[pc] = sId;
        pc++;
      }
      sId++;
    }
    
    if (s.pairTable == null) {
      if (debugPrints) { println("pairTable array is null!"); }
      errorMessage = "pair_table argument is missing from the load file!";
      return new SimulationState();
    }
    else if (s.pairTable.length != s.lenTot) {
      if (debugPrints) { println("pairTable length does not equal lenTot!"); }
      errorMessage = "Length of pair_table does not match the number of residues!";
      return new SimulationState();
    }
    checkPairTable(s.pairTable, s.lenTot);
    if (!"".equals(errorMessage)) { return new SimulationState(); }
    
    s.initTables();
    
    for (Map.Entry<String, String> nc : s.nonCanonicals.entrySet()) {
      String[] saveKey = nc.getKey().split(" ");
      int id1 = int(saveKey[0]);
      int id2 = int(saveKey[1]);
      if (id1 == id2) {
        if (debugPrints) {
          println("In save file, non canonical pair \"" + nc + "\" describes a base paired with itself, which is impossible");
          warningMessage += "In save file, non canonical pair \"" + nc + "\" describes a base paired with itself, which is impossible";
        }
      } else if (id1 < 0 || id1 >= s.lenTot || id2 < 0 || id2 >= s.lenTot) {
        println("In save file, non canonical pair \"" + nc + "\" is out of bounds!");
        warningMessage += "In save file, non canonical pair \"" + nc + "\" is out of bounds!";
      }
    }
    
    if (coordinates == null) {
      if (debugPrints) { println("coordinates array is null!"); }
      errorMessage = "coordinates of bases are missing from the load file!";
      return new SimulationState();
    }
    else if (coordinates.length != s.lenTot) {
      if (debugPrints) { println("Number of coordinate arguments does not match the number of residues!"); }
      errorMessage = "Number of coordinate arguments does not match the number of residues!";
      return new SimulationState();
    }
    
    if (spheres.displayOrder.size() > s.lenTot) {
      if (debugPrints) { println("displayOrder greater than lenTot!"); }
      errorMessage = "display_order is longer than the number of residues!";
      return new SimulationState();
    }
    else if (spheres.displayOrder.size() < s.lenTot) {
      if (debugPrints) { println("displayOrder less than lenTot!"); }
      errorMessage = "display_order is shorter than the number of residues!";
      return new SimulationState();
    }
    boolean[] numChecks = new boolean[s.lenTot];
    for (int i : spheres.displayOrder) {
      if (i < 0 || i >= s.lenTot) {
        if (debugPrints) { println("Value " + i + " in displayOrder is out of bounds!"); }
        errorMessage = "Value " + i + " in display_order is out of bounds!";
        return new SimulationState();
      }
      else if (numChecks[i]) {
        if (debugPrints) { println("Repeated value in displayOrder. A residue cannot have two positions in the order!"); }
        errorMessage = "Repeated value in display_order. A residue cannot have two positions in the order!";
        return new SimulationState();
      }
      
      numChecks[i] = true;
    }
    
    numChecks = new boolean[s.lenTot];
    for (Iterator<Integer> it = spheres.relaxedIds.keySet().iterator(); it.hasNext(); ) {
      int relaxedId = it.next();
      if (relaxedId < 0 || relaxedId >= s.lenTot) {
        if (debugPrints) { println("Relaxed ids are out of bounds!"); }
        errorMessage = "Relaxed ids are out of bounds!";
        return new SimulationState();
      }
      else if (numChecks[relaxedId]) {
        if (debugPrints) { println("Relaxed id " + relaxedId + " is repeated!"); }
        errorMessage = "Relaxed id " + relaxedId + " is repeated!";
        return new SimulationState();
      }
      
      numChecks[relaxedId] = true;
    }

    //movieMode = false;
    
    
    simState.spheres = s.initDefSphereList();
    SphereList newSpheres = simState.getSpheres();
    newSpheres.relaxedIds = spheres.relaxedIds;
    newSpheres.displayOrder = spheres.displayOrder;
    
    for (int i = 0; i < s.lenTot; ++i) {
      double x, y, z;
      String[] words = coordinates[i].split(",");
      if (words.length < 2) {
        if (debugPrints) { println("Not enough values for coordinate " + i + "!!!"); }
        // WARN MESSAGE?
        x = 0;
        y = 0;
        z = 0;
        
      } else {
        if ( words[0].matches("[-+]?\\d*\\.?\\d*") && words[1].matches("[-+]?\\d*\\.?\\d*") && (words.length == 2 || words[2].matches("[-+]?\\d*\\.?\\d*")) ) {
          x = stringToDouble(words[0]);
          y = stringToDouble(words[1]);
          
          if (words.length >= 3) {
            z = stringToDouble(words[2]);
          } else {
            z = 0;
          }
        }
        else {
          if (debugPrints) { println("Coordinate " + i + " has non-number values!"); }
          errorMessage = "Coordinate " + i + " has non-number values!";
          return new SimulationState();
        }
      }
      newSpheres.get(i).x = x;
      newSpheres.get(i).y = y;
      newSpheres.get(i).z = z;
    }
    
    if (colorData.size() == s.lenTot) {
      if (debugPrints) { println("Loading color data"); }
      for (int i = 0; i < s.lenTot; i++) {
        double colorNumber = colorData.get(i);
        newSpheres.get(i).colorNum = (float) colorNumber;
        newSpheres.get(i).baseColor = lerpColor(tempds.startGradient, tempds.endGradient, (float)colorNumber);
      }
    }
    
    ff = tempff;
    ds = tempds;
    
    if (noRigidLoopsEntry) {
      rigidLoops = false;
    }
    
    checkBoxes.get("Outlines").checked = ds.outlineMode;
    checkBoxes.get("Labels").checked = ds.labelMode;
    checkBoxes.get("Rigid Loops").checked = rigidLoops;
    checkBoxes.get("Rigid Hairpins").checked = rigidHairpins;
    checkBoxes.get("Sim. Selected Only").checked = !simulateAllMode;
    
    textDisplay = "SAVE FILE LOADED";
    
    return simState;
  }
  
  
  void checkPairTable(Short[] pairTable, int lenTot) {
    boolean[] pairCheck = new boolean[lenTot];
    for (short i = 0; i < pairTable.length; i++) {
      short pairId = pairTable[i];
      if (pairId < -1 || pairId >= lenTot) {
        if (debugPrints) { println("Residue " + (i+1) + " paired with an out-of-bounds number"); }
        errorMessage = "Residue " + (i+1) + " paired with an out of bounds number!";
        return;
      }
      
      if (pairId != -1) {
        if (pairId == i) {
          if (debugPrints) { println("Residue " + i + " paired with itself!"); }
          errorMessage = "Residue " + (i+1) + " cannot be paired with itself!";
          return;
        }
        else if (pairCheck[pairId]) {
          if (debugPrints) { println("Residue " + (pairId+1) + " is paired with multiple others!"); }
          errorMessage = "Residue " + (pairId+1) + " is paired with multiple other residues!";
          return;
        }
        else if (pairTable[pairId] != i) {
          if (debugPrints) {
            println("Residue " + (i+1) + " is paired with residue " + (pairId+1) + ", but residue " + (pairId+1) + " is paired with " + (pairTable[pairId]+1));
          }
          errorMessage = "Residue " + (i+1) + " is paired with residue " + (pairId+1) + ", but residue " + (pairId+1) + " is paired with " + (pairTable[pairId]+1);
          return;
        }
        
        pairCheck[pairId] = true;
      }
    }
  }
  
  /** Determine optimal order to place strands. */
  int[] orderStrands() {
    ArrayList<InterHelix> helices = new ArrayList<InterHelix>(); // stores interStrand helices
    int startStrandId = -1, endStrandId = -1;
    boolean foundHelix = false;
    int curHelSize = 0; // Number of pairs in last helix (default 1)
    int curStartStrand = -1, curEndStrand = -1;

    for (int i = 0, l = foldSteps.size(); i < l; ++i) {
      int start = foldSteps.get(i)[0];
      for (int j = 0; j < starts.length; j++) {
        if (start >= starts[j] && (j >= starts.length-1 || start < starts[j+1])) {
          startStrandId = j;
          break;
        }
      }
      int end = foldSteps.get(i)[1];
      for (int j = 0; j < starts.length; j++) {
        if (end >= starts[j] && (j >= starts.length-1 || end < starts[j+1])) {
          endStrandId = j;
          break;
        }
      }

      if (foundHelix) {
        InterHelix lastHelix = helices.get(helices.size() - 1);
        boolean sameHelix = (lastHelix.end1.origId + curHelSize == start && lastHelix.end2.origId - 1 == end);
        boolean connectedStrands = (curStartStrand == startStrandId || curEndStrand == endStrandId);
        if (sameHelix && connectedStrands) {
          curHelSize++;
          curStartStrand = startStrandId;
          curEndStrand = endStrandId;
          lastHelix.end2.origId = end;
          lastHelix.end2.strandId = endStrandId;
          lastHelix.end2.relId = end - starts[endStrandId];
        } else {
          foundHelix = false;
        }
      }
      if (!foundHelix) {
        if (startStrandId >= 0 && endStrandId >= 1 && startStrandId < endStrandId) { // Found an inter-strand helix
          helices.add( new InterHelix(startStrandId, start - starts[startStrandId], start, endStrandId, end - starts[endStrandId], end) );
          foundHelix = true;
          curHelSize = 1;
          curStartStrand = startStrandId;
          curEndStrand = endStrandId;
        }
      }
    }

    ArrayList<Strand> strands = new ArrayList<Strand>();
    for (int i = 0; i < sequences.length; i++) {
      int strandLength = 0;
      int strandStart = starts[i];
      if (i < sequences.length-1) {
        strandLength = starts[i+1] - strandStart;
      } else {
        strandLength = lenTot - strandStart;
      }
      strands.add(new Strand(i, strandLength));
    }

    // Add helix ends to respective strands
    for (InterHelix h : helices) {
      strands.get(h.end1.strandId).addEnd(h.end1);
      strands.get(h.end2.strandId).addEnd(h.end2);
    }

    for (int j = 0; j < sequences.length - 1; j++) {
      // Find next helix to place by finding the smallest possible distance to place.
      int smallestDistance = 32765;
      int bestHelixId = -1;
      boolean forward = true;
      for (int i = 0; i < helices.size(); i++) {
        InterHelix h = helices.get(i);
        if (h.end1.strandId == h.end2.strandId) {
          continue;
        }
        int strandLength1 = strands.get(h.end1.strandId).length;
        int strandLength2 = strands.get(h.end2.strandId).length;
        //assert(strandLength1 > h.end1.relId);
        //assert(strandLength2 > h.end2.relId);
        int distanceSum;
        distanceSum = (strandLength1 - h.end1.relId) + h.end2.relId + 1;
        if (distanceSum < smallestDistance) {
          smallestDistance = distanceSum;
          bestHelixId = i;
          forward = true;
        }
        distanceSum = (strandLength2 - h.end2.relId) + h.end1.relId + 1;
        if (distanceSum < smallestDistance) {
          smallestDistance = distanceSum;
          bestHelixId = i;
          forward = false;
        }
      }
      if (bestHelixId == -1) { // no more interstrand helices
        break;
      }  
      InterHelix bestHelix = helices.get(bestHelixId);
      Strand strand1 = strands.get(bestHelix.end1.strandId);
      Strand strand2 = strands.get(bestHelix.end2.strandId);
      if (!forward) {
        Strand temp = strand1;
        strand1 = strand2;
        strand2 = temp;
      }
      
      // Concatenate strands (Shift strand2's helices' relative ids and place them into strand1)
      for (HelixEnd end : strand2.helixEnds) {
        end.relId += strand1.length;
        end.strandId = strand1.id;
        strand1.helixEnds.add(end);
      }
      strand1.length += strand2.length;
      for (int id : strand2.strandIds) {
        strand1.strandIds.add(id);
      }
      
      // Clear strand 2
      strand2.length = 0;
      strand2.helixEnds.clear();
      strand2.strandIds.clear();
    }
    
    ArrayList<Integer> strandOrder = new ArrayList<Integer>();
    for (Strand s : strands) {
      strandOrder.addAll(s.strandIds);
    }
    int[] strandOrderOutput = new int[strandOrder.size()];
    if (debugPrints) { print("\nOrder:"); }
    for (int i = 0; i < strandOrder.size(); i++) {
      strandOrderOutput[i] = strandOrder.get(i);
      if (debugPrints) { print(" " + strandOrderOutput[i]); }
    }
    if (debugPrints) { println(""); }
    
    return strandOrderOutput;
  }
  
  
  /** Outputs an array storing the place in the new order (based on the strand order calculated from orderStrands) at each index cooresponding to the original id of each residue. */
  int[] initOrder(SimulationSetup s) {
    int[] residueOrderIds = new int[s.lenTot]; // Each index is the original id, stored value is order in index
    int orderId = 0;
    for (int i : s.strandOrder) {
      int strandLength = 0;
      int strandStart = s.starts[i];
      if (i < s.sequences.length-1) {
        strandLength = s.starts[i+1] - strandStart;
      } else {
        strandLength = s.lenTot - strandStart;
      }
      for (int j = strandStart; j < strandStart + strandLength; j++) {
        residueOrderIds[j] = orderId;
        orderId++;
      }
    }
    
    return residueOrderIds;
  }
  
  
  void generatePairTableOptimal() {
    int[] residueOrderSequence = new int[lenTot]; // Reverse of residueOrderIds: each index is the place in the new order, each stored value is the original id
    for (int i = 0; i < lenTot; i++) {
      residueOrderSequence[ residueOrder[i] ] = i;
    }
    
    pairTableOptimal = new ArrayList<Short>(pairTable.length); // Stores the position in the new order that each newly ordered residue is paired with
    
    for (int i = 0; i < pairTable.length; i++) { // Runs through the new order, building the new pairTableOptimal
      int pairedId = pairTable[residueOrderSequence[i]]; // The original ID of the residue paired with the residue in the new order at position i
      
      short pairedOrderedId; // The position in the new order of the residue paired with the residue in the new order at position i
      
      if (pairedId != -1) {
        pairedOrderedId = (short)residueOrder[ pairedId ];
      } else {
        pairedOrderedId = -1;
      }
      
      pairTableOptimal.add( pairedOrderedId );
    }
    
    ArrayList<Helix> helices = new ArrayList<Helix>();
    int end2 = -1;
    int size = 0;
    boolean inHelix = false;
    short pairedW;
    
    // Form helices
    for (int i = 0; i < pairTableOptimal.size(); ) {
      pairedW = pairTableOptimal.get(i);
      if (!inHelix) {
        if (i < pairedW) { // "i < pairedW" ensures each bond is only counted once and that pairedW is not -1 (meaning unpaired). WILL THIS BREAK IF A RESIDUE IS BONDED TO ITSELF?
          inHelix = true;
          size = 1;
          end2 = pairedW;
          helices.add(new Helix(i, end2, i));
        }
        i++;
      } else {
        end2--;
        if (end2 == pairedW) {
          i++;
          size++;
        } else { // end of helix
          helices.get(helices.size() - 1).setSize(size);
          inHelix = false;
        }
      }
    }
    
    for (int i = 0, l = helices.size(); i < l; i++) {
      Helix h1 = helices.get(i);
      for (int j = i+1; j < l; j++) {
        Helix h2 = helices.get(j);
        if ( (h1.ends[0] < h2.ends[0] && h2.ends[0] < h1.ends[1]) ^
             (h1.ends[0] < h2.ends[1] && h2.ends[1] < h1.ends[1]) ) {
          h1.crossedW.put(h2.id, h2);
          h2.crossedW.put(h1.id, h1);
          h1.numCrossed += h2.size;
          h2.numCrossed += h1.size;
        }
      }
    }
    
    ArrayList<Helix> ignoredHelices = new ArrayList<Helix>(); // The helices whose basepairs will be ignored in the radial placement algorithm, so the input is pseudoknot free.
    
    Helix mostCrossed;
    boolean finished = false;
    for (int i = 0; i < helices.size()-1; i++) {
      if (finished) { break; }
      mostCrossed = helices.get(0); // arbitrary start
      for (int j = 1; j < helices.size(); j++) {
        Helix nextH = helices.get(j);
        if (nextH.numCrossed > mostCrossed.numCrossed) {
          mostCrossed = nextH;
        }
      }
      
      if (mostCrossed.numCrossed < 0) {
        if (debugPrints) { println("NUM CROSSED OUT OF BOUNDS IN generatePairTableOptimal!"); }
      }
      
      if (mostCrossed.numCrossed == 0) { // No more pseudoknots!
        finished = true;
        break;
      }
      
      ignoredHelices.add(mostCrossed);
      
      for (int id : mostCrossed.crossedW.keySet()) {
        Helix h = mostCrossed.crossedW.get(id);
        h.crossedW.remove(mostCrossed.id);
        h.numCrossed -= mostCrossed.size;
      }
      
      mostCrossed.numCrossed = 0;
      mostCrossed.crossedW.clear();
    }
    
    
    // Modify pairTableOptimal to take out the basepairs in the ignored helices
    for (Helix h : ignoredHelices) {
      for (int i = 0; i < h.size; i++) {
        pairTableOptimal.set(h.ends[0]+i, (short) (-1));
        pairTableOptimal.set(h.ends[1]-i, (short) (-1));
      }
    }
  }
  
  void setupDebug() {
    println("\n*** SETUP ***");
    println("lenTot: " + lenTot);
    println("seqTot: " + seqTot);
    print("Starts:");
    for (int s : starts) {
      print(" " + s);
    }
    println("");
    
    print("Pair Table:");
    for (int i = 0; i < pairTable.length; ++i) {
      print(" " + pairTable[i]);
    }
    println("\n");
    
    println("foldSteps:");
    for (int[] pair : foldSteps) {
      println(pair[0] + " " + pair[1]);
    }
  }
  
  
  /** Sets up initial state as a circular diagram. */
  SimulationState createCircleState(SphereList spheres) {
    if (debugPrints) { println("\nStarting SimulationSetup.createCircleState..."); }
    
    SimulationState simState = new SimulationState();
    simState.sim = this; // sets the SimulationSetup object held in the SimulationState to be this SimulationSetup
    
    // Sets the spheres in a circle
    double circleReduction = 10.0;
    double startCircleRadius = (height / 2.0) - (ds.startPadding + (height / circleReduction));
    if (startCircleRadius < height / 8.0) {
      startCircleRadius = height / 8.0;
    }
    ff.radius = (Math.PI * startCircleRadius) / lenTot;
    if (ff.radius < ff.minRadius) { ff.radius = ff.minRadius; }
    if (ff.radius > ff.minRadius + ff.radiusScale) { ff.radius = ff.minRadius + ff.radiusScale; }
    
    ff.backboneDist = ff.radius * 2.5;
    if (ff.backboneDist < ff.minBackboneDist) { ff.backboneDist = ff.minBackboneDist; }
    if (ff.backboneDist > (ff.minBackboneDist + ff.backboneDistScale)) { ff.backboneDist = ff.minBackboneDist + ff.backboneDistScale; }
    
    ff.bpDist = ff.radius * 3.75;
    if (ff.bpDist < ff.minBpDist) { ff.bpDist = ff.minBpDist; }
    if (ff.bpDist > (ff.minBpDist + ff.bpDistScale)) { ff.bpDist = ff.minBpDist + ff.bpDistScale; }
    
    
    ds.padding = ds.startPadding;
    startCircleRadius += height / circleReduction;
    spheres.displayOrder.clear();
    spheres.relaxedIds.clear();
    for (int i = 0; i < lenTot; ++i) {
      double ang = residueOrder[i] * 2.0 * Math.PI / (double)(lenTot);
      spheres.get(i).x = (width/2) + Math.cos(ang) * startCircleRadius;
      spheres.get(i).y = (height/2) + Math.sin(ang) * startCircleRadius;
      spheres.displayOrder.add(i);
    }
    
    simulationMode = false;
    //movieMode = false;
    rigidHelices = true;
    
    simState.spheres = spheres;
    
    if (debugPrints) { println("Finished SimulationSetup.createCircleState!"); }
    return simState;
  }
  
  
  /** Sets up the state based on a radial algorithm. */
  SimulationState createRadialState(SphereList spheres) {
    if (debugPrints) { println("\nStarting SimulationSetup.createRadialState..."); }
    
    SimulationState simState = new SimulationState();
    simState.sim = this; // sets the SimulationSetup object held in the SimulationState to be this SimulationSetup
    
    ArrayList<Double> X = new ArrayList<Double>();
    ArrayList<Double> Y = new ArrayList<Double>();
    
    Radial radial = new Radial();
    radial.xyCoordinates(pairTableOptimal, X, Y); // Takes in a pseudoknot free pairTable of shorts.
    if (!"".equals(radialAlgorithmException)) {
      if (debugPrints) {
        println("Caught RadialAlgorithmException in createRadialState!");
        println(radialAlgorithmException);
      }
      errorMessage = radialAlgorithmException;
      return new SimulationState();
    }
    
    if (X.size() == Y.size()) {
      //println("Both coordinate arrays sized at " + X.size());
    } else {
      if (debugPrints) { println("Coordinate arrays of unequal size!"); }
      errorMessage = "Coordinate arrays of unequal size";
      return new SimulationState();
    }
    
    
    // Determine initial layout based on screen size and radial structure size
    double minX = X.get(0), maxX = X.get(0);
    double minY = Y.get(0), maxY = Y.get(0);
    
    for (int i = 0; i < X.size(); i++) {
      double x = X.get(i);
      double y = Y.get(i);
      
      if (x < minX) {
        minX = x;
      } else if (x > maxX) {
        maxX = x;
      }
      
      if (y < minY) {
        minY = y;
      } else if (y > maxY) {
        maxY = y;
      }
    }
    
    double xDif = maxX - minX;
    double yDif = maxY - minY;
    
    double xRatio = xDif / width;
    double yRatio = yDif / height;
    double scaleFactor;
    
    if (xRatio > yRatio) {
      scaleFactor = .85 / xRatio; // Change .9 to reflect start buffer?
    } else {
      scaleFactor = .85 / yRatio;
    }
    
    double xShift = (width / 2.0) - (scaleFactor * ((xDif/2.0) + minX));
    double yShift = (height / 2.0) - (scaleFactor * ((yDif/2.0) + minY));
    
    for (int i = 0; i < X.size(); i++) {
      X.set(i, (X.get(i) * scaleFactor) + xShift);
      Y.set(i, (Y.get(i) * scaleFactor) + yShift);
    }
    
    
    // Calculate starting backbone length
    double avgDist = 0;
    int numDists = 0;
    
    int numCounted = int(X.size() / 4);
    if (numCounted < 50) { numCounted = 50; }
    if (numCounted > X.size()) { numCounted = X.size(); }
    
    for (int i = 1; i < numCounted; i++) {
      if (i >= X.size()) { break; }
      numDists++;
      avgDist += normFunction(X.get(i) - X.get(i-1), Y.get(i) - Y.get(i-1), 0);
    }
    
    if (numDists != 0) {
      avgDist /= numDists;
    } else {
      avgDist = ff.minBackboneDist;
    }
    
    if (avgDist < 8) {
      ff.backboneDist = avgDist * .35;
    } else {
      ff.backboneDist = avgDist * .95;
    }
    
    if (ff.backboneDist < ff.minBackboneDist) { ff.backboneDist = ff.minBackboneDist; }
    if (ff.backboneDist > (ff.minBackboneDist + ff.backboneDistScale)) { ff.backboneDist = ff.minBackboneDist + ff.backboneDistScale; }
    
    ff.radius = avgDist / 3;
    if (ff.radius < ff.minRadius) { ff.radius = ff.minRadius; }
    if (ff.radius > ff.minRadius + ff.radiusScale) { ff.radius = ff.minRadius + ff.radiusScale; }
    
    
    // Calculate starting bond length
    avgDist = 0;
    numDists = 0;
    for (int i = 0; i < pairTableOptimal.size(); i++) {
      int pairedWId = pairTableOptimal.get(i);
      if (pairedWId != -1 && i < pairedWId) {
        numDists++;
        avgDist += normFunction( X.get(i) - X.get(pairedWId), Y.get(i) - Y.get(pairedWId), 0 );
      }
    }
    
    if (numDists != 0) {
      avgDist /= numDists;
    } else {
      avgDist = ff.minBpDist;
    }
    
    ff.bpDist = avgDist;
    if (ff.bpDist < ff.minBpDist) { ff.bpDist = ff.minBpDist; }
    if (ff.bpDist > (ff.minBpDist + ff.bpDistScale)) { ff.bpDist = ff.minBpDist + ff.bpDistScale; }
    
    ds.padding = ds.startPadding;
    spheres.displayOrder.clear();
    spheres.relaxedIds.clear();
    for (int i = 0; i < lenTot; ++i) {
      int pairTableId = residueOrder[i];
      spheres.get(i).x = X.get(pairTableId);
      spheres.get(i).y = Y.get(pairTableId);
      spheres.displayOrder.add(i);
    }
    
    simulationMode = false;
    //movieMode = false;
    rigidHelices = true;
    
    textDisplay = "RADIAL LAYOUT";
    
    simState.spheres = spheres;
    
    if (debugPrints) { println("Finished SimulationSetup.createRadialState!"); }
    return simState;
  }
}