Skip to content

Commit

Permalink
Formatted code using black
Browse files Browse the repository at this point in the history
  • Loading branch information
LaurentRDC committed Nov 9, 2018
1 parent c4ba02f commit 8d43aac
Show file tree
Hide file tree
Showing 18 changed files with 12,898 additions and 11,450 deletions.
6 changes: 4 additions & 2 deletions CONTRIBUTING.rst
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,11 @@ To set up `crystals` for local development:

Now you can make your changes locally.

4. When you're done making changes, run all the checks, doc builder and spell checker with the `unittest` standard module::
4. When you're done making changes, build the package components as follows::

python -m unittest discover
python -m unittest discover # run all tests
black .\crystals # run code formatter
python setup.py build_sphinx # build documentation

5. Commit your changes and push your branch to GitHub::

Expand Down
4 changes: 4 additions & 0 deletions README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@ crystals
:alt: Supported Python versions
:target: https://pypi.python.org/pypi/crystals

.. image:: https://img.shields.io/badge/code%20style-black-000000.svg
:alt: Code formatting style
:target: https://github.com/ambv/black

``crystals`` is a library of data structure and algorithms to manipulate abstract crystals in a Pythonic way. ``crystals`` helps with reading crystallographic
files (like .cif and .pdb), provides access to atomic positions, scattering utilities, and allows for symmetry determination. Although ``crystals`` can be used on its own,
it was made to be integrated into larger projects (like `scikit-ued <https://github.com/LaurentRDC/scikit-ued>`_).
Expand Down
10 changes: 5 additions & 5 deletions crystals/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@
"""
This package allows for manipulation and modelling of atomic structures in crystalline form.
"""
__author__ = 'Laurent P. René de Cotret'
__email__ = 'laurent.renedecotret@mail.mcgill.ca'
__license__ = 'BSD3'
__version__ = '1.0.0dev'
__author__ = "Laurent P. René de Cotret"
__email__ = "laurent.renedecotret@mail.mcgill.ca"
__license__ = "BSD3"
__version__ = "1.0.0dev"

from .atom import Atom
from .atom import frac_coords
from .atom import real_coords
Expand Down
180 changes: 101 additions & 79 deletions crystals/affine.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,9 @@

import numpy as np

#standard basis
e1,e2,e3 = np.eye(3)
# standard basis
e1, e2, e3 = np.eye(3)


def affine_map(array):
"""
Expand All @@ -30,18 +31,21 @@ def affine_map(array):
------
ValueError : If the transformation matrix is neither 3x3 or 4x4
"""
if array.shape == (4,4): #Already the right shape
if array.shape == (4, 4): # Already the right shape
return array
elif array.shape == (3,3):
extended_matrix = np.zeros(shape = (4,4), dtype = array.dtype)
elif array.shape == (3, 3):
extended_matrix = np.zeros(shape=(4, 4), dtype=array.dtype)
extended_matrix[-1, -1] = 1
extended_matrix[:3,:3] = array
extended_matrix[:3, :3] = array
return extended_matrix
else:
raise ValueError('Array shape not 3x3 or 4x4, and thus is not a transformation matrix.')
raise ValueError(
"Array shape not 3x3 or 4x4, and thus is not a transformation matrix."
)


def transform(matrix, array):
"""
"""
Applies a matrix transform on an array.
Parameters
Expand All @@ -61,19 +65,24 @@ def transform(matrix, array):
------
ValueError : If the transformation matrix is neither 3x3 or 4x4
"""
if matrix.shape not in [(3,3), (4,4)]:
raise ValueError('Input matrix is neither a 3x3 or 4x4 matrix, but \
rather of shape {}.'.format(matrix.shape))

matrix = affine_map(matrix)
# Case of a vector (e.g. position vector):
if array.ndim == 1:
extended_vector = np.array([0,0,0,1], dtype = array.dtype)
extended_vector[:3] = array
return np.dot(matrix, extended_vector)[:3]
else:
array = affine_map(array)
return np.dot(matrix, array)
if matrix.shape not in [(3, 3), (4, 4)]:
raise ValueError(
"Input matrix is neither a 3x3 or 4x4 matrix, but \
rather of shape {}.".format(
matrix.shape
)
)

matrix = affine_map(matrix)
# Case of a vector (e.g. position vector):
if array.ndim == 1:
extended_vector = np.array([0, 0, 0, 1], dtype=array.dtype)
extended_vector[:3] = array
return np.dot(matrix, extended_vector)[:3]
else:
array = affine_map(array)
return np.dot(matrix, array)


def translation_matrix(direction):
"""
Expand All @@ -92,8 +101,9 @@ def translation_matrix(direction):
matrix[:3, 3] = np.asarray(direction)[:3]
return matrix

