Feature/issue 123 complex numbers with vars

make/libraries

@@ -78,8 +78,8 @@ $(BOOST_LIB)/mpi.so: $(BOOST)/user-config.jam
 $(BOOST_LIB)/libboost_serialization.so: $(BOOST_LIB)/mpi.so
 
 $(BOOST_LIB)/libboost_serialization.dylib: $(BOOST_LIB)/mpi.so
-	install_name_tool -add_rpath "$(BOOST_LIB_ABS)" "$(BOOST_LIB)/libboost_serialization.dylib"
-	install_name_tool -id @rpath/libboost_serialization.dylib "$(BOOST_LIB)/libboost_serialization.dylib"
+	install_name_tool -add_rpath "$(BOOST_LIB_ABS)" $(BOOST_LIB)/libboost_serialization.dylib
+	install_name_tool -id @rpath/libboost_serialization.dylib $(BOOST_LIB)/libboost_serialization.dylib
 
 $(BOOST_LIB)/libboost_mpi.so: $(BOOST_LIB)/mpi.so
 

stan/math/fwd/core/fvar.hpp

@@ -4,8 +4,8 @@
 #include <stan/math/prim/scal/meta/likely.hpp>
 #include <stan/math/prim/scal/fun/is_nan.hpp>
 #include <stan/math/fwd/scal/meta/ad_promotable.hpp>
-#include <boost/utility/enable_if.hpp>
 #include <ostream>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -97,13 +97,10 @@ struct fvar {
    * @tparam V type of value (must be assignable to the value and
    *   tangent type T)
    * @param[in] v value
-   * @param[in] dummy value given by default with enable-if
-   *   metaprogramming
    */
