text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2018.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg$
// $Authors: Stephan Aiche, Chris Bielow$
// --------------------------------------------------------------------------

#include <OpenMS/SIMULATION/RawMSSignalSimulation.h>
#include <OpenMS/FORMAT/TextFile.h>
#include <OpenMS/FORMAT/SVOutStream.h>
#include <OpenMS/SYSTEM/File.h>

#include <boost/random/uniform_real.hpp>
#include <boost/random/poisson_distribution.hpp>
#include <boost/random/exponential_distribution.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/math/distributions.hpp>

#ifdef _OPENMP
#include <omp.h>

#endif

namespace OpenMS
{

  /**
   * TODO: review baseline and noise code
   */
  RawMSSignalSimulation::RawMSSignalSimulation(SimTypes::MutableSimRandomNumberGeneratorPtr rng) :
    DefaultParamHandler("RawSignalSimulation"),
    ProgressLogger(),
    mz_error_mean_(),
    mz_error_stddev_(),
    intensity_scale_(),
    intensity_scale_stddev_(),
    res_model_(RES_CONSTANT),
    res_base_(0),
    rnd_gen_(rng),
    contaminants_(),
    contaminants_loaded_(false)
  {
    setDefaultParams_();
    updateMembers_();
  }

  RawMSSignalSimulation::RawMSSignalSimulation() :
    DefaultParamHandler("RawSignalSimulation"),
    mz_error_mean_(),
    mz_error_stddev_(),
    intensity_scale_(),
    intensity_scale_stddev_(),
    res_model_(RES_CONSTANT),
    res_base_(0),
    rnd_gen_(),
    contaminants_(),
    contaminants_loaded_(false)
  {
    setDefaultParams_();
    updateMembers_();
  }

  RawMSSignalSimulation::RawMSSignalSimulation(const RawMSSignalSimulation& source) :
    DefaultParamHandler(source),
    ProgressLogger(source),
    mz_error_mean_(source.mz_error_mean_),
    mz_error_stddev_(source.mz_error_stddev_),
    intensity_scale_(source.intensity_scale_),
    intensity_scale_stddev_(source.intensity_scale_stddev_),
    res_model_(source.res_model_),
    res_base_(source.res_base_),
    contaminants_(),
    contaminants_loaded_(false)
  {
    setParameters(source.getParameters());
    rnd_gen_ = source.rnd_gen_;
    updateMembers_();
  }

  RawMSSignalSimulation& RawMSSignalSimulation::operator=(const RawMSSignalSimulation& source)
  {
    setParameters(source.getParameters());
    rnd_gen_ = source.rnd_gen_;

    mz_error_mean_ = source.mz_error_mean_;
    mz_error_stddev_ = source.mz_error_stddev_;

    intensity_scale_ = source.intensity_scale_;
    intensity_scale_stddev_ = source.intensity_scale_stddev_;

    res_model_ = source.res_model_;
    res_base_ = source.res_base_;

    contaminants_ = source.contaminants_;
    contaminants_loaded_ = source.contaminants_loaded_;

    updateMembers_();
    return *this;
  }

  RawMSSignalSimulation::~RawMSSignalSimulation()
  {
  }

  void RawMSSignalSimulation::setDefaultParams_()
  {
    defaults_.setValue("enabled", "true", "Enable RAW signal simulation? (select 'false' if you only need feature-maps)");
    defaults_.setValidStrings("enabled", ListUtils::create<String>("true,false"));

    defaults_.setValue("ionization_type", "ESI", "Type of Ionization (MALDI or ESI)");
    defaults_.setValidStrings("ionization_type", ListUtils::create<String>("MALDI,ESI"));

    // peak and instrument parameter
    defaults_.setValue("resolution:value", 50000, "Instrument resolution at 400 Th.");
    defaults_.setValue("resolution:type", "linear", "How does resolution change with increasing m/z?! QTOFs usually show 'constant' behavior, FTs have linear degradation, and on Orbitraps the resolution decreases with square root of mass.");
    defaults_.setValidStrings("resolution:type", ListUtils::create<String>("constant,linear,sqrt"));

    defaults_.setValue("peak_shape", "Gaussian", "Peak Shape used around each isotope peak (be aware that the area under the curve is constant for both types, but the maximal height will differ (~ 2:3 = Lorentz:Gaussian) due to the wider base of the Lorentzian.");
    defaults_.setValidStrings("peak_shape", ListUtils::create<String>("Gaussian,Lorentzian"));

    // baseline
    defaults_.setValue("baseline:scaling", 0.0, "Scale of baseline. Set to 0 to disable simulation of baseline.");
    defaults_.setMinFloat("baseline:scaling", 0.0);
    defaults_.setValue("baseline:shape", 0.5, "The baseline is modeled by an exponential probability density function (pdf) with f(x) = shape*e^(- shape*x)");
    defaults_.setMinFloat("baseline:shape", 0.0);
    defaults_.setSectionDescription("baseline", "Baseline modeling for MALDI ionization");

    // mz sampling rate
    //       e.g. http://www.adronsystems.com/faqs.htm#rate states 8 points per peak on low-res instruments --> ~4 points at FWHM
    defaults_.setValue("mz:sampling_points", 3, "Number of raw data points per FWHM of the peak.");
    defaults_.setMinInt("mz:sampling_points", 2);

    // contaminants:
    defaults_.setValue("contaminants:file", "examples/simulation/contaminants.csv", "Contaminants file with sum formula and absolute RT interval. See 'OpenMS/examples/simulation/contaminants.txt' for details.");

    // VARIATION

    // m/z error
    // todo: also plan for affine trafo (as in RT shift?)
    defaults_.setValue("variation:mz:error_stddev", 0.0, "Standard deviation for m/z errors. Set to 0 to disable simulation of m/z errors.");
    defaults_.setValue("variation:mz:error_mean", 0.0, "Average systematic m/z error (Da)");

    defaults_.setValue("variation:intensity:scale", 100.0, "Constant scale factor of the feature intensity. Set to 1.0 to get the real intensity values provided in the FASTA file.");
    defaults_.setMinFloat("variation:intensity:scale", 0.0);
    defaults_.setValue("variation:intensity:scale_stddev", 0.0, "Standard deviation of peak intensity (relative to the scaled peak height). Set to 0 to get simple rescaled intensities.");
    defaults_.setMinFloat("variation:intensity:scale_stddev", 0.0);

    defaults_.setSectionDescription("variation:mz", "Shifts in mass to charge dimension of the simulated signals.");
    defaults_.setSectionDescription("variation:intensity", "Variations in intensity to model randomness in feature intensity.");
    defaults_.setSectionDescription("variation", "Random components that simulate biological and technical variations of the simulated data.");

    // NOISE

    // shot noise
    defaults_.setValue("noise:shot:rate", 0.0, "Poisson rate of shot noise per unit m/z. Set this to 0 to disable simulation of shot noise.");
    defaults_.setMinFloat("noise:shot:rate", 0.0);
    defaults_.setValue("noise:shot:intensity-mean", 1.0, "Shot noise intensity mean (exponentially distributed with given mean).");
    defaults_.setSectionDescription("noise:shot", "Parameters of Poisson and Exponential for shot noise modeling (set :rate OR :mean = 0 to disable).");

    // white noise
    defaults_.setValue("noise:white:mean", 0.0, "Mean value of white noise being added to each measured signal.");
    defaults_.setValue("noise:white:stddev", 0.0, "Standard deviation of white noise being added to each measured signal.");
    defaults_.setSectionDescription("noise:white", "Parameters of Gaussian distribution for white noise modeling (set :mean AND :stddev = 0 to disable).");

    // detector noise
    defaults_.setValue("noise:detector:mean", 0.0, "Mean value of the detector noise being added to the complete measurement.");
    defaults_.setValue("noise:detector:stddev", 0.0, "Standard deviation of the detector noise being added to the complete measurement.");
    defaults_.setSectionDescription("noise:detector", "Parameters of Gaussian distribution for detector noise modeling (set :mean AND :stddev = 0 to disable).");

    defaults_.setSectionDescription("noise", "Parameters modeling noise in mass spectrometry measurements.");

    defaultsToParam_();
  }

  double RawMSSignalSimulation::getResolution_(const double query_mz, const double resolution, const RESOLUTIONMODEL model) const
  {
    switch (model)
    {
    case RES_CONSTANT:
      return resolution;

    case RES_LINEAR:
      return resolution * (400 / query_mz);

    case RES_SQRT:
      return resolution * (std::sqrt(400.0) / sqrt(query_mz));

    default:
      throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Unknown RESOLUTIONMODEL encountered!");
    }
  }

  void RawMSSignalSimulation::updateMembers_()
  {
    res_base_ = (double) param_.getValue("resolution:value");
    String model = param_.getValue("resolution:type");
    if (model == "constant")
      res_model_ = RES_CONSTANT;
    else if (model == "linear")
      res_model_ = RES_LINEAR;
    else if (model == "sqrt")
      res_model_ = RES_SQRT;
    else
      throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Resolution:type given in parameters is unknown");

    sampling_points_per_FWHM_ = (Int) param_.getValue("mz:sampling_points") - 1;

    mz_error_mean_    = param_.getValue("variation:mz:error_mean");
    mz_error_stddev_  = param_.getValue("variation:mz:error_stddev");

    intensity_scale_ = param_.getValue("variation:intensity:scale");
    intensity_scale_stddev_ = param_.getValue("variation:intensity:scale_stddev");

    contaminants_loaded_ = false;
  }

