special_functions.hpp

Go to the documentation of this file.

#ifndef __STAN__AGRAD__AGRAD_SPECIAL_FUNCTIONS_HPP__
#define __STAN__AGRAD__AGRAD_SPECIAL_FUNCTIONS_HPP__

#include <stan/agrad/error_handling.hpp>
#include <stan/math/special_functions.hpp>

#include <boost/math/special_functions/acosh.hpp>
#include <boost/math/special_functions/asinh.hpp>
#include <boost/math/special_functions/atanh.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/hypot.hpp>

namespace stan {

  namespace agrad {
    
    namespace {
      
      class lgamma_vari : public op_v_vari {
      public:
        lgamma_vari(double value, vari* avi) :
          op_v_vari(value, avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * boost::math::digamma(avi_->val_);
        }
      };

      class tgamma_vari : public op_v_vari {
      public:
        tgamma_vari(vari* avi) :
          op_v_vari(boost::math::tgamma(avi->val_), avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * val_ * boost::math::digamma(avi_->val_);
        }
      };

      class log1p_vari : public op_v_vari {
      public:
        log1p_vari(vari* avi) :
          op_v_vari(stan::math::log1p(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (1 + avi_->val_);
        }
      };

      class log1m_vari : public op_v_vari {
      public:
        log1m_vari(vari* avi) :
          op_v_vari(stan::math::log1p(-avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (avi_->val_ - 1);
        }
      };

      class binary_log_loss_1_vari : public op_v_vari {
      public:
        binary_log_loss_1_vari(vari* avi) :
          op_v_vari(-std::log(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ / avi_->val_;
        }
      };

      class binary_log_loss_0_vari : public op_v_vari {
      public:
        binary_log_loss_0_vari(vari* avi) :
          op_v_vari(-stan::math::log1p(-avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (1.0 - avi_->val_);
        }
      };

      class fdim_vv_vari : public op_vv_vari {
      public:
        fdim_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ - bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_;
          bvi_->adj_ -= adj_;
        }
      };

      class fdim_vd_vari : public op_v_vari {
      public:
        fdim_vd_vari(vari* avi, double b) :
          op_v_vari(avi->val_ - b, avi) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };

      class fdim_dv_vari : public op_v_vari {
      public:
        fdim_dv_vari(double a, vari* bvi) :
          op_v_vari(a - bvi->val_, bvi) {
        }
        void chain() {
          // avi_ is bvi argument to constructor
          avi_->adj_ -= adj_;
        }
      };

      class fma_vvv_vari : public op_vvv_vari {
      public:
        fma_vvv_vari(vari* avi, vari* bvi, vari* cvi) :
          op_vvv_vari(avi->val_ * bvi->val_ + cvi->val_,
                      avi,bvi,cvi) {
        }
        void chain() {
          avi_->adj_ += adj_ * bvi_->val_;
          bvi_->adj_ += adj_ * avi_->val_;
          cvi_->adj_ += adj_;
        }
      };

      class fma_vvd_vari : public op_vv_vari {
      public:
        fma_vvd_vari(vari* avi, vari* bvi, double c) :
          op_vv_vari(avi->val_ * bvi->val_ + c,
                     avi,bvi) {
        }
        void chain() {
          avi_->adj_ += adj_ * bvi_->val_;
          bvi_->adj_ += adj_ * avi_->val_;
        }
      };

      class fma_vdv_vari : public op_vdv_vari {
      public:
        fma_vdv_vari(vari* avi, double b, vari* cvi) :
          op_vdv_vari(avi->val_ * b + cvi->val_,
                      avi,b,cvi) {
        }
        void chain() {
          avi_->adj_ += adj_ * bd_;
          cvi_->adj_ += adj_;
        }
      };

      class fma_vdd_vari : public op_vd_vari {
      public:
        fma_vdd_vari(vari* avi, double b, double c) : 
          op_vd_vari(avi->val_ * b + c,
                     avi,b) {
        }
        void chain() {
          avi_->adj_ += adj_ * bd_;
        }
      };

      class fma_ddv_vari : public op_v_vari {
      public:
        fma_ddv_vari(double a, double b, vari* cvi) :
          op_v_vari(a * b + cvi->val_, 
                    cvi) {
        }
        void chain() {
          // avi_ is cvi from constructor
          avi_->adj_ += adj_;
        }
      };

