special_functions.hpp

Go to the documentation of this file.

#ifndef __STAN__MATH__SPECIAL_FUNCTIONS_HPP__
#define __STAN__MATH__SPECIAL_FUNCTIONS_HPP__

#include <stdexcept>


#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/throw_exception.hpp>

#include <stan/math/constants.hpp>
#include <stan/math/error_handling.hpp>

namespace stan {

  namespace math {

    // C99 

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    exp2(T y) {
      using std::pow;
      return pow(2.0,y);
    }

    template <typename T1, typename T2>
    inline typename boost::math::tools::promote_args<T1, T2>::type
    fdim(T1 a, T2 b) {
      return (a > b) ? (a - b) : 0.0;
    } 
    
    template <typename T1, typename T2, typename T3>
    inline typename boost::math::tools::promote_args<T1,T2,T3>::type
    fma(T1 a, T2 b, T3 c) {
      return (a * b) + c;
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    log2(T a) {
      using std::log;
      const static double LOG2 = std::log(2.0);
      return log(a) / LOG2;
    }

    // OTHER BASIC FUNCTIONS

    template <typename T>
    unsigned int int_step(T y) {
      return y > 0;
    }

    template <typename T>
    inline int step(T y) {
      return y < 0.0 ? 0 : 1;
    }


    // PROBABILITY-RELATED FUNCTIONS
    
    template <typename T1, typename T2>
    inline typename boost::math::tools::promote_args<T1,T2>::type
    lbeta(T1 a, T2 b) {
      return lgamma(a)
        + lgamma(b)
        - lgamma(a + b);
    }

    template <typename T_N, typename T_n>
    inline typename boost::math::tools::promote_args<T_N, T_n>::type
    binomial_coefficient_log(T_N N, T_n n) {
      using std::log;

      const double cutoff = 1000;
      if ((N < cutoff) || (N - n < cutoff)) {
        return lgamma(N + 1.0) - lgamma(n + 1.0) - lgamma(N - n + 1.0);
      } else {
        return n * log(N - n) + (N + 0.5) * log(N/(N-n))
          + 1/(12*N) - n - 1/(12*(N-n)) - lgamma(n + 1.0);
      }
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    inv_logit(T a) {
      using std::exp;
      return 1.0 / (1.0 + exp(-a));
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    logit(T a) {
      using std::log;
      return log(a / (1.0 - a));
    }



    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    Phi(T x) {
      static const double INV_SQRT_TWO = 1.0 / std::sqrt(2.0);
      return 0.5 * (1.0 + boost::math::erf(INV_SQRT_TWO * x));
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    inv_cloglog(T x) {
      using std::exp;
      return exp(-exp(x));
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    binary_log_loss(int y, T y_hat) {
      using std::log;
      return -log(y ? y_hat : (1.0 - y_hat));
    }

    // hide helper for now; could use Eigen here
    namespace {
      template <typename Vector, typename Scalar>
      Scalar maximum(const Vector& x) {
        if(x.size() == 0)
          BOOST_THROW_EXCEPTION(std::invalid_argument ("x must have at least one element"));
        Scalar max_x(x[0]);
        for (typename Vector::size_type i = 1; i < x.size(); ++i)
          if (x[i] < max_x)
            max_x = x[i];
        return max_x;
      }
    }


    template <typename Vector, typename Scalar>
    void softmax(const Vector& x, Vector& simplex) {
      using std::exp;
      if(x.size() != simplex.size()) 
        BOOST_THROW_EXCEPTION(std::invalid_argument ("x.size() != simplex.size()"));
      Scalar sum(0.0); 
      Scalar max_x = maximum<Vector,Scalar>(x);
      for (typename Vector::size_type i = 0; i < x.size(); ++i)
        sum += (simplex[i] = exp(x[i]-max_x));
      for (typename Vector::size_type i = 0; i < x.size(); ++i)
        simplex[i] /= sum;
    }

    template <typename Vector>
    void inverse_softmax(const Vector& simplex, Vector& y) {
      using std::log;
      if(simplex.size() != y.size())
        BOOST_THROW_EXCEPTION(std::invalid_argument ("simplex.size() != y.size()"));
      for (size_t i = 0; i < simplex.size(); ++i)
        y[i] = log(simplex[i]);
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    log1p(T x) {
      using std::log;
      if (x < -1.0)
        BOOST_THROW_EXCEPTION(std::domain_error ("x can not be less than -1"));

      if (x > 1e-9 || x < -1e-9)
        return log(1.0 + x);     // direct, if distant from 1
      else if (x > 1e-16 || x < -1e-16)
        return x - 0.5 * x * x;  // 2nd order Taylor, if close to 1
      else
        return x;                // 1st order Taylor, if very close to 1
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    log1m(T x) {
      return log1p(-x);
    }

    namespace {
      const double LOG_PI_OVER_FOUR = std::log(boost::math::constants::pi<double>()) / 4.0;
    }

    template <typename T>
    inline typename boost::math::tools::promote_args<T>::type
    lmgamma(unsigned int k, T x) {
      typename boost::math::tools::promote_args<T>::type result 
        = k * (k - 1) * LOG_PI_OVER_FOUR;
      for (unsigned int j = 1; j <= k; ++j)
        result += lgamma(x + (1.0 - j) / 2.0);
      return result;
    }
      

    inline double if_else(bool c, double y_true, double y_false) {
      return c ? y_true : y_false;
    }

    template <typename T>
    inline T square(T x) {
      return x * x;
    }
    
    template <typename T_a, typename T_b>
    inline typename boost::math::tools::promote_args<T_a,T_b>::type
    multiply_log(T_a a, T_b b) {
      using std::log;
      if (b == 0.0 && a == 0.0)
        return 0.0;
      return a * log(b);
    }

    inline double log1p_exp(const double& a) {
      using std::exp;
      // like log_sum_exp below with b=0.0
      if (a > 0.0)
        return a + log1p(exp(-a));
      return log1p(exp(a));
    }

    inline double log_sum_exp(const double& a, const double& b) {
      using std::exp;
      if (a > b)
        return a + log1p(exp(b - a));
      return b + log1p(exp(a - b));
    }

    inline double ibeta(const double& a,
                        const double& b,
                        const double& x) {
      return boost::math::ibeta(a, b, x);
    }



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

    template <typename T1, typename T2, typename T3>
    inline double simple_var(double v, 
                             const T1& /*y1*/, const T1& /*dy1*/, 
                             const T2& /*y2*/, const T2& /*dy2*/,
                             const T3& /*y3*/, const T3& /*dy3*/) {
      return v;
    }

    template <typename T>
    inline double value_of(T x) {
      return static_cast<double>(x);
    }

    template <>
    inline double value_of<double>(double x) {
      return x; 
    }

    // CONSTANTS

    inline double pi() {
      return boost::math::constants::pi<double>();
    }

    inline double e() {
      return E;
    }

    inline double sqrt2() {
      return SQRT_2;
    }

    inline double log2() {
      return LOG_2;
    }

    inline double log10() {
      return LOG_10;
    }

    inline double positive_infinity() {
      return INFTY;
    }

    inline double negative_infinity() {
      return NEGATIVE_INFTY;
    }

    inline double not_a_number() {
      return NOT_A_NUMBER;
    }

    inline double epsilon() {
      return EPSILON;
    }

    inline double negative_epsilon() {
      return NEGATIVE_EPSILON;
    }

  }


}

#endif