13835820221

Committed 13 Mar 2025 01:29PM UTC coverage: 10.557% (-0.08%) from 10.638%

Build # 13835820221

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #64 from CCPBioSim/55-weighted-torque-bug

Weighted Torque Bug

Run Details

0 of 7 new or added lines in 1 file covered. (0.0%)

1 existing line in 1 file now uncovered.

55 of 521 relevant lines covered (10.56%)

0.32 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/CodeEntropy/GeometricFunctions.py

import numpy as np


def get_beads(data_container, level):
    """
    Function to define beads depending on the level in the hierarchy.

    Input
    -----
    data_container : the MDAnalysis universe
    level : the heirarchy level (polymer, residue, or united atom)

    Output
    ------
    list_of_beads : the relevent beads
    """

    if level == "polymer":
        list_of_beads = []
        atom_group = "all"
        list_of_beads.append(data_container.select_atoms(atom_group))

    if level == "residue":
        list_of_beads = []
        num_residues = len(data_container.residues)
        for residue in range(num_residues):
            atom_group = "resindex " + str(residue)
            list_of_beads.append(data_container.select_atoms(atom_group))

    if level == "united_atom":
        list_of_beads = []
        heavy_atoms = data_container.select_atoms("not name H*")
        for atom in heavy_atoms:
            atom_group = (
                "index "
                + str(atom.index)
                + " or (name H* and bonded index "
                + str(atom.index)
                + ")"
            )
            list_of_beads.append(data_container.select_atoms(atom_group))

    return list_of_beads


# END


def get_axes(data_container, level, index=0):
    """
    Function to set the translational and rotational axes.
    The translational axes are based on the principal axes of the unit one level larger
    than the level we are interested in (except for the polymer level where there is no
    larger unit). The rotational axes use the covalent links between residues or atoms
    where possible to define the axes, or if the unit is not bonded to others of the
    same level the prinicpal axes of the unit are used.

    Input
    -----
    data_container : the information about the molecule and trajectory
    level : the level (united atom, residue, or polymer) of interest
    index : residue index (integer)

    Output
    ------
    trans_axes : translational axes
    rot_axes : rotational axes
    """
    index = int(index)

    if level == "polymer":
        # for polymer use principle axis for both translation and rotation
        trans_axes = data_container.atoms.principal_axes()
        rot_axes = data_container.atoms.principal_axes()

    if level == "residue":
        # Translation
        # for residues use principal axes of whole molecule for translation
        trans_axes = data_container.atoms.principal_axes()

        # Rotation
        # find bonds between atoms in residue of interest and other residues
        # we are assuming bonds only exist between adjacent residues
        # (linear chains of residues)
        # TODO refine selection so that it will work for branched polymers
        index_prev = index - 1
        index_next = index + 1
        atom_set = data_container.select_atoms(
            f"(resindex {index_prev} or resindex {index_next}) and bonded resid {index}"
        )
        residue = data_container.select_atoms(f"resindex {index}")

        if len(atom_set) == 0:
            # if no bonds to other residues use pricipal axes of residue
            rot_axes = residue.atoms.principal_axes()

        else:
            # set center of rotation to center of mass of the residue
            center = residue.atoms.center_of_mass()

            # get vector for average position of bonded atoms
            vector = get_avg_pos(atom_set, center)

            # use spherical coordinates function to get rotational axes
            rot_axes = get_sphCoord_axes(vector)

    if level == "united_atom":
        # Translation
        # for united atoms use principal axes of residue for translation
        trans_axes = data_container.residues.principal_axes()

        # Rotation
        # for united atoms use heavy atoms bonded to the heavy atom
        atom_set = data_container.select_atoms(f"not name H* and bonded index {index}")

        # center at position of heavy atom
        atom_group = data_container.select_atoms(f"index {index}")
        center = atom_group.positions[0]

        # get vector for average position of hydrogens
        vector = get_avg_pos(atom_set, center)

        # use spherical coordinates function to get rotational axes
        rot_axes = get_sphCoord_axes(vector)

    return trans_axes, rot_axes


# END


