text:
// Boost.Polygon library voronoi_visualizer.cpp file

//          Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// See http://www.boost.org for updates, documentation, and revision history.

#include <iostream>
#include <vector>

#include <QtOpenGL/QGLWidget>
#include <QtGui/QtGui>

#include <boost/polygon/polygon.hpp>
#include <boost/polygon/voronoi.hpp>
using namespace boost::polygon;

#include "voronoi_visual_utils.hpp"

class GLWidget : public QGLWidget {
  Q_OBJECT

 public:
  explicit GLWidget(QMainWindow* parent = NULL) :
      QGLWidget(QGLFormat(QGL::SampleBuffers), parent),
      primary_edges_only_(false),
      internal_edges_only_(false) {
    startTimer(40);
  }

  QSize sizeHint() const {
    return QSize(600, 600);
  }

  void build(const QString& file_path) {
    // Clear all containers.
    clear();

    // Read data.
    read_data(file_path);

    // No data, don't proceed.
    if (!brect_initialized_) {
      return;
    }

    // Construct bounding rectangle.
    construct_brect();

    // Construct voronoi diagram.
    construct_voronoi(
        point_data_.begin(), point_data_.end(),
        segment_data_.begin(), segment_data_.end(),
        &vd_);

    // Color exterior edges.
    for (const_edge_iterator it = vd_.edges().begin();
         it != vd_.edges().end(); ++it) {
      if (!it->is_finite()) {
        color_exterior(&(*it));
      }
    }

    // Update view port.
    update_view_port();
  }

  void show_primary_edges_only() {
    primary_edges_only_ ^= true;
  }

  void show_internal_edges_only() {
    internal_edges_only_ ^= true;
  }

 protected:
  void initializeGL() {
    glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glEnable(GL_BLEND);
    glEnable(GL_POINT_SMOOTH);
  }

  void paintGL() {
    qglClearColor(QColor::fromRgb(255, 255, 255));
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    draw_points();
    draw_segments();
    draw_vertices();
    draw_edges();
  }

  void resizeGL(int width, int height) {
    int side = qMin(width, height);
    glViewport((width - side) / 2, (height - side) / 2, side, side);
  }

  void timerEvent(QTimerEvent* e) {
    update();
  }

 private:
  typedef double coordinate_type;
  typedef point_data<coordinate_type> point_type;
  typedef segment_data<coordinate_type> segment_type;
  typedef rectangle_data<coordinate_type> rect_type;
  typedef voronoi_builder<int> VB;
  typedef voronoi_diagram<coordinate_type> VD;
  typedef VD::cell_type cell_type;
  typedef VD::cell_type::source_index_type source_index_type;
  typedef VD::cell_type::source_category_type source_category_type;
  typedef VD::edge_type edge_type;
  typedef VD::cell_container_type cell_container_type;
  typedef VD::cell_container_type vertex_container_type;
  typedef VD::edge_container_type edge_container_type;
  typedef VD::const_cell_iterator const_cell_iterator;
  typedef VD::const_vertex_iterator const_vertex_iterator;
  typedef VD::const_edge_iterator const_edge_iterator;

  static const std::size_t EXTERNAL_COLOR = 1;

  void clear() {
    brect_initialized_ = false;
    point_data_.clear();
    segment_data_.clear();
    vd_.clear();
  }

  void read_data(const QString& file_path) {
    QFile data(file_path);
    if (!data.open(QFile::ReadOnly)) {
      QMessageBox::warning(
          this, tr("Voronoi Visualizer"),
          tr("Disable to open file ") + file_path);
    }
    QTextStream in_stream(&data);
    std::size_t num_points, num_segments;
    int x1, y1, x2, y2;
    in_stream >> num_points;
    for (std::size_t i = 0; i < num_points; ++i) {
      in_stream >> x1 >> y1;
      point_type p(x1, y1);
      update_brect(p);
      point_data_.push_back(p);
    }
    in_stream >> num_segments;
    for (std::size_t i = 0; i < num_segments; ++i) {
      in_stream >> x1 >> y1 >> x2 >> y2;
      point_type lp(x1, y1);
      point_type hp(x2, y2);
      update_brect(lp);
      update_brect(hp);
      segment_data_.push_back(segment_type(lp, hp));
    }
    in_stream.flush();
  }

  void update_brect(const point_type& point) {
    if (brect_initialized_) {
      encompass(brect_, point);
    } else {
      set_points(brect_, point, point);
      brect_initialized_ = true;
    }
  }

  void construct_brect() {
    double side = (std::max)(xh(brect_) - xl(brect_), yh(brect_) - yl(brect_));
    center(shift_, brect_);
    set_points(brect_, shift_, shift_);
    bloat(brect_, side * 1.2);
  }

  void color_exterior(const VD::edge_type* edge) {
    if (edge->color() == EXTERNAL_COLOR) {
      return;
    }
    edge->color(EXTERNAL_COLOR);
    edge->twin()->color(EXTERNAL_COLOR);
    const VD::vertex_type* v = edge->vertex1();
    if (v == NULL || !edge->is_primary()) {
      return;
    }
    v->color(EXTERNAL_COLOR);
    const VD::edge_type* e = v->incident_edge();
    do {
      color_exterior(e);
      e = e->rot_next();
    } while (e != v->incident_edge());
  }

  void update_view_port() {
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    rect_type view_rect = brect_;
    deconvolve(view_rect, shift_);
    glOrtho(xl(view_rect), xh(view_rect),
            yl(view_rect), yh(view_rect),
            -1.0, 1.0);
    glMatrixMode(GL_MODELVIEW);
  }

  void draw_points() {
    // Draw input points and endpoints of the input segments.
    glColor3f(0.0f, 0.5f, 1.0f);
    glPointSize(9);
    glBegin(GL_POINTS);
    for (std::size_t i = 0; i < point_data_.size(); ++i) {
      point_type point = point_data_[i];
      deconvolve(point, shift_);
      glVertex2f(point.x(), point.y());
    }
    for (std::size_t i = 0; i < segment_data_.size(); ++i) {
      point_type lp = low(segment_data_[i]);
      lp = deconvolve(lp, shift_);
      glVertex2f(lp.x(), lp.y());
      point_type hp = high(segment_data_[i]);
      hp = deconvolve(hp, shift_);
      glVertex2f(hp.x(), hp.y());
    }
    glEnd();
  }

  void draw_segments() {
    // Draw input segments.
    glColor3f(0.0f, 0.5f, 1.0f);
    glLineWidth(2.7f);
    glBegin(GL_LINES);
    for (std::size_t i = 0; i < segment_data_.size(); ++i) {
      point_type lp = low(segment_data_[i]);
      lp = deconvolve(lp, shift_);
      glVertex2f(lp.x(), lp.y());
      point_type hp = high(segment_data_[i]);
      hp = deconvolve(hp, shift_);
      glVertex2f(hp.x(), hp.y());
    }
    glEnd();
  }

  void draw_vertices() {
    // Draw voronoi vertices.
    glColor3f(0.0f, 0.0f, 0.0f);
    glPointSize(6);
    glBegin(GL_POINTS);
    for (const_vertex_iterator it = vd_.vertices().begin();
         it != vd_.vertices().end(); ++it) {
      if (internal_edges_only_ && (it->color() == EXTERNAL_COLOR)) {
        continue;
      }
      point_type vertex(it->x(), it->y());
      vertex = deconvolve(vertex, shift_);
      glVertex2f(vertex.x(), vertex.y());
    }
    glEnd();
  }
  void draw_edges() {
    // Draw voronoi edges.
    glColor3f(0.0f, 0.0f, 0.0f);
    glLineWidth(1.7f);
    for (const_edge_iterator it = vd_.edges().begin();
         it != vd_.edges().end(); ++it) {
      if (primary_edges_only_ && !it->is_primary()) {
        continue;
      }
      if (internal_edges_only_ && (it->color() == EXTERNAL_COLOR)) {
        continue;
      }
      std::vector<point_type> samples;
      if (!it->is_finite()) {
        clip_infinite_edge(*it, &samples);
      } else {
        point_type vertex0(it->vertex0()->x(), it->vertex0()->y());
        samples.push_back(vertex0);
        point_type vertex1(it->vertex1()->x(), it->vertex1()->y());
        samples.push_back(vertex1);
        if (it->is_curved()) {
          sample_curved_edge(*it, &samples);
        }
      }
      glBegin(GL_LINE_STRIP);
      for (std::size_t i = 0; i < samples.size(); ++i) {
        point_type vertex = deconvolve(samples[i], shift_);
        glVertex2f(vertex.x(), vertex.y());
      }
      glEnd();
    }
  }

  void clip_infinite_edge(
      const edge_type& edge, std::vector<point_type>* clipped_edge) {
    const cell_type& cell1 = *edge.cell();
    const cell_type& cell2 = *edge.twin()->cell();
    point_type origin, direction;
    // Infinite edges could not be created by two segment sites.
    if (cell1.contains_point() && cell2.contains_point()) {
      point_type p1 = retrieve_point(cell1);
      point_type p2 = retrieve_point(cell2);
      origin.x((p1.x() + p2.x()) * 0.5);
      origin.y((p1.y() + p2.y()) * 0.5);
      direction.x(p1.y() - p2.y());
      direction.y(p2.x() - p1.x());
    } else {
      origin = cell1.contains_segment() ?
          retrieve_point(cell2) :
          retrieve_point(cell1);
      segment_type segment = cell1.contains_segment() ?
          retrieve_segment(cell1) :
          retrieve_segment(cell2);
      coordinate_type dx = high(segment).x() - low(segment).x();
      coordinate_type dy = high(segment).y() - low(segment).y();
      if ((low(segment) == origin) ^ cell1.contains_point()) {
        direction.x(dy);
        direction.y(-dx);
      } else {
        direction.x(-dy);
        direction.y(dx);
      }
    }
    coordinate_type side = xh(brect_) - xl(brect_);
    coordinate_type koef =
        side / (std::max)(fabs(direction.x()), fabs(direction.y()));
    if (edge.vertex0() == NULL) {
      clipped_edge->push_back(point_type(
          origin.x() - direction.x() * koef,
          origin.y() - direction.y() * koef));
    } else {
      clipped_edge->push_back(
          point_type(edge.vertex0()->x(), edge.vertex0()->y()));
    }
    if (edge.vertex1() == NULL) {
      clipped_edge->push_back(point_type(
          origin.x() + direction.x() * koef,
          origin.y() + direction.y() * koef));
    } else {
      clipped_edge->push_back(
          point_type(edge.vertex1()->x(), edge.vertex1()->y()));
    }
  }

  void sample_curved_edge(
      const edge_type& edge,
      std::vector<point_type>* sampled_edge) {
    coordinate_type max_dist = 1E-3 * (xh(brect_) - xl(brect_));
    point_type point = edge.cell()->contains_point() ?
        retrieve_point(*edge.cell()) :
        retrieve_point(*edge.twin()->cell());
    segment_type segment = edge.cell()->contains_point() ?
        retrieve_segment(*edge.twin()->cell()) :
        retrieve_segment(*edge.cell());
    voronoi_visual_utils<coordinate_type>::discretize(
        point, segment, max_dist, sampled_edge);
  }

