Utilisation de dérivés comme fonctions dans CppAD

Question

Je suis en train de modifier l'exemple ici:

# 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;
}

Je suis en train de le modifier de la façon suivante:

Soit G le Jacobien jac , qui est calculé dans cet exemple, dans la ligne:

jac = f.Jacobian(x);

et, comme dans l'exemple, laissez - X être les variables indépendantes. Je voudrais construire une nouvelle fonction, H, qui est une fonction de l' jac, c'est à dire H(jacobian(X)) = quelque chose, de telle sorte que H est autodifferentiable. Un exemple peut être H(X) = jacobian( jacobian(X)[0]), c'est à dire le jacobien du premier élément de jacobian(X) w.r.t X (la dérivée seconde de toutes sortes).

Le problème est qu' jac comme écrit ici est de type Vector, ce qui est un type paramétré sur un raw double, pas un AD<double>. À ma connaissance, cela signifie que la sortie n'est pas autodifferentiable.

Je suis à la recherche de quelques conseils sur si il est possible d'utiliser le Jacobien dans une grande opération, et de prendre le Jacobien de cette grande opération (qui ne ressemble à aucun opérateur arithmétique) ou si cela n'est pas possible.

EDIT: Cela a été mis en place pour une prime à la fois, mais je suis en train de mettre de nouveau pour voir si il y a une meilleure solution, car je pense que c'est important. Pour être un peu plus clair, les éléments qui la "bonne" réponse de besoins:

a) Un moyen de calcul de l'ordre arbitraire de produits dérivés.

b) Une façon intelligente de ne pas avoir à spécifier l'ordre des dérivés de l'a priori. Si le maximum de l'ordre de dérivation doit être connue au moment de la compilation, de l'ordre de la dérivée ne peut pas être déterminé de manière algorithmique. En outre, la spécification d'un très grand ordre que dans le courant de la réponse donnée conduira à la mémoire des questions de répartition et, j'imagine, des problèmes de performances.

c) en faisant Abstraction de la création de modèles de dérivés de l'ordre l'écart de l'utilisateur final. Ceci est important, car il peut être difficile de garder une trace de l'ordre de dérivés de besoin. C'est probablement quelque chose qui vient "gratuitement" si b) est résolu.

Si quelqu'un peut se fissurer, ce serait une formidable contribution et un outil très utile de l'opération.

Answer 1

Si vous souhaitez imbriquer des fonctions, vous devez imbriquer les AD<> ainsi. Vous pouvez imbriquer des Jacobiens que d'autres fonctions, voir par exemple l'extrait de code ci-dessous, qui est le calcul de la double dérivés par l'imbrication Jacobien

#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;
}

Comme une note, depuis le C++14 ou 11 la mise en œuvre de l'expression des modèles et de la différentiation automatique est devenu plus facile et peut être fait avec beaucoup moins d'effort, comme le montre par exemple dans cette vidéo vers la fin https://www.youtube.com/watch?v=cC9MtflQ_nI (désolé pour la mauvaise qualité). Si j'avais à mettre en œuvre assez simple symbolique des opérations, je voudrais commencer à partir de zéro avec C++ moderne: vous pouvez écrire du code plus simple, et vous obtenez des erreurs que vous pouvez comprendre facilement.

Edit: Généraliser l'exemple de construire un ordre arbitraire des produits peut être un modèle de métaprogrammation exercice. L'extrait de code ci-dessous montre qu'il est possible à l'aide du modèle de la récursivité

#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;
} 

Answer 2

Il y a une nouvelle fonctionnalité dans CppAD qui élimine le besoin d'AD <AD>, voir https://coin-or.github.io/CppAD/doc/base2ad.cpp.htm