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.

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

ENH: spatial.transform.Rotation [GSOC2018] #8945

Merged
merged 105 commits into from Jun 30, 2018

Conversation

Projects
None yet
5 participants
@adbugger
Copy link
Contributor

adbugger commented Jun 19, 2018

This pull request implements, documents, and tests the following features to the Rotation class under a new sub-package scipy.spatial.transform:

  • Initializers from and converters to quaternions, direction cosine matrices, euler angles and rotation vector representations of rotations
  • Application of a rotation to vectors
  • Inversion and composition of rotations
  • Indexing into a Rotation object to extract specific rotation(s)
@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 20, 2018

Mentioning @nmayorov , @larsoner since I can't add reviewers.

One test case related to Euler angles is failing by returning nan instead of a value. Curiously, all tests pass on my local setup and on TravisCI and the failure only occurs on appveyor. Since it doesn't appear to be a platform issue, how do I start fixing it?

@ilayn

This comment has been minimized.

Copy link
Member

ilayn commented Jun 20, 2018

Sometimes we get a one-off error for no apparent reason. I've restarted the Appveyor tests and added the reviewers. Let's see how it goes.

@ilayn ilayn added this to the 1.2.0 milestone Jun 20, 2018

@larsoner
Copy link
Member

larsoner left a comment

The failure looks real and specific to 3.6 on Windows (x64):

https://ci.appveyor.com/project/scipy/scipy/build/1.0.3155/job/g6hk4ln5ne249xew/tests

Knowing the input values, any ideas why this could happen? Are some numbers close to numerical limits that would produce nan (e.g., division by zero)?