  point_type retrieve_point(const cell_type& cell) {
    source_index_type index = cell.source_index();
    source_category_type category = cell.source_category();
    if (category == SOURCE_CATEGORY_SINGLE_POINT) {
      return point_data_[index];
    }
    index -= point_data_.size();
    if (category == SOURCE_CATEGORY_SEGMENT_START_POINT) {
      return low(segment_data_[index]);
    } else {
      return high(segment_data_[index]);
    }
  }

  segment_type retrieve_segment(const cell_type& cell) {
    source_index_type index = cell.source_index() - point_data_.size();
    return segment_data_[index];
  }

  point_type shift_;
  std::vector<point_type> point_data_;
  std::vector<segment_type> segment_data_;
  rect_type brect_;
  VB vb_;
  VD vd_;
  bool brect_initialized_;
  bool primary_edges_only_;
  bool internal_edges_only_;
};

class MainWindow : public QWidget {
  Q_OBJECT

 public:
  MainWindow() {
    glWidget_ = new GLWidget();
    file_dir_ = QDir(QDir::currentPath(), tr("*.txt"));
    file_name_ = tr("");

    QHBoxLayout* centralLayout = new QHBoxLayout;
    centralLayout->addWidget(glWidget_);
    centralLayout->addLayout(create_file_layout());
    setLayout(centralLayout);

    update_file_list();
    setWindowTitle(tr("Voronoi Visualizer"));
    layout()->setSizeConstraint(QLayout::SetFixedSize);
  }

 private slots:
  void primary_edges_only() {
    glWidget_->show_primary_edges_only();
  }

  void internal_edges_only() {
    glWidget_->show_internal_edges_only();
  }

  void browse() {
    QString new_path = QFileDialog::getExistingDirectory(
        0, tr("Choose Directory"), file_dir_.absolutePath());
    if (new_path.isEmpty()) {
      return;
    }
    file_dir_.setPath(new_path);
    update_file_list();
  }

  void build() {
    file_name_ = file_list_->currentItem()->text();
    QString file_path = file_dir_.filePath(file_name_);
    message_label_->setText("Building...");
    glWidget_->build(file_path);
    message_label_->setText("Double click the item to build voronoi diagram:");
    setWindowTitle(tr("Voronoi Visualizer - ") + file_path);
  }

  void print_scr() {
    if (!file_name_.isEmpty()) {
      QImage screenshot = glWidget_->grabFrameBuffer(true);
      QString output_file = file_dir_.absolutePath() + tr("/") +
          file_name_.left(file_name_.indexOf('.')) + tr(".png");
      screenshot.save(output_file, 0, -1);
    }
  }

 private:
  QGridLayout* create_file_layout() {
    QGridLayout* file_layout = new QGridLayout;

    message_label_ = new QLabel("Double click item to build voronoi diagram:");

    file_list_ = new QListWidget();
    file_list_->connect(file_list_,
                        SIGNAL(itemDoubleClicked(QListWidgetItem*)),
                        this,
                        SLOT(build()));

    QCheckBox* primary_checkbox = new QCheckBox("Show primary edges only.");
    connect(primary_checkbox, SIGNAL(clicked()),
        this, SLOT(primary_edges_only()));

    QCheckBox* internal_checkbox = new QCheckBox("Show internal edges only.");
    connect(internal_checkbox, SIGNAL(clicked()),
        this, SLOT(internal_edges_only()));

    QPushButton* browse_button =
        new QPushButton(tr("Browse Input Directory"));
    connect(browse_button, SIGNAL(clicked()), this, SLOT(browse()));
    browse_button->setMinimumHeight(50);

    QPushButton* print_scr_button = new QPushButton(tr("Make Screenshot"));
    connect(print_scr_button, SIGNAL(clicked()), this, SLOT(print_scr()));
    print_scr_button->setMinimumHeight(50);

    file_layout->addWidget(message_label_, 0, 0);
    file_layout->addWidget(file_list_, 1, 0);
    file_layout->addWidget(primary_checkbox, 2, 0);
    file_layout->addWidget(internal_checkbox, 3, 0);
    file_layout->addWidget(browse_button, 4, 0);
    file_layout->addWidget(print_scr_button, 5, 0);

    return file_layout;
  }

  void update_file_list() {
    QFileInfoList list = file_dir_.entryInfoList();
    file_list_->clear();
    if (file_dir_.count() == 0) {
      return;
    }
    QFileInfoList::const_iterator it;
    for (it = list.begin(); it != list.end(); it++) {
      file_list_->addItem(it->fileName());
    }
    file_list_->setCurrentRow(0);
  }

  QDir file_dir_;
  QString file_name_;
  GLWidget* glWidget_;
  QListWidget* file_list_;
  QLabel* message_label_;
};

int main(int argc, char* argv[]) {
  QApplication app(argc, argv);
  MainWindow window;
  window.show();
  return app.exec();
}

#include "voronoi_visualizer.moc"


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2021.
//
// 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: $
// --------------------------------------------------------------------------
//

#include <OpenMS/TRANSFORMATIONS/RAW2PEAK/TwoDOptimization.h>

#include <Eigen/Core>
#include <unsupported/Eigen/NonLinearOptimization>
#include <boost/math/special_functions/acosh.hpp> // TODO replace with std::acosh

using Eigen::VectorXd;

namespace OpenMS
{

  TwoDOptimization::TwoDOptimization() :
    DefaultParamHandler("TwoDOptimization")
  {
    // 2D optimization parameters
    defaults_.setValue("penalties:position", 0.0, "If the position changes more than 0.2Da during the fitting it can be penalized");
    defaults_.setValue("penalties:height", 1.0, "penalty term for the fitting of the intensity:" \
                                                "If it gets negative during the fitting it can be penalized.");
    defaults_.setValue("penalties:left_width", 0.0, "penalty term for the fitting of the left width:" \
                                                    "If the left width gets too broad or negative during the fitting it can be penalized.");
    defaults_.setValue("penalties:right_width", 0.0, "penalty term for the fitting of the right width:" \
                                                     "If the right width gets too broad or negative during the fitting it can be penalized.");
    defaults_.setValue("2d:tolerance_mz", 2.2, "mz tolerance for cluster construction", {"advanced"});
    defaults_.setValue("2d:max_peak_distance", 1.2, "maximal peak distance in mz in a cluster", {"advanced"});
    defaults_.setValue("iterations", 10, "maximal number of iterations for the fitting step");

    defaultsToParam_();
    updateMembers_();
  }

  TwoDOptimization::TwoDOptimization(const TwoDOptimization & opt) :
    DefaultParamHandler(opt)
  {
    updateMembers_();
  }

  TwoDOptimization & TwoDOptimization::operator=(const TwoDOptimization & opt)
  {
    if (&opt == this)
    {
      return *this;
    }
    DefaultParamHandler::operator=(opt);
    updateMembers_();

    return *this;
  }

  void TwoDOptimization::findMatchingPeaks_(std::multimap<double, IsotopeCluster>::iterator & it, PeakMap & ms_exp)
  {
    IsotopeCluster::ChargedIndexSet::const_iterator iter = it->second.peaks.begin();
    for (; iter != it->second.peaks.end(); ++iter)
    {

      double mz = (ms_exp[iter->first][iter->second]).getMZ();
      mz *= 10;
      matching_peaks_[(Int)(mz + 0.5)].push_back(PeakIndex(iter->first, iter->second));
    }

#ifdef DEBUG_2D
    std::map<Int, PeakIndex>::iterator it2 = matching_peaks_.begin();
    for (; it2 != matching_peaks_.end(); ++it2)
    {
      std::cout << it2->first << " has " << it2->second.size() << " elements:" << std::endl;
      for (Size i = 0; i < it2->second.size(); ++i)
      {
        std::cout << it2->second[i]->getPeak(ms_exp).getMZ() << "\t";
      }
      std::cout << std::endl;
    }
#endif

  }

  // Finds the neighbour of the peak denoted by @p current_mz in the previous scan
  std::vector<double>::iterator TwoDOptimization::searchInScan_(std::vector<double>::iterator scan_begin,
                                                                    std::vector<double>::iterator scan_end,
                                                                    double current_mz)
  {

    // perform binary search to find the neighbour in rt dimension
    //  lower_bound finds the peak with m/z current_mz or the next larger peak if this peak does not exist.
    std::vector<double>::iterator insert_iter = lower_bound(scan_begin, scan_end, current_mz);

    // the peak found by lower_bound does not have to be the closest one, therefore we have
    // to check both neighbours
    if (insert_iter == scan_end)   // we are at the and have only one choice
    {
      return --insert_iter;
    }
    else
    {
      // if the found peak is at the beginning of the spectrum,
      // there is not much we can do.
      if (insert_iter == scan_begin)
      {
        return insert_iter;
      }
      else           // see if the next smaller one fits better
      {
        double delta_mz = fabs(*insert_iter - current_mz);
        --insert_iter;

        if (fabs(*insert_iter - current_mz) < delta_mz)
        {
          return insert_iter;                       // peak to the left is closer (in m/z dimension)
        }
        else
        {
          return ++insert_iter;                          // peak to the right is closer
        }
      }
    }

  } // end of searchInScan_

  void TwoDOptimization::updateMembers_()
  {
    penalties_.height = (double)param_.getValue("penalties:height");
    penalties_.pos = (double)param_.getValue("penalties:position");
    penalties_.lWidth = (double)param_.getValue("penalties:left_width");
    penalties_.rWidth = (double)param_.getValue("penalties:right_width");
    max_peak_distance_ = (double)param_.getValue("2d:max_peak_distance");
    tolerance_mz_ = (double)param_.getValue("2d:tolerance_mz");
    max_iteration_ = (UInt)param_.getValue("iterations");

  }

