ProductModel.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: $
// --------------------------------------------------------------------------
 
 
#pragma once
 
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/BaseModel.h>
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/ModelDescription.h>
#include <OpenMS/KERNEL/Peak2D.h>
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/IsotopeModel.h>
#include <OpenMS/CHEMISTRY/EmpiricalFormula.h>
 
namespace OpenMS
{
 
  template <UInt D>
  class ProductModel;
 
  template <>
  class OPENMS_DLLAPI ProductModel<2>:
    public BaseModel<2>
  {
public:
 
    enum {D = 2};
 
    typedef double IntensityType;
    typedef DPosition<D> PositionType;
    typedef BaseModel<D>::SamplesType SamplesType;
 
    ProductModel() :
      BaseModel<D>(),
      distributions_(D, nullptr)
    {
      this->setName(this->getProductName());
 
      //Register model info
      for (UInt dim = 0; dim < D; ++dim)
      {
        String name = Peak2D::shortDimensionName(dim);
        this->subsections_.push_back(name);
        this->defaults_.setValue(name, "GaussModel", "Name of the model used for this dimension");
      }
 
      //defaults
      this->defaults_.setValue("intensity_scaling", 1.0, "Scaling factor used to adjust the model distribution to the intensities of the data");
      this->defaultsToParam_();
    }
 
    ProductModel(const ProductModel & source) :
      BaseModel<D>(source),
      distributions_(D, nullptr),
      scale_(source.scale_)
    {
      for (UInt dim = 0; dim < D; ++dim)
      {
        // clone source model
        if (source.distributions_[dim])
        {
          ModelDescription<1> desc(source.distributions_[dim]);
          setModel(dim, desc.createModel());
        }
      }
      updateMembers_();
    }
 
    ~ProductModel() override
    {
      for (Size dim = 0; dim < D; ++dim)
      {
        delete distributions_[dim];
      }
    }
 
    virtual ProductModel & operator=(const ProductModel & source)
    {
      if (&source == this) return *this;
 
      BaseModel<D>::operator=(source);
      scale_ = source.scale_;
 
      for (UInt dim = 0; dim < D; ++dim)
      {
        if (source.distributions_[dim])
        {
          // clone source model
          ModelDescription<1> desc(source.distributions_[dim]);
          setModel(dim, desc.createModel());
        }
        else
        {
          distributions_[dim] = nullptr;
        }
      }
      updateMembers_();
 
      return *this;
    }
 
    IntensityType getIntensity(const PositionType & pos) const override
    {
      IntensityType intens(scale_);
      for (UInt dim = 0; dim < D; ++dim)
      {
        if (distributions_[dim] == nullptr)
        {
          throw Exception::BaseException(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("ProductModel: model for dimension ") + dim + " not set.", "");
        }
        intens *= distributions_[dim]->getIntensity(pos[dim]);
      }
      return intens;
    }
 
    static BaseModel<D> * create()
    {
      return new ProductModel<D>();
    }
 
    static const String getProductName()
    {
      return String("ProductModel") + D + "D";
    }
 
    ProductModel & setModel(UInt dim, BaseModel<1> * dist)
    {
      OPENMS_PRECONDITION(dim < D, "ProductModel<D>:getModel(Position): index overflow!");
      if (dist == nullptr || dist == distributions_[dim])
      {
        return *this;
      }
 
      delete distributions_[dim];
      distributions_[dim] = dist;
 
      // Update model info
      String name = Peak2D::shortDimensionName(dim);
      this->param_.removeAll(name + ':');
      this->param_.insert(name + ':', distributions_[dim]->getParameters());
      this->param_.setValue(name, distributions_[dim]->getName());
 
      return *this;
    }
 
    BaseModel<1> * getModel(UInt dim) const
    {
      OPENMS_PRECONDITION(dim < D, "ProductModel<D>:getModel(Position): index overflow!");
      return distributions_[dim];
    }
 
    IntensityType getScale() const
    {
      return scale_;
    }
 
    void setScale(IntensityType scale)
    {
      this->setCutOff(this->getCutOff() / scale_);  // remove scaling from cutoff
      scale_ = scale;
      this->param_.setValue("intensity_scaling", scale);
      this->setCutOff(this->getCutOff() * scale_);  // scale cutoff
    }
 
    void getSamples(SamplesType & cont) const override
    {
      cont.clear();
      typedef BaseModel<1>::SamplesType Samples1D;
      std::vector<Samples1D> samples(D);
      // get samples for each dimension
      for (Size dim = 0; dim < D; ++dim)
      {
        distributions_[dim]->getSamples(samples[dim]);
      }
 
      BaseModel<D>::PeakType peak;
      std::vector<UInt> i(D, 0);  // index vector
 
      while (i[D - 1] < samples[D - 1].size())
      {
        for (UInt dim = 0; dim < D; ++dim)
        {
          peak.getPosition()[dim] = samples[dim][i[dim]].getPosition()[0];
        }
        fillIntensity(peak);
        cont.push_back(peak);
 
        ++i[0];
        for (Size dim = 0; dim < D - 1; ++dim)
        {
          if (i[dim] >= samples[dim].size())
          {
            i[dim] = 0;
            ++i[dim + 1];
          }
        }
      }
      return;
    }
 
protected:
    void updateMembers_() override
    {
      BaseModel<D>::updateMembers_();
      scale_ = (double)(this->param_.getValue("intensity_scaling"));
      for (UInt dim = 0; dim < D; ++dim)
      {
        String name = Peak2D::shortDimensionName(dim);
        if (this->param_.exists(name))
        {
          delete distributions_[dim];
          distributions_[dim] = Factory<BaseModel<1> >::create(this->param_.getValue(name).toString());
          Param copy = this->param_.copy(name + ":", true);
          distributions_[dim]->setParameters(copy);
          if (distributions_[dim]->getName().hasSubstring("IsotopeModel"))
          {
            static_cast<IsotopeModel *>(distributions_[dim])->setSamples(static_cast<IsotopeModel *>(distributions_[dim])->getFormula());
          }
        }
      }
    }
 
    std::vector<BaseModel<1> *> distributions_;
    IntensityType scale_;
  };
}