Why do we have two random rotation methods?

[tjlane]

@dermen is one better? We should choose the best and keep that one. I guess even if best is just fastest. Having two is confusing!

https://github.com/tjlane/thor/blob/master/src/python/structure.py#L190
https://github.com/tjlane/thor/blob/master/src/python/structure.py#L227

[tjlane]

@dermen ping ping ping

[tjlane]

I'm just being annoying until you respond. We'll see if I can out-annoy Daniel to get your attention first :D.

[dermen]

got this just now!

[dermen]

the method thor.xray.structure.rand_rotate_molecule2 uses the method thor.xray.math2.rand_rot which I believe improves speed, so I would stick with that one..

Check the documentation in math2.rand_rot for a reference.

[tjlane]

@dermen are there some test results you have that show it being faster AND correct? I think we should probably delete one of the two methods and stick with the better one. You agree?

[dermen]

Yea, I did both of these, probably in old odin issues.

both passed @rmcgibbo 's gist tests and rand_rotate_molecule2 was slightly faster

[dermen]

To check for speed I did

In [14]: xyz = mdtraj.load_pdb('12mr.pdb').xyz[0]

In [15]: from thor.structure import rand_rotate_molecule

In [16]: from thor.structure import rand_rotate_molecule2

In [17]: %timeit( rand_rotate_molecule(xyz) )
10 loops, best of 3: 23.6 ms per loop

In [18]: %timeit( rand_rotate_molecule2(xyz) )
1000 loops, best of 3: 828 us per loop

To check for correctness I used the test code the @rmcgibbo created a while back

first file: test_rotations.py:

import itertools
import os

import numpy as np
from scipy.stats import chi2 as chi2distribution
try:
    import matplotlib.pyplot as pp
    HAVE_PYPLOT = True
except ImportError:
    HAVE_PYPLOT = False

from thor.math2 import rand_rot



def test_RandRot():
    "Test the RandRot() function"
    for t in _test_random_rotation(rand_rot):
        yield t


def _test_random_rotation(rotation_matrix_factory, n_tests=200, n_bins=20):
    """Main test driver. Takes as input a function that generates random
    rotation matricies and tries rotating a bunch of vectors. In polar coordinates,
    we check via a chi2 test whether the polar angles and azimuthal angles
    appear to be distributed as they should.
    """

    v = np.array([1, 0, 0])
    rotated = np.empty((n_tests, 3))

    for i in range(n_tests):
        R = rotation_matrix_factory()
        rotated[i] = np.dot(v, R)
    r = np.sqrt(np.sum(np.square(rotated), axis=1))

    # polar angle theta (0, pi)
    theta = np.arccos(rotated[:, 2] / r)

    # azimuthal angle phi [0, 2pi)
    phi = np.arctan2(rotated[:,1], rotated[:,0])
    phi[phi < 0.0] += 2 * np.pi

    np.testing.assert_array_almost_equal(r, np.ones_like(r),
        err_msg='the length of each of the vectors should still be one')

    def assert_phi_chi2():
        # do a chi-squared test that the phi angles are uniform from [0, 2*pi)
        observed_counts, bin_edges = np.histogram(phi, bins=n_bins, range=(0, 2*np.pi))
        expected_counts = n_tests * np.diff(bin_edges) / (2*np.pi)
        assert_chisquared(observed_counts, expected_counts, bin_edges, title='phi (azimuthal angle)')

    def assert_theta_chi2():
        # marginalized over phi, the number of points at each polar angle theta
        # should follow a sin curve. The differential surface area element
        # is r^2 \sin \theta d\theta d\phi, so integrating over \phi gives us
        # dA(theta) = 2*pi * r**2 * sin(theta) * dtheta    
        # after binning, then we should have -2*pi*cos(b)
        observed_counts, bin_edges = np.histogram(theta, bins=n_bins, range=(0, np.pi))
        # the number of expected counts in each bin comes from integrating dA(theta)
        # over the theta bin, and then properly normalizing
        bin_areas = np.array([-2*np.pi*np.cos(bin_edges[i+1]) + 2*np.pi*np.cos(bin_edges[i]) for i in range(n_bins)])
        assert np.sum(bin_areas) - 4*np.pi < 1e-10
        expected_counts = n_tests * bin_areas / (4*np.pi)

        assert_chisquared(observed_counts, expected_counts, bin_edges, alpha=0.05, title='theta (polar angle)')

    yield assert_phi_chi2
    yield assert_theta_chi2


