functionsanitychecks.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include "findlinesegment.hh"
 
#include <iostream>
 
#include <dune/common/float_cmp.hh>
 
#include <spdlog/spdlog.h>
namespace Ikarus::utils {
  namespace Impl {
    double drawResultAndReturnSlope(std::string&& functionName, const std::function<double(double)>& ftfunc, bool draw,
                                    int slopeOfReference);
  }  // namespace Impl
  struct CheckFlags {
    bool draw                        = true;
    bool writeSlopeStatementIfFailed = true;
    double tolerance                 = 1e-2;
  };
 
  template <typename NonlinearOperator, typename UpdateType = typename NonlinearOperator::template ParameterValue<0>>
  bool checkGradient(
      NonlinearOperator& nonLinOp, CheckFlags checkFlags = CheckFlags(),
      std::function<void(typename NonlinearOperator::template ParameterValue<0>&, const UpdateType&)> p_updateFunction
      = [](typename NonlinearOperator::template ParameterValue<0>& a, const UpdateType& b) { a += b; }) {
    auto& x         = nonLinOp.firstParameter();
    const auto xOld = x;
    UpdateType b;
    if constexpr (not std::is_floating_point_v<UpdateType>) {
      b.resizeLike(nonLinOp.derivative());
      b.setRandom();
      b /= b.norm();
    } else
      b = 1;
 
    nonLinOp.updateAll();
    const auto e = nonLinOp.value();
 
    double gradfv;
    if constexpr (not std::is_floating_point_v<UpdateType>)
      gradfv = nonLinOp.derivative().dot(b);
    else
      gradfv = nonLinOp.derivative() * b;
 
    auto ftfunc = [&](auto t) {
      p_updateFunction(x, t * b);
      nonLinOp.template update<0>();
      auto value = std::abs(nonLinOp.value() - e - t * gradfv);
      x          = xOld;
      return value;
    };
 
    const double slope = Impl::drawResultAndReturnSlope("Gradient", ftfunc, checkFlags.draw, 2);
 
    const bool checkPassed = Dune::FloatCmp::le(2.0, slope, checkFlags.tolerance);
 
    if (checkFlags.writeSlopeStatementIfFailed and not checkPassed) {
      spdlog::info("Gradient check:");
      spdlog::info("The slope should be 2. It seems to be {}.", slope);
      if (checkPassed)
        spdlog::info("We consider this as sufficient.");
      else
        spdlog::info("The gradient seems wrong.");
    }
 
    nonLinOp.updateAll();
    return checkPassed;
  }
 
  template <typename NonlinearOperator, typename UpdateType = typename NonlinearOperator::template ParameterValue<0>>
  bool checkJacobian(
      NonlinearOperator& nonLinOp, CheckFlags checkFlags = CheckFlags(),
      std::function<void(typename NonlinearOperator::template ParameterValue<0>&, const UpdateType&)> p_updateFunction
      = [](typename NonlinearOperator::template ParameterValue<0>& a, const UpdateType& b) { a += b; }) {
    auto& x         = nonLinOp.firstParameter();
    const auto xOld = x;
    UpdateType b;
    b.resizeLike(nonLinOp.derivative().row(0).transpose());
    b.setRandom();
    b /= b.norm();
 
    nonLinOp.updateAll();
    const auto e = nonLinOp.value();
 
    const auto jacofv = (nonLinOp.derivative() * b).eval();
 
    auto ftfunc = [&](auto t) {
      p_updateFunction(x, t * b);
      nonLinOp.template update<0>();
      auto value = (nonLinOp.value() - e - t * jacofv).norm();
      x          = xOld;
      return value;
    };
 
    const double slope = Impl::drawResultAndReturnSlope("Jacobian", ftfunc, checkFlags.draw, 2);
 
    const bool checkPassed = Dune::FloatCmp::le(2.0, slope, checkFlags.tolerance);
 
    if (checkFlags.writeSlopeStatementIfFailed and not checkPassed) {
      spdlog::info("Jacobian check:");
      spdlog::info("The slope should be 2. It seems to be {}.", slope);
      if (checkPassed)
        spdlog::info("We consider this as sufficient.");
      else
        spdlog::info("The Jacobian seems wrong.");
    }
    nonLinOp.updateAll();
    return checkPassed;
  }
 
  template <typename NonlinearOperator, typename UpdateType = typename NonlinearOperator::template ParameterValue<0>>
  bool checkHessian(
      NonlinearOperator& nonLinOp, CheckFlags checkFlags = CheckFlags(),
      std::function<void(typename NonlinearOperator::template ParameterValue<0>&, const UpdateType&)> p_updateFunction
      = [](typename NonlinearOperator::template ParameterValue<0>& a, const UpdateType& b) { a += b; }) {
    auto& x         = nonLinOp.firstParameter();
    const auto xOld = x;
    UpdateType b;
    if constexpr (not std::is_floating_point_v<UpdateType>) {
      b.resizeLike(nonLinOp.derivative());
      b.setRandom();
      b /= b.norm();
    } else
      b = 1;
 
    nonLinOp.updateAll();
    const auto e = nonLinOp.value();
 
    double gradfv, vhessv;
    if constexpr (not std::is_floating_point_v<UpdateType>) {
      gradfv = nonLinOp.derivative().dot(b);
      vhessv = (nonLinOp.secondDerivative() * b).dot(b);
    } else {
      gradfv = nonLinOp.derivative() * b;
      vhessv = nonLinOp.secondDerivative() * b * b;
    }
 
    auto ftfunc = [&](auto t) {
      p_updateFunction(x, t * b);
      nonLinOp.template update<0>();
      auto value = std::abs(nonLinOp.value() - e - t * gradfv - 0.5 * t * t * vhessv);
      x          = xOld;
      return value;
    };
 
    const double slope = Impl::drawResultAndReturnSlope("Hessian", ftfunc, checkFlags.draw, 3);
 
    const bool checkPassed = Dune::FloatCmp::le(3.0, slope, checkFlags.tolerance);
 
    if (checkFlags.writeSlopeStatementIfFailed and not checkPassed) {
      spdlog::info("Hessian check:");
      spdlog::info("The slope should be 3. It seems to be {}.", slope);
      if (checkPassed)
        spdlog::info("We consider this as sufficient.");
      else
        spdlog::info("The Hessian seems wrong.");
    }
    nonLinOp.updateAll();
    return checkPassed;
  }
}  // namespace Ikarus::utils