Skip to content

Extend templation in GenVector to better support vectorised Vc types. #394

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vgvassilev merged 54 commits into from
Mar 15, 2017
Merged

Extend templation in GenVector to better support vectorised Vc types. #394

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
54 commits
Select commit Hold shift + click to select a range
76a8f9d
Template the GenVector Plane3D, Transform3D and Translation3D classes…
cjones051073 Feb 25, 2017
7e9c2b7
Use the 'using namespace std' trick with Cartesian3D and Plane3D to f…
cjones051073 Feb 25, 2017
b15d62a
Extend the 'using namespace std' trick to a few more classes
cjones051073 Feb 25, 2017
2b2fc17
Change class to typename
cjones051073 Feb 25, 2017
2e74f0d
Explicitly cast some constants to the scalar type (needed for Vc types)
cjones051073 Feb 25, 2017
c3175c6
Remove empty implementation files
cjones051073 Feb 25, 2017
9bdfb8e
Template the Boost classes
cjones051073 Feb 25, 2017
c97245c
Remove empty implementation files for Boost classes
cjones051073 Feb 25, 2017
f12d368
Clean up some indentations
cjones051073 Feb 25, 2017
bbbbbd0
Add some missing template types to the Boost classes
cjones051073 Feb 25, 2017
9a0ee77
Remove empty Boost implementation file
cjones051073 Feb 25, 2017
371d50d
Small cleanups
cjones051073 Feb 25, 2017
dcb0b30
Clean up some initialisations
cjones051073 Feb 25, 2017
29e7e30
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
14e73b4
Add Scalar/Vector Plane3D::Normalise() methods using SFINAE
cjones051073 Feb 27, 2017
27bbb57
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
196a00b
Roll back some changes for now, as they are causing problems in the t…
cjones051073 Feb 27, 2017
0ad5828
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
1928c75
Add specific vectorised implementations of a few methods
cjones051073 Feb 28, 2017
5ccdaea
Extend ostream operator for Displacement and Position vectors to supp…
cjones051073 Feb 28, 2017
ef1aea9
Fix the vector abs method (not fabs) and Vector Unit() for vectorised…
cjones051073 Feb 28, 2017
6cac16a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 28, 2017
68ba28a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 1, 2017
939b7cb
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 4, 2017
d9b810b
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 6, 2017
776ce1a
Add override qualifier
cjones051073 Mar 6, 2017
fbe74d2
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 10, 2017
12b43b0
Fix initialisation bug
cjones051073 Mar 10, 2017
8367bdf
remove space
cjones051073 Mar 10, 2017
0f404ba
Some small changes following review comments
cjones051073 Mar 10, 2017
65a701d
remove spaces
cjones051073 Mar 10, 2017
4ba8a1b
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 10, 2017
ff94a5a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 11, 2017
b8ca822
Fix formatting using clang-format
cjones051073 Mar 11, 2017
b32cb6a
some more 'using namespace std' updates
cjones051073 Mar 11, 2017
febe783
Fix more clang format issues
cjones051073 Mar 11, 2017
a628a2a
fix (hopefully) one more clang format issue
cjones051073 Mar 11, 2017
89b2d26
clang format is picky...
cjones051073 Mar 11, 2017
6bf0368
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 11, 2017
bd96bd0
revert changes in VectorUtils.h
cjones051073 Mar 13, 2017
dea8d05
Add new test for GenVector with Vc types
cjones051073 Mar 13, 2017
a18221a
extend test to 96 photons
cjones051073 Mar 13, 2017
eb0b8ed
reorder headers
cjones051073 Mar 13, 2017
7da1945
remove comment that number of test photons must be multiple of 8, as …
cjones051073 Mar 13, 2017
a271eab
Add <cmath> header includes
cjones051073 Mar 13, 2017
decbe11
revert back to using std:: as it seems to cause issues with gcc build…
cjones051073 Mar 14, 2017
58c6315
Include Vc before ROOT in test application. Needed for std:: support...
cjones051073 Mar 14, 2017
17f42a3
Add some more missing std::
cjones051073 Mar 14, 2017
434396d
More math updates
cjones051073 Mar 14, 2017
b145899
Change enable_if implementation in GenVector
cjones051073 Mar 14, 2017
ad8d3a5
update comment
cjones051073 Mar 14, 2017
b018a03
Extend the GenVector Vc test with a timing test that asserts that the…
cjones051073 Mar 14, 2017
bb1e386
add fabs on timing diff and extend Vc genvector expected diff to 25%
cjones051073 Mar 14, 2017
1786383
do not return a failure from the Vc speed test, as test could be run …
cjones051073 Mar 14, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Extend the 'using namespace std' trick to a few more classes
  • Loading branch information
