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.Map;
  import java.util.HashMap;
  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;
  import java.io.File;
  
  public class AnimateNetCDFAroundMass {
    class XYZ {
      Double time;
      Double mz;
      Double intensity;
      public XYZ(Double time, Double mz, Double intensity) {
        this.time = time;
        this.mz = mz;
        this.intensity = intensity;
      }
    }
  
    private List<List<XYZ>> getSpectraInWindowAll(List<List<XYZ>> spectra, Double time, Double tWin, 
        Double mz, Double mzWin) {
      List<List<XYZ>> windowOnMultiple = new ArrayList<>();
      for (List<XYZ> one : spectra)
        windowOnMultiple.add(getSpectraInWindow(one, time, tWin, mz, mzWin));
      return windowOnMultiple;
    }
  
    final static int gridSize = 100;
  
    private List<XYZ> getSpectraInWindow(List<XYZ> spectra, Double time, Double tWin, Double mz, Double mzWin) {
  
      // tWin and mzWin are half windows around time and mz, respectively
      // Range of mz is ( mz - mzWin, mz + mzWin )
      // Range of time is ( time - timeWin, time + timeWin )
  
      double mzLow = mz - mzWin;
      double mzHigh = mz + mzWin;
      double timeLow = time - tWin;
      double timeHigh = time + tWin;
  
      Map<Pair<Long, Long>, List<Double>> gridVals = new HashMap<>();
      // initialize the grid to empty intensity data points 
      for (long i = 0; i <= gridSize; i++) {
        for (long j = 0; j <= gridSize; j++) {
          gridVals.put(Pair.of(i, j), new ArrayList<>());
        }
      }
  
      // what is the step size between x -> x+1 and y -> y+1
      double timeStep = (2 * tWin) / gridSize;
      double mzStep = (2 * mzWin) / gridSize;
  
      // for each x,y,z find the grid position it lies in
      // (if it indeed is within the grid window we care)
      // and add that to the data points we see within
      for (XYZ xyz : spectra) {
        // ignore if outside the window
        if (xyz.mz < mzLow || xyz.mz > mzHigh)
          continue;
        if (xyz.time < timeLow || xyz.time > timeHigh)
          continue;
  
        // if inside the window find the grid position
        long timeGridPos = Math.round((xyz.time - timeLow) / timeStep);
        long mzGridPos = Math.round((xyz.mz - mzLow) / mzStep);
        Pair<Long, Long> xy = Pair.of(timeGridPos, mzGridPos);
  
        gridVals.get(xy).add(xyz.intensity);
      }
  
      // average out all the z values that appear at each x,y
      Map<Pair<Long, Long>, Double> gridAvg = new HashMap<>();
      for (Pair<Long, Long> xy : gridVals.keySet()) {
        Double avg = 0.0;
        List<Double> intensities = gridVals.get(xy);
        for (Double intensity : intensities)
          avg += intensity;
        if (!intensities.isEmpty()) 
          avg /= intensities.size();
        gridAvg.put(xy, avg);
      }
  
      // convert back to list of x,y,z coordinates
      List<XYZ> grid = new ArrayList<>();
      for (long i = 0; i <= gridSize; i++) {
        for (long j = 0; j <= gridSize; j++) {
          Double t = timeLow + i * timeStep;
          Double m = mzLow + j * mzStep;
          Double it = gridAvg.get(Pair.of(i, j));
          grid.add(new XYZ(t, m, it));
        }
      }
  
      return grid;
    }
  
    private List<XYZ> getSpectra(Iterator<LCMSSpectrum> spectraIt, Double time, Double tWin, Double mz, 
        Double mzWin) {
      double mzLow = mz - mzWin;
      double mzHigh = mz + mzWin;
      double timeLow = time - tWin;
      double timeHigh = time + tWin;
  
      List<XYZ> spectra = new ArrayList<>();
  
      while (spectraIt.hasNext()) {
        LCMSSpectrum timepoint = spectraIt.next();
  
        if (timepoint.getTimeVal() < timeLow || timepoint.getTimeVal() > timeHigh)
          continue;
  
        // get all (mz, intensity) at this timepoint
        for (Pair<Double, Double> mz_int : timepoint.getIntensities()) {
          double mzHere = mz_int.getLeft();
          double intensity = mz_int.getRight();
  
          if (mzHere < mzLow || mzHere > mzHigh)
            continue;
  
          spectra.add(new XYZ(timepoint.getTimeVal(), mzHere, intensity));
        }
      }
  
      return spectra;
    }
  
    public List<List<XYZ>> getSpectra(String[] fnames, Double time, Double tWin, Double mz, 
        Double mzWin) throws Exception {
      List<List<XYZ>> extracted = new ArrayList<>();
      LCMSParser parser = new LCMSNetCDFParser();
  
      for (String fname : fnames) {
        Iterator<LCMSSpectrum> iter = parser.getIterator(fname);
        extracted.add(getSpectra(iter, time, tWin, mz, mzWin));
      }
  
      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 < 7 || !areNCFiles(Arrays.copyOfRange(args, 5, args.length))) {
        throw new RuntimeException("Needs: \n" + 
            "(1) mass value, e.g., 132.0772 \n" +
            "(2) time value, e.g., 39.2, (seconds), \n" +
            "(3) minimum Mz Precision, 0.04 \n" +
            "(4) max z axis, e.g., 20000 \n" +
            "(5) prefix for .data and rendered .pdf \n" +
            "(6..) 2 or more NetCDF .nc files"
            );
      }
  
      Double mz = Double.parseDouble(args[0]);
      Double time = Double.parseDouble(args[1]);
      Double minMzPrecision = Double.parseDouble(args[2]);
      Double maxZAxis = Double.parseDouble(args[3]);
      String outPrefix = args[4];
  
      // the mz values go from 50-950, we start with a big window and exponentially narrow down
      double mzWin = 100;
      // time values go from 0-450, we start with a big window and exponentially narrow down
      double timeWin = 50;
  
      // the factor by which to zoom in every step (has to be >1, a value of 2 is good)
      double factor = 1.2;
  
      // the animation frame count
      int frame = 1;
  
      AnimateNetCDFAroundMass c = new AnimateNetCDFAroundMass();
      String[] netCDFFnames = Arrays.copyOfRange(args, 5, args.length);
      List<List<XYZ>> spectra = c.getSpectra(netCDFFnames, time, timeWin, mz, mzWin);
  
      for (List<XYZ> s : spectra) {
        System.out.format("%d xyz datapoints in (initial narrowed) spectra\n", s.size());
      }
  
      String[] labels = new String[netCDFFnames.length];
      for (int i=0; i<labels.length; i++)
        labels[i] = "Dataset: " + i;
      // you could set labels to netCDFFnames to get precise labels on the graphs
  
      Gnuplotter plotter = new Gnuplotter();
      String fmt = "png";
  
      List<String> outImgFiles = new ArrayList<>(), outDataFiles = new ArrayList<>();
      while (mzWin > minMzPrecision) {
  
        // exponentially narrow windows down
        mzWin /= factor;
        timeWin /= factor;
  
        List<List<XYZ>> windowedSpectra = c.getSpectraInWindowAll(spectra, time, timeWin, mz, mzWin);
  
        String frameid = String.format("%03d", frame);
        String outPDF = outPrefix + frameid + "." + fmt;
        String outDATA = outPrefix + frameid + ".data";
        outImgFiles.add(outPDF); 
        outDataFiles.add(outDATA);
        frame++;
  
        // Write data output to outfile
        PrintStream out = new PrintStream(new FileOutputStream(outDATA));
  
        // print out the spectra to outDATA
        for (List<XYZ> windowOfSpectra : windowedSpectra) {
          for (XYZ xyz : windowOfSpectra) {
            out.format("%.4f\t%.4f\t%.4f\n", xyz.time, xyz.mz, xyz.intensity);
            out.flush();
          }
          // delimit this dataset from the rest
          out.print("\n\n");
        }
  
        // close the .data
        out.close();
  
        // render outDATA to outPDF using gnuplot
        plotter.plotMulti3D(outDATA, outPDF, fmt, labels, maxZAxis);
      }
  
      String outImgs = outPrefix + "*." + fmt;
      plotter.makeAnimatedGIF(outImgs, outPrefix + ".gif");
      // all the frames are now in the animated gif, remove the intermediate files
      for (String f: outDataFiles) 
        new File(f).delete();
      for (String f: outImgFiles) 
        new File(f).delete();
    }
  }
  