Home | History | Annotate | Download | only in eigen2 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra. Eigen itself is part of the KDE project.
      3 //
      4 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro (at) gmail.com>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #include "sparse.h"
     11 
     12 template<typename SetterType,typename DenseType, typename Scalar, int Options>
     13 bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
     14 {
     15   typedef SparseMatrix<Scalar,Options> SparseType;
     16   {
     17     sm.setZero();
     18     SetterType w(sm);
     19     std::vector<Vector2i> remaining = nonzeroCoords;
     20     while(!remaining.empty())
     21     {
     22       int i = ei_random<int>(0,remaining.size()-1);
     23       w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
     24       remaining[i] = remaining.back();
     25       remaining.pop_back();
     26     }
     27   }
     28   return sm.isApprox(ref);
     29 }
     30 
     31 template<typename SetterType,typename DenseType, typename T>
     32 bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
     33 {
     34   sm.setZero();
     35   std::vector<Vector2i> remaining = nonzeroCoords;
     36   while(!remaining.empty())
     37   {
     38     int i = ei_random<int>(0,remaining.size()-1);
     39     sm.coeffRef(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
     40     remaining[i] = remaining.back();
     41     remaining.pop_back();
     42   }
     43   return sm.isApprox(ref);
     44 }
     45 
     46 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
     47 {
     48   const int rows = ref.rows();
     49   const int cols = ref.cols();
     50   typedef typename SparseMatrixType::Scalar Scalar;
     51   enum { Flags = SparseMatrixType::Flags };
     52 
     53   double density = std::max(8./(rows*cols), 0.01);
     54   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
     55   typedef Matrix<Scalar,Dynamic,1> DenseVector;
     56   Scalar eps = 1e-6;
     57 
     58   SparseMatrixType m(rows, cols);
     59   DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
     60   DenseVector vec1 = DenseVector::Random(rows);
     61   Scalar s1 = ei_random<Scalar>();
     62 
     63   std::vector<Vector2i> zeroCoords;
     64   std::vector<Vector2i> nonzeroCoords;
     65   initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
     66 
     67   if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
     68     return;
     69 
     70   // test coeff and coeffRef
     71   for (int i=0; i<(int)zeroCoords.size(); ++i)
     72   {
     73     VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
     74     if(ei_is_same_type<SparseMatrixType,SparseMatrix<Scalar,Flags> >::ret)
     75       VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 );
     76   }
     77   VERIFY_IS_APPROX(m, refMat);
     78 
     79   m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
     80   refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
     81 
     82   VERIFY_IS_APPROX(m, refMat);
     83   /*
     84   // test InnerIterators and Block expressions
     85   for (int t=0; t<10; ++t)
     86   {
     87     int j = ei_random<int>(0,cols-1);
     88     int i = ei_random<int>(0,rows-1);
     89     int w = ei_random<int>(1,cols-j-1);
     90     int h = ei_random<int>(1,rows-i-1);
     91 
     92 //     VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
     93     for(int c=0; c<w; c++)
     94     {
     95       VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
     96       for(int r=0; r<h; r++)
     97       {
     98 //         VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
     99       }
    100     }
    101 //     for(int r=0; r<h; r++)
    102 //     {
    103 //       VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
    104 //       for(int c=0; c<w; c++)
    105 //       {
    106 //         VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
    107 //       }
    108 //     }
    109   }
    110 
    111   for(int c=0; c<cols; c++)
    112   {
    113     VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c));
    114     VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c));
    115   }
    116 
    117   for(int r=0; r<rows; r++)
    118   {
    119     VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
    120     VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
    121   }
    122   */
    123 
    124   // test SparseSetters
    125   // coherent setter
    126   // TODO extend the MatrixSetter
    127 //   {
    128 //     m.setZero();
    129 //     VERIFY_IS_NOT_APPROX(m, refMat);
    130 //     SparseSetter<SparseMatrixType, FullyCoherentAccessPattern> w(m);
    131 //     for (int i=0; i<nonzeroCoords.size(); ++i)
    132 //     {
    133 //       w->coeffRef(nonzeroCoords[i].x(),nonzeroCoords[i].y()) = refMat.coeff(nonzeroCoords[i].x(),nonzeroCoords[i].y());
    134 //     }
    135 //   }
    136 //   VERIFY_IS_APPROX(m, refMat);
    137 
    138   // random setter
    139 //   {
    140 //     m.setZero();
    141 //     VERIFY_IS_NOT_APPROX(m, refMat);
    142 //     SparseSetter<SparseMatrixType, RandomAccessPattern> w(m);
    143 //     std::vector<Vector2i> remaining = nonzeroCoords;
    144 //     while(!remaining.empty())
    145 //     {
    146 //       int i = ei_random<int>(0,remaining.size()-1);
    147 //       w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y());
    148 //       remaining[i] = remaining.back();
    149 //       remaining.pop_back();
    150 //     }
    151 //   }
    152 //   VERIFY_IS_APPROX(m, refMat);
    153 
    154     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) ));
    155     #ifdef EIGEN_UNORDERED_MAP_SUPPORT
    156     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) ));
    157     #endif
    158     #ifdef _DENSE_HASH_MAP_H_
    159     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) ));
    160     #endif
    161     #ifdef _SPARSE_HASH_MAP_H_
    162     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) ));
    163     #endif
    164 
    165     // test fillrand
    166     {
    167       DenseMatrix m1(rows,cols);
    168       m1.setZero();
    169       SparseMatrixType m2(rows,cols);
    170       m2.startFill();
    171       for (int j=0; j<cols; ++j)
    172       {
    173         for (int k=0; k<rows/2; ++k)
    174         {
    175           int i = ei_random<int>(0,rows-1);
    176           if (m1.coeff(i,j)==Scalar(0))
    177             m2.fillrand(i,j) = m1(i,j) = ei_random<Scalar>();
    178         }
    179       }
    180       m2.endFill();
    181       VERIFY_IS_APPROX(m2,m1);
    182     }
    183 
    184   // test RandomSetter
    185   /*{
    186     SparseMatrixType m1(rows,cols), m2(rows,cols);
    187     DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
    188     initSparse<Scalar>(density, refM1, m1);
    189     {
    190       Eigen::RandomSetter<SparseMatrixType > setter(m2);
    191       for (int j=0; j<m1.outerSize(); ++j)
    192         for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i)
    193           setter(i.index(), j) = i.value();
    194     }
    195     VERIFY_IS_APPROX(m1, m2);
    196   }*/
    197 //   std::cerr << m.transpose() << "\n\n"  << refMat.transpose() << "\n\n";
    198 //   VERIFY_IS_APPROX(m, refMat);
    199 
    200   // test basic computations
    201   {
    202     DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
    203     DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
    204     DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
    205     DenseMatrix refM4 = DenseMatrix::Zero(rows, rows);
    206     SparseMatrixType m1(rows, rows);
    207     SparseMatrixType m2(rows, rows);
    208     SparseMatrixType m3(rows, rows);
    209     SparseMatrixType m4(rows, rows);
    210     initSparse<Scalar>(density, refM1, m1);
    211     initSparse<Scalar>(density, refM2, m2);
    212     initSparse<Scalar>(density, refM3, m3);
    213     initSparse<Scalar>(density, refM4, m4);
    214 
    215     VERIFY_IS_APPROX(m1+m2, refM1+refM2);
    216     VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
    217     VERIFY_IS_APPROX(m3.cwise()*(m1+m2), refM3.cwise()*(refM1+refM2));
    218     VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
    219 
    220     VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
    221     VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
    222 
    223     VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
    224     VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
    225 
    226     VERIFY_IS_APPROX(m1.col(0).eigen2_dot(refM2.row(0)), refM1.col(0).eigen2_dot(refM2.row(0)));
    227 
    228     refM4.setRandom();
    229     // sparse cwise* dense
    230     VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4);
    231 //     VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
    232   }
    233 
    234   // test innerVector()
    235   {
    236     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
    237     SparseMatrixType m2(rows, rows);
    238     initSparse<Scalar>(density, refMat2, m2);
    239     int j0 = ei_random(0,rows-1);
    240     int j1 = ei_random(0,rows-1);
    241     VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
    242     VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
    243     //m2.innerVector(j0) = 2*m2.innerVector(j1);
    244     //refMat2.col(j0) = 2*refMat2.col(j1);
    245     //VERIFY_IS_APPROX(m2, refMat2);
    246   }
    247 
    248   // test innerVectors()
    249   {
    250     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
    251     SparseMatrixType m2(rows, rows);
    252     initSparse<Scalar>(density, refMat2, m2);
    253     int j0 = ei_random(0,rows-2);
    254     int j1 = ei_random(0,rows-2);
    255     int n0 = ei_random<int>(1,rows-std::max(j0,j1));
    256     VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
    257     VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
    258                      refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
    259     //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
    260     //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0);
    261   }
    262 
    263   // test transpose
    264   {
    265     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
    266     SparseMatrixType m2(rows, rows);
    267     initSparse<Scalar>(density, refMat2, m2);
    268     VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
    269     VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
    270   }
    271 
    272   // test prune
    273   {
    274     SparseMatrixType m2(rows, rows);
    275     DenseMatrix refM2(rows, rows);
    276     refM2.setZero();
    277     int countFalseNonZero = 0;
    278     int countTrueNonZero = 0;
    279     m2.startFill();
    280     for (int j=0; j<m2.outerSize(); ++j)
    281       for (int i=0; i<m2.innerSize(); ++i)
    282       {
    283         float x = ei_random<float>(0,1);
    284         if (x<0.1)
    285         {
    286           // do nothing
    287         }
    288         else if (x<0.5)
    289         {
    290           countFalseNonZero++;
    291           m2.fill(i,j) = Scalar(0);
    292         }
    293         else
    294         {
    295           countTrueNonZero++;
    296           m2.fill(i,j) = refM2(i,j) = Scalar(1);
    297         }
    298       }
    299     m2.endFill();
    300     VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
    301     VERIFY_IS_APPROX(m2, refM2);
    302     m2.prune(1);
    303     VERIFY(countTrueNonZero==m2.nonZeros());
    304     VERIFY_IS_APPROX(m2, refM2);
    305   }
    306 }
    307 
    308 void test_eigen2_sparse_basic()
    309 {
    310   for(int i = 0; i < g_repeat; i++) {
    311     CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(8, 8)) );
    312     CALL_SUBTEST_2( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) );
    313     CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(33, 33)) );
    314 
    315     CALL_SUBTEST_3( sparse_basic(DynamicSparseMatrix<double>(8, 8)) );
    316   }
    317 }
    318