  void RawMSSignalSimulation::loadContaminants()
  {
    // contaminants:
    String contaminants_file = param_.getValue("contaminants:file");

    if (contaminants_file.trim().size() != 0)
    {
      if (!File::readable(contaminants_file)) // look in OPENMS_DATA_PATH
      {
        contaminants_file = File::find(contaminants_file);
      }
      if (!File::readable(contaminants_file))
        throw Exception::FileNotReadable(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, contaminants_file);
      // read & parse file:
      TextFile tf(contaminants_file, true);
      contaminants_.clear();
      const UInt COLS_EXPECTED = 8;
      Size line_number = 1;
      for (TextFile::ConstIterator tf_it = tf.begin(); tf_it != tf.end(); ++tf_it, ++line_number)
      {
        if (tf_it->empty() || tf_it->hasPrefix("#")) continue; // skip comments

        String line = *tf_it;
        StringList cols;
        line.removeWhitespaces().split(',', cols, true);
        if (cols.size() != COLS_EXPECTED)
        {
          throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, line, "Expected " + String(COLS_EXPECTED) + " components, got " + String(cols.size()));
        }
        ContaminantInfo c;
        c.name = cols[0];
        try
        {
          c.sf = EmpiricalFormula(cols[1]);
          if (c.sf.getCharge() != 0)
          {
            throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, cols[1], "Line " + String(line_number) + " in " + contaminants_file + " contains forbidden charged sum formulas. Charges must be specified in another column. Remove all '+' or '-'!");
          }
        }
        catch (...)
        {
          throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, cols[1], "Could not parse line " + String(line_number) + " in " + contaminants_file + ".");
        }
        try
        {
          c.rt_start = cols[2].toDouble();
          c.rt_end = cols[3].toDouble();
          c.intensity = cols[4].toDouble();
          c.q = cols[5].toInt();
        }
        catch (...)
        {
          throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, line, "Could not parse line " + String(line_number) + " in " + contaminants_file + ".");
        }
        if (cols[6].toUpper() == "REC")
        {
          c.shape = RT_RECTANGULAR;
        }
        else if (cols[6].toUpper() == "GAUSS")
        {
          c.shape = RT_GAUSSIAN;
        }
        else
        {
          throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, line, "Unknown shape type: " + cols[6] + " in line " + String(line_number) + " of '" + contaminants_file + "'");
        }

        if (cols[7].toUpper() == "ESI")
        {
          c.im = IM_ESI;
        }
        else if (cols[7].toUpper() == "MALDI")
        {
          c.im = IM_MALDI;
        }
        else if (cols[7].toUpper() == "ALL")
        {
          c.im = IM_ALL;
        }
        else
        {
          throw Exception::ParseError(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, line, "Unknown ionization type: " + cols[7] + " in line " + String(line_number) + " of '" + contaminants_file + "'");
        }

        contaminants_.push_back(c);
      }
    }
    contaminants_loaded_ = true;
  }

  void RawMSSignalSimulation::generateRawSignals(SimTypes::FeatureMapSim& features, SimTypes::MSSimExperiment& experiment, SimTypes::MSSimExperiment& experiment_ct, SimTypes::FeatureMapSim& c_map)
  {
    OPENMS_LOG_INFO << "Raw MS1 Simulation ... ";
    // TODO: check if signal intensities scale linear with actual abundance, e.g. DOI: 10.1021/ac0202280 for NanoFlow-ESI

    // we rely on the same size of Raw and Peak Map
    if (experiment.size() != experiment_ct.size())
    {
      throw Exception::InvalidSize(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, experiment_ct.size());
    }

    if (param_.getValue("enabled") == "false")
    {
      OPENMS_LOG_INFO << "disabled" << std::endl;
      return;
    }
    else
    {
      OPENMS_LOG_INFO << "started" << std::endl;
    }

    // retrieve mz boundary parameters from experiment:
    SimTypes::SimCoordinateType minimal_mz_measurement_limit = experiment[0].getInstrumentSettings().getScanWindows()[0].begin;
    SimTypes::SimCoordinateType maximal_mz_measurement_limit = experiment[0].getInstrumentSettings().getScanWindows()[0].end;

    // grid is constant over scans, so we compute it only once
    getSamplingGrid_(grid_, minimal_mz_measurement_limit, maximal_mz_measurement_limit, 5); // every 5 Da we adjust the sampling width by local FWHM

    OPENMS_LOG_INFO << "  Simulating signal for " << features.size() << " features ..." << std::endl;

    this->startProgress(0, features.size(), "RawMSSignal");

    Size progress(0);
    // we have a bit of code duplication here but this eases the
    // parallelization step
    if (experiment.size() == 1) // MS only
    {
      for (FeatureMap::iterator feature_it = features.begin();
           feature_it != features.end();
           ++feature_it, ++progress)
      {
        add1DSignal_(*feature_it, experiment, experiment_ct);
        this->setProgress(progress);
      }
    }
    else // LC/MS
    {
      std::vector<SimTypes::MSSimExperiment*> experiments; // pointer to experiment(s)
      experiments.push_back(&experiment); // the master thread gets the original (just a reference, no copying here)

      std::vector<SimTypes::MSSimExperiment*> experiments_ct; // pointer to experiment(s)
      experiments_ct.push_back(&experiment_ct); // the master thread gets the original (just a reference, no copying here)

#ifdef _OPENMP
      // prepare random numbers for the different threads
      // each possible thread gets his own set of random
      // numbers
      Size thread_count = omp_get_max_threads();

      threaded_random_numbers_.resize(thread_count);
      threaded_random_numbers_index_.resize(thread_count);
      experiments.reserve(thread_count); // !reserve!
      experiments_ct.reserve(thread_count); // !reserve!
      std::vector<SimTypes::MSSimExperiment> experiments_tmp(thread_count - 1); // holds MSExperiments for slave threads
      std::vector<SimTypes::MSSimExperiment> experiments_ct_tmp(thread_count - 1); // holds MSExperiments (centroided) for slave threads

      for (Size i = 0; i < thread_count; ++i)
      {
        threaded_random_numbers_[i].resize(THREADED_RANDOM_NUMBER_POOL_SIZE_);
        threaded_random_numbers_index_[i] = THREADED_RANDOM_NUMBER_POOL_SIZE_;
      }

      if (thread_count > 1)
      {
        // prepare a temporary experiment to store the results
        SimTypes::MSSimExperiment e_tmp = experiment;
        SimTypes::MSSimExperiment e_ct_tmp = experiment_ct;
        // remove actual data
        for (Size i = 0; i < e_tmp.size(); ++i)
        {
          e_tmp[i].clear(false);
          e_ct_tmp[i].clear(false);
        }
        // each slave thread gets a copy
        for (Size i = 1; i < thread_count; ++i)
        {
          experiments_tmp[i - 1] = e_tmp;
          experiments_ct_tmp[i - 1] = e_ct_tmp;
          // assign it to the list of experiments (this is no real copy, but a reference!)
          experiments.push_back(&(experiments_tmp[i - 1]));
          experiments_ct.push_back(&(experiments_ct_tmp[i - 1]));
        }
      }
#else
      Size thread_count = 1;
#endif

      Size compress_size_intermediate = 20000 / thread_count; // compress map every X features, (10.000 feature are ~ 2 GB at 0.002 sampling rate)
      Size compress_count = 0; // feature count (for each thread)

#ifdef _OPENMP
#pragma omp parallel for firstprivate(compress_count)
#endif
      for (SignedSize f = 0; f < (SignedSize)features.size(); ++f)
      {
#ifdef _OPENMP // update experiment index if necessary
        const int current_thread = omp_get_thread_num();
#else
        const int current_thread(0);
#endif
        add2DSignal_(features[f], *(experiments[current_thread]), *(experiments_ct[current_thread]));

        // progresslogger, only master thread sets progress (no barrier here)
#ifdef _OPENMP
#pragma omp atomic
#endif
        ++progress;
        if (current_thread == 0)
        {
          this->setProgress(progress);
        }

        // intermediate compress to avoid memory problems
        ++compress_count;
        if (compress_count > compress_size_intermediate)
        {
          compress_count = 0;
          compressSignals_(*(experiments[current_thread]));
        }
      } // ! raw signal sim

#ifdef _OPENMP // merge back other experiments
      for (Size i = 1; i < experiments.size(); ++i)
      {
        // copy peak data from temporal experiment
        for (Size scan = 0; scan < experiment.size(); ++scan)
        {
          if ((*experiments[i])[scan].empty())
            continue; // we do not care if the spectrum wasn't touched at all
          // append all points from temp to org
          experiment[scan].insert(experiment[scan].end(), (*experiments[i])[scan].begin(), (*experiments[i])[scan].end());
          // delete from child experiment to save memory (otherwise the merge would double it!)
          (*experiments[i])[scan].clear(false);

          // peak GT ( small, so no need to compress)
          experiment_ct[scan].insert(experiment_ct[scan].end(), (*experiments_ct[i])[scan].begin(), (*experiments_ct[i])[scan].end());

        }
      }
#endif

    } // ! 1D or 2D

    this->endProgress();

    // finally sort generated data
    experiment.sortSpectra(true);
    experiment.updateRanges();

    // build contaminant feature map & add raw signal
    if (experiment.size() > 1) // LC/MS only currently
    {
      createContaminants_(c_map, experiment, experiment_ct);
    }

    if ((String)param_.getValue("ionization_type") == "MALDI")
    {
      addBaseLine_(experiment, minimal_mz_measurement_limit);
    }
    addShotNoise_(experiment, minimal_mz_measurement_limit, maximal_mz_measurement_limit);
    compressSignals_(experiment);

    // add white noise to the simulated data
    addWhiteNoise_(experiment);

    // add detector noise the simulated data
    addDetectorNoise_(experiment);
  }

  double RawMSSignalSimulation::getPeakWidth_(const double mz, const bool is_gaussian) const
  {
    double mz_local = std::max(mz, 400.0); // at least assume m/z=400, as otherwise FWHM might get ridiculously small
    // convert from resolution @ current m/z --> FWHM
    double fwhm = mz_local / getResolution_(mz_local, res_base_, res_model_);
    // Approximation for Gaussian-shaped signals,
    // i.e. sqrt(2*ln(2))*2 = 2.35482
    // , relating FWHM to Gaussian width
    if (is_gaussian)
      fwhm /= 2.35482;
    else
    {
    } // for Lorentzian, we do nothing as the scale parameter is exactly the FWHM
    return fwhm;
  }

  void RawMSSignalSimulation::add1DSignal_(Feature& active_feature, SimTypes::MSSimExperiment& experiment, SimTypes::MSSimExperiment& experiment_ct)
  {
    SimTypes::SimIntensityType scale = getFeatureScaledIntensity_(active_feature.getIntensity(), 100.0);

    SimTypes::SimChargeType q = active_feature.getCharge();
    EmpiricalFormula ef = active_feature.getPeptideIdentifications()[0].getHits()[0].getSequence().getFormula();
    ef += EmpiricalFormula(active_feature.getMetaValue("charge_adducts")); // adducts
    ef -= EmpiricalFormula(String("H") + String(q));
    ef.setCharge(q); // effectively subtract q electrons

    Param p1;
    p1.setValue("statistics:mean", ef.getAverageWeight() / q);
    p1.setValue("interpolation_step", 0.001);
    p1.setValue("isotope:mode:mode", param_.getValue("peak_shape"));
    p1.setValue("intensity_scaling", 0.001 * scale); // this removes the problem of to big isotope-model values
    p1.setValue("charge", q);
    double fwhm;
    if (param_.getValue("peak_shape") == "Gaussian")
    {
      fwhm = getPeakWidth_(active_feature.getMZ(), true);
      p1.setValue("isotope:mode:GaussianSD", fwhm);
    }
    else
    {
      fwhm = getPeakWidth_(active_feature.getMZ(), false);
      p1.setValue("isotope:mode:LorentzFWHM", fwhm);
    }

    IsotopeModel isomodel;
    isomodel.setParameters(p1);
    isomodel.setSamples(ef);

    SimTypes::SimCoordinateType mz_start = isomodel.getInterpolation().supportMin();
    SimTypes::SimCoordinateType mz_end = isomodel.getInterpolation().supportMax();

    samplePeptideModel1D_(isomodel, mz_start, mz_end, experiment, experiment_ct, active_feature);
  }

  void RawMSSignalSimulation::add2DSignal_(Feature& active_feature, SimTypes::MSSimExperiment& experiment, SimTypes::MSSimExperiment& experiment_ct)
  {
    SimTypes::SimIntensityType scale = getFeatureScaledIntensity_(active_feature.getIntensity(), 1.0);

    SimTypes::SimChargeType q = active_feature.getCharge();
    EmpiricalFormula ef;
    if (active_feature.metaValueExists("sum_formula"))
    {
      ef = EmpiricalFormula(active_feature.getMetaValue("sum_formula"));
    }
    else
    {
      ef = EmpiricalFormula(active_feature.getPeptideIdentifications()[0].getHits()[0].getSequence().getFormula());
    }
    ef += EmpiricalFormula(active_feature.getMetaValue("charge_adducts")); // adducts
    ef -= EmpiricalFormula(String("H") + String(q));
    ef.setCharge(q); // effectively subtract q electrons

    Param p1;
    p1.setValue("statistics:mean", ef.getAverageWeight() / q);
    p1.setValue("interpolation_step", 0.001);
    p1.setValue("isotope:mode:mode", param_.getValue("peak_shape"));
    p1.setValue("intensity_scaling", 0.001); // this removes the problem of to big isotope-model values
    p1.setValue("charge", q);
    double fwhm;
    if (param_.getValue("peak_shape") == "Gaussian")
    {
      fwhm = getPeakWidth_(active_feature.getMZ(), true);
      p1.setValue("isotope:mode:GaussianSD", fwhm);
    }
    else
    {
      fwhm = getPeakWidth_(active_feature.getMZ(), false);
      p1.setValue("isotope:mode:LorentzFWHM", fwhm);
    }

    IsotopeModel* isomodel = new IsotopeModel();
    isomodel->setParameters(p1); // this needs to come BEFORE setSamples() - otherwise the default setSamples() is called here!
    isomodel->setSamples(ef); // this already includes adducts

    if (experiment.size() < 2)
    {
      throw Exception::InvalidSize(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, experiment.size());
    }
    double rt_sampling_rate = experiment[1].getRT() - experiment[0].getRT();
    EGHModel* elutionmodel = new EGHModel();
    chooseElutionProfile_(elutionmodel, active_feature, 1.0, rt_sampling_rate, experiment);
    ProductModel<2> pm;
    pm.setModel(0, elutionmodel); // new'ed models will be deleted by the pm! no need to delete them manually
    pm.setModel(1, isomodel); // new'ed models will be deleted by the pm! no need to delete them manually
    pm.setScale(scale); // scale

    // start and end points of the sampling
    SimTypes::SimCoordinateType rt_start(elutionmodel->getInterpolation().supportMin());
    SimTypes::SimCoordinateType rt_end(elutionmodel->getInterpolation().supportMax());
    if (active_feature.metaValueExists("RT_width_start") && active_feature.metaValueExists("RT_width_end")) // this is a contaminant with sampling restrictions
    {
      rt_start = active_feature.getMetaValue("RT_width_start");
      rt_end = active_feature.getMetaValue("RT_width_end");
    }
    SimTypes::SimCoordinateType mz_start(isomodel->getInterpolation().supportMin());
    SimTypes::SimCoordinateType mz_end(isomodel->getInterpolation().supportMax());

    // add peptide to GLOBAL MS map
    // add CH and new intensity to feature
    samplePeptideModel2D_(pm, mz_start, mz_end, rt_start, rt_end, experiment, experiment_ct, active_feature);
  }

  void RawMSSignalSimulation::samplePeptideModel1D_(const IsotopeModel& pm,
                                                    const SimTypes::SimCoordinateType mz_start,
                                                    const SimTypes::SimCoordinateType mz_end,
                                                    SimTypes::MSSimExperiment& experiment,
                                                    SimTypes::MSSimExperiment& experiment_ct,
                                                    Feature& active_feature)
  {
    SimTypes::SimIntensityType intensity_sum = 0.0;

    //OPENMS_LOG_DEBUG << "Sampling at [mz] " << mz_start << ":" << mz_end << std::endl;

    SimTypes::SimPointType point;

    // centroided GT
    for (IsotopeDistribution::const_iterator iter = pm.getIsotopeDistribution().begin();
         iter != pm.getIsotopeDistribution().end(); ++iter)
    {
      point.setMZ(iter->getMZ());
      point.setIntensity(iter->getIntensity());

      if (point.getIntensity() <= 0.0)
        continue;

      experiment_ct[0].push_back(point);
    }

    std::vector<SimTypes::SimCoordinateType>::const_iterator it_grid = lower_bound(grid_.begin(), grid_.end(), mz_start);
    boost::normal_distribution<double> ndist(mz_error_mean_, mz_error_stddev_);
    for (; it_grid != grid_.end() && (*it_grid) < mz_end; ++it_grid)
    {
      point.setMZ(*it_grid);
      point.setIntensity(pm.getIntensity(DPosition<1>(*it_grid)));

      if (point.getIntensity() <= 0.0)
        continue;

      // add Gaussian distributed m/z error
      double mz_err = ndist(rnd_gen_->getTechnicalRng());
      point.setMZ(fabs(point.getMZ() + mz_err));

      intensity_sum += point.getIntensity();
      experiment[0].push_back(point);
    }
    active_feature.setIntensity(intensity_sum);
  }

  void RawMSSignalSimulation::samplePeptideModel2D_(const ProductModel<2>& pm,
                                                    const SimTypes::SimCoordinateType mz_start,
                                                    const SimTypes::SimCoordinateType mz_end,
                                                    SimTypes::SimCoordinateType rt_start,
                                                    SimTypes::SimCoordinateType rt_end,
                                                    SimTypes::MSSimExperiment& experiment,
                                                    SimTypes::MSSimExperiment& experiment_ct,
                                                    Feature& active_feature)
  {
    if (rt_start <= 0)
      rt_start = 0;

    SimTypes::MSSimExperiment::iterator exp_start = experiment.RTBegin(rt_start);
    SimTypes::MSSimExperiment::iterator exp_ct_start = experiment_ct.RTBegin(rt_start);

    if (exp_start == experiment.end())
    {
      throw Exception::InvalidSize(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, 0);
    }

    SimTypes::SimIntensityType intensity_sum(0.0);

#ifdef OPENMS_ASSERTIONS
    Int end_scan = std::numeric_limits<Int>::min(); // only used in Debug build
#endif

    IsotopeModel* isomodel = static_cast<IsotopeModel*>(pm.getModel(1));
    IsotopeDistribution iso_dist = isomodel->getIsotopeDistribution();
    SimTypes::SimCoordinateType mz_mono = active_feature.getMZ();
    SimTypes::SimCoordinateType iso_peakdist = isomodel->getParameters().getValue("isotope:distance");
    Int q = active_feature.getCharge();

    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // Sample the model ...
    SimTypes::SimCoordinateType rt(0);
    SimTypes::MSSimExperiment::iterator exp_iter = exp_start;
    SimTypes::MSSimExperiment::iterator exp_ct_iter = exp_ct_start;
    for (; rt < rt_end && exp_iter != experiment.end(); ++exp_iter, ++exp_ct_iter)
    {
      rt = exp_iter->getRT();
      double distortion = double(exp_iter->getMetaValue("distortion"));
      double rt_intensity = ((EGHModel*)pm.getModel(0))->getIntensity(rt);

      // centroided GT
      Size iso_pos(0);
      SimTypes::SimPointType point;
      for (IsotopeDistribution::const_iterator iter = iso_dist.begin(); iter != iso_dist.end(); ++iter, ++iso_pos)
      {
        point.setMZ(mz_mono + (iso_pos * iso_peakdist / q));
        point.setIntensity(iter->getIntensity() * rt_intensity * distortion);

        if (point.getIntensity() <= 0.0)
          continue;

        exp_ct_iter->push_back(point);
      }

      // RAW signal (sample it on the grid)
      std::vector<SimTypes::SimCoordinateType>::const_iterator it_grid = lower_bound(grid_.begin(), grid_.end(), mz_start);
      for (; it_grid != grid_.end() && (*it_grid) < mz_end; ++it_grid)
      {
        ProductModel<2>::IntensityType intensity = pm.getIntensity(DPosition<2>(rt, *it_grid)) * distortion;
        if (intensity <= 0.0)
          continue; // intensity cutoff (below that we don't want to see a signal)

        point.setMZ(*it_grid);
        point.setIntensity(intensity);

        //OPENMS_LOG_ERROR << "Sampling " << rt << " , " << mz << " -> " << point.getIntensity() << std::endl;

        // add Gaussian distributed m/z error
#ifdef _OPENMP
        int CURRENT_THREAD = omp_get_thread_num();
        // check if we need to refill the random number pool for this thread
        if (threaded_random_numbers_index_[CURRENT_THREAD] == THREADED_RANDOM_NUMBER_POOL_SIZE_)
        {
          if (mz_error_stddev_ != 0.0)
          {
#pragma omp critical(generate_random_number_for_thread)
            {
              boost::normal_distribution<double> ndist(mz_error_mean_, mz_error_stddev_);
              for (Size i = 0; i < THREADED_RANDOM_NUMBER_POOL_SIZE_; ++i)
              {
                threaded_random_numbers_[CURRENT_THREAD][i] = ndist(rnd_gen_->getTechnicalRng());
              }
            }
          }
          else
          {
            // we do not need to care about concurrency here
            fill(threaded_random_numbers_[CURRENT_THREAD].begin(), threaded_random_numbers_[CURRENT_THREAD].end(), mz_error_mean_);
          }
          // reset index for this thread to first position
          threaded_random_numbers_index_[CURRENT_THREAD] = 0;
        }

        const double mz_err = threaded_random_numbers_[CURRENT_THREAD][threaded_random_numbers_index_[CURRENT_THREAD]++];
#else
        // we can use the normal Gaussian ran-gen if we do not use OPENMP
        boost::normal_distribution<double> ndist(mz_error_mean_, mz_error_stddev_);
        const double mz_err = ndist(rnd_gen_->getTechnicalRng());
#endif
        point.setMZ(std::fabs(point.getMZ() + mz_err));
        exp_iter->push_back(point);

        intensity_sum += point.getIntensity();
      }
      //update last scan affected
#ifdef OPENMS_ASSERTIONS
      end_scan = exp_iter - experiment.begin();
#endif
    }

    OPENMS_POSTCONDITION(end_scan != std::numeric_limits<Int>::min(), "RawMSSignalSimulation::samplePeptideModel2D_(): setting RT bounds failed!");

    // new intensity is AREA==SUM of all peaks
    active_feature.setIntensity(intensity_sum);

    // -------------------------
    // --- store convex hull ---
    // -------------------------
    active_feature.getConvexHulls().clear();

    // use isotope model (to determine mass traces)

    DoubleList isotope_intensities;
    for (IsotopeDistribution::iterator iter = iso_dist.begin();
         iter != iso_dist.end(); ++iter)
    {
      const SimTypes::SimCoordinateType mz = mz_mono + double(iter->getMZ() - iso_dist.begin()->getMZ()) / q; // this is only an approximated trace' m/z position (as we do assume 1Da space between them)

      SimTypes::SimCoordinateType rt_min =  std::numeric_limits<SimTypes::SimCoordinateType>::max();
      SimTypes::SimCoordinateType rt_max = -std::numeric_limits<SimTypes::SimCoordinateType>::max();
      bool has_data = false;

      // for each trace, sample the model again and see how far it extends
      SimTypes::SimCoordinateType rt(0);
      for (exp_iter = exp_start; rt < rt_end && exp_iter != experiment.end(); ++exp_iter)
      {
        rt = exp_iter->getRT();
        double distortion = double(exp_iter->getMetaValue("distortion"));
        ProductModel<2>::IntensityType intensity = pm.getIntensity(DPosition<2>(rt, mz)) * distortion;
        if (intensity <= 0.0)
          continue; // intensity cutoff (below that we don't want to see a signal)

        // update min&max
        if (rt_min > rt)
          rt_min = rt;
        if (rt_max < rt)
          rt_max = rt;
        has_data = true;
      }
      if (!has_data)
        continue;

      // add four edge points of mass trace
      ConvexHull2D hull;
      std::vector<DPosition<2> > points;
      points.push_back(DPosition<2>(rt_min, mz - 0.001));
      points.push_back(DPosition<2>(rt_min, mz + 0.001));
      points.push_back(DPosition<2>(rt_max, mz - 0.001));
      points.push_back(DPosition<2>(rt_max, mz + 0.001));
      hull.addPoints(points);
      active_feature.getConvexHulls().push_back(hull);

      isotope_intensities.push_back(iter->getIntensity());
    }

    active_feature.setMetaValue("isotope_intensities", isotope_intensities);

  }

  void RawMSSignalSimulation::chooseElutionProfile_(EGHModel* const elutionmodel, Feature& feature, const double scale, const double rt_sampling_rate, const SimTypes::MSSimExperiment& experiment)
  {
    SimTypes::SimCoordinateType f_rt = feature.getRT();

    Param p;
    // WARNING: step used to be 'rt_sampling_rate / 3.0', but distortion is not part of RT sim, and thus only
    //          modeled 1:1
    p.setValue("interpolation_step", rt_sampling_rate);
    p.setValue("statistics:variance", 1.0);
    p.setValue("statistics:mean", f_rt);

    p.setValue("egh:height", scale);
    p.setValue("egh:retention", f_rt);

    if (feature.metaValueExists("RT_width_gaussian")) // this is for contaminants only (we want the gaussian distribution width A+B at 5% of maximal height)
    {
      p.setValue("egh:alpha", 0.05);
      p.setValue("egh:A", double(feature.getMetaValue("RT_width_gaussian")) / 2.0 * 0.9); // make width a little smaller as this is only the 5% height cutoff
      p.setValue("egh:B", double(feature.getMetaValue("RT_width_gaussian")) / 2.0 * 0.9);
    }
    else if (feature.metaValueExists("RT_egh_variance") && feature.metaValueExists("RT_egh_tau"))
    {
      // for CE we want wider profiles with higher MT
      double width_factor(1); // default for HPLC
      if (feature.metaValueExists("RT_CE_width_factor"))
        width_factor = feature.getMetaValue("RT_CE_width_factor");

      p.setValue("egh:guess_parameter", "false");
      p.setValue("egh:tau", (double) feature.getMetaValue("RT_egh_tau"));
      p.setValue("egh:sigma_square", ((double) feature.getMetaValue("RT_egh_variance")) * width_factor);
    }
    else
    {
      throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Elution profile shape cannot be created. Wrong meta-values!", "");
    }

    elutionmodel->setParameters(p); // does the calculation
    //----------------------------------------------------------------------

    SimTypes::SimCoordinateType rt_em_start = elutionmodel->getInterpolation().supportMin();
    SimTypes::SimCoordinateType rt_em_end = elutionmodel->getInterpolation().supportMax();

    // find scan in experiment at which our elution starts
    SimTypes::MSSimExperiment::ConstIterator exp_it = experiment.RTBegin(rt_em_start);
    if (exp_it == experiment.end())
      --exp_it; // we need the last valid RT below, so .end() is not useful

    DoubleList elution_intensities;
    DoubleList elution_bounds;
    elution_bounds.resize(4); // store min and max RT (in seconds and index terms)
    elution_bounds[0] = std::distance(experiment.begin(), exp_it);
    elution_bounds[1] = exp_it->getRT();
    elution_bounds[2] = elution_bounds[0];
    elution_bounds[3] = elution_bounds[1];

    for (; (exp_it != experiment.end()) && (exp_it->getRT() <= rt_em_end); ++exp_it) // .. and disturb values by (an already smoothed) distortion diced in RTSimulation
    {
      double intensity = (double) exp_it->getMetaValue("distortion") * elutionmodel->getInterpolation().value(exp_it->getRT());
      // store elution profile in feature MetaValue
      elution_intensities.push_back(intensity);
      elution_bounds[2] = std::distance(experiment.begin(), exp_it);
      elution_bounds[3] = exp_it->getRT();
    }
    // set elution profile details in feature -> used for MS^E precursor selection in tandem MS later
    feature.setMetaValue("elution_profile_intensities", elution_intensities);
    feature.setMetaValue("elution_profile_bounds", elution_bounds);
  }

  void RawMSSignalSimulation::createContaminants_(SimTypes::FeatureMapSim& c_map, SimTypes::MSSimExperiment& exp, SimTypes::MSSimExperiment& exp_ct)
  {
    if (exp.size() == 1)
    {
      throw Exception::NotImplemented(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); // not implemented for 1D yet
    }

    if (!contaminants_loaded_)
      loadContaminants();

    IONIZATIONMETHOD this_im = (String)param_.getValue("ionization_type") == "ESI" ? IM_ESI : IM_MALDI;
    c_map.clear(true);

    Size out_of_range_RT(0), out_of_range_MZ(0);
    SimTypes::SimCoordinateType minimal_mz_measurement_limit = exp[0].getInstrumentSettings().getScanWindows()[0].begin;
    SimTypes::SimCoordinateType maximal_mz_measurement_limit = exp[0].getInstrumentSettings().getScanWindows()[0].end;

    for (Size i = 0; i < contaminants_.size(); ++i)
    {
      if (contaminants_[i].im != IM_ALL && contaminants_[i].im != this_im)
        continue;

      if (exp.getMinRT() > contaminants_[i].rt_end || contaminants_[i].rt_start > exp.getMaxRT())
      {
        ++out_of_range_RT;
        continue;
      }
      // ... create contaminants...
      SimTypes::FeatureMapSim::FeatureType feature;
      feature.setRT((contaminants_[i].rt_end + contaminants_[i].rt_start) / 2);
      feature.setMZ((contaminants_[i].sf.getMonoWeight() / contaminants_[i].q) + Constants::PROTON_MASS_U); // m/z (incl. protons)
      if (!(minimal_mz_measurement_limit < feature.getMZ() && feature.getMZ() < maximal_mz_measurement_limit))
      {
        ++out_of_range_MZ;
        continue;
      }
      feature.setIntensity(contaminants_[i].intensity);
      if (contaminants_[i].shape == RT_RECTANGULAR)
      {
        feature.setMetaValue("RT_width_gaussian", 1e6);
        feature.setMetaValue("RT_width_start", contaminants_[i].rt_start);
        feature.setMetaValue("RT_width_end", contaminants_[i].rt_end);
      }
      else
      {
        feature.setMetaValue("RT_width_gaussian", contaminants_[i].rt_end - contaminants_[i].rt_start);
      }
      feature.setMetaValue("sum_formula", contaminants_[i].sf.toString()); // formula without adducts or charges
      feature.setCharge(contaminants_[i].q);
      feature.setMetaValue("charge_adducts", "H" + String(contaminants_[i].q)); // adducts separately
      add2DSignal_(feature, exp, exp_ct);
      c_map.push_back(feature);
    }

    c_map.applyMemberFunction(&UniqueIdInterface::ensureUniqueId);
    OPENMS_LOG_INFO << "Contaminants out-of-RT-range: " << out_of_range_RT << " / " << contaminants_.size() << std::endl;
    OPENMS_LOG_INFO << "Contaminants out-of-MZ-range: " << out_of_range_MZ << " / " << contaminants_.size() << std::endl;

  }

  void RawMSSignalSimulation::addShotNoise_(SimTypes::MSSimExperiment& experiment, SimTypes::SimCoordinateType minimal_mz_measurement_limit, SimTypes::SimCoordinateType maximal_mz_measurement_limit)
  {
    const SimTypes::SimCoordinateType window_size = 100.0;
    SimTypes::SimCoordinateType mz_lw = minimal_mz_measurement_limit;
    SimTypes::SimCoordinateType mz_up = window_size + minimal_mz_measurement_limit;

    // we model the amount of (background) noise as Poisson process
    // i.e. the number of noise data points per unit m/z interval follows a Poisson
    // distribution. Noise intensity is assumed to be exponentially-distributed.
    double rate = param_.getValue("noise:shot:rate");
    double intensity_mean = param_.getValue("noise:shot:intensity-mean");

    // avoid sampling 0 values
    if (rate == 0.0 || intensity_mean == 0.0)
      return;

    // we distribute the rate in 100 Th windows
    double scaled_rate = rate * window_size;
    SimTypes::SimPointType shot_noise_peak;

    //distributions to sample from
    boost::random::poisson_distribution<UInt, double> pdist(scaled_rate);
    boost::uniform_real<SimTypes::SimCoordinateType> udist(mz_lw, mz_up);
    boost::random::exponential_distribution<SimTypes::SimCoordinateType> edist(intensity_mean);

    OPENMS_LOG_INFO << "Adding shot noise to spectra ..." << std::endl;
    Size num_intervals = std::ceil((maximal_mz_measurement_limit - minimal_mz_measurement_limit) / window_size);

    for (SimTypes::MSSimExperiment::Iterator spectrum_it = experiment.begin(); spectrum_it != experiment.end(); ++spectrum_it)
    {
      for (Size j = 0; j < num_intervals; ++j)
      {
        UInt counts = pdist(rnd_gen_->getTechnicalRng());
        for (UInt c = 0; c < counts; ++c)
        {
          SimTypes::SimCoordinateType mz = udist(rnd_gen_->getTechnicalRng());
          SimTypes::SimCoordinateType intensity = edist(rnd_gen_->getTechnicalRng());

          // we only add points if they have an intensity>0 and are inside of the measurement range
          if (mz < maximal_mz_measurement_limit)
          {
            shot_noise_peak.setIntensity(intensity);
            shot_noise_peak.setMZ(mz);
            spectrum_it->push_back(shot_noise_peak);
          }
        }

        mz_up += window_size;
      }
    } // end of each scan

    experiment.updateRanges();

  }

  void RawMSSignalSimulation::addBaseLine_(SimTypes::MSSimExperiment& experiment, SimTypes::SimCoordinateType minimal_mz_measurement_limit)
  {
    double scale = param_.getValue("baseline:scaling");
    double shape = param_.getValue("baseline:shape");

    if (scale == 0.0)
      return;

    // TODO: switch to iterator
    for (Size i = 0; i < experiment.size(); ++i)
    {
      for (Size j = 0; j < experiment[i].size(); ++j)
      {
        SimTypes::SimCoordinateType x = (experiment[i][j].getMZ() - minimal_mz_measurement_limit);
        //if (x >= 1000.0) continue; // speed-up TODO: revise this ..

        boost::math::exponential_distribution<double> ed(shape);
        double bx = boost::math::pdf(ed, x);
        bx *= scale;
        experiment[i][j].setIntensity(experiment[i][j].getIntensity() + bx);
      }
    }
  }

  void RawMSSignalSimulation::addWhiteNoise_(SimTypes::MSSimExperiment& experiment)
  {
    OPENMS_LOG_INFO << "Adding white noise to spectra ..." << std::endl;

    // get white noise parameters
    double white_noise_mean = param_.getValue("noise:white:mean");
    double white_noise_stddev = param_.getValue("noise:white:stddev");

    if (white_noise_mean == 0.0 && white_noise_stddev == 0.0)
    {
      return;
    }

    boost::normal_distribution<SimTypes::SimIntensityType> ndist(white_noise_mean, white_noise_stddev);

    for (SimTypes::MSSimExperiment::iterator spectrum_it = experiment.begin(); spectrum_it != experiment.end(); ++spectrum_it)
    {
      SimTypes::MSSimExperiment::SpectrumType new_spec = (*spectrum_it);
      new_spec.clear(false);

      for (SimTypes::MSSimExperiment::SpectrumType::iterator peak_it = (*spectrum_it).begin(); peak_it != (*spectrum_it).end(); ++peak_it)
      {
        SimTypes::SimIntensityType intensity = peak_it->getIntensity() + ndist(rnd_gen_->getTechnicalRng());
        if (intensity > 0.0)
        {
          peak_it->setIntensity(intensity);
          new_spec.push_back(*peak_it);
        }
      }

      *spectrum_it = new_spec;
    }
  }

  void RawMSSignalSimulation::addDetectorNoise_(SimTypes::MSSimExperiment& experiment)
  {
    OPENMS_LOG_INFO << "Adding detector noise to spectra ..." << std::endl;

    // get white noise parameters
    double detector_noise_mean = param_.getValue("noise:detector:mean");
    double detector_noise_stddev = param_.getValue("noise:detector:stddev");

    if (detector_noise_mean == 0.0 && detector_noise_stddev == 0.0)
    {
      OPENMS_LOG_INFO << "Detector noise was disabled." << std::endl;
      return;
    }

    boost::normal_distribution<SimTypes::SimIntensityType> ndist(detector_noise_mean, detector_noise_stddev);
    for (SimTypes::MSSimExperiment::iterator spectrum_it = experiment.begin(); spectrum_it != experiment.end(); ++spectrum_it)
    {
      SimTypes::MSSimExperiment::SpectrumType new_spec = (*spectrum_it);
      new_spec.clear(false);

      std::vector<SimTypes::SimCoordinateType>::iterator grid_it = grid_.begin();
      SimTypes::MSSimExperiment::SpectrumType::iterator peak_it = spectrum_it->begin();
      for (; grid_it != grid_.end(); ++grid_it)
      {
        // if peak is in grid
        if (peak_it != spectrum_it->end() && *grid_it == peak_it->getMZ())
        {
          SimTypes::SimIntensityType intensity = peak_it->getIntensity() + ndist(rnd_gen_->getTechnicalRng());
          if (intensity > 0.0)
          {
            peak_it->setIntensity(intensity);
            new_spec.push_back(*peak_it);
          }
          ++peak_it;
        }
        else // we have no point here, generate one if noise is above 0
        {
          SimTypes::SimIntensityType intensity = ndist(rnd_gen_->getTechnicalRng());
          if (intensity > 0.0)
          {
            SimTypes::MSSimExperiment::SpectrumType::PeakType noise_peak;
            noise_peak.setMZ(*grid_it);
            noise_peak.setIntensity(intensity);
            new_spec.push_back(noise_peak);
          }
        }
      }

      *spectrum_it = new_spec;
    }

  }

  void RawMSSignalSimulation::getSamplingGrid_(std::vector<SimTypes::SimCoordinateType>& grid, const SimTypes::SimCoordinateType mz_min, const SimTypes::SimCoordinateType mz_max, const Int step_Da)
  {
    if (fabs(mz_max - mz_min) < step_Da)
    {
      throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Sampling grid seems very small. This cannot be computed!");
    }
    grid.clear();
    SimTypes::SimCoordinateType mz = mz_min;
    double sampling_rate(0);
    while (mz <= mz_max)
    {
      SimTypes::SimCoordinateType fwhm = getPeakWidth_(mz, param_.getValue("peak_shape") == "Gaussian");
      sampling_rate = (fwhm / sampling_points_per_FWHM_);
      SimTypes::SimCoordinateType mz_cell_end = std::min(mz + step_Da, mz_max);
      for (; mz <= mz_cell_end; mz += sampling_rate)
      {
        grid.push_back(mz);
      }
    }
    grid.push_back(mz + sampling_rate); // one more point after mz_max, for binary search later
    return;
  }

  // TODO: add instrument specific sampling technique
  void RawMSSignalSimulation::compressSignals_(SimTypes::MSSimExperiment& experiment)
  {
    if (experiment.size() < 1 || experiment[0].getInstrumentSettings().getScanWindows().size() < 1)
    {
      throw Exception::IllegalSelfOperation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
    }
    SimTypes::SimCoordinateType min_mz = experiment[0].getInstrumentSettings().getScanWindows()[0].begin;
    SimTypes::SimCoordinateType max_mz = experiment[0].getInstrumentSettings().getScanWindows()[0].end;

    if (min_mz >= max_mz)
    {
      OPENMS_LOG_WARN << "No data to compress." << std::endl;
      return;
    }

    typedef std::vector<SimTypes::SimCoordinateType>::iterator GridTypeIt;
    std::vector<SimTypes::SimCoordinateType> grid;
    getSamplingGrid_(grid, min_mz, max_mz, 5); // every 5 Da we adjust the sampling width by local FWHM

    if (grid.size() < 3)
    {
      OPENMS_LOG_WARN << "Data spacing is weird - either you selected a very small interval or a very low resolution - or both. Not compressing." << std::endl;
      return;
    }

    Size point_count_before(0), point_count_after(0);
    SimTypes::SimPointType p;
    for (Size i = 0; i < experiment.size(); ++i)
    {
      if (experiment[i].size() <= 1)
        continue;

      if (experiment[i].isSorted() == false) // this should be true - however we check
      {
        experiment[i].sortByPosition();
      }

      // copy Spectrum and remove Peaks ..
      SimTypes::MSSimExperiment::SpectrumType cont = experiment[i];
      cont.clear(false);

      GridTypeIt grid_pos = grid.begin();
      GridTypeIt grid_pos_next(grid_pos + 1);

      double int_sum(0);
      bool break_scan(false);
      // match points to closest grid point
      for (Size j = 0; j < experiment[i].size(); ++j)
      {
        Size advance_by_binary_search = 3;
        while (fabs((*grid_pos_next) - experiment[i][j].getMZ()) < fabs((*grid_pos) - experiment[i][j].getMZ()))
        {
          if (int_sum > 0) // we collected some points before --> save them
          {
            p.setIntensity(int_sum);
            p.setMZ(*grid_pos);
            cont.push_back(p);
            int_sum = 0; // reset
          }

          if (--advance_by_binary_search == 0)
          {
            // advance using binary search
            grid_pos_next = std::lower_bound(grid_pos, grid.end(), experiment[i][j].getMZ());
            grid_pos = grid_pos_next - 1; // this should always work, since we ran at least 3 steps forward before
            advance_by_binary_search = 10; // just so we do not run into here again
          }
          else
          {
            // advance to next grid element
            ++grid_pos;
            ++grid_pos_next;
          }

          if (grid_pos_next == grid.end())
          {
            break_scan = true;
            break;
          }
        }
        if (break_scan)
          break; // skip remaining points of the scan (we reached the end of the grid)

        int_sum += experiment[i][j].getIntensity();

      } // end of scan

      if (int_sum > 0) // don't forget the last one
      {
        p.setIntensity(int_sum);
        p.setMZ(*grid_pos);
        cont.push_back(p);
      }

      point_count_before += experiment[i].size(); // stats
      experiment[i] = cont;
      point_count_after += experiment[i].size();
    }

    if (point_count_before != 0)
    {
      OPENMS_LOG_INFO << "Compressed data to grid ... " <<  point_count_before << " --> " << point_count_after << " (" << (point_count_after * 100 / point_count_before) << "%)\n";
    }
    else
    {
      OPENMS_LOG_INFO << "Not enough points in map .. did not compress!\n";
    }

    return;
  }

  SimTypes::SimIntensityType RawMSSignalSimulation::getFeatureScaledIntensity_(const SimTypes::SimIntensityType feature_intensity, const SimTypes::SimIntensityType natural_scaling_factor)
  {
    SimTypes::SimIntensityType intensity = feature_intensity * natural_scaling_factor * intensity_scale_;

    // add some noise
    // TODO: German comment
    // TODO: variables model f??r den intensit??ts-einfluss
    // e.g. sqrt(intensity) || ln(intensity)
    boost::normal_distribution<SimTypes::SimIntensityType> ndist(0, intensity_scale_stddev_ * intensity);
    intensity += ndist(rnd_gen_->getTechnicalRng());

    return intensity;
  }

}