-  template <typename V>
-  fvar(const V& v,
-       typename boost::enable_if_c<ad_promotable<V, T>::value>::type* dummy = 0)
-      : val_(v), d_(0.0) {
+  template <typename V,
+            typename std::enable_if_t<ad_promotable<V, T>::value>* = nullptr>
+  fvar(const V& v) : val_(v), d_(0.0) {  // NOLINT(runtime/explicit)
     if (unlikely(is_nan(v)))
       d_ = v;
   }

stan/math/fwd/core/operator_division.hpp

@@ -46,6 +46,7 @@ inline fvar<T> operator/(double x1, const fvar<T>& x2) {
   // TODO(carpenter): store x1 / x2.val_ and reuse
   return fvar<T>(x1 / x2.val_, -x1 * x2.d_ / (x2.val_ * x2.val_));
 }
+
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/fwd/cplx.hpp

@@ -0,0 +1,7 @@
+#ifndef STAN_MATH_FWD_CPLX_HPP
+#define STAN_MATH_FWD_CPLX_HPP
+
+#include <stan/math/fwd/cplx/fvar.hpp>
+#include <stan/math/fwd/cplx/operator_division.hpp>
+
+#endif

stan/math/fwd/cplx/fvar.hpp

@@ -0,0 +1,21 @@
+#ifndef STAN_MATH_FWD_CPLX_FVAR_HPP
+#define STAN_MATH_FWD_CPLX_FVAR_HPP
+
+#include <stan/math/fwd/scal/meta/ad_promotable.hpp>
+#include <stan/math/fwd/scal/meta/is_fvar.hpp>
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/prim/cplx/complex.hpp>
+
+namespace std {
+
+/**
+ *  std::complex<fvar> inherits from stan::math::complex. Details
+ *  are provided there.
+ */
+template <class T>
+struct complex<stan::math::fvar<T>> : stan::math::complex<stan::math::fvar<T>> {
+  using stan::math::complex<stan::math::fvar<T>>::complex;
+};
+
+}  // namespace std
+#endif

stan/math/fwd/cplx/operator_division.hpp

@@ -0,0 +1,29 @@
+#ifndef STAN_MATH_FWD_CPLX_OPERATOR_DIVISION_HPP
+#define STAN_MATH_FWD_CPLX_OPERATOR_DIVISION_HPP
+
+#include <stan/math/fwd/cplx/fvar.hpp>
+#include <stan/math/fwd/scal/fun/pow.hpp>
+#include <stan/math/prim/cplx/operator_division.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the result of dividing the first argument by the second.
+ * This overload exists because libc++ uses logb and scalbn for
+ * complex division, which aren't defined for fvars.
+ *
+ * @tparam T type of complex fvar value and tangent
+ * @param t first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+template <class T>
+inline std::complex<fvar<T>> operator/(std::complex<fvar<T>> const& t,
+                                       std::complex<fvar<T>> const& u) {
+  return stan::math::operator_division(t, u);  // no recursion
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp

@@ -18,7 +18,7 @@ namespace math {
  * autodiff variable.
  */
 template <typename F, typename T>
-struct apply_scalar_unary<F, fvar<T> > {
+struct apply_scalar_unary<F, fvar<T>> {
   /**
    * Function return type, which is same as the argument type for
    * the function, <code>fvar&lt;T&gt;</code>.

stan/math/fwd/scal.hpp

@@ -1,6 +1,8 @@
 #ifndef STAN_MATH_FWD_SCAL_HPP
 #define STAN_MATH_FWD_SCAL_HPP
 
+#include <stan/math/fwd/cplx.hpp>
+
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/scal/meta/ad_promotable.hpp>
 #include <stan/math/fwd/scal/meta/is_fvar.hpp>
@@ -79,6 +81,7 @@
 #include <stan/math/fwd/scal/fun/Phi_approx.hpp>
 #include <stan/math/fwd/scal/fun/pow.hpp>
 #include <stan/math/fwd/scal/fun/primitive_value.hpp>
+#include <stan/math/fwd/scal/fun/proj.hpp>
 #include <stan/math/fwd/scal/fun/rising_factorial.hpp>
 #include <stan/math/fwd/scal/fun/round.hpp>
 #include <stan/math/fwd/scal/fun/sin.hpp>

stan/math/fwd/scal/fun/abs.hpp

@@ -1,7 +1,9 @@
 #ifndef STAN_MATH_FWD_SCAL_FUN_ABS_HPP
 #define STAN_MATH_FWD_SCAL_FUN_ABS_HPP
 
+#include <stan/math/fwd/cplx.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/scal/fun/hypot.hpp>
 #include <stan/math/fwd/scal/fun/value_of.hpp>
 #include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/abs.hpp>
@@ -25,4 +27,12 @@ inline fvar<T> abs(const fvar<T>& x) {
 
 }  // namespace math
 }  // namespace stan
+
+namespace std {
+// constrained complex overload to forward zeroing fvar to std::abs
+template <class T>
+inline stan::math::fvar<T> abs(std::complex<stan::math::fvar<T>> const& t) {
+  return abs(std::complex<stan::math::zeroing<stan::math::fvar<T>>>(t));
+}
+}  // namespace std
 #endif

stan/math/fwd/scal/fun/is_inf.hpp

@@ -20,6 +20,12 @@ inline int is_inf(const fvar<T>& x) {
   return is_inf(x.val());
 }
 
+// forwarding for ADL
+template <class T>
+inline auto isinf(const fvar<T>& a) {
+  return is_inf(a);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/fwd/scal/fun/is_nan.hpp

@@ -20,6 +20,12 @@ inline int is_nan(const fvar<T>& x) {
   return is_nan(x.val());
 }
 
+// forwarding for ADL
+template <class T>
+inline auto isnan(const fvar<T>& a) {
+  return is_nan(a);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/fwd/scal/fun/proj.hpp

@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_FWD_SCAL_FUN_PROJ_HPP
+#define STAN_MATH_FWD_SCAL_FUN_PROJ_HPP
+
+#include <stan/math/fwd/cplx.hpp>
+#include <stan/math/fwd/scal/fun/is_inf.hpp>
+#include <stan/math/prim/scal/fun/proj.hpp>
+
+namespace std {
+
+/**
+ * Return a projection of a complex number onto
+ * the Riemann sphere.
+ *
+ * This overloads an erroneous definition in
+ * libstdc++ for general types.
+ *
+ * @tparam T auto diff variable type
+ * @param[in] t Variable input.
+ * @return projection onto Riemann sphere
+ */
+template <class T>
+inline std::complex<stan::math::fvar<T>> proj(
+    std::complex<stan::math::fvar<T>> const& t) {
+  return stan::math::proj(t);
+}
+
+}  // namespace std
+#endif

stan/math/prim/cplx.hpp

@@ -0,0 +1,11 @@
+#ifndef STAN_MATH_PRIM_CPLX_HPP
+#define STAN_MATH_PRIM_CPLX_HPP
+
+#include <stan/math/prim/cplx/complex.hpp>
+#include <stan/math/prim/cplx/operator_addition.hpp>
+#include <stan/math/prim/cplx/operator_division.hpp>
+#include <stan/math/prim/cplx/operator_multiplication.hpp>
+#include <stan/math/prim/cplx/operator_subtraction.hpp>
+#include <stan/math/prim/cplx/zeroing.hpp>
+
+#endif

stan/math/prim/cplx/complex.hpp

@@ -0,0 +1,57 @@
+#ifndef STAN_MATH_PRIM_CPLX_COMPLEX_HPP
+#define STAN_MATH_PRIM_CPLX_COMPLEX_HPP
+
+#include <stan/math/prim/cplx/zeroing.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ *  Complex class that forwards the interface of std::complex, brining
+ *  stan::math's namespace into ADL for functions called on data types
+ *  that inherit from this class, as well as allowing template
+ *  specializations of std::complex to indirectly work with different
+ *  underlying data types. For example, std::complex<var> specializes
+ *  the std::complex<T> template to inherit from
+ *  stan::math::complex<T>, which in turn inherits from
+ *  std::complex<zeroing<T>>. This causes stan::math to be pulled into
+ *  ADL on operations involving std::complex<var>, to handle, e.g.
+ *  division operations that the base std::complex doesn't have
+ *  defined usably for var, and, since stan::math::zeroing<var> zero
+ *  initializes itself, will also correctly work with the remaining
+ *  algorithms in std::complex<T> that require the expression T() to
+ *  be 0.
+ */
+template <class T>
+struct complex : std::complex<zeroing<T>> {
+  using std::complex<zeroing<T>>::complex;
+  // NOLINTNEXTLINE(runtime/explicit)
+  complex(const std::complex<zeroing<T>>& t)
+      : std::complex<zeroing<T>>(t) {}  ///< allow promotion
+  /**
+   * Without this converting ctor, copy initialization of
+   * std::complex<zeroing<(f)var>>> from an (f)var fails, because the
+   * default implementation in the std::complex template requires the
+   * argument type to match the template type. This is the same reason
+   * why std::complex<double> can't be multiplied by an int.
+   *
+   * tparam R type of real argument. Won't fire on non-arithmetic types.
+   * tparam I type of imaginary argument.
+   * param[in] real, the pure real component of the complex number.
+   * param[in] imag, the pure imaginary component of the complex number
+   */
+  template <class R = double, class I = R,
+            std::enable_if_t<is_arith_like<R>::value>* = nullptr>
+  // NOLINTNEXTLINE(runtime/explicit)
+  complex(R const& real = 0, I const& imag = 0)
+      : std::complex<zeroing<T>>(zeroing<T>(real), zeroing<T>(imag)) {}
+  // downcast
+  operator std::complex<T>() const {
+    return *static_cast<std::complex<T> const*>(this);
+  }
+};
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/cplx/operator_addition.hpp

@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_CPLX_OPERATOR_ADDITION_HPP
+#define STAN_MATH_PRIM_CPLX_OPERATOR_ADDITION_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the sum of the specified arguments
+ *
+ * @tparam T type of complex first argument
+ * @tparam U type of second argument
+ * @param t complex first argument
+ * @param u second argument
+ * @return sum
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && (is_complex<U>::value
+                               || is_arith_like<U>::value)>* = nullptr>
+inline auto operator+(std::complex<T> const& t, U const& u) {
+  return complex_promote<U>(t) += u;
+}
+
+/**
+ * Return the sum of the specified arguments
+ *
+ * @tparam T type of first argument
+ * @tparam U type of complex second argument
+ * @param t first argument
+ * @param u complex second argument
+ * @return sum
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && is_arith_like<T>::value>* = nullptr>
+inline auto operator+(T const& t, std::complex<U> const& u) {
+  return complex_promote<U>(t) += u;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif