Merge branch 'master' into feature/run_all_examples

examples/at_focused_annular_array_3D/at_focused_annular_array_3D.py

@@ -22,32 +22,32 @@
 verbose: bool = False
 
 # medium parameters
-c0: float            = 1500.0  # sound speed [m/s]
-rho0: float          = 1000.0  # density [kg/m^3]
+c0: float = 1500.0  # sound speed [m/s]
+rho0: float = 1000.0  # density [kg/m^3]
 
 # source parameters
-source_f0            = 1.0e6                                    # source frequency [Hz]
-source_roc           = 30e-3                                    # bowl radius of curvature [m]
-source_amp           = np.array([0.5e6, 1e6, 0.75e6])           # source pressure [Pa]
-source_phase         = np.deg2rad(np.array([0.0, 10.0, 20.0]))  # source phase [radians]
+source_f0 = 1.0e6  # source frequency [Hz]
+source_roc = 30e-3  # bowl radius of curvature [m]
+source_amp = np.array([0.5e6, 1e6, 0.75e6])  # source pressure [Pa]
+source_phase = np.deg2rad(np.array([0.0, 10.0, 20.0]))  # source phase [radians]
 
 # aperture diameters of the elements given an inner, outer pairs [m]
 diameters = np.array([[0.0, 5.0], [10.0, 15.0], [20.0, 25.0]]) * 1e-3
 diameters = diameters.tolist()
 
 # grid parameters
-axial_size: float   = 40e-3  # total grid size in the axial dimension [m]
+axial_size: float = 40e-3  # total grid size in the axial dimension [m]
 lateral_size: float = 45e-3  # total grid size in the lateral dimension [m]
 
 # computational parameters
-ppw: int             = 3      # number of points per wavelength
-t_end: float         = 40e-6  # total compute time [s] (this must be long enough to reach steady state)
-record_periods: int  = 1      # number of periods to record
-cfl: float           = 0.5    # CFL number
-source_x_offset: int = 20     # grid points to offset the source
-bli_tolerance: float = 0.01   # tolerance for truncation of the off-grid source points
-upsampling_rate: int = 10     # density of integration points relative to grid
-verbose_level: int   = 0      # verbosity of k-wave executable
+ppw: int = 3  # number of points per wavelength
+t_end: float = 40e-6  # total compute time [s] (this must be long enough to reach steady state)
+record_periods: int = 1  # number of periods to record
+cfl: float = 0.5  # CFL number
+source_x_offset: int = 20  # grid points to offset the source
+bli_tolerance: float = 0.01  # tolerance for truncation of the off-grid source points
+upsampling_rate: int = 10  # density of integration points relative to grid
+verbose_level: int = 0  # verbosity of k-wave executable
 
 # =========================================================================
 # RUN SIMULATION