"""
from __future__ import division, print_function, absolute_import

from .rotation import *

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

Better to be explicit and do from .rotation import Rotation


from .rotation import *

__all__ = [s for s in dir() if not s.startswith('_')]

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

Here you might as well just do __all__ = ['Rotation']. The * import above will import many things without leading _ that we do not want in __all__, such as np, re, etc. And even if you change that import, then division, print_function, absolute_import will still be in __all__ (and the code here will have the problem that if someone else does such a * import in the future they will have this problem).

Although the explicit list has some redundancy, it does not suffer from these problems.

@property
def n(self):
"""Number of rotations contained in object."""
return self._quat.shape[0]

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

I would just make this def __len__ and users do len(rots) instead of rots.n

if self._single:
return self._quat[0]
else:
return self._quat

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

This should probably return a .copy() of these quantities, and state that a copy is returned. Otherwise users could accidentally change the returned arrays in-place, and the Rotation object will be corrupted.

This comment has been minimized.

@adbugger

adbugger Jun 20, 2018

Author Contributor

Given that users do not (should not?) know how we store the rotation internally, should we say that a copy is being returned?

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

In that case, no don't mention it, but do the .copy() here

"""Initialize class from rotation vector.
A rotation vector is a 3 dimensional vector which is co-directional to
the axis of rotaion and whose norm gives the angle of rotation (in

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

rotaion -> rotation

result = result[0]
return self.__class__(result, normalized=True, copy=False)

def apply(self, vectors, inverse=False):

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

I'm not sure inverse is needed. One can just as easily do rots.inv().apply(vecs) if they want this behavior

This comment has been minimized.

@nmayorov

nmayorov Jun 23, 2018

Contributor

How do we resolve this suggestions? I'm slightly in favor for keeping inverse, because this is how I'm used to think about rotation application --- one objects does rotations in both directions. In other words, we don't need to have a rotation and its inverse to do the job.

can either be specified with shape `(3, )` or `(1, 3)`.
inverse : boolean, optional
If `inverse` is `True` then the inverse of the rotation(s) is
applied to the input vectors. Default is `False`.

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

Missing Returns section

"""Returns the inverse of the current rotation.
This function returns a new `Rotation` instance containing the inverse
of all the rotations in the current instance.

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

Missing Returns section. Most methods appear to be missing it

Specifies which rotation(s) to extract.
"""
# __init__ now copies by default
return self.__class__(self._quat[indexer], normalized=True)

This comment has been minimized.

@larsoner

larsoner Jun 20, 2018

Member

this should have a .copy() because if indexer makes a view (e.g., is a slice) then the new object will not own its own data, i.e.:

>>> x = np.zeros(3)
>>> y = np.asarray(x, float)
>>> y[2] = 1
>>> x
array([0., 0., 1.])

It would be good to add a tiny test case that exposes this problem, and then fix it (true TDD!)

This comment has been minimized.

@adbugger

adbugger Jun 20, 2018

Author Contributor

We've already taken care of this in the __init__ function. Even if self._quat[indexer] returns a view (which it does in most cases if my understanding of python indexing and slicing is correct), the self.__class__() calls __init__ where we make a copy of the input quaternions before setting the private _quat variable.

This has been verified by the following test case, which I will include in the test_rotation file:

>>> r = R.from_quaternion(np.random.uniform(size=(3,4)))
>>> r._quat
array([[0.61921279, 0.49144496, 0.592531  , 0.15480437],
       [0.47882848, 0.60457664, 0.24404027, 0.58792407],
       [0.10501252, 0.1217286 , 0.25670971, 0.95302395]])
>>> s = r[0:2]
>>> s._quat
array([[0.61921279, 0.49144496, 0.592531  , 0.15480437],
       [0.47882848, 0.60457664, 0.24404027, 0.58792407]])
>>> r._quat[0] = [0,0,0,1]
>>> r._quat
array([[0.        , 0.        , 0.        , 1.        ],
       [0.47882848, 0.60457664, 0.24404027, 0.58792407],
       [0.10501252, 0.1217286 , 0.25670971, 0.95302395]])
>>> s._quat
array([[0.61921279, 0.49144496, 0.592531  , 0.15480437],
       [0.47882848, 0.60457664, 0.24404027, 0.58792407]])

The new s object still contains its initial quaternions and is not modified when r is.

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 20, 2018

The following changes have been done:

  • Overall documentation and style changes
  • Better ordering of methods in the "methods" section
  • Ensured that the user cannot accidentally corrupt the object (Added a test case specifically to check for this)
  • Addition of a __len__ function and documentation for it.

Will figure out the test failure next.


# Step 4
angles = np.empty((num_rotations, 3))
angles[:, 1] = np.arccos(dcm_transformed[:, 2, 2])

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

This place causes nans, sure about 99%. You need to ensure that dcm_transformed[:, 2, 2] doesn't exceed 1 by absolute value. Or you could assign these angle to 0 or pi depending on the sign. Probably the first approach is cleaner.

@ilayn

This comment has been minimized.

Copy link
Member

ilayn commented Jun 21, 2018

You can also use rebasing and squashing your local commits if they are of similar nature. For reference https://help.github.com/articles/about-git-rebase/

@nmayorov
Copy link
Contributor

nmayorov left a comment

Couldn't check the whole code, but you should get main ideas.

Rotations in 3 dimensions
-------------------------
This package supports `Rotation` transforms with the following representations

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

I suggest to remove this section all together for now. So keep it as

Rotations in 3 dimensions
-------------------------
.. autosummary::
      ...

import numpy as np
import scipy.linalg
import re

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Perhaps put re and warning first.

import warnings


AXIS_TO_IND = {'x': 0, 'y': 1, 'z': 2}

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Use underscore for consistency?

This comment has been minimized.

@adbugger

adbugger Jun 21, 2018

Author Contributor

We removed underscore in favor of upper case names in a previous suggestion. Should I change it back?

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

I suggest underscore and upper case as a private constant.


def _make_elementary_quat(axis, angles):
num_rotations = angles.shape[0]
quat = np.zeros((num_rotations, 4))

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Perhaps just np.zeros((angle.shape[0], 4))?



class Rotation(object):
"""Rotation transform in 3 dimensions.

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

I guess it is usually said as "Rotation transformation" or you can just put "Rotation in 3 dimensions".

Returns
-------
output : `numpy.ndarray`, shape (4,) or (N, 4)

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

I think I would prefer to have more specific names.

Rotations in 3 dimensions can be represented using unit norm
quaternions [1]_. The mapping from quaternions to rotations is
two-to-one, i.e. a quaternion and its additive inverse represent the

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Additive inverse? Write something simpler, like "If you change the sign of each component you will get the same rotation." (But better than that, couldn't formulate it quickly).

else:
return self._quat.copy()

def as_dcm(self):

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

I think it makes sense to order methods exactly like in "Methods" section.

def as_dcm(self):
"""Represent rotations as direction cosine matrices.
3D rotations can be represented using direction cosine matrices, which

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Some reference is useful.

def from_rotvec(cls, rotvec):
"""Initialize class from rotation vector.
A rotation vector is a 3 dimensional vector which is co-directional to

This comment has been minimized.

@nmayorov

nmayorov Jun 21, 2018

Contributor

Reference is useful.

# Ensure less than unit norm
positive_unity = dcm_transformed[:, 2, 2] > 1
negative_unity = dcm_transformed[:, 2, 2] < -1
dcm_transformed[positive_unity, 2, 2] = 1.0

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

I think it makes sense to write it as 1 or 1.0 everywhere, I prefer 1.

@nmayorov nmayorov changed the title Initial Pull Request for Rotation class ENH: spatial.transform.Rotation [GSOC2018] Jun 22, 2018

@nmayorov
Copy link
Contributor

nmayorov left a comment

Besides some comments I made, I don't see any other significant problems. So consider them and after that we can wait for couple of days and merge.

Also go through everything again and fix any flaws you will find.

References
----------
.. [1] F. Landis Markley, `Unit Quaternion from Rotation Matrix

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

Super minor, but maybe reorder [1] and [2].

is_single = False
dcm = np.asarray(dcm, dtype=float)

if dcm.ndim not in [2, 3] or (dcm.shape[-2], dcm.shape[-1]) != (3, 3):

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

dcm.shape[-2:]

raise ValueError("Expected axis specification to be a non-empty "
"string of upto 3 characters, got {}".format(seq))

intrinsic = (re.compile('^[XYZ]{1,3}$').match(seq) is not None)

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

I remember Eric commented that you can do it simpler (without compile). Can you check that?

return cls(quat[0] if is_single else quat, normalized=True, copy=False)

def as_quat(self):
"""Represent rotations as quaternions.

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

I suggest to consider whether we use plural or singular forms in similar places. Probably the plural form is more natural.

"""Represent rotations as direction cosine matrices.
3D rotations can be represented using direction cosine matrices, which
are 3 x 3 real orthogonal matrices with eigenvalue equal to +1 [1]_.

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

You meant determinant?

rotation `p[i]` is composed with the corresponding rotation
`q[i]` and `output` contains `N` rotations.
"""
if not(self._quat.shape[0] == 1 or other._quat.shape[0] == 1 or

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

Can use len here. Check other places too.

"""Invert this rotation.
Composition of a rotation with its inverse results in an identity
rotation, or no rotation.

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

In an "identity transformation"?

Returns
-------
rotation : `Rotation` instance

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

inverse?

Returns
-------
rotation : `Rotation` instance
`output` contains

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

Just "Contains"?

Indexing within a rotation is supported since multiple rotation transforms
can be stored within a single `Rotation` instance.
Initialization using the `from_...` classmethods such as `from_quat` is

This comment has been minimized.

@nmayorov

nmayorov Jun 22, 2018

Contributor

I think we should say To create `Rotation` objects use `from_...` classmethods, `__init__` is not supposed to be used directly.


@classmethod
def from_quat(cls, quat, normalized=False):
"""Initialize Rotation from quaternions.

This comment has been minimized.

@nmayorov

nmayorov Jun 23, 2018

Contributor

I think you can omit "Rotation" in similar places. That is "Initialize from quaternions".

return cls(quat[0] if is_single else quat, normalized=True, copy=False)

def as_quat(self):
"""Represent rotations as quaternions.

This comment has been minimized.

@nmayorov

nmayorov Jun 23, 2018

Contributor

Maybe shorten it to "Represent as quanterions" (and in all similar places).

return rotvec

def as_euler(self, seq, degrees=False):
"""Compute Euler angles for rotations with specified axis sequence.

This comment has been minimized.

@nmayorov

nmayorov Jun 23, 2018

Contributor

Maybe follow the format you used in other as_ methods, i.e. "Represent [rotations, Rotation] as Euler angles."

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Jun 23, 2018

Reviews from everyone else would be appreciated. You could check the overall language and style of documentation.

I plan to merge it on Tuesday or Wednesday. I don't think it is too soon, because Aditya continues to work on the project and will improve small things in Rotation class as he progresses.

@rgommers
Copy link
Member

rgommers left a comment

Overall it looks good to me from a read through (didn't test). A couple of comments:

  • this is a new public submodule, so should be added to doc/API.rst.txt
  • Examples sections in the class docstring and method docstrings are still missing. It would be good to have an example for each of the basic operations, and also for things like storing and working with multiple rotations in a single Rotation instance.
@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 25, 2018

Added examples to each method depicting usage with single and multiple rotations. I think I've covered everything in the methods so I wasn't sure what to add to the class docstring. Take a look and let me know what seems appropriate.

@rgommers
Copy link
Member

rgommers left a comment

Examples LGTM, thanks @adbugger. Just a couple of minor comments

Examples
--------
>>> import numpy as np

This comment has been minimized.

@rgommers

rgommers Jun 25, 2018

Member

this import is never necessary in examples, it's always implied

>>> r.as_dcm().shape
(2, 3, 3)
If input matrices are not special orthogonal, then a special orthogonal

This comment has been minimized.

@rgommers

rgommers Jun 25, 2018

Member

I know what orthogonal is, but "special orthogonal" may deserve an explanation. what is it?

>>> np.linalg.det(dcm)
1.0000000000000002
It is also possible to have a stack containng a single rotation:

This comment has been minimized.

@rgommers

rgommers Jun 25, 2018

Member

typo in containing

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Jun 25, 2018

My thought about examples.

  1. They show that our design is quite robust and versatile. I like it.
  2. I think that you focused too much on formal things, like that we can initialize from a single object or stack. I think this is not so hard of a concept to understand. Instead I would prefer to see more insightful examples when user can say "yes it looks familiar, it makes sense".
  3. Maybe we can write a single example section for the whole class. A scheme I had in mind: initialize a single rotation from each representation --- I suggest to use a rotation around a single axis (use different ones, identity it also fine to start with) and show conversion to all other formats. This way a user could gain an insights how different representations are related. Mention that we can use stack initialization as well. Then create a stack rotation(s) and show how we can apply, compose, inverse and index them.

I understand that the current approach has its advantages as well. And this is just my suggestion to consider, I think it can be more interesting.

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 26, 2018

I have added an examples section in the class docstring. I think the class documentation serves as a good overview into the class and I suggest that we let the individual methods' examples remain as they are. They serve as a formal list of all possible input and output formats each method supports and are a nice addition.

A `Rotation` instance can be indexed and sliced as if it were a single
1D array or list:
>>> r = R.from_euler('zyx', [

This comment has been minimized.

@nmayorov

nmayorov Jun 26, 2018

Contributor

Looks like you can reuse the previously defined rotation.

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Jun 26, 2018

Really liked the new examples section. And totally don't mind keeping individual examples as well.

What happened with all these errors from some optimize code?

@larsoner

This comment has been minimized.

Copy link
Member

larsoner commented Jun 26, 2018

They are fixed in #8975 which needs review

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 29, 2018

It looks like the optimize failures have gone away. @nmayorov can you restart that one failed job on travis-ci? I want to see if it's a one-off worker failure or a real failure.

@ilayn

This comment has been minimized.

Copy link
Member

ilayn commented Jun 29, 2018

Done

adbugger added some commits Jun 18, 2018

STY: Changes for clairty
- Reorder imports
- Reorder functions to match order in methods section
- Remove unhelpful comments
- Shorten ..._quaternion to ..._quat.
DOC: Changed function docstrings.
More specific output names. Added references for rotvecs and dcms.
Simplified two-to-one mapping explanation in as_quat.
BUG: Fix nan error in degenerate euler axes
For certain gimbal locked cases, it appears that certain entries in
`dcm_transformed` are greater than unit norm. That results in nan values
when arccos function us used.

@adbugger adbugger force-pushed the adbugger:rotation branch from cf5158b to 7f8b9c9 Jun 30, 2018

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jun 30, 2018

@nmayorov , I have rebased and shortened commit history here. If this looks good, I'll start doing the same for the other PRs.

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Jun 30, 2018

I'm merging it now. It is in a decent state at least and we will improve it if we spot something during the end of GSoC.

Congratulations @adbugger for a solid contribution.

@nmayorov nmayorov merged commit 8835355 into scipy:master Jun 30, 2018

5 checks passed

ci/circleci: build_docs Your tests passed on CircleCI!
Details
ci/circleci: pypy3 Your tests passed on CircleCI!
Details
codecov/project 76.59% (target 1%)
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment