// File: test_timsxicextractor.cpp
// Created by: Olivier Langella
// Created on: 3/2/2021
//
/*******************************************************************************
 * Copyright (c) 2021 Olivier Langella
 *<olivier.langella@universite-paris-saclay.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     PAPPSOms++ 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.
 *
 *     PAPPSOms++ 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 PAPPSOms++.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

// Run using this command line:
// tests/catch2-only-tests "Creating mass spectrum via TIMS XIC extraction from
// timsdata" -s

#include <catch2/catch_test_macros.hpp>
#include <catch2/benchmark/catch_benchmark.hpp>


#include <QDebug>

#include <pappsomspp/core/mzrange.h>
#include <pappsomspp/core/msfile/msfileaccessor.h>
#include <pappsomspp/core/msrun/private/timsmsrunreader.h>
#include <pappsomspp/core/msrun/private/timsmsrunreaderms2.h>
#include <pappsomspp/core/pappsoexception.h>
#include <pappsomspp/core/processing/filters/filtersuitestring.h>
#include <pappsomspp/core/xicextractor/msrunxicextractorfactory.h>
#include <pappsomspp/core/processing/uimonitor/uimonitortext.h>
#include <pappsomspp/core/vendors/tims/timsddaprecursors.h>
#include "tests/tests-config.h"
#include "../common.h"
#include <time.h>


using namespace pappso;

TEST_CASE("Creating mass spectrum via TIMS XIC extraction from timsdata",
          "[massSpecTimsXicExtractor]")
{
  qSetMessagePattern(QString("%{file}@%{line}, %{function}(): %{message}"));

#if USE_PAPPSO_TREE == 1

  QString file_path_name;

  // file_path_name =
  //"/home/rusconi/devel/dataForMzml/bruker/20210126_HeLa.d/"
  //"1-26-2021_1_QC_HeLa10ng_826.d/analysis.tdf";

  file_path_name =
    "/gorgone/pappso/jouy/raw/2021_Tims_TOF/20211124_HeLa/"
    "11-25-2021_1_HeLa200ng_2321.d/analysis.tdf";

  // file_path_name =
  //"/home/rusconi/devel/dataForMzml/bruker/2021_Tims_TOF/20211124_HeLa/"
  //"11-25-2021_1_HeLa200ng_2321.d/analysis.tdf";

  qDebug() << "The file to use as a test base is: " << file_path_name;
  // When not in debug mode.
  std::cout << __FILE__ << ":" << __LINE__
            << " The file to use as a test base is: "
            << file_path_name.toStdString() << std::endl;

  pappso::MsFileAccessor accessor(file_path_name, "a1");

  qDebug() << "Setting the preferred file reader type to tims2.";

  accessor.setPreferredFileReaderType(Enums::MsDataFormat::brukerTims,
                                      Enums::FileReaderType::tims_ms2);

  pappso::MsRunReaderSPtr p_msreader =
    accessor.msRunReaderSPtr(accessor.getMsRunIds().front());

  REQUIRE(p_msreader != nullptr);

  REQUIRE(accessor.getFileReaderType() == Enums::FileReaderType::tims_ms2);

  pappso::TimsMsRunReaderMs2 *tims_reader =
    dynamic_cast<pappso::TimsMsRunReaderMs2 *>(p_msreader.get());
  REQUIRE(tims_reader != nullptr);

  pappso::TimsDataSp tims_data = tims_reader->getTimsDataSPtr();

  // std::vector<std::size_t> precursor_list = {132929};
  std::vector<std::size_t> precursor_list = {31693};

  // Now craft the mz list. We have an acquisition m/z range from 100 to 1700.
  std::vector<pappso_double> mz_list;

  // Seed the lower m/z value of the m/z range and push it back to the mz_list.
  double mz = 100;
  mz_list.push_back(mz);

  // We'll need a precision.
  PrecisionPtr precision_p = PrecisionFactory::getPpmInstance(30);

  // Compute the last mz value accounting for the delta (we do not want to
  // over-trip the 1700 hard limit.
  double last_mz = 1700;
  last_mz -= 2 * precision_p->delta(last_mz);

  while(mz < last_mz)
    {
      mz += 2 * precision_p->delta(mz);
      mz_list.push_back(mz);
    }

  qDebug() << "The mz_list now has " << mz_list.size() << " m/z values.";
  std::cout << "The mz_list now has " << mz_list.size() << " m/z values.";

  std::vector<pappso::XicCoordSPtr> xic_list;
  std::vector<pappso::XicCoordSPtr> xic_coord_list;

  clock_t start = clock();

  for(auto precursor_id : precursor_list)
    {
      XicCoordSPtr xic_coord_tims_struct_sp =
        tims_data.get()
          ->getTimsDdaPrecursorsPtr()
          ->getXicCoordTimsFromPrecursorId(precursor_id,
                                           PrecisionFactory::getPpmInstance(30))
          .initializeAndClone();

      xic_list.push_back(xic_coord_tims_struct_sp);
    }


  for(pappso::pappso_double mz : mz_list)
    {
      pappso::XicCoordSPtr new_xic_coord =
        xic_list[0].get()->initializeAndClone();
      new_xic_coord.get()->mzRange =
        pappso::MzRange(mz, PrecisionFactory::getPpmInstance(30.0));
      xic_coord_list.push_back(new_xic_coord);
    }

  MsRunXicExtractorInterfaceSp xic_extractor =
    MsRunXicExtractorFactory::getInstance().buildMsRunXicExtractorSp(
      p_msreader);
  xic_extractor.get()->setXicExtractMethod(Enums::XicExtractMethod::sum);

  QTextStream outputStream(stdout, QIODevice::WriteOnly);
  UiMonitorText monitor(outputStream);
  INFO("monitor.setStatus");
  monitor.setStatus("Now starting the actual XIC extraction");

  INFO("extractXicCoordSPtrList start");

  // DO THE EXTRACTION!!
  xic_extractor.get()->extractXicCoordSPtrList(monitor, xic_coord_list);

  INFO("extractXicCoordSPtrList stop");

  qInfo() << QString("Time taken: %1\n")
               .arg((double)(clock() - start) / CLOCKS_PER_SEC);

  // At this point we have tens of thousands of XICs performed each time over
  // the whole set of frames of the acquisition.

  // A XIC is a total ion current chromatogram for a given m/z value over the
  // whole rt range. The XIC is a Trace, with x:rt y:int. And there is such a
  // Trace for each m/z value in the mz_list above.

  // If we want a full mass spectral integration over the whole acquisiion, all
  // we need to do is iterate in the xic_coord_list, get the m/z value for each
  // item in this list and sum all the y values of the XIC Trace.

  Trace full_mass_spectrum_trace;

  for(auto xic_coord : xic_coord_list)
    {
      double mz    = xic_coord->mzRange.getMz();
      double sum_y = xic_coord.get()->xicSptr.get()->sumY();

      full_mass_spectrum_trace.push_back(DataPoint(mz, sum_y));
    }

  std::cout << "The full mass spectrum for the acquisition:\n"
            << full_mass_spectrum_trace.toString().toStdString();

#elif USE_PAPPSO_TREE == 1

  std::cout << std::endl << "..:: NO test TIMS TDF parsing ::.." << std::endl;

#endif
}