      class inv_logit_vari : public op_v_vari {
      public:
        inv_logit_vari(vari* avi) :
          op_v_vari(math::inv_logit(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ +=  adj_ * val_ * (1.0 - val_);
        }
      };

      class acosh_vari : public op_v_vari {
      public:
        acosh_vari(vari* avi) :
          op_v_vari(boost::math::acosh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / std::sqrt(avi_->val_ * avi_->val_ - 1.0);
        }
      };

      class asinh_vari : public op_v_vari {
      public:
        asinh_vari(vari* avi) :
          op_v_vari(boost::math::asinh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / std::sqrt(avi_->val_ * avi_->val_ + 1.0);
        }
      };

      class atanh_vari : public op_v_vari {
      public:
        atanh_vari(vari* avi) :
          op_v_vari(boost::math::atanh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (1.0 - avi_->val_ * avi_->val_);
        }
      };

      const double TWO_OVER_SQRT_PI = 2.0 / std::sqrt(boost::math::constants::pi<double>());

      class erf_vari : public op_v_vari {
      public:
        erf_vari(vari* avi) :
          op_v_vari(boost::math::erf(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * TWO_OVER_SQRT_PI * std::exp(- avi_->val_ * avi_->val_);
        }
      };

      const double NEG_TWO_OVER_SQRT_PI = - TWO_OVER_SQRT_PI;

      class erfc_vari : public op_v_vari {
      public:
        erfc_vari(vari* avi) :
          op_v_vari(boost::math::erfc(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * NEG_TWO_OVER_SQRT_PI * std::exp(- avi_->val_ * avi_->val_);
        }
      };

      const double LOG_2 = std::log(2.0);

      class exp2_vari : public op_v_vari {
      public:
        exp2_vari(vari* avi) :
          op_v_vari(std::pow(2.0,avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * val_ * LOG_2;
        }
      };

      class expm1_vari : public op_v_vari {
      public:
        expm1_vari(vari* avi) :
          op_v_vari(std::exp(avi->val_) - 1.0,avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * val_;
        }
      };

      class hypot_vv_vari : public op_vv_vari {
      public:
        hypot_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(boost::math::hypot(avi->val_,bvi->val_),
                     avi,bvi) {
        }
        void chain() {
          avi_->adj_ += adj_ * avi_->val_ / val_;
          bvi_->adj_ += adj_ * bvi_->val_ / val_;
        }
      };

      class hypot_vd_vari : public op_v_vari {
      public:
        hypot_vd_vari(vari* avi, double b) :
          op_v_vari(boost::math::hypot(avi->val_,b),
                    avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * avi_->val_ / val_;
        }
      };

      const double LOG2 = std::log(2.0);

      class log2_vari : public op_v_vari {
      public:
        log2_vari(vari* avi) :
          op_v_vari(stan::math::log2(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (LOG2 * avi_->val_); 
        }
      };

      class cbrt_vari : public op_v_vari {
      public:
        cbrt_vari(vari* avi) :
          op_v_vari(boost::math::cbrt(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (3.0 * val_ * val_);
        }
      };

      // derivative 0 almost everywhere
      class round_vari : public vari {
      public:
        round_vari(vari* avi) :
          vari(boost::math::round(avi->val_)) {
        }
      };

      // derivative 0 almost everywhere
      class trunc_vari : public vari {
      public:
        trunc_vari(vari* avi) :
          vari(boost::math::trunc(avi->val_)) { 
        }
      };

      class inv_cloglog_vari : public op_v_vari {
      public:
        inv_cloglog_vari(vari* avi) :
          op_v_vari(stan::math::inv_cloglog(avi->val_), avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ * std::exp(avi_->val_ - std::exp(avi_->val_));
        }
      };

      class Phi_vari : public op_v_vari {
      public:
        Phi_vari(vari* avi) :
          op_v_vari(stan::math::Phi(avi->val_), avi) {
        }
        void chain() {
          static const double NEG_HALF = -0.5;
          static const double INV_SQRT_TWO_PI 
            = 1.0 / std::sqrt(2.0 * boost::math::constants::pi<double>());
          avi_->adj_ += adj_ * INV_SQRT_TWO_PI * std::exp(NEG_HALF * avi_->val_ * avi_->val_);
        }
      };

