Использование производных как функций в 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) по сравнению с X (своего рода вторая производная).
Проблема в том, что jac, как здесь написано, относится к типу Vector, который является параметризованным типом на необработанном double, а не на AD<double>. Насколько мне известно, это означает, что выходные данные не поддаются автоматической дифференциации.
Мне нужен совет о том, можно ли использовать якобиан в более крупной операции и взять якобиан этой более крупной операции (в отличие от любого арифметического оператора), или если это невозможно.
Обновлено: Это однажды было объявлено за награду, но я снова поднимаю его, чтобы увидеть, есть ли лучшее решение, потому что я думаю, что это важно. Чтобы быть немного более ясным, элементы, которые необходимы для "правильного" ответа:
а) Средство вычисления производных произвольного порядка.
б) Интеллектуальный способ не указывать порядок производных априори. Если производная максимального порядка должна быть известна во время компиляции, порядок производной не может быть определен алгоритмически. Кроме того, указание чрезвычайно большого порядка, как в данном ответе, приведет к проблемам с распределением памяти и, как я полагаю, к проблемам с производительностью.
c) Абстрагирование шаблонов производного ордера от конечного пользователя. Это важно, потому что может быть сложно отследить порядок необходимых производных. Вероятно, это то, что приходит «бесплатно», если решено б).
Если кто-нибудь сможет это взломать, это будет отличный вклад и чрезвычайно полезная операция.
Пока что похоже, что несколько уровней могут быть правильным ответом (например, в coin-or.github.io/CppAD/doc/mul_level.cpp.htm), но мне не ясно, как применять это произвольное количество раз, что мне нужно.
Ответы 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://thewikihow.com/video_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;
}
Спасибо за пример. Это почти работает - за исключением того, что вам нужно знать количество уровней дифференциации, которые вы хотите априори. Есть ли способ взять производные n-го порядка без предварительного указания n?
Я добавил еще один фрагмент, в котором вычисляется произвольный порядок вывода
Спасибо. Я еще не проверял это (надеюсь, завтра), но я все равно назначил вам награду, так как я не смогу проверить ее до тех пор, пока она не закончится.
Не могли бы вы более подробно объяснить, что делает каждый из этих компонентов?
Я считаю, что в этой инструкции действительно есть ошибки, так как для order = 2 она выводит "f", вычисленное с помощью CppAD = 0 ", что неверно; должно быть 2.0.
да ладно, была ошибка. Я называл Independent (X) и Y [0] = f (X) в неправильном порядке. Я починил это. Я добавлю немного дополнительных пояснений ...
Спасибо за быстрое исправление. Еще один вопрос - неясно, как можно использовать этот код для передачи значения для зависимого вектора. Как бы это сделать?
В CppAD появилась новая функция, которая устраняет необходимость в AD <AD>, см. https://coin-or.github.io/CppAD/doc/base2ad.cpp.htm
Несмотря на то, что это технически отвечает на вопрос, было бы намного лучше, если бы он включал какой-либо соответствующий код в сам ответ.
Я добавил запрошенный код. Извини за это!