  int TwoDOptimization::TwoDOptFunctor::operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec)
  {

    // x contains the parameters to be optimized.
    // In our case these are the weighted average mz-position and left and right width for all
    // matching peaks in all considered scans. Additionally the intensities of all peaks are stored.
    // Params might contain any additional parameters. We handle these using class members
    // instead.
    // The vector f is supposed to contain the result when we return from this function.
    double computed_signal, experimental_signal, step, last_position;
    double p_height, p_position, p_width;
    Int count = 0;
    Int counter_posf = 0;
    const std::vector<std::pair<SignedSize, SignedSize> > & signal2D = m_data->signal2D;
    std::multimap<double, IsotopeCluster>::iterator iso_map_iter = m_data->iso_map_iter;
    Size total_nr_peaks = m_data->total_nr_peaks;
    const std::map<Int, std::vector<PeakIndex> > & matching_peaks = m_data->matching_peaks;
    const PeakMap & picked_peaks = m_data->picked_peaks;
    PeakMap::ConstIterator raw_data_first = m_data->raw_data_first;
    const OptimizationFunctions::PenaltyFactorsIntensity & penalties = m_data->penalties;

    Size num_scans = signal2D.size() / 2;
    IsotopeCluster::ChargedIndexSet::iterator peak_iter = iso_map_iter->second.peaks.begin();
    fvec.setZero();

    //iterate over all scans
    for (Size current_scan = 0; current_scan < num_scans; ++current_scan)
    {
      Size curr_scan_idx = current_scan + iso_map_iter->second.peaks.begin()->first;
      double current_position = ((raw_data_first
                           + signal2D[2 * current_scan].first)->begin()
                          + signal2D[2 * current_scan].second)->getMZ();
      //iterate over all points of the signal
      for (Int current_point = 1;
           current_point +  signal2D[2 * current_scan].second
           <= signal2D[2 * current_scan + 1].second;
           ++current_point)
      {
        last_position = current_position;

        computed_signal   = 0.;
        current_position  = ((raw_data_first
                              + signal2D[2 * current_scan].first)->begin()
                             + signal2D[2 * current_scan].second + current_point)->getMZ();
        experimental_signal = ((raw_data_first
                                + signal2D[2 * current_scan].first)->begin()
                               + signal2D[2 * current_scan].second + current_point)->getIntensity();
        step = current_position - last_position;
#ifdef DEBUG_2D
        std::cout << "experimental signal rt " << (raw_data_first
                                                   + signal2D[2 * current_scan].first)->getRT()
                  << "\tmz " << ((raw_data_first
                        + signal2D[2 * current_scan].first)->begin()
                       + signal2D[2 * current_scan].second + current_point)->getMZ()
        << "\tint " << ((raw_data_first
                         + signal2D[2 * current_scan].first)->begin()
                        + signal2D[2 * current_scan].second + current_point)->getIntensity() << std::endl;
#endif

        // size_t current_peak = 0;
        peak_iter = iso_map_iter->second.peaks.begin();
        while (peak_iter != iso_map_iter->second.peaks.end() && peak_iter->first != curr_scan_idx)
        {
          ++peak_iter;
        }
        //iterate over all peaks of the current scan
        while (peak_iter != iso_map_iter->second.peaks.end() && peak_iter->first == curr_scan_idx)
        {
          Int peak_idx = distance(iso_map_iter->second.peaks.begin(), peak_iter);
          double mz_in_hash = ((picked_peaks[peak_iter->first]).begin() + peak_iter->second)->getMZ() * 10;
          std::map<Int, std::vector<PeakIndex> >::const_iterator  m_spec_iter = matching_peaks.begin();
          Int map_idx = 0;
          while (m_spec_iter->first != (Int)(mz_in_hash + 0.5))
          {
            ++map_idx;
            ++m_spec_iter;
          }
          //Store the current parameters for this peak
          p_position    = x(total_nr_peaks + 3 * map_idx);
          p_height      = x(peak_idx);
          p_width       = (current_position <= p_position) ?
                          x(total_nr_peaks + 3 * map_idx + 1) :
                          x(total_nr_peaks + 3 * map_idx + 2);
          ++count;

          //is it a Lorentz or a Sech - Peak?
          if ((PeakShape::Type)(Int) Math::round((picked_peaks[peak_iter->first]).getFloatDataArrays()[5][peak_iter->second]) == PeakShape::LORENTZ_PEAK)
          {
#ifdef DEBUG_2D
            std::cout << "p_height " << p_height << "\tp_position " << p_position << "\tcurrent_position "
                      << current_position << std::endl;
#endif
            computed_signal += p_height / (1. + pow(p_width * (current_position - p_position), 2));
          }
          // if it's a Sech - Peak
          else
          {
#ifdef DEBUG_2D
            std::cout << "p_height " << p_height << "\tp_position " << p_position << "\tcurrent_position "
                      << current_position << std::endl;
#endif
            computed_signal += p_height / pow(cosh(p_width * (current_position - p_position)), 2);
          }
          // ++current_peak;
          ++peak_iter;

        }                        // end while
#ifdef DEBUG_2D
        std::cout << "computed vs experimental signal: " << computed_signal << "\t"
                  << experimental_signal << std::endl;
#endif
        fvec(counter_posf) = step * (computed_signal - experimental_signal);
        ++counter_posf;
      }

    }

    // penalties : especially negative heights have to be penalised
    double penalty = 0.;

    //iterate over all peaks again to compute the penalties
    // first look at all positions and width parameters
    UInt peak = 0, current_peak = 0;
    std::map<Int, std::vector<PeakIndex> >::const_iterator map_iter = matching_peaks.begin();
    for (; map_iter != matching_peaks.end(); ++map_iter)
    {
      std::vector<PeakIndex>::const_iterator vec_iter = map_iter->second.begin();
      double old_position = 0, old_width_l = 0, old_width_r = 0;
      double weight = 0;
      for (; vec_iter != map_iter->second.end(); ++vec_iter)
      {
        // intensity at peak position is stored in the meta data array, not as intensity of the peak
        double old_height = picked_peaks[vec_iter->spectrum].getFloatDataArrays()[1][vec_iter->peak];
        weight += old_height;
        old_position += (vec_iter)->getPeak(picked_peaks).getMZ() * old_height;
        old_width_l += picked_peaks[vec_iter->spectrum].getFloatDataArrays()[3][vec_iter->peak] * old_height;
        old_width_r += picked_peaks[vec_iter->spectrum].getFloatDataArrays()[4][vec_iter->peak] * old_height;

        double p_height = x(peak);
        ++peak;

        if (p_height < 1)
        {
          penalty += 1000000 * penalties.height * pow(fabs(p_height - old_height), 2);
        }

      }
      old_position /= weight;
      old_width_l /= weight;
      old_width_r /= weight;

      double p_position   = x(total_nr_peaks + 3 * current_peak);
      double p_width_l    = x(total_nr_peaks + 3 * current_peak + 1);
      double p_width_r    = x(total_nr_peaks + 3 * current_peak + 2);
      if (p_width_l < 0)
      {
        penalty += 1e7 * penalties.lWidth * pow(fabs(p_width_l - old_width_l), 2);
      }
      else if (p_width_l < 1)
      {
        penalty += 1000 * penalties.lWidth * pow(fabs(p_width_l - old_width_l), 2);
      }
      if (p_width_r < 0)
      {
        penalty += 1e7 * penalties.rWidth * pow(fabs(p_width_r - old_width_r), 2);
      }
      else if (p_width_r < 1)
      {
        penalty += 1000 * penalties.rWidth * pow(fabs(p_width_r - old_width_r), 2);
      }
      if (p_position < 0)
      {
        penalty += 100 * penalties.pos * pow(p_position - old_position, 2);
      }
      if (fabs(old_width_r - p_width_r) > 1)
      {
        penalty += 1000 * penalties.rWidth * pow(old_width_r - p_width_r, 2);
      }
      if (fabs(old_width_l - p_width_l) > 1)
      {
        penalty += 1000 * penalties.lWidth * pow(old_width_l - p_width_l, 2);
      }
      if (fabs(old_position - p_position) > 0.2)
      {
        penalty += 1000 * penalties.pos * pow(p_position - old_position, 2);
      }

      ++current_peak;
    }

    fvec(fvec.size()-1) = penalty;

    return 0;
  }

  // compute Jacobian matrix for the different parameters
  int TwoDOptimization::TwoDOptFunctor::df(const Eigen::VectorXd &x, Eigen::MatrixXd &J)
  {
    // For the conventions on x and params c.f. the commentary in residual()
    //
    // The matrix J is supposed to contain the result when we return from this function.
    // Note: expects the Jacobian as follows:
    //          - each row corresponds to one data point
    //          - each column corresponds to one parameter

    double last_position, step;
    double p_height, p_position, p_width;
    double diff, denom_inv, ddl_left, ddl_right, ddx0, sinh_term;
    Int count = 0;
    Int counter_posf = 0;

    const std::vector<std::pair<SignedSize, SignedSize> > & signal2D = m_data->signal2D;
    std::multimap<double, IsotopeCluster>::iterator iso_map_iter = m_data->iso_map_iter;
    Size total_nr_peaks = m_data->total_nr_peaks;
    const std::map<Int, std::vector<PeakIndex> > & matching_peaks = m_data->matching_peaks;
    std::vector<double> ov_weight(matching_peaks.size(), 0);
    const PeakMap & picked_peaks = m_data->picked_peaks;
    PeakMap::ConstIterator raw_data_first = m_data->raw_data_first;
    const OptimizationFunctions::PenaltyFactorsIntensity & penalties = m_data->penalties;
//          std::vector<double> &positions=static_cast<TwoDOptimization::Data*> (params) ->positions;
//          std::vector<double> &signal=static_cast<TwoDOptimization::Data*> (params) ->signal;
    IsotopeCluster::ChargedIndexSet::iterator peak_iter = iso_map_iter->second.peaks.begin();
    Size num_scans = signal2D.size() / 2;
    //iterate over all scans
    for (Size current_scan = 0; current_scan < num_scans; ++current_scan)
    {
      Size curr_scan_idx = current_scan + iso_map_iter->second.peaks.begin()->first;
      double current_position = ((raw_data_first
                           + signal2D[2 * current_scan].first)->begin()
                          + signal2D[2 * current_scan].second)->getMZ();
      // iterate over all points of the signal
      for (Int current_point = 1;
           current_point +  signal2D[2 * current_scan].second
           <= signal2D[2 * current_scan + 1].second;
           ++current_point)
      {
        last_position = current_position;
        current_position  = ((raw_data_first
                              + signal2D[2 * current_scan].first)->begin()
                             + signal2D[2 * current_scan].second + current_point)->getMZ();
//       double experimental_signal = ((raw_data_first + signal2D[2 * current_scan].first)->begin()
//                               + signal2D[2 * current_scan].second + current_point)->getIntensity();

        step = current_position - last_position;

#ifdef DEBUG_2D
        std::cout << "experimental signal rt " << (raw_data_first
                                                   + signal2D[2 * current_scan].first)->getRT()
                  << "\tmz " << ((raw_data_first
                        + signal2D[2 * current_scan].first)->begin()
                       + signal2D[2 * current_scan].second + current_point)->getMZ()
        << "\tint " << ((raw_data_first
                         + signal2D[2 * current_scan].first)->begin()
                        + signal2D[2 * current_scan].second + current_point)->getIntensity() << std::endl;

#endif

        // size_t current_peak = 0;
        peak_iter = iso_map_iter->second.peaks.begin();
        while (peak_iter != iso_map_iter->second.peaks.end() && peak_iter->first != curr_scan_idx)
        {
          ++peak_iter;
        }
        //iterate over all peaks of the current scan
        while (peak_iter != iso_map_iter->second.peaks.end() && peak_iter->first == curr_scan_idx)
        {
          Int peak_idx = distance(iso_map_iter->second.peaks.begin(), peak_iter);
          double mz_in_hash = ((picked_peaks[peak_iter->first]).begin() + peak_iter->second)->getMZ() * 10;
          std::map<Int, std::vector<PeakIndex> >::const_iterator  m_spec_iter =  matching_peaks.begin();
          Int map_idx = 0;
          while (m_spec_iter->first != (Int)(mz_in_hash + 0.5))
          {
            ++map_idx;
            ++m_spec_iter;
          }
          // if the current peak is in the reference scan take all parameters from the vector x
#ifdef DEBUG_2D
          std::cout << "ref_scan : " << current_peak << "\t "
                    << x(3 * current_peak) << "\t" << x(3 * current_peak + 1)
                    << "\t" << x(3 * current_peak + 2) << std::endl;
#endif
          // Store the current parameters for this peak
          p_position    = x(total_nr_peaks + 3 * map_idx);
          p_height      = x(peak_idx);
          p_width       = (current_position <= p_position) ?
                          x(total_nr_peaks + 3 * map_idx + 1) :
                          x(total_nr_peaks + 3 * map_idx + 2);
          ++count;
          double weight = step * picked_peaks[peak_iter->first].getFloatDataArrays()[1][peak_iter->second];
          ov_weight[map_idx] += weight;
          double ddx0_old = J(counter_posf, total_nr_peaks + 3 * map_idx);
          double ddl_left_old = J(counter_posf, total_nr_peaks + 3 * map_idx + 1);
          double ddl_right_old = J(counter_posf, total_nr_peaks + 3 * map_idx + 2);
          //is it a Lorentz or a Sech - Peak?

          if ((PeakShape::Type)(Int) Math::round((picked_peaks[peak_iter->first]).getFloatDataArrays()[5][peak_iter->second]) == PeakShape::LORENTZ_PEAK)
          {
            diff      = current_position - p_position;
            // partial derivative with respect to the height,...
            denom_inv = 1. / (1. + pow(p_width * diff, 2));
            // left width,...
            ddl_left  = (current_position <= p_position)
                        ? -2 * p_height * pow(diff, 2) * p_width * pow(denom_inv, 2) :
                          0;
            // right width ...
            ddl_right = (current_position  > p_position)
                        ? -2 * p_height * pow(diff, 2) * p_width * pow(denom_inv, 2) :
                          0;

            // and position
            ddx0 = 2 * p_height * pow(p_width, 2) * diff * pow(denom_inv, 2);

            J(counter_posf, total_nr_peaks + 3 * map_idx) = ddx0 * weight  + ddx0_old;
            J(counter_posf, peak_idx) = step * denom_inv;
            J(counter_posf, total_nr_peaks + 3 * map_idx + 1) = ddl_left * weight + ddl_left_old;
            J(counter_posf, total_nr_peaks + 3 * map_idx + 2) = ddl_right * weight + ddl_right_old;
          }
          // if it's a Sech - Peak
          else
          {
            diff      = current_position - p_position;
            denom_inv = 1. / cosh(p_width * diff);
            // The remaining computations are not stable if denom_inv == 0. In that case, we are far away from the peak
            // and can assume that all derivatives vanish
            sinh_term = (fabs(denom_inv) < 1e-6) ? 0.0 : sinh(p_width * diff);

            ddl_left  = (current_position <= p_position)
                        ? -2 * p_height * sinh_term * diff * pow(denom_inv, 3) :
                          0;
            ddl_right = (current_position  > p_position)
                        ? -2 * p_height * sinh_term * diff * pow(denom_inv, 3) :
                          0;

            ddx0      = 2 * p_height * p_width * sinh_term * pow(denom_inv, 3);

            J(counter_posf, total_nr_peaks + 3 * map_idx) = ddx0 * weight + ddx0_old;
            J(counter_posf, peak_idx) = step * pow(denom_inv, 2);
            J(counter_posf, total_nr_peaks + 3 * map_idx  + 1) = ddl_left * weight + ddl_left_old;
            J(counter_posf, total_nr_peaks + 3 * map_idx  + 2) = ddl_right * weight + ddl_right_old;
          }
          // ++current_peak;
          ++peak_iter;

        }                        // end while

        ++counter_posf;
      }

    }

    for (Size cluster = 0; cluster < matching_peaks.size(); ++cluster)
    {
      for (int j = 0; j < J.rows() - 1; ++j)
      {
        J(j, total_nr_peaks + 3 * cluster)
              = J(j, total_nr_peaks + 3 * cluster) / ov_weight[cluster];
        J(j, total_nr_peaks + 3 * cluster + 1)
              = J(j, total_nr_peaks + 3 * cluster + 1) / ov_weight[cluster];
        J(j, total_nr_peaks + 3 * cluster + 2)
              = J(j, total_nr_peaks + 3 * cluster + 2) / ov_weight[cluster];
      }
    }

    //iterate over all peaks again to compute the penalties
    // first look at all positions and width parameters
    UInt peak = 0, current_peak = 0;
    std::map<Int, std::vector<PeakIndex> >::const_iterator map_iter = matching_peaks.begin();

    for (; map_iter != matching_peaks.end(); ++map_iter)
    {
      std::vector<PeakIndex>::const_iterator vec_iter   = map_iter->second.begin();
      double old_position = 0, old_width_l = 0, old_width_r = 0;
      double weight = 0;
      double penalty_h = 0, penalty_l = 0, penalty_r = 0, penalty_p = 0;
      for (; vec_iter != map_iter->second.end(); ++vec_iter)
      {
        double old_height = picked_peaks[vec_iter->spectrum].getFloatDataArrays()[1][vec_iter->peak];
        weight += old_height;
        old_position += (vec_iter)->getPeak(picked_peaks).getMZ() * old_height;
        old_width_l += picked_peaks[vec_iter->spectrum].getFloatDataArrays()[3][vec_iter->peak] * old_height;
        old_width_r += picked_peaks[vec_iter->spectrum].getFloatDataArrays()[4][vec_iter->peak] * old_height;

        double p_height     = x(peak);

        double penalty_height = 2. * penalties.height * fabs(p_height - old_height);
        if (p_height < 1)
        {
          penalty_h += 1000000 * penalty_height;
        }
        J(counter_posf, peak) = penalty_h;
        ++peak;
      }
      old_position /= weight;
      old_width_l /= weight;
      old_width_r /= weight;

      // std::cout << old_position << "vs. ";
      double p_position   = x(total_nr_peaks + 3 * current_peak);
      double p_width_l    = x(total_nr_peaks + 3 * current_peak + 1);
      double p_width_r    = x(total_nr_peaks + 3 * current_peak + 2);
      double penalty_lwidth = 2. * penalties.lWidth * fabs(p_width_l - old_width_l);
      double penalty_rwidth = 2. * penalties.rWidth * fabs(p_width_r - old_width_r);
      double penalty_pos    = 2. * penalties.pos * fabs(p_position - old_position);
      //std::cout << p_position<<std::endl;
#ifdef DEBUG_2D
      std::cout << "penalty_lwidth " << penalty_lwidth << "penalty_rwidth " << penalty_rwidth
                << "penalty_pos " << penalty_pos << std::endl;
#endif
      if (p_width_l < 0)
      {
        penalty_l += 1e7 * penalty_lwidth;
      }
      else if (p_width_l < 1)
      {
        penalty_l += 2000 * penalties.lWidth * (fabs(p_width_l - old_width_l));
      }
      if (p_width_r < 0)
      {
        penalty_r += 1e7 * penalty_rwidth;
      }
      else if (p_width_r < 1)
      {
        penalty_r += 2000 * penalties.rWidth * (fabs(p_width_r - old_width_r));
      }
      if (p_position < 0)
      {
        penalty_p += 200 * penalty_pos;
      }

      if (fabs(old_position - p_position) > 0.2)
      {
        penalty_p += 2000 * penalties.pos * fabs(p_position - old_position);
      }
      if (fabs(old_width_r - p_width_r) > 1)
      {
        penalty_r += 1000 * penalty_rwidth;
      }
      if (fabs(old_width_l - p_width_l) > 1)
      {
        penalty_l += 1000 * penalty_lwidth;
      }

      J(counter_posf, total_nr_peaks + 3 * current_peak + 1) = penalty_l;
      J(counter_posf, total_nr_peaks + 3 * current_peak + 2) = penalty_r;
      J(counter_posf, total_nr_peaks + 3 * current_peak) = penalty_p;

      ++current_peak;
    }
    return 0;
  }

  void TwoDOptimization::optimize(InputSpectrumIterator first, InputSpectrumIterator last, PeakMap& ms_exp, bool real2D)
  {
    //#define DEBUG_2D
    //check if the input maps have the same number of spectra
    if ((UInt)distance(first, last) != ms_exp.size())
    {
      throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Error in Two2Optimization: Raw and peak map do not have the same number of spectra");
    }
    //do nothing if there are no scans
    if (ms_exp.empty())
    {
      return;
    }
    //check if required meta data arrays are present (for each scan)
    for (Size i = 0; i < ms_exp.size(); ++i)
    {
      //check if enough meta data arrays are present
      if (ms_exp[i].getFloatDataArrays().size() < 6)
      {
        throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Error in Two2Optimization: Not enough meta data arrays present (1:area, 5:shape, 3:left width, 4:right width)");
      }
      bool area = ms_exp[i].getFloatDataArrays()[1].getName() == "maximumIntensity";
      bool wleft = ms_exp[i].getFloatDataArrays()[3].getName() == "leftWidth";
      bool wright = ms_exp[i].getFloatDataArrays()[4].getName() == "rightWidth";
      bool shape = ms_exp[i].getFloatDataArrays()[5].getName() == "peakShape";

      if (!area || !wleft || !wright || !shape)
      {
        throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Error in Two2Optimization: One or several meta data arrays missing (1:intensity, 5:shape, 3:left width, 4:right width)");
      }
    }
    real_2D_ = real2D;
    typedef MSSpectrum SpectrumType;

    typename PeakMap::Iterator ms_exp_it = ms_exp.begin();
    typename PeakMap::Iterator ms_exp_it_end = ms_exp.end();
    if (ms_exp.empty())
    {
      return;
    }
    // stores the monoisotopic peaks of isotopic clusters
    std::vector<double> iso_last_scan;
    std::vector<double> iso_curr_scan;
    std::vector<std::multimap<double, IsotopeCluster>::iterator> clusters_last_scan;
    std::vector<std::multimap<double, IsotopeCluster>::iterator> clusters_curr_scan;
    std::multimap<double, IsotopeCluster>::iterator cluster_iter;
    double current_rt = ms_exp_it->getRT(), last_rt  = 0;

    // retrieve values for accepted peaks distances
    max_peak_distance_ = param_.getValue("2d:max_peak_distance");
    double tolerance_mz = param_.getValue("2d:tolerance_mz");

    UInt current_charge     = 0; // charge state of the current isotopic cluster
    double mz_in_hash   = 0; // used as reference to the current isotopic peak

    // sweep through scans
    for (UInt curr_scan = 0; ms_exp_it + curr_scan != ms_exp_it_end; ++curr_scan)
    {
      Size nr_peaks_in_scan = (ms_exp_it + curr_scan)->size();
      if (nr_peaks_in_scan == 0)
        continue;

      //last_rt = current_rt;
      current_rt = (ms_exp_it + curr_scan)->getRT();
      typename PeakMap::SpectrumType::Iterator peak_it  = (ms_exp_it + curr_scan)->begin();

      // copy cluster information of least scan
      iso_last_scan = iso_curr_scan;
      iso_curr_scan.clear();
      clusters_last_scan = clusters_curr_scan;
      clusters_curr_scan.clear();

#ifdef DEBUG_2D
      std::cout << "Next scan with rt: " << current_rt << std::endl;
      std::cout << "Next scan, rt = " << current_rt << " last_rt: " << last_rt << std::endl;
      std::cout << "---------------------------------------------------------------------------" << std::endl;
#endif
      MSSpectrum s;
      s.setRT(current_rt);
      // check if there were scans in between
      if (last_rt == 0 || // are we in the first scan
          ((lower_bound(first, last, s, typename SpectrumType::RTLess()) - 1)->getRT() == last_rt))
      {

        for (UInt curr_peak = 0; curr_peak < (ms_exp_it + curr_scan)->size() - 1; ++curr_peak)
        {

          // store the m/z of the current peak
          double curr_mz         = (peak_it + curr_peak)->getMZ();
          double dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - curr_mz;

          if (dist2nextpeak <= max_peak_distance_) // one single peak without neighbors isn't optimized
          {
#ifdef DEBUG_2D
            std::cout << "Isotopic pattern found ! " << std::endl;
            std::cout << "We are at: " << (peak_it + curr_peak)->getMZ()  << " " << curr_mz << std::endl;
#endif
            if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
            {
              std::sort(iso_last_scan.begin(), iso_last_scan.end());
              // there were some isotopic clusters in the last scan...
              std::vector<double>::iterator it =
                searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);

              double delta_mz = fabs(*it - curr_mz);
              //std::cout << delta_mz << " "<< tolerance_mz << std::endl;
              if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
              {
                mz_in_hash = curr_mz; // update current hash key

                // create new isotopic cluster
// #ifdef DEBUG_2D
//                                                      std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
// #endif
                IsotopeCluster new_cluster;
                new_cluster.peaks.charge  = current_charge;
                new_cluster.scans.push_back(curr_scan);
                cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));

              }
              else
              {
// //#ifdef DEBUG_2D
//                                                      std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
// //#endif
                cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];

                // check whether this scan is already contained
                if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
                    == cluster_iter->second.scans.end())
                {
                  cluster_iter->second.scans.push_back(curr_scan);
                }

//                                                      //#ifdef DEBUG_2D
//                                                      std::cout << "Cluster with " << cluster_iter->second.peaks.size()
//                                                                          << " peaks retrieved." << std::endl;
//                                                      //#endif
              }

            }
            else                             // last scan did not contain any isotopic cluster
            {
//                                              //#ifdef DEBUG_2D
//                                              std::cout << "Last scan was empty => creating new cluster." << std::endl;
//                                              std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
//                                              //#endif

              mz_in_hash = curr_mz; // update current hash key

              // create new isotopic cluster
              IsotopeCluster new_cluster;
              new_cluster.peaks.charge  = current_charge;
              new_cluster.scans.push_back(curr_scan);
              cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));

            }

