Goal: Learn how to add your own aggregated types to CoDiPack.
Prerequisite: Tutorial 2 - Reverse mode AD
Function: Norm of two scaled vector.
template<typename VectorType>
ActiveNumber func(ActiveNumber const& s1, VectorType const& v1, ActiveNumber const& s2, VectorType const& v2) {
return (s1 * v1 + s2 * v2).norm();
}
Full code:
using Number = double;
using Vector4 = Eigen::Matrix<Number, 4, 1>;
using Vector4Transpose = Eigen::Matrix<Number, 1, 4>;
template<>
}
template<typename T>
return s * v;
}
return v;
}
return s;
}
};
template<typename ArgS, typename ArgV>
}
template<typename T>
return v1 + v2;
}
return 1.0;
}
return 1.0;
}
};
template<typename ArgV1, typename ArgV2>
}
template<typename T>
return v.norm();
}
return v / result;
}
};
template<typename Impl>
public:
return codi::ComputeExpression<Number, VectorNormOperation, Impl>(this->
cast());
}
private:
return static_cast<Impl const&
>(*this);
}
};
using Base = codi::AggregatedActiveType<Vector4, ActiveNumber, ActiveVector4>;
template<typename Arg1, typename Arg2, typename Arg3>
codi::ExpressionInterface<Number, Arg1> const& arg1,
codi::ExpressionInterface<Number, Arg2> const& arg2,
codi::ExpressionInterface<Number, Arg3> const& arg3,
codi::ExpressionInterface<Number, Arg3>
const& arg4) :
Base() {
}
}
};
template<typename VectorType>
ActiveNumber func(ActiveNumber const& s1, VectorType const& v1, ActiveNumber const& s2, VectorType const& v2) {
return (s1 * v1 + s2 * v2).norm();
}
template<typename VectorType>
void test() {
VectorType v1(ActiveNumber(1.0), ActiveNumber(2.0), ActiveNumber(4.0), ActiveNumber(8.0));
VectorType v2(ActiveNumber(0.1), ActiveNumber(0.2), ActiveNumber(0.4), ActiveNumber(0.8));
ActiveNumber s1(5.0);
ActiveNumber s2(0.5);
for(int i = 0; i < 4; i += 1) {
}
ActiveNumber res = func(s1, v1, s2, v2);
std::cout <<
"d f/d s1 = " << s1.
getGradient() << std::endl;
std::cout <<
"d f/d s2 = " << s2.
getGradient() << std::endl;
for(int i = 0; i < 4; i += 1) {
std::cout << "d f/d v1[" << i << "] = " << v1[i].getGradient() << std::endl;
}
for(int i = 0; i < 4; i += 1) {
std::cout << "d f/d v2[" << i << "] = " << v2[i].getGradient() << std::endl;
}
}
int main() {
using Vector4WithActiveType = Eigen::Matrix<ActiveNumber, 4, 1>;
std::cout << "Running example with 'Vector4WithActiveType' vector type. No specialization are used for Eigen vector." << std::endl;
test<Vector4WithActiveType>();
std::cout << std::endl;
std::cout << "Running example with 'ActiveVector4' vector type. The specializations are used for Eigen vector." << std::endl;
test<ActiveVector4>();
return 0;
}
Introduction
An aggregated type, in the sense of CoDiPack, is a structure that can be expressed by a set of double values. E.g. std::complex<double> can be represented by two double values. If aggregated types are added to the CoDiPack expression tree, then memory and runtime of differentiated programs can be reduced.
Defining the basic types and the active vector type
In this tutorial we want to add specialized expressions for the Eigen Vector4 type. Therefore, we decalare first the basic types as:
using Number = double;
using Vector4 = Eigen::Matrix<Number, 4, 1>;
using Vector4Transpose = Eigen::Matrix<Number, 1, 4>;
It is now possible to declare the vector type that integrates into the CoDiPack expressions. The simplest way is to extend from codi::AggregatedActiveType. This base class implements all the necessary interfaces and functions for the compatibility with the CoDiPack expressions. In order to do this, codi::AggregatedActiveType needs some information about the used type, which is in our case Vector4. This is done by the specialization of codi::RealTraits::AggregatedTypeTraits. Here, the array constructor, the array access methods, and the number of elements for the type are defined. If the type can be described as an array of floating point values, then the default implementation codi::RealTraits::ArrayAggregatedTypeTraitsBase can be used. In addition we need to specify a transposed operation for the type. All together this is done with:
using Base = codi::AggregatedActiveType<Vector4, ActiveNumber, ActiveVector4>;
template<typename Arg1, typename Arg2, typename Arg3>
codi::ExpressionInterface<Number, Arg1> const& arg1,
codi::ExpressionInterface<Number, Arg2> const& arg2,
codi::ExpressionInterface<Number, Arg3> const& arg3,
codi::ExpressionInterface<Number, Arg3>
const& arg4) :
Base() {
}
}
};
Currently we only define a constructor and the array access for the ActiveVector4.
Definition of operations
For our function example, we need two operations the scalar multiplication and the vector addition. Both can be defined by implementing the codi::BinaryOperation interface. This interface defines the primal operation and the gradient computations. For the real valued case, the operation implementations can always assume that the Jacobian has the same type as the primal result. For the vector operations, this is no longer the case. Therfore, we also need to declare the Jacobian types in the operation implementations:
template<typename T>
return s * v;
}
return v;
}
return s;
}
};
template<typename ArgS, typename ArgV>
}
template<typename T>
return v1 + v2;
}
return 1.0;
}
return 1.0;
}
};
template<typename ArgV1, typename ArgV2>
}
Definition of member operations
If the base type contains member operations, then these need to be specified to. These operations could be specified in the implementation of ActiveVector4, but then they would not be availalbe on expressions. For example, it would not be possible to write (v1 + v2).norm(), since (v1 + v2) has the type codi::ComputeExpression. It is therefore necessary to inject the member operations into all expressions. This can be done by a specialization of codi::ExpressionMemberOperations for Vector4:
template<typename T>
return v.norm();
}
return v / result;
}
};
template<typename Impl>
public:
return codi::ComputeExpression<Number, VectorNormOperation, Impl>(this->
cast());
}
private:
return static_cast<Impl const&
>(*this);
}
};
Advantages and further hints
With the finished implementation we can now call the function defined above. We will do this once with an unspecialized type (using Vector4WithActiveType = Eigen::Matrix<ActiveNumber, 4, 1>) and once with the ActiveVector4. In the test we measure the tape memory that is recorded during the function evaluation. The result is:
Running example with 'Vector4WithActiveType' vector type. No specialization are used for Eigen vector.
-------------------------------------
CoDi Tape Statistics ( JacobianLinearTape )
-------------------------------------
Total memory used : 648.00 B
-------------------------------------
Statement entries
-------------------------------------
Total number : 20
-------------------------------------
Jacobian entries
-------------------------------------
Total number : 39
-------------------------------------
Running example with 'ActiveVector4' vector type. The specializations are used for Eigen vector.
-------------------------------------
CoDi Tape Statistics ( JacobianLinearTape )
-------------------------------------
Total memory used : 129.00 B
-------------------------------------
-------------------------------------
Statement entries
-------------------------------------
Total number : 1
-------------------------------------
Jacobian entries
-------------------------------------
Total number : 10
-------------------------------------
With the specialization only 129 bytes are used which is a reduction of about 80%. Since we integrated the ActiveVector4 into the CoDiPack expressions, the tape no longer records all intermediate steps in the function evaluation. The statement in the function is seen by CoDiPack as one big expression which reduces the number of statment and Jacobian entries.
Currently, all uses of Vector4 in an application would need to be exchanged with ActiveVector4. A declaration like Eigen::Matrix<ActiveNumber, 4, 1> would not automatically use the specialization. An automatic use can be achieved by a direct specialization of Eigen::Matrix for a CoDiPack type. E.g:
template<>
struct Eigen::Matrix<ActiveNumber, 4, 1, Eigen::AutoAlign, 4, 1> :
...
};
Represents a concrete aggregated lvalue int the CoDiPack expression tree.
Definition aggregatedActiveType.hpp:164
With this specialization both runs of the test function have the same reduced memory consumption.
Notes on implementation
The current implementation of aggregated types is tailored to the implementation of std::complex. The current expression framework of CoDiPack has problems in handling operations where the operators have not the commutativity property. It is also problematic to handle operations where the AD reverse mode can not be expressed as a Jacobian times a bar value. These cases can be handled, but their implementation is quite cumbersome. We will extend the expression framework in the near future such that these cases can be handled in an improved way.
1.13.2