Использование производных в качестве функций в CppAD
Я пытаюсь изменить пример здесь:
# include <cppad/cppad.hpp>
namespace { // ---------------------------------------------------------
// define the template function JacobianCases<Vector> in empty namespace
template <typename Vector>
bool JacobianCases()
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
using CppAD::exp;
using CppAD::sin;
using CppAD::cos;
// domain space vector
size_t n = 2;
CPPAD_TESTVECTOR(AD<double>) X(n);
X[0] = 1.;
X[1] = 2.;
// declare independent variables and starting recording
CppAD::Independent(X);
// a calculation between the domain and range values
AD<double> Square = X[0] * X[0];
// range space vector
size_t m = 3;
CPPAD_TESTVECTOR(AD<double>) Y(m);
Y[0] = Square * exp( X[1] );
Y[1] = Square * sin( X[1] );
Y[2] = Square * cos( X[1] );
// create f: X -> Y and stop tape recording
CppAD::ADFun<double> f(X, Y);
// new value for the independent variable vector
Vector x(n);
x[0] = 2.;
x[1] = 1.;
// compute the derivative at this x
Vector jac( m * n );
jac = f.Jacobian(x);
/*
F'(x) = [ 2 * x[0] * exp(x[1]) , x[0] * x[0] * exp(x[1]) ]
[ 2 * x[0] * sin(x[1]) , x[0] * x[0] * cos(x[1]) ]
[ 2 * x[0] * cos(x[1]) , -x[0] * x[0] * sin(x[i]) ]
*/
ok &= NearEqual( 2.*x[0]*exp(x[1]), jac[0*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*sin(x[1]), jac[1*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*cos(x[1]), jac[2*n+0], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *exp(x[1]), jac[0*n+1], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *cos(x[1]), jac[1*n+1], eps99, eps99);
ok &= NearEqual(-x[0] * x[0] *sin(x[1]), jac[2*n+1], eps99, eps99);
return ok;
}
} // End empty namespace
# include <vector>
# include <valarray>
bool Jacobian(void)
{ bool ok = true;
// Run with Vector equal to three different cases
// all of which are Simple Vectors with elements of type double.
ok &= JacobianCases< CppAD::vector <double> >();
ok &= JacobianCases< std::vector <double> >();
ok &= JacobianCases< std::valarray <double> >();
return ok;
}
Я пытаюсь изменить его следующим образом:
Пусть G якобиан jac который рассчитывается в этом примере, в строке:
jac = f.Jacobian(x);
и, как в примере, пусть X быть независимыми переменными. Я хотел бы построить новую функцию, H, который является функцией jacт.е. H(jacobian(X)) = что-то такое, что H автодифференцируемо. Примером может быть H(X) = jacobian( jacobian(X)[0])Якобиан первого элемента jacobian(X) WRT X (вторая производная сортов).
Проблема в том, что jac как написано здесь, имеет тип Vector, который является параметризованным типом на сыром doubleне AD<double>, Насколько мне известно, это означает, что выходные данные не являются автоматически дифференцируемыми.
Я ищу несколько советов о том, можно ли использовать якобиан в более крупной операции и взять якобиан из этой более крупной операции (в отличие от любого арифметического оператора) или если это невозможно.
РЕДАКТИРОВАТЬ: Это было выставлено за награду один раз, но я поднимаю его снова, чтобы посмотреть, есть ли лучшее решение, потому что я думаю, что это важно. Чтобы быть немного более понятным, элементы, которые нужны для "правильного" ответа:
а) Средства для вычисления производных произвольного порядка.
б) разумный способ не указывать порядок производных априори. Если производная максимального порядка должна быть известна во время компиляции, порядок производной не может быть определен алгоритмически. Кроме того, указание чрезвычайно большого порядка в текущем ответе приведет к проблемам с выделением памяти и, как мне кажется, к проблемам с производительностью.
в) Абстрагирование шаблонов производного заказа от конечного пользователя. Это важно, потому что может быть трудно отслеживать порядок необходимых производных. Это, вероятно, то, что приходит "бесплатно", если б) решена.
Если кто-нибудь сможет это взломать, это будет потрясающий вклад и чрезвычайно полезная операция.
2 ответа
Если вы хотите вложить функции, вы должны также вложить AD<>. Вы можете вкладывать якобианы в другие функции, например, смотрите фрагмент кода ниже, который вычисляет двойную производную путем вложения якобиана
#include <cstring>
#include <iostream> // standard input/output
#include <vector> // standard vector
#include <cppad/cppad.hpp> // the CppAD package http://www.coin-or.org/CppAD/
// main program
int main(void)
{ using CppAD::AD; // use AD as abbreviation for CppAD::AD
using std::vector; // use vector as abbreviation for std::vector
size_t i; // a temporary index
// domain space vector
auto Square = [](auto t){return t*t;};
vector< AD<AD<double>> > X(1); // vector of domain space variables
// declare independent variables and start recording operation sequence
CppAD::Independent(X);
// range space vector
vector< AD<AD<double>> > Y(1); // vector of ranges space variables
Y[0] = Square(X[0]); // value during recording of operations
// store operation sequence in f: X -> Y and stop recording
CppAD::ADFun<AD<double>> f(X, Y);
// compute derivative using operation sequence stored in f
vector<AD<double>> jac(1); // Jacobian of f (m by n matrix)
vector<AD<double>> x(1); // domain space vector
CppAD::Independent(x);
jac = f.Jacobian(x); // Jacobian for operation sequence
CppAD::ADFun<double> f2(x, jac);
vector<double> result(1);
vector<double> x_res(1);
x_res[0]=15.;
result=f2.Jacobian(x_res);
// print the results
std::cout << "f'' computed by CppAD = " << result[0] << std::endl;
}
Как примечание, поскольку C++14 или 11, внедряющие шаблоны выражений и автоматическое разграничение, стали проще и могут быть выполнены с гораздо меньшими усилиями, как показано, например, в этом видео ближе к концу https://www.youtube.com/watch?v=cC9MtflQ_nI (извините за плохое качество). Если бы мне пришлось реализовывать достаточно простые символические операции, я бы начал с нуля с современным C++: вы можете писать более простой код и получать ошибки, которые вы можете легко понять.
Редактировать: обобщение примера для построения производных произвольного порядка может быть шаблоном метапрограммирования. Приведенный ниже фрагмент кода показывает, что это возможно при использовании рекурсии шаблона.
#include <cstring>
#include <iostream>
#include <vector>
#include <cppad/cppad.hpp>
using CppAD::AD;
using std::vector;
template<typename T>
struct remove_ad{
using type=T;
};
template<typename T>
struct remove_ad<AD<T>>{
using type=T;
};
template<int N>
struct derivative{
using type = AD<typename derivative<N-1>::type >;
static constexpr int order = N;
};
template<>
struct derivative<0>{
using type = double;
static constexpr int order = 0;
};
template<typename T>
struct Jac{
using value_type = typename remove_ad<typename T::type>::type;
template<typename P, typename Q>
auto operator()(P & X, Q & Y){
CppAD::ADFun<value_type> f(X, Y);
vector<value_type> jac(1);
vector<value_type> x(1);
CppAD::Independent(x);
jac = f.Jacobian(x);
return Jac<derivative<T::order-1>>{}(x, jac);
}
};
template<>
struct Jac<derivative<1>>{
using value_type = derivative<0>::type;
template<typename P, typename Q>
auto operator()(P & x, Q & jac){
CppAD::ADFun<value_type> f2(x, jac);
vector<value_type> res(1);
vector<value_type> x_res(1);
x_res[0]=15.;
return f2.Jacobian(x_res);
}
};
int main(void)
{
constexpr int order=4;
auto Square = [](auto t){return t*t;};
vector< typename derivative<order>::type > X(1);
vector< typename derivative<order>::type > Y(1);
CppAD::Independent(X);
Y[0] = Square(X[0]);
auto result = Jac<derivative<order>>{}(X, Y);
std::cout << "f'' computed by CppAD = " << result[0] << std::endl;
}
В CppAD появилась новая функция, которая устраняет необходимость в AD