//                                      //#ifdef DEBUG_2D
//                                      std::cout << "Storing found peak in current isotopic cluster" << std::endl;
//                                      //#endif

            cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));

            iso_curr_scan.push_back(mz_in_hash);
            clusters_curr_scan.push_back(cluster_iter);
            ++curr_peak;

            cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
            iso_curr_scan.push_back((peak_it + curr_peak)->getMZ());
            clusters_curr_scan.push_back(cluster_iter);

            // check distance to next peak
            if ((curr_peak + 1) >= nr_peaks_in_scan)
              break;
            dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() -  (peak_it + curr_peak)->getMZ();

            // loop until end of isotopic pattern in this scan
            while (dist2nextpeak <= max_peak_distance_
                  &&  curr_peak < (nr_peaks_in_scan - 1))
            {
              cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak + 1)); // save peak in cluster
              iso_curr_scan.push_back((peak_it + curr_peak + 1)->getMZ());
              clusters_curr_scan.push_back(cluster_iter);
              // std::cout << "new entered: "<<(peak_it+curr_peak+1)->getMZ()<<" im while"<<std::endl;
              ++curr_peak;
              if (curr_peak >= nr_peaks_in_scan - 1)
                break;
              dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() -  (peak_it + curr_peak)->getMZ(); // get distance to next peak

            } // end while(...)

          } // end of if (dist2nextpeak <= max_peak_distance_)
          else
          {
            if (!iso_last_scan.empty()) // Did we find any isotopic cluster in the last scan?
            {
              std::sort(iso_last_scan.begin(), iso_last_scan.end());
              // there were some isotopic clusters in the last scan...
              std::vector<double>::iterator it =
                searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);

              double delta_mz = fabs(*it - curr_mz);
              // std::cout << delta_mz << " "<< tolerance_mz << std::endl;
              if (delta_mz > tolerance_mz) // check if first peak of last cluster is close enough
              {
                mz_in_hash = curr_mz; // update current hash key

                // create new isotopic cluster
//                                                      //#ifdef DEBUG_2D
//                                                      std::cout << "Last peak cluster too far, creating new cluster at "<<curr_mz << std::endl;
//                                                      //#endif
                IsotopeCluster new_cluster;
                new_cluster.peaks.charge  = current_charge;
                new_cluster.scans.push_back(curr_scan);
                cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));

              }
              else
              {
//                                                      //#ifdef DEBUG_2D
//                                                      std::cout << "Found neighbouring peak with distance (m/z) " << delta_mz << std::endl;
//                                                      //#endif
                cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];

                // check whether this scan is already contained
                if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
                    == cluster_iter->second.scans.end())
                {
                  cluster_iter->second.scans.push_back(curr_scan);
                }