      inline double calculate_chain(const double& x, const double& val) {
        return std::exp(x - val); // works out to inv_logit(x)
      }

      double log_sum_exp_as_double(const std::vector<var>& x) {
        using std::numeric_limits;
        using std::exp;
        using std::log;
        double max = -numeric_limits<double>::infinity();
        for (size_t i = 0; i < x.size(); ++i) 
          if (x[i] > max) 
            max = x[i].val();
        double sum = 0.0;
        for (size_t i = 0; i < x.size(); ++i) 
          if (x[i] != -numeric_limits<double>::infinity()) 
            sum += exp(x[i].val() - max);
        return max + log(sum);
      }

      class log1p_exp_v_vari : public op_v_vari {
      public:
        log1p_exp_v_vari(vari* avi) :
          op_v_vari(stan::math::log1p_exp(avi->val_),
                    avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
        }
      };      
      class log_sum_exp_vv_vari : public op_vv_vari {
      public:
        log_sum_exp_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(stan::math::log_sum_exp(avi->val_, bvi->val_),
                     avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
          bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_);
        }
      };
      class log_sum_exp_vd_vari : public op_vd_vari {
      public:
        log_sum_exp_vd_vari(vari* avi, double b) :
          op_vd_vari(stan::math::log_sum_exp(avi->val_, b),
                     avi, b) {
        }
        void chain() {
          avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
        }
      };
      class log_sum_exp_dv_vari : public op_dv_vari {
      public:
        log_sum_exp_dv_vari(double a, vari* bvi) :
          op_dv_vari(stan::math::log_sum_exp(a, bvi->val_),
                     a, bvi) {
        }
        void chain() {
          bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_);
        }
      };
      class log_sum_exp_vector_vari : public op_vector_vari {
      public:
        log_sum_exp_vector_vari(const std::vector<var>& x) :
          op_vector_vari(log_sum_exp_as_double(x), x) {
        }
        void chain() {
          for (size_t i = 0; i < size_; ++i) {
            vis_[i]->adj_ += adj_ * calculate_chain(vis_[i]->val_, val_);
          }
        }
      };

      class square_vari : public op_v_vari {
      public:
        square_vari(vari* avi) :
          op_v_vari(avi->val_ * avi->val_,avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * 2.0 * avi_->val_;
        }
      };

