// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#include "main.h"
#include <contrib/Eigen/MatrixFunctions>

// Returns a matrix with eigenvalues clustered around 0, 1 and 2.
template <typename MatrixType>
MatrixType randomMatrixWithRealEivals(const Index size) {
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  MatrixType diag = MatrixType::Zero(size, size);
  for (Index i = 0; i < size; ++i) {
    diag(i, i) =
        Scalar(RealScalar(internal::random<int>(0, 2))) + internal::random<Scalar>() * Scalar(RealScalar(0.01));
  }
  MatrixType A = MatrixType::Random(size, size);
  HouseholderQR<MatrixType> QRofA(A);
  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}

template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
struct randomMatrixWithImagEivals {
  // Returns a matrix with eigenvalues clustered around 0 and +/- i.
  static MatrixType run(const Index size);
};

// Partial specialization for real matrices
template <typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 0> {
  static MatrixType run(const Index size) {
    typedef typename MatrixType::Scalar Scalar;
    MatrixType diag = MatrixType::Zero(size, size);
    Index i = 0;
    while (i < size) {
      Index randomInt = internal::random<Index>(-1, 1);
      if (randomInt == 0 || i == size - 1) {
        diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
        ++i;
      } else {
        Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
        diag(i, i + 1) = alpha;
        diag(i + 1, i) = -alpha;
        i += 2;
      }
    }
    MatrixType A = MatrixType::Random(size, size);
    HouseholderQR<MatrixType> QRofA(A);
    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  }
};

// Partial specialization for complex matrices
template <typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 1> {
  static MatrixType run(const Index size) {
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    const Scalar imagUnit(0, 1);
    MatrixType diag = MatrixType::Zero(size, size);
    for (Index i = 0; i < size; ++i) {
      diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit +
                   internal::random<Scalar>() * Scalar(RealScalar(0.01));
    }
    MatrixType A = MatrixType::Random(size, size);
    HouseholderQR<MatrixType> QRofA(A);
    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  }
};

template <typename MatrixType>
void testMatrixExponential(const MatrixType& A) {
  typedef typename internal::traits<MatrixType>::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef std::complex<RealScalar> ComplexScalar;

  VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
}

template <typename MatrixType>
void testMatrixLogarithm(const MatrixType& A) {
  typedef typename internal::traits<MatrixType>::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;

  MatrixType scaledA;
  RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
  if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
    scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
  else
    scaledA = A;

  // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
  MatrixType expA = scaledA.exp();
  MatrixType logExpA = expA.log();
  // The round trip loses accuracy proportional to the magnitude of the intermediate exp(A), which is unbounded
  // relative to A itself once A is close to zero.
  VERIFY_IS_APPROX_SCALED(logExpA, scaledA, max_abs_coeff(expA));
}

template <typename MatrixType>
void testHyperbolicFunctions(const MatrixType& A) {
  // exp(A) and exp(-A) can cancel, so compare against the magnitude of the terms rather than of the result.
  const MatrixType expA = A.exp();
  const MatrixType expmA = (-A).exp();
  const typename NumTraits<typename MatrixType::Scalar>::Real scale =
      numext::maxi(max_abs_coeff(expA), max_abs_coeff(expmA));
  VERIFY_IS_APPROX_SCALED(A.sinh(), (expA - expmA) / 2, scale);
  VERIFY_IS_APPROX_SCALED(A.cosh(), (expA + expmA) / 2, scale);
}

template <typename MatrixType>
void testGonioFunctions(const MatrixType& A) {
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef std::complex<RealScalar> ComplexScalar;
  typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, MatrixType::Options>
      ComplexMatrix;

  ComplexScalar imagUnit(0, 1);
  ComplexScalar two(2, 0);

  ComplexMatrix Ac = A.template cast<ComplexScalar>();

  ComplexMatrix exp_iA = (imagUnit * Ac).exp();
  ComplexMatrix exp_miA = (-imagUnit * Ac).exp();

  // As in testHyperbolicFunctions, the two exponentials can cancel.
  const RealScalar scale = numext::maxi(max_abs_coeff(exp_iA), max_abs_coeff(exp_miA));

  ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
  VERIFY_IS_APPROX_SCALED(sinAc, (exp_iA - exp_miA) / (two * imagUnit), scale);
  ComplexMatrix result(A.rows(), A.cols());
  MatrixBase<ComplexMatrix>& destination = result;
  Ac.sin().evalTo(destination);
  VERIFY_IS_APPROX_SCALED(result, (exp_iA - exp_miA) / (two * imagUnit), scale);

  ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
  VERIFY_IS_APPROX_SCALED(cosAc, (exp_iA + exp_miA) / 2, scale);
}