//                                                      //#ifdef DEBUG_2D
//                                                      std::cout << "Cluster with " << cluster_iter->second.peaks.size()
//                                                                          << " peaks retrieved." << std::endl;
//                                                      //#endif
              }

            }
            else                             // last scan did not contain any isotopic cluster
            {
//                                              //#ifdef DEBUG_2D
//                                              std::cout << "Last scan was empty => creating new cluster." << std::endl;
//                                              std::cout << "Creating new cluster at m/z: " << curr_mz << std::endl;
//                                              //#endif

              mz_in_hash = curr_mz; // update current hash key

              // create new isotopic cluster
              IsotopeCluster new_cluster;
              new_cluster.peaks.charge  = current_charge;
              new_cluster.scans.push_back(curr_scan);
              cluster_iter = iso_map_.insert(std::pair<double, IsotopeCluster>(mz_in_hash, new_cluster));

            }

//                                      //#ifdef DEBUG_2D
//                                      std::cout << "Storing found peak in current isotopic cluster" << std::endl;
//                                      //#endif

            cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));

            iso_curr_scan.push_back(mz_in_hash);
            clusters_curr_scan.push_back(cluster_iter);

          }

          current_charge = 0; // reset charge
        } // end for (...)
      }
      last_rt = current_rt;
    }
    curr_region_ = iso_map_.begin();
#ifdef DEBUG_2D
    std::cout << iso_map_.size() << " isotopic clusters were found ! " << std::endl;
#endif

    if (real_2D_)
      optimizeRegions_(first, last, ms_exp);
    else
      optimizeRegionsScanwise_(first, last, ms_exp);
    //#undef DEBUG_2D
  }

  void TwoDOptimization::optimizeRegions_(InputSpectrumIterator& first,
                                          InputSpectrumIterator& last,
                                          PeakMap& ms_exp)
  {
    Int counter = 0;
    // go through the clusters
    for (std::multimap<double, IsotopeCluster>::iterator it = iso_map_.begin();
         it != iso_map_.end();
         ++it)
    {
#ifdef DEBUG_2D
      std::cout << "element: " << counter << std::endl;
      std::cout << "mz: " << it->first << std::endl << "rts: ";
//              for(Size i=0;i<it->second.scans.size();++i) std::cout << it->second.scans[i] << "\n";
      std::cout << std::endl << "peaks: ";
      IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
      for (; iter != it->second.peaks.end(); ++iter)
        std::cout << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;

//for(Size i=0;i<it->second.peaks.size();++i) std::cout << ms_exp[it->first].getRT() << " "<<(ms_exp[it->first][it->second]).getMZ()<<std::endl;
      std::cout << std::endl << std::endl;

#endif

      // prepare for optimization:
      // determine the matching peaks
      matching_peaks_.clear();
      findMatchingPeaks_(it, ms_exp);
      TwoDOptimization::Data twoD_data;
      twoD_data.penalties = penalties_;
      twoD_data.matching_peaks = matching_peaks_;
      // and the endpoints of each isotope pattern in the cluster
      getRegionEndpoints_(ms_exp, first, last, counter, 400, twoD_data);

      // peaks have to be stored globally
      twoD_data.iso_map_iter = it;

      twoD_data.picked_peaks = ms_exp;
      twoD_data.raw_data_first =  first;

      Size nr_diff_peaks = matching_peaks_.size();
      twoD_data.total_nr_peaks = it->second.peaks.size();

      Size nr_parameters = nr_diff_peaks * 3 + twoD_data.total_nr_peaks;

      // initialize and set parameters for optimization
      Eigen::VectorXd x_init (nr_parameters);
      x_init.setZero();

      std::map<Int, std::vector<PeakIndex> >::iterator m_peaks_it = twoD_data.matching_peaks.begin();
      Int peak_counter = 0;
      Int diff_peak_counter = 0;
      // go through the matching peaks
      for (; m_peaks_it != twoD_data.matching_peaks.end(); ++m_peaks_it)
      {
        double av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0, height;
        std::vector<PeakIndex>::iterator iter_iter = (m_peaks_it)->second.begin();
        for (; iter_iter != m_peaks_it->second.end(); ++iter_iter)
        {
          height = ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[1][(iter_iter)->peak]; //(iter_iter)->getPeak(ms_exp).getIntensity();
          avr_height += height;
          av_mz += (iter_iter)->getPeak(ms_exp).getMZ() * height;
          av_lw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[3][(iter_iter)->peak] * height; //left width
          av_rw +=    ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[4][(iter_iter)->peak] * height; //right width
          x_init(peak_counter) = height;
          ++peak_counter;
        }
        x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter) = av_mz / avr_height;
        x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter + 1) = av_lw / avr_height;
        x_init(twoD_data.total_nr_peaks + 3 * diff_peak_counter + 2) = av_rw / avr_height;
        ++diff_peak_counter;
      }

#ifdef DEBUG_2D
      std::cout << "----------------------------\n\nstart_value: " << std::endl;
      for (Size k = 0; k < start_value->size; ++k)
      {
          std::cout << x_init(k) << std::endl;
      }
#endif
      Int num_positions = 0;
      for (Size i = 0; i < twoD_data.signal2D.size(); i += 2)
      {
        num_positions += (twoD_data.signal2D[i + 1].second - twoD_data.signal2D[i].second + 1);
#ifdef DEBUG_2D
        std::cout << twoD_data.signal2D[i + 1].second << " - " << twoD_data.signal2D[i].second << " +1 " << std::endl;
#endif

      }
#ifdef DEBUG_2D
      std::cout << "num_positions : " << num_positions << std::endl;