      class multiply_log_vv_vari : public op_vv_vari {
      public:
        multiply_log_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(stan::math::multiply_log(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          using std::log;
          avi_->adj_ += adj_ * log(bvi_->val_);
          if (bvi_->val_ == 0.0 && avi_->val_ == 0)
            bvi_->adj_ += adj_ * std::numeric_limits<double>::infinity();
          else
            bvi_->adj_ += adj_ * avi_->val_ / bvi_->val_;
        }
      };
      class multiply_log_vd_vari : public op_vd_vari {
      public:
        multiply_log_vd_vari(vari* avi, double b) :
          op_vd_vari(stan::math::multiply_log(avi->val_,b),avi,b) {
        }
        void chain() {
          using std::log;
          avi_->adj_ += adj_ * log(bd_);
        }
      };
      class multiply_log_dv_vari : public op_dv_vari {
      public:
        multiply_log_dv_vari(double a, vari* bvi) :
          op_dv_vari(stan::math::multiply_log(a,bvi->val_),a,bvi) {
        }
        void chain() {
          if (bvi_->val_ == 0.0 && ad_ == 0.0)
            bvi_->adj_ += adj_ * std::numeric_limits<double>::infinity();
          else
            bvi_->adj_ += adj_ * ad_ / bvi_->val_;
        }
      };
      namespace {
        double ibeta_hypergeometric_helper(double a, double b, double z, double precision=1e-8, double max_steps=1000) {
          double val=0;
          double diff=1;
          double k=0;
          double a_2 = a*a;
          double bprod = 1;
          while (std::abs(diff) > precision && (++k < max_steps) && !std::isnan(diff)) {
            val += diff;
            bprod *= b+k-1.0;
            diff = a_2*std::pow(a+k,-2)*bprod*std::pow(z,k)/boost::math::tgamma(k+1);
          }
          return val;
        }
      }
      class ibeta_vvv_vari : public op_vvv_vari {
      public:
        ibeta_vvv_vari(vari* avi, vari* bvi, vari* xvi) :
          op_vvv_vari(stan::math::ibeta(avi->val_,bvi->val_,xvi->val_),avi,bvi,xvi) {
        }
        void chain() {
          double a = avi_->val_;
          double b = bvi_->val_;
          double c = cvi_->val_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          avi_->adj_ += adj_ *
            (log(c) - digamma(a) + digamma(a+b)) * val_ - 
            tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
          bvi_->adj_ += adj_ * 
            (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
             + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
          cvi_->adj_ += adj_ * 
            boost::math::ibeta_derivative(a,b,c);
        }
      };
      class ibeta_vvd_vari : public op_vvd_vari {
      public:
        ibeta_vvd_vari(vari* avi, vari* bvi, double x) :
          op_vvd_vari(stan::math::ibeta(avi->val_,bvi->val_,x),avi,bvi,x) {
        }
        void chain() {
          double a = avi_->val_;
          double b = bvi_->val_;
          double c = cd_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          avi_->adj_ += adj_ *
            (log(c) - digamma(a) + digamma(a+b)) * val_ - 
            tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
          bvi_->adj_ += adj_ * 
            (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
             + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
        }
      };
      class ibeta_vdv_vari : public op_vdv_vari {
      public:
        ibeta_vdv_vari(vari* avi, double b, vari* xvi) :
          op_vdv_vari(stan::math::ibeta(avi->val_,b,xvi->val_),avi,b,xvi) {
        }
        void chain() {
          double a = avi_->val_;
          double b = bd_;
          double c = cvi_->val_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          avi_->adj_ += adj_ *
            (log(c) - digamma(a) + digamma(a+b)) * val_ - 
            tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
          cvi_->adj_ += adj_ * 
            boost::math::ibeta_derivative(a,b,c);
        }
      };
      class ibeta_vdd_vari : public op_vdd_vari {
      public:
        ibeta_vdd_vari(vari* avi, double b, double x) :
          op_vdd_vari(stan::math::ibeta(avi->val_,b,x),avi,b,x) {
        }
        void chain() {
          double a = avi_->val_;
          double b = bd_;
          double c = cd_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          avi_->adj_ += adj_ *
            (log(c) - digamma(a) + digamma(a+b)) * val_ - 
            tgamma(a)*tgamma(a+b)/tgamma(b) * pow(c,a) / tgamma(1+a) / tgamma(1+a) * ibeta_hypergeometric_helper(a, 1-b, c);
        }
      };
      class ibeta_dvv_vari : public op_dvv_vari {
      public:
        ibeta_dvv_vari(double a, vari* bvi, vari* xvi) :
          op_dvv_vari(stan::math::ibeta(a,bvi->val_,xvi->val_),a,bvi,xvi) {
        }
        void chain() {
          double a = ad_;
          double b = bvi_->val_;
          double c = cvi_->val_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          bvi_->adj_ += adj_ * 
            (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
             + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
          cvi_->adj_ += adj_ * 
            boost::math::ibeta_derivative(a,b,c);
        }
      };
      class ibeta_dvd_vari : public op_dvd_vari {
      public:
        ibeta_dvd_vari(double a, vari* bvi, double x) :
          op_dvd_vari(stan::math::ibeta(a,bvi->val_,x),a,bvi,x) {
        }
        void chain() {
          double a = ad_;
          double b = bvi_->val_;
          double c = cd_;

          using std::sin;
          using std::pow;
          using std::log;
          using boost::math::constants::pi;
          using boost::math::tgamma;
          using boost::math::digamma;
          using boost::math::ibeta;
          using stan::agrad::ibeta_hypergeometric_helper;
          bvi_->adj_ += adj_ * 
            (tgamma(b)*tgamma(a+b)/tgamma(a)*pow(1-c,b) * ibeta_hypergeometric_helper(b,1-a,1-c)/tgamma(b+1)/tgamma(b+1)
             + ibeta(b, a, 1-c) * (digamma(b) - digamma(a+b) - log(1-c)));
        }
      };
      class ibeta_ddv_vari : public op_ddv_vari {
      public:
        ibeta_ddv_vari(double a, double b, vari* xvi) :
          op_ddv_vari(stan::math::ibeta(a,b,xvi->val_),a,b,xvi) {
        }
        void chain() {
          double a = ad_;
          double b = bd_;
          double c = cvi_->val_;

          cvi_->adj_ += adj_ * 
            boost::math::ibeta_derivative(a,b,c);
        }
      };
    }