@cjones051073
cjones051073 committed Feb 25, 2017
commit b15d62aebad03d5c5d4e4219c919a9ca29ebff0d
24 changes: 14 additions & 10 deletions math/genvector/inc/Math/GenVector/Cartesian2D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,8 +89,9 @@ public :
Scalar X() const { return fX;}
Scalar Y() const { return fY;}
Scalar Mag2() const { return fX*fX + fY*fY; }
Scalar R() const { return std::sqrt( Mag2());}
Scalar Phi() const { return (fX==0 && fY==0) ? 0.0 : atan2(fY,fX);}
Scalar R() const { using namespace std; return sqrt( Mag2()); }
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm, what is the advantage of the using namespace std; pattern, instead of using return std::sqrt...?

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is required to support Vc types. For Vc types you need to use the versions of sqrt etc. shipped as part of that library. Using std::sqrt does not allow this. The namespace trick means you use the std:: versions when appropriate, but also allows other implementations when not.

Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Huh, thanks for explaining! Neat and subtle at the same time. Btw, there is an inherited extra space after the open paren.

Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw, it's @lmoneta's call but I'd prefer to have a short comment saying what the pattern does, or maybe a #define R__enable_vc_types using namespace std; This would make this easier to read by an unarmed eye.

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dislike preprocessor directives unless absolutely required, and I don't think here adding one really helps. To my eye the pattern as is is clear, but as you say its up to @lmoneta . I would be OK with adding some sort of comment somewhere ?

Copy link
Member

@pcanal pcanal Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a comment would be as clear and more concise than a ifdef

Copy link
Member

@lmoneta lmoneta Mar 13, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi,

I don't understand this. Are you sure it is needed for Vc ? In the past I have seen Vc replacing automatically its vectorised function implementations in the std namespace

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 13, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am sure I ran into some issue that lead me to doing this... Also note the templation is in principle there for other types than Vc. This change in principle supports types that do not extend the math functions under std::, as well as those that do.

Scalar Phi() const { using namespace std;
return (fX==Scalar(0) && fY==Scalar(0)) ? Scalar(0) : atan2(fY,fX);}

/**
set the x coordinate value keeping y constant
Expand Down Expand Up @@ -124,8 +125,9 @@ public :
rotate by an angle
*/
void Rotate(Scalar angle) {
Scalar s = std::sin(angle);
Scalar c = std::cos(angle);
using namespace std;
const Scalar s = sin(angle);
const Scalar c = cos(angle);
SetCoordinates( c*fX - s*fY, s*fX + c * fY );
}

Expand Down Expand Up @@ -161,10 +163,11 @@ public :
template <class T2>
explicit Cartesian2D( const Polar2D<T2> & v )
{
Scalar r = v.R(); // re-using this instead of calling v.X() and v.Y()
using namespace std;
const Scalar r = v.R(); // re-using this instead of calling v.X() and v.Y()
// is the speed improvement
fX = r * std::cos(v.Phi());
fY = r * std::sin(v.Phi());
fX = r * cos(v.Phi());
fY = r * sin(v.Phi());
}
// Technical note: This works even though only Polar2Dfwd.h is
// included (and in fact, including Polar2D.h would cause circularity
Expand All @@ -174,9 +177,10 @@ public :
template <class T2>
Cartesian2D & operator = (const Polar2D<T2> & v)
{
Scalar r = v.R();
fX = r * std::cos(v.Phi());
fY = r * std::sin(v.Phi());
using namespace std;
const Scalar r = v.R();
fX = r * cos(v.Phi());
fY = r * sin(v.Phi());
return *this;
}

Expand Down
2 changes: 1 addition & 1 deletion math/genvector/inc/Math/GenVector/Cartesian3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ public :

// pseudorapidity
Scalar Eta() const {
return Impl::Eta_FromRhoZ ( Rho(),fZ);
return Impl::Eta_FromRhoZ( Rho(), fZ );
}

