Skip to content

Commit

Permalink
Added Google-style docstrings to all functions+classes in algorithms
Browse files Browse the repository at this point in the history
  • Loading branch information
PrinceWalnut committed Dec 13, 2023
1 parent b11ac43 commit d2133e6
Show file tree
Hide file tree
Showing 5 changed files with 327 additions and 209 deletions.
141 changes: 54 additions & 87 deletions src/laue_dials/algorithms/diffgeo.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,19 +9,14 @@
def get_UB_matrices(A):
"""
Convert the reciprocal space indexing matrix, A into the product of
an orthogonal matrix U and a lower triangular matrix B
an orthogonal matrix U and a lower triangular matrix B.
Parameters
----------
A : np.array
A 3x3 indexing matrix such that the scattering vector `S1-S0=Q=A@h`
Parameters:
A (np.array): A 3x3 indexing matrix such that the scattering vector `S1-S0=Q=A@h`.
Returns
-------
U : np.array
An orthogonal 3x3 matrix
B : np.array
A 3x3 lower triangular matrix
Returns:
U (np.array): An orthogonal 3x3 matrix with the crystal U matrix.
B (np.array): A 3x3 lower triangular matrix for the crystal B matrix.
"""
from scipy.linalg import rq

Expand All @@ -34,17 +29,13 @@ def get_UB_matrices(A):

def normalize(A):
"""
Normalize the last dimension of an array by dividing by its L2 norm
Normalize the last dimension of an array by dividing by its L2 norm.
Parameters
----------
A : np.array
A non-normal matrix
Parameters:
A (np.array): A non-normal matrix.
Returns
-------
A : np.array
A normalized matrix
Returns:
A (np.array): A normalized matrix.
"""
normalized_A = A / np.linalg.norm(A, axis=-1)[..., None]