    using stan::math::check_not_nan;
    using stan::math::check_greater_or_equal;

    inline var acosh(const stan::agrad::var& a) {
      return var(new acosh_vari(a.vi_));
    }

    inline var asinh(const stan::agrad::var& a) {
      return var(new asinh_vari(a.vi_));
    }

    inline var atanh(const stan::agrad::var& a) {
      return var(new atanh_vari(a.vi_));
    }

    inline var erf(const stan::agrad::var& a) {
      return var(new erf_vari(a.vi_));
    }

    inline var erfc(const stan::agrad::var& a) {
      return var(new erfc_vari(a.vi_));
    }

    inline var exp2(const stan::agrad::var& a) {
      return var(new exp2_vari(a.vi_));
    }

    inline var expm1(const stan::agrad::var& a) {
      return var(new expm1_vari(a.vi_));
    }

    inline var lgamma(const stan::agrad::var& a) {
      double lgamma_a = boost::math::lgamma(a.val());
      return var(new lgamma_vari(lgamma_a, a.vi_));
    }

    inline var log1p(const stan::agrad::var& a) {
      return var(new log1p_vari(a.vi_));
    }

    inline var log1m(const stan::agrad::var& a) {
      return var(new log1m_vari(a.vi_));
    }

    inline var fma(const stan::agrad::var& a,
                   const stan::agrad::var& b,
                   const stan::agrad::var& c) {
      return var(new fma_vvv_vari(a.vi_,b.vi_,c.vi_));
    }

    inline var fma(const stan::agrad::var& a,
                   const stan::agrad::var& b,
                   const double& c) {
      return var(new fma_vvd_vari(a.vi_,b.vi_,c));
    }

    inline var fma(const stan::agrad::var& a,
                   const double& b,
                   const stan::agrad::var& c) {
      return var(new fma_vdv_vari(a.vi_,b,c.vi_));
    }

    inline var fma(const stan::agrad::var& a,
                   const double& b, 
                   const double& c) {
      return var(new fma_vdd_vari(a.vi_,b,c));
    }

    inline var fma(const double& a,
                   const stan::agrad::var& b,
                   const double& c) {
      return var(new fma_vdd_vari(b.vi_,a,c));
    }

    inline var fma(const double& a,
                   const double& b,
                   const stan::agrad::var& c) {
      return var(new fma_ddv_vari(a,b,c.vi_));
    }

    inline var fma(const double& a,
                   const stan::agrad::var& b,
                   const stan::agrad::var& c) {
      return var(new fma_vdv_vari(b.vi_,a,c.vi_)); // a-b symmetry
    }


    inline var fmax(const stan::agrad::var& a,
                    const stan::agrad::var& b) {
      return a.vi_->val_ > b.vi_->val_ ? a : b;
    }

    inline var fmax(const stan::agrad::var& a,
                    const double& b) {
      return a.vi_->val_ >= b ? a : var(b);
    }

    inline var fmax(const double& a,
                    const stan::agrad::var& b) {
      return a > b.vi_->val_ ? var(a) : b;
    }

    inline var fmin(const stan::agrad::var& a,
                    const stan::agrad::var& b) {
      return a.vi_->val_ < b.vi_->val_ ? a : b;
    }

    inline var fmin(const stan::agrad::var& a,
                    const double& b) {
      return a.vi_->val_ <= b ? a : var(b);
    }

    inline var fmin(const double& a,
                    const stan::agrad::var& b) {
      return a < b.vi_->val_ ? var(a) : b;
    }

    inline var hypot(const stan::agrad::var& a,
                     const stan::agrad::var& b) {
      return var(new hypot_vv_vari(a.vi_,b.vi_));
    }

    inline var hypot(const stan::agrad::var& a,
                     const double& b) {
      return var(new hypot_vd_vari(a.vi_,b));
    }

    inline var hypot(const double& a,
                     const stan::agrad::var& b) {
      return var(new hypot_vd_vari(b.vi_,a));
    }

