agrad.hpp

Go to the documentation of this file.

#ifndef __STAN__AGRAD__AGRAD_HPP__
#define __STAN__AGRAD__AGRAD_HPP__

#include <cmath>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <vector>
#include <ostream>
#include <iostream>

#include <stan/memory/stack_alloc.hpp>

namespace stan {

  namespace agrad {

    class chainable;
    class vari;
    class var;

    // FIXME: manage all this as a single singleton (thread local)
    extern std::vector<chainable*> var_stack_; 
    extern memory::stack_alloc memalloc_;

    static void recover_memory();

    static void grad(chainable* vi);

    class chainable {

    public:

      chainable() { }

      ~chainable() { 
        // handled automatically
      }

      virtual void chain() {
      }

      virtual void init_dependent() {
      }

      virtual void set_zero_adjoint() {
      }
      

      static inline void* operator new(size_t nbytes) {
        return memalloc_.alloc(nbytes);
      }

    };

    class vari : public chainable {
    private:
      friend class var;

    public: 

      const double val_;

      double adj_;
      
      vari(const double x): 
        val_(x),
        adj_(0.0) {
        var_stack_.push_back(this);
      }

      ~vari() { 
        throw std::logic_error("vari destruction handled automatically");
      }

      virtual void init_dependent() {
        adj_ = 1.0;   // droot/droot = 1
      }

      virtual void set_zero_adjoint() {
        adj_ = 0.0;
      }

      friend std::ostream& operator<<(std::ostream& os, const vari* v) {
        return os << v << "  " << v->val_ << " : " << v->adj_;
      }

    };

    inline void print_stack(std::ostream& o) {
      o << "STACK, size=" << var_stack_.size() << std::endl;
      for (size_t i = 0; i < var_stack_.size(); ++i)
        o << i 
          << "  " << var_stack_[i]
          << "  " << (static_cast<vari*>(var_stack_[i]))->val_
          << " : " << (static_cast<vari*>(var_stack_[i]))->adj_
          << std::endl;
    }


    class var {
    public:

      // FIXME: doc what this is for
      typedef double Scalar;

      vari * vi_;

      bool is_uninitialized() {
        return (vi_ == static_cast<vari*>(0U));
      }

      explicit var(vari* vi)
      : vi_(vi) 
      {  }

       var() 
       : vi_(static_cast<vari*>(0U))
       { }

      var(bool b) :
        vi_(new vari(static_cast<double>(b))) {
      }

      var(char c) :
        vi_(new vari(static_cast<double>(c))) {
      }

      var(short n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned short n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(long int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned long int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned long long n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(long long n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(float x) :
        vi_(new vari(static_cast<double>(x))) {
      }

      var(double x) :
        vi_(new vari(x)) {
      }

      var(long double x) :
        vi_(new vari(static_cast<double>(x))) {
      }

      inline double val() const {
        return vi_->val_;
      }

      inline double adj() const {
        return vi_->adj_;
      }

      void grad(std::vector<var>& x,
                std::vector<double>& g) {
        stan::agrad::grad(vi_);
        g.resize(x.size());
        for (size_t i = 0; i < x.size(); ++i) 
          g[i] = x[i].vi_->adj_;
        recover_memory();
      }

      void grad() {
        stan::agrad::grad(vi_);
        recover_memory();
      }

      // POINTER OVERRIDES
      
      inline vari& operator*() {
        return *vi_;
      }

      inline vari* operator->() {
        return vi_;
      }

      // COMPOUND ASSIGNMENT OPERATORS
    
      inline var& operator+=(const var& b); 

      inline var& operator+=(const double b);

      inline var& operator-=(const var& b);

      inline var& operator-=(const double b);

      inline var& operator*=(const var& b);

      inline var& operator*=(const double b);

      inline var& operator/=(const var& b);

      inline var& operator/=(const double b);

      friend std::ostream& operator<<(std::ostream& os, const var& v) {
        return os << v.val() << ':' << v.adj();
      }
    };

    namespace {

      class op_v_vari : public vari {
      protected:
        vari* avi_;
      public:
        op_v_vari(double f, vari* avi) :
          vari(f),
          avi_(avi) {
        }
      };

      class op_vv_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
      public:
        op_vv_vari(double f, vari* avi, vari* bvi):
          vari(f),
          avi_(avi),
          bvi_(bvi) {
        }
      };

      class op_vd_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
      public:
        op_vd_vari(double f, vari* avi, double b) :
          vari(f),
          avi_(avi),
          bd_(b) {
        }
      };

      class op_dv_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
      public:
        op_dv_vari(double f, double a, vari* bvi) :
          vari(f),
          ad_(a),
          bvi_(bvi) {
        }
      };