def get_avg_pos(atom_set, center):
    """
    Function to get the average position of a set of atoms.

    Input
    -----
    atoms : MDAnalysis atom group
    center : position for center of rotation

    Output
    ------
    avg_position : three dimensional vector
    """
    # start with an empty vector
    avg_position = np.zeros((3))

    # get number of atoms
    number_atoms = len(atom_set.names)

    if number_atoms != 0:
        # sum positions for all atoms in the given set
        for atom_index in range(number_atoms):
            atom_position = atom_set.atoms[atom_index].position

            avg_position += atom_position

        avg_position /= number_atoms  # divide by number of atoms to get average

    else:
        # if no atoms in set the unit has no bonds to restrict its rotational motion,
        # so we can use a random vector to get the spherical coordinates axes
        avg_position = np.random.random(3)

    # transform the average position to a coordinate system with the origin at center
    avg_position = avg_position - center

    return avg_position


# END


def get_sphCoord_axes(arg_r):
    """
    For a given vector in space, treat it is a radial vector rooted at 0,0,0 and
    derive a curvilinear coordinate system according to the rules of polar spherical
    coordinates
    """

    x2y2 = arg_r[0] ** 2 + arg_r[1] ** 2
    r2 = x2y2 + arg_r[2] ** 2

    # Check for division by zero
    if r2 == 0.0:
        raise ValueError("r2 is zero, cannot compute spherical coordinates.")

    if x2y2 == 0.0:
        raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")

    # Check for non-negative values inside the square root
    if x2y2 / r2 < 0:
        raise ValueError(
            f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
            f"Cannot take square root."
        )

    if x2y2 < 0:
        raise ValueError(
            f"Negative value encountered for sin_phi and cos_phi calculation: {x2y2}. "
            f"Cannot take square root."
        )

    if x2y2 != 0.0:
        sin_theta = np.sqrt(x2y2 / r2)
        cos_theta = arg_r[2] / np.sqrt(r2)

        sin_phi = arg_r[1] / np.sqrt(x2y2)
        cos_phi = arg_r[0] / np.sqrt(x2y2)

    else:
        sin_theta = 0.0
        cos_theta = 1

        sin_phi = 0.0
        cos_phi = 1

    # if abs(sin_theta) > 1 or abs(sin_phi) > 1:
    #     print('Bad sine : T {} , P {}'.format(sin_theta, sin_phi))

    # cos_theta = np.sqrt(1 - sin_theta*sin_theta)
    # cos_phi = np.sqrt(1 - sin_phi*sin_phi)

    # print('{} {} {}'.format(*arg_r))
    # print('Sin T : {}, cos T : {}'.format(sin_theta, cos_theta))
    # print('Sin P : {}, cos P : {}'.format(sin_phi, cos_phi))

    spherical_basis = np.zeros((3, 3))

    # r^
    spherical_basis[0, :] = np.asarray(
        [sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]
    )

    # Theta^
    spherical_basis[1, :] = np.asarray(
        [cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta]
    )

    # Phi^
    spherical_basis[2, :] = np.asarray([-sin_phi, cos_phi, 0.0])

    return spherical_basis


# END


def get_weighted_forces(
    data_container, bead, trans_axes, highest_level, force_partitioning=0.5
):
    """
    Function to calculate the mass weighted forces for a given bead.

    Input
    -----
    bead : the part of the system to be considered
    trans_axes : the axes relative to which the forces are located

    Output
    ------
    weighted_force : the mass weighted sum of the forces in the bead
    """

    forces_trans = np.zeros((3,))

    # Sum forces from all atoms in the bead
    for atom in bead.atoms:
        # update local forces in translational axes
        forces_local = np.matmul(trans_axes, data_container.atoms[atom.index].force)
        forces_trans += forces_local

    if highest_level:
        # multiply by the force_partitioning parameter to avoid double counting
        # of the forces on weakly correlated atoms
        # the default value of force_partitioning is 0.5 (dividing by two)
        forces_trans = force_partitioning * forces_trans

    # divide the sum of forces by the mass of the bead to get the weighted forces
    mass = bead.total_mass()

    # Check that mass is positive to avoid division by 0 or negative values inside sqrt
    if mass <= 0:
        raise ValueError(
            f"Invalid mass value: {mass}. Mass must be positive to compute the square "
            f"root."
        )

    weighted_force = forces_trans / np.sqrt(mass)

    return weighted_force


