November 2008 - IMP-dev

add support for triplets of particles
by Daniel Russel 17 Nov '08

17 Nov '08

The patch adds a new type ParticleTriplet. This patch is needed for some upcoming changes to the containers. It also, as a side effect to clean up header inclusion order, moves the container macros to internal/container_macros.h (the macros need various headers to be included and so aren't as simple as the rest of the macros in macros.h).

1 0

Trac?
by Keren Lasker 17 Nov '08

17 Nov '08

can we use Trac ( or an alternative solution) to visualize the IMP SVN tracking. http://trac.edgewall.org/ ??

2 3

Fwd: [EMBED-dev] 2D geometry
by Javier Ángel Velázquez Muriel 15 Nov '08

15 Nov '08

---------- Forwarded message ---------- From: Javier Ángel Velázquez-Muriel <javi.velazquez(a)gmail.com> Date: 2008/11/15 Subject: Re: [EMBED-dev] 2D geometry To: List for EMBED development <embed-dev(a)salilab.org> Cc: List for IMP development <imp-dev(a)salilab.org> This is what I have so far. I expect to work on Matrix2D and Matrix3D (volume) soon as derived classes from MultiArray 2008/11/13 Keren Lasker <kerenl(a)salilab.org> > Frido/Javi > How about … [View More]

1 0

Re: [IMP-dev] [EMBED-dev] 2D geometry
by Keren Lasker 13 Nov '08

13 Nov '08

Frido/Javi How about Monday 10am or 11am SF time ? Javi please send the classes headers before we meet so we can save some unnecessary iterations trying to understand one another - communicating via skype is hard enough ;) On Nov 13, 2008, at 3:07 AM, Friedrich Foerster wrote: > yes, would be a good idea. > javi, if you plan on being compatible with spider/xmipp i urge you to > modify the 3d embed structure. it definitely makes no sense to treat > 2d and 3d images … [View More]differently. typically, all operations (rotations, > translations, fourier, etc) are 2d and 3d in these packages as well. > otherwise, embed, imp, or wherever you put it will end up in a mess > with all kinds of different classes. > if you want to change the 3d embed data structure for being compatible > with other packages you got my vote by all means. > > best > > frido > > > On Nov 13, 2008, at 8:47 AM, Keren Lasker wrote: > >> yes Javi - see my last email. >> The geometry manipulation classes are not related to data storage. >> Anyway - it would be nice to have a discussion between you, Frido and >> me on data storage - would you mind having a three way skype call >> soon? >> On Nov 12, 2008, at 7:43 PM, Javier Ángel Velázquez Muriel wrote: >> >>> Keren, the more I see your Matrix3D class the more realize that my >>> code is going to serve a different purpose. The matrix you define in >>> Matrix3D is a 3x3x3, but what I want is a NxLxM, something like the >>> map that you had for EMBed, and in my case, what I understood by >>> matrix 2D was a NxM, not a 2x2, for example. Vectors here don't >>> apply, >>> but still is going to be very helpful to have a templated vector. >>> Again, I am building the templated array. >>> Javi >>> >>> 2008/11/12 Keren Lasker <kerenl(a)salilab.org>: >>>> Javi, >>>> >>>> For your 2D geometry classes (Vector2D, Matrix2D ...) - why not >>>> just using >>>> the corresponding 3D classes and keeping Z fixed? >>>> I think it is not too bad in performance and will save you a lot of >>>> coding/debugging time. >>>> >>>> Keren. >>>> >> >> _______________________________________________ >> EMBED-dev mailing list >> EMBED-dev(a)salilab.org >> https://salilab.org/mailman/listinfo/embed-dev > > -- > > Friedrich Foerster > Planegger Str 65 > D-81241 Muenchen > Germany > > frifoe(a)gmail.com > > Tel: +49 89 81305541 > > _______________________________________________ > EMBED-dev mailing list > EMBED-dev(a)salilab.org > https://salilab.org/mailman/listinfo/embed-dev [View Less]

1 0

2D geometry
by Keren Lasker 13 Nov '08

13 Nov '08

Javi, For your 2D geometry classes (Vector2D, Matrix2D ...) - why not just using the corresponding 3D classes and keeping Z fixed? I think it is not too bad in performance and will save you a lot of coding/debugging time. Keren.

5 9

Re: [IMP-dev] [IMP-commits] r854 - trunk/modules/misc/include
by Daniel Russel 13 Nov '08

13 Nov '08

