#include <chrono>
#include <iostream>

#include <Eigen/Dense>

#include "fung/fung.hh"
#include "fung/examples/biomechanics/adipose_tissue_sommer_holzapfel.hh"
#include "fung/examples/biomechanics/muscle_tissue_martins.hh"

using namespace std::chrono;
using std::cout;
using std::endl;

template <class M, class Function>
void runTest(Function& f)
{
  auto numberOfEvaluations = 10'000'000u;

  M a = FunG::LinearAlgebra::unitMatrix<M>();
  M da(a);
  da *= 2;
  auto val = 0.;
  auto scale = 1.01;

  auto startTime = high_resolution_clock::now();
  for(auto i=0u; i<numberOfEvaluations; ++i)
  {
    (i%2==0) ? a*=scale : a*=1./scale; // manipulate a such that the compiler can't optimize the loop away
    f.update(a);
  }
  cout << "update: " << duration_cast<milliseconds>(high_resolution_clock::now() - startTime).count() << "ms\n";

  startTime = high_resolution_clock::now();
  for(auto i=0u; i<numberOfEvaluations; ++i)
  {
    (i%2==0) ? da *= scale : da *= 1./scale;
    val = f.d1(da);
  }
  cout << "d1: " << duration_cast<milliseconds>(high_resolution_clock::now() - startTime).count() << "ms\n";

  startTime = high_resolution_clock::now();
  for(auto i=0u; i<numberOfEvaluations; ++i)
  {
    (i%2==0) ? da *= scale : da *= 1./scale;
    val = f.d2(da,da);
  }
  cout << "d2: " << duration_cast<milliseconds>(high_resolution_clock::now() - startTime).count() << "ms\n";

  startTime = high_resolution_clock::now();
  for(auto i=0u; i<numberOfEvaluations; ++i)
  {
    (i%2==0) ? da *= scale : da *= 1./scale;
    val = f.d3(da,da,da);
  }
  cout << "d3: " << duration_cast<milliseconds>(high_resolution_clock::now() - startTime).count() << "ms\n";
  cout << endl;
}

int main()
{
  using FunG::LN;
  using FunG::Pow;
  constexpr int dim = 3;
  using M = Eigen::Matrix<double,dim,dim>;

  // material parameters
  auto d0 = 1., d1 = 1.;
  // fiberTensor: structural tensor, here associated with v=(1,0,0)
  // I: initial deformation, unit matrix corresponds to undeformed initial state
  auto fiberTensor = FunG::zero<M>(), I = FunG::zero<M>();
  I(0,0) = I(1,1) = I(2,2) = fiberTensor(0,0) = 1.;

  // models for adipose and muscle tissue
  auto compressibleAdipose = FunG::compressibleAdiposeTissue_SommerHolzapfel<Pow<2>,LN>(d0,d1,fiberTensor,I);
  auto compressibleMuscle = FunG::compressibleMuscleTissue_Martins<Pow<2>,LN>(d0,d1,fiberTensor,I);

  cout << "compressible adipose (Sommer,Holzapfel et al)" << endl;
  runTest<M>(compressibleAdipose);

  cout << "compressible muscle (Martins et al)" << endl;
  runTest<M>(compressibleMuscle);
}
