Update record_reader.py to remove warnings (#211)

* Update record_reader.py to remove warnings

* Update beamform.py

* Update simulation_options.py

* clean up zero division

* clean up minor warnings

* add N assertion

* fix element-wise string comparison warning

* fix num_radial

* try something for get_win

* add back d_radial

* comment out divide by zero assertation

* check for power exponent of 1

* test new assertion

* revert changes

* modify assertion

* remove fixes that don't help

* assert Z is greater than 0

* catch zero distance case

* attemt to remove warning with another assertion

* Revert fix. discuss matlab behaviour

* treat Z as an array

* assert N > 1

* Try to solve divide by zero warnings in signals

* Add TODO

* revert assertion

* update value checking

* update assertion errors

* change assertion the validation

---------

Co-authored-by: Farid Yagubbayli <[email protected]>

kwave/options/simulation_options.py

@@ -285,7 +285,7 @@ def option_factory(kgrid: "kWaveGrid", options: SimulationOptions):
                 options.plot_scale = [-1, 1]
 
         # replace defaults with user defined values if provided and check inputs
-        if (val := options.pml_alpha) is not None and options.pml_alpha != 'auto':
+        if (val := options.pml_alpha) is not None and not isinstance(options.pml_alpha, str):
             # check input is correct size
             val = np.atleast_1d(val)
             if val.size > kgrid.dim:

kwave/reconstruction/beamform.py

@@ -30,7 +30,8 @@ def focus(kgrid, input_signal, source_mask, focus_position, sound_speed):
          input_signal_mat:  matrix of time series following the source points
     """
 
-    assert kgrid.t_array != 'auto', "kgrid.t_array must be defined."
+    assert not isinstance(kgrid.t_array, str), "kgrid.t_array must be a numeric array."
+
     if isinstance(sound_speed, int):
         sound_speed = float(sound_speed)
 

kwave/utils/atten_comp.py

@@ -46,7 +46,7 @@ def atten_comp(
         dt: time step [s]
         c: sound speed [m/s]
         alpha_0: power law absorption prefactor [dB/(MHz^y cm)]
-        y: power law absorption exponent [0 < y < 3, y ~= 1]
+        y: power law absorption exponent [0 < y < 3, y != 1]
         display_updates: Boolean controlling whether command line updates
 %       and compute time are printed to the command line
         distribution: default TF distribution
@@ -75,6 +75,9 @@ def atten_comp(
     # extract signal characteristics
     N, num_signals = signal.shape
 
+    if y == 1:
+        raise ValueError("A power exponent [y] of 1 is not valid.")
+
     # convert absorption coefficient to nepers
     alpha_0 = db2neper(alpha_0, y)
 
@@ -252,15 +255,22 @@ def findClosest(arr: NDArray[Shape["Dim1"], Float], value: float):
     dist_vec = c * dt * (np.arange(N) - t0)
     dist_vec[dist_vec < 0] = 0
 
-    # create the time variant filter
+    # Check if f_array and dist_vec are valid
+    assert f_array is not None and len(f_array) > 0, "f_array must have non-zero length."
+    assert dist_vec is not None and len(dist_vec) > 0, "dist_vec must have non-zero length."
+
+    # Create the time variant filter
     f_mat, dist_mat = np.meshgrid(f_array, dist_vec)
+
+    assert y != 1, "A power exponent [y] of 1 is not valid." # this is a duplicate assertion to attempt to remove a warning
+
+    # Add conditionals or use np.where to manage zero and NaN
     part_1 = (2 * np.pi * np.abs(f_mat)) ** y
     part_2 = 1j * np.tan(np.pi * y / 2)
     part_3 = (2 * np.pi * f_mat)
     part_4 = (2 * np.pi * np.abs(f_mat)) ** (y - 1)
-    tv_filter = alpha_0 * dist_mat * (
-            part_1 - part_2 * part_3 * part_4
-    )
+
+    tv_filter = alpha_0 * dist_mat * (part_1 - part_2 * part_3 * part_4)
 
     # convert cutoff frequency to a window size
     N_win_array = np.floor((cutoff_freq_array / f_array[-1]) * N) - 1

kwave/utils/mapgen.py

@@ -64,7 +64,11 @@ def make_spiral_circle_points(num_points: int, radius: float) -> np.ndarray:
         return p0
 
     def make_concentric_circle_points(num_points: int, radius: float) -> Tuple[np.ndarray, int]:
+
+        assert num_points >= 1, "The number of points must be greater or equal to 1." 
+
         num_radial = int(np.ceil(np.sqrt(num_points / np.pi)))
+
         try:
             d_radial = radius / (num_radial - 1)
         except ZeroDivisionError:
@@ -217,7 +221,7 @@ def make_cart_bowl(
     bowl = R @ p0 + b
 
     # plot results
-    if plot_bowl:
+    if plot_bowl is True:
         # select suitable axis scaling factor
         _, scale, prefix, unit = scale_SI(np.max(bowl))
 
@@ -580,6 +584,9 @@ def make_cart_sphere(
     r = np.sqrt(1 - (y ** 2))
     phi = k * inc
 
+    if num_points <= 0:
+        raise ValueError("num_points must be greater than 0")
+
     # create the sphere
     sphere = radius * np.concatenate([np.cos(phi) * r[np.newaxis, :], y[np.newaxis, :], np.sin(phi) * r[np.newaxis, :]])
 
@@ -2831,6 +2838,8 @@ def calculate_axial_pressure() -> Tuple[
     def calculate_lateral_pressure() -> NDArray[Shape["N"], Float]:
         # calculate magnitude of the lateral pressure at the geometric focus
         Z = k * lateral_positions * diameter / (2 * radius)
+        # TODO: this should work
+        # assert np.all(Z) > 0, 'Z must be greater than 0'
         lateral_pressure = 2. * density * sound_speed * velocity * k * h * scipy.special.jv(1, Z) / Z
 
         # replace origin with limit
@@ -3098,7 +3107,7 @@ def make_cart_spherical_segment(
     segment = R @ p0 + b
 
     # plot results
-    if plot_bowl:
+    if plot_bowl is True:
         _, scale, prefix, unit = scale_SI(np.max(segment))
 
         # create the figure

kwave/utils/math.py

@@ -410,8 +410,16 @@ def compute_linear_transform(pos1, pos2, offset=None):
     # Compute vector pointing from pos1 to pos2
     beam_vec = pos2 - pos1
 
+    magnitude = np.linalg.norm(beam_vec)
+
+    # TODO: it seems matlab behaviour is to return nans when positions are the same.
+    #       we should open an issue and change our behaviour once matlab is fixed.
+    # if magnitude == 0:
+    #     # "pos1 and pos2 are the same"
+    #     return np.eye(3), 0 if offset is None else offset
+
     # Normalise to give unit beam vector
-    beam_vec = beam_vec / np.linalg.norm(beam_vec)
+    beam_vec = beam_vec /  magnitude
 
     # Canonical normalised beam_vec (canonical pos1 is [0, 0, 1])
     beam_vec0 = np.array([0, 0, -1])

kwave/utils/signals.py

@@ -146,8 +146,12 @@ def cosine_series(n: int, N: int, coeffs: List[float]) -> np.ndarray:
 
     # create the window
     if N.size == 1:
+        # TODO: what should this behaviour be if N is a list of ints? make windows of multiple lengths?
         n = np.arange(0, N)
 
+        # TODO: find failure cases in test suite when N is zero.
+        # assert np.all(N) > 1, 'Signal length N must be greater than 1'
+
         if type_ == 'Bartlett':
             win = (2 / (N - 1) * ((N - 1) / 2 - abs(n - (N - 1) / 2))).T
         elif type_ == 'Bartlett-Hanning':

tests/matlab_test_data_collectors/python_testers/utils/record_reader.py

@@ -1,12 +1,10 @@
 import numpy as np
-
 from scipy.io import loadmat
 
 
 class TestRecordReader(object):
     # Will make `pytest` to ignore this class as a test class
     __test__ = False
-
     def __init__(self, record_filename):
         recorded_data = loadmat(record_filename, simplify_cells=True)
         self._records = recorded_data