Automatic differentiation using C++ exploiting meta programming (part 2)

A readable & functioning version of the code can be found at github.

As was pointed out, the naive application of the class shown in part 1 of this series would cause unnecessary operations, using statically zero derivative values in order to create other statically zero derivatives, which in turn...

It would constitute a gigantic waste of CPU power.

What is needed is a method to attach information to an instance of the DualNumber class about which independent variables it is carrying. This attaching must necessarily be done via template meta programming, because we want that anything that can be decided at compile time, will be decided at compile time and avoid burdening the CPU at runtime with unnecessary code, which would anyway yield the same result at every execution.

1. Independent variables are identified by enumeration values.
2. The class DualNumber must be a template class, and the template argument must contain the information about which derivatives -- wrt which independent variables -- are contained.
3. To avoid unnecessarily complicating the compiler's job there must be only a single type of this class for a given set of independent variables -- and also because we want to exploit performance advantages of operations between identical types. This can be realized by that the passed set of independent variables must be sorted.
4. The template argument must be used to determine the size of the contained array of derivatives.
5. Operations (e.g., multiplication) between two different DualNumber types A & B , must return another DualNumber type C, which necessarily carries more derivatives than either of the argument types (sizeof(C) > max(sizeof(A), sizeof(B))).

This functionality can be realized by using boost::mp11. E.g. a possible template argument type identifying a DualNumber type containing derivatives wrt independent variables 1 and 3:

typedef boost::mp11::mp_list<
        boost::mp11::mp_size_t<1>,
        boost::mp11::mp_size_t<3>
> oneAndThree;        

In order to extract the size of this type, the boost::mp11::mp_size<> meta function can be used:

static constexpr std::size_t SIZE = boost::mp11::mp_size<oneAndThree>::value;
typedef std::array<double, SIZE> ARRAY;        

In order to realize the 5th point above, we need a meta function to merge two such argument types like oneAndThree. This can be realized by boost::mp11::mp_set_union<>. But this meta function alone would not result in a sorted set. So in addition boost::mp11::mp_sort<> must be used:

template<typename L, typename R>
struct merge
{       typedef boost::mp11::mp_set_union<L, R> merged;
        typedef boost::mp11::mp_sort<merged, boost::mp11::mp_less> type;
};        

With this we have everything needed to create the declaration for any operator using two different DualNumber types as input -- using the multiplication operator as an example and by exploiting above introduced meta function called merge:

/// within template<typename SET> struct DualNumber;
template<typename SET1>
typename merge<SET, SET1>::type operator*(
        const DualNumber<SET1>&_r
) const;        

Now to the implementation of this multiplication operator -- there are two ways. One via generated constexpr tables and one via boost::mp11::mp_for_each<SET>(functor). I would like to demonstrate the latter:

/// within template<typename SET> struct DualNumber;
template<typename SET1>
struct multiplication
{	typedef typename merge<SET, SET1>::type TYPE;
	DualNumber<TYPE>&m_rResult;
	const DualNumber<SET>&m_rL;
	const DualNumber<SET1>&m_rR;
	multiplication(
		DualNumber<TYPE>&_rResult,
		const DualNumber<SET>&_rL,
		const DualNumber<SET1>&_rR
	)
		:m_rResult(_rResult),
		m_rL(_rL),
		m_rR(_rR)
	{
	}
		/// derivative is contained in both arguments
	template<typename TPOS, typename LPOS, typename RPOS>
	void operator()(const TPOS&, const LPOS&, const RPOS&) const
	{	m_rResult.m_adDer[TPOS::value] =
			m_rL.m_adDer[LPOS::value]*m_rR.m_dValue
				+ m_rR.m_adDer[RPOS::value]*m_rL.m_dValue;
	}
		/// derivative is not contained in the left argument
	template<typename TPOS, typename RPOS>
	void operator()(
		const TPOS&,
		const boost::mp11::mp_size_t<DualNumber<SET>::SIZE>&,
		const RPOS&
	) const
	{	m_rResult.m_adDer[TPOS::value] =
			m_rR.m_adDer[RPOS::value]*m_rL.m_dValue;
	}
		/// derivative is not contained in the right argument
	template<typename TPOS, typename LPOS>
	void operator()(
		const TPOS&,
		const LPOS&,
		const boost::mp11::mp_size_t<DualNumber<SET1>::SIZE>&
	) const
	{	m_rResult.m_adDer[TPOS::value] =
			m_rL.m_adDer[LPOS::value]*m_rR.m_dValue;
	}
		/// entry point called for a particular independent variable in the TARGET
	template<typename INDEP>
	void operator()(const INDEP&) const
	{	(*this)(
			boost::mp11::mp_find<TYPE, INDEP>(),
			boost::mp11::mp_find<SET, INDEP>(),
			boost::mp11::mp_find<SET1, INDEP>()
		);
	}
};
	/// multiplication operator for different DualNumber types
template<typename SET1>
auto operator*(
	const DualNumber<SET1>&_r
) const
{	typedef typename merge<SET, SET1>::type TYPE;
	DualNumber<TYPE> s;
	s.m_dValue = m_dValue*_r.m_dValue;
	boost::mp11::mp_for_each<TYPE>(multiplication<SET1>(s, *this, _r));
	return s;
}        

The above code is using boost::mp11::mp_find<> to translate an enumeration value representing an independent variable into a position (or boost::mp11::mp_size<> in case of not found).

A working example can be found at github.

A fully developed implementation of this idea can be found at https://github.com/ExcessPhase/ctaylor/blob/master/cjacobian.h.