    inline var log2(const stan::agrad::var& a) {
      return var(new log2_vari(a.vi_));
    }

    inline var cbrt(const stan::agrad::var& a) {
      return var(new cbrt_vari(a.vi_));
    }

    inline var round(const stan::agrad::var& a) {
      return var(new round_vari(a.vi_));
    }

    inline var trunc(const stan::agrad::var& a) {
      return var(new trunc_vari(a.vi_));
    }


    inline var fdim(const stan::agrad::var& a,
                    const stan::agrad::var& b) {
      if (a.vi_->val_ > b.vi_->val_)
        return var(new fdim_vv_vari(a.vi_,b.vi_));
      else
        return var(new vari(0.0));
    }

    inline var fdim(const double& a,
                    const stan::agrad::var& b) {
      return a > b.vi_->val_
        ? var(new fdim_dv_vari(a,b.vi_))
        : var(new vari(0.0));
    }

    inline var fdim(const stan::agrad::var& a,
                    const double& b) {
      return a.vi_->val_ > b
        ? var(new fdim_vd_vari(a.vi_,b))
        : var(new vari(0.0));
    }

    inline var tgamma(const stan::agrad::var& a) {
      return var(new tgamma_vari(a.vi_));
    }

    inline var step(const stan::agrad::var& a) {
      return var(new vari(a.vi_->val_ < 0.0 ? 0.0 : 1.0));
    }

    inline var inv_cloglog(const stan::agrad::var& a) {
      return var(new inv_cloglog_vari(a.vi_));
    }

    inline var Phi(const stan::agrad::var& a) {
      return var(new Phi_vari(a.vi_));
    }

    inline var inv_logit(const stan::agrad::var& a) {
      return var(new inv_logit_vari(a.vi_));
    }

    inline var log_loss(const int& y, 
                        const stan::agrad::var& y_hat) {
      return y == 0  
        ? var(new binary_log_loss_0_vari(y_hat.vi_))
        : var(new binary_log_loss_1_vari(y_hat.vi_));
    }
    
    inline var log1p_exp(const stan::agrad::var& a) {
      return var(new log1p_exp_v_vari(a.vi_));
    }

    inline var log_sum_exp(const stan::agrad::var& a,
                           const stan::agrad::var& b) {
      return var(new log_sum_exp_vv_vari(a.vi_, b.vi_));
    }
    inline var log_sum_exp(const stan::agrad::var& a,
                           const double& b) {
      return var(new log_sum_exp_vd_vari(a.vi_, b));
    }
    inline var log_sum_exp(const double& a,
                           const stan::agrad::var& b) {
      return var(new log_sum_exp_dv_vari(a, b.vi_));
    }
    inline var log_sum_exp(const std::vector<var>& x) {
      return var(new log_sum_exp_vector_vari(x));
    }

    inline var square(const var& x) {
      return var(new square_vari(x.vi_));
    }

    // OTHER FUNCTIONS: stan/math/special_functions.hpp implementations
    inline var multiply_log(const var& a, const var& b) {
      return var(new multiply_log_vv_vari(a.vi_,b.vi_));
    }
    inline var multiply_log(const var& a, const double b) {
      return var(new multiply_log_vd_vari(a.vi_,b));
    }
    inline var multiply_log(const double a, const var& b) {
      if (a == 1.0)
        return log(b);
      return var(new multiply_log_dv_vari(a,b.vi_));
    }


    inline var if_else(bool c, const var& y_true, const var&y_false) {
      return c ? y_true : y_false;
    }
    inline var if_else(bool c, double y_true, const var& y_false) {
      if (c) 
        return var(y_true);
      else 
        return y_false;
    }
    inline var if_else(bool c, const var& y_true, const double y_false) {
      if (c) 
        return y_true;
      else 
        return var(y_false);
    }

    inline var ibeta(const var& a,
                     const var& b,
                     const var& x) {
      return var(new ibeta_vvv_vari(a.vi_, b.vi_, x.vi_));
    }

    inline double value_of(const agrad::var& v) {
      return v.vi_->val_;
    }

    inline int as_bool(const agrad::var& v) {
      return 0.0 != v.vi_->val_;
    }

  } // namespace math

} // namespace stan


#endif