eigen/test/nullary.cpp

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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/.

#include "main.h"

template <typename MatrixType>
bool equalsIdentity(const MatrixType& A) {
  bool offDiagOK = true;
  for (Index i = 0; i < A.rows(); ++i) {
    for (Index j = i + 1; j < A.cols(); ++j) {
      offDiagOK = offDiagOK && numext::is_exactly_zero(A(i, j));
    }
  }
  for (Index i = 0; i < A.rows(); ++i) {
    for (Index j = 0; j < (std::min)(i, A.cols()); ++j) {
      offDiagOK = offDiagOK && numext::is_exactly_zero(A(i, j));
    }
  }

  bool diagOK = (A.diagonal().array() == 1).all();
  return offDiagOK && diagOK;
}

template <typename VectorType>
void check_extremity_accuracy(const VectorType& v, const typename VectorType::Scalar& low,
                              const typename VectorType::Scalar& high) {
  typedef typename VectorType::Scalar Scalar;
  typedef typename VectorType::RealScalar RealScalar;

  RealScalar prec = internal::is_same<RealScalar, float>::value ? NumTraits<RealScalar>::dummy_precision() * 10
                                                                : NumTraits<RealScalar>::dummy_precision() / 10;
  Index size = v.size();

  if (size < 20) return;

  for (int i = 0; i < size; ++i) {
    if (i < 5 || i > size - 6) {
      Scalar ref =
          (low * RealScalar(size - i - 1)) / RealScalar(size - 1) + (high * RealScalar(i)) / RealScalar(size - 1);
      if (std::abs(ref) > 1) {
        if (!internal::isApprox(v(i), ref, prec))
          std::cout << v(i) << " != " << ref << "  ; relative error: " << std::abs((v(i) - ref) / ref)
                    << "  ; required precision: " << prec << "  ; range: " << low << "," << high << "  ; i: " << i
                    << "\n";
        VERIFY(internal::isApprox(
            v(i),
            (low * RealScalar(size - i - 1)) / RealScalar(size - 1) + (high * RealScalar(i)) / RealScalar(size - 1),
            prec));
      }
    }
  }
}