ext: cpp

max_stars_repo_name: alexeigurevich/OpenMS

max_issues_repo_name: alexeigurevich/OpenMS

max_forks_repo_name: alexeigurevich/OpenMS


text:
#ifdef USE_EIGEN_BACKEND

/** Eigen3 backend.
 *
 * Only supports float32 computation with NHWC memory layout (at runtime and as input).
 */

// CR lpuchallafiore: Add multi-threading support (see "Evaluating with a Thread Pool" in the Eigen Tensor docs).

#include "../neuralnet/nninterface.h"

#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>
#include <zstr/src/zstr.hpp>

#include "../neuralnet/desc.h"
#include "../neuralnet/modelversion.h"
#include "../neuralnet/nninputs.h"
#include "../neuralnet/nneval.h"

using namespace std;
using Eigen::Tensor;
using Eigen::TensorMap;

//Eigen doesn't seem to have a way to make a const tensor map out of a const float* ??
//So we have to cast away qualifiers to build it.
#pragma GCC diagnostic ignored "-Wcast-qual"

// Eigen tensors are stored in column-major order, so an NHWC memory layout is given by Tensor<4>(C,W,H,N).

#define SCALAR float
#define TENSOR2 Tensor<SCALAR, 2>
#define TENSOR3 Tensor<SCALAR, 3>
#define TENSOR4 Tensor<SCALAR, 4>
#define TENSORMAP2 TensorMap<Tensor<SCALAR, 2>>
#define TENSORMAP3 TensorMap<Tensor<SCALAR, 3>>
#define TENSORMAP4 TensorMap<Tensor<SCALAR, 4>>

