#include <cmath>

#include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/Math.hpp>

#include "Tests.hpp"

using namespace AMDiS;
using namespace Dune;

// -----------------------------------------------------------------------------
void test0()
{
  using V = FieldVector<double, 3>;
  V a{1.0, 2.0, 3.0};
  V b{2.0, 3.0, 4.0};

  AMDIS_TEST_EQ( a+b, (V{3.0, 5.0, 7.0}) );
  AMDIS_TEST_EQ( a-b, (V{-1.0, -1.0, -1.0}) );

  AMDIS_TEST_EQ( a*2,   (V{2.0, 4.0, 6.0}) );
  AMDIS_TEST_EQ( 2.0*a, (V{2.0, 4.0, 6.0}) );

  AMDIS_TEST_EQ( b/2.0, (V{1.0, 1.5, 2.0}) );

  AMDIS_TEST_EQ( a + 2*b, (V{5.0, 8.0, 11.0}) );
}

// -----------------------------------------------------------------------------
void test1()
{
  using V = FieldVector<double, 3>;
  V a{1.0, 2.0, 3.0};
  V b{2.0, 3.0, 4.0};
  V c{3.0, -5.0, 4.0};
  V d{-1.0, -2.0, -3.0};

  AMDIS_TEST_EQ( unary_dot(a), 14.0 );
  AMDIS_TEST_EQ( dot(a, b), 20.0 );
  AMDIS_TEST_EQ( dot(a, b), dot(b, a) );

  AMDIS_TEST_APPROX( one_norm(a), 6.0 );
  AMDIS_TEST_APPROX( two_norm(a), std::sqrt(14.0) );
  AMDIS_TEST_APPROX( two_norm(a), p_norm<2>(a) );
  AMDIS_TEST_APPROX( p_norm<3>(a), std::pow(36.0, 1.0/3) );
  AMDIS_TEST_APPROX( distance(a, b), std::sqrt(3.0) );

  AMDIS_TEST_EQ( cross(a, b), (V{-1.0, 2.0, -1.0}) );

  AMDIS_TEST_EQ( sum(c),  2.0 );
  AMDIS_TEST_EQ( min(b),  2.0 );
  AMDIS_TEST_EQ( min(c), -5.0 );
  AMDIS_TEST_EQ( max(c),  4.0 );
  AMDIS_TEST_EQ( max(d), -1.0 );
  AMDIS_TEST_EQ( abs_min(c), 3.0 );
  AMDIS_TEST_EQ( abs_max(c), 5.0 );

  using V3 = FieldVector<int, 3>;
  V3 f{2, 3, 4};
  V3 g{3, -5, 4};
  V3 h{-1, -2, -3};
  AMDIS_TEST_EQ( sum(g),  2 );
  AMDIS_TEST_EQ( min(f),  2 );
  AMDIS_TEST_EQ( min(g), -5 );
  AMDIS_TEST_EQ( max(g),  4 );
  AMDIS_TEST_EQ( max(h), -1 );
  AMDIS_TEST_EQ( abs_min(g), 3 );
  AMDIS_TEST_EQ( abs_max(g), 5 );
}

// -----------------------------------------------------------------------------
void test2()
{
  using V = FieldVector<double, 3>;
  using M = FieldMatrix<double, 3, 3>;

  M A{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} };
  M B{ {1.0, 1.0, 1.0}, {2.0, 1.0, 0.0}, {3.0, 1.0, -1.0} };
  M C{ {0.5, 2.0, -1.0}, {1.0, -1.0, 1.0}, {-0.25, 0.5, 4.0} };
  M D{ {-0.25, 1.0, -1.0}, {2.0, -2.0, 8.0}, {-0.125, -0.125, 2.0} };

  AMDIS_TEST_EQ( A+B,   (M{ {2.0, 1.0, 1.0}, { 2.0, 2.0, 0.0}, { 3.0, 1.0, 0.0} }) );
  AMDIS_TEST_EQ( A-B,   (M{ {0.0,-1.0,-1.0}, {-2.0, 0.0, 0.0}, {-3.0,-1.0, 2.0} }) );

  AMDIS_TEST_EQ( A*2,   (M{ {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0} }) );
  AMDIS_TEST_EQ( 2.0*A, (M{ {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0} }) );

  AMDIS_TEST_EQ( A/2.0, (M{ {0.5, 0.0, 0.0}, {0.0, 0.5, 0.0}, {0.0, 0.0, 0.5} }) );


  V a{1.0, 2.0, 3.0};
  V b{2.0, 3.0, 4.0};

  // AMDIS_TEST_EQ( outer(a, b), (M{ {2.0, 3.0, 4.0}, {4.0, 6.0, 8.0}, {6.0, 9.0, 12.0} }) );
  // AMDIS_TEST_EQ( diagonal(a), (M{ {1.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 3.0} }) );

  // AMDIS_TEST_EQ( diagonal(A), (V{1.0, 1.0, 1.0}) );

  // some tests for the determinant calculation
  // ------------------------------------------
  AMDIS_TEST_EQ( det(A), 1.0 );
  AMDIS_TEST_EQ( det(B), 0.0 );

  AMDIS_TEST_EQ( det(multiplies(C,D)), det(C)*det(D) );
  AMDIS_TEST_EQ( det(trans(C)), det(C) );
  AMDIS_TEST_EQ( det(2.0*D), Math::pow<3>(2.0)*det(D) );

  // some tests for the trace
  // ------------------------
  // AMDIS_TEST_EQ( trace(D), trace(trans(D)) );
  // AMDIS_TEST_EQ( trace(2.0*C + 0.5*D), 2.0*trace(C) + 0.5*trace(D) );
  // AMDIS_TEST_EQ( trace(C*D), trace(D*C) );
  // AMDIS_TEST_EQ( trace(C), sum(diagonal(C)) );

  // some tests for the norm
  // -----------------------
  // AMDIS_TEST_APPROX( sqr(frobenius_norm(B)), trace(hermitian(B)*B) );
  // AMDIS_TEST( frobenius_norm(C*D) <= frobenius_norm(B)*frobenius_norm(C) );

  // test the inverse
  // ----------------
  // AMDIS_TEST_EQ( inv(C), adj(C)/det(C) );
  AMDIS_TEST_EQ( inv(A), A );
  AMDIS_TEST_EQ( inv(D), M({ {4.0/5, 1.0/2, -8.0/5}, {4.0/3, 1.0/6, 0.0}, {2.0/15, 1.0/24, 2.0/5} }) );

  V sol;
  solve(A,sol,a);
  AMDIS_TEST_EQ( sol, a );
}