def assert_chisquared(observed, expected, bin_edges=None, alpha=0.05, title=''):
    """Assert that the "observed" counts are close enough the "expected"
    counts with a chi2 test at the given confidence level. If the test fails,
    we'll save a histogram to disk.
    """


    chi2_value = np.sum((observed - expected)**2 / expected)

    n_dof = len(observed) - 1 # number of degrees of freedom for the chi2 test
    pval = 1 - chi2distribution.cdf(chi2_value, n_dof)

    print 'observed'
    print observed
    print 'expected'
    print expected
    print 'pval', pval

    if pval < alpha:
        if HAVE_PYPLOT:
            pp.clf()
            pp.title(title)
            pp.bar(bin_edges[0:-1], observed, width=bin_edges[1] - bin_edges[0],
                   alpha=0.5, label='observed', color='red')
            pp.bar(bin_edges[0:-1], expected, width=bin_edges[1] - bin_edges[0],
                   alpha=0.5, label='expected', color='blue')
            pp.legend()
            for i in itertools.count():
                path = 'hist-%d.png' % i
                if not os.path.exists(path):
                    break
            pp.savefig(path)
        raise ValueError('p=%f (<%f), we reject null hypothesis that the distribution '
            'matches the expected one. saved histogram as %s' % (pval, alpha, path))

and the second file rotations.py:

import numpy as np

# Portions of this code are from transformations.py, by Christoph Gohlke
# available at http://www.lfd.uci.edu/~gohlke/code/transformations.py.html, and
# released under the following 3-clause BSD terms:

# Copyright (c) 2006-2013, Christoph Gohlke
# Copyright (c) 2006-2013, The Regents of the University of California
# Produced at the Laboratory for Fluorescence Dynamics
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright
#   notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
#   notice, this list of conditions and the following disclaimer in the
#   documentation and/or other materials provided with the distribution.
# * Neither the name of the copyright holders nor the names of any
#   contributors may be used to endorse or promote products derived
#   from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

# Other portions of this code are from Pyrr, which is released under a 2-clause
# MIT license:

# Pyrr is developed by Twisted Pair Development (twistedpairdevelopment.wordpress.com).
# Pyrr is released under the MIT license (a very relaxed licence), but it is encouraged that any modifications are submitted back to the master for inclusion.
# Created by Adam Griffiths.
# Copyright (c) 2012, Twisted Pair Development.
# All rights reserved.
# twistedpairdevelopment.wordpress.com
# @twistedpairdev
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met: 
# 
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 
# 
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# 
# The views and conclusions contained in the software and documentation are those
# of the authors and should not be interpreted as representing official policies, 
# either expressed or implied, of the FreeBSD Project.

##############################################################################
# @dermen's random rotation matrix code
##############################################################################

from thor.math2 import rand_rot



##############################################################################
# random rotation matrix code adapted from Christoph Gohlke's transformation.py
# and Pyrr
##############################################################################

def random_rotation_matrix(rand=None):
    """Return uniform random rotation matrix.

    Parameters
    ----------
    rand: array like
        Three independent random variables that are uniformly distributed
        between 0 and 1 for each returned quaternion.

    Example
    ------
    >>> R = random_rotation_matrix()
    >>> np.allclose(np.dot(R.T, R), np.identity(3))
    True

    """
    return create_from_quaternion(random_quaternion(rand))


def random_quaternion(rand=None):
    """Return uniform random unit quaternion.

    Parameters
    ----------
    rand: array like or None
        Three independent random variables that are uniformly distributed
        between 0 and 1.

    Examples
    --------
    >>> q = random_quaternion()
    >>> np.allclose(1, vector_norm(q))
    True
    >>> q = random_quaternion(np.random.random(3))
    >>> len(q.shape), q.shape[0]==4
    (1, True)

    Notes
    -----
    This code is adapted from  Christoph Gohlke's transformation.py
        ..[2] http://www.lfd.uci.edu/~gohlke/code/transformations.py.html
    """
    if rand is None:
        rand = np.random.rand(3)
    else:
        assert len(rand) == 3
    r1 = np.sqrt(1.0 - rand[0])
    r2 = np.sqrt(rand[0])
    pi2 = np.pi * 2.0
    t1 = pi2 * rand[1]
    t2 = pi2 * rand[2]
    return np.array([np.cos(t2)*r2, np.sin(t1)*r1,
                        np.cos(t1)*r1, np.sin(t2)*r2])


def create_from_quaternion(quat):
    """Creates a matrix with the same rotation as a quaternion.

    Parameters
    ----------
    quat : np.ndarray, shape(4,),
        The quaternion to create the matrix from.

    Returns
    -------
    rotation : np.ndarray
        A matrix with shape (3,3) with the quaternion's rotation.

    Notes
    -----
    This code is adapted from pyrr
        ..[1] https://github.com/adamlwgriffiths/Pyrr/blob/master/pyrr/matrix33.py
    """
    x = quat[ 0 ]
    y = quat[ 1 ]
    z = quat[ 2 ]
    w = quat[ 3 ]

    y2 = y**2
    x2 = x**2
    z2 = z**2
    xy = x * y
    xz = x * z
    yz = y * z
    wx = w * x
    wy = w * y
    wz = w * z

    return np.array(
        [
            # m1
            [
                # m11 = 1.0 - 2.0 * (q.y * q.y + q.z * q.z)
                1.0 - 2.0 * (y2 + z2),
                # m12 = 2.0 * (q.x * q.y + q.w * q.z)
                2.0 * (xy + wz),
                # m13 = 2.0 * (q.x * q.z - q.w * q.y)
                2.0 * (xz - wy)
                ],
            # m2
            [
                # m21 = 2.0 * (q.x * q.y - q.w * q.z)
                2.0 * (xy - wz),
                # m22 = 1.0 - 2.0 * (q.x * q.x + q.z * q.z)
                1.0 - 2.0 * (x2 + z2),
                # m23 = 2.0 * (q.y * q.z + q.w * q.x)
                2.0 * (yz + wx)
                ],
            # m3
            [
                # m31 = 2.0 * (q.x * q.z + q.w * q.y)
                2.0 * (xz + wy),
                # m32 = 2.0 * (q.y * q.z - q.w * q.x)
                2.0 * (yz - wx),
                # m33 = 1.0 - 2.0 * (q.x * q.x + q.y * q.y)
                1.0 - 2.0 * (x2 + y2)
                ]
            ]
        )