      class op_vvv_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
        vari* cvi_;
      public:
        op_vvv_vari(double f, vari* avi, vari* bvi, vari* cvi) :
          vari(f),
          avi_(avi),
          bvi_(bvi),
          cvi_(cvi) {
        }
      };

      class op_vvd_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
        double cd_;
      public:
        op_vvd_vari(double f, vari* avi, vari* bvi, double c) :
          vari(f),
          avi_(avi),
          bvi_(bvi),
          cd_(c) {
        }
      };

      class op_vdv_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
        vari* cvi_;
      public:
        op_vdv_vari(double f, vari* avi, double b, vari* cvi) :
          vari(f),
          avi_(avi),
          bd_(b), 
          cvi_(cvi) {
        }
      };

      class op_vdd_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
        double cd_;
      public:
        op_vdd_vari(double f, vari* avi, double b, double c) :
          vari(f),
          avi_(avi),
          bd_(b), 
          cd_(c) {
        }
      };

      class op_dvv_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
        vari* cvi_;
      public:
        op_dvv_vari(double f, double a, vari* bvi, vari* cvi) :
          vari(f),
          ad_(a),
          bvi_(bvi),
          cvi_(cvi) {
        }
      };

      class op_dvd_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
        double cd_;
      public:
        op_dvd_vari(double f, double a, vari* bvi, double c) :
          vari(f),
          ad_(a),
          bvi_(bvi),
          cd_(c) {
        }
      };

      class op_ddv_vari : public vari {
      protected:
        double ad_;
        double bd_;
        vari* cvi_;
      public:
        op_ddv_vari(double f, double a, double b, vari* cvi) :
          vari(f),
          ad_(a),
          bd_(b),
          cvi_(cvi) {
        }
      };

      // FIXME: memory leak -- copy vector to local memory
      class op_vector_vari : public vari {
      protected:
        const size_t size_;
        vari** vis_;
      public:
        op_vector_vari(double f, const std::vector<stan::agrad::var>& vs) :
          vari(f),
          size_(vs.size()) {
          vis_ = (vari**) operator new(sizeof(vari*[vs.size()]));
          for (size_t i = 0; i < vs.size(); ++i)
            vis_[i] = vs[i].vi_;
        }
        vari* operator[](size_t n) const {
          return vis_[n];
        }
        size_t size() {
          return size_;
        }
      };

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


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

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

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

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

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

      class subtract_dv_vari : public op_dv_vari {
      public:
        subtract_dv_vari(double a, vari* bvi) :
          op_dv_vari(a - bvi->val_, a, bvi) {
        }
        void chain() {
          bvi_->adj_ -= adj_;
        }
      };

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

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

      // (a/b)' = a' * (1 / b) - b' * (a / [b * b])
      class divide_vv_vari : public op_vv_vari {
      public:
        divide_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ / bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_ / bvi_->val_;
          bvi_->adj_ -= adj_ * avi_->val_ / (bvi_->val_ * bvi_->val_);
        }
      };

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

      class divide_dv_vari : public op_dv_vari {
      public:
        divide_dv_vari(double a, vari* bvi) :
          op_dv_vari(a / bvi->val_, a, bvi) {
        }
        void chain() {
          bvi_->adj_ -= adj_ * ad_ / (bvi_->val_ * bvi_->val_);
        }
      };

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

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

      double LOG_10 = std::log(10.0);
    
      class log10_vari : public op_v_vari {
      public:
        const double exp_val_;
        log10_vari(vari* avi) :
          op_v_vari(std::log10(avi->val_),avi),
          exp_val_(avi->val_) {
        }
        void chain() {
          avi_->adj_ += adj_ / (LOG_10 * exp_val_);
        }
      };

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

      class pow_vv_vari : public op_vv_vari {
      public:
        pow_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::pow(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          if (avi_->val_ == 0.0) return; // partials zero, avoids 0 & log(0)
          avi_->adj_ += adj_ * bvi_->val_ * val_ / avi_->val_;
          bvi_->adj_ += adj_ * std::log(avi_->val_) * val_;
        }
      };

      class pow_vd_vari : public op_vd_vari {
      public:
        pow_vd_vari(vari* avi, double b) :
          op_vd_vari(std::pow(avi->val_,b),avi,b) {
        }
        void chain() {
          if (avi_->val_ == 0.0) return; // partials zero, avoids 0 & log(0)
          avi_->adj_ += adj_ * bd_ * val_ / avi_->val_;
        }
      };

      class pow_dv_vari : public op_dv_vari {
      public:
        pow_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::pow(a,bvi->val_),a,bvi) {
        }
        void chain() {
          if (ad_ == 0.0) return; // partials zero, avoids 0 & log(0)
          bvi_->adj_ += adj_ * std::log(ad_) * val_;
        }
      };

      class cos_vari : public op_v_vari {
      public:
        cos_vari(vari* avi) :
          op_v_vari(std::cos(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ * std::sin(avi_->val_);
        }
      };

      class sin_vari : public op_v_vari {
      public:
        sin_vari(vari* avi) :
          op_v_vari(std::sin(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::cos(avi_->val_);
        }
      };

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

      class acos_vari : public op_v_vari {
      public:
        acos_vari(vari* avi) :
          op_v_vari(std::acos(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ / std::sqrt(1.0 - (avi_->val_ * avi_->val_));
        }
      };

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

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

      class atan2_vv_vari : public op_vv_vari {
      public:
        atan2_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::atan2(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          double a_sq_plus_b_sq = (avi_->val_ * avi_->val_) + (bvi_->val_ * bvi_->val_);
          avi_->adj_ += bvi_->val_ / a_sq_plus_b_sq;
          bvi_->adj_ -= avi_->val_ / a_sq_plus_b_sq;
        }
      };

      class atan2_vd_vari : public op_vd_vari {
      public:
        atan2_vd_vari(vari* avi, double b) :
          op_vd_vari(std::atan2(avi->val_,b),avi,b) {
        }
        void chain() {
          double a_sq_plus_b_sq = (avi_->val_ * avi_->val_) + (bd_ * bd_);
          avi_->adj_ += bd_ / a_sq_plus_b_sq;
        }
      };

      class atan2_dv_vari : public op_dv_vari {
      public:
        atan2_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::atan2(a,bvi->val_),a,bvi) {
        }
        void chain() {
          double a_sq_plus_b_sq = (ad_ * ad_) + (bvi_->val_ * bvi_->val_);
          bvi_->adj_ -= ad_ / a_sq_plus_b_sq;
        }
      };

      class cosh_vari : public op_v_vari {
      public:
        cosh_vari(vari* avi) :
          op_v_vari(std::cosh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::sinh(avi_->val_);
        }
      };

      class sinh_vari : public op_v_vari {
      public:
        sinh_vari(vari* avi) :
          op_v_vari(std::sinh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::cosh(avi_->val_);
        }
      };

      class tanh_vari : public op_v_vari {
      public:
        tanh_vari(vari* avi) :
          op_v_vari(std::tanh(avi->val_),avi) {
        }
        void chain() {
          double cosh = std::cosh(avi_->val_);
          avi_->adj_ += adj_ / (cosh * cosh);
        }
      };


      class floor_vari : public vari {
      public:
        floor_vari(vari* avi) :
          vari(std::floor(avi->val_)) {
        }
      };

      class ceil_vari : public vari {
      public:
        ceil_vari(vari* avi) :
          vari(std::ceil(avi->val_)) {
        }
      };

      class fmod_vv_vari : public op_vv_vari {
      public:
        fmod_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::fmod(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          avi_->adj_ += adj_;
          bvi_->adj_ -= adj_ * static_cast<int>(avi_->val_ / bvi_->val_);
        }
      };

      class fmod_vd_vari : public op_v_vari {
      public:
        fmod_vd_vari(vari* avi, double b) :
          op_v_vari(std::fmod(avi->val_,b),avi) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };
      
      class fmod_dv_vari : public op_dv_vari {
      public:
        fmod_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::fmod(a,bvi->val_),a,bvi) {
        }
        void chain() {
          int d = static_cast<int>(ad_ / bvi_->val_);
          bvi_->adj_ -= adj_ * d;
        }
      };
    }


    // COMPARISON OPERATORS

    inline bool operator==(const var& a, const var& b) {
      return a.val() == b.val();
    }

    inline bool operator==(const var& a, const double b) {
      return a.val() == b;
    }
  
    inline bool operator==(const double a, const var& b) {
      return a == b.val();
    }

    inline bool operator!=(const var& a, const var& b) {
      return a.val() != b.val();
    }

    inline bool operator!=(const var& a, const double b) {
      return a.val() != b;
    }

    inline bool operator!=(const double a, const var& b) {
      return a != b.val();
    }

    inline bool operator<(const var& a, const var& b) {
      return a.val() < b.val();
    }

    inline bool operator<(const var& a, const double b) {
      return a.val() < b;
    }

    inline bool operator<(const double a, const var& b) {
      return a < b.val();
    }

    inline bool operator>(const var& a, const var& b) {
      return a.val() > b.val();
    }

    inline bool operator>(const var& a, const double b) {
      return a.val() > b;
    }

    inline bool operator>(const double a, const var& b) {
      return a > b.val();
    }

    inline bool operator<=(const var& a, const var& b) {
      return a.val() <= b.val();
    }

    inline bool operator<=(const var& a, const double b) {
      return a.val() <= b;
    }

    inline bool operator<=(const double a, const var& b) {
      return a <= b.val();
    }

    inline bool operator>=(const var& a, const var& b) {
      return a.val() >= b.val();
    }

    inline bool operator>=(const var& a, const double b) {
      return a.val() >= b;
    }

    inline bool operator>=(const double a, const var& b) {
      return a >= b.val();
    }

    // LOGICAL OPERATORS

    inline bool operator!(const var& a) {
      return !a.val();
    }


    // ARITHMETIC OPERATORS

    inline var operator+(const var& a) {
      return a;
    }

    inline var operator-(const var& a) {
      return var(new neg_vari(a.vi_));
    }

    inline var operator+(const var& a, const var& b) {    
      return var(new add_vv_vari(a.vi_,b.vi_));
    }


    inline var operator+(const var& a, const double b) {
      if (b == 0.0)
        return a;
      return var(new add_vd_vari(a.vi_,b));
    }

    inline var operator+(const double a, const var& b) {
      if (a == 0.0)
        return b;
      return var(new add_vd_vari(b.vi_,a)); // by symmetry
    }

    inline var operator-(const var& a, const var& b) {
      return var(new subtract_vv_vari(a.vi_,b.vi_));
    }

    inline var operator-(const var& a, const double b) {
      if (b == 0.0)
        return a;
      return var(new subtract_vd_vari(a.vi_,b));
    }

    inline var operator-(const double a, const var& b) {
      return var(new subtract_dv_vari(a,b.vi_));
    }

    inline var operator*(const var& a, const var& b) {
      return var(new multiply_vv_vari(a.vi_,b.vi_));
    }

    inline var operator*(const var& a, const double b) {
      if (b == 1.0)
        return a;
      return var(new multiply_vd_vari(a.vi_,b));
    }

    inline var operator*(const double a, const var& b) {
      if (a == 1.0)
        return b;
      return var(new multiply_vd_vari(b.vi_,a)); // by symmetry
    }

    inline var operator/(const var& a, const var& b) {
      return var(new divide_vv_vari(a.vi_,b.vi_));
    }

    inline var operator/(const var& a, const double b) {
      if (b == 1.0)
        return a;
      return var(new divide_vd_vari(a.vi_,b));
    }

    inline var operator/(const double a, const var& b) {
      return var(new divide_dv_vari(a,b.vi_));
    }

    inline var& operator++(var& a) {
      a.vi_ = new increment_vari(a.vi_);
      return a;
    }

    inline var operator++(var& a, int dummy) {
      var temp(a);
      a.vi_ = new increment_vari(a.vi_);
      return temp;
    }

    inline var& operator--(var& a) {
      a.vi_ = new decrement_vari(a.vi_);
      return a;
    }

    inline var operator--(var& a, int dummy) {
      var temp(a);
      a.vi_ = new decrement_vari(a.vi_);
      return temp;
    }

    // CMATH EXP AND LOG
  
    inline var exp(const var& a) {
      return var(new exp_vari(a.vi_));
    }

    inline var log(const var& a) {
      return var(new log_vari(a.vi_));
    }

    inline var log10(const var& a) {
      return var(new log10_vari(a.vi_));
    }


    // POWER FUNCTIONS

    inline var sqrt(const var& a) {
      return var(new sqrt_vari(a.vi_));
    }

    inline var pow(const var& base, const var& exponent) {
      return var(new pow_vv_vari(base.vi_,exponent.vi_));
    }
  
    inline var pow(const var& base, const double exponent) {
      if (exponent == 0.5)
        return sqrt(base);
      if (exponent == 1.0)
        return base;
      if (exponent == 2.0)
        return base * base; // FIXME: square() functionality from special_functions
      return var(new pow_vd_vari(base.vi_,exponent));
    }

    inline var pow(const double base, const var& exponent) {
      return var(new pow_dv_vari(base,exponent.vi_));
    }


    // TRIG FUNCTIONS

    inline var cos(const var& a) {
      return var(new cos_vari(a.vi_));
    }

    inline var sin(const var& a) {
      return var(new sin_vari(a.vi_));
    }

    inline var tan(const var& a) {
      return var(new tan_vari(a.vi_));
    }

    inline var acos(const var& a) {
      return var(new acos_vari(a.vi_));
    }

    inline var asin(const var& a) {
      return var(new asin_vari(a.vi_));
    }

    inline var atan(const var& a) {
      return var(new atan_vari(a.vi_));
    }

    inline var atan2(const var& a, const var& b) {
      return var(new atan2_vv_vari(a.vi_,b.vi_));
    }

    inline var atan2(const var& a, const double b) {
      return var(new atan2_vd_vari(a.vi_,b));
    }

    inline var atan2(const double a, const var& b) {
      return var(new atan2_dv_vari(a,b.vi_));
    }

    // HYPERBOLIC FUNCTIONS

    inline var cosh(const var& a) {
      return var(new cosh_vari(a.vi_));
    }

    inline var sinh(const var& a) {
      return var(new sinh_vari(a.vi_));
    }

    inline var tanh(const var& a) {
      return var(new tanh_vari(a.vi_));
    }
  

    // ROUNDING FUNCTIONS

    inline var fabs(const var& a) {
      // cut-and-paste from abs()
      if (a.val() > 0.0)
        return a;
      if (a.val() < 0.0)
        return var(new neg_vari(a.vi_));
      // FIXME:  is this right?  breaks connection to a
      return var(new vari(0.0));
    }

    inline var floor(const var& a) {
      return var(new floor_vari(a.vi_));
    }

    inline var ceil(const var& a) {
      return var(new ceil_vari(a.vi_));
    }

    inline var fmod(const var& a, const var& b) {
      return var(new fmod_vv_vari(a.vi_,b.vi_));
    }
  
    inline var fmod(const var& a, const double b) {
      return var(new fmod_vd_vari(a.vi_,b));
    }

    inline var fmod(const double a, const var& b) {
      return var(new fmod_dv_vari(a,b.vi_));
    }


    // STD LIB FUNCTIONS

    inline var abs(const var& a) {   
      // cut-and-paste from fabs()
      if (a.val() > 0.0)
        return a;
      if (a.val() < 0.0)
        return var(new neg_vari(a.vi_));
      // FIXME:  same as fabs() -- is this right?
      return var(new vari(0.0));
    }


    static void recover_memory() {
      var_stack_.resize(0);
      memalloc_.recover_all();
    }

    static void free_memory() {
      memalloc_.free_all();
    }

    static void grad(chainable* vi) {
      // old with subtle *2 bug
      // std::vector<chainable*>::iterator begin = var_stack_.begin();
      // std::vector<chainable*>::iterator it = var_stack_.end();
      // for (; (it >= begin) && (*it != vi); --it) ;

      std::vector<chainable*>::iterator begin = var_stack_.begin();
      std::vector<chainable*>::iterator it  = var_stack_.end();
      if (begin == it) return; // nothing on stack
      for (--it; (it >= begin) && (*it != vi); --it) ;

      vi->init_dependent(); 
      // propagate derivates for remaining vars
      for (; it >= begin; --it)
        (*it)->chain();
    }

    static void set_zero_all_adjoints() {
      for (size_t i = 0; i < var_stack_.size(); ++i)
        var_stack_[i]->set_zero_adjoint();
    }

    inline void jacobian(std::vector<var>& dependents,
                         std::vector<var>& independents,
                         std::vector<std::vector<double> >& jacobian) {
      jacobian.resize(dependents.size());
      for (size_t i = 0; i < dependents.size(); ++i) {
        jacobian[i].resize(independents.size());
        if (i > 0) 
          set_zero_all_adjoints();
        jacobian.push_back(std::vector<double>(0));
        grad(dependents[i].vi_);
        for (size_t j = 0; j < independents.size(); ++j)
          jacobian[i][j] = independents[j].adj();
      }
    }


    inline var& var::operator+=(const var& b) {
      vi_ = new add_vv_vari(vi_,b.vi_);
      return *this;
    }

    inline var& var::operator+=(const double b) {
      if (b == 0.0)
        return *this;
      vi_ = new add_vd_vari(vi_,b);
      return *this;
    }

    inline var& var::operator-=(const var& b) {
      vi_ = new subtract_vv_vari(vi_,b.vi_);
      return *this;
    }

    inline var& var::operator-=(const double b) {
      if (b == 0.0)
        return *this;
      vi_ = new subtract_vd_vari(vi_,b);
      return *this;
    }

    inline var& var::operator*=(const var& b) {
      vi_ = new multiply_vv_vari(vi_,b.vi_);
      return *this;
    }

    inline var& var::operator*=(const double b) {
      if (b == 1.0)
        return *this;
      vi_ = new multiply_vd_vari(vi_,b);
      return *this;
    }

    inline var& var::operator/=(const var& b) {
        vi_ = new divide_vv_vari(vi_,b.vi_);
        return *this;
      }

    inline var& var::operator/=(const double b) {
        if (b == 1.0)
          return *this;
        vi_ = new divide_vd_vari(vi_,b);
        return *this;
      }
  }

}