def change_of_basis(basis1, basis2 = (e1,e2,e3)):
"""

def change_of_basis(basis1, basis2=(e1, e2, e3)):
"""
Returns the matrix that goes from one basis to the other.
Parameters
Expand All @@ -109,18 +119,19 @@ def change_of_basis(basis1, basis2 = (e1,e2,e3)):
Change-of-basis matrix that, applied to `basis`, will
return `basis2`.
"""
# Calculate the transform that goes from basis 1 to standard basis
basis1 = [np.asarray(vector).reshape(3,1) for vector in basis1]
basis1_to_standard = np.hstack(tuple(basis1))
# Calculate the transform that goes from basis 1 to standard basis
basis1 = [np.asarray(vector).reshape(3, 1) for vector in basis1]
basis1_to_standard = np.hstack(tuple(basis1))

# Calculate the transform that goes from standard basis to basis2
basis2 = [np.asarray(vector).reshape(3,1) for vector in basis2]
standard_to_basis2 = np.linalg.inv(np.hstack(tuple(basis2)))
# Calculate the transform that goes from standard basis to basis2
basis2 = [np.asarray(vector).reshape(3, 1) for vector in basis2]
standard_to_basis2 = np.linalg.inv(np.hstack(tuple(basis2)))

return np.dot(standard_to_basis2, basis1_to_standard)

return np.dot(standard_to_basis2, basis1_to_standard)

def is_basis(basis):
"""
"""
Returns true if the set of vectors forms a basis. This is done by checking
whether basis vectors are independent via an eigenvalue calculation.
Expand All @@ -133,7 +144,8 @@ def is_basis(basis):
out : bool
Whether or not the basis is valid.
"""
return (0 not in np.linalg.eigvals(np.asarray(basis)))
return 0 not in np.linalg.eigvals(np.asarray(basis))


def is_rotation_matrix(matrix):
"""
Expand All @@ -151,17 +163,18 @@ def is_rotation_matrix(matrix):
result : bool
If True, input could be a rotation matrix.
"""
# TODO: is this necessary? should a composite transformation
# of translation and rotation return True?
#if matrix.shape == (4,4):
# TODO: is this necessary? should a composite transformation
# of translation and rotation return True?
# if matrix.shape == (4,4):
# matrix = matrix[:3,:3]

is_orthogonal = np.allclose(np.linalg.inv(matrix), np.transpose(matrix))
unit_determinant = np.allclose(abs(np.linalg.det(matrix)), 1)
return is_orthogonal and unit_determinant

