autodiff.h

Go to the documentation of this file.

// clang-format off
 
#ifndef __AUTODIFF_H
#define __AUTODIFF_H
 
#ifndef EIGEN_DONT_PARALLELIZE
#define EIGEN_DONT_PARALLELIZE
#endif
 
#include <Eigen/Core>
#include <cmath>
#include <stdexcept>
#include <limits>
 
struct DiffScalarBase
{
    // ======================================================================
    // ======================================================================
 
    static inline void setVariableCount(size_t value)
    {
        m_variableCount = value;
    }
 
    static inline size_t getVariableCount()
    {
        return m_variableCount;
    }
 
    // ======================================================================
 
    // #ifdef WIN32
    // static __declspec(thread) size_t m_variableCount;
    // #else
    // static __thread size_t m_variableCount;
    // #endif
    static thread_local size_t m_variableCount;
};
 
// #ifdef WIN32
// #define DECLARE_DIFFSCALAR_BASE()
//  __declspec(thread) size_t DiffScalarBase::m_variableCount = 0
// #else
// #define DECLARE_DIFFSCALAR_BASE()
//  __thread size_t DiffScalarBase::m_variableCount = 0
// #endif
 
#define DECLARE_DIFFSCALAR_BASE() \
    thread_local size_t DiffScalarBase::m_variableCount = 0
 
template <typename _Scalar, typename _Gradient = Eigen::Matrix<_Scalar, Eigen::Dynamic, 1>>
struct DScalar1 : public DiffScalarBase
{
public:
    typedef _Scalar Scalar;
    typedef _Gradient Gradient;
    typedef Eigen::Matrix<DScalar1, 2, 1> DVector2;
    typedef Eigen::Matrix<DScalar1, 3, 1> DVector3;
 
    // ======================================================================
    // ======================================================================
 
    explicit DScalar1(Scalar value_ = (Scalar)0) : value(value_)
    {
        size_t variableCount = getVariableCount();
        grad.resize(variableCount);
        grad.setZero();
    }
 
    DScalar1(size_t index, const Scalar &value_)
        : value(value_)
    {
        size_t variableCount = getVariableCount();
        grad.resize(variableCount);
        grad.setZero();
        grad(index) = 1;
    }
 
    DScalar1(Scalar value_, const Gradient &grad_)
        : value(value_), grad(grad_) {}
 
    DScalar1(const DScalar1 &s)
        : value(s.value), grad(s.grad) {}
 
    inline const Scalar &getValue() const { return value; }
    inline const Gradient &getGradient() const { return grad; }
 
    // ======================================================================
    // ======================================================================
    friend DScalar1 operator+(const DScalar1 &lhs, const DScalar1 &rhs)
    {
        return DScalar1(lhs.value + rhs.value, lhs.grad + rhs.grad);
    }
 
    friend DScalar1 operator+(const DScalar1 &lhs, const Scalar &rhs)
    {
        return DScalar1(lhs.value + rhs, lhs.grad);
    }
 
    friend DScalar1 operator+(const Scalar &lhs, const DScalar1 &rhs)
    {
        return DScalar1(rhs.value + lhs, rhs.grad);
    }
 
    inline DScalar1 &operator+=(const DScalar1 &s)
    {
        value += s.value;
        grad += s.grad;
        return *this;
    }
 