namespace std {

  template<> 
  struct numeric_limits<stan::agrad::var> {
    static const bool is_specialized = true;
    static stan::agrad::var min() { return numeric_limits<double>::min(); }
    static stan::agrad::var max() { return numeric_limits<double>::max(); }
    static const int digits = numeric_limits<double>::digits;
    static const int digits10 = numeric_limits<double>::digits10;
    static const bool is_signed = numeric_limits<double>::is_signed;
    static const bool is_integer = numeric_limits<double>::is_integer;
    static const bool is_exact = numeric_limits<double>::is_exact;
    static const int radix = numeric_limits<double>::radix;
    static stan::agrad::var epsilon() { return numeric_limits<double>::epsilon(); }
    static stan::agrad::var round_error() { return numeric_limits<double>::round_error(); }

    static const int  min_exponent = numeric_limits<double>::min_exponent;
    static const int  min_exponent10 = numeric_limits<double>::min_exponent10;
    static const int  max_exponent = numeric_limits<double>::max_exponent;
    static const int  max_exponent10 = numeric_limits<double>::max_exponent10;

    static const bool has_infinity = numeric_limits<double>::has_infinity;
    static const bool has_quiet_NaN = numeric_limits<double>::has_quiet_NaN;
    static const bool has_signaling_NaN = numeric_limits<double>::has_signaling_NaN;
    static const float_denorm_style has_denorm = numeric_limits<double>::has_denorm;
    static const bool has_denorm_loss = numeric_limits<double>::has_denorm_loss;
    static stan::agrad::var infinity() { return numeric_limits<double>::infinity(); }
    static stan::agrad::var quiet_NaN() { return numeric_limits<double>::quiet_NaN(); }
    static stan::agrad::var signaling_NaN() { return numeric_limits<double>::signaling_NaN(); }
    static stan::agrad::var denorm_min() { return numeric_limits<double>::denorm_min(); }

    static const bool is_iec559 = numeric_limits<double>::is_iec559;
    static const bool is_bounded = numeric_limits<double>::is_bounded;
    static const bool is_modulo = numeric_limits<double>::is_modulo;

    static const bool traps = numeric_limits<double>::traps;
    static const bool tinyness_before = numeric_limits<double>::tinyness_before;
    static const float_round_style round_style = numeric_limits<double>::round_style;
  };

  inline int isnan(const stan::agrad::var& a) {
    return isnan(a.val());
  }

  inline int isinf(const stan::agrad::var& a) {
    return isinf(a.val());
  }

}

#endif