fitExpModGaussian.hpp

Go to the documentation of this file.

/** \file fitExpModGaussian.hpp
 * \author Jared R. Males
 * \brief Tools for fitting the GammaDistribution distribution to data.
 * \ingroup fitting_files
 *
 */
 
//***********************************************************************//
// Copyright 2023 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef fitExpModGaussian_hpp
#define fitExpModGaussian_hpp
 
#include "levmarInterface.hpp"
#include "../func/expModGaussian.hpp"
 
namespace mx
{
namespace math
{
namespace fit
{
 
template<typename realT>
struct array2FitExpModGaussian;
 
 
/** \defgroup expModGauss_peak_fit Exponetially Modified Gaussian Distribution
  * \brief Fitting the Exponetially Modified Gaussian Distribution to data.
  * 
  * The Exponetially Modified Gaussian Distribution is fit to data.
  * 
  * \ingroup peak_fit 
  */
 
///Class to manage fitting the Exponetially Modified Gaussian Distribution  to data via the \ref levmarInterface
/** In addition to the requirements on fitterT specified by \ref levmarInterface
  * this class also requires this definition in fitterT
  * \code
  * static const int nparams = 3;
  * \endcode
  * where the number 3 is replaced by the number of parameters that fitterT expects to fit.
  *
  * \tparam fitterT a type meeting the above requirements.
  *
  * \ingroup expModGauss_peak_fit
  *
  */
template<typename fitterT>
class fitExpModGaussian : public levmarInterface<fitterT>
{
   
public:
   
   typedef typename fitterT::realT realT;
 
   static const int nparams = fitterT::nparams;
 
   array2FitExpModGaussian<realT> arr;
   
   void initialize()
   {
      this->allocate_params(nparams);
      this->adata = &arr;      
   }
   
   fitExpModGaussian()
   {
      initialize();
   }
      
   ~fitExpModGaussian()
   {
   }
   
   ///Set the initial guess.
   void setGuess( realT G0,    ///< [in] the location parameter
                  realT mu,     ///< [in] the shape parameter
                  realT sigma, ///< [in] the scale parameter
                  realT lambda  ///< [in] the denominator or 1/peak-scale
                )
   {
      static_assert( nparams==4 , "fitExpModGaussian: Wrong setGuess called for no location parameter.");
 
      this->p[3] = G0; //this is p[2] to make it easy to leave out
      this->p[0] = mu;
      this->p[1] = sigma;
      this->p[2] = lambda;
   }
 
   ///Set the initial guess.
   void setGuess( realT mu,   ///< [in] the location parameter
                  realT sigma,    ///< [in] the shape parameter
                  realT lambda ///< [in] the scale parameter
                )
   {
      static_assert( nparams==3 , "fitExpModGaussian: Wrong setGuess called for no location parameter.");
      
      this->p[0] = mu; //this is p[2] to make it easy to leave out
      this->p[1] = sigma;
      this->p[2] = lambda;
   }
 
   void setArray(realT *data, int n)
   {
      arr.data = data;
      arr.n = n;
 
      this->n = n;
      
   }
   
   void G0( realT nG0 )
   {
      arr.G0 = nG0;
      if(nparams == 4)
      {
         this->p[3] = nG0;  //this is p[2] to make it easy to leave out
      }
   }
   
   void mu( realT nmu )
   {
      arr.mu = nmu;
      this->p[0] = nmu;
   }
   
   void sigma( realT ns )
   {
      arr.sigma = ns;
      this->p[1] = ns;
   }
 
   void lambda( realT nl)
   {
      arr.lambda = nl;
      this->p[2] = nl;
   }
 
   int fit()
   {
      return levmarInterface<fitterT>::fit();      
   }
      
   realT G0()
   {
      if(nparams == 4)
      {
         return this->p[3];
      }
      else
      {
         return 0;
      }
   }
   
   realT mu()
   {
      return this->p[0];
   }
   
   realT sigma()
   {
      return this->p[1];
   }
   
   realT lambda()
   {
      return this->p[2];
   }
 
   
};
 
 
///Wrapper for a native array to pass to \ref levmarInterface, with Exponentially Modified Gaussian details.
/** \ingroup gammaDist_peak_fit
  */
template<typename realT>
struct array2FitExpModGaussian
{
   realT * data {nullptr}; ///< Pointer to the array
   size_t n {0}; ///< dimension of the array
   
   realT G0 {0}; ///< the location parameter.
   realT mu {0}; ///< the shape parameter
   realT sigma {0}; ///< the scale parameter
   realT lambda {0}; ///< the denominator or 1/peak-scale
};
 
///\ref levmarInterface fitter structure for the Exponentially Modified Gaussian with fixed constant level.
/**
  *
  * \ingroup gammaDist_peak_fit
  *
  */
template<typename _realT>
struct fitExpModGaussian_3param_fitter
{
   typedef _realT realT;
 
   static const int nparams = 3;
 
   static void func(realT *p, realT *hx, int m, int n, void *adata)
   {
      array2FitExpModGaussian<realT> * arr = (array2FitExpModGaussian<realT> *) adata;
 
      for(int i=0;i<arr->n; i++)
      {
         hx[i] = func::expModGaussian<realT>(i, p[0], p[1], p[2]) - arr->data[i];
      }
   }
};
 
///\ref levmarInterface fitter structure for the Exponentially Modified Gaussian with arbitrary constant level.
/**
  *
  * \ingroup gammaDist_peak_fit
  *
  */
template<typename _realT>
struct fitExpModGaussian_4param_fitter
{
   typedef _realT realT;
 
   static const int nparams = 4;
 
   static void func(realT *p, realT *hx, int m, int n, void *adata)
   {
      array2FitExpModGaussian<realT> * arr = (array2FitExpModGaussian<realT> *) adata;
 
      for(int i=0;i<arr->n; i++)
      {
         hx[i] = p[3] + func::expModGaussian<realT>(i, p[0], p[1], p[2]) - arr->data[i];
      }
   }
};
 
 
 
} //namespace fit
} //namespace math
} //namespace mx
 
#endif //fitExpModGaussian_hpp