#define CONSTTENSOR2 const Tensor<SCALAR, 2>
#define CONSTTENSOR3 const Tensor<SCALAR, 3>
#define CONSTTENSOR4 const Tensor<SCALAR, 4>
#define CONSTTENSORMAP2 const TensorMap<Tensor<SCALAR, 2>>
#define CONSTTENSORMAP3 const TensorMap<Tensor<SCALAR, 3>>
#define CONSTTENSORMAP4 const TensorMap<Tensor<SCALAR, 4>>

// Debugging -----------------------------------------------------------------------------------------------------------
// #define DEBUG true

template <typename T>
void printTensorShape(const string& name, const T* t) {
  auto d = t->dimensions();
  cout << name << " rank=" << d.size() << " - (";
  for (int i = 0; i < d.size(); i++) {
    cout << d[i] << ",";
  }
  cout << ")" << endl;
}

#if DEBUG
#define DSHAPE(n, x) printTensorShape(n,x)
#define DTENSOR(n, x) cout << n << *x << endl
#else
#define DSHAPE(n, x)
#define DTENSOR(n, x)
#endif

// LoadedModel / ModelDesc ---------------------------------------------------------------------------------------------

struct LoadedModel {
  ModelDesc modelDesc;

  LoadedModel(const string& fileName) {
    ModelDesc::loadFromFileMaybeGZipped(fileName,modelDesc);
  }