I don't like having two equivalent constructors from thevsame data or using a vector to hold angles. What I had meant with my suggestion is initialization from an axis and an amount of rotation about that axis as that is a somewhat less problematic basis than Euler angles (and one I use :-) On Nov 12, 2008, at 2:31 PM, Notification of IMP commits <imp-commits(a)salilab.org > wrote: > Author: kerenl(a)SALILAB.ORG > Date: 2008-11-12 14:31:59 -0800 (Wed, 12 Nov 2008) > … [View More]New Revision: 854 > > Modified: > trunk/modules/misc/include/Rotation3D.h > Log: > 1. Change protected to private > 2. Add constructors from Vector3D and an angle > 3. update documentation > > Modified: trunk/modules/misc/include/Rotation3D.h > =================================================================== > --- trunk/modules/misc/include/Rotation3D.h 2008-11-12 21:15:36 > UTC (rev 853) > +++ trunk/modules/misc/include/Rotation3D.h 2008-11-12 22:31:59 > UTC (rev 854) > @@ -19,26 +19,28 @@ > > class IMPMISCEXPORT Rotation3D { > public: > - //! Initialize a rotation in x-y-z order from three angles > + Rotation3D(){ > + } > + //! Initialize a rotation in x-y-z order from three Euler angles > /** \param[in] xr Rotation around the X axis > \param[in] yr Rotation around the Y axis > \param[in] zr Rotation around the Z axis > */ > Rotation3D(Float xr, Float yr, Float zr) { > - Float cx = cos(xr); Float cy = cos(yr); Float cz = cos(zr); > - Float sx = sin(xr); Float sy = sin(yr); Float sz = sin(zr); > - Float m00 = cz*cy; > - Float m11 = -sy*sx*sz + cx*cz; > - Float m22 = cy*cx; > - quat_[0] = sqrt(1+m00+m11+m22)/2.0; > - quat_[1] = sqrt(1+m00-m11-m22)/2.0; > - quat_[2] = sqrt(1-m00+m11-m22)/2.0; > - quat_[3] = sqrt(1-m00-m11+m22)/2.0; > - if (cy*sx + sy*cx*sz + sx*cz < 0.0) quat_[1] = -quat_[1]; > - if (sz*sx - sy*cx*cz - sy < 0.0) quat_[2] = -quat_[2]; > - if (sz*cy + sy*sx*cz + sz*cx < 0.0) quat_[3] = -quat_[3]; > + init_angles(xr,yr,zr); > } > - > + //! Initialize a rotation in x-y-z order from three identical > Euler angles > + /** \param[in] e_angle Rotation around first the X axis, Y axis > and Z axis > + */ > + Rotation3D(Float e_angle){ > + init_angles(e_angle, e_angle, e_angle); > + } > + //! Initialize a rotation in x-y-z order from three Euler angles > + /** \param[in] v A vector that holds three Euler angles (x-y-z > order) > + */ > + Rotation3D(const Vector3D &v){ > + init_angles(v[0],v[1],v[2]); > + } > Matrix3D get_matrix() const { > const Float a = quat_[0]; > const Float b = quat_[1]; > @@ -76,7 +78,22 @@ > return atan2(matrix21(), matrix11()); > } > > -protected: > +private: > + void init_angles(Float xr, Float yr, Float zr) { > + Float cx = cos(xr); Float cy = cos(yr); Float cz = cos(zr); > + Float sx = sin(xr); Float sy = sin(yr); Float sz = sin(zr); > + Float m00 = cz*cy; > + Float m11 = -sy*sx*sz + cx*cz; > + Float m22 = cy*cx; > + quat_[0] = sqrt(1+m00+m11+m22)/2.0; > + quat_[1] = sqrt(1+m00-m11-m22)/2.0; > + quat_[2] = sqrt(1-m00+m11-m22)/2.0; > + quat_[3] = sqrt(1-m00-m11+m22)/2.0; > + if (cy*sx + sy*cx*sz + sx*cz < 0.0) quat_[1] = -quat_[1]; > + if (sz*sx - sy*cx*cz - sy < 0.0) quat_[2] = -quat_[2]; > + if (sz*cy + sy*sx*cz + sz*cx < 0.0) quat_[3] = -quat_[3]; > + } > + > Float matrix11() const { > return sqr(quat_[0]) + sqr(quat_[1]) - sqr(quat_[2]) - > sqr(quat_[3]); > } > > _______________________________________________ > IMP-commits mailing list > IMP-commits(a)salilab.org > https://salilab.org/mailman/listinfo/imp-commits [View Less]

4 6

Re: [IMP-dev] IMP-dev Digest, Vol 13, Issue 8
by Javier Ángel Velázquez Muriel 13 Nov '08

13 Nov '08

