check_grad_ode.hpp

View page source

C++ Check Gradient of an_ode

Syntax

      # include "check_grad_ode.hpp"
      ok = check_grad_ode( grad_ode )

Prototype

template <class Gradient>
bool check_grad_ode( Gradient& grad_ode )

grad_ode

Is a cpp_fun_obj object that computes the gradient for the ref:cpp_an_ode-name algorithm.

ok

is true (false) if the gradient passes (fails) the test.

Gradient

We use \(r\) to denote the range space component of the an_ode algorithm that we are computing the gradient for.

\[y_r (t) = \frac{ t^{r+1} }{ (r+1)! } \prod{j=0}^r x_j\]
\[\begin{split}\frac{ \partial y_r (t) }{ \partial x_j } = \begin{cases} y_r (t) / x_j & \text{if} \; j \leq r \\ 0 & \text{otherwise} \end{cases}\end{split}\]

Source Code

# include <cmpad/uniform_01.hpp>
# include <cmpad/near_equal.hpp>

// BEGIN PROTOTYPE
template <class Gradient>
bool check_grad_ode( Gradient& grad_ode )
// END PROTOTYPE
{  //
   // ok
   bool ok = true;
   //
   // rel_error
   double rel_error = 500. * std::numeric_limits<double>::epsilon();
   //
   // near_equal
   using cmpad::near_equal;
   //
   // n_arg
   size_t n_arg = 4;
   //
   // time_setup
   for(bool time_setup : { true, false } )
   {  //
      // option
      cmpad::option_t option;
      option.n_arg      = n_arg;
      option.n_other    = 10;
      option.time_setup = time_setup;
      //
      // grad_ode
      grad_ode.setup(option);
      //
      // x
      // note that x[i] != 0.0 so can divide by it
      cmpad::vector<double> x(n_arg);
      cmpad::uniform_01(x);
      for(size_t i = 0; i < n_arg; ++i)
         x[i] += 1.0;
      //
      // g
      cmpad::vector<double> g = grad_ode(x);
      //
      // r
      size_t r = n_arg - 1;
      //
      // y_r
      double tf  = 2.0;
      double y_r = x[0] * tf;
      for(size_t j = 1; j <= r; ++j)
      {  y_r = y_r * x[j] * tf / double(j+1);
      }
      //
      // ok
      for(size_t j = 0; j <= r; ++j)
         ok &= cmpad::near_equal( g[j], y_r / x[j], rel_error );
      for(size_t j = r+1; j < n_arg; ++j)
         ok &= g[j] == 0.0;
   }
   return ok;
}