sequencialDeconvolution.cpp

// The "more usual approach":
// - smoothing
// - SNIP (baseline removal)
// - l1-deconvolution
//
#include <csv.h>
#include <cxxopts.hpp>

#include <fstream>
#include <sstream>

#include "seqDeconv.hpp"
#include "sgFilter.hpp"
#include "snip.hpp"

namespace JointDeconv
{
};  // TODO: to remove

using namespace JointDeconv;

int main(int argc, char* argv[])
{
  //////////////////////
  // Argument parsing //
  //////////////////////

  const auto to_string = [](const auto& any) {
    std::stringstream converter;
    converter << any;
    return std::string(converter.str());
  };

  std::string command_line;
  for (int i = 0; i < argc; i++)
  {
    command_line += argv[i];
    command_line += ' ';
  }
  cxxopts::Options options(argv[0],
                           "A more classical approach "
                           "smoothing->SNIP->deconvolution (for comparison "
                           "purpose only), contact vincent.picaud@cea.fr");

  // Default parameter values
  //
  double sigma_left = 10;
  double sigma_right = 10;
  Size_t snip_halfWindowSize = 20;
  double peakMinHeight = 0.01;
  double lambda_1 = 0.1;
  double lambda_2 = 0.00001;
  double mu = 500;
  double eps = 1e-4;
  Size_t max_iter = 5000;

  bool gnuplot = false;

  options.add_options()
      //
      ("i,input",
       "Input file (two columns X,Y)",
       cxxopts::value<std::string>(),
       "FILE")
      //
      ("o,output",
       "Output file",
       cxxopts::value<std::string>()->default_value("$(FILE).out"),
       "OUTPUT FILE")
      //
      ("snip",
       "Snip half window size (>0), used to compute baseline",
       cxxopts::value<Size_t>(snip_halfWindowSize),
       to_string(snip_halfWindowSize))
      //
      ("sigma_left",
       "Peak shape factor (>0)",
       cxxopts::value<double>(sigma_left),
       to_string(sigma_left))
      //
      ("sigma_right",
       "Peak shape factor (>0)",
       cxxopts::value<double>(sigma_right),
       to_string(sigma_right))
      //
      ("peakMinHeight",
       "Minimal height to accept peak (>=0)",
       cxxopts::value<double>(peakMinHeight),
       to_string(peakMinHeight))
      //
      ("lambda_1",
       "lambda_1 penalty term (>=0)",
       cxxopts::value<double>(lambda_1),
       to_string(lambda_1))
      //
      ("lambda_2",
       "lambda_2 penalty term (>=0)",
       cxxopts::value<double>(lambda_2),
       to_string(lambda_2))
      //
      ("mu", "mu penalty term (>0)", cxxopts::value<double>(mu), to_string(mu))
      //
      ("eps", "eps goal (>=0)", cxxopts::value<double>(eps), to_string(eps))
      //
      ("max_iter",
       "maximum number of iterations (>0)",
       cxxopts::value<Size_t>(max_iter),
       to_string(max_iter))
      //
      ("p,gnuplot", "Gnuplot script", cxxopts::value<bool>(gnuplot))
      //
      ("help", "Print help");

  options.parse_positional("input");

  // CAVEAT: after options.parse(argc, argv); argc is set=1
  // -> we have to check its value before
  //
  bool show_help = argc == 1;

  options.parse(argc, argv);

  show_help |= options.count("help") > 0;

  if (show_help)
  {
    std::cout << options.help() << std::endl;
    exit(0);
  }

  // Sanity check
  //
  if (options.count("input") == 0)
  {
    std::cout << "#Error: missing input FILE" << std::endl;
    return EXIT_FAILURE;
  }
  if (sigma_left <= 0)
  {
    std::cout << "#Error: sigma_left= " << sigma_left
              << " is not a positive number" << std::endl;
    return EXIT_FAILURE;
  }
  if (sigma_right <= 0)
  {
    std::cout << "#Error: sigma_right= " << sigma_right
              << " is not a positive number" << std::endl;
    return EXIT_FAILURE;
  }
  if (peakMinHeight < 0)
  {
    std::cout << "#Error: peakMinHeight= " << peakMinHeight
              << " is not a nonegative number" << std::endl;
    return EXIT_FAILURE;
  }
  if (lambda_1 < 0)
  {
    std::cout << "#Error: lambda_1= " << lambda_1
              << " is not a nonnegative number" << std::endl;
    return EXIT_FAILURE;
  }
  if (lambda_2 < 0)
  {
    std::cout << "#Error: lambda_2= " << lambda_2
              << " is not a nonnegative number" << std::endl;
    return EXIT_FAILURE;
  }
  if (mu <= 0)
  {
    std::cout << "#Error: mu= " << mu << " is not a positive number"
              << std::endl;
    return EXIT_FAILURE;
  }
  if (eps < 0)
  {
    std::cout << "#Error: eps= " << eps << " is not a nonnegative number"
              << std::endl;
    return EXIT_FAILURE;
  }
  if (max_iter <= 0)
  {
    std::cout << "#Error: max_iter= " << max_iter << " is not a positive number"
              << std::endl;
    return EXIT_FAILURE;
  }

  // Generates output filename
  //
  const std::string output_filename((options.count("output") == 0)
                                        ? options["input"].as<std::string>() +
                                              ".out"
                                        : options["output"].as<std::string>());

  ///////////////////
  // Read CSV file //
  ///////////////////

  std::vector<std::pair<double, double>> data;
  try
  {
    io::CSVReader<2,
                  io::trim_chars<' ', '\t'>,
                  io::no_quote_escape<','>,
                  io::throw_on_overflow,
                  io::single_line_comment<'#'>>
        in(options["input"].as<std::string>());
    double X_i, Y_i;

    while (in.read_row(X_i, Y_i))
    {
      data.push_back(std::make_pair(X_i, Y_i));
    }
  }
  catch (std::exception& e)
  {
    std::cerr << "#Error: " << e.what() << std::endl;
    return EXIT_FAILURE;
  }

  ///////////////////////////////////////
  // Call the deconvolution subroutine //
  ///////////////////////////////////////

  const Size_t n = data.size();

  if (n <= 2)
  {
    std::cerr << "#Error: spectrum must have at least 2 points" << std::endl;
    return EXIT_FAILURE;
  }

  Vector x(n);
  Vector y(n);
  Vector y_smoothed(n);
  Vector deconvolvedPeak(n);
  Vector convolvedPeak(n);
  Vector baseline(n);

  for (Index_t i = 0; i < n; ++i)
  {
    x[i] = data[i].first;
    y[i] = data[i].second;
  }

  data.clear();

  SeqDeconv_InputParameters inputParameters;
  inputParameters.lambda_1 = lambda_1;
  inputParameters.lambda_2 = lambda_2;
  inputParameters.solver_inputParameters.iter_max = max_iter;
  inputParameters.solver_inputParameters.eps_goal = eps;

  // Smooth
  sgFilter(y, y_smoothed);
  // SNIP
  snip(y_smoothed, baseline, snip_halfWindowSize);
  // Deconv
  y -= baseline;
  seqDeconv_GaussianPeaks(x,
                          y,
                          peakMinHeight,
                          sigma_left,
                          sigma_right,
                          deconvolvedPeak,
                          convolvedPeak,
                          inputParameters);
  y += baseline;

  //////////////////
  // Write output //
  //////////////////

  std::ofstream output(output_filename);
  output << "# Generated by: " << command_line << std::endl;
  for (Index_t i = 0; i < n; ++i)
  {
    output << x[i] << ", " << y[i] << ", " << deconvolvedPeak[i] << ", "
           << baseline[i] << ", " << convolvedPeak[i] << ", "
           << baseline[i] + convolvedPeak[i] << std::endl;
  }

  output.close();

  // Gnuplot
  // load "test.csv.out.gnuplot"
  if (gnuplot)
  {
    std::ofstream output_gnuplot(output_filename + ".gnuplot");

    output_gnuplot << "set datafile separator ','" << std::endl;
    output_gnuplot << "plot \"" << output_filename
                   << "\" u 1 : 2 w l t \"Y raw\" lc \"blue\"" << std::endl;
    output_gnuplot << "replot \"" << output_filename
                   << "\" u 1 : 6 w l t \"reconstructed\" lw 2 lc \"grey30\""
                   << std::endl;
    output_gnuplot << "replot \"" << output_filename
                   << "\" u 1 : 4 w l t \"baseline\" lw 2 lc \"grey70\""
                   << std::endl;
    output_gnuplot << "replot \"" << output_filename
                   << "\" u 1 : 5 w l t \"convolved peaks\" lc \"green\""
                   << std::endl;
    output_gnuplot << "replot \"" << output_filename
                   << "\" u 1 : 3 w i t \"peaks\" lw 2 lc \"red\"" << std::endl;
    output_gnuplot << "set terminal eps" << std::endl;
    output_gnuplot << "set output \"" << output_filename + ".eps"
                   << "\"" << std::endl;
    output_gnuplot << "replot" << std::endl;
    output_gnuplot << "set terminal png" << std::endl;
    output_gnuplot << "set output \"" << output_filename + ".png"
                   << "\"" << std::endl;
    output_gnuplot << "replot" << std::endl;
    output_gnuplot << "set terminal qt" << std::endl;
    output_gnuplot << "set output" << std::endl;
  }

  return EXIT_SUCCESS;
}