I don't like the convention x-y-z for the Euler angles and we absolutely need to be consistent around IMP about them or enforce clear specification. The most used one in EM by far is ZYZ and the one that I am going to use. Javi > > > 2008/11/12 <imp-dev-request(a)salilab.org>: >> Send IMP-dev mailing list submissions to >> imp-dev(a)salilab.org >> >> To subscribe or unsubscribe via the World Wide Web, visit >> https://salilab.org/… [View More]mailman/listinfo/imp-dev >> or, via email, send a message with subject or body 'help' to >> imp-dev-request(a)salilab.org >> >> You can reach the person managing the list at >> imp-dev-owner(a)salilab.org >> >> When replying, please edit your Subject line so it is more specific >> than "Re: Contents of IMP-dev digest..." >> >> >> Today's Topics: >> >> 1. Re: [IMP-commits] r854 - trunk/modules/misc/include >> (Dina Schneidman) >> 2. Re: [IMP-commits] r854 - trunk/modules/misc/include >> (Daniel Russel) >> 3. Re: 2D geometry ( Javier ?ngel Vel?zquez Muriel ) >> 4. Re: 2D geometry ( Javier ?ngel Vel?zquez Muriel ) >> 5. Re: 2D geometry (Daniel Russel) >> 6. Re: 2D geometry (Daniel Russel) >> 7. Re: 2D geometry (Keren Lasker) >> >> >> ---------------------------------------------------------------------- >> >> Message: 1 >> Date: Wed, 12 Nov 2008 15:59:13 -0800 >> From: "Dina Schneidman" <duhovka(a)gmail.com> >> Subject: Re: [IMP-dev] [IMP-commits] r854 - trunk/modules/misc/include >> To: "List for IMP development" <imp-dev(a)salilab.org> >> Message-ID: >> <55d63ce50811121559p4995f334kf9e9c1247f21d35f(a)mail.gmail.com> >> Content-Type: text/plain; charset="iso-8859-1" >> >> I think these are Euler angles. >> http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/inde… >> Daniel, I agree that conceptually it is wrong to store them in Vector3D, but >> it is very convenient in practice. >> >> On Wed, Nov 12, 2008 at 10:33 PM, Keren Lasker <kerenl(a)salilab.org> wrote: >> >>> Daniel - the angles are not Euler angles but just three angles around >>> X/Y/Z axis. >>> Would you like the option to rotate around another axis, ,something >>> like: >>> (1,1,0),45 - meaning rotating around the diagonal of the XY plane in >>> 45 degrees? >>> On Nov 12, 2008, at 5:43 PM, Daniel Russel wrote: >>> >>> > I don't like having two equivalent constructors from thevsame data or >>> > using a vector to hold angles. >>> > >>> > What I had meant with my suggestion is initialization from an axis and >>> > an amount of rotation about that axis as that is a somewhat less >>> > problematic basis than Euler angles (and one I use :-) >>> > >>> > >>> > >>> > On Nov 12, 2008, at 2:31 PM, Notification of IMP commits < >>> imp-commits(a)salilab.org >>> >> wrote: >>> > >>> >> Author: kerenl(a)SALILAB.ORG >>> >> Date: 2008-11-12 14:31:59 -0800 (Wed, 12 Nov 2008) >>> >> New Revision: 854 >>> >> >>> >> Modified: >>> >> trunk/modules/misc/include/Rotation3D.h >>> >> Log: >>> >> 1. Change protected to private >>> >> 2. Add constructors from Vector3D and an angle >>> >> 3. update documentation >>> >> >>> >> Modified: trunk/modules/misc/include/Rotation3D.h >>> >> =================================================================== >>> >> --- trunk/modules/misc/include/Rotation3D.h 2008-11-12 21:15:36 >>> >> UTC (rev 853) >>> >> +++ trunk/modules/misc/include/Rotation3D.h 2008-11-12 22:31:59 >>> >> UTC (rev 854) >>> >> @@ -19,26 +19,28 @@ >>> >> >>> >> class IMPMISCEXPORT Rotation3D { >>> >> public: >>> >> - //! Initialize a rotation in x-y-z order from three angles >>> >> + Rotation3D(){ >>> >> + } >>> >> + //! Initialize a rotation in x-y-z order from three Euler angles >>> >> /** \param[in] xr Rotation around the X axis >>> >> \param[in] yr Rotation around the Y axis >>> >> \param[in] zr Rotation around the Z axis >>> >> */ >>> >> Rotation3D(Float xr, Float yr, Float zr) { >>> >> - Float cx = cos(xr); Float cy = cos(yr); Float cz = cos(zr); >>> >> - Float sx = sin(xr); Float sy = sin(yr); Float sz = sin(zr); >>> >> - Float m00 = cz*cy; >>> >> - Float m11 = -sy*sx*sz + cx*cz; >>> >> - Float m22 = cy*cx; >>> >> - quat_[0] = sqrt(1+m00+m11+m22)/2.0; >>> >> - quat_[1] = sqrt(1+m00-m11-m22)/2.0; >>> >> - quat_[2] = sqrt(1-m00+m11-m22)/2.0; >>> >> - quat_[3] = sqrt(1-m00-m11+m22)/2.0; >>> >> - if (cy*sx + sy*cx*sz + sx*cz < 0.0) quat_[1] = -quat_[1]; >>> >> - if (sz*sx - sy*cx*cz - sy < 0.0) quat_[2] = -quat_[2]; >>> >> - if (sz*cy + sy*sx*cz + sz*cx < 0.0) quat_[3] = -quat_[3]; >>> >> + init_angles(xr,yr,zr); >>> >> } >>> >> - >>> >> + //! Initialize a rotation in x-y-z order from three identical >>> >> Euler angles >>> >> + /** \param[in] e_angle Rotation around first the X axis, Y axis >>> >> and Z axis >>> >> + */ >>> >> + Rotation3D(Float e_angle){ >>> >> + init_angles(e_angle, e_angle, e_angle); >>> >> + } >>> >> + //! Initialize a rotation in x-y-z order from three Euler angles >>> >> + /** \param[in] v A vector that holds three Euler angles (x-y-z >>> >> order) >>> >> + */ >>> >> + Rotation3D(const Vector3D &v){ >>> >> + init_angles(v[0],v[1],v[2]); >>> >> + } >>> >> Matrix3D get_matrix() const { >>> >> const Float a = quat_[0]; >>> >> const Float b = quat_[1]; >>> >> @@ -76,7 +78,22 @@ >>> >> return atan2(matrix21(), matrix11()); >>> >> } >>> >> >>> >> -protected: >>> >> +private: >>> >> + void init_angles(Float xr, Float yr, Float zr) { >>> >> + Float cx = cos(xr); Float cy = cos(yr); Float cz = cos(zr); >>> >> + Float sx = sin(xr); Float sy = sin(yr); Float sz = sin(zr); >>> >> + Float m00 = cz*cy; >>> >> + Float m11 = -sy*sx*sz + cx*cz; >>> >> + Float m22 = cy*cx; >>> >> + quat_[0] = sqrt(1+m00+m11+m22)/2.0; >>> >> + quat_[1] = sqrt(1+m00-m11-m22)/2.0; >>> >> + quat_[2] = sqrt(1-m00+m11-m22)/2.0; >>> >> + quat_[3] = sqrt(1-m00-m11+m22)/2.0; >>> >> + if (cy*sx + sy*cx*sz + sx*cz < 0.0) quat_[1] = -quat_[1]; >>> >> + if (sz*sx - sy*cx*cz - sy < 0.0) quat_[2] = -quat_[2]; >>> >> + if (sz*cy + sy*sx*cz + sz*cx < 0.0) quat_[3] = -quat_[3]; >>> >> + } >>> >> + >>> >> Float matrix11() const { >>> >> return sqr(quat_[0]) + sqr(quat_[1]) - sqr(quat_[2]) - >>> >> sqr(quat_[3]); >>> >> } >>> >> >>> >> _______________________________________________ >>> >> IMP-commits mailing list >>> >> IMP-commits(a)salilab.org >>> >> https://salilab.org/mailman/listinfo/imp-commits >>> > _______________________________________________ >>> > IMP-dev mailing list >>> > IMP-dev(a)salilab.org >>> > https://salilab.org/mailman/listinfo/imp-dev >>> >>> _______________________________________________ >>> IMP-dev mailing list >>> IMP-dev(a)salilab.org >>> https://salilab.org/mailman/listinfo/imp-dev >>> >> [View Less]

4 4

Re: [IMP-dev] IMP-dev Digest, Vol 13, Issue 9
by Javier Ángel Velázquez Muriel 12 Nov '08

12 Nov '08

Actually, on second thought, what I am working on is not necessary a Matrix in 2D or 3D as Keren has wrote in her class. I am working in a template array able to store different types of data and then access to them in a 1D, 2D or 3D fashion. It is the way images and volumes are treated in Xmipp, and it is very efficient. Your point of view, guys, is more of storing points, and therefore compatible with what I am doing. J

3 4

Re: [IMP-dev] Build system failure, 11/12/2008
by Keren Lasker 12 Nov '08

12 Nov '08

hi, It seems that the problem was with my recent updates, for some reason Transformation3D was not recognized as a module of misc. I'll fix it today as well as addressing Daniel's comments. Keren. On Nov 12, 2008, at 10:00 AM, build-system(a)salilab.org wrote: > Build system failed on 11/12/2008. > Failed components: MODELLER, IMP. > Please see https://salilab.org/internal/imp/tests.html for full > details. > _______________________________________________ > IMP-dev … [View More]