Expand All @@ -58,19 +49,13 @@ def hkl2ray(hkl, wavelength=None):
Convert a miller index to the shortest member of its central ray.
Optionally, adjust its wavelength accordingly.
Parameters
----------
hkl : array
`n x 3` array of miller indices. the dtype must be interpretable as an integer
wavelength : array (optional)
length `n` array of wavelengths corresponding to each miller index
Parameters:
hkl (np.array): `n x 3` array of miller indices. The dtype must be interpretable as an integer.
wavelength (Optional[np.array]): Length `n` array of wavelengths corresponding to each miller index.
Returns
-------
reduced_hkl : array
The miller index of the shortest vector on the same central ray as hkl
reduced_wavelength : array (optional)
The wavelengths corresponding to reduced_hkl
Returns:
reduced_hkl (np.array): The miller index of the shortest vector on the same central ray as the original hkl.
reduced_wavelength (Optional[np.array]): The wavelengths corresponding to reduced_hkl.
"""
gcd = np.gcd.reduce(hkl.astype(int), axis=-1)
if wavelength is not None:
Expand All @@ -81,19 +66,15 @@ def hkl2ray(hkl, wavelength=None):

def is_ray_equivalent(hkl1, hkl2):
"""
Test for equivalency between two miller indices in a Laue experiment. Returns a boolean array for each of the `n` hkls in `hkl{1,2}`.
Test for equivalency between two miller indices in a Laue experiment.
Returns a boolean array for each of the `n` hkls in `hkl{1,2}`.
Parameters
----------
hkl1 : array
`n x 3` array of miller indices. the dtype must be interpretable as an integer
hkl2 : array
`n x 3` array of miller indices. the dtype must be interpretable as an integer
Parameters:
hkl1 (np.array): `n x 3` array of miller indices. The dtype must be interpretable as an integer.
hkl2 (np.array): `n x 3` array of miller indices. The dtype must be interpretable as an integer.
Returns
-------
equal_hkl : array
Boolean array for hkl equivalency by index from original arrays
Returns:
equal_hkl (np.array): Boolean array for hkl equivalency by index from original arrays.
"""
return np.all(np.isclose(hkl2ray(hkl1), hkl2ray(hkl2)), axis=1)

Expand All @@ -102,21 +83,14 @@ def align_hkls(reference, target, spacegroup, anomalous=True):
"""
Use the point group operators in `spacegroup` to align target hkls to reference.
Parameters
----------
reference : array
n x 3 array of miller indices
target : array
n x 3 array of miller indices
spacegroup : gemmi.Spacegroup
The space group of reference/target.
anomalous : bool(optional)
If true, test Friedel symmetries too.
Returns
-------
aligned : array
n x 3 array of miller indices equivalent to target
Parameters:
reference (np.array): n x 3 array of miller indices.
target (np.array): n x 3 array of miller indices.
spacegroup (gemmi.Spacegroup): The space group of reference/target.
anomalous (bool, optional): If true, test Friedel symmetries too.
Returns:
aligned (np.array): n x 3 array of miller indices equivalent to target.
"""
aligned = target
cc = -1.0
Expand All @@ -140,19 +114,14 @@ def align_hkls(reference, target, spacegroup, anomalous=True):

def orthogonalization(a, b, c, alpha, beta, gamma):
"""
Compute the orthogonalization matrix from cell params
Compute the orthogonalization matrix from cell params.
Parameters
----------
a,b,c : floats
Floats representing the magnitudes of the three unit cell axes
alpha, beta, gamma : floats
Floats representing the three unit cell angles
Parameters:
a, b, c (float): Floats representing the magnitudes of the three unit cell axes.
alpha, beta, gamma (float): Floats representing the three unit cell angles.
Returns
-------
orthogonalization_matrix : array
3 x 3 orthogonalization matrix for unit cell
Returns:
orthogonalization_matrix (np.array): 3 x 3 orthogonalization matrix for unit cell.
"""
alpha, beta, gamma = (
np.pi * alpha / 180.0,
Expand Down Expand Up @@ -182,17 +151,16 @@ def orthogonalization(a, b, c, alpha, beta, gamma):

def mat_to_rot_xyz(R, deg=True):
"""
Decompose a rotation matrix into euler angles
Decompose a rotation matrix into euler angles.
Parameters
----------
R : array
3 x 3 rotation matrix
Parameters:
R (np.array): 3 x 3 rotation matrix.
deg (Optional[bool]): If True, angles are returned in degrees.
Returns
-------
rot_x, rot_y, rot_z : floats
Euler angles to generate rotation matrix
Returns:
rot_x (float): Euler angle around x-axis to generate rotation matrix.
rot_y (float): Euler angle around y-axis to generate rotation matrix.
rot_z (float): Euler angle around z-axis to generate rotation matrix.
"""
if R[2, 0] < 1:
if R[2, 0] > -1:
Expand All @@ -216,17 +184,16 @@ def mat_to_rot_xyz(R, deg=True):

def rot_xyz_to_mat(rot_x, rot_y, rot_z, deg=True):
"""
Convert euler angles into a rotation matrix
Convert euler angles into a rotation matrix.
Parameters
----------
rot_x, rot_y, rot_z : floats
Euler angles to generate rotation matrix
Parameters:
rot_x (float): Euler angle around x-axis to generate rotation matrix.
rot_y (float): Euler angle around y-axis to generate rotation matrix.
rot_z (float): Euler angle around z-axis to generate rotation matrix.
deg (Optional[bool]): If True, input angles are assumed to be in degrees.
Returns
-------
R : array
3 x 3 rotation matrix
Returns:
R (np.array): 3 x 3 rotation matrix.
"""
if deg:
rot_x = np.deg2rad(rot_x)
Expand Down
97 changes: 97 additions & 0 deletions src/laue_dials/algorithms/integration.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,21 @@

class Profile:
def __init__(
"""
Initialize a Profile instance.
Parameters:
x (np.array): x-coordinates.
y (np.array): y-coordinates.
counts (np.array): Counts associated with each coordinate.
cen_x (Optional[np.array]): Center x-coordinate.
cen_y (Optional[np.array]): Center y-coordinate.
fg_cutoff (float): Foreground cutoff value.
bg_cutoff (float): Background cutoff value.
minfrac (float): Minimum fraction of pixels in a profile.
eps (float): Small positive epsilon value to remove negative intensities.
frac_step_size (float): Fractional step size for fitting iterations.
"""
self,
x,
y,
Expand Down Expand Up @@ -49,6 +64,19 @@ def __init__(
self.update_mask()

def set_profile_params(self, scale, slope, intercept, cen_x, cen_y):
"""
Set profile parameters.
Parameters:
scale (np.array): Scaling factor.
slope (np.array): Slope of the profile.
intercept (float): Intercept value.
cen_x (np.array): Center x-coordinate.
cen_y (np.array): Center y-coordinate.
Returns:
self (Profile): Updated Profile instance.
"""
self.scale = scale
self.slope = slope
self.intercept = intercept
Expand All @@ -61,6 +89,12 @@ def set_profile_params(self, scale, slope, intercept, cen_x, cen_y):
return self

def fit(self, nsteps=5):
"""
Fit the profile.
Parameters:
nsteps (int): Number of fitting steps.
"""
self.success = True
for i in range(nsteps):
try:
Expand All @@ -73,20 +107,52 @@ def fit(self, nsteps=5):

@classmethod
def from_dataframe(cls, df, x_key="dx", y_key="dy", count_key="counts", **kwargs):
"""
Create a Profile instance from a DataFrame.
Parameters:
df (pd.DataFrame): Input DataFrame.
x_key (str): x-coordinate column name.
y_key (str): y-coordinate column name.
count_key (str): Counts column name.
kwargs: Additional keyword arguments.
Returns:
Profile: Created Profile instance.
"""
x = df[x_key].astype("float32")
y = df[y_key].astype("float32")
counts = df[count_key].astype("float32")
return cls(x, y, counts, **kwargs)

@property
def difference_vectors(self):
"""
Get the difference vectors.
Returns:
np.array: Difference vectors.
"""
return np.column_stack((self.x - self.cen_x, self.y - self.cen_y))

@property
def background(self):
"""
Get the background.
Returns:
np.array: Background values.
"""
return self.difference_vectors @ self.slope + self.intercept

def update_mask(self, fg_cutoff=None, bg_cutoff=None):
"""
Update the mask.
Parameters:
fg_cutoff (Optional[float]): Foreground cutoff value.
bg_cutoff (Optional[float]): Background cutoff value.
"""
if fg_cutoff is None:
fg_cutoff = self.fg_cutoff
if bg_cutoff is None:
Expand Down Expand Up @@ -124,6 +190,12 @@ def update_mask(self, fg_cutoff=None, bg_cutoff=None):
return

def update_background_plane(self, alpha=0.9):
"""
Update the background plane.
Parameters:
alpha (float): Alpha value for updating slope and intercept.
"""
y = self.counts
X = np.pad(self.difference_vectors, [[0, 0], [0, 1]], constant_values=1.0)
weights = np.reciprocal(self.counts) # Inverse variance poisson because #stats
Expand All @@ -139,6 +211,12 @@ def update_background_plane(self, alpha=0.9):
self.intercept = (1.0 - alpha) * self.intercept + alpha * intercept

def update_profile(self, alpha=0.9):
"""
Update the profile.
Parameters:
alpha (float): Alpha value for updating scale factor and profile centroid.
"""
weights = self.counts - self.background
weights = np.maximum(
weights, self.eps
Expand All @@ -158,6 +236,7 @@ def update_profile(self, alpha=0.9):
self.cen_y = (1.0 - alpha) * self.cen_y + alpha * cen_y

def integrate(self):
"""Integrate the profile."""
bg = self.background

self.I = ((self.counts - bg) * self.fg_mask).sum()
Expand All @@ -166,6 +245,14 @@ def integrate(self):

class SegmentedImage:
def __init__(self, pixels, centroids, radius=20):
"""
Initialize a SegmentedImage instance.
Parameters:
pixels (np.array): Pixel values.
centroids (np.array): Centroids of the pixels.
radius (int): Radius value.
"""
super().__init__()
self.radius = radius
pixels = np.array(pixels).astype("float32")
Expand Down Expand Up @@ -217,6 +304,16 @@ def __init__(self, pixels, centroids, radius=20):
self.distance_vectors = distance_vectors

def integrate(self, isigi_cutoff=2.0, knn=5):
"""
Integrate the SegmentedImage.
Parameters:
isigi_cutoff (float): I/SIGI cutoff value.
knn (int): Number of nearest neighbors.
Returns:
strong_idx (np.array): Strong indices.
"""
# Implement k nearest-neighbors for weak spots to average profile params from nearby strong spots
strong_idx = np.array(
[i for i, p in enumerate(self.profiles) if p.I / p.SigI > isigi_cutoff]
Expand Down
Loading

0 comments on commit d2133e6

Please sign in to comment.