template <typename VectorType>
void testVectorType(const VectorType& base) {
  typedef typename VectorType::Scalar Scalar;
  typedef typename VectorType::RealScalar RealScalar;

  const Index size = base.size();

  Scalar high = internal::random<Scalar>(-500, 500);
  Scalar low = (size == 1 ? high : internal::random<Scalar>(-500, 500));
  if (numext::real(low) > numext::real(high)) std::swap(low, high);

  // check low==high
  if (internal::random<float>(0.f, 1.f) < 0.05f) low = high;
  // check abs(low) >> abs(high)
  else if (size > 2 && std::numeric_limits<RealScalar>::max_exponent10 > 0 && internal::random<float>(0.f, 1.f) < 0.1f)
    low = -internal::random<Scalar>(1, 2) *
          RealScalar(std::pow(RealScalar(10), std::numeric_limits<RealScalar>::max_exponent10 / 2));

  const Scalar step = ((size == 1) ? 1 : (high - low) / RealScalar(size - 1));

  // check whether the result yields what we expect it to do
  VectorType m(base), o(base);
  m.setLinSpaced(size, low, high);
  o.setEqualSpaced(size, low, step);

  if (!NumTraits<Scalar>::IsInteger) {
    VectorType n(size);
    for (int i = 0; i < size; ++i) n(i) = low + RealScalar(i) * step;
    VERIFY_IS_APPROX(m, n);
    VERIFY_IS_APPROX(n, o);

    CALL_SUBTEST(check_extremity_accuracy(m, low, high));
  }

  RealScalar range_length = numext::real(high - low);
  if ((!NumTraits<Scalar>::IsInteger) || (range_length >= size && (Index(range_length) % (size - 1)) == 0) ||
      (Index(range_length + 1) < size && (size % Index(range_length + 1)) == 0)) {
    VectorType n(size);
    if ((!NumTraits<Scalar>::IsInteger) || (range_length >= size))
      for (int i = 0; i < size; ++i) n(i) = size == 1 ? low : (low + ((high - low) * Scalar(i)) / RealScalar(size - 1));
    else
      for (int i = 0; i < size; ++i)
        n(i) = size == 1 ? low : low + Scalar((double(range_length + 1) * double(i)) / double(size));
    VERIFY_IS_APPROX(m, n);

    // random access version
    m = VectorType::LinSpaced(size, low, high);
    VERIFY_IS_APPROX(m, n);
    VERIFY(internal::isApprox(m(m.size() - 1), high));
    VERIFY(size == 1 || internal::isApprox(m(0), low));
    VERIFY_IS_EQUAL(m(m.size() - 1), high);
    if (!NumTraits<Scalar>::IsInteger) CALL_SUBTEST(check_extremity_accuracy(m, low, high));
  }

  VERIFY(numext::real(m(m.size() - 1)) <= numext::real(high));
  VERIFY((m.array().real() <= numext::real(high)).all());
  VERIFY((m.array().real() >= numext::real(low)).all());

  VERIFY(numext::real(m(m.size() - 1)) >= numext::real(low));
  if (size >= 1) {
    VERIFY(internal::isApprox(m(0), low));
    VERIFY_IS_EQUAL(m(0), low);
  }

  // check whether everything works with row and col major vectors
  Matrix<Scalar, Dynamic, 1> row_vector(size);
  Matrix<Scalar, 1, Dynamic> col_vector(size);
  row_vector.setLinSpaced(size, low, high);
  col_vector.setLinSpaced(size, low, high);
  // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit
  // when computing the squared sum in isApprox, thus the 2x factor.
  VERIFY(row_vector.isApprox(col_vector.transpose(), RealScalar(2) * NumTraits<Scalar>::epsilon()));

  Matrix<Scalar, Dynamic, 1> size_changer(size + 50);
  size_changer.setLinSpaced(size, low, high);
  VERIFY(size_changer.size() == size);

  typedef Matrix<Scalar, 1, 1> ScalarMatrix;
  ScalarMatrix scalar;
  scalar.setLinSpaced(1, low, high);
  VERIFY_IS_APPROX(scalar, ScalarMatrix::Constant(high));
  VERIFY_IS_APPROX(ScalarMatrix::LinSpaced(1, low, high), ScalarMatrix::Constant(high));

  // regression test for bug 526 (linear vectorized transversal)
  if (size > 1 && (!NumTraits<Scalar>::IsInteger)) {
    m.tail(size - 1).setLinSpaced(low, high);
    VERIFY_IS_APPROX(m(size - 1), high);
  }

  // regression test for bug 1383 (LinSpaced with empty size/range)
  {
    Index n0 = VectorType::SizeAtCompileTime == Dynamic ? 0 : VectorType::SizeAtCompileTime;
    low = internal::random<Scalar>();
    m = VectorType::LinSpaced(n0, low, low - RealScalar(1));
    VERIFY(m.size() == n0);

    if (VectorType::SizeAtCompileTime == Dynamic) {
      VERIFY_IS_EQUAL(VectorType::LinSpaced(n0, 0, Scalar(RealScalar(n0 - 1))).sum(), Scalar(0));
      VERIFY_IS_EQUAL(VectorType::LinSpaced(n0, low, low - RealScalar(1)).sum(), Scalar(0));
    }

    m.setLinSpaced(n0, 0, Scalar(RealScalar(n0 - 1)));
    VERIFY(m.size() == n0);
    m.setLinSpaced(n0, low, low - RealScalar(1));
    VERIFY(m.size() == n0);

    // empty range only:
    VERIFY_IS_APPROX(VectorType::LinSpaced(size, low, low), VectorType::Constant(size, low));
    m.setLinSpaced(size, low, low);
    VERIFY_IS_APPROX(m, VectorType::Constant(size, low));

    if (NumTraits<Scalar>::IsInteger) {
      VERIFY_IS_APPROX(VectorType::LinSpaced(size, low, low + Scalar(RealScalar(size - 1))),
                       VectorType::LinSpaced(size, low + Scalar(RealScalar(size - 1)), low).reverse());

      if (VectorType::SizeAtCompileTime == Dynamic) {
        // Check negative multiplicator path:
        for (Index k = 1; k < 5; ++k)
          VERIFY_IS_APPROX(VectorType::LinSpaced(size, low, low + Scalar(RealScalar((size - 1) * k))),
                           VectorType::LinSpaced(size, low + Scalar(RealScalar((size - 1) * k)), low).reverse());
        // Check negative divisor path:
        for (Index k = 1; k < 5; ++k)
          VERIFY_IS_APPROX(VectorType::LinSpaced(size * k, low, low + Scalar(RealScalar(size - 1))),
                           VectorType::LinSpaced(size * k, low + Scalar(RealScalar(size - 1)), low).reverse());
      }
    }
  }

  // test setUnit()
  if (m.size() > 0) {
    for (Index k = 0; k < 10; ++k) {
      Index i = internal::random<Index>(0, m.size() - 1);
      m.setUnit(i);
      VERIFY_IS_APPROX(m, VectorType::Unit(m.size(), i));
    }
    if (VectorType::SizeAtCompileTime == Dynamic) {
      Index i = internal::random<Index>(0, 2 * m.size() - 1);
      m.setUnit(2 * m.size(), i);
      VERIFY_IS_APPROX(m, VectorType::Unit(m.size(), i));
    }
  }
}