2 2

VectorD
by Daniel Russel 12 Nov '08

12 Nov '08

As part of the discussed reorg of vector and matrix stuff, here is the patch to switch to a d-dimensional implementation for vectors. In order to do nice error checking it adds a preprocessor macro to the _wrap compilation which is used to disable the compile time checks (since SWIG instantiates all methods, for obvious reasons) and replace them with runtime checks. For some reason I need a typedef for Vector3D in each of the module .i files. I don't see why, but I can't get things to … [View More]work without it. Index: kernel/test/misc/test_vector3d.py =================================================================== --- kernel/test/misc/test_vector3d.py (revision 855) +++ kernel/test/misc/test_vector3d.py (working copy) @@ -94,6 +94,11 @@ for i in range(3): self.assertEqual(prod[i], expected_prod[i]) self.assertEqual(v1[i], expected_prod[i]) + def test_generators(self): + """Check the Vector3D generators""" + # test calling since it is a bit non-trivial in SWIG + v= IMP.random_vector_in_unit_sphere() + v= IMP.random_vector_in_sphere(IMP.Vector3D(0,0,0), 1) if __name__ == '__main__': unittest.main() Index: kernel/include/Vector3D.h =================================================================== --- kernel/include/Vector3D.h (revision 855) +++ kernel/include/Vector3D.h (working copy) @@ -8,7 +8,7 @@ #ifndef IMP_VECTOR_3D_H #define IMP_VECTOR_3D_H -#include "IMP_config.h" +#include "VectorD.h" #include "base_types.h" #include "macros.h" #include "exception.h" @@ -17,199 +17,8 @@ IMP_BEGIN_NAMESPACE -//! Simple 3D vector class -/** \ingroup helper - */ -class IMPDLLEXPORT Vector3D -{ - bool is_default() const {return false;} -public: - // public for swig - typedef Vector3D This; +typedef VectorD<3> Vector3D; - //! Initialize the vector from separate x,y,z values. - Vector3D(Float x, Float y, Float z) { - vec_[0] = x; - vec_[1] = y; - vec_[2] = z; - } - - //! Default constructor - Vector3D() {} - - IMP_COMPARISONS_3(vec_[0], vec_[1], vec_[2]); - - //! \return A single component of this vector (0-2). - Float operator[](unsigned int i) const { - IMP_assert(i < 3, "Invalid component of vector requested"); - return vec_[i]; - } - - //! \return A single component of this vector (0-2). - Float& operator[](unsigned int i) { - IMP_assert(i < 3, "Invalid component of vector requested"); - return vec_[i]; - } - - //! \return the scalar product of two vectors. - /** \param[in] vec2 The other vector to use in the product. - */ - Float scalar_product(const Vector3D &vec2) const { - return vec_[0] * vec2.vec_[0] + vec_[1] * vec2.vec_[1] - + vec_[2] * vec2.vec_[2]; - } - - //! scalar product - Float operator*(const Vector3D &o) const { - return scalar_product(o); - } - - //! product with scalar - Vector3D operator*(Float s) const { - return Vector3D(operator[](0) * s, - operator[](1) * s, - operator[](2) * s); - } - - //! divide by a scalar - Vector3D operator/(Float s) const { - return Vector3D(operator[](0) / s, - operator[](1) / s, - operator[](2) / s); - } - - //! negation - Vector3D operator-() const { - return Vector3D(-operator[](0), - -operator[](1), - -operator[](2)); - } - - //! \return the vector product of two vectors. - /** \param[in] vec2 The other vector to use in the product. - */ - Vector3D vector_product(const Vector3D &vec2) const { - return Vector3D(vec_[1] * vec2.vec_[2] - vec_[2] * vec2.vec_[1], - vec_[2] * vec2.vec_[0] - vec_[0] * vec2.vec_[2], - vec_[0] * vec2.vec_[1] - vec_[1] * vec2.vec_[0]); - } - - //! \return The square of the magnitude of this vector. - Float get_squared_magnitude() const { - return vec_[0] * vec_[0] + vec_[1] * vec_[1] + vec_[2] * vec_[2]; - } - - //! \return The magnitude of this vector. - Float get_magnitude() const { - return std::sqrt(get_squared_magnitude()); - } - - //! \return This vector normalized to unit length. - Vector3D get_unit_vector() const { - Float mag = get_magnitude(); - // avoid division by zero - mag = std::max(mag, static_cast<Float>(1e-12)); - return Vector3D(vec_[0] / mag, vec_[1] / mag, vec_[2] / mag); - } - - //! \return Difference between two vectors. - Vector3D operator-(const Vector3D &o) const { - return Vector3D(operator[](0)-o[0], - operator[](1)-o[1], - operator[](2)-o[2]); - } - - //! \return Sum of two vectors. - Vector3D operator+(const Vector3D &o) const { - return Vector3D(operator[](0) + o[0], - operator[](1) + o[1], - operator[](2) + o[2]); - } - - //! Accumulate the vector - Vector3D& operator+=(const Vector3D &o) { - vec_[0] += o[0]; - vec_[1] += o[1]; - vec_[2] += o[2]; - return *this; - } - - //! Rescale the vector - Vector3D& operator/=(Float f) { - vec_[0] /= f; - vec_[1] /= f; - vec_[2] /= f; - return *this; - } - - //! Rescale the vector - Vector3D& operator*=(Float f) { - vec_[0] *= f; - vec_[1] *= f; - vec_[2] *= f; - return *this; - } - - void show(std::ostream &out=std::cout) const { - out << "(" << operator[](0) << ", " << operator[](1) << ", " - << operator[](2) << ")"; - } - -private: - Float vec_[3]; -}; - -IMP_OUTPUT_OPERATOR(Vector3D); - - -//! product with scalar -inline Vector3D operator*(Float s, const Vector3D &o) { - return Vector3D(o[0]*s, - o[1]*s, - o[2]*s); -} - -//! Generate a random vector in a box with uniform density -IMPDLLEXPORT Vector3D -random_vector_in_box(const Vector3D &lb=Vector3D(0,0,0), - const Vector3D &ub=Vector3D(10,10,10)); - -//! Generate a random vector in a sphere with uniform density -IMPDLLEXPORT Vector3D -random_vector_in_sphere(const Vector3D &center=Vector3D(0,0,0), - Float radius=1); - -//! Generate a random vector on a sphere with uniform density -IMPDLLEXPORT Vector3D -random_vector_on_sphere(const Vector3D &center=Vector3D(0,0,0), - Float radius=1); - -struct SpacesIO -{ - const Vector3D &v_; - SpacesIO(const Vector3D &v): v_(v){} -}; - - -inline std::ostream &operator<<(std::ostream &out, const SpacesIO &s) -{ - out << s.v_[0] << " " << s.v_[1] << " " << s.v_[2]; - return out; -} - -struct CommasIO -{ - const Vector3D &v_; - CommasIO(const Vector3D &v): v_(v){} -}; - - -inline std::ostream &operator<<(std::ostream &out, const CommasIO &s) -{ - out << s.v_[0] << ", " << s.v_[1] << ", " << s.v_[2]; - return out; -} - IMP_END_NAMESPACE #endif /* IMP_VECTOR_3D_H */ Index: kernel/include/SConscript =================================================================== --- kernel/include/SConscript (revision 855) +++ kernel/include/SConscript (working copy) @@ -4,7 +4,7 @@ 'Particle.h', 'ScoreState.h', 'OptimizerState.h', 'IMP_config.h', 'log.h', 'DerivativeAccumulator.h', 'Key.h', 'utility.h', 'Restraint.h', 'Optimizer.h', - 'DecoratorBase.h', 'Vector3D.h', + 'DecoratorBase.h', 'Vector3D.h', 'VectorD.h', 'UnaryFunction.h', 'PairScore.h', 'SingletonScore.h', 'macros.h', 'TripletScore.h', 'exception.h', 'VersionInfo.h', 'Object.h', 'Pointer.h', 'RefCountedObject.h', 'ParticleRefiner.h'] Index: kernel/include/VectorD.h =================================================================== --- kernel/include/VectorD.h (revision 0) +++ kernel/include/VectorD.h (revision 0) @@ -0,0 +1,365 @@ +/** + * \File VectorD.h \brief Simple D vector class. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#ifndef IMP_VECTOR_D_H +#define IMP_VECTOR_D_H + +#include "IMP_config.h" +#include "base_types.h" +#include "macros.h" +#include "exception.h" +#include "random.h" +#include "internal/constants.h" + +#include <boost/random/uniform_real.hpp> + +#include <cmath> + +IMP_BEGIN_NAMESPACE + +//! Simple D vector class +/** \ingroup helper + */ +template <unsigned int D> +class IMPDLLEXPORT VectorD +{ + bool is_default() const {return false;} +public: + // public for swig + typedef VectorD<D> This; + + //! Initialize the 3-vector from separate x,y,z values. + VectorD(Float x, Float y) { +#ifdef SWIG_WRAPPER + IMP_check(D==2, "Need " << D << " to construct a " + << D << "-vector.", ValueException); +#else + BOOST_STATIC_ASSERT(D==2); +#endif + vec_[0] = x; + vec_[1] = y; + } + + //! Initialize the 2-vector from separate x,y values. + VectorD(Float x, Float y, Float z) { +#ifdef SWIG_WRAPPER + IMP_check(D==3, "Need " << D << " to construct a " + << D << "-vector.", ValueException); +#else + BOOST_STATIC_ASSERT(D==3); +#endif + vec_[0] = x; + vec_[1] = y; + vec_[2] = z; + } + + //! Initialize the 2-vector from separate x,y values. + VectorD(Float w, Float x, Float y, Float z) { +#ifdef SWIG_WRAPPER + IMP_check(D==4, "Need " << D << " to construct a " + << D << "-vector.", ValueException); +#else + BOOST_STATIC_ASSERT(D==4); +#endif + vec_[0] = w; + vec_[1] = x; + vec_[2] = y; + vec_[2] = z; + } + + //! Default constructor + VectorD() {} + + IMP_COMPARISONS; + + //! \return A single component of this vector (0-D). + Float operator[](unsigned int i) const { + IMP_assert(i < D, "Invalid component of vector requested"); + return vec_[i]; + } + + //! \return A single component of this vector (0-D). + Float& operator[](unsigned int i) { + IMP_assert(i < D, "Invalid component of vector requested"); + return vec_[i]; + } + + //! \return the scalar product of two vectors. + /** \param[in] vec2 The other vector to use in the product. + */ + Float scalar_product(const This &o) const { + Float ret=0; + for (unsigned int i=0; i< D; ++i) { + ret += vec_[i]* o.vec_[i]; + } + return ret; + } + + //! scalar product + Float operator*(const This &o) const { + return scalar_product(o); + } + + //! product with scalar + VectorD operator*(Float s) const { + This ret; + for (unsigned int i=0; i< D; ++i) { + ret.vec_[i] = vec_[i] * s; + } + return ret; + } + + //! divide by a scalar + VectorD operator/(Float s) const { + This ret; + for (unsigned int i=0; i< D; ++i) { + ret.vec_[i] = vec_[i] / s; + } + return ret; + } + + //! negation + VectorD operator-() const { + This ret; + for (unsigned int i=0; i< D; ++i) { + ret.vec_[i] = -vec_[i]; + } + return ret; + } + + //! \return the vector product of two vectors. + /** \param[in] vec2 The other vector to use in the product. + */ + VectorD vector_product(const VectorD &vec2) const { + BOOST_STATIC_ASSERT(D==3); + return VectorD(vec_[1] * vec2.vec_[2] - vec_[2] * vec2.vec_[1], + vec_[2] * vec2.vec_[0] - vec_[0] * vec2.vec_[2], + vec_[0] * vec2.vec_[1] - vec_[1] * vec2.vec_[0]); + } + + //! \return The square of the magnitude of this vector. + Float get_squared_magnitude() const { + return scalar_product(*this); + } + + //! \return The magnitude of this vector. + Float get_magnitude() const { + return std::sqrt(get_squared_magnitude()); + } + + //! \return This vector normalized to unit length. + VectorD get_unit_vector() const { + Float mag = get_magnitude(); + // avoid division by zero + mag = std::max(mag, static_cast<Float>(1e-12)); + return operator/(mag); + } + + //! \return Difference between two vectors. + VectorD operator-(const VectorD &o) const { + This ret; + for (unsigned int i=0; i< D; ++i) { + ret.vec_[i] = vec_[i] - o.vec_[i]; + } + return ret; + } + + //! \return Sum of two vectors. + VectorD operator+(const VectorD &o) const { + This ret; + for (unsigned int i=0; i< D; ++i) { + ret.vec_[i] = vec_[i] + o.vec_[i]; + } + return ret; + } + + //! Accumulate the vector + VectorD& operator+=(const VectorD &o) { + for (unsigned int i=0; i< D; ++i) { + vec_[i] += o[i]; + } + return *this; + } + + //! Rescale the vector + VectorD& operator/=(Float f) { + for (unsigned int i=0; i< D; ++i) { + vec_[i] /= f; + } + return *this; + } + + //! Rescale the vector + VectorD& operator*=(Float f) { + for (unsigned int i=0; i< D; ++i) { + vec_[i] *= f; + } + return *this; + } + + void show(std::ostream &out=std::cout, std::string delim=", ") const { + out << "("; + for (unsigned int i=0; i< D; ++i) { + out << vec_[i]; + if (i != D-1) { + out << delim; + } + } + out << ")"; + } + +private: + int compare(const This &o) const { + for (unsigned int i=0; i< D; ++i) { + if (vec_[i] < o.vec_[i]) return -1; + else if (vec_[i] > o.vec_[i]) return 1; + } + return 0; + } + + Float vec_[D]; +}; + +template <unsigned int D> +std::ostream &operator<<(std::ostream &out, const VectorD<D> &v) { + v.show(out); + return out; +} + + +//! product with scalar +template <unsigned int D> +inline VectorD<D> operator*(Float s, const VectorD<D> &o) { + return o*s; +} + +//! create a constant vector +/** This is not the right name. + */ +template <unsigned int D> +VectorD<D> constant_vector(Float s) { + VectorD<D> ret; + for (unsigned int i= 0; i < D; ++i) { + ret[i]=s; + } + return ret; +} + + +//! Generate a random vector in a box with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_in_box(const VectorD<D> &lb, + const VectorD<D> &ub) { + VectorD<D> ret; + for (unsigned int i=0; i< D; ++i) { + IMP_check(lb[i] < ub[i], "Box for randomize must be non-empty", + ValueException); + ::boost::uniform_real<> rand(lb[i], ub[i]); + ret[i]=rand(random_number_generator); + } + return ret; +} + +//! Generate a random vector in a box with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_in_unit_box() { + return random_vector_in_box(VectorD<D>(0,0,0), + VectorD<D>(1,1,1)); +} + +//! Generate a random vector in a sphere with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_in_sphere(const VectorD<D> &center, + Float radius){ + IMP_check(radius > 0, "Radius in randomize must be postive", + ValueException); + VectorD<D> rad= constant_vector<D>(radius); + VectorD<D> min= center - rad; + VectorD<D> max= center + rad; + Float norm; + VectorD<D> ret; + // \todo This algorithm could be more efficient. + do { + ret=random_vector_in_box(min, max); + norm= (center- ret).get_magnitude(); + } while (norm > radius); + return ret; +} + +//! Generate a random vector in a unit sphere with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_in_unit_sphere(){ + return random_vector_in_sphere(VectorD<D>(0,0,0), 1); +} + +//! Generate a random vector on a sphere with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_on_sphere(const VectorD<D> &center, + Float radius) { + // could be made general + BOOST_STATIC_ASSERT(D==3); + IMP_check(radius > 0, "Radius in randomize must be postive", + ValueException); + ::boost::uniform_real<> rand(-1,1); + VectorD<D> up; + up[2]= rand(random_number_generator); + ::boost::uniform_real<> trand(0, 2*IMP::internal::PI); + Float theta= trand(random_number_generator); + // radius of circle + Float r= std::sqrt(1-square(up[2])); + up[0]= std::sin(theta)*r; + up[1]= std::cos(theta)*r; + IMP_assert(std::abs(up.get_magnitude() -1) < .1, + "Error generating unit vector on sphere"); + IMP_LOG(VERBOSE, "Random vector on sphere is " << up << std::endl); + return center+ up*radius; +} + +//! Generate a random vector on a sphere with uniform density +template <unsigned int D> +IMPDLLEXPORT VectorD<D> +random_vector_on_unit_sphere() { + return random_vector_on_sphere(VectorD<D>(0,0,0), 1); +} + + +template <unsigned int D> +struct SpacesIO +{ + const VectorD<D> &v_; + SpacesIO(const VectorD<D> &v): v_(v){} +}; + +template <unsigned int D> +inline std::ostream &operator<<(std::ostream &out, const SpacesIO<D> &s) +{ + s.v_.show(out, " "); + return out; +} + +template <unsigned int D> +struct CommasIO +{ + const VectorD<D> &v_; + CommasIO(const VectorD<D> &v): v_(v){} +}; + +template <unsigned int D> +inline std::ostream &operator<<(std::ostream &out, const CommasIO<D> &s) +{ + s.v_.show(out, ", "); + return out; +} + +IMP_END_NAMESPACE + +#endif /* IMP_VECTOR_D_H */ Index: kernel/include/macros.h =================================================================== --- kernel/include/macros.h (revision 855) +++ kernel/include/macros.h (working copy) @@ -8,6 +8,39 @@ #ifndef IMP_MACROS_H #define IMP_MACROS_H +//! Implement comparison in a class using a compare function +/** The macro requires that This be defined as the type of the current class. + The compare function should take a const This & and return -1, 0, 1 as + appropriate. + */ +#define IMP_COMPARISONS \ + /** */ bool operator==(const This &o) const { \ + return (compare(o) == 0); \ + } \ + /** */ bool operator!=(const This &o) const { \ + return (compare(o) != 0); \ + } \ + /** */ bool operator<(const This &o) const { \ + IMP_assert(!is_default() && !o.is_default(), \ + "Ordering with uninitialized index is undefined"); \ + return (compare(o) <0); \ + } \ + /** */ bool operator>(const This &o) const { \ + IMP_assert(!is_default() && !o.is_default(), \ + "Ordering with uninitialized index is undefined"); \ + return (compare(o) > 0); \ + } \ + /** */ bool operator>=(const This &o) const { \ + IMP_assert(!is_default() && !o.is_default(), \ + "Ordering with uninitialized index is undefined"); \ + return !(compare(o) < 0); \ + } \ + /** */ bool operator<=(const This &o) const { \ + IMP_assert(!is_default() && !o.is_default(), \ + "Ordering with uninitialized index is undefined"); \ + return !(compare(o) > 0); \ + } + //! Implement comparison in a class using field as the variable to compare /** The macro requires that This be defined as the type of the current class. */ Index: kernel/include/internal/SConscript =================================================================== --- kernel/include/internal/SConscript (revision 855) +++ kernel/include/internal/SConscript (working copy) @@ -1,7 +1,15 @@ -files = ['AttributeTable.h', 'Vector.h', - 'ref_counting.h', 'ObjectContainer.h', - 'kernel_version_info.h', 'constants.h', 'units.h', - 'utility.h', 'Unit.h', 'ExponentialNumber.h'] +files = [ + 'AttributeTable.h', + 'ExponentialNumber.h', + 'ObjectContainer.h', + 'Unit.h', + 'Vector.h', + 'constants.h', + 'kernel_version_info.h', + 'ref_counting.h', + 'units.h', + 'utility.h', + ] -files = [File(f) for f in files] +files = [File(x) for x in files] Return('files') Index: kernel/src/SConscript =================================================================== --- kernel/src/SConscript (revision 855) +++ kernel/src/SConscript (working copy) @@ -12,7 +12,7 @@ 'Particle.cpp', 'ScoreState.cpp', 'Object.cpp', 'OptimizerState.cpp', 'Log.cpp', 'Restraint.cpp', 'Optimizer.cpp', 'random.cpp', 'Key.cpp', 'exception.cpp', 'ParticleRefiner.cpp', - 'Vector3D.cpp', 'UnaryFunction.cpp', 'PairScore.cpp', + 'UnaryFunction.cpp', 'PairScore.cpp', 'SingletonScore.cpp', 'TripletScore.cpp' ] + internal_files Index: kernel/src/Vector3D.cpp =================================================================== --- kernel/src/Vector3D.cpp (revision 855) +++ kernel/src/Vector3D.cpp (working copy) @@ -1,65 +0,0 @@ -/** - * \file Vector3D.cpp \brief Simple 3D vector class. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - * - */ - -#include "IMP/Vector3D.h" -#include "IMP/random.h" -#include "IMP/internal/constants.h" -#include "IMP/utility.h" - -#include <boost/random/uniform_real.hpp> - -IMP_BEGIN_NAMESPACE - -Vector3D random_vector_in_box(const Vector3D &min, const Vector3D &max) -{ - Vector3D ret; - for (unsigned int i=0; i< 3; ++i) { - IMP_check(min[i] < max[i], "Box for randomize must be non-empty", - ValueException); - ::boost::uniform_real<> rand(min[i], max[i]); - ret[i]=rand(random_number_generator); - } - return ret; -} - - -Vector3D random_vector_in_sphere(const Vector3D &center, Float radius) -{ - IMP_check(radius > 0, "Radius in randomize must be postive", - ValueException); - Vector3D min(center[0]-radius, center[1]-radius, center[2]-radius); - Vector3D max(center[0]+radius, center[1]+radius, center[2]+radius); - float norm; - Vector3D ret; - // \todo This algorithm could be more efficient. - do { - ret=random_vector_in_box(min, max); - norm= (center- ret).get_magnitude(); - } while (norm > radius); - return ret; -} - -Vector3D random_vector_on_sphere(const Vector3D &center, Float radius) -{ - IMP_check(radius > 0, "Radius in randomize must be postive", - ValueException); - ::boost::uniform_real<> rand(-1,1); - Vector3D up; - up[2]= rand(random_number_generator); - ::boost::uniform_real<> trand(0, 2*internal::PI); - Float theta= trand(random_number_generator); - // radius of circle - Float r= std::sqrt(1-square(up[2])); - up[0]= std::sin(theta)*r; - up[1]= std::cos(theta)*r; - IMP_assert(std::abs(up.get_magnitude() -1) < .1, - "Error generating unit vector on sphere"); - IMP_LOG(VERBOSE, "Random vector on sphere is " << up << std::endl); - return center+ up*radius; -} - -IMP_END_NAMESPACE Index: kernel/src/internal/SConscript =================================================================== --- kernel/src/internal/SConscript (revision 855) +++ kernel/src/internal/SConscript (working copy) @@ -1,6 +1,7 @@ -Import('env') +files = [ + 'constants.cpp', + 'kernel_version_info.cpp', + ] -files = ['kernel_version_info.cpp', 'constants.cpp'] - files = [File(x) for x in files] Return('files') Index: kernel/pyext/Vector3D.i =================================================================== --- kernel/pyext/Vector3D.i (revision 855) +++ kernel/pyext/Vector3D.i (working copy) @@ -1,30 +1,31 @@ + /* Provide our own implementations for some operators */ -%ignore IMP::Vector3D::operator[]; -%ignore IMP::Vector3D::operator+=; -%ignore IMP::Vector3D::operator*=; -%ignore IMP::Vector3D::operator/=; +%ignore IMP::VectorD::operator[]; +%ignore IMP::VectorD::operator+=; +%ignore IMP::VectorD::operator*=; +%ignore IMP::VectorD::operator/=; /* Make sure that we return the original Python object from C++ inplace operators (not a new Python proxy around the same C++ object) */ namespace IMP { - %feature("shadow") Vector3D::__iadd__(const Vector3D &) %{ + %feature("shadow") VectorD::__iadd__(const VectorD &) %{ def __iadd__(self, *args): $action(self, *args) return self %} - %feature("shadow") Vector3D::__imul__(Float) %{ + %feature("shadow") VectorD::__imul__(Float) %{ def __imul__(self, *args): $action(self, *args) return self %} - %feature("shadow") Vector3D::__idiv__(Float) %{ + %feature("shadow") VectorD::__idiv__(Float) %{ def __idiv__(self, *args): $action(self, *args) return self %} } -%extend IMP::Vector3D { +%extend IMP::VectorD { Float __getitem__(unsigned int index) const { return self->operator[](index); } @@ -33,9 +34,19 @@ } /* Ignore C++ return value from inplace operators, so that SWIG does not generate a new SWIG wrapper for the return value (see above). */ - void __iadd__(const Vector3D &o) { self->operator+=(o); } + void __iadd__(const VectorD &o) { self->operator+=(o); } void __imul__(Float f) { self->operator*=(f); } void __idiv__(Float f) { self->operator/=(f); } }; -%include "IMP/Vector3D.h" +%include "IMP/VectorD.h" + +namespace IMP { + %template(Vector3D) VectorD<3>; + %template(random_vector_on_sphere) random_vector_on_sphere<3>; + %template(random_vector_in_sphere) random_vector_in_sphere<3>; + %template(random_vector_in_box) random_vector_in_box<3>; + %template(random_vector_on_unit_sphere) random_vector_on_unit_sphere<3>; + %template(random_vector_in_unit_sphere) random_vector_in_unit_sphere<3>; + %template(random_vector_in_unit_box) random_vector_in_unit_box<3>; +} \ No newline at end of file Index: kernel/pyext/SConscript =================================================================== --- kernel/pyext/SConscript (revision 855) +++ kernel/pyext/SConscript (working copy) @@ -5,6 +5,7 @@ # Get a modified build environment suitable for building Python extensions: e = get_pyext_environment(env, 'IMP', cplusplus=True) e.Append(CPPPATH=['#/build/include', e['BOOST_CPPPATH']]) +e.Append(CPPDEFINES=['SWIG_WRAPPER']) e.Append(LIBPATH=['#/build/lib'], LIBS='imp') # Build the Python extension from SWIG interface file: Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 855) +++ kernel/pyext/IMP.i (working copy) @@ -97,6 +97,7 @@ %feature("director") IMP::ParticleRefiner; %include "IMP/base_types.h" +%include "Vector3D.i" %include "IMP/Object.h" %include "IMP/RefCountedObject.h" %include "IMP/Index.h" @@ -113,7 +114,6 @@ %include "IMP/SingletonScore.h" %include "IMP/TripletScore.h" %include "Particle.i" -%include "Vector3D.i" %include "IMP/DecoratorBase.h" %include "IMP/Optimizer.h" @@ -137,4 +137,8 @@ %template(Floats) ::std::vector<Float>; %template(Strings) ::std::vector<String>; %template(Ints) ::std::vector<Int>; + /*%template(Vector3D) VectorD<3>; + %template(random_vector_on_sphere) random_vector_on_sphere<3>; + %template(random_vector_in_sphere) random_vector_in_sphere<3>; + %template(random_vector_in_box) random_vector_in_box<3>;*/ } Index: modules/core/test/states/test_close_pairs_score_state.py =================================================================== --- modules/core/test/states/test_close_pairs_score_state.py (revision 855) +++ modules/core/test/states/test_close_pairs_score_state.py (working copy) @@ -59,7 +59,7 @@ for p in ps: d= IMP.core.XYZRDecorator.cast(p) - d.set_coordinates(IMP.random_vector_in_box()) + d.set_coordinates(IMP.random_vector_in_unit_box()) self._compare_lists(m, cpss) print "changing radius" cpss.set_radius_key(IMP.core.XYZDecorator.get_xyz_keys()[0]) Index: modules/core/test/states/test_close_bipartite_pairs_score_state.py =================================================================== --- modules/core/test/states/test_close_bipartite_pairs_score_state.py (revision 855) +++ modules/core/test/states/test_close_bipartite_pairs_score_state.py (working copy) @@ -64,12 +64,12 @@ for p in ps0: d= IMP.core.XYZRDecorator.cast(p) - d.set_coordinates(IMP.random_vector_in_box()) + d.set_coordinates(IMP.random_vector_in_unit_box()) self._compare_lists(m, cpss) for p in ps1: d= IMP.core.XYZRDecorator.cast(p) - d.set_coordinates(IMP.random_vector_in_box()) + d.set_coordinates(IMP.random_vector_in_unit_box()) self._compare_lists(m, cpss) print "changing radius" Index: modules/core/test/states/test_maximum_change.py =================================================================== --- modules/core/test/states/test_maximum_change.py (revision 855) +++ modules/core/test/states/test_maximum_change.py (working copy) @@ -25,7 +25,7 @@ print "perturbing" for i in range(0,len(ps)): d= IMP.core.XYZDecorator(ps[i]) - v= IMP.random_vector_in_sphere() + v= IMP.random_vector_in_unit_sphere() print v[0] mmax=max(mmax, float(v[0]), float(v[1]), float(v[2])) d.set_coordinates(d.get_coordinates()+v) Index: modules/core/pyext/core.i =================================================================== --- modules/core/pyext/core.i (revision 855) +++ modules/core/pyext/core.i (working copy) @@ -57,6 +57,10 @@ %import "kernel/pyext/IMP_keys.i" namespace IMP { + typedef VectorD<3> Vector3D; +} + +namespace IMP { namespace core { %template(AtomTypeBase) ::IMP::KeyBase<IMP_ATOM_TYPE_INDEX>; %template(ResidueTypeBase) ::IMP::KeyBase<IMP_RESIDUE_TYPE_INDEX>; Index: modules/misc/pyext/misc.i =================================================================== --- modules/misc/pyext/misc.i (revision 857) +++ modules/misc/pyext/misc.i (working copy) @@ -17,6 +17,10 @@ /* Get definitions of kernel base classes (but do not wrap) */ %import "kernel/pyext/IMP.i" +namespace IMP { + typedef VectorD<3> Vector3D; +} + namespace IMP::misc { IMP_OWN_FIRST_CONSTRUCTOR(BondCoverPairScore) IMP_OWN_FIRST_SECOND_CONSTRUCTOR(RefineOncePairScore) [View Less]

1 0