#include <cassert>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fufem/interval.hh>

#include <dune/tnnmg/problem-classes/bisection.hh>

#include <dune/tectonic/mydirectionalconvexfunction.hh>

#include "compute_state_dieterich.hh"

// psi(beta) = V/L beta + e^(-beta)
class DieterichNonlinearity {
public:
  typedef Dune::FieldVector<double, 1> VectorType;
  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;

  DieterichNonlinearity(double _VoL) : VoL(_VoL) {}

  void directionalSubDiff(VectorType const &u, VectorType const &v,
                          Interval<double> &D) const {
    D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
  }

  void directionalDomain(VectorType const &, VectorType const &,
                         Interval<double> &dom) const {
    dom[0] = -std::numeric_limits<double>::max();
    dom[1] = std::numeric_limits<double>::max();
  }

private:
  double const VoL;
};

double state_update_dieterich_bisection(double tau, double VoL,
                                        double old_state) {
  DieterichNonlinearity::VectorType const start(0);
  DieterichNonlinearity::VectorType const direction(1);
  DieterichNonlinearity const phi(VoL);
  MyDirectionalConvexFunction<DieterichNonlinearity> const J(
      1.0 / tau, old_state / tau, phi, start, direction);

  int bisectionsteps = 0;
  Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
}

double state_update_dieterich(double tau, double VoL, double old_state) {
  double const ret = state_update_dieterich_bisection(tau, VoL, old_state);
  /* We have

     ret - old_state + tau * VoL = tau*exp(-ret)

     or

     log((ret - old_state)/tau + VoL) = -ret
  */
  assert(std::min(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))),
                  std::abs(std::log((ret - old_state) / tau + VoL) + ret)) <
         1e-12);
  return ret;
}