#endif

      TwoDOptFunctor functor (nr_parameters, std::max(num_positions + 1, (Int)(nr_parameters)), &twoD_data);
      Eigen::LevenbergMarquardt<TwoDOptFunctor> lmSolver (functor);
      Eigen::LevenbergMarquardtSpace::Status status = lmSolver.minimize(x_init);

      //the states are poorly documented. after checking the source, we believe that
      //all states except NotStarted, Running and ImproperInputParameters are good
      //termination states.
      if (status <= Eigen::LevenbergMarquardtSpace::ImproperInputParameters)
      {
          throw Exception::UnableToFit(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "UnableToFit-TwoDOptimization:", "Could not fit the data: Error " + String(status));
      }

      Int peak_idx = 0;
      std::map<Int, std::vector<PeakIndex> >::iterator itv
        = twoD_data.matching_peaks.begin();
      for (; itv != twoD_data.matching_peaks.end(); ++itv)
      {
        Int i = distance(twoD_data.matching_peaks.begin(), itv);
        for (Size j = 0; j < itv->second.size(); ++j)
        {

#ifdef DEBUG_2D
          std::cout << "pos: " << itv->second[j].getPeak(ms_exp).getMZ() << "\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak] //itv->second[j].getPeak(ms_exp).getIntensity()
                    << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
                    << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";

#endif
          double mz = x_init(twoD_data.total_nr_peaks + 3 * i);
          ms_exp[itv->second[j].spectrum][itv->second[j].peak].setMZ(mz);
          double height = x_init(peak_idx);
          ms_exp[itv->second[j].spectrum].getFloatDataArrays()[1][itv->second[j].peak] = height;
          double left_width = x_init(twoD_data.total_nr_peaks + 3 * i + 1);
          ms_exp[itv->second[j].spectrum].getFloatDataArrays()[3][itv->second[j].peak] = left_width;
          double right_width = x_init(twoD_data.total_nr_peaks + 3 * i + 2);

          ms_exp[itv->second[j].spectrum].getFloatDataArrays()[4][itv->second[j].peak] = right_width;
          // calculate area
          if ((PeakShape::Type)(Int)ms_exp[itv->second[j].spectrum].getFloatDataArrays()[5][itv->second[j].peak] == PeakShape::LORENTZ_PEAK)
          {
            double x_left_endpoint = mz - 1 / left_width* sqrt(height / 1 - 1);
            double x_right_endpoint = mz + 1 / right_width* sqrt(height / 1 - 1);
            double area_left = -height / left_width* atan(left_width * (x_left_endpoint - mz));
            double area_right = -height / right_width* atan(right_width * (mz - x_right_endpoint));
            ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
          }
          else         // it's a sech peak
          {
            double x_left_endpoint = mz - 1 / left_width* boost::math::acosh(sqrt(height / 0.001));
            double x_right_endpoint = mz + 1 / right_width* boost::math::acosh(sqrt(height / 0.001));
            double area_left = -height / left_width * (sinh(left_width * (mz - x_left_endpoint)) / cosh(left_width * (mz - x_left_endpoint)));
            double area_right = -height / right_width * (sinh(right_width * (mz - x_right_endpoint)) / cosh(right_width * (mz - x_right_endpoint)));
            ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
          }

#ifdef DEBUG_2D
          std::cout << "pos: " << itv->second[j].getPeak(ms_exp).getMZ() << "\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak] //itv->second[j].getPeak(ms_exp).getIntensity()
                    << "\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
                    << "\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] << "\n";
#endif

          ++peak_idx;

        }
      }

      ++counter;
    } // end for
    //#undef DEBUG_2D
  }

  void TwoDOptimization::optimizeRegionsScanwise_(InputSpectrumIterator& first,
                                                  InputSpectrumIterator& last,
                                                  PeakMap& ms_exp)
  {
    Int counter = 0;
    TwoDOptimization::Data d;
    d.picked_peaks = ms_exp;
    d.raw_data_first =  first;

    struct OpenMS::OptimizationFunctions::PenaltyFactors penalties;

    ParamValue pv = param_.getValue("penalties:position");
    if (pv.isEmpty() || pv.toString().empty())
      penalties.pos = 0.;
    else
      penalties.pos = (float)pv;

    pv = param_.getValue("penalties:left_width");
    if (pv.isEmpty() || pv.toString().empty())
      penalties.lWidth = 1.;
    else
      penalties.lWidth = (float)pv;

    pv = param_.getValue("penalties:right_width");
    if (pv.isEmpty() || pv.toString().empty())
      penalties.rWidth = 1.;
    else
      penalties.rWidth = (float)pv;
#ifdef DEBUG_2D
    std::cout << penalties.pos << " "
              << penalties.rWidth << " "
              << penalties.lWidth << std::endl;
#endif
//      PeakMap::const_iterator help = first;
//      // std::cout << "\n\n\n\n---------------------------------------------------------------";
//      while(help!=last)
//          {
//              // std::cout<<help->getRT()<<std::endl;
//              ++help;
//          }
    // std::cout << "---------------------------------------------------------------\n\n\n\n";

    UInt max_iteration;
    pv = param_.getValue("iterations");
    if (pv.isEmpty() || pv.toString().empty())
      max_iteration = 15;
    else
      max_iteration = (UInt)pv;

    std::vector<PeakShape> peak_shapes;

    // go through the clusters
    for (std::multimap<double, IsotopeCluster>::iterator it = iso_map_.begin();
         it != iso_map_.end();
         ++it)
    {
      d.iso_map_iter = it;
#ifdef DEBUG_2D
      std::cerr << "element: " << counter << std::endl;
      std::cerr << "mz: " << it->first << std::endl << "rts: ";
      for (Size i = 0; i < it->second.scans.size(); ++i)
        std::cerr << it->second.scans[i] << "\n";
      std::cerr << std::endl << "peaks: ";
      IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
      for (; iter != it->second.peaks.end(); ++iter)
        std::cerr << ms_exp[iter->first].getRT() << " " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
      //for(Size i=0;i<it->second.peaks_.size();++i) std::cout << ms_exp[it->first].getRT() << " "<<(ms_exp[it->first][it->second]).getMZ()<<std::endl;
      std::cerr << std::endl << std::endl;

#endif
      // prepare for optimization:
      // determine the matching peaks
      // and the endpoints of each isotope pattern in the cluster

      getRegionEndpoints_(ms_exp, first, last, counter, 400, d);
      OptimizePick::Data data;

      Size idx = 0;
      for (Size i = 0; i < d.signal2D.size() / 2; ++i)
      {
        data.positions.clear();
        data.signal.clear();

        PeakMap::SpectrumType::const_iterator ms_it =
          (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i].second;
        Int size = distance(ms_it, (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second);
        data.positions.reserve(size);
        data.signal.reserve(size);

        while (ms_it != (d.raw_data_first + d.signal2D[2 * i].first)->begin() + d.signal2D[2 * i + 1].second)
        {
          data.positions.push_back(ms_it->getMZ());
          data.signal.push_back(ms_it->getIntensity());
          ++ms_it;
        }

        IsotopeCluster::IndexPair pair;
        pair.first =  d.iso_map_iter->second.peaks.begin()->first + idx;

        IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
                                                                        d.iso_map_iter->second.peaks.end(),
                                                                        pair, [](auto& left, auto& right){return left.first < right.first;});

        // find the last entry with this rt-value
        ++pair.first;
        IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(d.iso_map_iter->second.peaks.begin(),
                                                                         d.iso_map_iter->second.peaks.end(),
                                                                         pair, [](auto& left, auto& right){return left.first < right.first;});

        while (set_iter != set_iter2)
        {
          const Size peak_index = set_iter->second;
          const MSSpectrum& spec = ms_exp[set_iter->first];
          PeakShape shape(spec.getFloatDataArrays()[1][peak_index], //intensity
                          spec[peak_index].getMZ(),
                          spec.getFloatDataArrays()[3][peak_index], //left width
                          spec.getFloatDataArrays()[4][peak_index], //right width
                          spec[peak_index].getIntensity(), //area is stored in peak intensity
                          PeakShape::Type(Int(spec.getFloatDataArrays()[5][peak_index]))); //shape
          peak_shapes.push_back(shape);
          ++set_iter;
        }
#ifdef DEBUG_2D
        std::cout << "rt "
                  << (d.raw_data_first + d.signal2D[2 * i].first)->getRT()
                  << "\n";
#endif
        OptimizePick opt(penalties, max_iteration);
#ifdef DEBUG_2D
        std::cout << "vorher\n";

        for (Size p = 0; p < peak_shapes.size(); ++p)
        {
          std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
                    << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width  << std::endl;
        }
#endif
        opt.optimize(peak_shapes, data);
#ifdef DEBUG_2D
        std::cout << "nachher\n";
        for (Size p = 0; p < peak_shapes.size(); ++p)
        {
          std::cout << peak_shapes[p].mz_position << "\t" << peak_shapes[p].height
                    << "\t" << peak_shapes[p].left_width << "\t" << peak_shapes[p].right_width  << std::endl;
        }
#endif
        std::sort(peak_shapes.begin(), peak_shapes.end(), PeakShape::PositionLess());
        pair.first =  d.iso_map_iter->second.peaks.begin()->first + idx;

        set_iter = lower_bound(d.iso_map_iter->second.peaks.begin(),
                               d.iso_map_iter->second.peaks.end(),
                               pair, [](auto& left, auto& right){return left.first < right.first;});
        Size p = 0;
        while (p < peak_shapes.size())
        {
          MSSpectrum& spec = ms_exp[set_iter->first];
          spec[set_iter->second].setMZ(peak_shapes[p].mz_position);
          spec.getFloatDataArrays()[3][set_iter->second] = peak_shapes[p].left_width;
          spec.getFloatDataArrays()[4][set_iter->second] = peak_shapes[p].right_width;
          spec.getFloatDataArrays()[1][set_iter->second] = peak_shapes[p].height; // maximum intensity
          // calculate area
          if (peak_shapes[p].type == PeakShape::LORENTZ_PEAK)
          {
            PeakShape& ps = peak_shapes[p];
            double x_left_endpoint = ps.mz_position - 1 / ps.left_width* sqrt(ps.height / 1 - 1);
            double x_right_endpoint = ps.mz_position + 1 / ps.right_width* sqrt(ps.height / 1 - 1);
            double area_left = -ps.height / ps.left_width* atan(ps.left_width * (x_left_endpoint - ps.mz_position));
            double area_right = -ps.height / ps.right_width* atan(ps.right_width * (ps.mz_position - x_right_endpoint));
            spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
          }
          else        //It's a Sech - Peak
          {
            PeakShape& ps = peak_shapes[p];
            double x_left_endpoint = ps.mz_position - 1 / ps.left_width* boost::math::acosh(sqrt(ps.height / 0.001));
            double x_right_endpoint = ps.mz_position + 1 / ps.right_width* boost::math::acosh(sqrt(ps.height / 0.001));
            double area_left = ps.height / ps.left_width * (sinh(ps.left_width * (ps.mz_position - x_left_endpoint)) / cosh(ps.left_width * (ps.mz_position - x_left_endpoint)));
            double area_right = -ps.height / ps.right_width * (sinh(ps.right_width * (ps.mz_position - x_right_endpoint)) / cosh(ps.right_width * (ps.mz_position - x_right_endpoint)));
            spec[set_iter->second].setIntensity(area_left + area_right); // area is stored as peak intensity
          }
          ++set_iter;
          ++p;
        }
        ++idx;
        peak_shapes.clear();
      }

      ++counter;
    }
  }

  void TwoDOptimization::getRegionEndpoints_(PeakMap& exp,
                                             InputSpectrumIterator& first,
                                             InputSpectrumIterator& last,
                                             Size iso_map_idx,
                                             double noise_level,
                                             TwoDOptimization::Data& d)
  {
    d.signal2D.clear();
    typedef typename InputSpectrumIterator::value_type InputExperimentType;
    typedef typename InputExperimentType::value_type InputPeakType;
    typedef std::multimap<double, IsotopeCluster> MapType;

    double rt, first_peak_mz, last_peak_mz;

    typename PeakMap::SpectrumType spec;
    InputPeakType peak;

    MapType::iterator iso_map_iter = iso_map_.begin();
    for (Size i = 0; i < iso_map_idx; ++i)
      ++iso_map_iter;

#ifdef DEBUG2D
    std::cout << "rt begin: " << exp[iso_map_iter->second.scans[0]].getRT()
              << "\trt end: " << exp[iso_map_iter->second.scans[iso_map_iter->second.scans.size() - 1]].getRT()
              << " \t" << iso_map_iter->second.scans.size() << " scans"
              << std::endl;
#endif

    // get left and right endpoint for all scans in the current cluster
    for (Size i = 0; i < iso_map_iter->second.scans.size(); ++i)
    {
      typename PeakMap::iterator exp_it;

      // first the right scan through binary search
      rt = exp[iso_map_iter->second.scans[i]].getRT();
      spec.setRT(rt);
      InputSpectrumIterator iter = lower_bound(first, last, spec, MSSpectrum::RTLess());
      //  if(iter->getRT() != rt) --iter;
      exp_it = exp.RTBegin(rt);
#ifdef DEBUG2D
      std::cout << exp_it->getRT() << " vs " << iter->getRT() << std::endl;
#endif
      // now the right mz
      IsotopeCluster::IndexPair pair;
      pair.first =  iso_map_iter->second.peaks.begin()->first + i;
      // get iterator in peaks-set that points to the first peak in the current scan
      IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(iso_map_iter->second.peaks.begin(),
                                                                      iso_map_iter->second.peaks.end(),
                                                                      pair, [](auto& left, auto& right){return left.first < right.first;});

      // consider a bit more of the signal to the left
      first_peak_mz = (exp_it->begin() + set_iter->second)->getMZ() - 1;

      // find the last entry with this rt-value
      ++pair.first;
      IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(iso_map_iter->second.peaks.begin(),
                                                                       iso_map_iter->second.peaks.end(),
                                                                       pair, [](auto& left, auto& right){return left.first < right.first;});

      if (i == iso_map_iter->second.scans.size() - 1)
      {
        set_iter2 = iso_map_iter->second.peaks.end();
        --set_iter2;
      }
      else if (set_iter2 != iso_map_iter->second.peaks.begin())
        --set_iter2;

      last_peak_mz = (exp_it->begin() + set_iter2->second)->getMZ() + 1;

      //std::cout << rt<<": first peak mz "<<first_peak_mz << "\tlast peak mz "<<last_peak_mz <<std::endl;
      peak.setPosition(first_peak_mz);
      typename PeakMap::SpectrumType::const_iterator raw_data_iter
        = lower_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
      if (raw_data_iter != iter->begin())
      {
        --raw_data_iter;
      }
      double intensity = raw_data_iter->getIntensity();
      // while the intensity is falling go to the left
      while (raw_data_iter != iter->begin() && (raw_data_iter - 1)->getIntensity() < intensity &&
             (raw_data_iter - 1)->getIntensity() > noise_level)
      {
        --raw_data_iter;
        intensity = raw_data_iter->getIntensity();
      }
      ++raw_data_iter;
      IsotopeCluster::IndexPair left, right;
      left.first = distance(first, iter);
      left.second = raw_data_iter - iter->begin();
#ifdef DEBUG2D
      std::cout << "left: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
#endif
      // consider a bit more of the signal to the right
      peak.setPosition(last_peak_mz + 1);
      raw_data_iter
        = upper_bound(iter->begin(), iter->end(), peak, typename InputPeakType::PositionLess());
      if (raw_data_iter == iter->end())
        --raw_data_iter;
      intensity = raw_data_iter->getIntensity();
      // while the intensity is falling go to the right
      while (raw_data_iter + 1 != iter->end() && (raw_data_iter + 1)->getIntensity() < intensity)
      {
        ++raw_data_iter;
        intensity = raw_data_iter->getIntensity();
        if ((raw_data_iter + 1 != iter->end()) && (raw_data_iter + 1)->getIntensity() > noise_level)
          break;
      }
      right.first = left.first;
      right.second = raw_data_iter - iter->begin();
#ifdef DEBUG2D
      std::cout << "right: " << iter->getRT() << "\t" << raw_data_iter->getMZ() << std::endl;
#endif
      // region endpoints are stored in global vector
      d.signal2D.push_back(left);
      d.signal2D.push_back(right);
    }
#ifdef DEBUG2D
    std::cout << first_peak_mz << "\t" << last_peak_mz << std::endl;
#endif
  }

}

