check_derivatives.h

Go to the documentation of this file.

/* This file is part of the Tomographer project, which is distributed under the
 * terms of the MIT license.
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2015 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef CHECK_DERIVATIVES_H
#define CHECK_DERIVATIVES_H

#include <cstddef>
#include <iostream>

#include <Eigen/Core>

namespace Tomographer
{
namespace Tools
{

template<typename Der1, typename Der2, typename fnType, typename ErrorStream>
bool check_derivatives(const Eigen::ArrayBase<Der1> & derivatives, const Eigen::MatrixBase<Der2> & point,
                       fnType fn, std::size_t valdims,
                       typename Eigen::MatrixBase<Der2>::Scalar delta = 1e-6,
                       typename Eigen::MatrixBase<Der1>::Scalar tol = 1e-6,
                       ErrorStream & error_stream = std::cerr
                       )
{
  bool ok = true;

  const std::size_t xdims = derivatives.cols();
  eigen_assert(point.rows() == (int)xdims);
  eigen_assert(point.cols() == (int)1);
  eigen_assert(derivatives.rows() == (int)valdims);

  Eigen::VectorXd val0(valdims);
  Eigen::VectorXd dval1(valdims);
  Eigen::VectorXd dvalFromDer(valdims);

  Eigen::VectorXd pt2(point.rows());

  // calculate the base point
  fn(val0, point);

  std::size_t i;
  for (i = 0; i < xdims; ++i) {
    // numerically calculate the differences ...

    pt2 = point;
    pt2(i) += delta;

    fn(dval1, pt2);
    dval1 -= val0;

    dvalFromDer = delta * derivatives.matrix().col(i);//derivatives.matrix() * (delta*dir);

    double thediff = (dval1 - dvalFromDer).norm();

    if (thediff/delta > tol ) {
      // Error in the derivative

      // direction in which we probed, for the error_stream
      Eigen::VectorXd dir = Eigen::VectorXd::Zero(xdims);
      dir(i) = 1.0;

      ok = false;
      error_stream
        << "Error in derivative check: Derivative wrong in direction\n"
        << "dir = " << dir.transpose() << "   [basis vector #"<<i<<"]\n"
        << "\tpoint = \t" << point.transpose() << "\n"
        << "\tval0  = \t" << val0.transpose() << "\n"
        << "\tdval1 = \t" << dval1.transpose() << "\n"
        << "\tdvalFromDer = \t"<<dvalFromDer.transpose() << "\n"
        << "\tderivative in this direction =\n\t\t\t\t" << derivatives.transpose().block(i,0,1,valdims) << "\n"
        << "--> difference in p2-points: \t" << thediff << "\n"
        << "--> difference in derivatives: \t" << thediff/delta << "\n\n";
    }
  }

  return ok;
}



} // Tools
} // Tomographer



#endif