\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

cppad_det_minor.cpp

View page source

Cppad Speed: Gradient of Determinant by Minor Expansion

Specifications

See link_det_minor .

Implementation

# include <cppad/cppad.hpp>
# include <cppad/speed/det_by_minor.hpp>
# include <cppad/speed/uniform_01.hpp>

// Note that CppAD uses global_option["memory"] at the main program level
# include <map>
extern std::map<std::string, bool> global_option;
// see comments in main program for this external
extern size_t global_cppad_thread_alloc_inuse;

namespace {
   // typedefs
   typedef CppAD::AD<double>       a_double;
   typedef CppAD::vector<a_double> a_vector;
   //
   // setup
   void setup(
      // inputs
      size_t size             ,
      // outputs
      CppAD::ADFun<double>& f )
   {
      // object for computing determinant
      CppAD::det_by_minor<a_double>   a_det(size);
      //
      // number of independent variables
      size_t nx = size * size;
      //
      // choose a matrix
      CppAD::vector<double> matrix(nx);
      CppAD::uniform_01(nx, matrix);
      //
      // copy to independent variables
      a_vector   a_A(nx);
      for(size_t j = 0; j < nx; ++j)
         a_A[j] = matrix[j];
      //
      // declare independent variables for function computation
      bool record_compare   = false;
      size_t abort_op_index = 0;
      CppAD::Independent(a_A, abort_op_index, record_compare);
      //
      // AD computation of the determinant
      a_vector a_detA(1);
      a_detA[0] = a_det(a_A);
      //
      // f : A -> detA
      f.Dependent(a_A, a_detA);
      //
      // optimize
      std::string optimize_options =
      "no_conditional_skip no_compare_op no_print_for_op no_cumulative_sum_op";
      if( global_option["val_graph"] )
         optimize_options += " val_graph";
      if( global_option["optimize"] )
         f.optimize(optimize_options);
   }

}

bool link_det_minor(
   const std::string&         job      ,
   size_t                     size     ,
   size_t                     repeat   ,
   CppAD::vector<double>     &matrix   ,
   CppAD::vector<double>     &gradient )
{  global_cppad_thread_alloc_inuse = 0;

   // --------------------------------------------------------------------
   // check global options
   const char* valid[] = { "memory", "onetape", "optimize", "val_graph"};
   size_t n_valid = sizeof(valid) / sizeof(valid[0]);
   typedef std::map<std::string, bool>::iterator iterator;
   //
   for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
   {  if( itr->second )
      {  bool ok = false;
         for(size_t i = 0; i < n_valid; i++)
            ok |= itr->first == valid[i];
         if( ! ok )
            return false;
      }
   }
   // ---------------------------------------------------------------------
   //
   // AD function mapping matrix to determinant
   static CppAD::ADFun<double> static_f;
   //
   // size corresponding to static_f
   static size_t static_size = 0;
   //
   // number of independent variables
   size_t nx = size * size;
   //
   // vectors of reverse mode weights
   CppAD::vector<double> w(1);
   w[0] = 1.;
   //
   // onetape
   bool onetape = global_option["onetape"];
   // -----------------------------------------------------------------------
   if( job == "setup" )
   {  if( onetape )
      {  setup(size, static_f);
         static_size = size;
      }
      else
      {  static_size = 0;
      }
      return true;
   }
   if( job ==  "teardown" )
   {  static_f = CppAD::ADFun<double>();
      return true;
   }
   // -----------------------------------------------------------------------
   CPPAD_ASSERT_UNKNOWN( job == "run" );
   while(repeat--)
   {  if( onetape )
      {  // use if before assert to avoid warning
         if( size != static_size )
         {  CPPAD_ASSERT_UNKNOWN( size == static_size );
         }
      }
      else
      {  setup(size, static_f);
      }
      // get next matrix
      CppAD::uniform_01(nx, matrix);

      // evaluate the gradient
      static_f.Forward(0, matrix);
      gradient = static_f.Reverse(1, w);
   }
   size_t thread                   = CppAD::thread_alloc::thread_num();
   global_cppad_thread_alloc_inuse = CppAD::thread_alloc::inuse(thread);
   return true;
}