Put these in a directory and run nosetests to check that rand_rot is passable

mender@texas:/data/work/mender/gist_rotate$ /data/ana/epd/bin/nosetests 
..
----------------------------------------------------------------------
Ran 2 tests in 1.073s

OK

[dermen]

Assuming the test code is flawless maybe we should stick with rand_rotate_molecule2 ? and implement in the cuda-code ?

[tjlane]

K sounds good to me. Are you up for doin' it? I can it will just be a while.

[dermen]

Im not so cuda-savy but if I find time I will look into it. Hasn't been a bottleneck so far, but I do want to learn some cuda..

Do we have an official cuda-dev-dude on board with thor?

[tjlane]

We are the official dudes :)

You don't really need to know any CUDA, though, the function to modify is this one:
https://github.com/tjlane/thor/blob/master/src/scatter/_gpuscatter.cu#L73

[dermen]

I implemented rand_rotate_molecule2 in cpuscatter to replace the quaternion method and created a factor of 2-3 reduction in speed! Counter-intuitive.. I suppose the python version uses optimized matrix multiplication and hence the speedup when comparing python methods.. here is the new code anyway, maybe we can speed it up using some optimized matrix functions...

void rotate(  float &rand1, float &rand2, float &rand3, float &r1, float &r2, float &r3, float &ax, float &ay, float &az )
{

/*
rand1,rand2,rand3 are 3 random numbers between 0-1
r1,r2,r3 are the atomic positions before rotating
a1,a2,a3 are the atomic positions after rotating

 SOURCE
   ------
 http://www.google.com/url?sa=t&rct=j&q=uniform%20random%20rotation%20matrix&source=web&cd=5&ved=0CE8QFjAE&url=http%3A%2F%2Fciteseerx.ist.psu.edu%2Fviewdoc%2Fdownload%3Fdoi%3D10.1.1.53.1357%26rep%3Drep1%26type%3Dps&ei=lw2cUa2eIMKRiQKXnIDwBQ&usg=AFQjCNHViFXwwa8kv_tobzteWYM8EaKF-w&sig2=148RpesMoZvJmtse2oerjg&bvm=bv.46751780,d.cGE
*/

float x1 = rand1*PI*2;
float x2 = rand2*PI*2;
float x3 = rand3;

// z-rotation
float z11 = cos(x1);
float z12 = sin(x1);
float z13 = 0;
float z21 = -sin(x1);
float z22 = cos(x1);
float z23 = 0;
float z31 = 0;
float z32 = 0;
float z33 = 1;

// vector to move the pole
float v1 = cos(x2)*sqrt(x3);
float v2 = sin(x2)*sqrt(x3);
float v3 = sqrt(1-x3);

// identity minus the outer product of vector
float v11 = 1-2*v1*v1; 
float v12 = -2*v1*v2; 
float v13 = -2*v1*v3; 
float v21 = -2*v2*v1; 
float v22 = 1-2*v2*v2; 
float v23 = -2*v2*v3; 
float v31 = -2*v3*v1; 
float v32 = -2*v3*v2; 
float v33 = 1-2*v3*v3; 

// matrix product
float r11 = -v11*z11 - v12*z21 - v13*z31;
float r12 = -v11*z12 - v12*z22 - v13*z32;
float r13 = -v11*z13 - v12*z23 - v13*z33;
float r21 = -v21*z11 - v22*z21 - v23*z31;
float r22 = -v21*z12 - v22*z22 - v23*z32;
float r23 = -v21*z13 - v22*z23 - v23*z33;
float r31 = -v31*z11 - v32*z21 - v33*z31;
float r32 = -v31*z12 - v32*z22 - v33*z32;
float r33 = -v31*z13 - v32*z23 - v33*z33;

// rotate the atom
ax = r11*r1 + r12*r2 + r13*r3;
ay = r21*r1 + r22*r2 + r23*r3;
az = r31*r1 + r32*r2 + r33*r3;
}

[tjlane]

you can always call BLAS -- not sure how much faster that will be

[tjlane]

I resolved the original issue taking @dermen 's advice above. I think we can move the discussion of the CUDA/C++ scattering code to issue #5 since it seems to be different from the python version & pertains directly to the performance of the scattering code rather than the public python interface.