/usr/lib/python2.7/dist-packages/cogent/maths/geometry.py is in python-cogent 1.9-9.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 | #!/usr/bin/env python
"""Code for geometric operations, e.g. distances and center of mass."""
from __future__ import division
from numpy import array, take, sum, newaxis, sqrt, sqrt, sin, cos, pi, c_, \
vstack, dot, ones
__author__ = "Sandra Smit"
__copyright__ = "Copyright 2007-2016, The Cogent Project"
__credits__ = ["Sandra Smit", "Gavin Huttley", "Rob Knight", "Daniel McDonald",
"Marcin Cieslik"]
__license__ = "GPL"
__version__ = "1.9"
__maintainer__ = "Sandra Smit"
__email__ = "sandra.smit@colorado.edu"
__status__ = "Production"
def center_of_mass(coordinates, weights= -1):
"""Calculates the center of mass for a dataset.
coordinates, weights can be two things:
either: coordinates = array of coordinates, where one column contains
weights, weights = index of column that contains the weights
or: coordinates = array of coordinates, weights = array of weights
weights = -1 by default, because the simplest case is one dataset, where
the last column contains the weights.
If weights is given as a vector, it can be passed in as row or column.
"""
if isinstance(weights, int):
return center_of_mass_one_array(coordinates, weights)
else:
return center_of_mass_two_array(coordinates, weights)
def center_of_mass_one_array(data, weight_idx= -1):
"""Calculates the center of mass for a dataset
data should be an array of x1,...,xn,r coordinates, where r is the
weight of the point
"""
data = array(data)
coord_idx = range(data.shape[1])
del coord_idx[weight_idx]
coordinates = take(data, (coord_idx), 1)
weights = take(data, (weight_idx,), 1)
return sum(coordinates * weights, 0) / sum(weights, 0)
def center_of_mass_two_array(coordinates, weights):
"""Calculates the center of mass for a set of weighted coordinates
coordinates should be an array of coordinates
weights should be an array of weights. Should have same number of items
as the coordinates. Can be either row or column.
"""
coordinates = array(coordinates)
weights = array(weights)
try:
return sum(coordinates * weights, 0) / sum(weights, 0)
except ValueError:
weights = weights[:, newaxis]
return sum(coordinates * weights, 0) / sum(weights, 0)
def distance(first, second):
"""Calculates Euclideas distance between two vectors (or arrays).
WARNING: Vectors have to be the same dimension.
"""
return sqrt(sum(((first - second) ** 2).ravel()))
def sphere_points(n):
"""Calculates uniformly distributed points on a unit sphere using the
Golden Section Spiral algorithm.
Arguments:
-n: number of points
"""
points = []
inc = pi * (3 - sqrt(5))
offset = 2 / float(n)
for k in xrange(int(n)):
y = k * offset - 1 + (offset / 2)
r = sqrt(1 - y * y)
phi = k * inc
points.append([cos(phi) * r, y, sin(phi) * r])
return array(points)
def coords_to_symmetry(coords, fmx, omx, mxs, mode):
"""Applies symmetry transformation matrices on coordinates. This is used to
create a crystallographic unit cell or a biological molecule, requires
orthogonal coordinates, a fractionalization matrix (fmx),
an orthogonalization matrix (omx) and rotation matrices (mxs).
Returns all coordinates with included identity, which should be the first
matrix in mxs.
Arguments:
- coords: an array of orthogonal coordinates
- fmx: fractionalization matrix
- omx: orthogonalization matrix
- mxs: a sequence of 4x4 rotation matrices
- mode: if mode 'table' assumes rotation matrices operate on
fractional coordinates (like in crystallographic tables).
"""
all_coords = [coords] # the first matrix is identity
if mode == 'fractional': # working with fractional matrices
coords = dot(coords, fmx.transpose())
# add column of 1.
coords4 = c_[coords, array([ones(len(coords))]).transpose()]
for i in xrange(1, len(mxs)): # skip identity
rot_mx = mxs[i].transpose()
new_coords = dot(coords4, rot_mx)[:, :3] # rotate and translate, remove
if mode == 'fractional': # ones column
new_coords = dot(new_coords, omx.transpose()) # return to orthogonal
all_coords.append(new_coords)
# a vstack(arrays) with a following reshape is faster then
# the equivalent creation of a new array via array(arrays).
return vstack(all_coords).reshape((len(all_coords), coords.shape[0], 3))
def coords_to_crystal(coords, fmx, omx, n=1):
"""Applies primitive lattice translations to produce a crystal from the
contents of a unit cell.
Returns all coordinates with included zero translation (0, 0, 0).
Arguments:
- coords: an array of orthogonal coordinates
- fmx: fractionalization matrix
- omx: orthogonalization matrix
- n: number of layers of unit-cells == (2*n+1)^2 unit-cells
"""
rng = range(-n, n + 1) # a range like -2, -1, 0, 1, 2
fcoords = dot(coords, fmx.transpose()) # fractionalize
vectors = [(x, y, z) for x in rng for y in rng for z in rng]
# looking for the center Thickened cube numbers:
# a(n)=n*(n^2+(n-1)^2)+(n-1)*2*n*(n-1) ;)
all_coords = []
for primitive_vector in vectors:
all_coords.append(fcoords + primitive_vector)
# a vstack(arrays) with a following reshape is faster then
# the equivalent creation of a new array via array(arrays)
all_coords = vstack(all_coords).reshape((len(all_coords), \
coords.shape[0], coords.shape[1], 3))
all_coords = dot(all_coords, omx.transpose()) # orthogonalize
return all_coords
|