template <typename MatrixType>
void testMatrix(const MatrixType& A) {
  testMatrixExponential(A);
  testMatrixLogarithm(A);
  testHyperbolicFunctions(A);
  testGonioFunctions(A);
}

template <typename MatrixType>
void testMatrixType(const MatrixType& m) {
  // Matrices with clustered eigenvalue lead to different code paths
  // in MatrixFunction.h and are thus useful for testing.

  const Index size = m.rows();
  for (int i = 0; i < g_repeat; i++) {
    testMatrix(MatrixType::Random(size, size).eval());
    testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
    testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
  }
}

template <typename MatrixType>
void testMapRef(const MatrixType& A) {
  // Test if passing Ref and Map objects is possible
  // (Regression test for Bug #1796)
  Index size = A.rows();
  MatrixType X;
  X.setRandom(size, size);
  MatrixType Y(size, size);
  Ref<MatrixType> R(Y);
  Ref<const MatrixType> Rc(X);
  Map<MatrixType> M(Y.data(), size, size);
  Map<const MatrixType> Mc(X.data(), size, size);

  X = X * X;  // make sure sqrt is possible
  Y = X.sqrt();
  R = Rc.sqrt();
  M = Mc.sqrt();
  Y = X.exp();
  R = Rc.exp();
  M = Mc.exp();
  X = Y;  // make sure log is possible
  Y = X.log();
  R = Rc.log();
  M = Mc.log();

  Y = X.cos() + Rc.cos() + Mc.cos();
  Y = X.sin() + Rc.sin() + Mc.sin();

  Y = X.cosh() + Rc.cosh() + Mc.cosh();
  Y = X.sinh() + Rc.sinh() + Mc.sinh();
}

template <typename MatrixType>
void testMatrixLogarithmSingular() {
  // Regression test for bug #1613: log() of a singular matrix used to
  // hang because sqrt(0) = 0 is a fixed point of the inverse-scaling loop
  // in matrix_log_compute_big(). The implementation now returns NaN.
  MatrixType A = MatrixType::Zero(5, 5);
  VERIFY(A.log().array().isNaN().all());
  // An infinite entry of an atomic block also used to hang that loop.
  A.setIdentity();
  A(0, 4) = std::numeric_limits<typename MatrixType::RealScalar>::infinity();
  VERIFY(A.log().array().isNaN().all());
}

// The loop could also create an infinite entry: the first square root divides A(0, 1) by 2 sqrt(A(0, 0)) = 2e-10,
// which overflows, as log(A)(0, 1) = A(0, 1) / A(0, 0) does. The three eigenvalues form one atomic block.
void testMatrixLogarithmSquareRootOverflow() {
  MatrixXd A(3, 3);
  A << 1e-20, 1e300, 0, 0, 1e-20, 1, 0, 0, 1e-20;
  VERIFY(MatrixXd(A.log()).array().isNaN().all());
}

// An empty matrix has no eigenvalue clusters, but the start of the first block was written anyway.
template <typename MatrixType>
void testEmptyMatrix() {
  const MatrixType A(0, 0);
  VERIFY_IS_EQUAL(MatrixType(A.sin()).size(), 0);
  VERIFY_IS_EQUAL(MatrixType(A.log()).size(), 0);
}

EIGEN_DECLARE_TEST(matrix_function) {
  CALL_SUBTEST_1(testMatrixType(Matrix<float, 1, 1>()));
  CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
  CALL_SUBTEST_3(testMatrixType(MatrixXf(8, 8)));
  CALL_SUBTEST_4(testMatrixType(Matrix2d()));
  CALL_SUBTEST_5(testMatrixType(Matrix<double, 5, 5, RowMajor>()));
  CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
  CALL_SUBTEST_7(testMatrixType(MatrixXd(13, 13)));

  CALL_SUBTEST_1(testMapRef(Matrix<float, 1, 1>()));
  CALL_SUBTEST_2(testMapRef(Matrix3cf()));
  CALL_SUBTEST_3(testMapRef(MatrixXf(8, 8)));
  CALL_SUBTEST_7(testMapRef(MatrixXd(13, 13)));

  CALL_SUBTEST_3(testMatrixLogarithmSingular<MatrixXf>());
  CALL_SUBTEST_7(testMatrixLogarithmSingular<MatrixXd>());
  CALL_SUBTEST_7(testMatrixLogarithmSquareRootOverflow());
  CALL_SUBTEST_3(testEmptyMatrix<MatrixXf>());
  CALL_SUBTEST_7(testEmptyMatrix<MatrixXcd>());
}