@@ -161,7 +161,7 @@
 amp_on_axis = amp[:, Ny // 2]
 
 # define axis vectors for plotting
-x_vec = kgrid.x_vec[source_x_offset + 1 : -1, :] - kgrid.x_vec[source_x_offset]
+x_vec = kgrid.x_vec[source_x_offset + 1 :, :] - kgrid.x_vec[source_x_offset]
 y_vec = kgrid.y_vec
 
 # =========================================================================

examples/ivp_photoacoustic_waveforms/ivp_photoacoustic_waveforms.ipynb

@@ -37,6 +37,7 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
     "from kwave.data import Vector\n",
     "from kwave.kgrid import kWaveGrid\n",
     "from kwave.kmedium import kWaveMedium\n",
@@ -246,7 +247,6 @@
    "outputs": [],
    "source": [
     "t_sc, t_scale, t_prefix, _ = scale_SI(t_end)\n",
-    "import matplotlib.pyplot as plt\n",
     "_, ax1 = plt.subplots()\n",
     "ax1.plot(np.squeeze(kgrid2.t_array * t_scale), sensor_data_2D['p'] / np.max(np.abs(sensor_data_2D['p'])), 'r-')\n",
     "ax1.plot(np.squeeze(kgrid3.t_array * t_scale), sensor_data_3D[\"p\"] / np.max(np.abs(sensor_data_3D[\"p\"])), \"k-\", label=\"3D\")\n",

examples/na_controlling_the_pml/na_controlling_the_pml.ipynb

@@ -36,7 +36,6 @@
    "source": [
     "import os\n",
     "from copy import deepcopy\n",
-    "from tempfile import gettempdir\n",
     "from typing import Dict\n",
     "\n",
     "import numpy as np\n",

examples/us_bmode_linear_transducer/us_bmode_linear_transducer.ipynb

@@ -6,11 +6,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "try:\n",
-    "  import google.colab\n",
-    "  IN_COLAB = True\n",
-    "except:\n",
-    "  IN_COLAB = False"
+    "import os\n",
+    "import importlib.util\n",
+    "\n",
+    "def is_running_in_colab()->bool:\n",
+    "    return 'COLAB_RELEASE_TAG' in os.environ\n",
+    "\n",
+    "def is_library_installed(library_name):\n",
+    "    spec = importlib.util.find_spec(library_name)\n",
+    "    return spec is not None\n"
    ]
   },
   {
@@ -20,12 +24,11 @@
    "outputs": [],
    "source": [
     "%%capture\n",
-    "# install k-wave-python\n",
-    "try:\n",
-    "    import k_wave\n",
-    "except:\n",
-    "    !pip install k-wave-python\n",
-    "if  IN_COLAB:\n",
+    "\n",
+    "if not is_library_installed('kwave'):\n",
+    "    %pip install k-wave-python\n",
+    "\n",
+    "if is_running_in_colab():\n",
     "    !wget https://raw.githubusercontent.com/waltsims/k-wave-python/master/examples/us_bmode_linear_transducer/example_utils.py"
    ]
   },

examples/us_defining_transducer/us_defining_transducer.ipynb

@@ -21,17 +21,14 @@
     "\n",
     "from kwave.data import Vector\n",
     "from kwave.kgrid import kWaveGrid\n",
-    "from kwave.kmedium import kWaveMedium\n",
     "from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D\n",
     "from kwave.ksensor import kSensor\n",
     "from kwave.ktransducer import kWaveTransducerSimple, NotATransducer\n",
     "from kwave.kWaveSimulation import SimulationOptions\n",
     "\n",
     "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n",
-    "from kwave.utils.colormap import get_color_map\n",
     "from kwave.utils.dotdictionary import dotdict\n",
     "from kwave.utils.filters import spect\n",
-    "from kwave.utils.kwave_array import kWaveArray\n",
     "from kwave.utils.plot import voxel_plot\n",
     "from kwave.utils.signals import tone_burst"
    ]
@@ -239,7 +236,7 @@
     "axes[1].plot(f * 1e-6, as_3 / np.max(as_1.flatten()), 'r-', label = 'Sensor Position 3')\n",
     "axes[1].legend()\n",
     "axes[1].set_xlabel('Frequency [MHz]')\n",
-    "axes[1].set_ylabel('Normalised Amplitude Spectrum [au]');\n",
+    "axes[1].set_ylabel('Normalised Amplitude Spectrum [au]')\n",
     "f_max = medium.sound_speed / (2 * dx)\n",
     "axes[1].set_xlim([0, f_max * 1e-6])\n"
    ]

kwave/kWaveSimulation.py

@@ -160,8 +160,8 @@ def blank_sensor(self):
 
         """
         fields = ["p", "p_max", "p_min", "p_rms", "u", "u_non_staggered", "u_split_field", "u_max", "u_min", "u_rms", "I", "I_avg"]
-        if not any(self.record.is_set(fields)) and not self.time_rev:
-            return False
+        if not (isinstance(self.sensor, NotATransducer) or any(self.record.is_set(fields)) or self.time_rev):
+            return True
         return False
 
     @property
@@ -1311,7 +1311,7 @@ def create_sensor_variables(self) -> None:
         """
         # define the output variables and mask indices if using the sensor
         if self.use_sensor:
