content: /*
   * The MIT License (MIT)
   *
   * Copyright (c) 2011-2015 Broad Institute, Aiden Lab
   *
   * Permission is hereby granted, free of charge, to any person obtaining a copy
   * of this software and associated documentation files (the "Software"), to deal
   * in the Software without restriction, including without limitation the rights
   * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   * copies of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included in
   * all copies or substantial portions of the Software.
   *
   *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   *  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   *  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   *  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   *  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   *  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
   *  THE SOFTWARE.
   *
   * Based on PicardTools 
   * Juicer version 1.5
   */
  import java.text.NumberFormat;
  import java.text.ParseException;
  import java.util.Locale;
  import java.io.FileReader;
  import java.io.BufferedReader;
  import java.io.IOException;
  import java.io.File;
  class LibraryComplexity {
    public static void main(String[] args) {
      if (args.length != 2 && args.length != 3 && args.length != 1) {
        System.out.println("Usage: java LibraryComplexity <directory> <output file>");
        System.out.println("     : java LibraryComplexity <unique> <pcr> <opt>");
        System.exit(0);
      }
      NumberFormat nf = NumberFormat.getInstance(Locale.US);
  
      BufferedReader reader = null;
      long readPairs = 0;
      long uniqueReadPairs = 0;
      long opticalDups = 0;
      long totalReadPairs = 0;
      if (args.length == 2 || args.length==1) {
        try {
          File f = new File(args[0] + "/opt_dups.txt");
          if (f.exists()) {
            reader = new BufferedReader(new FileReader(f));
            while (reader.readLine() != null) opticalDups++;
            reader.close();
          }
          f = new File(args[0] + "/merged_nodups.txt");
          if (f.exists()) {
            reader = new BufferedReader(new FileReader(f));
            while (reader.readLine() != null) uniqueReadPairs++;
            reader.close();				
          }
          f = new File(args[0] + "/dups.txt");
          if (f.exists()) {
            reader = new BufferedReader(new FileReader(args[0] + "/dups.txt"));
            while (reader.readLine() != null) readPairs++;
            reader.close();				
            readPairs += uniqueReadPairs; 
          }
          String fname = "inter.txt";
          if (args.length == 2) fname = args[1];
          f = new File(args[0] + "/" + fname);
          if (f.exists()) {
            reader = new BufferedReader(new FileReader(args[0] + "/" + fname));
            String line = reader.readLine();
            boolean done = false;
            while (line != null && !done) {
              if (line.contains("Sequenced Read")) {
                String[] parts = line.split(":");
                try {
                  totalReadPairs = nf.parse(parts[1].trim()).longValue();
                }
                catch (ParseException e) {
                  totalReadPairs = 0;
                }
                done = true;
              }
              line = reader.readLine(); 
            }
            reader.close();
          }
        } catch (IOException error) {
          System.err.println("Problem counting lines in merged_nodups and dups");
          System.exit(1);
        }
      }
      else {
        try {
          uniqueReadPairs = Integer.valueOf(args[0]);
          readPairs = Integer.valueOf(args[1]);
          opticalDups = Integer.valueOf(args[2]);
        }
        catch (NumberFormatException error) {
          System.err.println("When called with three arguments, must be integers");
          System.exit(1);
        }
        readPairs += uniqueReadPairs;
      }
      NumberFormat decimalFormat = NumberFormat.getPercentInstance();
      decimalFormat.setMinimumFractionDigits(2);
      decimalFormat.setMaximumFractionDigits(2);
  
      System.out.print("Unique Reads: " + NumberFormat.getInstance().format(uniqueReadPairs) + " ");
      if (totalReadPairs > 0) {
        System.out.println("(" + decimalFormat.format(uniqueReadPairs/(double)totalReadPairs) + ")");
      }
      else {
        System.out.println();
      }
  
      System.out.print("PCR Duplicates: " + nf.format(readPairs - uniqueReadPairs) + " ");
      if (totalReadPairs > 0) {
        System.out.println("(" + decimalFormat.format((readPairs - uniqueReadPairs)/(double)totalReadPairs) + ")");
      }
      else {
        System.out.println();
      }
      System.out.print("Optical Duplicates: " + nf.format(opticalDups) + " ");
      if (totalReadPairs > 0) {
        System.out.println("(" + decimalFormat.format(opticalDups/(double)totalReadPairs) + ")");
      }
      else {
        System.out.println();
      }
      long result;
      try {
        result = estimateLibrarySize(readPairs, uniqueReadPairs);
      }
      catch (NullPointerException e) {
        result = 0;
      }
      System.out.println("Library Complexity Estimate: " + nf.format(result));
    }
  
    /**
     * Estimates the size of a library based on the number of paired end molecules 
     * observed and the number of unique pairs observed.
     * <br>
     * Based on the Lander-Waterman equation that states:<br>
     *     C/X = 1 - exp( -N/X )<br>
     * where<br>
     *     X = number of distinct molecules in library<br>
     *     N = number of read pairs<br>
     *     C = number of distinct fragments observed in read pairs<br>
     */
    public static Long estimateLibrarySize(final long readPairs, 
                                           final long uniqueReadPairs) {
      final long readPairDuplicates = readPairs - uniqueReadPairs;
      
      if (readPairs > 0 && readPairDuplicates > 0) {
  	    long n = readPairs;
  	    long c = uniqueReadPairs;
  	    
  	    double m = 1.0, M = 100.0;
  	    
  	    if (c >= n || f(m*c, c, n) < 0) {
          throw new IllegalStateException("Invalid values for pairs and unique pairs: " + n + ", " + c);
  	    }
        
  	    while( f(M*c, c, n) >= 0 ) {
          m = M;
          M *= 10.0;
  	    }
        
  	    double r = (m+M)/2.0;
  	    double u = f( r * c, c, n );
  	    int i = 0;
  	    while (u != 0 && i < 1000) {
          if (u > 0) m = r;
          else M = r;
          r = (m+M)/2.0;
          u = f( r * c, c, n );
          i++;
  	    }
  	    if (i == 1000) {
          System.err.println("Iterated 1000 times, returning estimate");
  	    }
  	    
  	    return (long) (c * (m+M)/2.0);
      }
      else {
  	    return null;
      }
    }
    
    /** Method that is used in the computation of estimated library size. */
    private static double f(double x, double c, double n) {
      return c/x - 1 + Math.exp(-n/x);
    }
  }    
  