content: /*************************************************************************
  *                                                                        *
  *  This file is part of the 20n/act project.                             *
  *  20n/act enables DNA prediction for synthetic biology/bioengineering.  *
  *  Copyright (C) 2017 20n Labs, Inc.                                     *
  *                                                                        *
  *  Please direct all queries to act@20n.com.                             *
  *                                                                        *
  *  This program is free software: you can redistribute it and/or modify  *
  *  it under the terms of the GNU General Public License as published by  *
  *  the Free Software Foundation, either version 3 of the License, or     *
  *  (at your option) any later version.                                   *
  *                                                                        *
  *  This program is distributed in the hope that it will be useful,       *
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  *  GNU General Public License for more details.                          *
  *                                                                        *
  *  You should have received a copy of the GNU General Public License     *
  *  along with this program.  If not, see <http://www.gnu.org/licenses/>. *
  *                                                                        *
  *************************************************************************/
  
  package com.act.lcms;
  
  import java.util.List;
  import java.util.Arrays;
  import java.util.ArrayList;
  import java.util.Iterator;
  import java.util.Collections;
  import java.io.PrintStream;
  import java.io.OutputStream;
  import java.io.FileOutputStream;
  import org.apache.commons.lang3.tuple.Pair;
  
  public class CompareTwoNetCDFAroundMass {
  
    static final double mzTolerance = 0.01;
    static final int maxNumMzTolerated = 3;
  
    private double extractMZ(double mzWanted, List<Pair<Double, Double>> intensities, double time) {
      double intensityFound = 0;
      int numWithinPrecision = 0;
      double mzLowRange = mzWanted - mzTolerance;
      double mzHighRange = mzWanted + mzTolerance;
      // we expect there to be pretty much only one intensity value in the precision
      // range we are looking at. But if a lot of masses show up then complain
      for (Pair<Double, Double> mz_int : intensities) {
        double mz = mz_int.getLeft();
        double intensity = mz_int.getRight();
  
  
        if (mz > mzLowRange && mz < mzHighRange) {
          intensityFound += intensity;
          numWithinPrecision++;
          if (intensity > 100000)
            System.out.format("time: %f\tmz: %f\t intensity: %f\n", time, mz, intensity);
        }
      }
  
      if (numWithinPrecision > maxNumMzTolerated) {
        System.out.format("Only expected %d, but found %d in the mz range [%f, %f]\n", maxNumMzTolerated, 
            numWithinPrecision, mzLowRange, mzHighRange);
      }
  
      return intensityFound;
    }
  
    private List<Pair<Double, Double>> extractMZSlice(double mz, Iterator<LCMSSpectrum> spectraIt, int numOut, String tagfname) {
      List<Pair<Double, Double>> mzSlice = new ArrayList<>();
  
      int pulled = 0;
      // iterate through first few, or all if -1 given
      while (spectraIt.hasNext() && (numOut == -1 || pulled < numOut)) {
        LCMSSpectrum timepoint = spectraIt.next();
  
        // get all (mz, intensity) at this timepoint
        List<Pair<Double, Double>> intensities = timepoint.getIntensities();
  
        // this time point is valid to look at if its max intensity is around
        // the mass we care about. So lets first get the max peak location
        double intensityForMz = extractMZ(mz, intensities, timepoint.getTimeVal());
  
        // the above is Pair(mz_extracted, intensity), where mz_extracted = mz
        // we now add the timepoint val and the intensity to the output
        mzSlice.add(Pair.of(timepoint.getTimeVal(), intensityForMz));
  
        pulled++;
      }
      System.out.format("Done: %s\n\n", tagfname);
  
      return mzSlice;
    }
  
    /**
     * Extracts a window around a particular m/z value within the spectra
     * @param file The netCDF file with the full LCMS (lockmass corrected) spectra
     * @param mz The m/z value +- 1 around which to extract data for
     * @param numTimepointsToExamine Since we stream the netCDF in, we can choose
     *                               to look at the first few timepoints or -1 for all
     */
    public List<List<Pair<Double, Double>>> getSpectraForMass(double mz, String[] fnames, int numTimepointsToExamine) 
      throws Exception {
      List<List<Pair<Double, Double>>> extracted = new ArrayList<>();
      LCMSParser parser = new LCMSNetCDFParser();
  
      for (String fname : fnames) {
        Iterator<LCMSSpectrum> iter = parser.getIterator(fname);
        List<Pair<Double, Double>> mzSlice = extractMZSlice(mz, iter, numTimepointsToExamine, fname);
        extracted.add(mzSlice);
      }
  
      return extracted;
    }
  
    private static boolean areNCFiles(String[] fnames) {
      for (String n : fnames) {
        System.out.println(".nc file = " + n);
        if (!n.endsWith(".nc"))
          return false;
      }
      return true;
    }
  
    public static void main(String[] args) throws Exception {
      if (args.length < 5 || !areNCFiles(Arrays.copyOfRange(args, 3, args.length))) {
        throw new RuntimeException("Needs: \n" + 
            "(1) mass value, e.g., 132.0772 for debugging, \n" +
            "(2) how many timepoints to process (-1 for all), \n" +
            "(3) prefix for .data and rendered .pdf \n" +
            "(4,5..) 2 or more NetCDF .nc files"
            );
      }
  
      String fmt = "pdf";
      Double mz = Double.parseDouble(args[0]);
      Integer numSpectraToProcess = Integer.parseInt(args[1]);
      String outPrefix = args[2];
      String outImg = outPrefix.equals("-") ? null : outPrefix + "." + fmt;
      String outData = outPrefix.equals("-") ? null : outPrefix + ".data";
  
      CompareTwoNetCDFAroundMass c = new CompareTwoNetCDFAroundMass();
      String[] netCDF_fnames = Arrays.copyOfRange(args, 3, args.length);
      List<List<Pair<Double, Double>>> spectra = c.getSpectraForMass(mz, netCDF_fnames, numSpectraToProcess);
  
      // Write data output to outfile
      PrintStream out = outData == null ? System.out : new PrintStream(new FileOutputStream(outData));
  
      // print out the spectra to outData
      for (List<Pair<Double, Double>> spectraInFile : spectra) {
        for (Pair<Double, Double> xy : spectraInFile) {
          out.format("%.4f\t%.4f\n", xy.getLeft(), xy.getRight());
          out.flush();
        }
        // delimit this dataset from the rest
        out.print("\n\n");
      }
      // find the ymax across all spectra, so that we can have a uniform y scale
      Double yrange = 0.0;
      for (List<Pair<Double, Double>> spectraInFile : spectra) {
        Double ymax = 0.0;
        for (Pair<Double, Double> xy : spectraInFile) {
          Double intensity = xy.getRight();
          if (ymax < intensity) ymax = intensity;
        }
        if (yrange < ymax) yrange = ymax;
      }
  
      if (outData != null) {
        // if outData is != null, then we have written to .data file
        // now render the .data to the corresponding .pdf file
  
        // first close the .data
        out.close();
  
        // render outData to outFILE using gnuplo
        Gnuplotter plotter = new Gnuplotter();
        plotter.plot2D(outData, outImg, netCDF_fnames, "time in seconds", yrange, "intensity", fmt);
      }
    }
  }
  