    inline DScalar1 &operator+=(const Scalar &v)
    {
        value += v;
        return *this;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    friend DScalar1 operator-(const DScalar1 &lhs, const DScalar1 &rhs)
    {
        return DScalar1(lhs.value - rhs.value, lhs.grad - rhs.grad);
    }
 
    friend DScalar1 operator-(const DScalar1 &lhs, const Scalar &rhs)
    {
        return DScalar1(lhs.value - rhs, lhs.grad);
    }
 
    friend DScalar1 operator-(const Scalar &lhs, const DScalar1 &rhs)
    {
        return DScalar1(lhs - rhs.value, -rhs.grad);
    }
 
    friend DScalar1 operator-(const DScalar1 &s)
    {
        return DScalar1(-s.value, -s.grad);
    }
 
    inline DScalar1 &operator-=(const DScalar1 &s)
    {
        value -= s.value;
        grad -= s.grad;
        return *this;
    }
 
    inline DScalar1 &operator-=(const Scalar &v)
    {
        value -= v;
        return *this;
    }
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
    friend DScalar1 operator/(const DScalar1 &lhs, const Scalar &rhs)
    {
        if (rhs == 0)
            throw std::runtime_error("DScalar1: Division by zero!");
        Scalar inv = 1.0f / rhs;
        return DScalar1(lhs.value * inv, lhs.grad * inv);
    }
 
    friend DScalar1 operator/(const Scalar &lhs, const DScalar1 &rhs)
    {
        return lhs * inverse(rhs);
    }
 
    friend DScalar1 operator/(const DScalar1 &lhs, const DScalar1 &rhs)
    {
        return lhs * inverse(rhs);
    }
 
    friend DScalar1 inverse(const DScalar1 &s)
    {
        Scalar valueSqr = s.value * s.value,
               invValueSqr = (Scalar)1 / valueSqr;
 
        // vn = 1/v, Dvn = -1/(v^2) Dv
        return DScalar1((Scalar)1 / s.value, s.grad * -invValueSqr);
    }
 
    inline DScalar1 &operator/=(const Scalar &v)
    {
        value /= v;
        grad /= v;
        return *this;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
    inline friend DScalar1 operator*(const DScalar1 &lhs, const Scalar &rhs)
    {
        return DScalar1(lhs.value * rhs, lhs.grad * rhs);
    }
 
    inline friend DScalar1 operator*(const Scalar &lhs, const DScalar1 &rhs)
    {
        return DScalar1(rhs.value * lhs, rhs.grad * lhs);
    }
 
    inline friend DScalar1 operator*(const DScalar1 &lhs, const DScalar1 &rhs)
    {
        // Product rule
        return DScalar1(lhs.value * rhs.value,
                        rhs.grad * lhs.value + lhs.grad * rhs.value);
    }
 
    inline DScalar1 &operator*=(const Scalar &v)
    {
        value *= v;
        grad *= v;
        return *this;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    friend DScalar1 abs(const DScalar1 &s)
    {
        return s.value < 0?  -s: s;
    }
 
    friend DScalar1 sqrt(const DScalar1 &s)
    {
        Scalar sqrtVal = std::sqrt(s.value),
               temp = (Scalar)1 / ((Scalar)2 * sqrtVal);
 
        // vn = sqrt(v)
        // Dvn = 1/(2 sqrt(v)) Dv
        return DScalar1(sqrtVal, s.grad * temp);
    }
 
    friend DScalar1 pow(const DScalar1 &s, const Scalar &a)
    {
        Scalar powVal = std::pow(s.value, a),
               temp = a * std::pow(s.value, a - 1);
        // vn = v ^ a, Dvn = a*v^(a-1) * Dv
        return DScalar1(powVal, s.grad * temp);
    }
 
    friend DScalar1 exp(const DScalar1 &s)
    {
        Scalar expVal = std::exp(s.value);
 
        // vn = exp(v), Dvn = exp(v) * Dv
        return DScalar1(expVal, s.grad * expVal);
    }
 
    friend DScalar1 log(const DScalar1 &s)
    {
        Scalar logVal = std::log(s.value);
 
        // vn = log(v), Dvn = Dv / v
        return DScalar1(logVal, s.grad / s.value);
    }
 
    friend DScalar1 sin(const DScalar1 &s)
    {
        // vn = sin(v), Dvn = cos(v) * Dv
        return DScalar1(std::sin(s.value), s.grad * std::cos(s.value));
    }
 
    friend DScalar1 cos(const DScalar1 &s)
    {
        // vn = cos(v), Dvn = -sin(v) * Dv
        return DScalar1(std::cos(s.value), s.grad * -std::sin(s.value));
    }
 
    friend DScalar1 acos(const DScalar1 &s)
    {
        if (std::abs(s.value) >= 1)
            throw std::runtime_error("acos: Expected a value in (-1, 1)");
 
        Scalar temp = -std::sqrt((Scalar)1 - s.value * s.value);
 
        // vn = acos(v), Dvn = -1/sqrt(1-v^2) * Dv
        return DScalar1(std::acos(s.value),
                        s.grad * ((Scalar)1 / temp));
    }
 
    friend DScalar1 asin(const DScalar1 &s)
    {
        if (std::abs(s.value) >= 1)
            throw std::runtime_error("asin: Expected a value in (-1, 1)");
 
        Scalar temp = std::sqrt((Scalar)1 - s.value * s.value);
 
        // vn = asin(v), Dvn = 1/sqrt(1-v^2) * Dv
        return DScalar1(std::asin(s.value),
                        s.grad * ((Scalar)1 / temp));
    }
 
    friend DScalar1 atan2(const DScalar1 &y, const DScalar1 &x)
    {
        Scalar denom = x.value * x.value + y.value * y.value;
 
        // vn = atan2(y, x), Dvn = (x*Dy - y*Dx) / (x^2 + y^2)
        return DScalar1(std::atan2(y.value, x.value),
                        y.grad * (x.value / denom) - x.grad * (y.value / denom));
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    inline void operator=(const DScalar1 &s)
    {
        value = s.value;
        grad = s.grad;
    }
    inline void operator=(const Scalar &v)
    {
        value = v;
        grad.setZero();
    }
    inline bool operator<(const DScalar1 &s) const { return value < s.value; }
    inline bool operator<=(const DScalar1 &s) const { return value <= s.value; }
    inline bool operator>(const DScalar1 &s) const { return value > s.value; }
    inline bool operator>=(const DScalar1 &s) const { return value >= s.value; }
    inline bool operator<(const Scalar &s) const { return value < s; }
    inline bool operator<=(const Scalar &s) const { return value <= s; }
    inline bool operator>(const Scalar &s) const { return value > s; }
    inline bool operator>=(const Scalar &s) const { return value >= s; }
    inline bool operator==(const Scalar &s) const { return value == s; }
    inline bool operator!=(const Scalar &s) const { return value != s; }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    static inline DVector2 vector(const Eigen::Matrix<Scalar, 2, 1> &v)
    {
        return DVector2(DScalar1(v.x()), DScalar1(v.y()));
    }
 
    static inline DVector3 vector(const Eigen::Matrix<Scalar, 3, 1> &v)
    {
        return DVector3(DScalar1(v.x()), DScalar1(v.y()), DScalar1(v.z()));
    }
 
#if defined(__MITSUBA_MITSUBA_H_) /* Mitsuba-specific */
    static inline DVector2 vector(const mitsuba::TVector2<Scalar> &v)
    {
        return DVector2(DScalar1(v.x), DScalar1(v.y));
    }
 
    static inline DVector2 vector(const mitsuba::TPoint2<Scalar> &p)
    {
        return DVector2(DScalar1(p.x), DScalar1(p.y));
    }
 
    static inline DVector3 vector(const mitsuba::TVector3<Scalar> &v)
    {
        return DVector3(DScalar1(v.x), DScalar1(v.y), DScalar1(v.z));
    }
 
    static inline DVector3 vector(const mitsuba::TPoint3<Scalar> &p)
    {
        return DVector3(DScalar1(p.x), DScalar1(p.y), DScalar1(p.z));
    }
#endif
 
    // ======================================================================
protected:
    Scalar value;
    Gradient grad;
};
 
template <typename Scalar, typename VecType>
std::ostream &operator<<(std::ostream &out, const DScalar1<Scalar, VecType> &s)
{
    out << "[" << s.getValue()
        << ", grad=" << s.getGradient().format(Eigen::IOFormat(4, 1, ", ", "; ", "", "", "[", "]"))
        << "]";
    return out;
}
 
template <typename _Scalar, typename _Gradient = Eigen::Matrix<_Scalar, Eigen::Dynamic, 1>,
          typename _Hessian = Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic>>
struct DScalar2 : public DiffScalarBase
{
public:
    typedef _Scalar Scalar;
    typedef _Gradient Gradient;
    typedef _Hessian Hessian;
    typedef Eigen::Matrix<DScalar2, 2, 1> DVector2;
    typedef Eigen::Matrix<DScalar2, 3, 1> DVector3;
 
    // ======================================================================
    // ======================================================================
 
    explicit DScalar2(Scalar value_ = (Scalar)0) : value(value_)
    {
        size_t variableCount = getVariableCount();
 
        grad.resize(variableCount);
        grad.setZero();
        hess.resize(variableCount, variableCount);
        hess.setZero();
    }
 
    DScalar2(size_t index, const Scalar &value_)
        : value(value_)
    {
        size_t variableCount = getVariableCount();
 
        grad.resize(variableCount);
        grad.setZero();
        grad(index) = 1;
        hess.resize(variableCount, variableCount);
        hess.setZero();
    }
 
    DScalar2(Scalar value_, const Gradient &grad_, const Hessian &hess_)
        : value(value_), grad(grad_), hess(hess_) {}
 
    DScalar2(const DScalar2 &s)
        : value(s.value), grad(s.grad), hess(s.hess) {}
 
    inline const Scalar &getValue() const { return value; }
    inline const Gradient &getGradient() const { return grad; }
    inline const Hessian &getHessian() const { return hess; }
 
    // ======================================================================
    // ======================================================================
    friend DScalar2 operator+(const DScalar2 &lhs, const DScalar2 &rhs)
    {
        return DScalar2(lhs.value + rhs.value,
                        lhs.grad + rhs.grad, lhs.hess + rhs.hess);
    }
 
    friend DScalar2 operator+(const DScalar2 &lhs, const Scalar &rhs)
    {
        return DScalar2(lhs.value + rhs, lhs.grad, lhs.hess);
    }
 
    friend DScalar2 operator+(const Scalar &lhs, const DScalar2 &rhs)
    {
        return DScalar2(rhs.value + lhs, rhs.grad, rhs.hess);
    }
 
    inline DScalar2 &operator+=(const DScalar2 &s)
    {
        value += s.value;
        grad += s.grad;
        hess += s.hess;
        return *this;
    }
 
    inline DScalar2 &operator+=(const Scalar &v)
    {
        value += v;
        return *this;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    friend DScalar2 operator-(const DScalar2 &lhs, const DScalar2 &rhs)
    {
        return DScalar2(lhs.value - rhs.value, lhs.grad - rhs.grad, lhs.hess - rhs.hess);
    }
 
    friend DScalar2 operator-(const DScalar2 &lhs, const Scalar &rhs)
    {
        return DScalar2(lhs.value - rhs, lhs.grad, lhs.hess);
    }
 
    friend DScalar2 operator-(const Scalar &lhs, const DScalar2 &rhs)
    {
        return DScalar2(lhs - rhs.value, -rhs.grad, -rhs.hess);
    }
 
    friend DScalar2 operator-(const DScalar2 &s)
    {
        return DScalar2(-s.value, -s.grad, -s.hess);
    }
 
    inline DScalar2 &operator-=(const DScalar2 &s)
    {
        value -= s.value;
        grad -= s.grad;
        hess -= s.hess;
        return *this;
    }
 
    inline DScalar2 &operator-=(const Scalar &v)
    {
        value -= v;
        return *this;
    }
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
    friend DScalar2 operator/(const DScalar2 &lhs, const Scalar &rhs)
    {
        if (rhs == 0)
            throw std::runtime_error("DScalar2: Division by zero!");
        Scalar inv = 1.0f / rhs;
        return DScalar2(lhs.value * inv, lhs.grad * inv, lhs.hess * inv);
    }
 
    friend DScalar2 operator/(const Scalar &lhs, const DScalar2 &rhs)
    {
        return lhs * inverse(rhs);
    }
 
    friend DScalar2 operator/(const DScalar2 &lhs, const DScalar2 &rhs)
    {
        return lhs * inverse(rhs);
    }
 
    friend DScalar2 inverse(const DScalar2 &s)
    {
        Scalar valueSqr = s.value * s.value,
               valueCub = valueSqr * s.value,
               invValueSqr = (Scalar)1 / valueSqr;
 
        // vn = 1/v
        DScalar2 result((Scalar)1 / s.value);
 
        // Dvn = -1/(v^2) Dv
        result.grad = s.grad * -invValueSqr;
 
        // D^2vn = -1/(v^2) D^2v + 2/(v^3) Dv Dv^T
        result.hess = s.hess * -invValueSqr;
        result.hess += s.grad * s.grad.transpose()
                       * ((Scalar)2 / valueCub);
 
        return result;
    }
 
    inline DScalar2 &operator/=(const Scalar &v)
    {
        value /= v;
        grad /= v;
        hess /= v;
        return *this;
    }
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
    friend DScalar2 operator*(const DScalar2 &lhs, const Scalar &rhs)
    {
        return DScalar2(lhs.value * rhs, lhs.grad * rhs, lhs.hess * rhs);
    }
 
    friend DScalar2 operator*(const Scalar &lhs, const DScalar2 &rhs)
    {
        return DScalar2(rhs.value * lhs, rhs.grad * lhs, rhs.hess * lhs);
    }
 
    friend DScalar2 operator*(const DScalar2 &lhs, const DScalar2 &rhs)
    {
        DScalar2 result(lhs.value * rhs.value);
 
        result.grad = rhs.grad * lhs.value + lhs.grad * rhs.value;
 
        // (i,j) = g*F_xixj + g*G_xixj + F_xi*G_xj + F_xj*G_xi
        result.hess = rhs.hess * lhs.value;
        result.hess += lhs.hess * rhs.value;
        result.hess += lhs.grad * rhs.grad.transpose();
        result.hess += rhs.grad * lhs.grad.transpose();
 
        return result;
    }
 
    inline DScalar2 &operator*=(const Scalar &v)
    {
        value *= v;
        grad *= v;
        hess *= v;
        return *this;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    friend DScalar2 abs(const DScalar2 &s)
    {
        return s.value < 0?  -s: s;
    }
 
    friend DScalar2 sqrt(const DScalar2 &s)
    {
        Scalar sqrtVal = std::sqrt(s.value),
               temp = (Scalar)1 / ((Scalar)2 * sqrtVal);
 
        // vn = sqrt(v)
        DScalar2 result(sqrtVal);
 
        // Dvn = 1/(2 sqrt(v)) Dv
        result.grad = s.grad * temp;
 
        // D^2vn = 1/(2 sqrt(v)) D^2v - 1/(4 v*sqrt(v)) Dv Dv^T
        result.hess = s.hess * temp;
        result.hess += s.grad * s.grad.transpose()
                       * (-(Scalar)1 / ((Scalar)4 * s.value * sqrtVal));
 
        return result;
    }
 
    friend DScalar2 pow(const DScalar2 &s, const Scalar &a)
    {
        Scalar powVal = std::pow(s.value, a),
               temp = a * std::pow(s.value, a - 1);
        // vn = v ^ a
        DScalar2 result(powVal);
 
        // Dvn = a*v^(a-1) * Dv
        result.grad = s.grad * temp;
 
        // D^2vn = a*v^(a-1) D^2v - 1/(4 v*sqrt(v)) Dv Dv^T
        result.hess = s.hess * temp;
        result.hess += s.grad * s.grad.transpose()
                       * (a * (a - 1) * std::pow(s.value, a - 2));
 
        return result;
    }
 
    friend DScalar2 exp(const DScalar2 &s)
    {
        Scalar expVal = std::exp(s.value);
 
        // vn = exp(v)
        DScalar2 result(expVal);
 
        // Dvn = exp(v) * Dv
        result.grad = s.grad * expVal;
 
        // D^2vn = exp(v) * Dv*Dv^T + exp(v) * D^2v
        result.hess = (s.grad * s.grad.transpose()
                       + s.hess)
                      * expVal;
 
        return result;
    }
 
    friend DScalar2 log(const DScalar2 &s)
    {
        Scalar logVal = std::log(s.value);
 
        // vn = log(v)
        DScalar2 result(logVal);
 
        // Dvn = Dv / v
        result.grad = s.grad / s.value;
 
        // D^2vn = (v*D^2v - Dv*Dv^T)/(v^2)
        result.hess = s.hess / s.value - (s.grad * s.grad.transpose() / (s.value * s.value));
 
        return result;
    }
 
    friend DScalar2 sin(const DScalar2 &s)
    {
        Scalar sinVal = std::sin(s.value),
               cosVal = std::cos(s.value);
 
        // vn = sin(v)
        DScalar2 result(sinVal);
 
        // Dvn = cos(v) * Dv
        result.grad = s.grad * cosVal;
 
        // D^2vn = -sin(v) * Dv*Dv^T + cos(v) * Dv^2
        result.hess = s.hess * cosVal;
        result.hess += s.grad * s.grad.transpose() * -sinVal;
 
        return result;
    }
 
    friend DScalar2 cos(const DScalar2 &s)
    {
        Scalar sinVal = std::sin(s.value),
               cosVal = std::cos(s.value);
        // vn = cos(v)
        DScalar2 result(cosVal);
 
        // Dvn = -sin(v) * Dv
        result.grad = s.grad * -sinVal;
 
        // D^2vn = -cos(v) * Dv*Dv^T - sin(v) * Dv^2
        result.hess = s.hess * -sinVal;
        result.hess += s.grad * s.grad.transpose() * -cosVal;
 
        return result;
    }
 
    friend DScalar2 acos(const DScalar2 &s)
    {
        if (std::abs(s.value) >= 1)
            throw std::runtime_error("acos: Expected a value in (-1, 1)");
 
        Scalar temp = -std::sqrt((Scalar)1 - s.value * s.value);
 
        // vn = acos(v)
        DScalar2 result(std::acos(s.value));
 
        // Dvn = -1/sqrt(1-v^2) * Dv
        result.grad = s.grad * ((Scalar)1 / temp);
 
        // D^2vn = -1/sqrt(1-v^2) * D^2v - v/[(1-v^2)^(3/2)] * Dv*Dv^T
        result.hess = s.hess * ((Scalar)1 / temp);
        result.hess += s.grad * s.grad.transpose()
                       * s.value / (temp * temp * temp);
 
        return result;
    }
 
    friend DScalar2 asin(const DScalar2 &s)
    {
        if (std::abs(s.value) >= 1)
            throw std::runtime_error("asin: Expected a value in (-1, 1)");
 
        Scalar temp = std::sqrt((Scalar)1 - s.value * s.value);
 
        // vn = asin(v)
        DScalar2 result(std::asin(s.value));
 
        // Dvn = 1/sqrt(1-v^2) * Dv
        result.grad = s.grad * ((Scalar)1 / temp);
 
        // D^2vn = 1/sqrt(1-v*v) * D^2v + v/[(1-v^2)^(3/2)] * Dv*Dv^T
        result.hess = s.hess * ((Scalar)1 / temp);
        result.hess += s.grad * s.grad.transpose()
                       * s.value / (temp * temp * temp);
 
        return result;
    }
 
    friend DScalar2 atan2(const DScalar2 &y, const DScalar2 &x)
    {
        // vn = atan2(y, x)
        DScalar2 result(std::atan2(y.value, x.value));
 
        // Dvn = (x*Dy - y*Dx) / (x^2 + y^2)
        Scalar denom = x.value * x.value + y.value * y.value,
               denomSqr = denom * denom;
        result.grad = y.grad * (x.value / denom)
                      - x.grad * (y.value / denom);
 
        // D^2vn = (Dy*Dx^T + xD^2y - Dx*Dy^T - yD^2x) / (x^2+y^2)
        //    - [(x*Dy - y*Dx) * (2*x*Dx + 2*y*Dy)^T] / (x^2+y^2)^2
        result.hess = (y.hess * x.value
                       + y.grad * x.grad.transpose()
                       - x.hess * y.value
                       - x.grad * y.grad.transpose())
                      / denom;
 
        result.hess -=
            (y.grad * (x.value / denomSqr) - x.grad * (y.value / denomSqr)) * (x.grad * ((Scalar)2 * x.value) + y.grad * ((Scalar)2 * y.value)).transpose();
 
        return result;
    }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
    inline void operator=(const DScalar2 &s)
    {
        value = s.value;
        grad = s.grad;
        hess = s.hess;
    }
    inline void operator=(const Scalar &v)
    {
        value = v;
        grad = Gradient::Zero(getVariableCount());
        hess = Hessian::Zero(getVariableCount(), getVariableCount());
    }
    inline bool operator<(const DScalar2 &s) const { return value < s.value; }
    inline bool operator<=(const DScalar2 &s) const { return value <= s.value; }
    inline bool operator>(const DScalar2 &s) const { return value > s.value; }
    inline bool operator>=(const DScalar2 &s) const { return value >= s.value; }
    inline bool operator!=(const DScalar2 &s) const { return value != s.value; }
 
    inline bool operator<(const Scalar &s) const { return value < s; }
    inline bool operator<=(const Scalar &s) const { return value <= s; }
    inline bool operator>(const Scalar &s) const { return value > s; }
    inline bool operator>=(const Scalar &s) const { return value >= s; }
    inline bool operator==(const Scalar &s) const { return value == s; }
    inline bool operator!=(const Scalar &s) const { return value != s; }
 
    // ======================================================================
 
    // ======================================================================
    // ======================================================================
 
#if defined(__MITSUBA_MITSUBA_H_) /* Mitsuba-specific */
    static inline DVector2 vector(const mitsuba::TVector2<Scalar> &v)
    {
        return DVector2(DScalar2(v.x), DScalar2(v.y));
    }
 
    static inline DVector2 vector(const mitsuba::TPoint2<Scalar> &p)
    {
        return DVector2(DScalar2(p.x), DScalar2(p.y));
    }
 
    static inline DVector3 vector(const mitsuba::TVector3<Scalar> &v)
    {
        return DVector3(DScalar2(v.x), DScalar2(v.y), DScalar2(v.z));
    }
 
    static inline DVector3 vector(const mitsuba::TPoint3<Scalar> &p)
    {
        return DVector3(DScalar2(p.x), DScalar2(p.y), DScalar2(p.z));
    }
#endif
 
    static inline DVector2 vector(const Eigen::Matrix<Scalar, 2, 1> &v)
    {
        return DVector2(DScalar2(v.x()), DScalar2(v.y()));
    }
 
    static inline DVector3 vector(const Eigen::Matrix<Scalar, 3, 1> &v)
    {
        return DVector3(DScalar2(v.x()), DScalar2(v.y()), DScalar2(v.z()));
    }
 
    // ======================================================================
protected:
    Scalar value;
    Gradient grad;
    Hessian hess;
};
 
template <typename Scalar, typename VecType, typename MatType>
std::ostream &operator<<(std::ostream &out, const DScalar2<Scalar, VecType, MatType> &s)
{
    out << "[" << s.getValue()
        << ", grad=" << s.getGradient().format(Eigen::IOFormat(4, 1, ", ", "; ", "", "", "[", "]"))
        << ", hess=" << s.getHessian().format(Eigen::IOFormat(4, 0, ", ", "; ", "", "", "[", "]"))
        << "]";
    return out;
}
 
// clang-format on
 
namespace std {
 
template <typename _Scalar, typename _Gradient>
class numeric_limits<DScalar1<_Scalar, _Gradient>>
{
public:
    static const bool is_specialized = std::numeric_limits<_Scalar>::is_specialized;
    static const bool is_signed = std::numeric_limits<_Scalar>::is_signed;
    static const bool is_integer = std::numeric_limits<_Scalar>::is_integer;
    static const bool is_exact = std::numeric_limits<_Scalar>::is_exact;
    static const int radix = std::numeric_limits<_Scalar>::radix;
    static const bool has_infinity = std::numeric_limits<_Scalar>::has_infinity;
    static const bool has_quiet_NaN = std::numeric_limits<_Scalar>::has_quiet_NaN;
    static const bool has_signaling_NaN = std::numeric_limits<_Scalar>::has_signaling_NaN;
    static const bool is_iec559 = std::numeric_limits<_Scalar>::is_iec559;
    static const bool is_bounded = std::numeric_limits<_Scalar>::is_bounded;
    static const bool is_modulo = std::numeric_limits<_Scalar>::is_modulo;
    static const bool traps = std::numeric_limits<_Scalar>::traps;
    static const bool tinyness_before = std::numeric_limits<_Scalar>::tinyness_before;
    static const std::float_denorm_style has_denorm = std::numeric_limits<_Scalar>::has_denorm;
    static const bool has_denorm_loss = std::numeric_limits<_Scalar>::has_denorm_loss;
    static const int min_exponent = std::numeric_limits<_Scalar>::min_exponent;
    static const int max_exponent = std::numeric_limits<_Scalar>::max_exponent;
    static const int min_exponent10 = std::numeric_limits<_Scalar>::min_exponent10;
    static const int max_exponent10 = std::numeric_limits<_Scalar>::max_exponent10;
    static const std::float_round_style round_style = std::numeric_limits<_Scalar>::round_style;
    static const int digits = std::numeric_limits<_Scalar>::digits;
    static const int digits10 = std::numeric_limits<_Scalar>::digits10;
    static const int max_digits10 = std::numeric_limits<_Scalar>::max_digits10;
 
    static constexpr _Scalar min() { return std::numeric_limits<_Scalar>::min(); }
 
    static constexpr _Scalar max() { return std::numeric_limits<_Scalar>::max(); }
 
    static constexpr _Scalar lowest() { return std::numeric_limits<_Scalar>::lowest(); }
 
    static constexpr _Scalar epsilon() { return std::numeric_limits<_Scalar>::epsilon(); }
 
    static constexpr _Scalar round_error() { return std::numeric_limits<_Scalar>::round_error(); }
 
    static constexpr _Scalar infinity() { return std::numeric_limits<_Scalar>::infinity(); }
 
    static constexpr _Scalar quiet_NaN() { return std::numeric_limits<_Scalar>::quiet_NaN(); }
 
    static constexpr _Scalar signaling_NaN()
    {
        return std::numeric_limits<_Scalar>::signaling_NaN();
    }
 
    static constexpr _Scalar denorm_min() { return std::numeric_limits<_Scalar>::denorm_min(); }
};
 
template <typename _Scalar, typename _Gradient, typename _Hessian>
class numeric_limits<DScalar2<_Scalar, _Gradient, _Hessian>>
{
public:
    static const bool is_specialized = std::numeric_limits<_Scalar>::is_specialized;
    static const bool is_signed = std::numeric_limits<_Scalar>::is_signed;
    static const bool is_integer = std::numeric_limits<_Scalar>::is_integer;
    static const bool is_exact = std::numeric_limits<_Scalar>::is_exact;
    static const int radix = std::numeric_limits<_Scalar>::radix;
    static const bool has_infinity = std::numeric_limits<_Scalar>::has_infinity;
    static const bool has_quiet_NaN = std::numeric_limits<_Scalar>::has_quiet_NaN;
    static const bool has_signaling_NaN = std::numeric_limits<_Scalar>::has_signaling_NaN;
    static const bool is_iec559 = std::numeric_limits<_Scalar>::is_iec559;
    static const bool is_bounded = std::numeric_limits<_Scalar>::is_bounded;
    static const bool is_modulo = std::numeric_limits<_Scalar>::is_modulo;
    static const bool traps = std::numeric_limits<_Scalar>::traps;
    static const bool tinyness_before = std::numeric_limits<_Scalar>::tinyness_before;
    static const std::float_denorm_style has_denorm = std::numeric_limits<_Scalar>::has_denorm;
    static const bool has_denorm_loss = std::numeric_limits<_Scalar>::has_denorm_loss;
    static const int min_exponent = std::numeric_limits<_Scalar>::min_exponent;
    static const int max_exponent = std::numeric_limits<_Scalar>::max_exponent;
    static const int min_exponent10 = std::numeric_limits<_Scalar>::min_exponent10;
    static const int max_exponent10 = std::numeric_limits<_Scalar>::max_exponent10;
    static const std::float_round_style round_style = std::numeric_limits<_Scalar>::round_style;
    static const int digits = std::numeric_limits<_Scalar>::digits;
    static const int digits10 = std::numeric_limits<_Scalar>::digits10;
    static const int max_digits10 = std::numeric_limits<_Scalar>::max_digits10;
 
    static constexpr _Scalar min() { return std::numeric_limits<_Scalar>::min(); }
 
    static constexpr _Scalar max() { return std::numeric_limits<_Scalar>::max(); }
 
    static constexpr _Scalar lowest() { return std::numeric_limits<_Scalar>::lowest(); }
 
    static constexpr _Scalar epsilon() { return std::numeric_limits<_Scalar>::epsilon(); }
 
    static constexpr _Scalar round_error() { return std::numeric_limits<_Scalar>::round_error(); }
 
    static constexpr _Scalar infinity() { return std::numeric_limits<_Scalar>::infinity(); }
 
    static constexpr _Scalar quiet_NaN() { return std::numeric_limits<_Scalar>::quiet_NaN(); }
 
    static constexpr _Scalar signaling_NaN()
    {
        return std::numeric_limits<_Scalar>::signaling_NaN();
    }
 
    static constexpr _Scalar denorm_min() { return std::numeric_limits<_Scalar>::denorm_min(); }
};
 
} // namespace std
 
#endif /* __AUTODIFF_H */