#ifndef __JMATH__JGAUSS__
#define __JMATH__JGAUSS__

#include <istream>
#include <ostream>
#include <iomanip>
#include <cmath>
#include <limits>

#include "JLang/JEquals.hh"
#include "JLang/JManip.hh"
#include "JMath/JMath.hh"
#include "JMath/JConstants.hh"

/**
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  using JLANG::JEquals;

  /**
   * Gauss model.
   */
  struct JGauss_t :
    public JMath<JGauss_t>,
    public JEquals<JGauss_t>
  {
    /**
     * Default constructor.
     */
    JGauss_t() :
      mean      (0.0),
      sigma     (0.0),
      signal    (0.0),
      background(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     * \param  signal        signal
     * \param  background    background
     */
    JGauss_t(const double mean,
	     const double sigma,
	     const double signal,
	     const double background) :
      mean      (mean),
      sigma     (sigma),
      signal    (signal),
      background(background)
    {}


    /**
     * Equality.
     *
     * \param  gauss         gauss
     * \param  eps           numerical precision
     * \return               true if gauss's identical; else false
     */
    bool equals(const JGauss_t& gauss,
                const double    eps = std::numeric_limits<double>::min()) const
    {
      return (fabs(mean       - gauss.mean)       <=  eps  &&
	      fabs(sigma      - gauss.sigma)      <=  eps  &&
	      fabs(signal     - gauss.signal)     <=  eps  &&
	      fabs(background - gauss.background) <=  eps);
    }


    /**
     * Add gauss.
     *
     * \param  gauss         gauss
     * \return               this gauss
     */
    JGauss_t& add(const JGauss_t& gauss)
    {
      mean       += gauss.mean;
      sigma      += gauss.sigma;
      signal     += gauss.signal;
      background += gauss.background;

      return *this;
    }


    /**
     * Subtract gauss.
     *
     * \param  gauss         gauss
     * \return               this gauss
     */
    JGauss_t& sub(const JGauss_t& gauss)
    {
      mean       -= gauss.mean;
      sigma      -= gauss.sigma;
      signal     -= gauss.signal;
      background -= gauss.background;

      return *this;
    }


    /**
     * Scale gauss.
     *
     * \param  factor        multiplication factor
     * \return               this gauss
     */
    JGauss_t& mul(const double factor)
    {
      mean       *= factor;
      sigma      *= factor;
      signal     *= factor;
      background *= factor;

      return *this;
    }


    /**
     * Write Gauss to input stream.
     *
     * \param  in            input stream
     * \param  gauss         gauss
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JGauss_t& gauss) 
    {
      return in >> gauss.mean >> gauss.sigma >> gauss.signal >> gauss.background;
    }


    /**
     * Write Gauss to output stream.
     *
     * \param  out           output stream
     * \param  gauss         gauss
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JGauss_t& gauss) 
    {
      using namespace std;

      return out << FIXED(7,3) << gauss.mean   << ' '
		 << FIXED(7,3) << gauss.sigma  << ' '
		 << FIXED(9,3) << gauss.signal << ' '
		 << FIXED(9,3) << gauss.background;
    }

    double mean;
    double sigma;
    double signal;
    double background;
  };


  /**
   * Gauss function object.
   *
   * Evaluates function, derivative and gradient values.
   */
  struct JGauss :
    public JGauss_t
  {
    /**
     * Type definition of fit parameter.
     */
    typedef double JGauss::*parameter_type;


    /**
     * Default constructor.
     */
    JGauss() :
      JGauss_t()
    {}


    /**
     * Copy constructor.
     *
     * \param  gauss         gauss
     */
    JGauss(const JGauss_t& gauss) :
      JGauss_t(gauss)
    {}


    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     * \param  signal        signal
     * \param  background    background
     */
    JGauss(const double mean,
	   const double sigma,
	   const double signal     = 1.0,
	   const double background = 0.0) :
      JGauss_t(mean, sigma, signal, background)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u  = (x - mean) / sigma;

      return signal * get(u) + background;
    }


    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u  = (x - mean) / sigma;

      return signal * get(u) * -u;
    }


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double operator()(const double x) const
    {
      return getValue(x);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    const JGauss_t& getGradient(const double x) const
    {
      const double w  = 1.0 / sigma;
      const double u  = (x - mean) * w;
      const double f0 =  get(u);
      const double fs =  signal * f0;

      gradient.mean       =  fs * (u)                   * w;             // d(f)/d(mean)
      gradient.sigma      =  fs * (u + 1.0) * (u - 1.0) * w;             // d(f)/d(sigma)
      gradient.signal     =  f0;                                         // d(f)/d(signal)
      gradient.background =  1.0;                                        // d(f)/d(background)

      return gradient;
    }

  private:
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      return exp(-0.5*u*u) / (sqrt(2.0*PI) * sigma);
    }

    mutable JGauss_t gradient;
  };
}

#endif