/**
Expand Down
19 changes: 12 additions & 7 deletions math/genvector/inc/Math/GenVector/Cylindrical3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,10 +106,12 @@ public :
{rho=fRho; zz=fZ; phi=fPhi;}

private:
inline static Scalar pi() { return M_PI; }
inline void Restrict() {
inline static Scalar pi() { return Scalar(M_PI); }
inline void Restrict()
{
using namespace std;
if ( fPhi <= -pi() || fPhi > pi() )
fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
fPhi = fPhi - floor( fPhi/(2*pi()) +.5 ) * 2*pi();
return;
}
public:
Expand All @@ -120,13 +122,16 @@ public :
Scalar Z() const { return fZ; }
Scalar Phi() const { return fPhi; }

Scalar X() const { return fRho*std::cos(fPhi); }
Scalar Y() const { return fRho*std::sin(fPhi); }
Scalar X() const { using namespace std; return fRho*cos(fPhi); }
Scalar Y() const { using namespace std; return fRho*sin(fPhi); }

Scalar Mag2() const { return fRho*fRho + fZ*fZ; }
Scalar R() const { return std::sqrt( Mag2()); }
Scalar R() const { using namespace std; return sqrt( Mag2()); }
Scalar Perp2() const { return fRho*fRho; }
Scalar Theta() const { return (fRho==0 && fZ==0 ) ? 0 : atan2(fRho,fZ); }
Scalar Theta() const {
using namespace std;
return (fRho==Scalar(0) && fZ==Scalar(0) ) ? Scalar(0) : atan2(fRho,fZ);
}

// pseudorapidity - use same implementation as in Cartesian3D
Scalar Eta() const {
Expand Down
7 changes: 4 additions & 3 deletions math/genvector/inc/Math/GenVector/Polar2D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,8 +97,8 @@ public :

Scalar R() const { return fR;}
Scalar Phi() const { return fPhi; }
Scalar X() const { return fR*std::cos(fPhi);}
Scalar Y() const { return fR*std::sin(fPhi);}
Scalar X() const { using namespace std; return fR*cos(fPhi);}
Scalar Y() const { using namespace std; return fR*sin(fPhi);}
Scalar Mag2() const { return fR*fR;}


Expand Down Expand Up @@ -134,8 +134,9 @@ public :
restrict abgle hi to be between -PI and PI
*/
inline void Restrict() {
using namespace std;
if ( fPhi <= -pi() || fPhi > pi() )
fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
fPhi = fPhi - floor( fPhi/(2*pi()) +.5 ) * 2*pi();
return;
}

Expand Down
44 changes: 25 additions & 19 deletions math/genvector/inc/Math/GenVector/PtEtaPhiE4D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,26 +137,28 @@ public :

// other coordinate representation

Scalar Px() const { return fPt*cos(fPhi);}
Scalar Px() const { using namespace std; return cos(fPhi);}
Scalar X () const { return Px(); }
Scalar Py() const { return fPt*sin(fPhi);}
Scalar Py() const { using namespace std; return sin(fPhi);}
Scalar Y () const { return Py(); }
Scalar Pz() const {
return fPt > 0 ? fPt*std::sinh(fEta) :
fEta == 0 ? 0 :
fEta > 0 ? fEta - etaMax<Scalar>() :
fEta + etaMax<Scalar>() ;
using namespace std;
return ( fPt > 0 ? fPt*sinh(fEta) :
fEta == 0 ? 0 :
fEta > 0 ? fEta - etaMax<Scalar>() :
fEta + etaMax<Scalar>() );
}
Scalar Z () const { return Pz(); }

/**
magnitude of momentum
*/
Scalar P() const {
return fPt > 0 ? fPt*std::cosh(fEta) :
fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>() :
fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() :
0 ;
using namespace std;
return ( fPt > 0 ? fPt*cosh(fEta) :
fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>() :
fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() :
0 );
}
Scalar R() const { return P(); }

Expand All @@ -175,13 +177,14 @@ public :
invariant mass
*/
Scalar M() const {
using namespace std;
Scalar mm = M2();
if (mm >= 0) {
return std::sqrt(mm);
return sqrt(mm);
} else {
GenVector::Throw ("PtEtaPhiE4D::M() - Tachyonic:\n"
" Pt and Eta give P such that P^2 > E^2, so the mass would be imaginary");
return -std::sqrt(-mm);
return -sqrt(-mm);
}
}
Scalar Mag() const { return M(); }
Expand All @@ -201,13 +204,14 @@ public :
transverse mass
*/
Scalar Mt() const {
using namespace std;
Scalar mm = Mt2();
if (mm >= 0) {
return std::sqrt(mm);
return sqrt(mm);
} else {
GenVector::Throw ("PtEtaPhiE4D::Mt() - Tachyonic:\n"
" Pt and Eta give Pz such that Pz^2 > E^2, so the mass would be imaginary");
return -std::sqrt(-mm);
return -sqrt(-mm);
}
}