  LoadedModel() = delete;
  LoadedModel(const LoadedModel&) = delete;
  LoadedModel& operator=(const LoadedModel&) = delete;
};

LoadedModel* NeuralNet::loadModelFile(const string& file) {
  LoadedModel* loadedModel = new LoadedModel(file);
  return loadedModel;
}

void NeuralNet::freeLoadedModel(LoadedModel* loadedModel) {
  delete loadedModel;
}

string NeuralNet::getModelName(const LoadedModel* loadedModel) {
  return loadedModel->modelDesc.name;
}

int NeuralNet::getModelVersion(const LoadedModel* loadedModel) {
  return loadedModel->modelDesc.version;
}

Rules NeuralNet::getSupportedRules(const LoadedModel* loadedModel, const Rules& desiredRules, bool& supported) {
  return loadedModel->modelDesc.getSupportedRules(desiredRules, supported);
}

//------------------------------------------------------------------------------

struct ComputeContext {
  int nnXLen;
  int nnYLen;
};

ComputeContext* NeuralNet::createComputeContext(
  const std::vector<int>& gpuIdxs,
  Logger* logger,
  int nnXLen,
  int nnYLen,
  const string& openCLTunerFile,
  const string& homeDataDirOverride,
  bool openCLReTunePerBoardSize,
  enabled_t useFP16Mode,
  enabled_t useNHWCMode,
  const LoadedModel* loadedModel
) {
  (void)gpuIdxs;
  (void)logger;
  (void)openCLTunerFile;
  (void)homeDataDirOverride;
  (void)openCLReTunePerBoardSize;
  (void)loadedModel;

  bool useFP16 = useFP16Mode == enabled_t::True ? true : false;
  bool useNHWC = useNHWCMode == enabled_t::False ? false : true;

  if(useFP16)
    throw StringError("Eigen backend: useFP16 = true not supported");
  if(!useNHWC)
    throw StringError("Eigen backend: useNHWC = false not supported");

  ComputeContext* context = new ComputeContext();
  context->nnXLen = nnXLen;
  context->nnYLen = nnYLen;
  return context;
}

void NeuralNet::freeComputeContext(ComputeContext* computeContext) {
  delete computeContext;
}

// Helpers --------------------------------------------------------------------------------------------------------------

static void computeMaskSum(CONSTTENSORMAP3* mask, float* maskSum) {
  for (int n = 0; n < mask->dimension(2); n++) {
    float s = 0.f;
    for (int h = 0; h < mask->dimension(1); h++) {
      for (int w = 0; w < mask->dimension(0); w++) {
        s += (*mask)(w, h, n);
      }
    }
    maskSum[n] = s;
  }
}

// in NxHxWxC, bias NxC
static void addNCBiasInplace(TENSORMAP4* in, CONSTTENSORMAP2* bias) {
  assert(in->dimension(0) == bias->dimension(0) && in->dimension(3) == bias->dimension(1));
  for (int n = 0; n < in->dimension(3); n++) {
    for (int h = 0; h < in->dimension(2); h++) {
      for (int w = 0; w < in->dimension(1); w++) {
        for (int c = 0; c < in->dimension(0); c++) {
          (*in)(c,w,h,n) += (*bias)(c,n);
        }
      }
    }
  }
}

static void poolRowsGPool(CONSTTENSORMAP4* in, TENSORMAP2* out, const float* maskSum) {
  for (int n = 0; n < in->dimension(3); n++) {
    for (int c = 0; c < in->dimension(0); c++) {
      float s = 0.f;
      float m = 0.f;
      for (int h = 0; h < in->dimension(2); h++) {
        for (int w = 0; w < in->dimension(1); w++) {
          float x = (*in)(c, w, h, n);
          s += x;
          m = max(m, x);
        }
      }
      float div = maskSum[n];
      float sqrtdiv = sqrt(div);
      float mean = s / div;
      (*out)(c, n) = mean;
      (*out)(c + in->dimension(0), n) = mean * (sqrtdiv - 14.f) * 0.1f;
      (*out)(c + 2*in->dimension(0), n) = m;
    }
  }
}

// // Given input [n,w,h,c] fills output of shape [n,c] with sum over c.
// static void poolRowsSum(CONSTTENSORMAP4* in, TENSORMAP2* out, float scaleSum) {
//   for (int n = 0; n < in->dimension(3); n++) {
//     for (int c = 0; c < in->dimension(0); c++) {
//       float s = 0.f;
//       for (int h = 0; h < in->dimension(2); h++) {
//         for (int w = 0; w < in->dimension(1); w++) {
//           float x = in(c, w, h, n);
//           s += x;
//         }
//       }
//       out(c, n) = s * scaleSum;
//     }
//   }
// }

static void poolRowsValueHead(CONSTTENSORMAP4* in, TENSORMAP2* out, const float* maskSum) {
  for (int n = 0; n < in->dimension(3); n++) {
    for (int c = 0; c < in->dimension(0); c++) {
      float s = 0.f;
      for (int h = 0; h < in->dimension(2); h++) {
        for (int w = 0; w < in->dimension(1); w++) {
          float x = (*in)(c, w, h, n);
          s += x;
        }
      }
      float div = maskSum[n];
      float sqrtdiv = sqrt(div);
      float mean = s / div;
      (*out)(c, n) = mean;
      (*out)(c + in->dimension(0), n) = mean * (sqrtdiv - 14.f) * 0.1f;
      (*out)(c + 2*in->dimension(0), n) = mean * ((sqrtdiv - 14.0f) * (sqrtdiv - 14.0f) * 0.01f - 0.1f);
    }
  }
}

// Layers --------------------------------------------------------------------------------------------------------------

// Convolution layer with zero-padding.
struct ConvLayer {
  string name;

  Eigen::array<pair<int, int>, 4> paddings;
  TENSOR4 kernel;
  int inChannels, outChannels;

  ConvLayer() = delete;
  ConvLayer(const ConvLayer&) = delete;
  ConvLayer& operator=(const ConvLayer&) = delete;

  ConvLayer(const ConvLayerDesc& desc) {
    name = desc.name;
    int convYSize = desc.convYSize;
    int convXSize = desc.convXSize;
    inChannels = desc.inChannels;
    outChannels = desc.outChannels;
    //Currently eigen impl doesn't support dilated convs
    int dilationY = desc.dilationY;
    int dilationX = desc.dilationX;
    int paddingX = (convXSize / 2) * dilationX;
    int paddingY = (convYSize / 2) * dilationY;

    if(dilationX != 1 || dilationY != 1)
      throw StringError("Eigen backend: Encountered convolution dilation factors other than 1, not supported");

    assert(convXSize % 2 == 1);
    assert(convYSize % 2 == 1);

    paddings[0] = make_pair(0, 0);                // C
    paddings[1] = make_pair(paddingX, paddingX);  // W
    paddings[2] = make_pair(paddingY, paddingY);  // H
    paddings[3] = make_pair(0, 0);                // N

    // CR-someday lpuchallafiore: optimize NHWC vs NCHW, etc.
    kernel = TensorMap<const Tensor<const SCALAR, 4>>(
      desc.weights.data(), convXSize, convYSize, inChannels, outChannels);
  }

  void apply(CONSTTENSORMAP4* input, TENSORMAP4* output, bool accumulate) const {
    auto padded = input->pad(paddings);
    assert(output->dimension(0) == outChannels);
    for(int n = 0; n < input->dimension(3); n++) {
      auto inN = padded.chip(n, 3);
      for(int oc = 0; oc < outChannels; oc++) {
        TENSOR2 sum(input->dimension(1), input->dimension(2));
        sum.setZero();

        for(int ic = 0; ic < inChannels; ic++) {
          Eigen::array<ptrdiff_t, 2> dims({0, 1});
          auto kChip = kernel.chip(oc, 3).chip(ic, 2);
          auto inNC = inN.chip(ic, 0);
          sum += inNC.convolve(kChip, dims);
        }

        if(accumulate)
          output->chip(n, 3).chip(oc, 0) += sum;
        else
          output->chip(n, 3).chip(oc, 0) = sum;
      }
    }
  }
};

//--------------------------------------------------------------

struct BatchNormLayer {
  string name;

  vector<float> mergedScale;
  vector<float> mergedBias;

  BatchNormLayer() = delete;
  BatchNormLayer(const BatchNormLayer&) = delete;
  BatchNormLayer& operator=(const BatchNormLayer&) = delete;

  BatchNormLayer(const BatchNormLayerDesc& desc) {
    name = desc.name;
    int numChannels = desc.numChannels;
    float epsilon = desc.epsilon;

    mergedScale.resize(numChannels);
    mergedBias.resize(numChannels);
    for(int c = 0; c < numChannels; c++) {
      mergedScale[c] = desc.scale[c] / sqrt(desc.variance[c] + epsilon);
      mergedBias[c] = desc.bias[c] - mergedScale[c] * desc.mean[c];
    }
  }

  // Mask should be in 'NHW' format (no "C" channel).
  void apply(
    bool applyRelu,
    CONSTTENSORMAP4* input,
    TENSORMAP4* output,
    CONSTTENSORMAP3* mask
  ) const {
    for(int c = 0; c < input->dimension(0); c++) {
      auto inC = input->chip(c, 0);
      auto x = inC * mergedScale[c] + mergedBias[c];
      auto z = TENSOR3(mask->dimension(0), mask->dimension(1), mask->dimension(2)).setZero();
      if(applyRelu)
        output->chip(c, 0) = (*mask == 1.f).select(x.cwiseMax(0.f), z);
      else
        output->chip(c, 0) = (*mask == 1.f).select(x, z);
    }
  }
};

//--------------------------------------------------------------

struct ActivationLayer {
  string name;

  ActivationLayer() = delete;
  ActivationLayer(const ActivationLayer&) = delete;
  ActivationLayer& operator=(const ActivationLayer&) = delete;

  ActivationLayer(const ActivationLayerDesc& desc) { name = desc.name; }

  template <int N>
  void apply(const Tensor<SCALAR, N>* input, Tensor<SCALAR, N>* output) const { *output = input->cwiseMax(0.f); }
  template <int N>
  void apply(const TensorMap<Tensor<SCALAR, N>>* input, TensorMap<Tensor<SCALAR, N>>* output) const { *output = input->cwiseMax(0.f); }
};

//--------------------------------------------------------------

struct MatMulLayer {
  string name;
  TENSOR2 weights;

  MatMulLayer() = delete;
  MatMulLayer(const MatMulLayer&) = delete;
  MatMulLayer& operator=(const MatMulLayer&) = delete;

  MatMulLayer(const MatMulLayerDesc& desc)
    : name(desc.name)
  {
    weights = TENSOR2(desc.outChannels, desc.inChannels);
    memcpy(weights.data(), desc.weights.data(), sizeof(SCALAR) * weights.size());
  }

  void apply(CONSTTENSORMAP2* in, TENSORMAP2* out) const {
    Eigen::array<Eigen::IndexPair<int>, 1> product_dims = { Eigen::IndexPair<int>(1, 0) };
    *out = weights.contract(*in, product_dims);
  }
};

struct MatBiasLayer {
  string name;
  std::vector<float> weights;

  MatBiasLayer() = delete;
  MatBiasLayer(const MatBiasLayer&) = delete;
  MatBiasLayer& operator=(const MatBiasLayer&) = delete;

  MatBiasLayer(const MatBiasLayerDesc& desc)
    : name(desc.name),
      weights(desc.weights) {}

  void apply(TENSORMAP2* mat) const {
    for(int n = 0; n < mat->dimension(1); n++) {
      for(int c = 0; c < mat->dimension(0); c++) {
        (*mat)(c, n) += weights[c];
      }
    }
  }
};

// Blocks
// --------------------------------------------------------------------------------------------------------------