# END


def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5):
    """
    Function to calculate the moment of inertia weighted torques for a given bead.

    This function computes torques in a rotated frame and then weights them using
    the moment of inertia tensor. To prevent numerical instability, it treats
    extremely small diagonal elements of the moment of inertia tensor as zero
    (since values below machine precision are effectively zero). This avoids
    unnecessary use of extended precision (e.g., float128).

    Additionally, if the computed torque is already zero, the function skips
    the division step, reducing unnecessary computations and potential errors.

    Parameters
    ----------
    data_container : object
        Contains atomic positions and forces.
    bead : object
        The part of the molecule to be considered.
    rot_axes : np.ndarray
        The axes relative to which the forces and coordinates are located.
    force_partitioning : float, optional
        Factor to adjust force contributions, default is 0.5.

    Returns
    -------
    np.ndarray
        The mass-weighted sum of the torques in the bead.
    """

    torques = np.zeros((3,))
    weighted_torque = np.zeros((3,))

    for atom in bead.atoms:

        # update local coordinates in rotational axes
        coords_rot = data_container.atoms[atom.index].position - bead.center_of_mass()
        coords_rot = np.matmul(rot_axes, coords_rot)
        # update local forces in rotational frame
        forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)

        # multiply by the force_partitioning parameter to avoid double counting
        # of the forces on weakly correlated atoms
        # the default value of force_partitioning is 0.5 (dividing by two)
        forces_rot = force_partitioning * forces_rot

        # define torques (cross product of coordinates and forces) in rotational axes
        torques_local = np.cross(coords_rot, forces_rot)
        torques += torques_local

    # divide by moment of inertia to get weighted torques
    # moment of inertia is a 3x3 tensor
    # the weighting is done in each dimension (x,y,z) using the diagonal elements of
    # the moment of inertia tensor
    moment_of_inertia = bead.moment_of_inertia()

    for dimension in range(3):
        # Skip calculation if torque is already zero
        if np.isclose(torques[dimension], 0):
            weighted_torque[dimension] = 0
            continue

        # Check for zero moment of inertia
        if np.isclose(moment_of_inertia[dimension, dimension], 0):
            raise ZeroDivisionError(
                f"Attempted to divide by zero moment of inertia in dimension "
                f"{dimension}."
            )

        # Check for negative moment of inertia
        if moment_of_inertia[dimension, dimension] < 0:
            raise ValueError(
                f"Negative value encountered for moment of inertia: "
                f"{moment_of_inertia[dimension, dimension]} "
                f"Cannot compute weighted torque."
            )

        # Compute weighted torque
        weighted_torque[dimension] = torques[dimension] / np.sqrt(
            moment_of_inertia[dimension, dimension]
        )

    return weighted_torque


# END


def create_submatrix(data_i, data_j, number_frames):
    """
    Function for making covariance matrices.

    Input
    -----
    data_i : values for bead i
    data_j : valuees for bead j

    Output
    ------
    submatrix : 3x3 matrix for the covariance between i and j
    """

    # Start with 3 by 3 matrix of zeros
    submatrix = np.zeros((3, 3))

    # For each frame calculate the outer product (cross product) of the data from the
    # two beads and add the result to the submatrix
    for frame in range(number_frames):
        outer_product_matrix = np.outer(data_i[frame], data_j[frame])
        submatrix = np.add(submatrix, outer_product_matrix)

    # Divide by the number of frames to get the average
    submatrix /= number_frames

    return submatrix


# END