ext: cpp

max_stars_repo_name: proteomicsunitcrg/OpenMS

max_issues_repo_name: proteomicsunitcrg/OpenMS

max_forks_repo_name: proteomicsunitcrg/OpenMS


text:
#ifdef __linux__

#include "generator.linux.hpp"

#include <boost/algorithm/string.hpp>
#include <boost/variant/apply_visitor.hpp>
#include <boost/variant/variant.hpp>

#include "blackhole/extensions/format.hpp"

namespace blackhole {
inline namespace v1 {
namespace datetime {
namespace {

struct visitor_t : public boost::static_visitor<> {
    fmt::MemoryWriter& stream;
    std::tm tm;
    std::uint64_t usec;
    char buffer[1024];

    visitor_t(fmt::MemoryWriter& stream, const std::tm& tm, std::uint64_t usec) :
        stream(stream),
        tm(tm),
        usec(usec)
    {}

    auto operator()(const literal_t& value) -> void {
        std::size_t ret = std::strftime(buffer, sizeof(buffer), value.value.c_str(), &tm);
        stream << fmt::StringRef(buffer, ret);
    }

    auto operator()(epoch_t) -> void {
        stream.write("{}", std::mktime(&tm));
    }

    auto operator()(usecond_t) -> void {
        stream.write("{:06d}", usec);
    }
};

}  // namespace

generator_t::generator_t(std::string pattern) {
    std::string literal;

    auto pos = std::begin(pattern);
    while (pos != std::end(pattern)) {
        if (boost::starts_with(boost::make_iterator_range(pos, std::end(pattern)), "%f")) {
            tokens.emplace_back(literal_t{std::move(literal)});
            tokens.emplace_back(usecond_t{});
            literal.clear();
            ++pos;
            ++pos;
        } else if (boost::starts_with(boost::make_iterator_range(pos, std::end(pattern)), "%s")) {
            tokens.emplace_back(literal_t{std::move(literal)});
            tokens.emplace_back(epoch_t{});
            literal.clear();
            ++pos;
            ++pos;
        } else {
            literal.push_back(*pos);
            ++pos;
        }
    }

    if (!literal.empty()) {
        tokens.emplace_back(literal_t{std::move(literal)});
    }
}

auto generator_t::operator()(writer_type& stream, const std::tm& tm, std::uint64_t usec) const -> void {
    visitor_t visitor(stream, tm, usec);

    for (auto& token : tokens) {
        boost::apply_visitor(visitor, token);
    }
}

auto make_generator(const std::string& pattern) -> generator_t {
    return generator_t(pattern);
}

} // namespace datetime
} // namespace v1
} // namespace blackhole

#endif

max_stars_repo_name: JakariaBlaine/blackhole

max_issues_repo_name: JakariaBlaine/blackhole

max_forks_repo_name: JakariaBlaine/blackhole


text:
/*
 [auto_generated]
 boost/numeric/odeint/external/openmp/openmp_state.hpp

 [begin_description]
 Wrappers for OpenMP.
 [end_description]

 Copyright 2013 Karsten Ahnert
 Copyright 2013 Mario Mulansky
 Copyright 2013 Pascal Germroth

 Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or
 copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_STATE_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_STATE_HPP_INCLUDED

#include <omp.h>
#include <vector>
#include <algorithm>
#include <boost/range/adaptor/sliced.hpp>
#include <boost/numeric/odeint/util/copy.hpp>
#include <boost/numeric/odeint/util/split.hpp>
#include <boost/numeric/odeint/util/resize.hpp>
#include <boost/numeric/odeint/external/openmp/openmp_nested_algebra.hpp>

namespace boost {
namespace numeric {
namespace odeint {

/** \brief A container that is split into distinct parts, for threading.
 * Just a wrapper for vector<vector<T>>, use `copy` for splitting/joining.
 */
template< class T >
struct openmp_state : public std::vector< std::vector< T > >
{
    openmp_state() {}

    openmp_state(size_t n, const std::vector<T>& val = std::vector<T>())
    : std::vector< std::vector< T > >(n, val) {}

    template<class InputIterator>
    openmp_state(InputIterator first, InputIterator last)
    : std::vector< std::vector< T > >(first, last) {}

    openmp_state(const std::vector< std::vector< T > > &orig)
    : std::vector< std::vector< T > >(orig) {}

};

template< class T >
struct is_resizeable< openmp_state< T > > : boost::true_type { };

template< class T >
struct same_size_impl< openmp_state< T > , openmp_state< T > >
{
    static bool same_size( const openmp_state< T > &x , const openmp_state< T > &y )
    {
        if( x.size() != y.size() ) return false;
        for( size_t i = 0 ; i != x.size() ; i++ )
            if( x[i].size() != y[i].size() ) return false;
        return true;
    }
};

template< class T >
struct resize_impl< openmp_state< T > , openmp_state< T > >
{
    static void resize( openmp_state< T > &x , const openmp_state< T > &y )
    {
        x.resize( y.size() );
#       pragma omp parallel for schedule(dynamic)
        for(size_t i = 0 ; i < x.size() ; i++)
            x[i].resize( y[i].size() );
    }
};

/** \brief Copy data between openmp_states of same size. */
template< class T >
struct copy_impl< openmp_state< T >, openmp_state< T > >
{
    static void copy( const openmp_state< T > &from, openmp_state< T > &to )
    {
#       pragma omp parallel for schedule(dynamic)
        for(size_t i = 0 ; i < from.size() ; i++)
            std::copy( from[i].begin() , from[i].end() , to.begin() );
    }
};

/** \brief Copy data from some container to an openmp_state and resize it.
 * Target container size will determine number of blocks to split into.
 * If it is empty, it will be resized to the maximum number of OpenMP threads.
 * SourceContainer must support `s::value_type`, `s::const_iterator`, `s.begin()`, `s.end()` and `s.size()`,
 * with Random Access Iterators; i.e. it must be a Random Access Container. */
template< class SourceContainer >
struct split_impl< SourceContainer, openmp_state< typename SourceContainer::value_type > >
{
    static void split( const SourceContainer &from, openmp_state< typename SourceContainer::value_type > &to )
    {
        if(to.size() == 0) to.resize( omp_get_max_threads() );
        const size_t part = from.size() / to.size();
#       pragma omp parallel for schedule(dynamic)
        for(size_t i = 0 ; i < to.size() ; i++) {
            typedef typename SourceContainer::const_iterator it_t;
            const it_t begin = from.begin() + i * part;
            it_t end = begin + part;
            // for cases where from.size() % to.size() > 0
            if(i + 1 == to.size() || end > from.end()) end = from.end();
            to[i].resize(end - begin);
            std::copy(begin, end, to[i].begin());
        }
    }
};

/** \brief Copy data from an openmp_state to some container and resize it.
 * TargetContainer must support `s::value_type`, `s::iterator`, `s.begin()` and `s.resize(n)`,
 * i.e. it must be a `std::vector`. */