template <typename MatrixType>
void testMatrixType(const MatrixType& m) {
  using std::abs;
  const Index rows = m.rows();
  const Index cols = m.cols();
  typedef typename MatrixType::Scalar Scalar;

  Scalar s1 = internal::random<Scalar>();

  MatrixType A;
  A.setIdentity(rows, cols);
  VERIFY(equalsIdentity(A));
  VERIFY(equalsIdentity(MatrixType::Identity(rows, cols)));

  A = MatrixType::Constant(rows, cols, s1);
  Index i = internal::random<Index>(0, rows - 1);
  Index j = internal::random<Index>(0, cols - 1);
  VERIFY_IS_APPROX(MatrixType::Constant(rows, cols, s1)(i, j), s1);
  VERIFY_IS_APPROX(MatrixType::Constant(rows, cols, s1).coeff(i, j), s1);
  VERIFY_IS_APPROX(A(i, j), s1);
}

template <int>
void bug79() {
  // Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
  VERIFY((MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() <
         std::numeric_limits<double>::epsilon());
}

template <int>
void bug1630() {
  Array4d x4 = Array4d::LinSpaced(0.0, 1.0);
  Array3d x3(Array4d::LinSpaced(0.0, 1.0).head(3));
  VERIFY_IS_APPROX(x4.head(3), x3);
}

template <int>
void nullary_overflow() {
  // Check possible overflow issue
  int n = 60000;
  ArrayXi a1(n), a2(n), a_ref(n);
  a1.setLinSpaced(n, 0, n - 1);
  a2.setEqualSpaced(n, 0, 1);
  for (int i = 0; i < n; ++i) a_ref(i) = i;
  VERIFY_IS_APPROX(a1, a_ref);
  VERIFY_IS_APPROX(a2, a_ref);
}

template <int>
void nullary_internal_logic() {
  // check some internal logic
  VERIFY((internal::has_nullary_operator<internal::scalar_constant_op<double> >::value));
  VERIFY((!internal::has_unary_operator<internal::scalar_constant_op<double> >::value));
  VERIFY((!internal::has_binary_operator<internal::scalar_constant_op<double> >::value));
  VERIFY((internal::functor_has_linear_access<internal::scalar_constant_op<double> >::ret));

  VERIFY((!internal::has_nullary_operator<internal::scalar_identity_op<double> >::value));
  VERIFY((!internal::has_unary_operator<internal::scalar_identity_op<double> >::value));
  VERIFY((internal::has_binary_operator<internal::scalar_identity_op<double> >::value));
  VERIFY((!internal::functor_has_linear_access<internal::scalar_identity_op<double> >::ret));

  VERIFY((!internal::has_nullary_operator<internal::linspaced_op<float> >::value));
  VERIFY((internal::has_unary_operator<internal::linspaced_op<float> >::value));
  VERIFY((!internal::has_binary_operator<internal::linspaced_op<float> >::value));
  VERIFY((internal::functor_has_linear_access<internal::linspaced_op<float> >::ret));

  // Regression unit test for an MSVC bug.
  // Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details.
  // See also traits<Ref>::match.
  {
    MatrixXf A = MatrixXf::Random(3, 3);
    Ref<const MatrixXf> R = 2.0 * A;
    VERIFY_IS_APPROX(R, A + A);

    Ref<const MatrixXf> R1 = MatrixXf::Random(3, 3) + A;

    VectorXi V = VectorXi::Random(3);
    Ref<const VectorXi> R2 = VectorXi::LinSpaced(3, 1, 3) + V;
    VERIFY_IS_APPROX(R2, V + Vector3i(1, 2, 3));

    VERIFY((internal::has_nullary_operator<internal::scalar_constant_op<float> >::value));
    VERIFY((!internal::has_unary_operator<internal::scalar_constant_op<float> >::value));
    VERIFY((!internal::has_binary_operator<internal::scalar_constant_op<float> >::value));
    VERIFY((internal::functor_has_linear_access<internal::scalar_constant_op<float> >::ret));

    VERIFY((!internal::has_nullary_operator<internal::linspaced_op<int> >::value));
    VERIFY((internal::has_unary_operator<internal::linspaced_op<int> >::value));
    VERIFY((!internal::has_binary_operator<internal::linspaced_op<int> >::value));
    VERIFY((internal::functor_has_linear_access<internal::linspaced_op<int> >::ret));
  }
}

// Test LinSpaced at vectorization boundary sizes.
// The packetOp in linspaced_op_impl uses mask/select logic to handle
// the last partial packet (when vector size is not a multiple of PacketSize).
// This exercises those boundaries with element-by-element verification.
template <typename Scalar>
void linspaced_boundary() {
  typedef typename NumTraits<Scalar>::Real RealScalar;
  const Index PS = internal::packet_traits<Scalar>::size;
  const Index sizes[] = {1, 2, 3, PS - 1, PS, PS + 1, 2 * PS - 1, 2 * PS, 2 * PS + 1, 4 * PS, 4 * PS + 1};
  typedef Matrix<Scalar, Dynamic, 1> Vec;

  for (int si = 0; si < 11; ++si) {
    Index n = sizes[si];
    if (n <= 0) continue;

    Scalar low(1), high(100);
    Vec v = Vec::LinSpaced(n, low, high);

    // With n==1, LinSpaced returns [high] by design.
    if (n == 1) {
      VERIFY_IS_EQUAL(v(0), high);
    } else {
      VERIFY_IS_EQUAL(v(0), low);
      VERIFY_IS_EQUAL(v(n - 1), high);

      // Verify monotonicity.
      for (Index k = 1; k < n; ++k) {
        VERIFY(numext::real(v(k)) >= numext::real(v(k - 1)));
      }

      // Verify against scalar reference computation.
      for (Index k = 0; k < n; ++k) {
        Scalar ref = Scalar(low + (high - low) * RealScalar(k) / RealScalar(n - 1));
        VERIFY_IS_APPROX(v(k), ref);
      }
    }
  }

  // Test the "flip" path: when |high| < |low|, the implementation uses
  // a reversed computation for better precision. Verify at packet boundaries.
  for (int si = 0; si < 11; ++si) {
    Index n = sizes[si];
    if (n <= 0 || n == 1) continue;  // skip n=1, flip irrelevant for single element

    Scalar low(1000), high(1);  // |high| < |low| triggers flip
    Vec v = Vec::LinSpaced(n, low, high);

    VERIFY_IS_EQUAL(v(0), low);
    VERIFY_IS_EQUAL(v(n - 1), high);

    // Verify monotonicity (decreasing).
    for (Index k = 1; k < n; ++k) {
      VERIFY(numext::real(v(k)) <= numext::real(v(k - 1)));
    }

    // Verify against scalar reference.
    for (Index k = 0; k < n; ++k) {
      Scalar ref = Scalar(low + (high - low) * RealScalar(k) / RealScalar(n - 1));
      VERIFY_IS_APPROX(v(k), ref);
    }
  }
}

// Test integer LinSpaced divisor path.
// When (abs(high - low) + 1) < num_steps, the integer LinSpaced uses
// a divisor-based formula instead of multiplication. This path is
// barely covered by existing tests which use random ranges.
template <int>
void linspaced_integer_divisor() {
  typedef Matrix<int, Dynamic, 1> VecI;

  // Case: num_steps much larger than range → triggers divisor path.
  // LinSpaced(12, 0, 5): 12 steps over range [0,5], so range+1=6, 6 < 12.
  {
    VecI v = VecI::LinSpaced(12, 0, 5);
    VERIFY_IS_EQUAL(v(0), 0);
    // All values must be in [0, 5].
    for (Index k = 0; k < 12; ++k) {
      VERIFY(v(k) >= 0 && v(k) <= 5);
    }
    // Must be non-decreasing.
    for (Index k = 1; k < 12; ++k) {
      VERIFY(v(k) >= v(k - 1));
    }
    // Each integer 0-5 should appear at least once.
    for (int val = 0; val <= 5; ++val) {
      bool found = false;
      for (Index k = 0; k < 12; ++k) {
        if (v(k) == val) {
          found = true;
          break;
        }
      }
      VERIFY(found);
    }
  }

  // Case: range exactly divides steps → each value should appear equally.
  // LinSpaced(20, 0, 3): range+1=4, 20%4==0, so each of 0,1,2,3 appears 5 times.
  {
    VecI v = VecI::LinSpaced(20, 0, 3);
    VERIFY_IS_EQUAL(v(0), 0);
    for (Index k = 0; k < 20; ++k) {
      VERIFY(v(k) >= 0 && v(k) <= 3);
    }
    for (Index k = 1; k < 20; ++k) {
      VERIFY(v(k) >= v(k - 1));
    }
  }

  // Reverse: LinSpaced(12, 5, 0) should be reverse of LinSpaced(12, 0, 5).
  {
    VecI fwd = VecI::LinSpaced(12, 0, 5);
    VecI rev = VecI::LinSpaced(12, 5, 0);
    VERIFY_IS_APPROX(fwd, rev.reverse());
  }

  // Single step: always returns high.
  {
    VecI v = VecI::LinSpaced(1, 3, 7);
    VERIFY_IS_EQUAL(v(0), 7);
  }
}

EIGEN_DECLARE_TEST(nullary) {
  CALL_SUBTEST_1(testMatrixType(Matrix2d()));
  CALL_SUBTEST_2(testMatrixType(MatrixXcf(internal::random<int>(1, 300), internal::random<int>(1, 300))));
  CALL_SUBTEST_3(testMatrixType(MatrixXf(internal::random<int>(1, 300), internal::random<int>(1, 300))));

  for (int i = 0; i < g_repeat * 10; i++) {
    CALL_SUBTEST_3(testVectorType(VectorXcd(internal::random<int>(1, 30000))));
    CALL_SUBTEST_4(testVectorType(VectorXd(internal::random<int>(1, 30000))));
    CALL_SUBTEST_5(testVectorType(Vector4d()));  // regression test for bug 232
    CALL_SUBTEST_6(testVectorType(Vector3d()));
    CALL_SUBTEST_7(testVectorType(VectorXf(internal::random<int>(1, 30000))));
    CALL_SUBTEST_8(testVectorType(Vector3f()));
    CALL_SUBTEST_8(testVectorType(Vector4f()));
    CALL_SUBTEST_8(testVectorType(Matrix<float, 8, 1>()));
    CALL_SUBTEST_8(testVectorType(Matrix<float, 1, 1>()));

    CALL_SUBTEST_9(testVectorType(VectorXi(internal::random<int>(1, 10))));
    CALL_SUBTEST_9(testVectorType(VectorXi(internal::random<int>(9, 300))));
    CALL_SUBTEST_9(testVectorType(Matrix<int, 1, 1>()));
  }

  CALL_SUBTEST_6(bug79<0>());
  CALL_SUBTEST_6(bug1630<0>());
  CALL_SUBTEST_9(nullary_overflow<0>());
  CALL_SUBTEST_10(nullary_internal_logic<0>());

  // LinSpaced at vectorization boundaries (deterministic, outside g_repeat).
  CALL_SUBTEST_11(linspaced_boundary<float>());
  CALL_SUBTEST_11(linspaced_boundary<double>());

  // Integer LinSpaced divisor path tests.
  CALL_SUBTEST_12(linspaced_integer_divisor<0>());
}