// test of scalar wrapper FieldVector<T,1> and FieldMatrix<T,1,1>
void test3()
{
  using V = FieldVector<double, 3>;
  using M = FieldMatrix<double, 3, 3>;

  V a{1.0, 2.0, 3.0};
  M A{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} };

  using VD = FieldVector<double, 1>;
  using MD = FieldMatrix<double, 1, 1>;

  VD vd1 = 1.0;
  MD md1 = 1.0;

  using VI = FieldVector<int, 1>;
  using MI = FieldMatrix<int, 1, 1>;

  VI vi1 = 1;
  MI mi1 = 1;

  // scale a vector
  AMDIS_TEST_EQ( 1*a, a );
  AMDIS_TEST_EQ( 1.0*a, a );
  AMDIS_TEST_EQ( a*1, a );
  AMDIS_TEST_EQ( a*1.0, a );
  AMDIS_TEST_EQ( vd1*a, a );
  AMDIS_TEST_EQ( a*vd1, a );
  AMDIS_TEST_EQ( vi1*a, a );
  AMDIS_TEST_EQ( a*vi1, a );
  AMDIS_TEST_EQ( md1*a, a );
  AMDIS_TEST_EQ( a*md1, a );
  AMDIS_TEST_EQ( mi1*a, a );
  AMDIS_TEST_EQ( a*mi1, a );
  AMDIS_TEST_EQ( a/1, a );
  AMDIS_TEST_EQ( a/1.0, a );
  AMDIS_TEST_EQ( a/vd1, a );
  AMDIS_TEST_EQ( a/vi1, a );
  AMDIS_TEST_EQ( a/md1, a );
  AMDIS_TEST_EQ( a/mi1, a );

  // scale a matrix
  AMDIS_TEST_EQ( 1*A, A );
  AMDIS_TEST_EQ( 1.0*A, A );
  AMDIS_TEST_EQ( A*1, A );
  AMDIS_TEST_EQ( A*1.0, A );
  AMDIS_TEST_EQ( vd1*A, A );
  AMDIS_TEST_EQ( A*vd1, A );
  AMDIS_TEST_EQ( vi1*A, A );
  AMDIS_TEST_EQ( A*vi1, A );
  AMDIS_TEST_EQ( md1*A, A );
  AMDIS_TEST_EQ( A*md1, A );
  AMDIS_TEST_EQ( mi1*A, A );
  AMDIS_TEST_EQ( A*mi1, A );
  AMDIS_TEST_EQ( A/1, A );
  AMDIS_TEST_EQ( A/1.0, A );
  AMDIS_TEST_EQ( A/vd1, A );
  AMDIS_TEST_EQ( A/vi1, A );
  AMDIS_TEST_EQ( A/md1, A );
  AMDIS_TEST_EQ( A/mi1, A );

  AMDIS_TEST_EQ( vd1*vd1, 1.0 );
  AMDIS_TEST_EQ( vd1*md1, 1.0 );
  AMDIS_TEST_EQ( md1*md1, 1.0 );
  AMDIS_TEST_EQ( md1*vd1, 1.0 );
  AMDIS_TEST_EQ( vd1*1.0, 1.0 );
  AMDIS_TEST_EQ( md1*1.0, 1.0 );
  AMDIS_TEST_EQ( 1.0*md1, 1.0 );
  AMDIS_TEST_EQ( 1.0*vd1, 1.0 );
  AMDIS_TEST_EQ( vi1*vi1, 1 );
  AMDIS_TEST_EQ( vi1*mi1, 1 );
  AMDIS_TEST_EQ( mi1*mi1, 1 );
  AMDIS_TEST_EQ( mi1*vi1, 1 );
}


int main()
{
  test0();
  test1();
  test2();
  test3();

  return report_errors();
}