#ifndef __JMATH__JMATHLIB2D__
#define __JMATH__JMATHLIB2D__

#include "JMath/JMathlib.hh"

/**
 * \file
 * Functional algebra in 2D.
 *
 * \author mdejong
 */

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

namespace JMATH {

  /**
   * Make 2D function of x from 1D function.
   */
  template<class JF1_t>
  struct JMake2X :
    public JMathlib< JMake2X<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JMake2X<JF1_t> >::operator();


    /**
     * Default constructor.
     */
    JMake2X()
    {}


    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JMake2X(const JF1_t& f1) :
      JF1_t(f1)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JMake2X(const Args& ...args) :
      JF1_t(args...)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y) const
    {
      return static_cast<const JF1_t&>(*this).getValue(x);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               gradient
     */
    const JMake2X& getGradient(const double x, const double y) const
    {
      static JMake2X gradient;

      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(x);

      return gradient;
    }
  };


  /**
   * Make 2D function of y from 1D function.
   */
  template<class JF1_t>
  struct JMake2Y :
    public JMathlib< JMake2Y<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JMake2Y<JF1_t> >::operator();


    /**
     * Default constructor.
     */
    JMake2Y()
    {}


    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JMake2Y(const JF1_t& f1) :
      JF1_t(f1)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JMake2Y(const Args& ...args) :
      JF1_t(args...)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y) const
    {
      return static_cast<const JF1_t&>(*this).getValue(y);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               gradient
     */
    const JMake2Y& getGradient(const double x, const double y) const
    {
      static JMake2Y gradient;

      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(y);

      return gradient;
    }
  };


  template<int ID_t, size_t N> using JPolynome2X = JMake2X< JPolynome<ID_t, N> >;                    //!< 2D polynomial function of x.
  template<int ID_t, size_t N> using JPolynome2Y = JMake2Y< JPolynome<ID_t, N> >;                    //!< 2D polynomial function of y.
  template<int ID_t, bool normalise = false> using JGauss2X = JMake2X< JGauss<ID_t, normalise> >;    //!< 2D Gauss function of x.
  template<int ID_t, bool normalise = false> using JGauss2Y = JMake2Y< JGauss<ID_t, normalise> >;    //!< 2D Gauss function of y.
  template<int ID_t, class JF1_t = _vF> using JPow2X = JMake2X< JPow<ID_t, JF1_t> >;                 //!< 2D power of function of x.
  template<int ID_t, class JF1_t = _vF> using JPow2Y = JMake2Y< JPow<ID_t, JF1_t> >;                 //!< 2D power of function of y.
  template<int N>  using JXn2X = JMake2X< JXn<N> >;                                                  //!< 2D fixed power of x.
  template<int N>  using JXn2Y = JMake2Y< JXn<N> >;                                                  //!< 2D fixed power of y.
  template<class JF1_t = _vF> using JSqrt2X = JMake2X< JSqrt<JF1_t> >;                               //!< 2D square root of function of x.
  template<class JF1_t = _vF> using JSqrt2Y = JMake2Y< JSqrt<JF1_t> >;                               //!< 2D square root of function of y.
  template<class JF1_t> using JSin2X = JMake2X< JSin<JF1_t> >;                                       //!< 2D sine of function of x.
  template<class JF1_t> using JSin2Y = JMake2Y< JSin<JF1_t> >;                                       //!< 2D sine of function of y.
  template<class JF1_t> using JCos2X = JMake2X< JCos<JF1_t> >;                                       //!< 2D cosine of function of x.
  template<class JF1_t> using JCos2Y = JMake2Y< JCos<JF1_t> >;                                       //!< 2D cosine of function of y.
  template<class JF1_t> using JExp2X = JMake2X< JExp<JF1_t> >;                                       //!< 2D exponent of function of x.
  template<class JF1_t> using JExp2Y = JMake2Y< JExp<JF1_t> >;                                       //!< 2D exponent of function of y.
  template<class JF1_t> using JLog2X = JMake2X< JLog<JF1_t> >;                                       //!< 2D logarithm of function of x.
  template<class JF1_t> using JLog2Y = JMake2Y< JLog<JF1_t> >;                                       //!< 2D logarithm of function of y.

  
  /**
   * 2D correlated Gauss function.
   */
  template<int ID_t, bool normalise = false>
  struct JGauss2D :
    public JMathlib< JGauss2D<ID_t, normalise> >,
    public JGauss<ID_t, normalise>
  {
    using JMathlib< JGauss2D<ID_t, normalise> >::operator();
    

    /**
     * Default constructor.
     */
    JGauss2D()
    {}


    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     */
    JGauss2D(const double mean,
             const double sigma) :
      JGauss<ID_t, normalise>(mean, sigma)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y) const
    {
      typedef JGauss<ID_t, normalise>  JGauss_t;

      return (static_cast<const JGauss_t&>(*this).getValue(x) *
              static_cast<const JGauss_t&>(*this).getValue(y));
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \return               gradient
     */
    const JGauss2D& getGradient(const double x, const double y) const
    {
      static JGauss2D gradient;

      typedef JGauss<ID_t, normalise>  JGauss_t;

      static_cast<JGauss_t&>(gradient)  = JGauss_t(static_cast<const JGauss_t&>(*this).getGradient(x)).mul(static_cast<const JGauss_t&>(*this).getValue(y));
      static_cast<JGauss_t&>(gradient) += JGauss_t(static_cast<const JGauss_t&>(*this).getGradient(y)).mul(static_cast<const JGauss_t&>(*this).getValue(x));

      return gradient;
    }
  };
}

#endif