Expand All @@ -218,7 +222,8 @@ public :
transverse energy
*/
Scalar Et() const {
return fE / std::cosh(fEta); // faster using eta
using namespace std;
return fE / cosh(fEta); // faster using eta
}

/**
Expand All @@ -230,8 +235,9 @@ public :
private:
inline static Scalar pi() { return M_PI; }
inline void Restrict() {
using namespace std;
if ( fPhi <= -pi() || fPhi > pi() )
fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
fPhi = fPhi - floor( fPhi/(2*pi()) +.5 ) * 2*pi();
return;
}
public:
Expand All @@ -240,9 +246,9 @@ public :
polar angle
*/
Scalar Theta() const {
if (fPt > 0) return 2* std::atan( exp( - fEta ) );
if (fEta >= 0) return 0;
return pi();
using namespace std;
return ( fPt > 0 ? Scalar(2) * atan( exp( - fEta ) ) :
fEta >= 0 ? Scalar(0) : pi() );
}

// --------- Set Coordinates of this system ---------------
Expand Down
23 changes: 14 additions & 9 deletions math/genvector/inc/Math/GenVector/PxPyPzE4D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,13 +148,14 @@ public :
invariant mass
*/
Scalar M() const {
using namespace std;
Scalar mm = M2();
if (mm >= 0) {
return std::sqrt(mm);
return sqrt(mm);
} else {
GenVector::Throw ("PxPyPzE4D::M() - Tachyonic:\n"
" P^2 > E^2 so the mass would be imaginary");
return -std::sqrt(-mm);
return -sqrt(-mm);
}
}
Scalar Mag() const { return M(); }
Expand All @@ -168,7 +169,7 @@ public :
/**
Transverse spatial component (P_perp or rho)
*/
Scalar Pt() const { return std::sqrt(Perp2());}
Scalar Pt() const { using namespace std; return std::sqrt(Perp2());}
Scalar Perp() const { return Pt();}
Scalar Rho() const { return Pt();}

Expand All @@ -181,13 +182,14 @@ public :
transverse mass
*/
Scalar Mt() const {
using namespace std;
Scalar mm = Mt2();
if (mm >= 0) {
return std::sqrt(mm);
if ( mm >= 0 ) {
return sqrt(mm);
} else {
GenVector::Throw ("PxPyPzE4D::Mt() - Tachyonic:\n"
" Pz^2 > E^2 so the transverse mass would be imaginary");
return -std::sqrt(-mm);
return -sqrt(-mm);
}
}

Expand All @@ -204,22 +206,25 @@ public :
transverse energy
*/
Scalar Et() const {
using namespace std;
Scalar etet = Et2();
return fT < 0.0 ? -std::sqrt(etet) : std::sqrt(etet);
return fT < 0.0 ? -sqrt(etet) : sqrt(etet);
}

/**
azimuthal angle
*/
Scalar Phi() const {
return (fX == 0.0 && fY == 0.0) ? 0 : std::atan2(fY,fX);
using namespace std;
return ( fX == Scalar(0) && fY == Scalar(0) ) ? Scalar(0) : atan2(fY,fX);
}

/**
polar angle
*/
Scalar Theta() const {
return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(),fZ);
using namespace std;
return (fX == Scalar(0) && fY == Scalar(0) && fZ == Scalar(0) ) ? Scalar(0) : atan2(Pt(),fZ);
}

/**
Expand Down
7 changes: 4 additions & 3 deletions math/genvector/inc/Math/GenVector/eta.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,16 @@ namespace ROOT {
*/
template<typename Scalar>
inline Scalar Eta_FromRhoZ(Scalar rho, Scalar z) {
using namespace std;
if (rho > 0) {

// value to control Taylor expansion of sqrt
static const Scalar big_z_scaled =
std::pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));
pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));

Scalar z_scaled = z/rho;
if (std::fabs(z_scaled) < big_z_scaled) {
return std::log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
if (fabs(z_scaled) < big_z_scaled) {
return log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
} else {
// apply correction using first order Taylor expansion of sqrt
return z>0 ? std::log(2.0*z_scaled + 0.5/z_scaled) : -std::log(-2.0*z_scaled);
Expand Down