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);
  m.setLinSpaced(size,low,high);

  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);

    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(n0-1)).sum(),Scalar(0));
      VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,low,low-RealScalar(1)).sum(),Scalar(0));
    }

    m.setLinSpaced(n0,0,Scalar(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(size-1)), VectorType::LinSpaced(size,low+Scalar(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((size-1)*k)), VectorType::LinSpaced(size,low+Scalar((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(size-1)), VectorType::LinSpaced(size*k,low+Scalar(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;
  typedef typename MatrixType::RealScalar RealScalar;

  Scalar s1;
  do {
    s1 = internal::random<Scalar>();
  } while(abs(s1)<RealScalar(1e-5) && (!NumTraits<Scalar>::IsInteger));

  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);
  a1.setLinSpaced(n, 0, n-1);
  for(int i=0; i<n; ++i)
    a2(i) = i;
  VERIFY_IS_APPROX(a1,a2);
}

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 a weird 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 ));
  }
}

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>() );
}