struct ResidualBlockIntf {
  virtual ~ResidualBlockIntf(){}

  virtual void apply(
    TENSORMAP4* trunk,
    TENSORMAP4* trunkScratch,
    TENSORMAP4* regularOut,
    TENSORMAP4* regularScratch,
    TENSORMAP4* midIn,
    TENSORMAP4* midScratch,
    TENSORMAP4* gpoolOut,
    TENSORMAP4* gpoolOut2,
    TENSORMAP2* gpoolConcat,
    TENSORMAP2* gpoolBias,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const = 0;
};

struct ResidualBlock final : public ResidualBlockIntf {
  string name;
  BatchNormLayer preBN;
  ConvLayer regularConv;
  BatchNormLayer midBN;
  ConvLayer finalConv;

  ResidualBlock() = delete;
  ResidualBlock(const ResidualBlock&) = delete;
  ResidualBlock& operator=(const ResidualBlock&) = delete;

  ~ResidualBlock(){}

  ResidualBlock(const ResidualBlockDesc& desc)
    : name(desc.name),
      preBN(desc.preBN),
      regularConv(desc.regularConv),
      midBN(desc.midBN),
      finalConv(desc.finalConv) {}

  void apply(
    TENSORMAP4* trunk,
    TENSORMAP4* trunkScratch,
    TENSORMAP4* regularOut,
    TENSORMAP4* regularScratch,
    TENSORMAP4* midIn,
    TENSORMAP4* midScratch,
    TENSORMAP4* gpoolOut,
    TENSORMAP4* gpoolOut2,
    TENSORMAP2* gpoolConcat,
    TENSORMAP2* gpoolBias,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const override {
    (void)regularOut;
    (void)regularScratch;
    (void)gpoolOut;
    (void)gpoolOut2;
    (void)gpoolConcat;
    (void)gpoolBias;
    (void)maskSum;
    const bool applyBNRelu = true;
    preBN.apply(applyBNRelu, trunk, trunkScratch, mask);
    regularConv.apply(trunkScratch, midIn, false);
    midBN.apply(applyBNRelu, midIn, midScratch, mask);
    finalConv.apply(midScratch, trunk, true);
  }
};

// // Given two tensors with shapes inA: [n, h, w, cA] and inB: [n, h, w, cB]
// // Copy them into a single tensor out: [n, h, w, cA + cB]
// TENSOR4 concatTensors(CONSTTENSOR4& a, CONSTTENSOR4& b) {
//   assert(a->dimension(1) == b->dimension(1) && a->dimension(2) == b->dimension(2) && a->dimension(3) == b->dimension(3));
//   TENSOR4 x = TENSOR4(/* C */ a->dimension(0) + b->dimension(0),
//                                           /* W */ a->dimension(1),
//                                           /* H */ a->dimension(2),
//                                           /* N */ a->dimension(3));
//   for (int n = 0; n < a->dimension(3); n++) {
//     for (int h = 0; h < a->dimension(2); h++) {
//       for (int w = 0; w < a->dimension(1); w++) {
//         int c = 0;
//         for (int ca = 0; a->dimension(0); ca++, c++) {
//           x(c,w,h,n) = a(ca,w,h,n);
//         }
//         for (int cb = 0; b->dimension(0); cb++, c++) {
//           x(c,w,h,n) = b(cb,w,h,n);
//         }
//       }
//     }
//   }
//   return x;
// }

struct GlobalPoolingResidualBlock final : public ResidualBlockIntf {
  string name;
  BatchNormLayer preBN;
  ActivationLayer preActivation;
  ConvLayer regularConv;
  ConvLayer gpoolConv;
  BatchNormLayer gpoolBN;
  ActivationLayer gpoolActivation;
  MatMulLayer gpoolToBiasMul;
  BatchNormLayer midBN;
  ActivationLayer midActivation;
  ConvLayer finalConv;

  GlobalPoolingResidualBlock() = delete;
  GlobalPoolingResidualBlock(const GlobalPoolingResidualBlock&) = delete;
  GlobalPoolingResidualBlock& operator=(const GlobalPoolingResidualBlock&) = delete;

  ~GlobalPoolingResidualBlock(){}

  GlobalPoolingResidualBlock(const GlobalPoolingResidualBlockDesc& desc)
    : name(desc.name),
      preBN(desc.preBN),
      preActivation(desc.preActivation),
      regularConv(desc.regularConv),
      gpoolConv(desc.gpoolConv),
      gpoolBN(desc.gpoolBN),
      gpoolActivation(desc.gpoolActivation),
      gpoolToBiasMul(desc.gpoolToBiasMul),
      midBN(desc.midBN),
      midActivation(desc.midActivation),
      finalConv(desc.finalConv) {}

  void apply(
    TENSORMAP4* trunk,
    TENSORMAP4* trunkScratch,
    TENSORMAP4* regularOut,
    TENSORMAP4* regularScratch,
    TENSORMAP4* midIn,
    TENSORMAP4* midScratch,
    TENSORMAP4* gpoolOut,
    TENSORMAP4* gpoolOut2,
    TENSORMAP2* gpoolConcat,
    TENSORMAP2* gpoolBias,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const override {
    (void)midIn;
    (void)midScratch;
    const bool applyBNRelu = true;
    DTENSOR("trunk", trunk);
    DTENSOR("mask", mask);
    preBN.apply(applyBNRelu, trunk, trunkScratch, mask);
    DTENSOR("trunkScratch", trunkScratch);
    regularConv.apply(trunkScratch, regularOut, false);
    DTENSOR("regularOut", regularOut);
    gpoolConv.apply(trunkScratch, gpoolOut, false);
    DTENSOR("gpoolOut", gpoolOut);
    gpoolBN.apply(applyBNRelu, gpoolOut, gpoolOut2, mask);
    DTENSOR("gpoolOut2", gpoolOut2);
    poolRowsGPool(gpoolOut2, gpoolConcat, maskSum);
    gpoolToBiasMul.apply(gpoolConcat, gpoolBias);
    addNCBiasInplace(regularOut, gpoolBias);
    midBN.apply(applyBNRelu, regularOut, regularScratch, mask);
    finalConv.apply(regularScratch, trunk, true);
    DSHAPE("trunk", trunk);
    DSHAPE("trunkScratch", trunkScratch);
    DSHAPE("regularOut", regularOut);
    DSHAPE("gpoolOut", gpoolOut);
    DSHAPE("gpoolOut2", gpoolOut2);
    DSHAPE("gpoolConcat", gpoolConcat);
    DSHAPE("gpoolBias", gpoolBias);
    DSHAPE("mask", mask);
  }
};

struct Trunk {
  string name;
  int version;
  int numBlocks;

  ConvLayer initialConv;
  MatMulLayer initialMatMul;
  vector<pair<int, ResidualBlockIntf*>> blocks;
  BatchNormLayer trunkTipBN;
  ActivationLayer trunkTipActivation;

  Trunk() = delete;
  Trunk(const Trunk&) = delete;
  Trunk& operator=(const Trunk&) = delete;

  Trunk(const TrunkDesc& desc)
    : name(desc.name),
      version(desc.version),
      numBlocks(desc.numBlocks),
      initialConv(desc.initialConv),
      initialMatMul(desc.initialMatMul),
      trunkTipBN(desc.trunkTipBN),
      trunkTipActivation(desc.trunkTipActivation)
  {
    for (int i = 0; i < numBlocks; ++i) {
      if (desc.blocks[i].first == ORDINARY_BLOCK_KIND) {
        ResidualBlockDesc* blockDesc = (ResidualBlockDesc*)desc.blocks[i].second;
        ResidualBlockIntf* block = new ResidualBlock(*blockDesc);
        blocks.push_back(make_pair(ORDINARY_BLOCK_KIND, block));
      }
      else if (desc.blocks[i].first == DILATED_BLOCK_KIND) {
        throw StringError("Eigen backend: Dilated residual blocks are not supported right now");
      }
      else if (desc.blocks[i].first == GLOBAL_POOLING_BLOCK_KIND) {
        GlobalPoolingResidualBlockDesc* blockDesc = (GlobalPoolingResidualBlockDesc*)desc.blocks[i].second;
        GlobalPoolingResidualBlock* block = new GlobalPoolingResidualBlock(*blockDesc);
        blocks.push_back(make_pair(GLOBAL_POOLING_BLOCK_KIND, block));
      }
      else {
        ASSERT_UNREACHABLE;
      }
    }
  }

  virtual ~Trunk() {
    for (auto p : blocks) {
      delete p.second;
    }
  }

  void apply(
    CONSTTENSORMAP4* input,
    CONSTTENSORMAP2* inputGlobal,
    TENSORMAP2* inputMatMulOut,
    TENSORMAP4* trunk,
    TENSORMAP4* trunkScratch,
    TENSORMAP4* regularOut,
    TENSORMAP4* regularScratch,
    TENSORMAP4* midIn,
    TENSORMAP4* midScratch,
    TENSORMAP4* gpoolOut,
    TENSORMAP4* gpoolOut2,
    TENSORMAP2* gpoolConcat,
    TENSORMAP2* gpoolBias,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const {

    initialConv.apply(input, trunkScratch, false);
    initialMatMul.apply(inputGlobal, inputMatMulOut);
    addNCBiasInplace(trunkScratch, inputMatMulOut);

    // apply blocks
    // Flip trunkBuf and trunkScratchBuf so that the result gets accumulated in trunkScratchBuf
    for (auto block : blocks) {
      block.second->apply(
        trunkScratch,
        trunk,
        regularOut,
        regularScratch,
        midIn,
        midScratch,
        gpoolOut,
        gpoolOut2,
        gpoolConcat,
        gpoolBias,
        mask,
        maskSum
      );
    }

    // And now with the final BN port it from trunkScratchBuf to trunkBuf.
    const bool applyBNRelu = true;
    trunkTipBN.apply(applyBNRelu, trunkScratch, trunk, mask);
  }
};

struct PolicyHead {
  string name;
  int version;

  ConvLayer p1Conv;
  ConvLayer g1Conv;
  BatchNormLayer g1BN;
  ActivationLayer g1Activation;
  MatMulLayer gpoolToBiasMul;
  BatchNormLayer p1BN;
  ActivationLayer p1Activation;
  ConvLayer p2Conv;
  MatMulLayer gpoolToPassMul;

  PolicyHead() = delete;
  PolicyHead(const PolicyHead&) = delete;
  PolicyHead& operator=(const PolicyHead&) = delete;

  PolicyHead(const PolicyHeadDesc& desc)
    : name(desc.name),
      version(desc.version),
      p1Conv(desc.p1Conv),
      g1Conv(desc.g1Conv),
      g1BN(desc.g1BN),
      g1Activation(desc.g1Activation),
      gpoolToBiasMul(desc.gpoolToBiasMul),
      p1BN(desc.p1BN),
      p1Activation(desc.p1Activation),
      p2Conv(desc.p2Conv),
      gpoolToPassMul(desc.gpoolToPassMul) {}

  void apply(
    CONSTTENSORMAP4* trunk,
    TENSORMAP4* p1Out,
    TENSORMAP4* p1Out2,
    TENSORMAP4* g1Out,
    TENSORMAP4* g1Out2,
    TENSORMAP2* g1Concat,
    TENSORMAP2* g1Bias,
    TENSORMAP2* policyPass,
    TENSORMAP4* policy,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const {
    const bool applyBNRelu = true;
    p1Conv.apply(trunk, p1Out, false);
    g1Conv.apply(trunk, g1Out, false);
    g1BN.apply(applyBNRelu, g1Out, g1Out2, mask);
    poolRowsGPool(g1Out2, g1Concat, maskSum);
    gpoolToBiasMul.apply(g1Concat, g1Bias);
    addNCBiasInplace(p1Out, g1Bias);
    p1BN.apply(true, p1Out, p1Out2, mask);
    p2Conv.apply(p1Out2, policy, false);
    gpoolToPassMul.apply(g1Concat, policyPass);
  }
};

struct ValueHead {
  string name;
  int version;

  ConvLayer v1Conv;
  BatchNormLayer v1BN;
  ActivationLayer v1Activation;
  MatMulLayer v2Mul;
  MatBiasLayer v2Bias;
  ActivationLayer v2Activation;
  MatMulLayer v3Mul;
  MatBiasLayer v3Bias;
  MatMulLayer sv3Mul;
  MatBiasLayer sv3Bias;
  ConvLayer vOwnershipConv;

  ValueHead() = delete;
  ValueHead(const ValueHead&) = delete;
  ValueHead& operator=(const ValueHead&) = delete;

  ValueHead(const ValueHeadDesc& desc)
    : name(desc.name),
      version(desc.version),
      v1Conv(desc.v1Conv),
      v1BN(desc.v1BN),
      v1Activation(desc.v1Activation),
      v2Mul(desc.v2Mul),
      v2Bias(desc.v2Bias),
      v2Activation(desc.v2Activation),
      v3Mul(desc.v3Mul),
      v3Bias(desc.v3Bias),
      sv3Mul(desc.sv3Mul),
      sv3Bias(desc.sv3Bias),
      vOwnershipConv(desc.vOwnershipConv) {}