-            if not self.blank_sensor or isinstance(self.options.save_to_disk, str):
+            if not self.blank_sensor or self.options.save_to_disk:
                 if self.cuboid_corners:
                     # create empty list of sensor indices
                     self.sensor_mask_index = []

kwave/kWaveSimulation_helper/display_simulation_params.py

@@ -13,7 +13,7 @@ def display_simulation_params(kgrid: kWaveGrid, medium: kWaveMedium, elastic_cod
     k_size = kgrid.size
 
     # display time step information
-    logging.log(logging.INFO, "  dt: ", f"{scale_SI(dt)[0]}s, t_end: {scale_SI(t_array_end)[0]}s, time steps:", Nt)
+    logging.log(logging.INFO, f"  dt: {scale_SI(dt)[0]} s, t_end: {scale_SI(t_array_end)[0]}s, time steps: {Nt}")
 
     c_min, c_min_comp, c_min_shear = get_min_sound_speed(medium, elastic_code)
 
@@ -58,40 +58,34 @@ def print_grid_size(kgrid, scale):
     grid_scale_str = " by ".join(map(str, grid_size_scale))
 
     # display grid size
-    logging.log(logging.INFO, f"  input grid size: {grid_size_str} grid points ({grid_scale_str} m)")
+    logging.log(logging.INFO, f"  input grid size: {grid_size_str} grid points ({grid_scale_str}m)")
 
 
 def print_max_supported_freq(kgrid, c_min):
     # display the grid size and maximum supported frequency
     k_max, k_max_all = kgrid.k_max, kgrid.k_max_all
 
+    def max_freq_str(kfreq):
+        scaled_num_str, _, _, _ = scale_SI(kfreq * c_min / (2 * np.pi))
+        return scaled_num_str
+
     if kgrid.dim == 1:
         # display maximum supported frequency
-        logging.log(logging.INFO, "  maximum supported frequency: ", scale_SI(k_max_all * c_min / (2 * np.pi))[0], "Hz")
+        logging.log(logging.INFO, f"  maximum supported frequency: {max_freq_str(k_max_all)}Hz")
 
     elif kgrid.dim == 2:
         # display maximum supported frequency
         if k_max.x == k_max.y:
-            logging.log(logging.INFO, "  maximum supported frequency: ", scale_SI(k_max_all * c_min / (2 * np.pi))[0], "Hz")
+            logging.log(logging.INFO, f"  maximum supported frequency: {max_freq_str(k_max_all)}Hz")
         else:
-            logging.log(
-                logging.INFO,
-                "  maximum supported frequency: ",
-                scale_SI(k_max.x * c_min / (2 * np.pi))[0],
-                "Hz by ",
-                scale_SI(k_max.y * c_min / (2 * np.pi))[0],
-                "Hz",
-            )
+            logging.log(logging.INFO, f"  maximum supported frequency: {max_freq_str(k_max.x)}Hz by {max_freq_str(k_max.y)}Hz")
 
     elif kgrid.dim == 3:
         # display maximum supported frequency
         if k_max.x == k_max.z and k_max.x == k_max.y:
-            logging.log(logging.INFO, "  maximum supported frequency: ", f"{scale_SI(k_max_all * c_min / (2*np.pi))[0]}Hz")
+            logging.log(logging.INFO, f"  maximum supported frequency: {max_freq_str(k_max_all)}Hz")
         else:
             logging.log(
                 logging.INFO,
-                "  maximum supported frequency: ",
-                f"{scale_SI(k_max.x * c_min / (2*np.pi))[0]}Hz by "
-                f"{scale_SI(k_max.y * c_min / (2*np.pi))[0]}Hz by "
-                f"{scale_SI(k_max.z * c_min / (2*np.pi))[0]}Hz",
+                f"  maximum supported frequency: {max_freq_str(k_max.x)}Hz by {max_freq_str(k_max.y)}Hz by {max_freq_str(k_max.z)}Hz",
             )

kwave/kWaveSimulation_helper/set_sound_speed_ref.py

@@ -36,7 +36,7 @@ def get_ordinary_sound_speed_ref(medium):
 
     """
     c_ref = _get_sound_speed_ref(medium.sound_speed_ref, medium.sound_speed)
-    logging.log(logging.INFO, "  reference sound speed: ", c_ref, "m/s")
+    logging.log(logging.INFO, f"  reference sound speed: {c_ref} m/s")
     return c_ref, None, None
 
 
@@ -51,7 +51,7 @@ def get_pstd_elastic_sound_speed_ref(medium: kWaveMedium):  # pragma: no cover
 
     """
     c_ref = _get_sound_speed_ref(medium.sound_speed_ref, medium.sound_speed_compression)
-    logging.log(logging.INFO, "  reference sound speed: ", c_ref, "m/s")
+    logging.log(logging.INFO, f"  reference sound speed: {c_ref} m/s")
     return c_ref, None, None
 
 
@@ -66,10 +66,10 @@ def get_kspace_elastic_sound_speed_ref(medium: kWaveMedium):  # pragma: no cover
 
     """
     c_ref_compression = _get_sound_speed_ref(medium.sound_speed_ref_compression, medium.sound_speed_compression)
-    logging.log(logging.INFO, "  reference sound speed (compression): ", c_ref_compression, "m/s")
+    logging.log(logging.INFO, f"  reference sound speed (compression): {c_ref_compression} m/s")
 
     c_ref_shear = _get_sound_speed_ref(medium.sound_speed_ref_shear, medium.sound_speed_shear)
-    logging.log(logging.INFO, "  reference sound speed (shear): ", c_ref_shear, "m/s")
+    logging.log(logging.INFO, f"  reference sound speed (shear): {c_ref_shear} m/s")
 
     return None, c_ref_compression, c_ref_shear
 
@@ -86,5 +86,5 @@ def _get_sound_speed_ref(reference, speed):
     else:
         c_ref = reductions["max"](speed)
 
-    logging.log(logging.INFO, "  reference sound speed: ", c_ref, "m/s")
+    logging.log(logging.INFO, f"  reference sound speed: {c_ref} m/s")
     return float(c_ref)

kwave/kmedium.py

@@ -145,7 +145,7 @@ def _check_absorbing_without_stokes(self) -> None:
         assert np.isscalar(self.alpha_power), "medium.alpha_power must be scalar."
 
         # check y is real and within 0 to 3
-        assert np.all(np.isreal(self.alpha_coeff)) and 0 < self.alpha_power < 3, "medium.alpha_power must be a real number between 0 and 3."
+        assert np.all(np.isreal(self.alpha_coeff)) and 0 <= self.alpha_power < 3, "medium.alpha_power must be a real number between 0 and 3."
 
         # display warning if y is close to 1 and the dispersion term has not been set to zero
         if self.alpha_mode != "no_dispersion":

kwave/utils/kwave_array.py

@@ -81,7 +81,7 @@ def __post_init__(self):
                 self.end_point = [self.end_point]
             self.end_point = np.array(self.end_point, dtype=float)
 
-        self.measure = round(float(self.measure), 8)
+        self.measure = float(self.measure)
 
 
 class kWaveArray(object):
@@ -586,7 +586,7 @@ def get_off_grid_points(self, kgrid, element_num, mask_only):
             )
 
             # keep points in the positive y domain
-            grid_weights = grid_weights[:, kgrid.Ny:]
+            grid_weights = grid_weights[:, kgrid.Ny :]
 
         else:
             # remove integration points which are outside grid

kwave/utils/matrix.py

@@ -143,7 +143,7 @@ def expand_matrix(
         opts["pad_width"] = exp_coeff
     elif len(matrix.shape) == 2:
         if n_coeff == 2:
-            opts["pad_width"] = exp_coeff
+            opts["pad_width"] = [(exp_coeff[0],), (exp_coeff[1],)]
         if n_coeff == 4:
             opts["pad_width"] = [(exp_coeff[0], exp_coeff[1]), (exp_coeff[2], exp_coeff[3])]
     elif len(matrix.shape) == 3:

pyproject.toml

@@ -24,14 +24,14 @@ classifiers = [
     "Programming Language :: Python :: 3",
 ]
 dependencies = [
-    "h5py==3.11.0",
+    "h5py==3.12.1",
     "scipy==1.13.1",
     "opencv-python==4.10.0.84",
-    "deepdiff==7.0.1",
-    "matplotlib==3.9.0",
+    "deepdiff==8.0.1",
+    "matplotlib==3.9.2",
     "numpy>=1.22.2,<1.27.0",
     "beartype==0.18.5",
-    "jaxtyping==0.2.31"
+    "jaxtyping==0.2.34"
 ]
 
 [project.urls]
@@ -44,14 +44,15 @@ Bug-tracker = "https://github.com/waltsims/k-wave-python/issues"
 test = ["pytest",
         "coverage==7.6.1",
         "phantominator",
+        "testfixtures==8.3.0",
         "requests==2.32.3"]
 example = ["gdown==5.2.0"]
-docs = [ "sphinx-mdinclude==0.6.1",
+docs = [ "sphinx-mdinclude==0.6.2",
     "sphinx-copybutton==0.5.2",
     "sphinx-tabs==3.4.5",
-    "sphinx-toolbox==3.6.0",
-    "furo==2024.5.6"]
-dev = ["pre-commit==3.7.1"]
+    "sphinx-toolbox==3.8.0",
+    "furo==2024.8.6"]
+dev = ["pre-commit==3.8.0"]
 
 [tool.hatch.version]
 path = "kwave/__init__.py"

tests/matlab_test_data_collectors/python_testers/utils/check_equality.py

@@ -1,5 +1,6 @@
 import functools
 import logging
+import dataclasses
 
 import numpy as np
 
@@ -13,6 +14,25 @@ def _getattr(obj, attr):
     return functools.reduce(_getattr, [obj] + attr.split("."))
 
 
+def check_element_equality(
+    actual_element: Element,
+    expected_element: Element,
+) -> bool:
+    for field in dataclasses.fields(expected_element):
+        expected = getattr(expected_element, field.name)
+        actual = getattr(actual_element, field.name)
+        if isinstance(expected, np.ndarray):
+            if not np.allclose(actual, expected):
+                return False
+        elif isinstance(expected, float):
+            if not np.isclose(actual, expected):
+                return False
+        else:
+            if actual != expected:
+                return False
+    return True
+
+
 def check_kgrid_equality(kgrid_object: kWaveArray, expected_kgrid_dict: dict):
     are_totally_equal = True
     for key, expected_value in expected_kgrid_dict.items():
@@ -107,14 +127,21 @@ def check_kwave_array_equality(kwave_array_object: kWaveArray, expected_kwave_ar
         actual_value = recursive_getattr(kwave_array_object, mapped_key, None)
 
         if key == "elements":
+            are_equal = True
             if isinstance(expected_value, dict):
                 expected_value = [Element(**expected_value)]
-            elif isinstance(expected_value, np.ndarray):
-                expected_value = expected_value.tolist()
+                for actual, expected in zip(actual_value, expected_value):
+                    are_equal &= check_element_equality(
+                        actual_element=actual,
+                        expected_element=expected,
+                    )
             elif isinstance(expected_value, list):
                 expected_value = [Element(**val) for val in expected_value]
-            # Hacky way to compare but works the best for now
-            are_equal = str(actual_value) == str(expected_value)
+                for actual, expected in zip(actual_value, expected_value):
+                    are_equal &= check_element_equality(
+                        actual_element=actual,
+                        expected_element=expected,
+                    )
         else:
             actual_value = np.squeeze(actual_value)
             expected_value = np.array(expected_value)