template< class TargetContainer >
struct unsplit_impl< openmp_state< typename TargetContainer::value_type >, TargetContainer >
{
    static void unsplit( const openmp_state< typename TargetContainer::value_type > &from , TargetContainer &to )
    {
        // resize target
        size_t total_size = 0;
        for(size_t i = 0 ; i < from.size() ; i++)
            total_size += from[i].size();
        to.resize( total_size );
        // copy parts
        typename TargetContainer::iterator out = to.begin();
        for(size_t i = 0 ; i < from.size() ; i++)
            out = std::copy(from[i].begin(), from[i].end(), out);
    }
};

/** \brief OpenMP-parallelized algebra.
 * For use with openmp_state.
 */
typedef openmp_nested_algebra< range_algebra > openmp_algebra;

/** \brief Use `openmp_algebra` for `openmp_state`. */
template< class T >
struct algebra_dispatcher< openmp_state< T > >
{
    typedef openmp_algebra algebra_type;
};

}
}
}

#endif

ext: hpp

max_stars_repo_name: rudylee/expo

max_issues_repo_name: rudylee/expo

max_forks_repo_name: rudylee/expo


text:
#include <iostream>
#include <numeric>

#include <boost/range/irange.hpp>

template<typename T = std::uint64_t>
constexpr auto get_code() {

	const T row = 3010;
	const T col = 3019;

	const T init_value = 20151125;
	const T multiplier = 252533;
	const T modulus    = 33554393;

	// the number at (row,col) represents the position
	// in the sequence of natural numbers that form
	// a triangle, as seen in the readme.txt
	const auto n = (row + col - 1);

	// this is the final number in the triangular arrangement
	// that the number at (row,col) is part of
	const auto triangular_number = ((n * (n + 1)) / 2);

	// target represents the end of the sequence we need to
	// go through in order to compute the value at (row,col).
	// this corresponds to the following mapping (see example
	// table in the readme.txt):
	// (x,y) => triangular(1,(x+y-1)) - x + 1
	const auto target = (triangular_number - row + 1);

	auto result = init_value;

	for(const auto i : boost::irange({1}, target)) {
		result *= multiplier;
		result %= modulus;
	}

	return result;
}

int main() {

	std::cout << get_code() << std::endl;

	return 0;
}

max_stars_repo_path: Day 25/main.cpp

max_issues_repo_path: Day 25/main.cpp

max_forks_repo_path: Day 25/main.cpp


text: #include <boost/numeric/odeint/util/is_resizeable.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschrÃ¤nkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
// This file is part of the Matrix Template Library
//
// See also license.mtl.txt in the distribution.

#include <iostream>
#include <vector>
#include <complex>
#include <boost/numeric/mtl/mtl.hpp>

using namespace std;
using namespace mtl;

int main(int, char**)
{
    float A[2][3]= {{2, 3, 6}, {1, 2, 0}};
    A[0][0]= 24.;

    std::complex<float> z(2, 5);

    std::vector<std::complex<float> > v(3, z);

    std::vector<float> v2(3, 2.5);

    dense_vector<double>   w(3), w2(9);
    w= 4, 5, 6;
    iota(w2);

    dense2D<double>  B(2, 3);
    B= 2;

    dense2D<double, mat::parameters<col_major> >  C(2, 3);
    C= 2, 3, 4,
       5, 6, 7;

    cout << "C is\n" << C;

    compressed2D<double>  D(2, 3);
    D= 2;

    compressed2D<double, mat::parameters<col_major> >  E(2, 3);
    E= 2, 3, 4,
       5, 6, 7;

    cout << "E is\n" << E;

    compressed2D<float> F;
    laplacian_setup(F, 3, 4);
    cout << "F is\n" << F;

    return 0;
}


text:
/*###############################################################################
#
# This software is distributed under the following ISC license:
#
# Copyright © 2017 BRAZILIAN SYNCHROTRON LIGHT SOURCE <sol@lnls.br>
#   Dawood Alnajjar <dawood.alnajjar@lnls.br>
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
# OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
# CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#
# Description:
# Analog output functions
#
###############################################################################*/

#include "Common.h"

#include <ctype.h>
#include <errno.h>
#include <stdexcept>
#include <time.h>
#include <boost/bimap.hpp>
#include "rt_var_handler.h"
#include "ANALOG.h"
#include "CrioLinux.h"
#include "utils.h"
/* boost namespace */
using boost::bimap;

/* typedefs */
typedef bimap< std::string, unsigned > bm_address_type;

/* ------------ Helpers -------------- */

static __inline__ int crioSetAO(struct crio_context *ctx, uint32_t address, float value);
static __inline__ int crioSetAOFixedPoint(struct crio_context *ctx, double value, struct fxp_ctx *fxp_data) ;

static __inline__ int crioSetAOFixedPoint(struct crio_context *ctx, double value, struct fxp_ctx *fxp_data) {
    uint64_t value_u64 = dbl_to_fxp(value, fxp_data);
    auto Res = NiFpga_WriteU64(NiFpga_Session(ctx->session), fxp_data->address, value_u64);
    if (NiFpga_IsError(Res)) throw (CrioLibException(E_VAR_ACCESS, "Cannot access address"));
    return 0;
}

static __inline__ int crioSetAO(struct crio_context *ctx, uint32_t address, float value) {
    auto Res = NiFpga_WriteSgl(NiFpga_Session(ctx->session), address, value);
    if (NiFpga_IsError(Res)) throw (CrioLibException(E_VAR_ACCESS, "Cannot access address"));
    return 0;
}

/* ------------ API FUNCTIONS -------------- */
int crioGetAOArraySize(struct crio_context *ctx, unsigned *size) {
    if (!ctx->session_open)
        throw (CrioLibException(E_SESSION_CLOSED , "[%s] Operation performed on closed session.", LIB_CRIO_LINUX ));
    *size = ctx->ao_count;
    return 0;
}

int crioSetAOItem(struct crio_context *ctx, const char *name, double value) {
    struct fxp_ctx * local_fxp_data;
    if (!ctx->session_open)
        throw (CrioLibException(E_SESSION_CLOSED , "[%s] Operation performed on closed session.", LIB_CRIO_LINUX ));
    try {
        /* Check if variable is an RT variable */
        if (is_rt_var(name) == true)
        {
            unsigned address = ctx->rt_variable_offsets[((bm_address_type *)ctx->rt_addresses)->left.at(name)];
            set_rt_val(ctx->shared_memory, address, value, name);
            if (ctx->debugCRIO)
            {
                printf ("RT AO name=%s, Value=%f, Index=%d\n" , name, value, address);
                fprintf (ctx->log, "RT AO name=%s, Value=%f, Index=%d\n" , name, value, address);
            }
        }
        else
        {
            /* Check if variable is an FPGA fixedpoint variable */
            if (is_fixed_point(name) == true)
            {
                local_fxp_data = &(((struct fxp_ctx * )ctx->fxps)[((bm_address_type *)ctx->ao_addresses)->left.at(name)]);
                crioSetAOFixedPoint(ctx, value, local_fxp_data);
                if (ctx->debugCRIO)
                {
                    printf ("FPGA FXP AO name=%s, Value=%f, Address=0x%05x, Integer word Length=%u, Word Length=%u, Sign=%d\n" , name, value, local_fxp_data->address, local_fxp_data->int_word_length, local_fxp_data->word_length, local_fxp_data->sign);
                    fprintf (ctx->log, "FPGA FXP AO name=%s, Value=%f, Address=0x%05x, Integer word Length=%u, Word Length=%u, Sign=%d\n" , name, value, local_fxp_data->address, local_fxp_data->int_word_length, local_fxp_data->word_length, local_fxp_data->sign);
                }
            }
            /* FPGA floating-point variable */
            else
            {
                unsigned address = ((bm_address_type *)ctx->ao_addresses)->left.at(name);
                crioSetAO(ctx, address, (float)value);
                if (ctx->debugCRIO)
                {
                    printf ("FPGA SGL AO name=%s, Value=%f, Address=0x%05x\n" , name, value, address);
                    fprintf (ctx->log, "FPGA SGL AO name=%s, Value=%f, Address=0x%05x\n" , name, value, address);
                }
            }
        }
    } catch (out_of_range) {
        throw (CrioLibException(E_OUT_OF_RANGE , "[%s] Property [%s]: Query returned null.", LIB_CRIO_LINUX , name ));
    } catch(CrioLibException &e) {
        throw (CrioLibException(e.errorcode, "[%s] Property [%s]: %s.", LIB_CRIO_LINUX , name, e.what()));
    }
    return 0;
}

ext: cc

max_stars_repo_name: lnls-sol/nheengatu-linuxlibs

max_issues_repo_name: lnls-sol/nheengatu-linuxlibs

max_forks_repo_name: lnls-sol/nheengatu-linuxlibs


text:
#include <iostream>
#include <string>
#include <boost/asio.hpp>
#include <boost/bind.hpp>
#include <boost/thread/thread.hpp>
#include <boost/array.hpp>

#define DURATION_MILLIS_COUNT(time_point) std::chrono::duration_cast<std::chrono::milliseconds> (time_point).count()

using namespace boost;
using namespace boost::asio;

bool isServerMode = false;
int port = -1;
std::string ipAddr = "";

void printHelp() {
    std::cout << "usage: mnc ([-l listen/server mode]) [destination multicast ip address] [port]" << std::endl;
}

bool processArgs(int argc, char *argv[]) {

    if(argc < 3) {
        std::cout << "not enough arguments" << std::endl;
        printHelp();
        return false;
    }

    char option;
    while ( (option = getopt(argc, argv,":l")) != -1 ) {
        switch(option) {
            case 'l':
                isServerMode = true;
                break;
            case '?':
                printf("Unknown option: %c\n", optopt);
                printHelp();
                return false;
                break;
        }
    }

    ipAddr = std::string(argv[argc - 2]);
    port =  std::atoi(argv[argc - 1]);

    if(port == -1) {
        std::cout << "missing port" << std::endl;
        printHelp();
        return false;
    }
    //if(ipAddr == "localhost") ipAddr = "127.0.0.1";
    if(ipAddr == "") {
        std::cout << "missing destination ip addr" << std::endl;
        printHelp();
        return false;
    }
    return true;
}

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

    const auto period = std::chrono::milliseconds(10);
    auto t0 = std::chrono::steady_clock::now();

    if(!processArgs(argc, argv)) return -1;

    io_service ios;
    ip::udp::endpoint ep(ip::address::from_string(ipAddr), port);

    auto receive_buffer = boost::array<char, 1024>();
    auto socket = ip::udp::socket(ios);

    socket.open(ep.protocol());
    socket.set_option(ip::udp::socket::reuse_address(true));
    socket.bind(ep);
    socket.set_option(ip::multicast::join_group(ep.address()));

    size_t packetSize;
    std::string packetString;

    while(true) {
        auto t = std::chrono::steady_clock::now();

//        std::cout << DURATION_MILLIS_COUNT(std::chrono::steady_clock::now() - t0) << std::endl;
//        std::this_thread::sleep_for(std::chrono::milliseconds(500));

        if(isServerMode) {
            packetSize = socket.receive_from(asio::buffer(receive_buffer), ep);
            packetString = std::string(receive_buffer.begin(), receive_buffer.begin() + packetSize);
            std::cout << packetString << std::flush;
        } else {
            std::getline(std::cin, packetString);
            socket.send_to(asio::buffer(packetString + "\n"), ep);
        }

        std::this_thread::sleep_until(t + period);
    }

    return 0;
}