def filter_zero_rows_columns(arg_matrix, verbose):
    """
    function for removing rows and columns that contain only zeros from a matrix

    Input
    -----
    arg_matrix : matrix

    Output
    ------
    arg_matrix : the reduced size matrix
    """

    # record the initial size
    init_shape = np.shape(arg_matrix)

    zero_indices = list(
        filter(
            lambda row: np.all(np.isclose(arg_matrix[row, :], 0.0)),
            np.arange(np.shape(arg_matrix)[0]),
        )
    )
    all_indices = np.ones((np.shape(arg_matrix)[0]), dtype=bool)
    all_indices[zero_indices] = False
    arg_matrix = arg_matrix[all_indices, :]

    all_indices = np.ones((np.shape(arg_matrix)[1]), dtype=bool)
    zero_indices = list(
        filter(
            lambda col: np.all(np.isclose(arg_matrix[:, col], 0.0)),
            np.arange(np.shape(arg_matrix)[1]),
        )
    )
    all_indices[zero_indices] = False
    arg_matrix = arg_matrix[:, all_indices]

    # get the final shape
    final_shape = np.shape(arg_matrix)

    if verbose and init_shape != final_shape:
        print(
            "A shape change has occured ({},{}) -> ({}, {})".format(
                *init_shape, *final_shape
            )
        )

    return arg_matrix


# END