def rotation_matrix(angle, axis = (0,0,1)):
"""

def rotation_matrix(angle, axis=(0, 0, 1)):
"""
Return matrix to rotate about axis defined by direction around the origin [0,0,0].
To combine rotation and translations, see http://www.euclideanspace.com/maths/geometry/affine/matrix4x4/index.htm
Expand All @@ -181,24 +194,29 @@ def rotation_matrix(angle, axis = (0,0,1)):
--------
translation_rotation_matrix
"""
sina, cosa = math.sin(angle), math.cos(angle)

# Make sure direction is a numpy vector of unit length
direction = np.asarray(axis)
direction = direction/np.linalg.norm(direction)

# rotation matrix around unit vector
R = np.diag([cosa, cosa, cosa])
R += np.outer(direction, direction) * (1.0 - cosa)
direction *= sina
R += np.array([[ 0.0, -direction[2], direction[1]],
[ direction[2], 0.0, -direction[0]],
[-direction[1], direction[0], 0.0]])

return R
sina, cosa = math.sin(angle), math.cos(angle)

# Make sure direction is a numpy vector of unit length
direction = np.asarray(axis)
direction = direction / np.linalg.norm(direction)

# rotation matrix around unit vector
R = np.diag([cosa, cosa, cosa])
R += np.outer(direction, direction) * (1.0 - cosa)
direction *= sina
R += np.array(
[
[0.0, -direction[2], direction[1]],
[direction[2], 0.0, -direction[0]],
[-direction[1], direction[0], 0.0],
]
)

return R


def translation_rotation_matrix(angle, axis, translation):
"""
"""
Returns a 4x4 matrix that includes a rotation and a translation.
Parameters
Expand All @@ -215,9 +233,10 @@ def translation_rotation_matrix(angle, axis, translation):
matrix : `~numpy.ndarray`, shape (4,4)
Affine transform matrix.
"""
rmat = affine_map(rotation_matrix(angle = angle, axis = axis))
rmat[:3, 3] = np.asarray(translation)
return rmat
rmat = affine_map(rotation_matrix(angle=angle, axis=axis))
rmat[:3, 3] = np.asarray(translation)
return rmat


def change_basis_mesh(xx, yy, zz, basis1, basis2):
"""
Expand All @@ -237,21 +256,24 @@ def change_basis_mesh(xx, yy, zz, basis1, basis2):
XX, YY, ZZ : `~numpy.ndarray`
"""
# Build coordinate array row-wise
changed = np.empty(shape = (3, xx.size), dtype = np.float)
linearized = np.empty(shape = (3, xx.size), dtype = np.float)
linearized[0,:] = xx.ravel()
linearized[1,:] = yy.ravel()
linearized[2,:] = zz.ravel()
changed = np.empty(shape=(3, xx.size), dtype=np.float)
linearized = np.empty(shape=(3, xx.size), dtype=np.float)
linearized[0, :] = xx.ravel()
linearized[1, :] = yy.ravel()
linearized[2, :] = zz.ravel()

# Change the basis at each row
COB = change_of_basis(basis1, basis2)
np.dot(COB, linearized, out = changed)
return (changed[0,:].reshape(xx.shape),
changed[1,:].reshape(yy.shape),
changed[2,:].reshape(zz.shape))
np.dot(COB, linearized, out=changed)
return (
changed[0, :].reshape(xx.shape),
changed[1, :].reshape(yy.shape),
changed[2, :].reshape(zz.shape),
)


def minimum_image_distance(xx, yy, zz, lattice):
"""
"""
Returns a periodic array according to the minimum image convention.
Parameters
Expand All @@ -266,17 +288,17 @@ def minimum_image_distance(xx, yy, zz, lattice):
r : `~numpy.ndarray`
Minimum image distance over the lattice
"""
COB = change_of_basis(np.eye(3), lattice)
linearized = np.empty(shape = (3, xx.size), dtype = np.float) # In the standard basis
ulinearized = np.empty_like(linearized) # In the unitcell basis
COB = change_of_basis(np.eye(3), lattice)
linearized = np.empty(shape=(3, xx.size), dtype=np.float) # In the standard basis
ulinearized = np.empty_like(linearized) # In the unitcell basis

linearized[0,:] = xx.ravel()
linearized[1,:] = yy.ravel()
linearized[2,:] = zz.ravel()
linearized[0, :] = xx.ravel()
linearized[1, :] = yy.ravel()
linearized[2, :] = zz.ravel()

# Go to unitcell basis, where the cell is cubic of side length 1
np.dot(COB, linearized, out = ulinearized)
ulinearized -= np.rint(ulinearized)
np.dot(np.linalg.inv(COB), ulinearized, out = linearized)
# Go to unitcell basis, where the cell is cubic of side length 1
np.dot(COB, linearized, out=ulinearized)
ulinearized -= np.rint(ulinearized)
np.dot(np.linalg.inv(COB), ulinearized, out=linearized)

return np.reshape(np.linalg.norm(linearized, axis = 0), xx.shape)
return np.reshape(np.linalg.norm(linearized, axis=0), xx.shape)

0 comments on commit 8d43aac

Please sign in to comment.