GsfTools/test/KullbackLeiblerDistance_t.cpp

0001 #include "TrackingTools/GsfTools/interface/KullbackLeiblerDistance.h"
0002 #include "TrackingTools/GsfTools/interface/DistanceBetweenComponents.h"
0003
0004 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
0005 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0006 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0007
0008 #include "DataFormats/GeometrySurface/interface/Surface.h"
0009 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0010
0011 #include "DataFormats/GeometrySurface/interface/Plane.h"
0012
0013 #include "FWCore/Utilities/interface/HRRealTime.h"
0014 #include <iostream>
0015 #include <vector>
0016
0017 bool isAligned(const void* data, long alignment) {
0018   // check that the alignment is a power of two
0019   assert((alignment & (alignment - 1)) == 0);
0020   return ((long)data & (alignment - 1)) == 0;
0021 }
0022
0023 void st() {}
0024 void en() {}
0025
0026 typedef DistanceBetweenComponents<5> Distance;
0027 typedef KullbackLeiblerDistance<5> KDistance;
0028 typedef SingleGaussianState<5> GS;
0029 typedef GS::Vector Vector;
0030 typedef GS::Matrix Matrix;
0031 typedef ROOT::Math::SMatrix<double, 6, 6, ROOT::Math::MatRepSym<double, 6> > Matrix6;
0032
0033 Distance const& distance() {
0034   static Distance* d = new KDistance;
0035   return *d;
0036 }
0037
0038 Matrix buildCovariance(float y) {
0039   // build a resonable covariance matrix as JIJ
0040
0041   Basic3DVector<float> axis(0.5, 1., 1);
0042
0043   Surface::RotationType rot(axis, 0.5 * M_PI);
0044
0045   Surface::PositionType pos(0., 0., 0.);
0046
0047   Plane plane(pos, rot);
0048   LocalTrajectoryParameters tp(1., 1., y, 0., 0., 1.);
0049
0050   JacobianLocalToCartesian jl2c(plane, tp);
0051   return ROOT::Math::SimilarityT(jl2c.jacobian(), Matrix6(ROOT::Math::SMatrixIdentity()));
0052   // return  ROOT::Math::Transpose(jl2c.jacobian())* jl2c.jacobian();
0053 }
0054
0055 int main(int argc, char* argv[]) {
0056   std::cout << "size of  SingleGaussianState<5>" << sizeof(GS) << std::endl;
0057   std::cout << "size of  Matrix" << sizeof(Matrix) << std::endl;
0058   std::cout << "size of  Vector" << sizeof(Vector) << std::endl;
0059
0060   Distance const& d = distance();
0061
0062   Matrix cov1 = buildCovariance(1.);
0063   Matrix cov2 = buildCovariance(2.);
0064
0065   GS* gs1 = new GS(Vector(1., 1., 1., 1., 1.), cov1);
0066   // GS gs1(Vector(1., 1.,1., 1.,1.),Matrix(ROOT::Math::SMatrixIdentity()));
0067
0068   GS* gs0 = new GS(Vector(1., 1., 1., 0., 0.), Matrix(ROOT::Math::SMatrixIdentity()));
0069   GS* gsP = new GS(Vector(1., 1., 1., 10., 10.), Matrix(ROOT::Math::SMatrixIdentity()));
0070
0071   std::cout << "GS " << ((isAligned(gs0, 16)) ? "a " : "n ") << std::endl;
0072   std::cout << "cov " << ((isAligned(&gs0->covariance(), 16)) ? "a " : "n ") << std::endl;
0073   std::cout << "mean " << ((isAligned(&gs0->mean(), 16)) ? "a " : "n ") << std::endl;
0074   std::cout << "weightM " << ((isAligned(&gs0->weightMatrix(), 16)) ? "a " : "n ") << std::endl;
0075
0076   // GS gs2(Vector(2., 2., 2., 2.,2.),cov);
0077   GS* gs2 = new GS(Vector(2., 2., 2., 2., 2.), cov2);
0078
0079   std::vector<GS> vgs(10000);
0080   vgs.front() = *gs1;
0081   vgs.back() = *gs2;
0082
0083   // make sure we load all code...
0084   edm::HRTimeType s0 = edm::hrRealTime();
0085   double res = d(*gs0, *gsP);
0086   edm::HRTimeType e0 = edm::hrRealTime();
0087   std::cout << e0 - s0 << std::endl;
0088
0089   double res2 = 0;
0090
0091   bool both = (argc < 2);
0092
0093   if (both || argv[1][0] == 'a') {
0094     st();
0095     edm::HRTimeType s = edm::hrRealTime();
0096     res2 = d(*gs1, *gs2);
0097     edm::HRTimeType e = edm::hrRealTime();
0098     en();
0099     std::cout << e - s << std::endl;
0100   }
0101   if (both || argv[1][0] == 'b') {
0102     st();
0103     edm::HRTimeType s = edm::hrRealTime();
0104     res2 = d(vgs.front(), vgs.back());
0105     edm::HRTimeType e = edm::hrRealTime();
0106     en();
0107     std::cout << e - s << std::endl;
0108   }
0109
0110   std::cout << res << " " << res2 << std::endl;
0111
0112   res = 0;
0113   st();
0114   edm::HRTimeType s = edm::hrRealTime();
0115   for (int i = 0; i < 100; ++i)
0116     for (auto const& s : vgs)
0117       res += d(vgs.front(), s);
0118   edm::HRTimeType e = edm::hrRealTime();
0119   en();
0120
0121   std::cout << e - s << std::endl;
0122
0123   std::cout << res << std::endl;
0124
0125   return 0;
0126 }