1	import numpy as np	×
2
3
4	def get_beads(data_container, level):	×
5	"""
6	Function to define beads depending on the level in the hierarchy.
7
8	Input
9	-----
10	data_container : the MDAnalysis universe
11	level : the heirarchy level (polymer, residue, or united atom)
12
13	Output
14	------
15	list_of_beads : the relevent beads
16	"""
17
18	if level == "polymer":	×
19	list_of_beads = []	×
20	atom_group = "all"	×
21	list_of_beads.append(data_container.select_atoms(atom_group))	×
22
23	if level == "residue":	×
24	list_of_beads = []	×
25	num_residues = len(data_container.residues)	×
26	for residue in range(num_residues):	×
27	atom_group = "resindex " + str(residue)	×
28	list_of_beads.append(data_container.select_atoms(atom_group))	×
29
30	if level == "united_atom":	×
31	list_of_beads = []	×
32	heavy_atoms = data_container.select_atoms("not name H*")	×
33	for atom in heavy_atoms:	×
34	atom_group = (	×
35	"index "
36	+ str(atom.index)
37	+ " or (name H* and bonded index "
38	+ str(atom.index)
39	+ ")"
40	)
41	list_of_beads.append(data_container.select_atoms(atom_group))	×
42
43	return list_of_beads	×
44
45
46	# END
47
48
49	def get_axes(data_container, level, index=0):	×
50	"""
51	Function to set the translational and rotational axes.
52	The translational axes are based on the principal axes of the unit one level larger
53	than the level we are interested in (except for the polymer level where there is no
54	larger unit). The rotational axes use the covalent links between residues or atoms
55	where possible to define the axes, or if the unit is not bonded to others of the
56	same level the prinicpal axes of the unit are used.
57
58	Input
59	-----
60	data_container : the information about the molecule and trajectory
61	level : the level (united atom, residue, or polymer) of interest
62	index : residue index (integer)
63
64	Output
65	------
66	trans_axes : translational axes
67	rot_axes : rotational axes
68	"""
69	index = int(index)	×
70
71	if level == "polymer":	×
72	# for polymer use principle axis for both translation and rotation
73	trans_axes = data_container.atoms.principal_axes()	×
74	rot_axes = data_container.atoms.principal_axes()	×
75
76	if level == "residue":	×
77	# Translation
78	# for residues use principal axes of whole molecule for translation
79	trans_axes = data_container.atoms.principal_axes()	×
80
81	# Rotation
82	# find bonds between atoms in residue of interest and other residues
83	# we are assuming bonds only exist between adjacent residues
84	# (linear chains of residues)
85	# TODO refine selection so that it will work for branched polymers
86	index_prev = index - 1	×
87	index_next = index + 1	×
88	atom_set = data_container.select_atoms(	×
89	f"(resindex {index_prev} or resindex {index_next}) and bonded resid {index}"
90	)
91	residue = data_container.select_atoms(f"resindex {index}")	×
92
93	if len(atom_set) == 0:	×
94	# if no bonds to other residues use pricipal axes of residue
95	rot_axes = residue.atoms.principal_axes()	×
96
97	else:
98	# set center of rotation to center of mass of the residue
99	center = residue.atoms.center_of_mass()	×
100
101	# get vector for average position of bonded atoms
102	vector = get_avg_pos(atom_set, center)	×
103
104	# use spherical coordinates function to get rotational axes
105	rot_axes = get_sphCoord_axes(vector)	×
106
107	if level == "united_atom":	×
108	# Translation
109	# for united atoms use principal axes of residue for translation
110	trans_axes = data_container.residues.principal_axes()	×
111
112	# Rotation
113	# for united atoms use heavy atoms bonded to the heavy atom
114	atom_set = data_container.select_atoms(f"not name H* and bonded index {index}")	×
115
116	# center at position of heavy atom
117	atom_group = data_container.select_atoms(f"index {index}")	×
118	center = atom_group.positions[0]	×
119
120	# get vector for average position of hydrogens
121	vector = get_avg_pos(atom_set, center)	×
122
123	# use spherical coordinates function to get rotational axes
124	rot_axes = get_sphCoord_axes(vector)	×
125
126	return trans_axes, rot_axes	×
127
128
129	# END
130
131
132	def get_avg_pos(atom_set, center):	×
133	"""
134	Function to get the average position of a set of atoms.
135
136	Input
137	-----
138	atoms : MDAnalysis atom group
139	center : position for center of rotation
140
141	Output
142	------
143	avg_position : three dimensional vector
144	"""
145	# start with an empty vector
146	avg_position = np.zeros((3))	×
147
148	# get number of atoms
149	number_atoms = len(atom_set.names)	×
150
151	if number_atoms != 0:	×
152	# sum positions for all atoms in the given set
153	for atom_index in range(number_atoms):	×
154	atom_position = atom_set.atoms[atom_index].position	×
155
156	avg_position += atom_position	×
157
158	avg_position /= number_atoms # divide by number of atoms to get average	×
159
160	else:
161	# if no atoms in set the unit has no bonds to restrict its rotational motion,
162	# so we can use a random vector to get the spherical coordinates axes
163	avg_position = np.random.random(3)	×
164
165	# transform the average position to a coordinate system with the origin at center
166	avg_position = avg_position - center	×
167
168	return avg_position	×
169
170
171	# END
172
173
174	def get_sphCoord_axes(arg_r):	×
175	"""
176	For a given vector in space, treat it is a radial vector rooted at 0,0,0 and
177	derive a curvilinear coordinate system according to the rules of polar spherical
178	coordinates
179	"""
180
181	x2y2 = arg_r[0] 2 + arg_r[1] 2	×
182	r2 = x2y2 + arg_r[2] ** 2	×
183
184	# Check for division by zero
185	if r2 == 0.0:	×
186	raise ValueError("r2 is zero, cannot compute spherical coordinates.")	×
187
188	if x2y2 == 0.0:	×
189	raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")	×
190
191	# Check for non-negative values inside the square root
192	if x2y2 / r2 < 0:	×
193	raise ValueError(	×
194	f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
195	f"Cannot take square root."
196	)
197
198	if x2y2 < 0:	×
199	raise ValueError(	×
200	f"Negative value encountered for sin_phi and cos_phi calculation: {x2y2}. "
201	f"Cannot take square root."
202	)
203
204	if x2y2 != 0.0:	×
205	sin_theta = np.sqrt(x2y2 / r2)	×
206	cos_theta = arg_r[2] / np.sqrt(r2)	×
207
208	sin_phi = arg_r[1] / np.sqrt(x2y2)	×
209	cos_phi = arg_r[0] / np.sqrt(x2y2)	×
210
211	else:
212	sin_theta = 0.0	×
213	cos_theta = 1	×
214
215	sin_phi = 0.0	×
216	cos_phi = 1	×
217
218	# if abs(sin_theta) > 1 or abs(sin_phi) > 1:
219	# print('Bad sine : T {} , P {}'.format(sin_theta, sin_phi))
220
221	# cos_theta = np.sqrt(1 - sin_theta*sin_theta)
222	# cos_phi = np.sqrt(1 - sin_phi*sin_phi)
223
224	# print('{} {} {}'.format(*arg_r))
225	# print('Sin T : {}, cos T : {}'.format(sin_theta, cos_theta))
226	# print('Sin P : {}, cos P : {}'.format(sin_phi, cos_phi))
227
228	spherical_basis = np.zeros((3, 3))	×
229
230	# r^
231	spherical_basis[0, :] = np.asarray(	×
232	[sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]
233	)
234
235	# Theta^
236	spherical_basis[1, :] = np.asarray(	×
237	[cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta]
238	)
239
240	# Phi^
241	spherical_basis[2, :] = np.asarray([-sin_phi, cos_phi, 0.0])	×
242
243	return spherical_basis	×
244
245
246	# END
247
248
249	def get_weighted_forces(	×
250	data_container, bead, trans_axes, highest_level, force_partitioning=0.5
251	):
252	"""
253	Function to calculate the mass weighted forces for a given bead.
254
255	Input
256	-----
257	bead : the part of the system to be considered
258	trans_axes : the axes relative to which the forces are located
259
260	Output
261	------
262	weighted_force : the mass weighted sum of the forces in the bead
263	"""
264
265	forces_trans = np.zeros((3,))	×
266
267	# Sum forces from all atoms in the bead
268	for atom in bead.atoms:	×
269	# update local forces in translational axes
270	forces_local = np.matmul(trans_axes, data_container.atoms[atom.index].force)	×
271	forces_trans += forces_local	×
272
273	if highest_level:	×
274	# multiply by the force_partitioning parameter to avoid double counting
275	# of the forces on weakly correlated atoms
276	# the default value of force_partitioning is 0.5 (dividing by two)
277	forces_trans = force_partitioning * forces_trans	×
278
279	# divide the sum of forces by the mass of the bead to get the weighted forces
280	mass = bead.total_mass()	×
281
282	# Check that mass is positive to avoid division by 0 or negative values inside sqrt
283	if mass <= 0:	×
284	raise ValueError(	×
285	f"Invalid mass value: {mass}. Mass must be positive to compute the square "
286	f"root."
287	)
288
289	weighted_force = forces_trans / np.sqrt(mass)	×
290
291	return weighted_force	×
292
293
294	# END
295
296
297	def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5):	×
298	"""
299	Function to calculate the moment of inertia weighted torques for a given bead.
300
301	This function computes torques in a rotated frame and then weights them using
302	the moment of inertia tensor. To prevent numerical instability, it treats
303	extremely small diagonal elements of the moment of inertia tensor as zero
304	(since values below machine precision are effectively zero). This avoids
305	unnecessary use of extended precision (e.g., float128).
306
307	Additionally, if the computed torque is already zero, the function skips
308	the division step, reducing unnecessary computations and potential errors.
309
310	Parameters
311	----------
312	data_container : object
313	Contains atomic positions and forces.
314	bead : object
315	The part of the molecule to be considered.
316	rot_axes : np.ndarray
317	The axes relative to which the forces and coordinates are located.
318	force_partitioning : float, optional
319	Factor to adjust force contributions, default is 0.5.
320
321	Returns
322	-------
323	np.ndarray
324	The mass-weighted sum of the torques in the bead.
325	"""
326
327	torques = np.zeros((3,))	×
328	weighted_torque = np.zeros((3,))	×
329
330	for atom in bead.atoms:	×
331
332	# update local coordinates in rotational axes
333	coords_rot = data_container.atoms[atom.index].position - bead.center_of_mass()	×
334	coords_rot = np.matmul(rot_axes, coords_rot)	×
335	# update local forces in rotational frame
336	forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)	×
337
338	# multiply by the force_partitioning parameter to avoid double counting
339	# of the forces on weakly correlated atoms
340	# the default value of force_partitioning is 0.5 (dividing by two)
341	forces_rot = force_partitioning * forces_rot	×
342
343	# define torques (cross product of coordinates and forces) in rotational axes
344	torques_local = np.cross(coords_rot, forces_rot)	×
345	torques += torques_local	×
346
347	# divide by moment of inertia to get weighted torques
348	# moment of inertia is a 3x3 tensor
349	# the weighting is done in each dimension (x,y,z) using the diagonal elements of
350	# the moment of inertia tensor
351	moment_of_inertia = bead.moment_of_inertia()	×
352
353	for dimension in range(3):	×
354	# Skip calculation if torque is already zero
NEW 355	if np.isclose(torques[dimension], 0):	×
NEW 356	weighted_torque[dimension] = 0	×
NEW 357	continue	×
358
359	# Check for zero moment of inertia
NEW 360	if np.isclose(moment_of_inertia[dimension, dimension], 0):	×
NEW 361	raise ZeroDivisionError(	×
362	f"Attempted to divide by zero moment of inertia in dimension "
363	f"{dimension}."
364	)
365
366	# Check for negative moment of inertia
NEW 367	if moment_of_inertia[dimension, dimension] < 0:	×
UNCOV 368	raise ValueError(	×
369	f"Negative value encountered for moment of inertia: "
370	f"{moment_of_inertia[dimension, dimension]} "
371	f"Cannot compute weighted torque."
372	)
373
374	# Compute weighted torque
NEW 375	weighted_torque[dimension] = torques[dimension] / np.sqrt(	×
376	moment_of_inertia[dimension, dimension]
377	)
378
379	return weighted_torque	×
380
381
382	# END
383
384
385	def create_submatrix(data_i, data_j, number_frames):	×
386	"""
387	Function for making covariance matrices.
388
389	Input
390	-----
391	data_i : values for bead i
392	data_j : valuees for bead j
393
394	Output
395	------
396	submatrix : 3x3 matrix for the covariance between i and j
397	"""
398
399	# Start with 3 by 3 matrix of zeros
400	submatrix = np.zeros((3, 3))	×
401
402	# For each frame calculate the outer product (cross product) of the data from the
403	# two beads and add the result to the submatrix
404	for frame in range(number_frames):	×
405	outer_product_matrix = np.outer(data_i[frame], data_j[frame])	×
406	submatrix = np.add(submatrix, outer_product_matrix)	×
407
408	# Divide by the number of frames to get the average
409	submatrix /= number_frames	×
410
411	return submatrix	×
412
413
414	# END
415
416
417	def filter_zero_rows_columns(arg_matrix, verbose):	×
418	"""
419	function for removing rows and columns that contain only zeros from a matrix
420
421	Input
422	-----
423	arg_matrix : matrix
424
425	Output
426	------
427	arg_matrix : the reduced size matrix
428	"""
429
430	# record the initial size
431	init_shape = np.shape(arg_matrix)	×
432
433	zero_indices = list(	×
434	filter(
435	lambda row: np.all(np.isclose(arg_matrix[row, :], 0.0)),
436	np.arange(np.shape(arg_matrix)[0]),
437	)
438	)
439	all_indices = np.ones((np.shape(arg_matrix)[0]), dtype=bool)	×
440	all_indices[zero_indices] = False	×
441	arg_matrix = arg_matrix[all_indices, :]	×
442
443	all_indices = np.ones((np.shape(arg_matrix)[1]), dtype=bool)	×
444	zero_indices = list(	×
445	filter(
446	lambda col: np.all(np.isclose(arg_matrix[:, col], 0.0)),
447	np.arange(np.shape(arg_matrix)[1]),
448	)
449	)
450	all_indices[zero_indices] = False	×
451	arg_matrix = arg_matrix[:, all_indices]	×
452
453	# get the final shape
454	final_shape = np.shape(arg_matrix)	×
455
456	if verbose and init_shape != final_shape:	×
457	print(	×
458	"A shape change has occured ({},{}) -> ({}, {})".format(
459	init_shape, final_shape
460	)
461	)
462
463	return arg_matrix	×
464
465
466	# END

CCPBioSim / CodeEntropy / 13835820221

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous