13769110208

Committed 10 Mar 2025 03:57PM UTC coverage: 1.774% (+0.05%) from 1.726%

Build # 13769110208

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #53 from CCPBioSim/7-nan-error-handling

Add error handling for NaN and invalid values in calculations:

Run Details

2 of 33 new or added lines in 3 files covered. (6.06%)

3 existing lines in 2 files now uncovered.

55 of 3101 relevant lines covered (1.77%)

0.05 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.

    Input
    -----
    bead : the part of the molecule to be considered
    rot_axes : the axes relative to which the forces and coordinates are located
    frame : the frame number from the trajectory

    Output
    ------
    weighted_torque : 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):
        # Check if the moment of inertia is valid for square root calculation
        inertia_value = moment_of_inertia[dimension, dimension]

        if np.isclose(inertia_value, 0):
            raise ValueError(
                f"Invalid moment of inertia value: {inertia_value}. "
                f"Cannot compute weighted torque."
            )

        # compute weighted torque if moment of inertia is valid
        weighted_torque[dimension] = torques[dimension] / np.sqrt(inertia_value)

    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
NEW 185	if r2 == 0.0:	×
NEW 186	raise ValueError("r2 is zero, cannot compute spherical coordinates.")	×
187
NEW 188	if x2y2 == 0.0:	×
NEW 189	raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")	×
190
191	# Check for non-negative values inside the square root
NEW 192	if x2y2 / r2 < 0:	×
NEW 193	raise ValueError(	×
194	f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
195	f"Cannot take square root."
196	)
197
NEW 198	if x2y2 < 0:	×
NEW 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
NEW 283	if mass <= 0:	×
NEW 284	raise ValueError(	×
285	f"Invalid mass value: {mass}. Mass must be positive to compute the square "
286	f"root."
287	)
288
UNCOV 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	Input
302	-----
303	bead : the part of the molecule to be considered
304	rot_axes : the axes relative to which the forces and coordinates are located
305	frame : the frame number from the trajectory
306
307	Output
308	------
309	weighted_torque : the mass weighted sum of the torques in the bead
310	"""
311
312	torques = np.zeros((3,))	×
313	weighted_torque = np.zeros((3,))	×
314
315	for atom in bead.atoms:	×
316
317	# update local coordinates in rotational axes
318	coords_rot = data_container.atoms[atom.index].position - bead.center_of_mass()	×
319	coords_rot = np.matmul(rot_axes, coords_rot)	×
320	# update local forces in rotational frame
321	forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)	×
322
323	# multiply by the force_partitioning parameter to avoid double counting
324	# of the forces on weakly correlated atoms
325	# the default value of force_partitioning is 0.5 (dividing by two)
326	forces_rot = force_partitioning * forces_rot	×
327
328	# define torques (cross product of coordinates and forces) in rotational axes
329	torques_local = np.cross(coords_rot, forces_rot)	×
330	torques += torques_local	×
331
332	# divide by moment of inertia to get weighted torques
333	# moment of inertia is a 3x3 tensor
334	# the weighting is done in each dimension (x,y,z) using the diagonal elements of
335	# the moment of inertia tensor
336	moment_of_inertia = bead.moment_of_inertia()	×
337
338	for dimension in range(3):	×
339	# Check if the moment of inertia is valid for square root calculation
NEW 340	inertia_value = moment_of_inertia[dimension, dimension]	×
341
NEW 342	if np.isclose(inertia_value, 0):	×
NEW 343	raise ValueError(	×
344	f"Invalid moment of inertia value: {inertia_value}. "
345	f"Cannot compute weighted torque."
346	)
347
348	# compute weighted torque if moment of inertia is valid
NEW 349	weighted_torque[dimension] = torques[dimension] / np.sqrt(inertia_value)	×
350
UNCOV 351	return weighted_torque	×
352
353
354	# END
355
356
357	def create_submatrix(data_i, data_j, number_frames):	×
358	"""
359	Function for making covariance matrices.
360
361	Input
362	-----
363	data_i : values for bead i
364	data_j : valuees for bead j
365
366	Output
367	------
368	submatrix : 3x3 matrix for the covariance between i and j
369	"""
370
371	# Start with 3 by 3 matrix of zeros
372	submatrix = np.zeros((3, 3))	×
373
374	# For each frame calculate the outer product (cross product) of the data from the
375	# two beads and add the result to the submatrix
376	for frame in range(number_frames):	×
377	outer_product_matrix = np.outer(data_i[frame], data_j[frame])	×
378	submatrix = np.add(submatrix, outer_product_matrix)	×
379
380	# Divide by the number of frames to get the average
381	submatrix /= number_frames	×
382
383	return submatrix	×
384
385
386	# END
387
388
389	def filter_zero_rows_columns(arg_matrix, verbose):	×
390	"""
391	function for removing rows and columns that contain only zeros from a matrix
392
393	Input
394	-----
395	arg_matrix : matrix
396
397	Output
398	------
399	arg_matrix : the reduced size matrix
400	"""
401
402	# record the initial size
403	init_shape = np.shape(arg_matrix)	×
404
405	zero_indices = list(	×
406	filter(
407	lambda row: np.all(np.isclose(arg_matrix[row, :], 0.0)),
408	np.arange(np.shape(arg_matrix)[0]),
409	)
410	)
411	all_indices = np.ones((np.shape(arg_matrix)[0]), dtype=bool)	×
412	all_indices[zero_indices] = False	×
413	arg_matrix = arg_matrix[all_indices, :]	×
414
415	all_indices = np.ones((np.shape(arg_matrix)[1]), dtype=bool)	×
416	zero_indices = list(	×
417	filter(
418	lambda col: np.all(np.isclose(arg_matrix[:, col], 0.0)),
419	np.arange(np.shape(arg_matrix)[1]),
420	)
421	)
422	all_indices[zero_indices] = False	×
423	arg_matrix = arg_matrix[:, all_indices]	×
424
425	# get the final shape
426	final_shape = np.shape(arg_matrix)	×
427
428	if verbose and init_shape != final_shape:	×
429	print(	×
430	"A shape change has occured ({},{}) -> ({}, {})".format(
431	init_shape, final_shape
432	)
433	)
434
435	return arg_matrix	×
436
437
438	# END

CCPBioSim / CodeEntropy / 13769110208

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