diff --git a/package/CHANGELOG b/package/CHANGELOG index a05f07222a2..d54a313609a 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -32,6 +32,9 @@ Enhancements * Added *Group.copy() methods returning an identical copy of the respective group (PR #1922) * Use a faster function to deduplicate indices (PR #1951) + * Calculations in *Group.center() are performed in double precision (#PR1936) + * Functions in lib.distances accept coordinate arrays of arbitrary dtype + (PR #1936) Fixes * rewind in the SingleFrameReader now reads the frame from the file (Issue #1929) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 816745d5a8b..68217167934 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -711,11 +711,19 @@ def center(self, weights, pbc=None, compound='group'): atoms = self.atoms + # enforce calculations in double precision: + dtype = np.float64 + if compound.lower() == 'group': if pbc: coords = atoms.pack_into_box(inplace=False) else: coords = atoms.positions + # promote coords or weights to dtype if required: + if weights is None: + coords = coords.astype(dtype, copy=False) + else: + weights = weights.astype(dtype, copy=False) return np.average(coords, weights=weights, axis=0) elif compound.lower() == 'residues': compound_indices = atoms.resindices @@ -728,14 +736,18 @@ def center(self, weights, pbc=None, compound='group'): " one of 'group', 'residues', or 'segments'." "".format(compound)) - # Sort positions and masses by compound index: + # Sort positions and weights by compound index and promote to dtype if + # required: sort_indices = np.argsort(compound_indices) compound_indices = compound_indices[sort_indices] coords = atoms.positions[sort_indices] - if weights is not None: + if weights is None: + coords = coords.astype(dtype, copy=False) + else: + weights = weights.astype(dtype, copy=False) weights = weights[sort_indices] # Allocate output array: - centers = np.zeros((n_compounds, 3), dtype=coords.dtype) + centers = np.zeros((n_compounds, 3), dtype=dtype) # Get sizes of compounds: unique_compound_indices, compound_sizes = np.unique(compound_indices, return_counts=True) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 11ba3c1e112..01bbeb868d0 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -227,27 +227,29 @@ def distance_array(reference, configuration, box=None, result=None, backend="ser Parameters ---------- - reference : numpy.array of numpy.float32 - Reference coordinate array. - configuration : numpy.array of numpy.float32 - Configuration coordinate array. - box : numpy.array or None + reference : numpy.ndarray + Reference coordinate array of shape ``(n, 3)`` (``dtype`` is arbitrary, + will be converted to ``dtype=numpy.float32`` internally) + configuration : numpy.ndarray + Configuration coordinate array of shape ``(m, 3)`` (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) + box : numpy.ndarray or None Dimensions of the cell; if provided, the minimum image convention is applied. The dimensions must be provided in the same format as returned by by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - result : numpy.array of numpy.float64, optional + result : numpy.ndarray(dtype=numpy.float64), optional Preallocated result array which must have the shape ``(len(ref), len(conf))`` and ``dtype=numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. [``None``] - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). Returns ------- - d : numpy.array + d : numpy.ndarray ``(len(reference),len(configuration))`` numpy array with the distances ``d[i,j]`` between reference coordinates `i` and configuration coordinates `j`. @@ -260,9 +262,11 @@ def distance_array(reference, configuration, box=None, result=None, backend="ser .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - ref = reference.copy('C') - conf = configuration.copy('C') + ref = reference.astype(np.float32, order='C', copy=True) + conf = configuration.astype(np.float32, order='C', copy=True) _check_array(conf, 'conf') _check_array(ref, 'ref') @@ -313,25 +317,25 @@ def self_distance_array(reference, box=None, result=None, backend="serial"): Parameters ---------- - reference : array - Reference coordinate array with ``N=len(ref)`` coordinates. - box : array or None + reference : numpy.ndarray + Reference coordinate array with ``N=len(ref)`` coordinates (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) + box : numpy.ndarray or None Dimensions of the cell; if provided, the minimum image convention is applied. The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - result : array, optional - Preallocated result array which must have the shape - ``(N*(N-1)/2,)`` and dtype ``numpy.float64``. Avoids creating - the array which saves time when the function is called - repeatedly. [``None``] - backend + result : numpy.ndarray(dtype=numpy.float64), optional + Preallocated result array which must have the shape ``(N*(N-1)/2,)`` and + dtype ``numpy.float64``. Avoids creating the array which saves time when + the function is called repeatedly. [``None``] + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). Returns ------- - d : array + d : numpy.ndarray ``N*(N-1)/2`` numpy 1D array with the distances dist[i,j] between ref coordinates i and j at position d[k]. Loop through d: @@ -345,13 +349,15 @@ def self_distance_array(reference, box=None, result=None, backend="serial"): Note ---- This method is slower than it could be because internally we need to make - copies of the coordinate arrays. + copies of the coordinate array. .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - ref = reference.copy('C') + ref = reference.astype(np.float32, order='C', copy=True) _check_array(ref, 'ref') @@ -400,13 +406,14 @@ def transform_RtoS(inputcoords, box, backend="serial"): Parameters ---------- inputcoords : array - A (3,) coordinate or (n x 3) array of coordinates, type ``np.float32``. + A (3,) coordinate or (n x 3) array of coordinates (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) box : array The unitcell dimesions for this system. The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). @@ -417,8 +424,10 @@ def transform_RtoS(inputcoords, box, backend="serial"): .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - coords = inputcoords.copy('C') + coords = inputcoords.astype(np.float32, order='C', copy=True) is_1d = False # True if only one vector coordinate if len(coords.shape) == 1: @@ -459,13 +468,14 @@ def transform_StoR(inputcoords, box, backend="serial"): Parameters ---------- inputcoords : array - A (3,) coordinate or (n x 3) array of coordinates, type ``np.float32``. + A (3,) coordinate or (n x 3) array of coordinates (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) box : array The unitcell dimesions for this system. The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). @@ -477,9 +487,10 @@ def transform_StoR(inputcoords, box, backend="serial"): .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - _check_array_dtype(inputcoords, 'S') - coords = inputcoords.copy('C') + coords = inputcoords.astype(np.float32, order='C', copy=True) is_1d = False # True if only one vector coordinate if len(coords.shape) == 1: @@ -527,9 +538,11 @@ def calc_bonds(coords1, coords2, box=None, result=None, backend="serial"): Parameters ---------- coords1 : array - An array of coordinates for one half of the bond. + An array of coordinates for one half of the bond (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) coords2 : array - An array of coordinates for the other half of bond + An array of coordinates for the other half of bond (``dtype`` is + arbitrary, will be converted to ``dtype=numpy.float32`` internally) box : array The unitcell dimesions for this system. The dimensions must be provided in the same format as returned @@ -539,7 +552,7 @@ def calc_bonds(coords1, coords2, box=None, result=None, backend="serial"): Preallocated result array which must be same length as coord arrays and ``dtype=numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. [None] - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). @@ -552,9 +565,11 @@ def calc_bonds(coords1, coords2, box=None, result=None, backend="serial"): .. versionadded:: 0.8 .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - atom1 = coords1.copy('C') - atom2 = coords2.copy('C') + atom1 = coords1.astype(np.float32, order='C', copy=True) + atom2 = coords2.astype(np.float32, order='C', copy=True) _check_array(atom1, 'atom1') _check_array(atom2, 'atom2') @@ -611,29 +626,32 @@ def calc_angles(coords1, coords2, coords3, box=None, result=None, backend="seria Parameters ---------- - coords1 : array - Coordinate array of one side of angles. - coords2 : array - Coordinate array of apex of angles. - coords3 : array - Coordinate array of other side of angles. - box : array + coords1 : numpy.ndarray + Coordinate array of one side of angles (``dtype`` is arbitrary, will be + converted to ``dtype=numpy.float32`` internally) + coords2 : numpy.ndarray + Coordinate array of apex of angles (``dtype`` is arbitrary, will be + converted to ``dtype=numpy.float32`` internally) + coords3 : numpy.ndarray + Coordinate array of other side of angles (``dtype`` is arbitrary, will be + converted to ``dtype=numpy.float32`` internally) + box : numpy.ndarray, optional The unitcell dimesions for this system. The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - result : array, optional + result : numpy.ndarray, optional Preallocated result array which must be same length as coord arrays and ``dtype=numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. [None] - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). Returns ------- - angles : array - A numpy.array of angles in radians. + angles : numpy.ndarray + An array of angles in radians. .. versionadded:: 0.8 @@ -641,10 +659,12 @@ def calc_angles(coords1, coords2, coords3, box=None, result=None, backend="seria Added optional box argument to account for periodic boundaries in calculation .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - atom1 = coords1.copy('C') - atom2 = coords2.copy('C') - atom3 = coords3.copy('C') + atom1 = coords1.astype(np.float32, order='C', copy=True) + atom2 = coords2.astype(np.float32, order='C', copy=True) + atom3 = coords3.astype(np.float32, order='C', copy=True) numatom = atom1.shape[0] _check_array(atom1, 'coords1') @@ -711,13 +731,17 @@ def calc_dihedrals(coords1, coords2, coords3, coords4, box=None, result=None, Parameters ---------- coords1 : array - Coordinate array of 1st atom in dihedrals. + Coordinate array of 1st atom in dihedrals (``dtype`` is arbitrary, will + be converted to ``dtype=numpy.float32`` internally) coords2 : array - Coordinate array of 2nd atom in dihedrals. + Coordinate array of 2nd atom in dihedrals (``dtype`` is arbitrary, will + be converted to ``dtype=numpy.float32`` internally) coords3 : array - Coordinate array of 3rd atom in dihedrals. + Coordinate array of 3rd atom in dihedrals (``dtype`` is arbitrary, will + be converted to ``dtype=numpy.float32`` internally) coords4 : array - Coordinate array of 4th atom in dihedrals. + Coordinate array of 4th atom in dihedrals (``dtype`` is arbitrary, will + be converted to ``dtype=numpy.float32`` internally) box : array The unitcell dimesions for this system. The dimensions must be provided in the same format as returned @@ -727,7 +751,7 @@ def calc_dihedrals(coords1, coords2, coords3, coords4, box=None, result=None, Preallocated result array which must be same length as coord arrays and ``dtype=numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. [None] - backend + backend : str Select the type of acceleration; "serial" is always available. Other possibilities are "OpenMP" (OpenMP). @@ -744,11 +768,13 @@ def calc_dihedrals(coords1, coords2, coords3, coords4, box=None, result=None, Renamed from calc_torsions to calc_dihedrals .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - atom1 = coords1.copy('C') - atom2 = coords2.copy('C') - atom3 = coords3.copy('C') - atom4 = coords4.copy('C') + atom1 = coords1.astype(np.float32, order='C', copy=True) + atom2 = coords2.astype(np.float32, order='C', copy=True) + atom3 = coords3.astype(np.float32, order='C', copy=True) + atom4 = coords4.astype(np.float32, order='C', copy=True) _check_array(atom1, 'atom1') _check_array(atom2, 'atom2') @@ -796,21 +822,22 @@ def apply_PBC(incoords, box, backend="serial"): Parameters ---------- - coords : array - Coordinate array (of type numpy.float32). + incoords : numpy.ndarray + Coordinate array of shape ``(n, 3)`` (``dtype`` is arbitrary, will be + converted to ``dtype=numpy.float32`` internally) box : array The unitcell dimesions for this system; can be either orthogonal or triclinic information. The dimensions must be provided in the same format as returned by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: ``[lx, ly, lz, alpha, beta, gamma]``. - backend - Select the type of acceleration; "serial" is always available. Other - possibilities are "OpenMP" (OpenMP). + backend : str + Select the type of acceleration; ``"serial"`` is always available. Other + possibilities are ``"OpenMP"`` (OpenMP). Returns ------- - newcoords : array + newcoords : numpy.ndarray(dtype=numpy.float32) Coordinates that are now all within the primary unit cell, as defined by box. @@ -818,8 +845,10 @@ def apply_PBC(incoords, box, backend="serial"): .. versionadded:: 0.8 .. versionchanged:: 0.13.0 Added *backend* keyword. + .. versionchanged:: 0.19.0 + Internal dtype conversion of input coordinates to ``numpy.float32``. """ - coords = incoords.copy('C') + coords = incoords.astype(np.float32, order='C', copy=True) _check_array(coords, 'coords') diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 71e74231239..9a6c5aff6c6 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -27,19 +27,13 @@ import MDAnalysis as mda -def test_transform_StoR_pass(): +@pytest.mark.parametrize('coord_dtype', (np.float32, np.float64)) +def test_transform_StoR_pass(coord_dtype): box = np.array([10, 7, 3, 45, 60, 90], dtype=np.float32) - s = np.array([[0.5, -0.1, 0.5]], dtype=np.float32) + s = np.array([[0.5, -0.1, 0.5]], dtype=coord_dtype) original_r = np.array([[ 5.75, 0.36066014, 0.75000012]], dtype=np.float32) test_r = mda.lib.distances.transform_StoR(s, box) assert_equal(original_r, test_r) - -def test_transform_StoR_fail(): - box = np.array([10, 7, 3, 45, 60, 90], dtype=np.float32) - s = np.array([[0.5, -0.1, 0.5]]) - - with pytest.raises(TypeError, match='S must be of type float32'): - r = mda.lib.distances.transform_StoR(s, box) diff --git a/testsuite/MDAnalysisTests/utils/test_distances.py b/testsuite/MDAnalysisTests/utils/test_distances.py index 5b29ca9b554..fa57a16d0ac 100644 --- a/testsuite/MDAnalysisTests/utils/test_distances.py +++ b/testsuite/MDAnalysisTests/utils/test_distances.py @@ -388,6 +388,10 @@ def positions(): d = np.array([[0., 0., 0.], [3., 3., 3.], [11., -11., 0.], [65., -65., 65.]], dtype=np.float32) return a, b, c, d + @staticmethod + def convert_position_dtype(a, b, c, d, dtype): + return a.astype(dtype), b.astype(dtype), c.astype(dtype), d.astype(dtype) + @staticmethod @pytest.fixture() def wrongtype(): @@ -402,8 +406,9 @@ def wronglength(): return np.array([[0., 0., 0.], [3., 3., 3.]], dtype=np.float32) - def test_bonds(self, positions, box, backend): - a, b, c, d = positions + @pytest.mark.parametrize('dtype', (np.float32, np.float64)) + def test_bonds(self, positions, box, backend, dtype): + a, b, c, d = self.convert_position_dtype(*positions, dtype=dtype) dists = MDAnalysis.lib.distances.calc_bonds(a, b, backend=backend) assert_equal(len(dists), 4, err_msg="calc_bonds results have wrong length") @@ -423,23 +428,6 @@ def test_bonds(self, positions, box, backend): assert_almost_equal(dists_pbc[3], 3.46410072, self.prec, err_msg="PBC check #w with box") - #Bad input checking - def test_bonds_wrongtype(self, positions, wrongtype, wronglength, backend): - a, b, c, d = positions - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_bonds(a, wrongtype, backend=backend) - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_bonds(wrongtype, b, backend=backend) - - with pytest.raises(ValueError): - MDAnalysis.lib.distances.calc_bonds(a, wronglength, backend=backend) - - with pytest.raises(ValueError): - MDAnalysis.lib.distances.calc_bonds(wronglength, b, backend= - backend) - - def test_bonds_badbox(self, positions, backend): a, b, c, d = positions badboxtype = np.array([10., 10., 10.], dtype=np.float64) @@ -469,8 +457,9 @@ def test_bonds_triclinic(self, positions, triclinic_box, backend): reference = np.array([0.0, 1.7320508, 1.4142136, 2.82842712]) assert_almost_equal(dists, reference, self.prec, err_msg="calc_bonds with triclinic box failed") - def test_angles(self, positions, backend): - a, b, c, d = positions + @pytest.mark.parametrize('dtype', (np.float32, np.float64)) + def test_angles(self, positions, backend, dtype): + a, b, c, d = self.convert_position_dtype(*positions, dtype=dtype) angles = MDAnalysis.lib.distances.calc_angles(a, b, c, backend=backend) # Check calculated values @@ -484,36 +473,6 @@ def test_angles(self, positions, backend): assert_almost_equal(angles[3], 0.098174833, self.prec, err_msg="Small angle failed in calc_angles") - # Check data type checks - def test_angles_wrongtype(self, positions, wrongtype, wronglength, backend): - a, b, c, d = positions - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_angles( - a, wrongtype, c, backend=backend) - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_angles() - wrongtype, b, c, backend = backend - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_angles( - a, b, wrongtype, backend=backend) - - with pytest.raises(ValueError): - MDAnalysis.lib.distances.calc_angles( - a, wronglength, c, - backend=backend) # try inputting arrays of different length - - with pytest.raises(ValueError): - MDAnalysis.lib.distances.calc_angles( - wronglength, b, c, - backend=backend) - - with pytest.raises(ValueError): - MDAnalysis.lib.distances.calc_angles( - a, b, wronglength, - backend=backend) - def test_angles_bad_result(self, positions, backend): a, b, c, d = positions badresult = np.zeros(len(a) - 1) @@ -521,8 +480,9 @@ def test_angles_bad_result(self, positions, backend): MDAnalysis.lib.distances.calc_angles( a, b, c, result=badresult, backend=backend) # Bad result array - def test_dihedrals(self, positions, backend): - a, b, c, d = positions + @pytest.mark.parametrize('dtype', (np.float32, np.float64)) + def test_dihedrals(self, positions, backend, dtype): + a, b, c, d = self.convert_position_dtype(*positions, dtype=dtype) dihedrals = MDAnalysis.lib.distances.calc_dihedrals(a, b, c, d, backend=backend) @@ -534,29 +494,6 @@ def test_dihedrals(self, positions, backend): assert_almost_equal(dihedrals[3], -0.50714064, self.prec, err_msg="arbitrary dihedral angle failed") - # Check data type checks - def test_dihedrals_wrongtype(self, positions, wrongtype, backend): - a, b, c, d = positions - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_dihedrals( - a, wrongtype, c, d, - backend=backend) # try inputting float64 values - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_dihedrals( - wrongtype, b, c, d, - backend=backend) - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_dihedrals( - a, b, wrongtype, d, - backend=backend) - - with pytest.raises(TypeError): - MDAnalysis.lib.distances.calc_dihedrals( - a, b, c, wrongtype, - backend=backend) - def test_dihedrals_wronglength(self, positions, wronglength, backend): a, b, c, d = positions with pytest.raises(ValueError):