#include <stan/math/prim.hpp>
#include <stan/math/rev/core.hpp>
#include <boost/math/special_functions/hypergeometric_pFq.hpp>

namespace model_gamma_model_namespace {
  inline auto gamma_icdf(const double& x, const double& alpha, const double& beta,
                           std::ostream* pstream__) {
    return boost::math::gamma_p_inv(alpha, x) / beta;
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const stan::math::var& x,
                                    const stan::math::var& alpha,
                                    const stan::math::var& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double x_val = x.val();
    double a_val = alpha.val();
    double b_val = beta.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(a_val, x_val);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / b_val;
    // Pre-compute some quantities for the gradients
    double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
    double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - a_val));
    double gamma_a_val = stan::math::tgamma(a_val);

    // Calculate derivative wrt 'x' parameter
    double x_d = (1 / b_val) * gamma_a_val * exp_gamma_p_inv_val * pow_gamma_p_inv_val;

    // Calculate derivative wrt to 'alpha' parameter
    double a_d = (1 / b_val) * exp_gamma_p_inv_val
                  * pow_gamma_p_inv_val
                  * (stan::math::square(gamma_a_val)
                     * std::pow(gamma_p_inv_val,a_val)
                        // Boost doesn't have a normalised Hypergeometric pFq, implement our own
                     *  (boost::math::hypergeometric_pFq({a_val,a_val},{a_val+1,a_val+1},-gamma_p_inv_val)
                          / stan::math::square(stan::math::tgamma(a_val+1)))
                     - x_val
                     * gamma_a_val
                     * std::log(gamma_p_inv_val)
                     + boost::math::polygamma(0, a_val)
                     * (gamma_a_val - boost::math::tgamma(a_val, gamma_p_inv_val)));

    // Calculate derivative wrt to 'beta' parameter
    double b_d = -gamma_p_inv_val / stan::math::square(b_val);

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[x, alpha, beta, x_d, a_d, b_d](auto& vi){
      x.adj() += vi.adj() * x_d;
      alpha.adj() += vi.adj() * a_d;
      beta.adj() += vi.adj() * b_d;
    });
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const double& x,
                                    const stan::math::var& alpha,
                                    const stan::math::var& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double x_val = x;
    double a_val = alpha.val();
    double b_val = beta.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(a_val, x_val);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / b_val;
    // Pre-compute some quantities for the gradients
    double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
    double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - a_val));
    double gamma_a_val = stan::math::tgamma(a_val);

    // Calculate derivative wrt to 'alpha' parameter
    double a_d = (1 / b_val) * exp_gamma_p_inv_val
                  * pow_gamma_p_inv_val
                  * (stan::math::square(gamma_a_val)
                     * std::pow(gamma_p_inv_val,a_val)
                        // Boost doesn't have a normalised Hypergeometric pFq, implement our own
                     *  (boost::math::hypergeometric_pFq({a_val,a_val},{a_val+1,a_val+1},-gamma_p_inv_val)
                          / stan::math::square(stan::math::tgamma(a_val+1)))
                     - x_val
                     * gamma_a_val
                     * std::log(gamma_p_inv_val)
                     + boost::math::polygamma(0, a_val)
                     * (gamma_a_val - boost::math::tgamma(a_val, gamma_p_inv_val)));

    // Calculate derivative wrt to 'beta' parameter
    double b_d = -gamma_p_inv_val / stan::math::square(b_val);

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[alpha, beta, a_d, b_d](auto& vi){
      alpha.adj() += vi.adj() * a_d;
      beta.adj() += vi.adj() * b_d;
    });
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const double& x,
                                    const double& alpha,
                                    const stan::math::var& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double b_val = beta.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(alpha, x);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / b_val;

    // Calculate derivative wrt to 'beta' parameter
    double b_d = -gamma_p_inv_val / stan::math::square(b_val);

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[beta, b_d](auto& vi){
      beta.adj() += vi.adj() * b_d;
    });
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const stan::math::var& x,
                                    const double& alpha,
                                    const stan::math::var& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double x_val = x.val();
    double b_val = beta.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(alpha, x_val);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / b_val;
    // Pre-compute some quantities for the gradients
    double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
    double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - alpha));
    double gamma_a_val = stan::math::tgamma(alpha);

    // Calculate derivative wrt 'x' parameter
    double x_d = (1 / b_val) * gamma_a_val * exp_gamma_p_inv_val * pow_gamma_p_inv_val;

    // Calculate derivative wrt to 'beta' parameter
    double b_d = -gamma_p_inv_val / stan::math::square(b_val);

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[x, beta, x_d, b_d](auto& vi){
      x.adj() += vi.adj() * x_d;
      beta.adj() += vi.adj() * b_d;
    });
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const stan::math::var& x,
                                    const double& alpha,
                                    const double& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double x_val = x.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(alpha, x_val);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / beta;
    // Pre-compute some quantities for the gradients
    double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
    double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - alpha));
    double gamma_a_val = stan::math::tgamma(alpha);

    // Calculate derivative wrt 'x' parameter
    double x_d = (1 / beta) * gamma_a_val * exp_gamma_p_inv_val * pow_gamma_p_inv_val;

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[x, x_d](auto& vi){
      x.adj() += vi.adj() * x_d;
    });
  }

  // Function to handle Stan parameters
  inline stan::math::var gamma_icdf(const stan::math::var& x,
                                    const stan::math::var& alpha,
                                    const double& beta,
                                    std::ostream* pstream__) {
    // Extract values (not derivatives) of inputs
    double x_val = x.val();
    double a_val = alpha.val();

    // Calculate inverse of the incomplete gamma function
    double gamma_p_inv_val = boost::math::gamma_p_inv(a_val, x_val);
    // Use to derive quantile of Gamma distribution
    double gamma_icdf = gamma_p_inv_val / beta;
    // Pre-compute some quantities for the gradients
    double exp_gamma_p_inv_val = std::exp(gamma_p_inv_val);
    double pow_gamma_p_inv_val = std::pow(gamma_p_inv_val,(1 - a_val));
    double gamma_a_val = stan::math::tgamma(a_val);

    // Calculate derivative wrt 'x' parameter
    double x_d = (1 / beta) * gamma_a_val * exp_gamma_p_inv_val * pow_gamma_p_inv_val;

    // Calculate derivative wrt to 'alpha' parameter
    double a_d = (1 / beta) * exp_gamma_p_inv_val
                  * pow_gamma_p_inv_val
                  * (stan::math::square(gamma_a_val)
                     * std::pow(gamma_p_inv_val,a_val)
                        // Boost doesn't have a normalised Hypergeometric pFq, implement our own
                     *  (boost::math::hypergeometric_pFq({a_val,a_val},{a_val+1,a_val+1},-gamma_p_inv_val)
                          / stan::math::square(stan::math::tgamma(a_val+1)))
                     - x_val
                     * gamma_a_val
                     * std::log(gamma_p_inv_val)
                     + boost::math::polygamma(0, a_val)
                     * (gamma_a_val - boost::math::tgamma(a_val, gamma_p_inv_val)));

    // Return new parameter with values of the function and gradients
    return stan::math::make_callback_var(gamma_icdf,[x, alpha, x_d, a_d](auto& vi){
      x.adj() += vi.adj() * x_d;
      alpha.adj() += vi.adj() * a_d;
    });
  }
}