  void apply(
    CONSTTENSORMAP4* trunk,
    TENSORMAP4* v1Out,
    TENSORMAP4* v1Out2,
    TENSORMAP2* v1Mean,
    TENSORMAP2* v2Out,
    TENSORMAP2* value,
    TENSORMAP2* scoreValue,
    TENSORMAP4* ownership,
    CONSTTENSORMAP3* mask,
    const float* maskSum
  ) const {
    bool applyBNRelu = true;
    v1Conv.apply(trunk, v1Out, false);
    v1BN.apply(applyBNRelu, v1Out, v1Out2, mask);
    poolRowsValueHead(v1Out2, v1Mean, maskSum);
    v2Mul.apply(v1Mean, v2Out);
    v2Bias.apply(v2Out);
    v2Activation.apply(v2Out, v2Out);
    v3Mul.apply(v2Out, value);
    v3Bias.apply(value);

    sv3Mul.apply(v2Out, scoreValue);
    sv3Bias.apply(scoreValue);

    vOwnershipConv.apply(v1Out2, ownership, false);
  }
};

// Model and Buffer I/O ------------------------------------------------------------------------------------------------

struct Model {
  string name;
  int version;
  int numInputChannels;
  int numInputGlobalChannels;
  int numValueChannels;
  int numScoreValueChannels;
  int numOwnershipChannels;

  Trunk trunk;
  PolicyHead policyHead;
  ValueHead valueHead;

  Model() = delete;
  Model(const Model&) = delete;
  Model& operator=(const Model&) = delete;

  Model(const ModelDesc& desc)
    : name(desc.name), version(desc.version), numInputChannels(desc.numInputChannels),
      numInputGlobalChannels(desc.numInputGlobalChannels),
      numValueChannels(desc.numValueChannels),
      numScoreValueChannels(desc.numScoreValueChannels),
      numOwnershipChannels(desc.numOwnershipChannels),
      trunk(desc.trunk),
      policyHead(desc.policyHead),
      valueHead(desc.valueHead) {}

  void apply(
    CONSTTENSORMAP4* input,
    CONSTTENSORMAP2* inputGlobal,
    TENSORMAP2* inputMatMulOut,
    TENSORMAP4* trunkBuf,
    TENSORMAP4* trunkScratch,
    TENSORMAP4* regularOut,
    TENSORMAP4* regularScratch,
    TENSORMAP4* midIn,
    TENSORMAP4* midScratch,
    TENSORMAP4* gpoolOut,
    TENSORMAP4* gpoolOut2,
    TENSORMAP2* gpoolConcat,
    TENSORMAP2* gpoolBias,

    TENSORMAP4* p1Out,
    TENSORMAP4* p1Out2,
    TENSORMAP4* g1Out,
    TENSORMAP4* g1Out2,
    TENSORMAP2* g1Concat,
    TENSORMAP2* g1Bias,
    TENSORMAP2* policyPass,
    TENSORMAP4* policy,

    TENSORMAP4* v1Out,
    TENSORMAP4* v1Out2,
    TENSORMAP2* v1Mean,
    TENSORMAP2* v2Out,
    TENSORMAP2* value,
    TENSORMAP2* scoreValue,
    TENSORMAP4* ownership,

    TENSORMAP3* mask,
    float* maskSum
  ) const {
    *mask = input->chip(0,0);
    computeMaskSum(mask,maskSum);

    trunk.apply(
      input,
      inputGlobal,
      inputMatMulOut,
      trunkBuf,
      trunkScratch,
      regularOut,
      regularScratch,
      midIn,
      midScratch,
      gpoolOut,
      gpoolOut2,
      gpoolConcat,
      gpoolBias,
      mask,
      maskSum
    );
    policyHead.apply(
      trunkBuf,
      p1Out,
      p1Out2,
      g1Out,
      g1Out2,
      g1Concat,
      g1Bias,
      policyPass,
      policy,
      mask,
      maskSum
    );
    valueHead.apply(
      trunkBuf,
      v1Out,
      v1Out2,
      v1Mean,
      v2Out,
      value,
      scoreValue,
      ownership,
      mask,
      maskSum
    );
  }
};

//--------------------------------------------------------------

struct Buffers {
  TENSOR2 inputMatMulOut;
  TENSOR4 trunk;
  TENSOR4 trunkScratch;
  TENSOR4 regularOut;
  TENSOR4 regularScratch;
  TENSOR4 midIn;
  TENSOR4 midScratch;
  TENSOR4 gpoolOut;
  TENSOR4 gpoolOut2;
  TENSOR2 gpoolConcat;
  TENSOR2 gpoolBias;

  TENSOR4 p1Out;
  TENSOR4 p1Out2;
  TENSOR4 g1Out;
  TENSOR4 g1Out2;
  TENSOR2 g1Concat;
  TENSOR2 g1Bias;
  TENSOR2 policyPass;
  TENSOR4 policy;

  TENSOR4 v1Out;
  TENSOR4 v1Out2;
  TENSOR2 v1Mean;
  TENSOR2 v2Out;
  TENSOR2 value;
  TENSOR2 scoreValue;
  TENSOR4 ownership;

  TENSOR3 mask;
  vector<float> maskSum;

  Buffers(
    const ModelDesc& desc,
    int maxBatchSize,
    int nnXLen,
    int nnYLen
  ) :
    inputMatMulOut(desc.trunk.trunkNumChannels, maxBatchSize),
    trunk(desc.trunk.trunkNumChannels, nnXLen, nnYLen, maxBatchSize),
    trunkScratch(desc.trunk.trunkNumChannels, nnXLen, nnYLen, maxBatchSize),
    regularOut(desc.trunk.regularNumChannels, nnXLen, nnYLen, maxBatchSize),
    regularScratch(desc.trunk.regularNumChannels, nnXLen, nnYLen, maxBatchSize),
    midIn(desc.trunk.midNumChannels, nnXLen, nnYLen, maxBatchSize),
    midScratch(desc.trunk.midNumChannels, nnXLen, nnYLen, maxBatchSize),
    gpoolOut(desc.trunk.gpoolNumChannels, nnXLen, nnYLen, maxBatchSize),
    gpoolOut2(desc.trunk.gpoolNumChannels, nnXLen, nnYLen, maxBatchSize),
    gpoolConcat(desc.trunk.gpoolNumChannels*3, maxBatchSize),
    gpoolBias(desc.trunk.regularNumChannels, maxBatchSize),

