Я пытаюсь изменить пример здесь :
# 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) относительно X (вторая производная от рода).
Проблема в том, что jac, как написано здесь, имеет типVector, который является параметризованным типом для необработанного double, а не AD<double>.Насколько мне известно, это означает, что выходные данные не являются авто-дифференцируемыми.
Я ищу несколько советов о том, можно ли использовать якобиан в более крупной операции и взять якобиан из этой более крупной операции (не такой как в любом другом случае).арифметический оператор) или, если это невозможно.
РЕДАКТИРОВАТЬ: Это было выставлено за награду один раз, но я поднимаю его снова, чтобы увидеть, есть ли лучшее решение, потому что я думаю, что это важно,Чтобы быть немного более понятным, элементы, в которых нуждается «правильный» ответ:
a) Средство вычисления производных произвольного порядка.
b) Разумный способ не указыватьпорядок производных априори.Если производная максимального порядка должна быть известна во время компиляции, порядок производной не может быть определен алгоритмически.Кроме того, указание чрезвычайно большого порядка, как в текущем ответе, приведет к проблемам с выделением памяти и, я думаю, к проблемам с производительностью.
c) Абстрагирование шаблонов производного порядка от конечного пользователя.Это важно, потому что может быть трудно отслеживать порядок необходимых производных.Это, вероятно, то, что приходит «бесплатно», если б) решено.
Если кто-нибудь сможет это взломать, это будет потрясающий вклад и чрезвычайно полезная операция.