Math/test/Similarity_t.cpp

0001 #define SMATRIX_USE_CONSTEXPR
0002 #include "Math/SMatrix.h"
0003 #include <cmath>
0004 #include <cstddef>
0005 #include <random>
0006
0007 typedef unsigned int IndexType;
0008 //typedef unsigned long long IndexType;
0009
0010 namespace ROOT {
0011
0012   namespace Math {
0013
0014     /**
0015        Force Expression evaluation from general to symmetric.
0016        To be used when is known (like in similarity products) that the result
0017        is symmetric
0018        Note this is function used in the simmilarity product: no check for temporary is
0019        done since in that case is not needed
0020    */
0021     struct AssignAsSym {
0022       /// assign a symmetric matrix from an expression
0023       template <class T, IndexType D, class A, class R>
0024       static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const Expr<A, T, D, D, R>& rhs) {
0025         // for(unsigned int i=0; i<D*D; ++i) lhs.fRep[i] = rhs.apply(i);
0026
0027         // in principle this will do only half+D evaluations (and is ~4 times faster than the above for a multiplication)
0028         for (IndexType i = 0; i < D; ++i)
0029           // storage of symmetric matrix is in lower block
0030           for (IndexType j = 0; j <= i; ++j)
0031             lhs(i, j) = rhs(i, j);
0032         // symmetrize
0033         for (IndexType i = 0; i != D - 1; ++i)
0034           for (IndexType j = i + 1; j != D; ++j)
0035             lhs(i, j) = lhs(j, i);
0036       }
0037       /// assign the "symmetric" matric  from a general matrix
0038       template <class T, IndexType D, class R>
0039       static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const SMatrix<T, D, D, R>& rhs) {
0040         for (IndexType i = 0; i != D * D; ++i)
0041           lhs.fRep[i] = rhs.apply(i);
0042         /*
0043        // useful only if we do not store the upper triangle
0044        for( IndexType i=0; i<D; ++i)
0045        // storage of symmetric matrix is in lower block
0046        for( IndexType j=0; j<=i; ++j)
0047        lhs(i,j) = rhs(i,j);
0048        for (IndexType i=0; i!=D-1; ++i)
0049      for (IndexType j=i+1; j!=D; ++j)
0050        lhs(i,j) = lhs(j,i);
0051        */
0052       }
0053     };  // struct AssignAsSym
0054
0055     template <class T, IndexType D1, IndexType D2, class R>
0056     inline SMatrix<T, D1, D1, MatRepSym<T, D1> > Similarity1(const SMatrix<T, D1, D2, R>& lhs,
0057                                                              const SMatrix<T, D2, D2, MatRepStd<T, D2> >& rhs) {
0058       SMatrix<T, D1, D2, MatRepStd<T, D1, D2> > tmp = lhs * rhs;
0059       typedef SMatrix<T, D1, D1, MatRepSym<T, D1> > SMatrixSym;
0060       SMatrixSym mret{SMatrixNoInit{}};
0061       AssignSym::Evaluate(mret, tmp * Transpose(lhs));
0062       return mret;
0063     }
0064
0065     template <class T, IndexType D1, IndexType D2>
0066     inline SMatrix<T, D1, D1, MatRepStd<T, D1> > Similarity2(const SMatrix<T, D1, D2, MatRepStd<T, D1, D2> >& lhs,
0067                                                              const SMatrix<T, D2, D2, MatRepStd<T, D2> >& rhs) {
0068       SMatrix<T, D1, D2, MatRepStd<T, D1, D2> > tmp = lhs * rhs;
0069       typedef SMatrix<T, D1, D1, MatRepStd<T, D1> > SMatrixSym;
0070       SMatrixSym mret = SMatrixNoInit();
0071       AssignAsSym::Evaluate(mret, tmp * Transpose(lhs));
0072       return mret;
0073     }
0074   }  // namespace Math
0075 }  // namespace ROOT
0076
0077 // U(i,k) * A(k,l) * U(j,l)
0078 template <typename T, IndexType D1, IndexType D2>
0079 inline void similarity(ROOT::Math::SMatrix<T, D1, D1, ROOT::Math::MatRepStd<T, D1> >& b,
0080                        ROOT::Math::SMatrix<T, D1, D2, ROOT::Math::MatRepStd<T, D1, D2> > const& u,
0081                        ROOT::Math::SMatrix<T, D2, D2, ROOT::Math::MatRepStd<T, D2> > const& a) {
0082   // brute force loop
0083   for (IndexType i = 0; i != D1; ++i)
0084     for (IndexType j = 0; j <= i; ++j)
0085       for (IndexType k = 0; k != D2; ++k)
0086         for (IndexType l = 0; l != D2; ++l)
0087           b(i, j) += u(i, k) * a(k, l) * u(j, l);
0088
0089   /*
0090   for (IndexType is=0; is<D1; is+=4) {
0091     IndexType ie=std::min(is+4,D1);
0092     //for (IndexType js=0; js<=is; js+=4) {
0093       for (IndexType ks=0; ks<D2; ks+=4) {
0094     IndexType ke=std::min(ks+4,D2);
0095     for (IndexType i=is; i!=ie; ++i)
0096       for (IndexType j=0; j<=i; ++j)
0097         //    for (IndexType j=js; j<=std::min(js+3,i); ++j)
0098         for (IndexType k=ks; k!=ke; ++k)
0099           for (IndexType l=0; l!=D2; ++l)
0100         b(i,j) += u(i,k)*a(k,l)*u(j,l);
0101       }
0102       //}
0103   }
0104   */
0105
0106   //    for (IndexType l=0; l<=k; ++l)
0107   //  b(i,j) += (u(i,k)*u(j,l)+u(i,l)*u(j,k))*a(k,l);
0108   /*
0109   T s[N];
0110   for (IndexType i=0; i!=N; ++i)
0111     for (IndexType j=0; j<=i; ++j) {
0112       for (IndexType k=0; k!=N; ++k) {
0113     s[k]=0;
0114     for (IndexType l=0; l!=N; ++l)
0115       s[k] += a(k,l)*u(j,l);
0116       }
0117       for (IndexType k=0; k!=N; ++k)
0118     b(i,j) += u(i,k)*s[k];
0119     }
0120   */
0121
0122   for (IndexType i = 0; i != D1 - 1; ++i)
0123     for (IndexType j = i + 1; j != D1; ++j)
0124       b(i, j) = b(j, i);
0125 }
0126
0127 template <typename M1, typename M2>
0128 double eps(M1 const& m1, M2 const& m2) {
0129   IndexType N = M1::kRows;
0130   double ret = 0.;
0131   for (IndexType i = 0; i != N; ++i)
0132     for (IndexType j = 0; j != N; ++j)
0133       ret = std::max(ret, std::abs(m1(i, j) - m2(i, j)));
0134   return ret;
0135 }
0136
0137 template <typename M>
0138 bool isSym(M const& m) {
0139   IndexType N = M::kRows;
0140   for (IndexType i = 0; i != N; ++i)
0141     for (IndexType j = 0; j <= i; ++j)
0142       if (m(i, j) != m(j, i))
0143         return false;
0144   return true;
0145 }
0146
0147 #include <iostream>
0148 #include "FWCore/Utilities/interface/HRRealTime.h"
0149
0150 namespace {
0151   std::mt19937 eng;
0152   std::uniform_real_distribution<double> rgen(-5., 5.);
0153
0154   inline double rr() { return rgen(eng); }
0155
0156   template <typename T, IndexType D1, IndexType D2>
0157   inline void fillRandom(ROOT::Math::SMatrix<T, D1, D2, ROOT::Math::MatRepStd<T, D1, D2> >& a) {
0158     for (IndexType i = 0; i != D1; ++i) {
0159       for (IndexType j = 0; j != D2; ++j)
0160         a(i, j) = rr();
0161     }
0162   }
0163   template <typename T, IndexType N>
0164   inline void fillRandomSym(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& a) {
0165     for (IndexType i = 0; i != N; ++i) {
0166       for (IndexType j = 0; j <= i; ++j)
0167         a(i, j) = rr();
0168     }
0169     for (IndexType i = 0; i != N - 1; ++i)
0170       for (IndexType j = i + 1; j != N; ++j)
0171         a(i, j) = a(j, i);
0172   }
0173 }  // namespace
0174
0175 bool ok = true;
0176
0177 template <typename T, IndexType D1, IndexType D2>
0178 void go(edm::HRTimeType& s1, edm::HRTimeType& s2, edm::HRTimeType& s3, edm::HRTimeType& s4, bool print) {
0179   typedef ROOT::Math::SMatrix<T, D1, D2, ROOT::Math::MatRepStd<T, D1, D2> > JMatrix;
0180   typedef ROOT::Math::SMatrix<T, D1, D1, ROOT::Math::MatRepStd<T, D1, D1> > Matrix1;
0181   typedef ROOT::Math::SMatrix<T, D1, D1, ROOT::Math::MatRepSym<T, D1> > SymMatrix1;
0182   typedef ROOT::Math::SMatrix<T, D2, D2, ROOT::Math::MatRepStd<T, D2, D2> > Matrix2;
0183   typedef ROOT::Math::SMatrix<T, D2, D2, ROOT::Math::MatRepSym<T, D2> > SymMatrix2;
0184
0185   JMatrix lh;
0186   fillRandom(lh);
0187   Matrix2 rh;
0188   fillRandomSym(rh);
0189
0190   SymMatrix2 srh;
0191   ROOT::Math::AssignSym::Evaluate(srh, rh);
0192
0193   SymMatrix1 res1;
0194   s1 = edm::hrRealTime();
0195   res1 = ROOT::Math::Similarity(lh, srh);
0196   s1 = edm::hrRealTime() - s1;
0197
0198   SymMatrix1 res2;
0199   s2 = edm::hrRealTime();
0200   res2 = ROOT::Math::Similarity1(lh, rh);
0201   s2 = edm::hrRealTime() - s2;
0202
0203   Matrix1 res3;
0204   s3 = edm::hrRealTime();
0205   res3 = ROOT::Math::Similarity2(lh, rh);
0206   s3 = edm::hrRealTime() - s3;
0207
0208   Matrix1 res4;
0209   s4 = edm::hrRealTime();
0210   similarity(res4, lh, rh);
0211   s4 = edm::hrRealTime() - s4;
0212
0213   if (print) {
0214     if (!isSym(rh))
0215       std::cout << " rh is not sym" << std::endl;
0216     if (!isSym(res1))
0217       std::cout << " res1 is not sym" << std::endl;
0218     if (!isSym(res2))
0219       std::cout << " res2 is not sym" << std::endl;
0220     if (!isSym(res3))
0221       std::cout << " res3 is not sym" << std::endl;
0222     if (!isSym(res4))
0223       std::cout << " res4 is not sym" << std::endl;
0224
0225     std::cout << D1 << "x" << D2 << std::endl;
0226     std::cout << "eps sim  " << eps(res1, res2) << std::endl;
0227     std::cout << "eps std  " << eps(res1, res3) << std::endl;
0228     std::cout << "eps loop  " << eps(res1, res4) << std::endl;
0229   }
0230
0231   ok &= isSym(rh) && isSym(res1) && isSym(res2) && isSym(res3) && isSym(res4);
0232 }
0233
0234 template <typename T, IndexType D1, IndexType D2>
0235 void loop(std::ostream& co) {
0236   ok = true;
0237
0238   int N = 100000;
0239   edm::HRTimeType t1 = 0;
0240   edm::HRTimeType t2 = 0;
0241   edm::HRTimeType t3 = 0;
0242   edm::HRTimeType t4 = 0;
0243   edm::HRTimeType s1 = 0, s2 = 0, s3 = 0, s4 = 0;
0244   for (int i = 0; i != N; ++i) {
0245     go<T, D1, D2>(s1, s2, s3, s4, false);
0246     t1 += s1;
0247     t2 += s2;
0248     t3 += s3;
0249     t4 += s4;
0250   }
0251
0252   std::cout << D1 << "x" << D2 << std::endl;
0253   std::cout << "root sim " << t1 / N << std::endl;
0254   std::cout << "sym  sim " << t2 / N << std::endl;
0255   std::cout << "std  sim " << t3 / N << std::endl;
0256   std::cout << "loop sim " << t4 / N << std::endl;
0257
0258   co << "|  " << t1 / N << "|  " << t2 / N << "|  " << t3 / N << "|  " << t4 / N << "|";
0259
0260   if (ok)
0261     std::cout << " OK " << std::endl;
0262 }
0263
0264 #include <sstream>
0265 int main() {
0266   edm::HRTimeType s1 = 0, s2 = 0, s3 = 0, s4 = 0;
0267
0268   go<double, 3, 3>(s1, s2, s3, s4, true);
0269   go<double, 5, 5>(s1, s2, s3, s4, true);
0270   go<double, 5, 15>(s1, s2, s3, s4, true);
0271   go<double, 15, 15>(s1, s2, s3, s4, true);
0272   std::cout << std::endl;
0273
0274   std::ostringstream co;
0275   co << "|  *3x3*  |||| |  *5x5*  |||| |  *5x15*  |||| |  *15x15*  ||||\n";
0276   co << "|  *root*  |  *sym*  |  *std*  |  *loop*  |";
0277   co << "|  *root*  |  *sym*  |  *std*  |  *loop*  |";
0278   co << "|  *root*  |  *sym*  |  *std*  |  *loop*  |";
0279   co << "|  *root*  |  *sym*  |  *std*  |  *loop*  |\n";
0280   loop<double, 3, 3>(co);
0281   loop<double, 5, 5>(co);
0282   loop<double, 5, 15>(co);
0283   loop<double, 15, 15>(co);
0284   std::cout << co.str() << std::endl;
0285
0286   return 0;
0287 }