    p1Out(desc.policyHead.p1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    p1Out2(desc.policyHead.p1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    g1Out(desc.policyHead.g1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    g1Out2(desc.policyHead.g1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    g1Concat(desc.policyHead.g1Conv.outChannels*3, maxBatchSize),
    g1Bias(desc.policyHead.gpoolToBiasMul.outChannels, maxBatchSize),
    policyPass(desc.policyHead.gpoolToPassMul.outChannels, maxBatchSize),
    policy(desc.policyHead.p2Conv.outChannels, nnXLen, nnYLen, maxBatchSize),

    v1Out(desc.valueHead.v1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    v1Out2(desc.valueHead.v1Conv.outChannels, nnXLen, nnYLen, maxBatchSize),
    v1Mean(desc.valueHead.v1Conv.outChannels*3, maxBatchSize),
    v2Out(desc.valueHead.v2Mul.outChannels, maxBatchSize),
    value(desc.valueHead.v3Mul.outChannels, maxBatchSize),
    scoreValue(desc.valueHead.sv3Mul.outChannels, maxBatchSize),
    ownership(desc.valueHead.vOwnershipConv.outChannels, nnXLen, nnYLen, maxBatchSize),

    mask(nnXLen, nnYLen, maxBatchSize),
    maskSum(maxBatchSize)
  {}
};

//--------------------------------------------------------------

struct InputBuffers {
  int maxBatchSize;

  size_t singleInputElts;
  size_t singleInputGlobalElts;

  size_t singlePolicyPassResultElts;
  size_t singlePolicyResultElts;
  size_t singleValueResultElts;
  size_t singleScoreValueResultElts;
  size_t singleOwnershipResultElts;

  std::vector<float> spatialInput;
  std::vector<float> globalInput;

  InputBuffers(const LoadedModel* loadedModel, int maxBatchSz, int nnXLen, int nnYLen) {
    const ModelDesc& m = loadedModel->modelDesc;

    int xSize = nnXLen;
    int ySize = nnYLen;

    maxBatchSize = maxBatchSz;
    singleInputElts = m.numInputChannels * xSize * ySize;
    singleInputGlobalElts = m.numInputGlobalChannels;

    singlePolicyPassResultElts = (size_t)(1);
    singlePolicyResultElts = (size_t)(xSize * ySize);
    singleValueResultElts = (size_t)m.numValueChannels;
    singleScoreValueResultElts = (size_t)m.numScoreValueChannels;
    singleOwnershipResultElts = (size_t)m.numOwnershipChannels * xSize * ySize;

    assert(NNModelVersion::getNumSpatialFeatures(m.version) == m.numInputChannels);
    assert(NNModelVersion::getNumGlobalFeatures(m.version) == m.numInputGlobalChannels);

    spatialInput = vector<float>(m.numInputChannels * xSize * ySize * maxBatchSize);
    globalInput = vector<float>(m.numInputGlobalChannels * maxBatchSize);
  }

  ~InputBuffers() { }

  InputBuffers() = delete;
  InputBuffers(const InputBuffers&) = delete;
  InputBuffers& operator=(const InputBuffers&) = delete;
};

InputBuffers* NeuralNet::createInputBuffers(const LoadedModel* loadedModel, int maxBatchSize, int nnXLen, int nnYLen) {
  return new InputBuffers(loadedModel, maxBatchSize, nnXLen, nnYLen);
}
void NeuralNet::freeInputBuffers(InputBuffers* inputBuffers) {
  delete inputBuffers;
}

// float* NeuralNet::getBatchEltSpatialInplace(InputBuffers* inputBuffers, int nIdx) {
//   assert(nIdx < inputBuffers->maxBatchSize);
//   return inputBuffers->spatialInput.data() + (inputBuffers->singleInputElts * nIdx);
// }

// float* NeuralNet::getBatchEltGlobalInplace(InputBuffers* inputBuffers, int rowIdx) {
//   assert(rowIdx < inputBuffers->maxBatchSize);
//   return inputBuffers->globalInput.data() + (inputBuffers->singleInputGlobalElts * rowIdx);
// }

// int NeuralNet::getBatchEltSpatialLen(const InputBuffers* inputBuffers) {
//   return inputBuffers->singleInputElts;
// }
// int NeuralNet::getBatchEltGlobalLen(const InputBuffers* inputBuffers) {
//   return inputBuffers->singleInputGlobalElts;
// }

// bool* NeuralNet::getSymmetriesInplace(InputBuffers* inputBuffers) {
//   return inputBuffers->symmetriesBuffer;
// }

// NeuralNet -----------------------------------------------------------------------------------------------------------

void NeuralNet::globalInitialize() {
  // no-op for cpu
}

void NeuralNet::globalCleanup() {
  // no-op for cpu
}

struct ComputeHandle {
  const ComputeContext* context;
  int maxBatchSize;
  bool inputsUseNHWC;
  Model model;
  Buffers buffers;

  ComputeHandle() = delete;
  ComputeHandle(const ComputeHandle&) = delete;
  ComputeHandle& operator=(const ComputeHandle&) = delete;

  ComputeHandle(const ComputeContext* ctx, const LoadedModel& loadedModel, int maxBSize, bool iNHWC)
    : context(ctx),
      maxBatchSize(maxBSize),
      inputsUseNHWC(iNHWC),
      model(loadedModel.modelDesc),
      buffers(loadedModel.modelDesc,maxBSize,ctx->nnXLen,ctx->nnYLen)
  {}
};

ComputeHandle* NeuralNet::createComputeHandle(
  ComputeContext* context,
  const LoadedModel* loadedModel,
  Logger* logger,
  int maxBatchSize,
  bool requireExactNNLen,
  bool inputsUseNHWC,
  int gpuIdxForThisThread
) {
  if(logger != NULL) {
    logger->write("Eigen (CPU) backend: Model version " + Global::intToString(loadedModel->modelDesc.version));
    logger->write("Eigen (CPU) backend: Model name: " + loadedModel->modelDesc.name);
  }

  (void)requireExactNNLen; //We don't bother with mask optimizations if we know exact sizes right now.
  (void)gpuIdxForThisThread; //Doesn't matter

  if(!inputsUseNHWC)
    throw StringError("Eigen backend: inputsUseNHWC = false unsupported");
  return new ComputeHandle(context, *loadedModel, maxBatchSize, inputsUseNHWC);
}

void NeuralNet::freeComputeHandle(ComputeHandle* gpuHandle) {
  delete gpuHandle;
}

void NeuralNet::getOutput(
  ComputeHandle* computeHandle,
  InputBuffers* inputBuffers,
  int numBatchEltsFilled,
  NNResultBuf** inputBufs,
  int symmetry,
  vector<NNOutput*>& outputs
) {
  assert(numBatchEltsFilled <= inputBuffers->maxBatchSize);
  assert(numBatchEltsFilled > 0);
  int batchSize = numBatchEltsFilled;
  int nnXLen = computeHandle->context->nnXLen;
  int nnYLen = computeHandle->context->nnYLen;
  int version = computeHandle->model.version;

  int numSpatialFeatures = NNModelVersion::getNumSpatialFeatures(version);
  int numGlobalFeatures = NNModelVersion::getNumGlobalFeatures(version);
  assert(numSpatialFeatures == computeHandle->model.numInputChannels);
  assert(numSpatialFeatures * nnXLen * nnYLen == inputBuffers->singleInputElts);
  assert(numGlobalFeatures == inputBuffers->singleInputGlobalElts);

  for(int nIdx = 0; nIdx<batchSize; nIdx++) {
    float* rowSpatialInput = inputBuffers->spatialInput.data() + (inputBuffers->singleInputElts * nIdx);
    float* rowGlobalInput = inputBuffers->globalInput.data() + (inputBuffers->singleInputGlobalElts * nIdx);

    const float* rowGlobal = inputBufs[nIdx]->rowGlobal;
    const float* rowSpatial = inputBufs[nIdx]->rowSpatial;
    std::copy(rowGlobal,rowGlobal+numGlobalFeatures,rowGlobalInput);
    SymmetryHelpers::copyInputsWithSymmetry(rowSpatial, rowSpatialInput, 1, nnYLen, nnXLen, numSpatialFeatures, computeHandle->inputsUseNHWC, symmetry);
  }

  Buffers& buffers = computeHandle->buffers;

  CONSTTENSORMAP4 input(inputBuffers->spatialInput.data(), numSpatialFeatures, nnXLen, nnYLen, batchSize);
  CONSTTENSORMAP2 inputGlobal(inputBuffers->globalInput.data(), numGlobalFeatures, batchSize);

#define MAP4(NAME) TENSORMAP4 NAME(buffers.NAME.data(), buffers.NAME.dimension(0), buffers.NAME.dimension(1), buffers.NAME.dimension(2), batchSize)
#define MAP3(NAME) TENSORMAP3 NAME(buffers.NAME.data(), buffers.NAME.dimension(0), buffers.NAME.dimension(1), batchSize)
#define MAP2(NAME) TENSORMAP2 NAME(buffers.NAME.data(), buffers.NAME.dimension(0), batchSize)

  MAP2(inputMatMulOut);
  MAP4(trunk);
  MAP4(trunkScratch);
  MAP4(regularOut);
  MAP4(regularScratch);
  MAP4(midIn);
  MAP4(midScratch);
  MAP4(gpoolOut);
  MAP4(gpoolOut2);
  MAP2(gpoolConcat);
  MAP2(gpoolBias);
  MAP4(p1Out);
  MAP4(p1Out2);
  MAP4(g1Out);
  MAP4(g1Out2);
  MAP2(g1Concat);
  MAP2(g1Bias);
  MAP2(policyPass);
  MAP4(policy);
  MAP4(v1Out);
  MAP4(v1Out2);
  MAP2(v1Mean);
  MAP2(v2Out);
  MAP2(value);
  MAP2(scoreValue);
  MAP4(ownership);
  MAP3(mask);
  vector<float>& maskSum = buffers.maskSum;

  computeMaskSum(&mask,maskSum.data());

  computeHandle->model.apply(
    &input,
    &inputGlobal,
    &inputMatMulOut,
    &trunk,
    &trunkScratch,
    &regularOut,
    &regularScratch,
    &midIn,
    &midScratch,
    &gpoolOut,
    &gpoolOut2,
    &gpoolConcat,
    &gpoolBias,
    &p1Out,
    &p1Out2,
    &g1Out,
    &g1Out2,
    &g1Concat,
    &g1Bias,
    &policyPass,
    &policy,
    &v1Out,
    &v1Out2,
    &v1Mean,
    &v2Out,
    &value,
    &scoreValue,
    &ownership,
    &mask,
    maskSum.data()
  );

  assert(outputs.size() == batchSize);

  float* policyData = policy.data();
  float* policyPassData = policyPass.data();
  float* valueData = value.data();
  float* scoreValueData = scoreValue.data();
  float* ownershipData = ownership.data();

  for(int row = 0; row < batchSize; row++) {
    NNOutput* output = outputs[row];
    assert(output->nnXLen == nnXLen);
    assert(output->nnYLen == nnYLen);

    const float* policySrcBuf = policyData + row * inputBuffers->singlePolicyResultElts;
    float* policyProbs = output->policyProbs;

    //These are not actually correct, the client does the postprocessing to turn them into
    //policy probabilities and white game outcome probabilities
    //Also we don't fill in the nnHash here either
    SymmetryHelpers::copyOutputsWithSymmetry(policySrcBuf, policyProbs, 1, nnYLen, nnXLen, symmetry);
    policyProbs[inputBuffers->singlePolicyResultElts] = policyPassData[row];

    int numValueChannels = computeHandle->model.numValueChannels;
    assert(numValueChannels == 3);
    output->whiteWinProb = valueData[row * numValueChannels];
    output->whiteLossProb = valueData[row * numValueChannels + 1];
    output->whiteNoResultProb = valueData[row * numValueChannels + 2];

    //As above, these are NOT actually from white's perspective, but rather the player to move.
    //As usual the client does the postprocessing.
    if(output->whiteOwnerMap != NULL) {
      const float* ownershipSrcBuf = ownershipData + row * nnXLen * nnYLen;
      assert(computeHandle->model.numOwnershipChannels == 1);
      SymmetryHelpers::copyOutputsWithSymmetry(ownershipSrcBuf, output->whiteOwnerMap, 1, nnYLen, nnXLen, symmetry);
    }

    if(version >= 8) {
      int numScoreValueChannels = computeHandle->model.numScoreValueChannels;
      assert(numScoreValueChannels == 4);
      output->whiteScoreMean = scoreValueData[row * numScoreValueChannels];
      output->whiteScoreMeanSq = scoreValueData[row * numScoreValueChannels + 1];
      output->whiteLead = scoreValueData[row * numScoreValueChannels + 2];
      output->varTimeLeft = scoreValueData[row * numScoreValueChannels + 3];
    }
    else if(version >= 4) {
      int numScoreValueChannels = computeHandle->model.numScoreValueChannels;
      assert(numScoreValueChannels == 2);
      output->whiteScoreMean = scoreValueData[row * numScoreValueChannels];
      output->whiteScoreMeanSq = scoreValueData[row * numScoreValueChannels + 1];
      output->whiteLead = output->whiteScoreMean;
      output->varTimeLeft = 0;
    }
    else if(version >= 3) {
      int numScoreValueChannels = computeHandle->model.numScoreValueChannels;
      assert(numScoreValueChannels == 1);
      output->whiteScoreMean = scoreValueData[row * numScoreValueChannels];
      //Version 3 neural nets don't have any second moment output, implicitly already folding it in, so we just use the mean squared
      output->whiteScoreMeanSq = output->whiteScoreMean * output->whiteScoreMean;
      output->whiteLead = output->whiteScoreMean;
      output->varTimeLeft = 0;
    }
    else {
      ASSERT_UNREACHABLE;
    }
  }
}

void NeuralNet::printDevices() {
}

// FOR TESTING ---------------------------------------------------------------------------------------------------------
bool NeuralNet::testEvaluateConv(
  const ConvLayerDesc* desc,
  int batchSize,
  int nnXLen,
  int nnYLen,
  bool useFP16,
  bool useNHWC,
  const std::vector<float>& inputBuffer,
  std::vector<float>& outputBuffer
) {
  if(!useNHWC || useFP16)
    return false;
  ConvLayer layer(*desc);
  TENSORMAP4 inTensor(
    (float*)inputBuffer.data(), desc->inChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 outTensorBuf(desc->outChannels, nnXLen, nnYLen, batchSize);
  TENSORMAP4 outTensor(outTensorBuf);

  layer.apply(&inTensor, &outTensor, false);

  outputBuffer.resize(outTensorBuf.size());
  memcpy(outputBuffer.data(), outTensorBuf.data(), sizeof(SCALAR) * outTensorBuf.size());
  return true;
}

// Mask should be in 'NHW' format (no "C" channel).
bool NeuralNet::testEvaluateBatchNorm(
  const BatchNormLayerDesc* desc,
  int batchSize,
  int nnXLen,
  int nnYLen,
  bool useFP16,
  bool useNHWC,
  const std::vector<float>& inputBuffer,
  const std::vector<float>& maskBuffer,
  std::vector<float>& outputBuffer
) {
  if(!useNHWC || useFP16)
    return false;
  BatchNormLayer layer(*desc);
  TENSORMAP4 inTensor((float*)inputBuffer.data(), desc->numChannels, nnXLen, nnYLen, batchSize);
  TENSORMAP3 mask((float*)maskBuffer.data(), nnXLen, nnYLen, batchSize);
  TENSOR4 outTensorBuf(desc->numChannels, nnXLen, nnYLen, batchSize);
  TENSORMAP4 outTensor(outTensorBuf);

  layer.apply(false, &inTensor, &outTensor, &mask);

  outputBuffer.resize(outTensorBuf.size());
  memcpy(outputBuffer.data(), outTensorBuf.data(), sizeof(SCALAR) * outTensorBuf.size());
  return true;
}

// CR lpuchallafiore: test evaluate activation layer.
// CR lpuchallafiore: test evaluate matmul layer.
// CR lpuchallafiore: test evaluate matbias layer.

bool NeuralNet::testEvaluateResidualBlock(
  const ResidualBlockDesc* desc,
  int batchSize,
  int nnXLen,
  int nnYLen,
  bool useFP16,
  bool useNHWC,
  const std::vector<float>& inputBuffer,
  const std::vector<float>& maskBuffer,
  std::vector<float>& outputBuffer
) {
  if(!useNHWC || useFP16)
    return false;
  ResidualBlock block(*desc);
  TENSORMAP4 inTensor((float*)inputBuffer.data(), desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSORMAP3 mask((float*)maskBuffer.data(), nnXLen, nnYLen, batchSize);

  TENSOR4 trunkBuf(desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 trunkScratchBuf(desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 midBuf(desc->finalConv.inChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 midScratchBuf(desc->finalConv.inChannels, nnXLen, nnYLen, batchSize);

  TENSORMAP4 trunk(trunkBuf);
  TENSORMAP4 trunkScratch(trunkScratchBuf);
  TENSORMAP4 mid(midBuf);
  TENSORMAP4 midScratch(midScratchBuf);

  trunk = inTensor;

  block.apply(
    &trunk,
    &trunkScratch,
    NULL,
    NULL,
    &mid,
    &midScratch,
    NULL,
    NULL,
    NULL,
    NULL,
    &mask,
    NULL
  );

  outputBuffer.resize(trunk.size());
  memcpy(outputBuffer.data(), trunk.data(), sizeof(SCALAR) * trunk.size());
  return true;
}

bool NeuralNet::testEvaluateGlobalPoolingResidualBlock(
  const GlobalPoolingResidualBlockDesc* desc,
  int batchSize,
  int nnXLen,
  int nnYLen,
  bool useFP16,
  bool useNHWC,
  const std::vector<float>& inputBuffer,
  const std::vector<float>& maskBuffer,
  std::vector<float>& outputBuffer) {
  if(!useNHWC || useFP16)
    return false;

  GlobalPoolingResidualBlock block(*desc);

  TENSORMAP4 inTensor((float*)inputBuffer.data(), desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSORMAP3 mask((float*)maskBuffer.data(), nnXLen, nnYLen, batchSize);

  TENSOR4 trunkBuf(desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 trunkScratchBuf(desc->preBN.numChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 regularOutBuf(desc->finalConv.inChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 regularScratchBuf(desc->finalConv.inChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 gpoolOutBuf(desc->gpoolConv.outChannels, nnXLen, nnYLen, batchSize);
  TENSOR4 gpoolOut2Buf(desc->gpoolConv.outChannels, nnXLen, nnYLen, batchSize);
  TENSOR2 gpoolConcatBuf(desc->gpoolConv.outChannels*3, batchSize);
  TENSOR2 gpoolBiasBuf(desc->gpoolToBiasMul.outChannels, batchSize);

  TENSORMAP4 trunk(trunkBuf);
  TENSORMAP4 trunkScratch(trunkScratchBuf);
  TENSORMAP4 regularOut(regularOutBuf);
  TENSORMAP4 regularScratch(regularScratchBuf);
  TENSORMAP4 gpoolOut(gpoolOutBuf);
  TENSORMAP4 gpoolOut2(gpoolOut2Buf);
  TENSORMAP2 gpoolConcat(gpoolConcatBuf);
  TENSORMAP2 gpoolBias(gpoolBiasBuf);

  std::vector<float> maskSum(batchSize);
  computeMaskSum(&mask,maskSum.data());

  trunk = inTensor;

  block.apply(
    &trunk,
    &trunkScratch,
    &regularOut,
    &regularScratch,
    NULL,
    NULL,
    &gpoolOut,
    &gpoolOut2,
    &gpoolConcat,
    &gpoolBias,
    &mask,
    maskSum.data()
  );

  outputBuffer.resize(trunk.size());
  memcpy(outputBuffer.data(), trunk.data(), sizeof(SCALAR) * trunk.size());

  return true;
}

#endif  // USE_EIGEN_BACKEND

// CR lpuchallafiore: test evaluate Trunk
// CR lpuchallafiore: test evaluate Policy Head
// CR lpuchallafiore: test evaluate Value Head

max_stars_repo_name: Harald-Han/KataGo

max_issues_repo_name: Harald-Han/KataGo

max_forks_repo_name: Harald-Han/KataGo


text: #include <boost/tti/has_static_member_function.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include <boost/bind.hpp>
#include <boost/ref.hpp>
#include <boost/parameter.hpp>
#include <string>
#include <functional>

BOOST_PARAMETER_NAME(s1)
BOOST_PARAMETER_NAME(s2)
BOOST_PARAMETER_NAME(s3)

template <typename ArgumentPack>
std::string f(ArgumentPack const& args)
{
    std::string const& s1 = args[_s1];
    std::string const& s2 = args[_s2];
    typename boost::parameter::binding<
        ArgumentPack, tag::s3, std::string
    >::type s3 = args[
        _s3 || boost::bind(
            std::plus<std::string>()
          , boost::ref(s1)
          , boost::ref(s2)
        )
    ];
    return s3;
}

#include <boost/core/lightweight_test.hpp>

int main()
{
    std::string x = f((_s1="hello,", _s2=" world", _s3="hi world"));
    BOOST_TEST_EQ(x, std::string("hi world"));
    return boost::report_errors();
}